DEFROST: A New Code for Simulating Preheating after Inflation
http://www.sfu.ca/physics/cosmology/defrost/
Abstract
At the end of inflation, dynamical instability can rapidly deposit the energy of homogeneous cold inflaton into excitations of other fields. This process, known as preheating, is rather violent, inhomogeneous and nonlinear, and has to be studied numerically. This paper presents a new code for simulating scalar field dynamics in expanding universe written for that purpose. Compared to available alternatives, it significantly improves both the speed and the accuracy of calculations, and is fully instrumented for 3D visualization. We reproduce previously published results on preheating in simple chaotic inflation models, and further investigate nonlinear dynamics of the inflaton decay. Surprisingly, we find that the fields do not want to thermalize quite the way one would think. Instead of directly reaching equilibrium, the evolution appears to be stuck in a rather simple but quite inhomogeneous state. In particular, onepoint distribution function of total energy density appears to be universal among various twofield preheating models, and is exceedingly well described by a lognormal distribution. It is tempting to attribute this state to scalar field turbulence.
pacs:
98.80.Cq, 05.10.aI Introduction
The idea of inflation is a cornerstone of the modern theory of the early universe. According to inflationary paradigm, universe at early times undergoes a period of rapid (quasiexponential) expansion, which wipes the initial state of the universe clean while seeding the primordial inhomogeneities with quantum fluctuations generated during expansion Linde:2005dd (); Linde:2005ht (). While universe is inflating, all of its energy sits in the homogeneous scalar field or condensate (known as inflaton), which is in a vacuumlike state with little entropy or particle excitations. But eventually inflation ends, and this energy has to be deposited into excitations of other matter fields, starting the thermal history of the universe with a hot big bang.
Decay of the inflaton can be very efficient if the fields experience dynamical instability at the end of inflation; such a stage became known as preheating. In most chaotic inflation models, oscillations of the inflaton field can cause instability via parametric resonance Dolgov:1989us (); Traschen:1990sw (); Kofman:1994rk (); Shtanov:1994ce (); Kofman:1995fi (). Although linear development of this instability can be understood analytically Kofman:1997yn (); Greene:1997fu (), it might be chaotic Podolsky:2002qv (), and one needs to resort to numerical simulations to investigate nonlinear dynamics that soon takes over Khlebnikov:1996mc (); Podolsky:2005bw (); Dufaux:2006ee (); Felder:2006cc (). In hybrid inflation models Linde:1993cn (), in addition to parametric resonance GarciaBellido:1997wm (), one also has a tachyonic instability associated with the symmetry breaking Felder:1998vq (), dynamics of which has been explored in Felder:2000hj (); Felder:2001kt (); GarciaBellido:2002aj ().
Nonequilibrium dynamics of preheating can lead to a multitude of interesting phenomena. Some of the topics discussed in the literature are formation of topological defects Tkachev:1998dc (); Battye:1998xe (), production of various particles (with applications to baryo and leptogenesis) GarciaBellido:1999sv (); GarciaBellido:2000dc (); GarciaBellido:2001cb (); GarciaBellido:2003wd (), possibility of primordial black hole formation Suyama:2004mz (); Suyama:2006sr (), generation of primordial magnetic fields DiazGil:2007dy (); DiazGil:2008tf (), and production of stochastic gravitational wave background Easther:2006gt (); Easther:2006vd (); GarciaBellido:2007dg (); GarciaBellido:2007af (); Dufaux:2007pt (); Caprini:2007xq (). Due to difficulties of dealing with nonlinear evolution equations, most of these studies rely on numerical simulations.
This paper describes a new code for simulating nonlinear scalar field dynamics in expanding universe developed to study preheating, called DEFROST, and the first results obtained with it. There are other codes available for this purpose, most notably LATTICEEASY by Gary Felder and Igor Tkachev Felder:2000hq (), and its parallel version CLUSTEREASY Felder:2007nz (). Through the use of more advanced algorithms and careful optimization, the new code significantly improves both accuracy and performance achievable in simulations of preheating. An important design goal has been the ease of visualization and analysis of the results, which is all too important for understanding dynamics of complex systems.
We present results on preheating in a simple twofield chaotic inflation model with massive inflaton decaying into another scalar field via quartic coupling Felder:2006cc (). Through our simulations, a new and somewhat simpler picture of the late stages of preheating emerges. After initial transient when instability develops, bubbles form and then breakup, the matter distribution soon arranges itself in a clumpy state which persists with little changes for a long time. Onepoint probability distribution function of total energy density in this state appears to be universally lognormal among various twofield preheating models. It is tempting to attribute this to relativistic turbulence Nordlund:1998wj (); Micha:2002ey (); Micha:2004bv (). We also see some evidence that the structure formed during preheating continues to grow in size on a much longer timescale. Somehow, this picture reminds one of large scale structure formation Bardeen:1985tr (); Bernardeau:1994aq (); Bond:1995yt ().
This paper is organized in the following way: In Section II, we introduce scalar field models of preheating, derive equations of motion, and discuss physical approximations we use. Section III describes the detailed implementation of numerical evolution algorithm. Initial conditions including quantum fluctuations of the fields produced during inflation are discussed in Section IV, with particular attention paid to implementation of Gaussian random field generator. We briefly recount the theory of preheating via broad parametric resonance at the end of chaotic inflation in Section V, and present our simulations in Section VI. We conclude by summarizing our main results in Section VII.
Ii Equations of Motion
As our baseline model of reheating we will take a system of scalar fields , minimally coupled to gravity and interacting through some (nonlinear) potential as described by the action
(1) 
The above action can be modified to describe more complicated geometrical quantities (like a vector field, for example) instead of a real scalar multiplet . Ultimately what we care about is only field equations of motion and their gravitational effects, not actual gauge symmetries. Variation of the action (1) with respect to the metric gives Einstein equation with a stressenergy tensor
(2) 
while variation with respect to the field gives equation of motion for the field
(3) 
Although in principle it is possible to solve the complete system of Einstein and scalar field equations numerically (see York:1971hw (); Shibata:1995we (); Baumgarte:1998te (); Pretorius:2004jg () for formulation and some approaches), in practice, it is not such an easy thing to do. Einstein equation solvers in dimensions are complex to implement, very expensive to run, and, despite marked improvement in the recent years Pretorius:2004jg (), still might have issues with numerical stability.
Fortunately for us, we do not have to solve the full Einstein equations. Although preheating is a rather violent process, and stressenergy tensor becomes very inhomogeneous due to nonlinear field dynamics, the smallness of gravitational coupling constant assures that gravitational backreaction of these inhomogeneities is rather small (at level in the simulations presented in this paper). Thus we are going to treat the scalar field evolution as if it was happening in a homogeneous flat FriedmanRobertsonWalker spacetime
(4) 
and calculate the inhomogeneous gravitational field due to matter distribution with stressenergy tensor (2) using linear perturbation theory Bardeen:1980kt (). We will ignore backreaction of metric perturbations on the scalar field evolution.
With these simplifying assumptions, the problem becomes much more tractable: we just need to solve a system of coupled scalar field equations of motion (3), which in spacetime with metric (4) become
(5) 
Here and later, dot denotes the derivative with respect to time , and spatial differential operators (gradient and Laplacian ) are taken with respect to comoving coordinates and threedimensional flat metric. The expansion rate and acceleration are determined by averaged Einstein equations
(6) 
where energy density and isotropic pressure are components of the stress energy tensor (2) given by
(7)  
(8) 
and averages are taken over the whole simulation volume.
To solve equation (5) numerically, one needs to know the Hubble parameter value . Rather than attempting to resolve the constraint equation (6) at every time step (which would result in an implicit evolution scheme), it is faster and more convenient to use the evolution equation
(9) 
to evaluate its value in the future. This is what is done in LATTICEEASY. However, we have one more trick up our sleeves: a disproportionately huge gain in numerical accuracy can be realized by evolving the Hubble length instead of the Hubble parameter by using
(10) 
For constant equation of state , Hubble parameter evolves as , while Hubble length evolves as . The latter variable has vanishing second (and higher) derivatives and correspondingly smaller truncation error when discretized to second order. As an added bonus, the spatial gradients (which are the single most expensive thing to calculate) cancel out when taken in combination, and do not enter evolution equation in the form (10).
Once the field equations of motion (5) are solved, we can evaluate all the components of stressenergy tensor (2), in particular energy density (7) and pressure (8), as well as calculate linear metric perturbations they create in a homogeneous spacetime (4). In this paper, we will focus on scalar perturbations, behaviour of which during preheating have not been widely studied yet. In the longitudinal gauge, perturbed metric can be written as
(11) 
The two scalar gravitational potentials and are equal in the absence of anisotropic stress. This is in general not the case for the scalar field models, and can be expected to hold only approximately and in the average sense. Both gravitational potentials and are nondynamical, being solutions of the Poissonlike equations. In particular, the equation for is
(12) 
Strictly speaking, what should stand in the right hand side is not the density , but the gaugeinvariant density variable , which has some extra terms in it Bardeen:1980kt (). At this stage of the code development, I will not make that distinction and simply ignore the missing terms (which usually works well on subhorizon scales), along with anisotropic stress contribution to potential . Full support of gaugeinvariant perturbations including tensor and vector modes is planned for the future code release.
Iii PDE Solver Implementation
Scalar field equations of motion (5) are coupled nonlinear partial differential equations, and have to be solved numerically. Fortunately, all the nonlinearity comes from the potential term only; the differential operator itself is simple, homogeneous and hyperbolic, which makes the numerical solution quite straightforward. For the solver, we adopt a secondorder accurate finite difference scheme based on leapfrog algorithm.
The scalar field values are discretized on threedimensional cubic grid in comoving coordinates with uniform spacing and periodic boundary conditions. Since evolution equations (5) are second order, values of the fields on two consecutive time slices are required to advance to the next one. We will denote the previous, current, and next time slices by indices dn, hr, and up correspondingly. The time derivatives of a quantity are discretized to second order as
(13) 
coefficient  cost  stability  

