{'fraction_non_alphanumeric': 0.030886157826649416, 'fraction_numerical': 0.007761966364812419, 'mean_word_length': 4.836353340883353, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
PSpectRe: A Pseudo-Spectral Code for (P)reheating

PSpectRe is a C++ program that uses Fourier-space pseudo-spectral methods to evolve interacting scalar fields in an expanding universe. PSpectRe is optimized for the analysis of parametric resonance in the post-inflationary universe, and provides an alternative to finite differencing codes, such as Defrost and LatticeEasy. PSpectRe has both second- (Velocity-Verlet) and fourth-order (Runge-Kutta) time integrators. Given the same number of spatial points and/or momentum modes, PSpectRe is not significantly slower than finite differencing codes, despite the need for multiple Fourier transforms at each timestep, and exhibits excellent energy conservation. Further, by computing the post-resonance equation of state, we show that in some circumstances PSpectRe obtains reliable results while using substantially fewer points than a finite differencing code. PSpectRe is designed to be easily extended to other problems in early-universe cosmology, including the generation of gravitational waves during phase transitions and pre-inflationary bubble collisions. Specific applications of this code will be pursued in future work.


Introduction
A key challenge for any inflationary model is to ensure that inflation ends and that the post-inflationary universe thermalizes, setting the stage for a conventional, hot-big-bang cosmology. A promising approach to this problem is provided by preheating and parametric resonance, in which the post-inflationary universe undergoes explosive particle production [1,2,3]. Preheating leads to a fascinating range of phenomenological consequences, possibly including superheavy dark matter [4], gravitational-wave production [5,6,7,8,9,10,11] and non-Gaussianity [12,13,14]. Resonance can occur in a broad swathe of cosmological models and, since it depends on the coupling between the inflaton sector and the "rest of physics," any observable consequences of resonance constitute a unique window into the inflationary epoch.
Given that resonance is an intrinsically nonlinear process, understanding it is a challenging task, and frequently requires numerical studies. Consequently, expanding on a theoretical framework initially developed by Khlebnikov and Tkachev [15], a number of codes have been developed to study scalar field dynamics in an expanding universe, including LatticeEasy [16] (not to mention the parallelized ClusterEasy [17], and GPU accelerated CUDAEASY [18]) and Defrost [19]. These programs all use finite differencing methods to compute spatial derivatives of the scalar fields. In this paper we describe PSpectRe, a Fourier-space spectral code for solving the Klein-Gordon equations of two scalar fields in an expanding universe with an arbitrary quartic interaction potential. Technically, this is a pseudo-spectral code, in that nonlinear terms in the potential and its derivatives are computed by first converting the fields into their position-space representation, performing the necessary multiplications, and then taking the inverse transform [20]. Since spatial derivatives are taken by multiplying the corresponding Fourier component by the wavenumber k i , pseudo-spectral methods are free of differencing noise 1 , and have excellent spatial resolution [20]. On the other hand, overhead from the Fourier transforms and inverse transforms exacts a cost in terms of computational efficiency. As currently implemented, PSpectRe allows both second-and fourth-order-accurate time integration steps. The latter is more computationally intensive, but delivers excellent energy conservation. PSpectRe is currently implemented to run on a single, shared memory machine, but threads efficiently across 16 (effective) cores. The rate limiting step is the fast Fourier transform (FFT) algorithm. It would be straightforward to port PSpectRe to a cluster environment, given that the FFT libraries we use exist in a fully parallelized form.
The goal of this paper is to introduce and describe PSpectRe, and to demonstrate that it provides a useful cross check on existing finite-differencing-based codes. Looking to the future, spectral codes have excellent accuracy and precision, and the investigation of preheating-generated observables such as gravitational waves and non-Gaussianity require accurate solutions of the background dynamics. The initial version of PSpectRe is equivalent to LatticeEasy and Defrost in that it solves the Klein-Gordon equation in a rigid, expanding background. However, Easther, Giblin and Lim have previously used a spectral algorithm to compute the gravitational-wave background generated during resonance [7,8]. This approach naturally lends itself to implementation within PSpectRe, and will be pursued in future work. Finally, in earlier work Easther and Parry used a pseudo-spectral code to study gravitational effects on resonance in a 1+1 dimensional scenario [21,22]. In terms of performance, a "stock" version of Defrost is a few times faster than PSpectRe on the same hardware (with a 256 3 grid), but the difference is not dramatic. However, we will compute the post-inflationary equation of state, and show that PSpectRe can achieve reliable results with a 64 3 grid in circumstances where Defrost requires 128 3 or even 256 3 points. We attribute this improvement to the lack of differencing noise in spectral codes.
The format of this paper is as follows. We first review the scalar field evolution equations, and the specific program variables used within PSpectRe, which are the same as those employed by LatticeEasy. We then describe PSpectRe's basic computational 1 By differencing noise we mean the error introduced when computing approximations to a derivative by subtracting neighboring-point values and then dividing by some multiple of the grid spacing. Since neighboring points often have similar values, subtracting them using finite-precision, floating-point arithmetic often yields a result with few significant digits. When the result is divided by some multiple of the grid spacing, the rounding error is greatly magnified, resulting in a significant source of numerical error. engine, and present results on energy conservation and the spatial distribution of energy following resonance for the well-studied scenario of an inflation, φ, with a quadratic potential, coupled to a massless field, χ, via a φ 2 χ 2 interaction term.

Equations of Motion
We consider a homogeneous and isotropic universe with a flat Friedmann-Robertson-Walker metric The scale factor a(t) obeys The universe contains an inflaton field φ(t, x, y, z) and a generic matter field χ(t, x, y, z) which obey the Klein-Gordon equation, PSpectRe supports the generic quartic potential (2.6) For this system of interacting scalar fields, the density and pressure are Consider scalar fields whose values and time derivatives are known at the points of an N × N × N rectangular lattice, which represents the spatial part of an FRW universe with periodic spatial boundary conditions. The same information is contained in the discrete spatial Fourier transform, and similarly for χ and X( k). Here, k labels grid points as wave vectors and r labels the grid points as spatial positions. The inverse transformation is The Φ( k) are, in general, complex-valued but, since φ and χ are real valued, φ( r) = φ( r) , so Φ( k) = Φ(− k) and the number of free parameters matches in both representations.

PSpectRe: Algorithm and Implementation
Working in momentum space, the field evolution is described by N 3 second-order, ordinary differential equations for the time dependence of each of Φ( k) and X( k), two further degrees of freedom for the scale factor (equation 2.3), subject to the overall Hamilton constraint, (2.2). Spatial derivatives do not appear, since ∂ i φ applied to equation (2.9) simply multiplies Φ( k) by k i , the usual result for the Fourier transform of a derivative. Consequently, while finite difference codes work by discretizing the spatial derivatives, (psuedo)-spectral methods simply replace ∇ 2 with k · k. These equations are nonlinear, however, via the interaction terms in the potential, and the "friction" terms in the Klein-Gordon equation.
To avoid confusion, we point out there is a second class of spectral methods for differential equations derived by discretizing the continuum equations and then expanding φ(x + δx) as a discrete Fourier series. This method is used, for example, to solve a linear differential equation by transforming the problem into a matrix equation for the Fourierseries coefficients [20]. While such methods can produce precision-limited solutions to the discretized equations, those solutions are, by definition, only approximate solutions to the continuum equations, since the field derivatives have been truncated to some finite order. The method used in PSpectRe is distinct from those which use the discretizefirst methodology, in that we make the transition to Fourier space before discretization. This has the benefit of representing the derivative operators exactly, but requires the following negligibility assumption. Namely, there is a multivalued relationship between the discrete-Fourier-series coefficients and the continuum momentum values: a discrete mode k corresponds not only to the continuum mode k, but also to the continuum mode k − 2πN L , where N is the total number of points per grid side and L is the physical length of each side of the simulated box. This is because contributions to those two modes are indistinguishable on a finite grid of regularly-spaced samples. PSpectRe uses the convention that the first N 2 + 1 Fourier-space components in any dimension represent the modes 0, . . . , πN L and the remaining N 2 − 1 points represent the modes − π(N −2) L through − 2π L . By simulating systems for which the modes π(N +2) L , PSpectRe's representation achieves high fidelity by exactly computing the solution for those that are modes explicitly represented. Importantly, the dispersion relation induced by the Laplacian operator is not distorted from its continuum behavior, as it would be if a discretize-first methodology were employed.
In the initial phases, resonance typically amplifies a small subset of Fourier modes, and in this scenario a Fourier basis is particularly natural. Conversely, one might worry that a pseudo-spectral code would be more likely to fail when applied to scenarios in which the field has substantial, localized peaks in position space which is possible at late times, or in scenarios with topological defects. In work currently in progress, we have used PSpectRe to analyze oscillon formation and evolution in three dimensions (c.f. [23]. As will be described in a forthcoming paper, our experience is that the pseudo-spectral code copes well with very sharply peaked field configurations, even at effective spatial resolutions where finite differencing schemes would fail.

Field Rescalings and the Velocity Terms
The time evolution can be handled by any suitable ordinary-differential-equation solver, and PSpectRe implements both a second-order (in time) Velocity-Verlet routine [24], and a fourth-order Runge-Kutta integrator [25]. To improve the stability of the integration, we remove the Hφ and Hχ terms from the equations of motion, making use of the same rescalings as LatticeEasy. 2 We quote them here for convenience If λ φ = 0 and λ χ is 0, we set where φ 0 is the initial value of the amplitude of the k = 0 mode for φ. In all cases, since the φ field is initially dominant. We also follow the LatticeEasy convention of using a prime to denote differentiation with respect to the rescaled program time, while dots denote differentiation with respect to physical time. Rescaled (or "program") variables are written with a subscript "pr".

Transformed Equations of Motion
To write the equations of motion for the Fourier coefficients of the fields, we take the Fourier transform of the Klein-Gordon equation. This results in the following equation for the amplitude of each momentum mode with wave vector k : where the hat above the partial derivative indicates a Fourier transform. The equation for X(k) is of precisely the same form.
In terms of rescaled variables, the equations of motion are The expression for ∂Vpr ∂Xpr is similar, with Φ pr and X pr exchanging roles. The variable rescalings are designed so that s − 2r + 3 ≡ 0, so the Φ pr dependence is eliminated.
The scale factor a is evolved using the same version of the Friedmann equations employed by LatticeEasy: where angle brackets denote averages over the simulation box. In terms of rescaled variables, this becomes

Nonlinear Terms
Spectral codes are simple to implement in linear systems, because the individual momentum modes do not mix. However, preheating and resonance are driven by the nonlinear couplings between modes, and properly accounting for these is the primary challenge faced by any spectral treatment of resonance, both in terms of algorithmic complexity and execution time.
We have to compute several different nonlinear terms. Firstly, to evolve the scale factor we need averages over the entire box for both quadratic (e.g. |∇φ| 2 ) expressions and the potential, which -in principle -can contain mixtures of arbitrary order in the fields. Writing a variable (e.g. φ) in terms of its Fourier expansion and integrating over the box generates delta functions, and we evaluate the corresponding terms as a sum of squares: Integrating the potential is trickier, since we now have cubic and higher-order terms. In principle, we could evaluate an M − th order term by writing Φ( k 1 )Φ( k 2 ) · · · Φ( k M ); performing the spatial integral gives a delta function which ensures that k 1 + k 2 +· · ·+ k M = 0, and we could then sum over all the corresponding momenta that satisfied this expression. However, the number of possible combinations is prohibitively expensive as M becomes large. Consequently, we shift the fields into the position space representation, compute the potential, and integrate over the (spatial) box using Simpson's rule. This approach is necessarily inaccurate, since the product of the Φ( k 1 )Φ( k 3 ) · · · Φ( k M ) will contain terms of frequency up to ∼ M k max , which will not be resolved on our spatial grid. In other words, if Φ contains a Fourier-space component e ikmaxx , then Φ M will contain a contribution e ikmaxx M = e i(M kmax)x . In practice, power in such off-grid modes is automatically aliased by the DFT to some mode represented on the grid. PSpectRe contains an option to "pad" the spatial grid for the computation of the energy integral, so that it contains (pN ) 3 points, where p is a (small) integer. The padding and unpadding algorithms are specified in the appendix. This slows the code, but dramatically improves instantaneous energy conservation, as we discuss below.
To check the energy conservation we need to compute the average density, ρ over the simulation volume. Quadratic terms are simply dealt with via Parseval's Theorem. In terms of rescaled variables, ρ is given by We use four applications per field of Parseval's theorem to compute every term in the energy density except for the potential: Finally, we need the Fourier components of the potential derivatives, ∂V ∂φ and ∂V ∂χ , which are generally nonlinear with respect to the fields, including terms proportional to χ 2 φ, φ 2 χ, φ 3 and χ 3 . We handle these terms by transforming the fields into position space, performing the multiplication (such as χ 2 φ), and then transforming the result back to Fourier space to proceed with the spectral evolution. This is necessarily time-consuming, and is the price we pay for evolving a nonlinear system of equations with a spectral algorithm.

Integrators
PSpectRe provides two user-selectable integrators: A standard fourth-order-in-time Runge-Kutta integrator and a second-order-in-time Velocity-Verlet integrator. The Runge-Kutta scheme is more expensive, since the nonlinear terms in the potential must be transformed four times per time step. The Velocity-Verlet scheme is less expensive, but only secondorder accurate in time. The scale factor is evolved along with the fields using the same integration scheme. In the case of the Velocity-Verlet scheme, this violates one of the scheme's assumptions, that the second derivative does not directly depend on the first derivative, sinceä explicitly depends onȧ. Sinceȧ has become small and slowly-varying by the time resonance sets in, the effect of this discrepancy is minor, but reduces the Verlet integrator to first order accuracy for the scalar factor.

Initial conditions
PSpectRe provides two routines for initializing the fields, which emulate the procedures of LatticeEasy and Defrost, respectively. In either case, the initial value for the k = 0 mode is provided by the user in physical units, and sets to the mean value of the corresponding field. In our experience, the post-resonance state of the systems we have simulated are qualitatively, if not quantitatively, unaffected by the choice of the initialconditions routine.

LatticeEasy Emulation
We refer the reader to the LatticeEasy documentation (section 6.3) for a complete discussion of the assumptions and approximations that underlie the following formulae. The initial values of all mode amplitudes except the k = 0 mode are assumed to be given by Rayleigh distributions with RMS values W k such that where (3.14) In terms of rescaled variables, the RMS values are where (3.16) Once we have computed these RMS amplitudes, we generate a uniform deviate Y and use it to compute a randomized quantity f kpr : Finally, we generate two random phases, θ l and θ r , for the left-moving and right-moving components of the momentum mode, and set the real and imaginary parts of its amplitude via Φ kpr = 1 √ 2 (f kpr )(e iθ l + e iθr ) . These initial values for the momentum-mode amplitudes must be multiplied by a factor of N 3 in order to be consistent with the Fourier transform conventions listed above. Since the fields are real-valued, we must ensure that these initial momentum mode amplitudes obey the appropriate conjugate symmetry.
For simplicity, the initial values of the time-derivatives of the k = 0 mode amplitudes are chosen to be zero for both fields. The initial values of the derivatives for the other modes are set according to the formula where H pr ≡ a a .

Defrost Emulation
Defrost uses a convolution-based method to initialize the fields and their derivatives.
The initialization process is viewed as constructing a sample from a Gaussian random field defined asφ where the mode amplitudes are uncorrelated random variables: and The effective mass is This is implemented by convolving white noise with a (spherically symmetric) kernel function sin kr kr .
Unfortunately, this kernel function is singular as r → 0, so it is regularized by introducing a Gaussian cut-off at some scale q below the Nyquist frequency Unlike Defrost, which uses a cut-off of 1 2 k Nyquist , PSpectRe uses a cut-off of 150 2 k Nyquist . The larger cut-off produces field-mode amplitudes closer to their physical values. The regularized kernel function, ξ(r), is then convolved with the {b k }: A similar procedure is used to initialize the velocity fields.

Units
The masses and coupling-constants appearing in the potential and the field amplitudes are specified in physical units normalized such that the Planck mass is 1. The size of the box and the program time are specified in terms of the (rescaled) program units. For a potential dominated by an inflaton mass term, this is equivalent to specifying these quantities in terms of the inflaton mass. Specifically: (3.27)

Validation
PSpectRe was validated by comparing to LatticeEasy and Defrost using the m 2 φ 2model, and also by ensuring that its own results are self-consistent as the timestep, inte- grator and box size are varied. We primarily focus upon the quadratic model which is the default model in Defrost, and all runs shown here have the same parameter values as those presented in [19]. We also tested our code against the quartic λφ 4 model which is implemented inside LatticeEasy. Formally, there is no conserved energy in the system, but we gauge the accuracy of the evolution by checking the accuracy with which the averaged Friedmann equation is satisfied, plotting a quantity that will vanish (up to ∼ δρ ρ , initially ∼ 10 −7 ) when energy is "perfectly" conserved. Figure 1 shows the energy as a function of the integration scheme. Initially the Runge-Kutta greatly outperforms the Verlet, but both degrade as resonance hits. We compute the average potential in our box by constructing V (φ, χ) is real space and performing a numerical integral. However, for this case the shortest momentum modes have a nontrivial amplitude and the resulting integral effectively under-samples the "fine structure" in the potential. To test this hypothesis we have added a "padding" feature to the code, where the N 3 momentum modes are Fourier transformed back into a (pN ) 3 spatial lattice, which effectively interpolates extra points into the numerical integral. The padding and unpadding algorithms are specified in the appendix. Figure 2 shows the same scenarios as those plotted in Figure 1 with a padding factor p = 2. In this case, the energy conservation for both integrators is improved by an order of magnitude. Finally, Figure 3 Figure 1, for Runge-Kutta integrators and a box size of 64. The red plot is unpadded, and the blue plot has a padding factor of 2: we show the latter simulation on its own in the lower panel. In this case there is an initial transient which does not depend on the algorithm, and the subsequent evolution is followed at better than 1 part in 10 6 . runs with and without padding. With padding we obtain an accuracy of the order of parts in 10 6 throughout the run. Increasing N (for constant box size, L) will often have a similar effect to using padded grid for the potential, if the higher momentum modes have a low amplitude throughout the simulation. However, using the padded grid increases the run time by a factor of approximately 2.5, whereas doubling N will increase runtime by an order of magnitude, since the higher resolution simulation tracks eight (= 2 3 ) times as many momentum modes as the lower resolution simulation.
More generally, our results here reflect several paradoxes associated with the discussion of energy conservation in a preheating code. Firstly, these systems are intrinsically nonlinear, and mode-mode couplings transfer power to momenta that are unresolved within our simulation. Consequently, a code that has perfect energy conservation while tracking a finite number of momentum modes is unphysical to the extent that it suppresses this transfer of power, and effectively imposes a low-pass filter on the dynamics. Our use of a padded grid avoids this problem, since we can accurately compute the potential-energy integral. Note that we evaluate the quadratic expressions for the momentum and gradients exactly without resorting to padding -writing these spatial integrals with the appropriate Fourier expansions creates a sequence of δ functions and our integrals are evaluated by summing the squared-amplitudes of the Fourier coefficients. Finally, this code (like most other preheating simulations, see [28,29] for exceptions) implicitly ignores any contribution from metric perturbations to the scalar field evolution. Despite these caveats, we have demonstrated that PSpectRe can conserve the Hamiltonian constraint for this system to high accuracy, if desired. In practice, PSpectRe's Verlet integrator is suitable in most situations, and offers savings in both runtime and memory use. On the other hand, if high accuracy is needed, the Runge-Kutta integrator provides better performance than setting the Verlet timestep to a very small value. We compare the energy conservation of PSpectRe and Defrost in Figure 4. We find that we can make the timestep in PSpectRe somewhat larger than Defrost's "default" value and obtain a comparable level of energy conservation. For this case Defrost ran approximately twice as fast as PSpectRe, with both codes saturating a 16 (effective) core, shared-memory machine. In practice, performance can degrade substantially when "additional" variables are computed, such as the Newtonian gravitational potential, or the user calls for frequent data output. Given that the Fast Fourier Transform is an O(3N 3 log N ) process we would expect PSpectRe to run more slowly than a comparable finite differencing code, which naively scales as O(N 3 ), but in practice the time penalty is not exorbitant.
While PSpectRe evolves the fields in momentum space, it can output its results in terms of the spatial values of fields, and of the various classes of energy within the model. As an example. Figure 5 shows the output of a quadratic potential with L = 3 at 256 3 , and we recover the spatial patterns described in [27]. Parametric resonance is typically associated with a fixed physical scale, which determines the length scales that the simulation must accurately resolve. Clearly, algorithms which allow this characteristic scale to be close to the smallest length (or highest momentum) are more memory-efficient than those which require a substantial amount of separation between these two scales. As noted in the introduction, spectral codes compute spatial derivatives using an algebraic identity, rather than by comparing values of the fields at adjacent points, and are thus not affected by differencing errors. Consequently, it is reasonable to hope that PSpectRe might need fewer points (or momentum modes) than a finite-differencing code applied to the same system. To investigate this, we extract the instantaneous equation of state w as a function of time [19], where w = p ρ . (4.4) In Figure 6 we show the computed value of w for the quadratic model with box size L = 5, using Defrost at N = 128 and N = 256. The two curves almost completely overlap, showing that the code is effectively in the continuum limit when N ≥ 128. In Figure 7 we plot the same quantify at smaller N , and find that at N = 64 and N = 32 the resulting values for w differ substantially. Conversely, with PSpectRe we find excellent agreement between N = 128 and N = 64, while the N = 32 case is still largely reliable. Consequently, it appears that pseudo-spectral code does allow us to accurately follow the dynamics of resonance with a lower cutoff length scale. From a practical perspective, this can either have a substantial effect on the runtime, given that at a constant timestep, performance naively scales as N 3 , or (alternatively) PSpectRe can simulate larger spatial volumes than a finite differencing-code, given a fixed amount of computational resources. Finally, while PSpectRe and Defrost both reach a "continuum limit" when the number of points  , whereas the right panel shows the results of three PSpectRe simulations, at 32 3 (red), 64 3 (blue) and 128 3 (black). Other than the number of points, the parameters are identical between these runs and those in Figure 6. Note the vertical scales are different between the two plots, and the time interval shown is subset of that used in Figure 6.
being tracked in the simulation is sufficiently large, the two codes do not overlap with each other as the detailed time-evolution of w depends very sensitively on the implementation of the initial conditions. However, we do recover the late time limiting value w ∼ 0.2 seen by Frolov.

Conclusion
We have demonstrated that PSpectRe yields a highly-accurate and computationallyefficient approach to the solution of scalar-field evolution during parametric resonance.
PSpectRe thus provides a useful contrast to finite-differencing algorithms such as those in Defrost and LatticeEasy.
The present discussion is intended primarily as a "proof of concept", demonstrating that a pseudo-spectral algorithm is a practical method for the analysis of preheating, despite the extra computational overhead introduced by the nonlinear terms in the evolution equations. In doing so, we have recovered a variety of existing results in the literature. We find that PSpectRe's energy conservation performance is excellent, and scales in wellunderstood ways as the grid size and timestep are varied. Moreover, we demonstrate that some energy is "lost" to the simulation when it is transferred to high frequency modes. While PSpectRe does not evolve these modes, it can account for them by "padding" the spatial integral used to compute the potential energy.
As a pseudo-spectral code, PSpectRe does not rely on comparing values of fields at adjacent lattice points to take spatial derivatives. Consequently, in some circumstances PSpectRe provides accurate results with smaller grids than those needed by finite differencing codes. When the grid sizes are equal, PSpectRe is still competitive with finitedifferencing codes in terms of runtime, but the repeated Fourier transforms do introduce a noticeable performance penalty. Reducing the grid size by a factor of 2 naively improves runtime by (slightly more than) a factor of 8, with fixed timestep, while maintaining the same effective resolution as a finite-differencing code. Consequently, by changing runtime settings, PSpectRe can be optimized for precision, speed, or the ability to run simulations over large spatial volumes.

A. Download and License
PSpectRe is publicly available 3 , and is made available under a BSD-style license. Its current features are on a par with those of Defrost and LatticeEasy, and we plan to extend its capabilities in future work. PSpectRe compiles with both the freely-available GNU C++ compiler and the Intel C++ compiler. The code calls either the FFTW library, version 3 [26] or Intel's Math Kernel Library to compute the Fourier transforms Our experience is that the Intel library outperformed FFTW on machines where both were available. If an extended precision FFTW library is available 4 , the user may select a mode in which all computation is done using extended floating-point precision. PSpectRe has been run on both Mac OS and Linux based systems. Finally, we warn that not all combinations of features and settings have been exhaustively tested.

B. Padding and Unpadding Algorithms
The conjugate symmetry allows the DFT of a grid of N × N × N real values to be stored using a Fourier-space grid of only N ×N × N 2 + 1 complex values. These values are further restricted such that they contain exactly N 3 degrees of freedom. The padding algorithm, which interpolates between a Fourier-space grid of size N × N × N 2 + 1 and one of size pN × pN × pN for some positive integer p, is specified in Algorithm 1. The unpadding algorithm, which reverses the procedure, is specified in Algorithm 2.