degeneracy  
standard  
isotropic A  
isotropic B  
isotropic C 
Discretization of spatial differential operators
(14) 
allows more freedom. Direct generalization of the second derivative discretization in equation (13) to three spatial dimensions leads to the oftenused secondorder accurate expression for Laplacian using six nearest neighbours of a point. However, this is not the only (or the best) choice. Truncation error for this scheme depends on direction, introducing anisotropic artifacts in the field evolution at short length scales. One can mitigate this unpleasant fact by increasing resolution, but smarter discretization scheme is a better solution. By using all 26 neighbours in a cube around a point, one can derive a family of discretizations of Laplacian operator which is secondorder accurate and fourthorder isotropic Patra:2005 (). For discretization to be isotropic, the coefficients in a linear combination approximating the Laplacian operator should only depend on the distance from the central point
(15) 
as illustrated in Figure 1. The values of coefficients for possible discretizations of the Laplacian operator are summarized in Table 1, along with their computational cost and stability properties. Isotropic discretizations and offer reduced computational cost, while isotropic discretization has the best accuracy and stability. As multiplications are cheap and additions are essentially free on modern CPUs, there is no reason not to use the best discretization scheme available. Thus, DEFROST is configured to use the isotropic discretization by default. That can be changed by uncommenting the other coefficient definitions in the source code, although profiling shows little gain in doing so for large grids (for which performance of DEFROST solver is apparently memorybandwidth dominated).
To calculate the energy density (7), we need to discretize the square of spatial gradients as well. It is very important for discretized energy to be conserved by discretized equations of motion. Otherwise, it will leak off the grid in the course of a long simulation, affecting overall accuracy or even giving incorrect results (for equation of state, for example). Simply squaring the first spatial derivative discretized like in equation (13), although secondorder accurate, is not conservative. Using discretized action approach Pen:1993nx (), one can show that the conservative secondorder accurate and fourthorder isotropic discretization of the gradientsquare operator is
(16) 
where coefficients are the same as in Laplacian (15). Evaluating this expression requires 30 multiplies per point for discretization , which is significantly more expensive than computing the Laplacian. However, this does not place much of a burden on the total runtime. As we mentioned before, we eliminated gradientsquare terms from evolution equations, so they only need to be calculated for output, which happens much less often.
Putting everything together, we end up with the following evolution scheme: the discretized field values are advanced to the next time slice using
(17) 
where is the discretized Laplace operator (15) and . The time step has to be small enough both to satisfy Courant stability condition ( listed in Table 1) and to resolve the period of the fastest oscillating field. All the coefficients in expression (17) except for the effective mass of the th field
(18) 
are constant on the grid and the same for all fields, and thus can be precomputed outside the evolution loop. Inside the same loop, kinetic and potential energy of the fields are accumulated to evaluate
(19) 
which is used to advance the Hubble length
(20) 
To avoid weak numerical instability associated with evenodd slice decoupling, we use instead of in the above equation, which then can be solved for either as an exact quadratic or iteratively (we use the latter in the code).
Equations (17) and (20) provide a complete recipe how to advance the field variables to the next time step. Once in a while (at user’s request) we would also like to calculate, analyze, and output for visualization some auxiliary variables like energy density and gravitational potential . Energy density (7) and pressure (8) are easy to calculate once we know the field values and their gradients (16). Finding the gravitational potential is a little less trivial, as we need to solve the Laplace equation (12) with periodic boundary conditions. The fastest way to do it is to use a fast Fourier transform (FFT). Applying FFT to the discretized Laplacian operator (15), we end up with an algebraic equation for in Fourier domain
(21) 
where is the integer Fourier mode wavenumber and polynomial follows from discretization (15)
(22) 
with . Here are once again the coefficients of discretization (15) with values listed in Table 1.
DEFROST implements the evolution scheme (17) in Fortran 90, fully taking advantage of capabilities of modern hardware and compilers by using both automatic vectorization (over field variables) and automatic parallelization (of the evolution loop). During code development, profiling showed that physical memory layout has an unexpectedly large impact on performance, which became apparent as the solver loop got optimized. Let us briefly discuss the storage model which was adopted after some investigation. The fields are sampled and stored in a large multidimensional array smp(N,0:p,0:p,0:p,3). The first index (minor in Fortran index ordering) enumerates the fields. The next three indices enumerate the three dimensions of the spatial grid, padded to elements for reasons discussed below. The last index (major in Fortran index ordering) enumerates the three time slices used in the evolution scheme. Rather than copying large amounts of data around, the allocation of indices for dn, hr, and up slices cycles every time step.
The layout of a single time slice is shown in Figure 2 (with one spatial dimension suppressed for clarity). The simulation volume is surrounded by a single cell wide boundary layer, introduced to implement the boundary conditions without conditional logic inside the evolution loop. It also allows for an easy transition to parallel cluster implementation using MPI, as required buffers are already in place. Somewhat counterintuitively, padding the grid by a few extra empty cells can significantly improve the runtime for large grids. The root cause for this phenomenon probably lies in some interaction between memory access pattern and hardware memory cache algorithm. To get the best FFT performance, the grid size is usually taken to be a power of two. While evaluating (17) on an unpadded array, memory would be accessed with a power of two stride, which could conceivably interfere with caching and prefetching done by the memory subsystem. Whatever the cause, padding the array so that its size is prime (or a product of a few large primes) can improve the runtime, so the user is advised to experiment.
Finally, a few words should be said about statistical estimators used to analyze the simulation data. The power spectral density (PSD) estimator is a fairly conventional one implemented using FFT. It employs antialiasing with fourthorder polynomial kernel when folding the spectrum into wavenumber bins to reduce sampling noise. The implementation of probability density function (PDF) estimator is less conventional, and does not use histogram binning at all. Instead, PDF is derived from cumulative density function (CDF) which is obtained by partially sorting the data cube into quantile brackets. Although more expensive and harder to parallelize, this approach is more robust and offers uniform sampling noise across the distribution.
Iv Initial Conditions
At the end of the inflation, most of the energy is still stored in the inflaton, and all the fields are homogeneous except for small quantum fluctuations. But the presence of these quantum fluctuations in the fields is essential to trigger the dynamical instability leading to preheating.
Initial conditions in the homogeneous field components depend on the model of inflation and are treated as an external input by DEFROST. They are straightforward to obtain by following the expanding FriedmanRobertsonWalker solution during inflation either analytically (if possible), or using any of the available numerical ordinary differential equation integrators. As inflationary trajectory is an attractor Belinsky:1985zd (), it is easy to find and no particular care is needed on where to start tracing it. The full threedimensional simulation of preheating should take over at some time near the end of the inflation, but as there are other factors at play (such as limited spatial dynamic range available to simulation), the exact moment is best decided on case by case basis.
To make a concrete example, Figure 3 shows expansion history of a typical chaotic inflation model with massive inflaton (discussed in more detail in the next section). As the Universe is inflating, its horizon size stays relatively constant, while the physical wavelength of comoving modes grows with the scale factor . The modes which were originally inside the horizon expand and leave the horizon during inflation. Eventually, inflation ends and the horizon size starts growing faster than the scale factor (for instance, during matter domination), at which point the modes begin to reenter the horizon. The moment in time when comoving modes stop leaving the horizon and begin reentry can be taken as the end of the inflation. This happens when
(23) 
or in terms of a slow roll parameter
(24) 
For simulations presented in this paper, we start exactly at the moment when inflation ends (24), and select the comoving box size of to cover all the length scales of interest, as illustrated in Figure 3.
As we already mentioned, it is crucial to include quantum field fluctuations in the mostlyhomogeneous initial conditions for the preheating instability to develop. The spectra of quantum field fluctuations are determined by effective masses of the fields involved. Let us briefly recount the standard derivation Linde:2005ht () to establish notation. Canonically normalized massive field with Lagrangian
(25) 
is quantized in flat Minkowski spacetime by promoting the field value and field momentum
(26) 
to quantum operators and obeying canonical commutation relation
(27) 
In flat spacetime, the field operator can be represented using Fourier mode decomposition as
(28) 
where the mode creationannihilation operators obey
(29) 
which directly follows from canonical commutation relation (27). The mode frequency of the massive field is related to its wavevector by a simple dispersion relation
(30) 
Using mode decomposition (28) and commutation relation (29), it is easy to show that the twopoint correlation functions of the field value and momentum are
(31)  
(32) 
The spectrum of the field fluctuations is simply .
One can repeat this procedure in the background of expanding homogeneous universe. In general, the time dependence of the modes will be different, but it turns out that quantum fluctuations of massive fields in de Sitter spacetime can be approximated quite well with above expressions if one simply replaces the field mass with an effective mass
(33) 
We will use this approximation in DEFROST. In addition, we will follow the established practice of treating quantum operators as Gaussian random variables. This is an assumption, but not totally unjustified one. Quantum modes essentially behave classically after they leave the horizon Polarski:1995jg (). Even the subhorizon modes (which are initially quantum) get large occupation numbers once the preheating instability kicks in, and can be treated classically Khlebnikov:1996mc (). Thus, we will initialize the field fluctuations as a Gaussian random field
(34) 
where the complex random operators and obey
(35) 
A lot of effort has gone into making the realization of random field initial conditions in DEFROST as statistically accurate as possible. The straightforward and often used way to generate a Gaussian random field on a discrete grid is to directly discretize equation (34) in Fourier space, assigning th mode a random Gaussian number with amplitude . Although simple, this procedure is spoiled by the finite grid size effects, and does not reproduce correct twopoint correlation functions in the real space Pen:1997up (). One ends up with a substantial lack of power on the scales comparable with the box size, which is not surprising if one considers that only a few longwavelength modes “fit” into the box, and naive discretization ignores all the power in the infrared part of the spectrum which should have been properly aliased into those few low modes.
A wealth of the literature is dedicated to the subject of generating Gaussian random fields of a given spectrum, particularly in the context of the body simulations of the large scale structure Fournier:1982 (); Efstathiou:1985re (); Hoffman:1991 (); Hoffman:1992 (); Salmon:1996 (); Pen:1997up (); Bertschinger:2001ng (); Plaszczynski:2005yp (). We draw on that experience, and in DEFROST for random field generator we adopt a method described in Salmon:1996 (); Pen:1997up (). Gaussian random field (34) with spectrum is realized by convolving white noise with a spherically symmetric kernel function
The kernel function can be evaluated analytically in terms of Bessel functions
(37) 
and has a power law ultraviolet divergence in the limit . For the purposes of discretization on a finite grid, we have to regularize this divergence, which we do by introducing a Gaussian cutoff at some scale below the Nyquist frequency
(38) 
The regularized kernel does not have a nice expression in terms of elementary functions, but is easy to evaluate numerically, of course. We use a onedimensional discrete sine transform (DST) on a substantially larger grid to calculate it (simply because DST is already provided by FFTW libraries we use), but other methods like quadrature integrators could be used as well. Once the sphericallysymmetric kernel is evaluated, it is sampled on the the threedimensional grid in real space, and the random field is initialized as a convolution
(39) 
The convolution is implemented using discrete FFT as
(40) 
where is the grid size, is the spacing of discrete wavemodes, and is a complex Gaussian random number generated using BoxMuller transformation
(41) 
from two real random numbers and uniformly distributed on a unit interval. The implementation of initial velocity generator is entirely analogous (and handled by the same procedure), and is not worth repeating here.
Finally, I have to point out that as of current writing, there is a bug in the random number generator implementation in LATTICEEASY. The formula (41) is modified there to produce two complex numbers using only three uniformly distributed real numbers. This results in correlated random numbers with quite nonGaussian distribution. Fortunately, the results reported so far do not seem to be too much affected by this problem, but as nongaussianity studies are becoming more prominent in the modern cosmology literature, some care must be taken in proper implementation of random number generators.
V Chaotic Inflation and Broad Parametric Resonance
So far, the discussion was rather general, as DEFROST is designed to be easily adoptable to study arbitrary models of preheating. This Section describes preheating model we selected for first simulations with DEFROST: chaotic inflation ending via broad parametric resonance Kofman:1994rk (). The development of linear instability in this model can be largely understood analytically Kofman:1997yn (), and its nonlinear dynamics has been widely studied numerically as well Podolsky:2005bw (); Felder:2006cc (). For all its simplicity, this model has very rich dynamics, and still holds surprises. Our very first simulations uncover new aspects of the evolution dynamics in this model, which are reported in the next Section.
In a minimal form, the model consists of two scalar fields: the massive inflaton and the massless decay product interacting via potential
(42) 
During inflation, the value of the inflaton is large, the field is overdamped by the large Hubble friction, and slowly rolls down its potential. As it reaches the value of around one in Planck units, the damping dips below critical, and homogeneous inflaton starts oscillating with decreasing amplitude
(43) 
Decay field is coupled to inflaton, and feels its oscillations through modulation of the effective mass; equation of motion for the Fourier mode with wavector is
(44) 
If the coupling is large enough, periodic modulation of the field mass leads to strong instability via parametric resonance. This can be understood analytically by applying general theory of differential equations with periodic coefficients Kofman:1997yn (). If one ignores the expansion and the Hubble drag term in equation (44), evolution for the Fourier mode is given by Mathieu equation
(45) 
where we have introduced dimensionless parameters
(46) 
and time variable to bring the equation into canonical form. According to Flouqet’s theorem, a general solution of Mathieu equation (45) is of the form , where is a periodic function with period . Floquet exponent depends on parameters and , and there is an elegant way to calculate its value Bateman:1955 (), which we (somewhat reluctantly) will omit here, and just quote the final result. Figure 4 shows the dependence of on parameters and as a density plot. For certain parameter values Floquet exponent has positive real part, leading to exponential instability of the solution; these unstable bands are marked on Figure 4. The value of for our problem (46) is restricted to lie above the red line in Figure 4, which corresponds to the homogeneous mode with . For sufficiently large coupling , large portion of available phase space volume is unstable, leading to fast development of instability. This regime is known as broad parametric resonance. The instability grows on a timescale comparable to (as Floquet exponent values are around ), and manifests itself after a few dozen of oscillations of the inflaton, which is very fast in cosmological terms. In the next Section, we describe the nonlinear field evolution after this instability develops.
Vi Numerical Results
Before we present our simulation results, a few words should be said about the units used throughout the code. As it is clear from the action (1), we prefer to work with dimensionless scalar fields, rather than canonically normalized ones. In this convention, one can simply think of the values of the scalar fields as measured in units of Planck mass . Note that we use reduced Planck mass rather than . The only other scale in the model (42) is the inflaton mass . The coupling constant (which has dimension of mass with our scalar field normalization), the frequency and wavevectors of the unstable modes, and everything else can be referred to it. It is convenient to scale all the quantities (except for field values) by the appropriate power of , making them dimensionless and suitable for numerical analysis. Thus we set in the code.
We simulate the preheating model (42) with inflaton mass , the value of coupling constant in a comoving box of size (the values chosen correspond to the ones used in Ref. Felder:2006cc ()). The grid size is taken to be , while the time step has to be reduced substantially below the Courant limit to resolve oscillations of the field (which is initially hundred times heavier than inflaton ). The simulation is started at the end of the inflation (24), when the value of the inflaton is and the Hubble constant is . The simulation box is initially about five Hubble lengths across. We let the code run until , which corresponds to time steps.
To get things going, we first reproduce the previously reported results Podolsky:2005bw () on expansion history of the preheating model (42). The left panel of Figure 5 shows evolution of horizon size during expansion, with growth corresponding to matterdominated expansion scaled out. The expansion history has a sharp break as instability develops, and energy gets deposited into relativistic inhomogeneous modes from a homogeneous oscillating inflaton (which behaves as a pressureless dust). This transition can be seen in terms of an effective equation of state parameter , the value of which is plotted in the right panel of Figure 5. It undergoes large amplitude oscillations (shown by thin pale line in Figure 5), but when averaged over a few periods (with a Kaiser window function), the underlying behavior is uncovered. The averaged equation of state (shown by thick red line) switches over from dustlike equation of state to a value slightly less than a quarter Podolsky:2005bw (), corresponding to a fairly relativistic fluid. (If evolved further, the residual homogeneous component in the inflaton will eventually come to dominate the evolution again, slowly lowering equation of state toward in the process.)
While we recover the results obtained by simulations with LATTICEEASY, the accuracy of the integrator used in DEFROST is significantly higher. The performance of integration scheme for expansion factor (20) is illustrated in Figure 6, which shows residual curvature
(47) 
which should be zero in the flat model we are evolving. As you can see, the constraint equation is satisfied to level (which is exceptionally good for a second order scheme with timestep), and the error does not accumulate with time. In fact, this error is mostly due to the fact that we neglected second order corrections to density from initial field fluctuations.
Having made sure that our code reproduces previously reported results, and after numerous checks of code integrity and accuracy, we move on to investigation of the field dynamics during preheating.
Evolution of field distributions and spectra for inflaton , decay field , total density and gravitational potential are presented in Figure 7. Left panel shows evolution of the median value (thick red line), along with 68% and 95% percentile brackets around it (which would correspond to and contours for a Gaussian distribution) shown by shaded outlines. The contour with the lightest shading spans the extremal values inside the simulation box, which serves to illustrate the extent of the tails of the distribution, although the exact percentile it corresponds to depends on the spatial resolution of the simulation. Dilution due to expansion has been scaled out to highlight the relative change of the distributions as evolution proceeds.
The evolution of the distribution of the decay field values (second row of Figure 7) clearly shows the onset of instability a little after , rapid spreading of the distribution due to exponential amplification of seed inhomogeneities, and selflimiting of the growth by nonlinear interactions when the scaled value of the decay field becomes of order unity. As the decay field perturbation grows and becomes nonlinear, it is drawing the energy from the zero mode of the inflaton (top row of Figure 7), reducing the amplitude of its oscillations, and eventually forcing the inflaton to become strongly inhomogeneous as well due to nonlinear backreaction. We should note here that although the amplitude of the coherent inflaton oscillation decays, it does not go away altogether, and the leftover homogeneous mode will eventually come to dominate the universe expansion Dufaux:2006ee (), as its equation of state is effectively that of pressureless dust, and it dilutes slower than inhomogeneous components, which have equation of state closer to relativistic one.
Of particular interest to us is the distribution of the total energy density, shown in the third row of Figure 7. It is clearly very inhomogeneous, with peak densities easily exceeding ten times the average. After a brief transient, it quickly settles to a nearly stationary distribution, which appears to be highly nonGaussian (and is in fact plotted on a logarithmic scale). We will come back to this point after we inspect the spatial picture of the energy density distribution inside the simulation box.
The bottom row of Figure 7 shows the evolution of the gravitational potential. Despite total energy density being highly inhomogeneous and having huge overdensities, the gravitational potential it produces is rather small (at level), and is further diluted away by the expansion with nearrelativistic equation of state. The maximal potential well depth inside the simulation box is , which is far too small to form any primordial black holes. This result is in line with observations of Suyama:2004mz (); Suyama:2006sr (), although now we have a more transparent diagnostic of black hole formation as we calculate the gravitational potential directly.
The right panel of Figure 7 shows evolution of the field spectra as the density plot in terms of time and wavenumber , with dilution of the fields due to expansion and overall powerlaw dependence on scaled out. The spectra of total energy density and gravitational potential show a very clear simultaneous peak soon after instability develops, which is sharply localized in both time and scale. Subsequent evolution of the two differs, however. The peak power in energy density is evolving toward higher and smaller scales, while the peak power in potential is evolving to lower , so the structure of potential wells is growing in spatial extent. The understanding of what’s going on will become more clear when we look at evolution in real space, which is what we are going to discuss next.
Figure 8 shows threedimensional volume renders of the contents of the simulation box. As the fields and oscillate rapidly in a standing wave pattern, quickly losing coherent phasing, the view of their values is messy and not very enlightening, and we will omit it here. Much more interesting is the picture of what is happening to the total energy density, which is an adiabatic invariant for the oscillating fields. The top row of Figure 8 shows density distribution inside the simulation box soon after onset of instability (at , left) and at the end of the simulation (at , right). Density field is shaded using logarithmic color map with linear transparency ramp applied, so that only the peaks of the density distribution are visible.
Immediately after the onset of instability, the density distribution in Figure 8 (top left) looks like smoke is filling the box. What you are seeing is actually the overdense bubble walls forming a threedimensional foamlike structure that fills the box. Its origin is easy to understand if one thinks about how seed inhomogeneities are amplified by instability. Broad parametric resonance amplifies wavemodes in a certain band, effectively serving as a lowpass filter (with a kernel that can be approximated analytically) and sets the characteristic size of the structure which grows out of the seed inhomogeneities. Original fluctuations are a Gaussian random field, which already has the structure of peaks, ridges, and valleys imprinted into it. The skeleton of this structure is essentially preserved unchanged as the growth of inhomogeneities due to instability increases density contrast. Once the density contrast becomes of order unity, nonlinear evolution takes over. This will happen to underdense regions first, with repulsive interaction term helping to evacuate the bubble interiors, and pushing the matter density into the bubble walls, thus forming the structure you see in Figure 8 (top left).
The evolution does not stop at forming bubbles, however. The repulsive interaction soon breaks the extended bubble walls into smaller localized blobs, moving more or less freely inside the simulation box (the animation of this process is available online). The final state is depicted in Figure 8 (top right), and persists with little change for a long, long time. This state seems quite distinct from thermal equilibrium, yet it is still longlived and statistically simple in a certain way, which is quite surprising. Even more surprisingly, the distribution of values of total energy density in units of quickly becomes statistically stationary and after a brief transient tends to a distribution with a probability density function shown in the lower right of Figure 8. It can be fitted with exceedingly high accuracy by a lognormal distribution
(48) 
with one free parameter ( or ), as the mean is unitnormalized by virtue of us scaling out the expansion. The corresponding median is . With statistical errors of PDF estimator being what they are, the apparent lognormality of a density distribution is undoubtedly not a mere coincidence, but must have a explanation rooted in scalar field dynamics. Moreover, further simulations of preheating models with different couplings and inflaton potentials seem to suggest that lognormal distribution of density is a universal feature of twofield preheating models. This observation presents a very interesting theoretical puzzle, which will be explored in detail elsewhere Frolov:preview ().
Although onepoint distribution of total energy density quickly becomes stationary as noted above, other quantities continue evolving on much longer time scales. The blobs continue to fragment, and their characteristic size slowly decreases with time. While this is obvious from visualizations, it can be further quantified by introducing the (physical) correlation length of the total energy density configuration
(49) 
The evolution of the comoving correlation length is plotted in the left panel of Figure 9. Initially it is very large (), as the density field is nearly homogeneous. As instability develops and structure forms, it abruptly drops to about , and then continues to decrease, but much more slowly. Although the graph clearly shows evolutionary trend in density correlation length , actual numbers should be taken with a grain of salt, as the density field does eventually become fragmented on a scale close to spatial grid resolution (which for grid we use would be reached around ).
While it is known that thermalization after preheating might take a long time Kofman:1995fi (), and there might be an intermediate scaling regime in the evolution of the fields Micha:2002ey (); Micha:2004bv (), the view of the process from the real space presented here is strikingly simple. One would think the field distributions will thermalize eventually, presumably forming a homogeneous fluidlike state with Maxwellian distribution of particle velocities. Thermalization does not happen in our simulations, which lack necessary quantum effects. Exactly how and when it does happen is an extremely interesting question, and the one which requires further study.
Finally, let us discuss gravitational potential , which is sourced by the evolving energy density distribution , and is shown in the middle row of Figure 8. To make the structure more visible, we have opted for a density plot on a threeslice through the simulation box rather than a volume rendering. The color map shows positive potential values (corresponding to underdense regions) as shades of red, and negative potential values (corresponding to overdense regions) as shades of blue, blending into white for zero potential value.
Immediately after the onset of instability, the gravitational potential in Figure 8 (middle left) clearly traces the foamlike structure of matter distribution. The isolated potential peaks (red) in the interior of the bubbles are separated by extended potential valleys (blue) created by overdense bubble walls. The gravitational potential configuration is asymmetric between positive and negative values, and is clearly nonGaussian. Subsequent evolution of the gravitational potential is rather interesting. As bubble walls break into smaller and smaller blobs, the structure of the gravitational potential does not follow suit. Instead, it begins to grow in spatial extent (the animation is available online). By the end of the simulation in Figure 8 (middle right), the size of the structure in the gravitational potential spans almost the whole box. The growth of structure can be quantified by introducing correlation length for gravitational potential the way we did for energy density
(50) 
Evolution of the comoving correlation length is plotted in the right panel of Figure 9, with pale line tracing its instantaneous value, and thick red line giving a running average (with Kaiser window) over few oscillations. While comoving correlation length grows overall, evolution shown in Figure 9 (right) still represents an initial transient, and the longtime asymptotic behaviour is not reached in the simulation reported here. Further investigation shows that correlation length continues to grow faster than comoving size, but not quite as fast as the horizon size. Eventually one might even have to worry about it outgrowing the finite simulation box size, but that would take a very long time, and is not reached by our simulations.
Vii Conclusions
This paper presents a new numerical code I developed for simulating preheating of the Universe after the end of the inflation, which I call DEFROST. It is small (about 600 lines of Fortran code), fast, easy to modify, and is fully instrumented for 3D visualizations (using LLNL’s VisIt, for example). The source code is available for download online at http://www.sfu.ca/physics/cosmology/defrost/ and is distributed under the terms of GNU Public License. While the main design goal of DEFROST has been the accuracy of the simulations, performance of the solver has also been significantly improved compared to LATTICEEASY Felder:2000hq (), which is the most mature and widely used reheating code publicly available today.
As a result of all the optimizations (and a bit of black magic), DEFROST outperforms LATTICEEASY by about a factor of four in raw PDE solver speed (for 2 fields on a grid in double precision on a dual Xeon 5160 machine) while using more accurate (and more expensive) discretization. If one takes into account the time spent on analysis of the results, the difference is even larger, as FFTW libraries used by DEFROST are vastly faster than FFT routines shipped with LATTICEEASY (especially on multiprocessor machines). The speedup offered by DEFROST is so significant that the studies done a few years ago on a big parallel cluster Podolsky:2005bw (); Dufaux:2006ee () using MPI version of LATTICEEASY Felder:2007nz () can now be carried out on a single fast workstation. The planned MPI version of DEFROST should be able to push the accessible simulation size over the barrier, provided the code scales well.
The code was tested on a number of chaotic inflation models which end via parametric resonance. In this paper, we report the simulations of the simplest twofield preheating model with massive inflaton and quartic coupling to decay field (42). We reproduce the previously published numerical results for this model Podolsky:2005bw (); Felder:2006cc () (and the ones for trilinear coupling Dufaux:2006ee (), which we will not discuss here, although our simulation data and results for that model are available online as well). We further investigate the dynamics of the scalar field evolution in these preheating models, taking advantage of advanced visualization and analysis capabilities DEFROST offers. In particular, we study the behaviour of energy density distribution and scalar gravitational potential during preheating, something which has not been looked at closely before. Our main science results are summarized by two observations, both novel and quite surprising.
First, the evolving scalar fields quickly end up in a simple state, which, although highly inhomogeneous, appears to have a certain universality to it. In this state, the onepoint distribution function of total energy density is nearly stationary (apart for the overall dilution due to expansion), and is described by a lognormal distribution for all twofield parametric resonance preheating models we tried so far, namely the ones described by interaction potentials

,

,

,

.
This is true even if distributions of field values or other correlators might still be evolving, and appears to be a very general statement about random scalar fields one encounters in preheating. It is tempting to attribute this state to scalar field turbulence Micha:2002ey (); Micha:2004bv (), especially since lognormal density distributions are known to occur in supersonic isothermal turbulence in hydrodynamics Nordlund:1998wj (). We do not see obvious signs of thermalization, even if the simulations are run for a time much longer than the dynamical timescale of the problem (the longest done so far for massive inflaton is corresponding to five folds since the end of inflation; this is limited mainly by my patience rather than the code stability).
Second, less general but still amusing, is the observation that the smallscale structure in the gravitational potential can grow faster than comoving box expands. It is not quite clear whether the reason for it happening is kinematic or dynamical in nature. As we neglected gravitational interactions in our simulations, the only thing that can cause the structure to grow is the interaction between scalar fields themselves. In our preheating model it is repulsive; yet the structure still grows! Although one might suspect that any inhomogeneity in gravitational potential on subhorizon scales would probably get washed away by subsequent evolution (and is too small to form primordial black holes), this effect still might have some interesting cosmological consequences.
All in all, we find that the picture of preheating dynamics is simpler in real space than what it looks like in particle representation. The final stage of preheating, with growing structure and lognormal density distribution, eerily reminds one of large scale structure formation in later cosmology (although of course it occurs on vastly smaller scales and is driven by completely different physics). Perhaps the analytical methods developed for the latter Bardeen:1985tr (); Bernardeau:1994aq (); Bond:1995yt () could be fruitfully applied to preheating as well. This is what we intend to explore next.
Acknowledgments
This work was supported by the Natural Sciences and Engineering Research Council of Canada under Discovery Grants program. Numerical computations were done in part on Sunnyvale cluster at Canadian Institute for Theoretical Astrophysics. The author is grateful for hospitality during his stays at CITA, where some portion of this work has been carried out, and would like to thank Lev Kofman, Andrei Linde, Gary Felder, and Mustafa Amin for helpful discussions of preheating, and Chris Matzner for bringing the subject of isothermal supersonic turbulence to my attention.
References
 (1) A. Linde, “Inflation and string cosmology,” eConf C040802, L024 (2004) [J. Phys. Conf. Ser. 24, 151 (2005 PTPSA,163,295322.2006)] [arXiv:hepth/0503195].
 (2) A. D. Linde, “Particle physics and inflationary cosmology,” Contemp. Concepts Phys. 5, 1362 (1990) [arXiv:hepth/0503203].
 (3) A. D. Dolgov and D. P. Kirilova, “Production of particles by a variable scalar field,” Sov. J. Nucl. Phys. 51, 172 (1990) [Yad. Fiz. 51, 273 (1990)].
 (4) J. H. Traschen and R. H. Brandenberger, “Particle production during outofequilibrium phase transitions,” Phys. Rev. D 42, 2491 (1990).
 (5) L. Kofman, A. D. Linde and A. A. Starobinsky, “Reheating after inflation,” Phys. Rev. Lett. 73, 3195 (1994) [arXiv:hepth/9405187].
 (6) Y. Shtanov, J. H. Traschen and R. H. Brandenberger, “Universe reheating after inflation,” Phys. Rev. D 51, 5438 (1995) [arXiv:hepph/9407247].
 (7) L. Kofman, A. D. Linde and A. A. Starobinsky, “Nonthermal phase transitions after inflation,” Phys. Rev. Lett. 76, 1011 (1996) [arXiv:hepth/9510119].
 (8) L. Kofman, A. D. Linde and A. A. Starobinsky, “Towards the theory of reheating after inflation,” Phys. Rev. D 56, 3258 (1997) [arXiv:hepph/9704452].
 (9) P. B. Greene, L. Kofman, A. D. Linde and A. A. Starobinsky, “Structure of resonance in preheating after inflation,” Phys. Rev. D 56, 6175 (1997) [arXiv:hepph/9705347].
 (10) D. I. Podolsky and A. A. Starobinsky, “Chaotic reheating,” Grav. Cosmol. Suppl. 8N1, 13 (2002) [arXiv:astroph/0204327].
 (11) S. Y. Khlebnikov and I. I. Tkachev, “Classical decay of inflaton,” Phys. Rev. Lett. 77, 219 (1996) [arXiv:hepph/9603378].
 (12) D. I. Podolsky, G. N. Felder, L. Kofman and M. Peloso, “Equation of state and beginning of thermalization after preheating,” Phys. Rev. D 73, 023501 (2006) [arXiv:hepph/0507096].
 (13) J. F. Dufaux, G. N. Felder, L. Kofman, M. Peloso and D. Podolsky, “Preheating with trilinear interactions: Tachyonic resonance,” JCAP 0607, 006 (2006) [arXiv:hepph/0602144].
 (14) G. N. Felder and L. Kofman, “Nonlinear inflaton fragmentation after preheating,” Phys. Rev. D 75, 043518 (2007) [arXiv:hepph/0606256].
 (15) A. D. Linde, “Hybrid inflation,” Phys. Rev. D 49, 748 (1994) [arXiv:astroph/9307002].
 (16) J. GarciaBellido and A. D. Linde, “Preheating in hybrid inflation,” Phys. Rev. D 57, 6075 (1998) [arXiv:hepph/9711360].
 (17) G. N. Felder, L. Kofman and A. D. Linde, “Instant preheating,” Phys. Rev. D 59, 123523 (1999) [arXiv:hepph/9812289].
 (18) G. N. Felder, J. GarciaBellido, P. B. Greene, L. Kofman, A. D. Linde and I. Tkachev, “Dynamics of symmetry breaking and tachyonic preheating,” Phys. Rev. Lett. 87, 011601 (2001) [arXiv:hepph/0012142].
 (19) G. N. Felder, L. Kofman and A. D. Linde, “Tachyonic instability and dynamics of spontaneous symmetry breaking,” Phys. Rev. D 64, 123517 (2001) [arXiv:hepth/0106179].
 (20) J. GarciaBellido, M. Garcia Perez and A. GonzalezArroyo, “Symmetry breaking and false vacuum decay after hybrid inflation,” Phys. Rev. D 67, 103501 (2003) [arXiv:hepph/0208228].
 (21) I. Tkachev, S. Khlebnikov, L. Kofman and A. D. Linde, “Cosmic strings from preheating,” Phys. Lett. B 440, 262 (1998) [arXiv:hepph/9805209].
 (22) R. A. Battye and J. Weller, “Cosmic structure formation in hybrid inflation models,” Phys. Rev. D 61, 043501 (2000) [arXiv:astroph/9810203].
 (23) J. GarciaBellido, D. Y. Grigoriev, A. Kusenko and M. E. Shaposhnikov, “Nonequilibrium electroweak baryogenesis from preheating after inflation,” Phys. Rev. D 60, 123504 (1999) [arXiv:hepph/9902449].
 (24) J. GarciaBellido, S. Mollerach and E. Roulet, “Fermion production during preheating after hybrid inflation,” JHEP 0002, 034 (2000) [arXiv:hepph/0002076].
 (25) J. GarciaBellido and E. Ruiz Morales, “Particle production from symmetry breaking after inflation,” Phys. Lett. B 536, 193 (2002) [arXiv:hepph/0109230].
 (26) J. GarciaBellido, M. GarciaPerez and A. GonzalezArroyo, “ChernSimons production during preheating in hybrid inflation models,” Phys. Rev. D 69, 023504 (2004) [arXiv:hepph/0304285].
 (27) T. Suyama, T. Tanaka, B. Bassett and H. Kudoh, “Are black holes overproduced during preheating?,” Phys. Rev. D 71, 063507 (2005) [arXiv:hepph/0410247].
 (28) T. Suyama, T. Tanaka, B. Bassett and H. Kudoh, “Black hole production in tachyonic preheating,” JCAP 0604, 001 (2006) [arXiv:hepph/0601108].
 (29) A. DiazGil, J. GarciaBellido, M. Garcia Perez and A. GonzalezArroyo, “Magnetic field production during preheating at the electroweak scale,” Phys. Rev. Lett. 100, 241301 (2008) [arXiv:0712.4263 [hepph]].
 (30) A. DiazGil, J. GarciaBellido, M. G. Perez and A. GonzalezArroyo, “Primordial magnetic fields from preheating at the electroweak scale,” JHEP 0807, 043 (2008) [arXiv:0805.4159 [hepph]].
 (31) R. Easther and E. A. Lim, “Stochastic gravitational wave production after inflation,” JCAP 0604, 010 (2006) [arXiv:astroph/0601617].
 (32) R. Easther, J. T. . Giblin and E. A. Lim, “Gravitational wave production at the end of inflation,” Phys. Rev. Lett. 99, 221301 (2007) [arXiv:astroph/0612294].
 (33) J. GarciaBellido and D. G. Figueroa, “A stochastic background of gravitational waves from hybrid preheating,” Phys. Rev. Lett. 98, 061302 (2007) [arXiv:astroph/0701014].
 (34) J. GarciaBellido, D. G. Figueroa and A. Sastre, “A gravitational wave background from reheating after hybrid inflation,” Phys. Rev. D 77, 043517 (2008) [arXiv:0707.0839 [hepph]].
 (35) J. F. Dufaux, A. Bergman, G. N. Felder, L. Kofman and J. P. Uzan, “Theory and numerics of gravitational waves from preheating after inflation,” Phys. Rev. D 76, 123517 (2007) [arXiv:0707.0875 [astroph]].
 (36) C. Caprini, R. Durrer and G. Servant, “Gravitational wave generation from bubble collisions in firstorder phase transitions: An analytic approach,” arXiv:0711.2593 [astroph].
 (37) G. N. Felder and I. Tkachev, “LATTICEEASY: A program for lattice simulations of scalar fields in an expanding universe,” arXiv:hepph/0011159.
 (38) G. N. Felder, “CLUSTEREASY: A program for simulating scalar field evolution on parallel computers,” arXiv:0712.0813 [hepph].
 (39) A. Nordlund and P. Padoan, “Density PDFs of supersonic turbulence,” arXiv:astroph/9810074.
 (40) R. Micha and I. I. Tkachev, “Relativistic turbulence: A long way from preheating to equilibrium,” Phys. Rev. Lett. 90, 121301 (2003) [arXiv:hepph/0210202].
 (41) R. Micha and I. I. Tkachev, “Turbulent thermalization,” Phys. Rev. D 70, 043538 (2004) [arXiv:hepph/0403101].
 (42) J. M. Bardeen, “Gauge invariant cosmological perturbations,” Phys. Rev. D 22, 1882 (1980).
 (43) J. M. Bardeen, J. R. Bond, N. Kaiser and A. S. Szalay, “The statistics of peaks of Gaussian random fields,” Astrophys. J. 304, 15 (1986).
 (44) F. Bernardeau and L. Kofman, “Properties of the cosmological density distribution function,” Astrophys. J. 443, 479 (1995) [arXiv:astroph/9403028].
 (45) J. R. Bond, L. Kofman and D. Pogosian, “How filaments are woven into the cosmic web,” Nature 380, 603 (1996) [arXiv:astroph/9512141].
 (46) J. W. York, “Gravitational degrees of freedom and the initialvalue problem,” Phys. Rev. Lett. 26, 1656 (1971).
 (47) M. Shibata and T. Nakamura, “Evolution of threedimensional gravitational waves: Harmonic slicing case,” Phys. Rev. D 52, 5428 (1995).
 (48) T. W. Baumgarte and S. L. Shapiro, “On the numerical integration of Einstein’s field equations,” Phys. Rev. D 59, 024007 (1999) [arXiv:grqc/9810065].
 (49) F. Pretorius, “Numerical relativity using a generalized harmonic decomposition,” Class. Quant. Grav. 22, 425 (2005) [arXiv:grqc/0407110].
 (50) M. Patra and M. Karttunen, “Stencils with isotropic discretization error for differential operators,” Num. Meth. for PDEs 22, 936 (2005).
 (51) U. L. Pen, D. N. Spergel and N. Turok, “Cosmic structure formation and microwave anisotropies from global field ordering,” Phys. Rev. D 49, 692 (1994).
 (52) V. A. Belinsky, I. M. Khalatnikov, L. P. Grishchuk and Y. B. Zeldovich, “Inflationary stages in cosmological models with a scalar field,” Phys. Lett. B 155, 232 (1985).
 (53) D. Polarski and A. A. Starobinsky, “Semiclassicality and decoherence of cosmological perturbations,” Class. Quant. Grav. 13, 377 (1996) [arXiv:grqc/9504030].
 (54) A. Fournier, D. Fussell and L. Carpenter, “Computer rendering of stochastic models,” Communications of the ACM 25, 371 (1982).
 (55) G. Efstathiou, M. Davis, C. S. Frenk and S. D. M. White, “Numerical techniques for large cosmological Nbody simulations,” Astrophys. J. Suppl. 57, 241 (1985).
 (56) Y. Hoffman and E. Ribak, “Constrained realizations of Gaussian fields  A simple algorithm,” Astrophys. J. 380, L5 (1991).
 (57) Y. Hoffman and E. Ribak, “Primordial Gaussian perturbation fields  Constrained realizations,” Astrophys. J. 384, 448 (1992).
 (58) J. Salmon, “Generation of correlated and constrained gaussian stochastic processes for Nbody simulations,” Astrophys. J. 460, 59 (1996).
 (59) U. L. Pen, “Generating cosmological gaussian random fields,” Astrophys. J. 490, L127 (1997) [arXiv:astroph/9709261].
 (60) E. Bertschinger, “Multiscale Gaussian random fields for cosmological simulations,” Astrophys. J. Suppl. 137, 1 (2001) [arXiv:astroph/0103301].
 (61) S. Plaszczynski [PLANCKHFI Collaboration], “Generating long streams of noise,” Fluct. Noise Lett. 7, R1 (2007) [arXiv:astroph/0510081].
 (62) H. Bateman, “Higher transcendental functions (Volume III),” McGrawHill (1955).
 (63) A. V. Frolov and L. Kofman, in preparation.