Cosmological Inflation in N-Dimensional Gaussian Random Fields with Algorithmic Data Compression

There is considerable interest in inflationary models with multiple inflaton fields. The inflaton field $\boldsymbol\phi$ that has been postulated to drive accelerating expansion in the very early universe has a corresponding potential $V$, the details of which are free parameters of the theory. We consider a natural hypothesis that $V$ ought to be maximally random. We realize this idea by defining $V$ as a Gaussian random field in some number $N$ of dimensions. Given a model that statistically determines the shape of $V$, we repeatedly evolve $\boldsymbol\phi$ under random potentials, cataloging a representative sample of trajectories associated with that model. On anthropic grounds, we impose a minimum with $V=0$ and only consider trajectories that reach that minimum. We simulate each path evolution stepwise through $\boldsymbol\phi$-space while simultaneously computing $V$ and its derivatives along the path via a Gaussian random process. When $N$ is large, this method significantly reduces computational load as compared to methods that generate the potential landscape all at once. Even so, the covariance matrix of constraints on $V$ can quickly become large and cumbersome. To solve this problem, we present data compression algorithms to prioritize the necessary information already simulated, then keep an arbitrarily large portion. With these optimizations, we simulate thousands of trajectories, extract matter, tensor, and isocurvature spectra from each, and then assemble statistical predictions of these quantities through repeated trials. We find that the Gaussian random potential is a highly versatile inflationary model with a rich volume of parameter space capable of reproducing modern observations.


INTRODUCTION
Characterization of plausible inflationary models is an extremely high priority in cosmology (Dodelson 2003).Since inflationary theory was developed in the late 1970s and 80s, the cosmological community has sought to constrain the parameter space of possible inflationary models (e.g., Langlois 2010).Some parts of this parameter space -in particular, some choices for the inflaton potential V -have been explored in exhaustive detail, but the vastness of the parameter space defies complete exploration.
Potentials with a large number of inflaton degrees of freedom N (Salopek et al. 1989) are of particular interest, as they may be a natural consequence of string theory (Baumann & McAllister 2015).Even if we lay string theory aside, it is extremely plausible to imagine that many additional fields may occur as we extend our model of particle physics to the vastly higher energies at which inflation is thought to have occurred.Fully exploring a theory with N ≫ 1 inflaton fields, with an unknown potential V (ϕ (0) , . . ., ϕ (N −1) ), is a computationally daunting task.
Rather than specifying a fixed potential a priori, it may be natural to consider models inspired by the "landscape" of string theory, in which V is taken to be a realization of a random process, often taken to be a Gaussian process with a Gaussian correlation function.One can imagine making predictions for this model by generat-ing many random realizations of the potential, simulating the dynamics of the inflaton fields, and using the result to make predictions for observable quantities (scalar and tensor power spectra, non-Gaussianity, etc.)In a high-dimensional space, simulating a realization throughout an entire box is computationally intractable.Fortunately, we need the simulation only along the inflaton path.
Other methods besides RMT have been used to explore multifield inflation in random potentials.Some involve simulating the potential in the vicinity of a single point (Masoumi et al. 2017;Bjorkmo & Marsh 2018).Other methods involve random walks (Tye et al. 2009), Bayesian networks (Price et al. 2016), and machine learning (Rudelius 2019).
In this paper, we present a method of simulating trajectories in a random landscape directly, without recourse to the approximations involved in techniques such as RMT.At each step along the numerical solution of the inflaton equations of motion, we simulate the potential, taking into account the information provided by the simulations at previous steps along the path.Similar methods are found in Blanco-Pillado et al. (2020).
At each step in our simulation, the values of V and its gradient are generated at the next point along the inflaton trajectory.These values are drawn from the conditional probability distribution that takes account of both the correlation function of the Gaussian random process from which the full function V is assumed to be drawn and all of the previous potential and gradient values that have been sampled.The result is equivalent to simulating the full potential function in advance and numerically solving for the trajectory but is much less computationally expensive.
In order to keep the computation from becoming intractable, we periodically "forget" information that no longer significantly influences the simulation.In principle, this forgetting step breaks the precise equivalence between our method and simulating the entire potential in advance, but we are able to quantify the degree of approximation involved and show that it is small.We illustrate our method by computing many simulated trajectories and using the results to predict distributions of various observable quantities, including the adiabatic, isocurvature, and tensor power spectra.In these illustrations, we consider relatively low values of N and focus on a particular portion of parameter space.With further code optimization and/or the use of more computation time, we expect to extend this work to larger numbers of dimensions and a wider range of input parameters.
A central feature of our approach is that we do not restrict our attention in advance to certain categories of trajectory (e.g., slow-roll slow turn or rapid turn).Rather, we simulate trajectories drawn at random from the probability distribution specified by the Gaussian random process for V , and let this determine the frequency with which these different types of trajectory appear.The one exception is that we impose a priori the constraint that inflation ends, by only considering trajectories that end at a minimum at which V = 0.This constraint is justified on anthropic grounds (e.g., Barrow & Tipler 1988;Efstathiou 1995;Susskind 2007): in the landscape models under consideration, most of the universe inflates forever or recollapses rapidly, but we live in one of the rare patches that did neither.
The paper is structured as follows: requisite multifield inflationary theory is presented in the remainder of Section 1, the generation of Gaussian random potentials is formalized in Section 2, our original code and data compression algorithms are detailed in Section 3, results are given in Section 4, and we conclude with Section 5.

Multi-field Inflation
As customary in inflationary theory, we suppose there existed some "inflaton" field ϕ(t, x) which permeated the early universe with time evolution governed by its associated potential function V (ϕ).We assume that the inflaton field was spatially homogeneous with small per-turbations, i.e., ϕ(t, x) ≈ ϕ(t) + δϕ(t, x) In general, ϕ = (ϕ (0) , ϕ (1) , ..., ϕ (N −1) ) is a vector of N components while V is scalar-valued.In addition to its existence, we assert that the inflaton potential dominates the energy density of the early universe so that its characteristics may determine the cosmic expansion history at early times.The Klein-Gordon equation for the evolution of where overdots signify derivatives with respect to cosmic time t, H is the Hubble parameter, and Planck units c = ℏ = 8πG = 1 are implemented throughout.We assume a flat field-space metric (e.g., Langlois 2010); generalization of our method to any given metric should be possible.
Equation ( 1) is written componentwise as for α = 0, . . ., N − 1.As the well-known analogy goes, this system of equations is equivalent to a ball rolling in the potential V with a friction force given by the Hubble parameter.
The Hubble parameter is related to the field by the Friedmann equation, Given the potential, the second-order differential equation (2) can be solved numerically with initial conditions ϕ i = ϕ(t i ) and φi = φ(t i ).Often, the initial rate of change is set to slow-roll equilibrium so that φi = 0.By (1), the approximation is Furthermore, if | φ| 2 is assumed to be negligibly small, then (3) reduces to 3H 2 ≈ V and the slow-roll approximation is φi For numerical and analytic purposes, it is often useful to measure time in terms of the expansion of the universe (e.g., Price et al. 2015).During inflation, the scale factor a(t) grows approximately exponentially, so a natural unit of time is the "number of e-folds" where t i is some arbitrary start time of inflation.Since H = ȧ/a, the left-hand side of (5) then simplifies to dϕ i /dN e and, under slow roll, some algebra transforms (2) into in which we have defined the slow-roll parameter ϵ as signaling the end of inflation at t e when ϵ ↑ 1.Note that H can be computed at any time step by Once numerical integration finishes at ϵ = 1, the values of ϕ, φ, V , ∇ ϕ V , ϵ, and H are known at discrete times between t i and t e .We refer to this set of information as the inflaton field trajectory.
1.2.Power Spectra Various observable quantities can be extracted from the trajectory, such as the power spectra of adiabatic curvature (scalar) fluctuations P S , tensor fluctuations P T , isocurvature fluctuations P iso and the non-Gaussianity of the CMB (Price et al. 2016).In this paper, we will focus on the power spectra.
Following Price et al. (2015), curvature perturbations evolve with the mode matrix ψ, an expansion of the comoving fluctuations δϕ in the inflaton field.For any given mode k, the mode matrix obeys (10) where the coupling matrix C is defined as Given conditions on ψ and dψ/dN e at t i , equation ( 10) can be solved numerically up through t e .At any given time, P S (k) is computed as where Ψ = φ • ψ/| φ|. and * is the conjugate transpose.
To be clear, the power spectrum, as with the mode matrix, is a function of both wavenumber and time.We are interested in the shape of P S (k) at the end of inflation, as this spectrum can be used to compare with observations. 1 We assume standard, rapid reheating at the end of inflation.We expect the power spectrum to resemble a power law of amplitude A S and index n S around the pivot scale k 0 = 0.05 Mpc −1 : For any reasonable inflationary model, P S can be evolved through time at several modes evenly log-spaced around k 0 via (10), then fit to (13) at t e .The best fit parameters A S and n S can then be used to gauge the  (Akrami et al. 2020).
viability and predictive power of the choice of inflaton potential and initial conditions.
Similarly, the power spectrum of isocurvature perturbations can be computed at any point in time as where Ψ α = ŝα • ψ for each unit vector ŝα orthogonal to φ.As with the curvature power spectrum, we sample P iso around k 0 , then fit to a power law, The isocurvature spectrum amplitude is typically expressed as the isocurvature-to-scalar ratio β iso (k 0 ) ≡ A iso /(A iso + A S ), so we will report extract and report this quantity from our own model.Following Langlois (2010), the power spectrum of tensor perturbations is easier to compute: The one-to-one and onto function t(k) yields the time at which a scale of magnitude k last exited the horizon, most easily expressed as its inverse, k(t) = a(t)H(t).Thus, in this definition, P T (k) gives by default the power at the end of inflation.As with the curvature spectrum, the tensor spectrum can be sampled at several wavenumbers around k 0 and fit to another anticipated power law, The tensor perturbation amplitude is often expressed in terms of the tensor-to-scalar ratio r, which is defined as the power law ratio for any given k, We will quote tensor-to-scalar ratios at the scale k t = 0.002 Mpc −1 .Constraints on A S , n S , r(k t ), and β iso (k 0 ) have been found by recent Planck missions and are summarized in Table 1 (Akrami et al. 2020).
When computing all of the power spectra described above, we follow the procedures used in the MultiMode-Code software (Price et al. 2015), and we have confirmed that the two methods give matching results in particular cases such as a quadratic potential.
The review in Section 1 was independent of the shape of inflaton potential, but a definition of V (ϕ) is the heart of every inflationary model.In this paper, we assume that V (ϕ) is a realization of an N -dimensional Gaussian random field with mean zero, coherence scale s, and energy scale V ⋆ (Tye et al. 2009).These three parameters (s, V ⋆ , N ) do not uniquely define a potential, but rather characterize the statistical properties of a family of random potentials.The coherence scale is a characteristic distance in ϕ-space below which values of V are strongly correlated.The inflationary energy scale V ⋆ is the standard deviation of V .In order to produce scalar fluctuations at levels that approximately match observations, V ⋆ should be of order 10 −5 E 4 Pl .N is the dimension of V , or the number of components of ϕ.

Generating Correlated Points
In order to simulate a trajectory, one could begin by simulating the potential V throughout an N -dimensional box of ϕ-space in which the trajectory is expected to lie, but for large N this is extremely inefficient.It is far more efficient to generate V only along the trajectory; potential values elsewhere are irrelevant and unnecessary for generating new values.We describe and implement this method of simulating V along the path of ϕ in the following sections.
Generating points on a Gaussian random field is a constrained random process; every new point is generated randomly constrained by every previously generated point.This is necessary to ensure that the potential has the desired correlation structure -in particular, that it is appropriately smooth on scales less than s.
Suppose that V is known at P individual "past" points in ϕ-space, and that those known values are grouped into a vector v P = (V (ϕ 0 ), V (ϕ 1 ), . . ., V (ϕ P −1 )).Now, we want to generate new potential values at F "future" points, constrained by the information in v P .First, we need to define the covariance matrix Γ of constraints on V as an (P + F ) × (P + F ) matrix with elements where i, j index v = (v P , v F ).This symmetric, positivedefinite matrix is subdivided into blocks corresponding to the type of information stored: The upper left block Γ PP contains covariances of past points, the offdiagonal blocks Γ PF = Γ T F P contain crosscovariances between past and future points, and Γ F F contains covariances of future points (Press et al. 2007).The vector v F is drawn from a Gaussian distribution (Hoffman & Ribak 1991;Rybicki & Press 1992) and covariance matrix where Γ C = L C L T C is the lower-triangular Cholesky decomposition and y is a vector of F independent standard normal values drawn from N (0, 1).
As an illustrative example, suppose we sample the potential at ϕ 0 to be V (ϕ 0 ) = aV ⋆ , and we seek to generate a second point, V (ϕ 1 ), constrained by V (ϕ 0 ).If the points are separated in ϕ-space by |ϕ 0 − ϕ 1 | = ∆ϕ, the covariance matrix is Since the number of future points is F = 1, the mean is a scalar, as is the conditional covariance matrix, The new point is then generated as , where y is drawn from a standard normal distribution.Figure 1 plots µ as a function of ∆ϕ with one and twosigma regions for a = 0, 1, and 2. As ∆ϕ grows, µ trends from V (ϕ 0 ) toward zero, the unconstrained mean, and the standard deviation increases from zero toward V ⋆ .In other words, the value of V 1 is highly constrained around V 0 when ϕ 1 is close to ϕ 0 , but the constraint weaken as the separation grows.
We also store derivatives of the potential and their covariances with other information.For example, if one element v i of v is a derivative with respect to ϕ (α) , ∂V /∂ϕ (α) , and another element v j is a value of V , then All correlations between first and second derivatives of V can be derived in a similar way, albeit with tedious algebra.
We emphasize that generating values of V and/or its derivatives along a trajectory by keeping track of the full conditional probability distribution at each point is precisely equivalent to generating an entire sample function V in advance and then sampling it at the given points.The reason is simply that the joint probability distribution of a set of n random variables is the same as the product of the appropriate conditional probabilities: p(x 1 , x 2 , . . ., x n ) = p(x 1 )p(x 2 |x 1 )p(x 3 |x 1 , x 2 ) . ... If we wanted to, we could go back at any later time and sample other quantities from the appropriate conditional probability distribution.For example, if we have simulated V and ∇V at a series of points, we could later go back and ) values as a function of ∆ϕ, following the in-text example.Given V (ϕ 0 ) = aV⋆ and seeking to generate V (ϕ 1 ) at some distance ∆ϕ = |ϕ 0 − ϕ 1 | away, the figure shows the mean µ (black line) and 1σ, 2σ intervals of the distribution from which V (ϕ 1 ) is drawn for a range of separations.When ∆ϕ is small, the points are close together and V 1 will be highly constrained around V 0 .However, the influence of V 0 on V 1 decays as ∆ϕ grows.

The Minimum Condition
Generic inflaton trajectories do not lead to universes like our own: the solution settles into a minimum of V that is generically quite far from zero, leading to a universe that either inflates forever or rapidly recollapses.When comparing with observation, we wish to look at models that are anthropically constrained (e.g., Barrow & Tipler 1988;Efstathiou 1995;Susskind 2007) to end in a minimum whose value is extremely close to zero.
We implment this anthropic constraint by imposing the following N 2 + N + 1 conditions on V before the trajectory is solved: These conditions ensure that the potential is zero, flat, and concave up at the origin.More precisely, we want the Hessian H to be positive-definite at the origin.Since there are no preferable locations in a Gaussian random field, we could have chosen any location for the minimum, but chose the origin for convenience.Constraints (32) and ( 33) are equivalent to diagonalizing H and forcing its diagonal elements to be positive.Physically, we reorient our axes to align with the principal directions of curvature.Figure 2 shows an example Gaussian random poten-tial subject to the minimum conditions at the origin.Smoothly connected valleys and peaks are visible in the landscape.Notice that there are some locations from which ϕ would not evolve toward the origin (e.g. the southwest and north, on the left panel), instead meandering away to some adjacent minimum.We discard these runs as unphysical universe on anthropic grounds since V e ̸ = 0.However, the likelihood that ϕ converges to the minimum at the origin increases as |ϕ i | decreases; in this case, |ϕ i | < s almost guarantees convergence.Note that at V = 0, our simulated minima are likely to be shallower than more-common minima at V < 0. It is likely, especially with more inflaton components, that at the concavity is only slightly greater than zero along at least one principal direction.For these common potentials, the minimum appears more like a saddle point, and ϕ "momentum" can easily carry the inflaton field through the shallow minimum and out the other side (also leading to a discarded run).

COMPUTATIONAL PROCESS
We developed the bushwhack code to solve the equations of motion for the inflaton field in an N -dimensional Gaussian random potential.Careful measures were implemented to simplify and compress data related to the random generation of V , the most computationally intensive feature.

The Program
Figure 3 shows the general progression of the bushwhack program, from the input of initial conditions to the output information about predicted observable quantities.We detail each step below.
Supply initial conditions.First, we characterize the potential by choosing its parameters, s, V ⋆ , and N .The initial ϕ (α) values are determined by randomly selecting a location near the origin at some user-inputted radius ϕ i ≡ |ϕ i |.This random choice of initial ϕ is necessary to collect a representative sample of trajectories.Given that we diagonalize the Hessian matrix at the origin (see (1) [s] Fig. 2.-A 2-dimensional Gaussian random potential V generated at nodes on a 40×40 grid evenly spaced in ϕ (0) , ϕ (1) .Left: a "birds-eye view" of the potential isocontours, colored by potential value.Right: the same potential function as a traditional 3-D projection, with V as the height.In both panels, the imposed minimum conditions are clearly visible.Note that we do not employ these techniques of generating potential values throughout an entire rectangular domain when solving the trajectory, but do so here for illustrative purposes.Fig. 3.-Flowchart for bushwhack code illustrating the chronology of computational processes involved in solving for the trajectory of ϕ and extracting observable quantities.The user supplies initial conditions, a minimum is seeded at the origin (an anthropic constraint), the inflaton field is evolved until a stopping criterion triggers, and curvature, tensor, and isocurvature spectra are extracted.
Section 2.2), forcing ϕ i = (ϕ i , 0, 0, . . ., 0), for instance, would yield paths which always start near a principal axis of curvature.The inflaton "velocity" φ need not be user-supplied, since its value will be automatically computed from V i and ∇ ϕ V i via equation ( 5).
Seed the minimum.Before solving for a trajectory, the code enforces conditions ( 30)-( 33) to allow for the possibility that the solution is physically plausible.In particular, condition (32) is imposed by randomly simulating N independent second derivative values, taking their absolute values to force them positive, and storing them as the diagonals of the Hessian of V at the origin.At this step, many arrays are initialized such as Γ, v, ϕ, and dϕ/dN e .The imposed constraints at the origin will remain the first N 2 + N + 1 elements of v throughout the simulation and constrain the rest of the potential via a corresponding upper-left block of Γ = ⟨vv T ⟩.
Simulate-solve cycle.Beginning at ϕ i , N + 1 data points are randomly generated, V i and ∇ ϕ V i .Those values are then substituted in equation ( 7) and the field is solved numerically through one time step to ϕ i+1 via a fourth-order Runge-Kutta solution algorithm.Then V i+1 and ∇ ϕ V i+1 are generated and so continues the "simulate-solve cycle".At any given time step, the potential is generated exactly along the path evolution of ϕ, constrained by all previously generated data stored in Γ.In the language of Section 2.1, past constraints define Γ PP , which gets larger and more burdensome as the simulation goes on, while Γ F F remains an (N + 1) × (N + 1) square representing the new potential and gradient values at the current position.
The most computationally expensive steps in this cycle are those involving Γ −1 PP (Equations 22 and 23).Because this matrix is symmetric positive definite, these are done via Cholesky decomposition and solving systems of equations with a triangular matrix, which is much faster and numerically stable.This decomposition can be updated incrementally to account for the new rows and columns added in each time step; it does not have to be done from scratch.
It is also worth noting that this matrix has a large condition number : many of the constraints are nearly redundant, leading to nearly null eigenvectors.To ensure numerical stability, one can add a small multiple of the identity to the matrix with negligible impact on subsequent potential generation.
Stopping conditions.Multiple stopping conditions are monitored at every time step to supervise the evolution of ϕ.The first slow-roll parameter, ϵ, turns out to be an excellent indicator of the end of inflation.ϵ(t i ) is a small fraction of unity, but ϵ ↑ 1 as ϕ accelerates into the minimum at the origin.This condition triggers before ϕ oscillates around the minimum in the potential, so no information about those oscillations are stored.However, potential minima in Gaussian random fields are expected to resemble N -dimensional quadratics, so the oscillations may be easily approximated without further simulation.Additionally, we define the "inwardness" parameter ι as the alignment between the slopes of V and the minimum at any given location, or where hats represent unit vectors.If ι < 0 at any point during the simulation, then V is sloping away from the origin and it is extremely unlikely that ϕ will converge to the V = 0 minimum.In the cases we examine below, we have verified that adopting this criterion results in rejecting a small fraction (≲ 1%) of trajectories that do in fact converge to the correct minimum, and hence makes a negligible difference to the final probability distributions of observables.However, as N increases, there is a higher probability that ϕ takes significant (> 90 • ) turns before converging to the minimum, so we may have to relax this condition in future work.As a last resort, V itself is monitored.If V < 0 at any point, then ϕ is in the process of converging to a negative (and, therefore, unphysical) minimum away from the origin.All three conditionsϵ ↑ 1, ι ↓ 0, and V ↓ 0 -serve to halt the simulate-solve cycle at the precise moment when either inflation ends or the evolution goes awry.The accumulated data is then bunched into an object and analyzed.Stored data includes the inflaton field, the potential, and stopping conditions at every location along the trajectory as well as the full covariance matrix.
Plausibility check.Once the equations of motion are solved, bushwhack makes an effort to determine whether the solution is a remotely plausible representative of our own universe.The verifying questions that the code asks are as follows: 1. Did the solution fail by diverging?In other words, was ι(ϕ) < 0 at any point?Even if the solution fails a plausibility check, the information is available for manual analysis.
Extract observables.If a trajectory is found to converge to the origin with enough inflation, scalar, tensor, and isocurvature power spectra are extracted via methods detailed in Section 1.2.To sample the curvature power spectrum, note from (11) that it is now required to know the second derivatives of V along the path, so those are simulated sparsely along the path using Γ.It is worth emphasizing here that even with the same set of potential parameters and initial conditions on the trajectory, infinitely many random potentials and trajectories are possible, each of which yield different power spectra.It is only through repeated trials that we compile the statistical predictions of observable quantities from a certain parameter set.Plots of data computed from one sample trajectory are presented in Figure 4.

Algorithmic Data Compression
Even when V is simulated only at points along the path evolution of ϕ, the covariance matrix Γ can become large and ill-conditioned.Since we add N + 1 constraints with each time step, the covariance matrix Γ after T time steps has over (N + 1) 2 T 2 elements and accumulates over (N + 1) 2 (2T + 1) new elements each step.Both populating the matrix and performing operations with it, such as Cholesky decompositions, can quickly become burdensome, so compiling a representative sample of solutions with high N without compression methods is computationally difficult.We present two methods of data compression that can be applied iteratively throughout a simulation to limit the size of Γ while maintaining reasonable levels of accuracy in predictions of spectral quantities.

"Unprincipled" Forgetting
The trajectory of ϕ evolves smoothly down the slopes of V ; at some time step T , the values nearest ϕ T are usually ϕ T −1 , ϕ T −2 , etc.By Equation ( 20), V T and its derivatives are thus most highly constrained by V T −1 , V T −2 , etc., and their derivatives.Similarly, if ϕ travels long distances during its evolution, then V T is unlikely to be strongly constrained by early data.How much of the information contained in Γ contributes a Left: the inflaton field meanders in the direction of ∇ ϕ V (hatches along the curve) from some initial ϕ i = s to the seeded minimum at the origin.Middle: slow roll parameter ε and inwardness ι, tracked to determine when to stop evolving ϕ.In this case, ϵ triggers when it accelerates to 1, a sign that the trajectory converged to the origin.The inwardness parameter detects that ϕ moved nearly perpendicular to the origin at first, but never triggers.Right: power-law fits to curvature (top) and tensor (bottom) power spectra.Analytic curve for P T (k) is shown in blue.Note that the approximation is only valid near k 0 = 0.05 Mpc −1 ≈ 1.31 × 10 −58 l −1 Pl .
non-negligible constraint on new data?It is possible that early, weak constraints can be deleted while retaining the recent, relevant core of Γ and therefore increasing the efficiency of solving for ϕ.We detail this forgetting method below, labeled "unprincipled" because the decision about which entries in v to forget is made in an ad hoc manner.Assume that ϕ is evolved through T − 1 timesteps and that Γ is large and ill-conditioned.The strength of previous constraints on V T can be quantified in the sum of the eigenvalues (i.e., the trace) of Γ C : lower values of Tr(Γ C ) correspond to stronger constraints (lower variance in generated points).Consider removing rows/columns of Γ corresponding to early data, V i , V 1 , V 2 etc., one at a time, leaving a modified covariance matrix Γ, and recomputing the constrained covariance matrix ΓC after each removal.We carried out this calculation with 50 converged trajectories in different random potentials.Figure 5 shows the strength of constraints from ΓC relative to those from Γ C computed from the full covariance matrix as a function of the amount of data deleted.As Γ gets smaller (moving from right to left on the plot), Tr( ΓC ) increases (higher on the vertical axis).In most cases, half of the rows in Γ can be deleted before Tr( ΓC ) changes by more than 10% of itself.However, the conditional covariance matrix does not directly tell us how deleting old data affects the accuracy of computed spectral quantities.
To jointly quantify the amount of time saved and the accuracy of the unprincipled forgetting method, we ran 50 simulations with the bushwhack code, parameters (s, V 2 ⋆ , N, ϕ i ) = (30, 5 × 10 −9 , 2, 0.8s), without any compression methods.For each run, we saved the total simulation time t sim , predicted value for A S , and random state, then reran the simulation (with the same random state) 5 times, each executing unprincipled forgetting steps once Γ reached some variable maximum size.For each rerun we, again, stored the compressed simulation time tsim and compressed prediction for the curvature spectral amplitude ÃS .The left panel of Figure 6 demonstrates that the fractional time taken tsim /t sim decreases as a power law with t sim .The time saved can achieve values of 1/5 with conservative compression and 1/50 with high compression while simultaneously predicting ÃS to be very nearly correct.The right panel of Figure 6 confirms that higher degrees of compression tend to produce less accurate predictions, but assures that these inaccuracies are small enough at moderately high levels.Given the speed-up, accuracy, and ease of implementation of unprincipled forgetting steps, we incorporate them permanently into the bushwhack code to execute at a maximum Γ size of 1000 2 .

"Principled" Forgetting
In the previous subsection, we described an method for increasing computational efficiency that involved compressing Γ under the assumption that earlier points weakly constrain newly generated data.This method was labeled "unprincipled" because it involved the reasonable, but ad hoc choice to delete points based on their distance from the point of interest.In this section, we Fig. 5.-Constraining power lost from the covariance matrix Γ as a function of kept information.As more information is removed from earlier in the simulation (moving left), the constraints on a new potential value become more relaxed (the trace of the conditional covariance matrix ΓC increases).However, the constraining power of Γ is not equally distributed in its rows: half of all rows/columns can be removed before Tr( ΓC ) increases by 10%.
describe a "principled" method of compressing Γ that allows for direct control over the portion of trace kept.We leave its implementation in random field inflationary trajectories to future work (see Section 5.1).
We seek to compress Γ so as to maximize its constraining power while minimizing its overall size.When v P gets large, we seek to find a compression matrix A such that ṽP = Av P , where ṽP is a smaller, compressed version of v P .We can transform blocks of Γ into compressed forms as As mentioned previously, lower values of Tr(Γ C ) correspond with tighter constraints on new potential generation.Thus, to maximize the constraining power of Γ, we need to reduce Tr( ΓC ) by maximizing the trace of its latter term, It is straightforward to show that τ only depends on the row space of A, so we are free to impose the additional constraint that the rows of A be orthonormal with respect to Γ PP .
This constrained maximization problem leads to where Λ is a matrix of Lagrange multipliers.We can once again take linear combinations of the rows of A to diagonalize Λ.Now, each row a of A is a solution to the generalized eigenvalue problem Since τ = i λ i , the optimal compression matrix will comprise some number of rows corresponding to the largest eigenvalues so that the size of ΓPP is much less than Γ PP and τ ≈ Tr(Γ F P Γ −1 PP Γ PF ).Significant compression may be possible with this method, as the density of data points required to resolve the path evolution of ϕ is often greater than what is necessary to describe the potential in a local domain.An overdensity of data points leads to many redundant constraints on new potential generation, and thus a high concentration of τ in the very largest eigenvalues.We leave proper implementation of principled forgetting to future work.

RESULTS
The bushwhack code (see Section 3.1) was parallelized and implemented with unprincipled forgetting (see Section 3.2.1)for execution on the University of Richmond Quark Supercomputing Cluster.As a preliminary investigation, we choose some reasonable parameter sets for the generation of V : The only variable between these sets is N , which may change predictions of observable quantities in nontrivial ways, though varying s and V ⋆ is worthwhile future work.
To gather a representative sample of trajectories from a given parameter set, we catalogue 1000 simulations with ϕ i chosen from each of 12 radial ϕ-shells of thickness 0.1s from 0.3s to 2.1s.That is, we choose the initial distance from the minimum to lie within a certain bin (for instance, ϕ i ∈ [0.8s, 0.9s]), then randomly choose the vector location ϕ i , all 1000 times in each bin.If ϕ i < 0.3s, the trajectory is very unlikely to accumulate enough inflation, and if ϕ i > 2.1s, exponentially unlikely to converge to the origin.

Weighting Predictions by Success Probabilities
We seek statistical predictions for observable quantities for every s, V ⋆ , N , and ϕ i .The simplest way to collect a representative sample of successful trajectories would be to choose ϕ i uniformly in the volume of a large sphere (say, with radius 3s) centered on the minimum at the origin, many times for many random potentials, and save only the successful runs.This would naturally represent all possible successful trajectories and the statistical predictions could be calculated by simple means and variances.However, in this approach, most random ϕ i would be large -toward the outskirts of the sphere -because most of the sphere's volume is at large radii.To obtain a sizeable number of successful trajectories with low ϕ i , the number of total solutions would need to be very high, and thus computationally expensive, especially in high-N potentials.
Instead, we subdivide the ϕ-volume into spherical shells of equal thickness (not equal volume) around the origin and catalog the same number of trajectories with ϕ i in each shell.This allows us to store an equal number of trajectories with ϕ i ≈ 0.4s as with ϕ i ≈ 2.0s.For instance, if parameter set (45) yields especially interesting predictions in the ϕ i ∈ [0.3s, 0.4s] bin, this method ensures there are plenty of simulations with this property Fig. 6.-Quantifying compression and accuracy of "unprincipled" forgetting.Left: A simulation without compression steps that takes a time t sim [seconds] to complete can finish in a fraction of the time after incorporating forgetting steps.That fraction tsim /t sim falls as a power law with t sim .Colors represent the maximum allowed size of Γ before a compression step is executed: lower maximum size means more frequent compression.Right: Even with high compression, simulations accurately reproduce uncompressed predictions.While higher compression is, indeed, less accurate, a moderate level of compression can be achieved while sustaining reasonably small error.46).While A S predictions exhibit significant tension with Planck data, slight variations in s and/or V⋆ may remedy these differences while preserving the agreements seen in n S , r, and β iso .
to examine.However, when computing overall predictions for spectral quantities associated with a particu- -The probability that a trajectory starting at an initial distance ϕ i from the origin will both converge and accumulate enough inflation.This function is necessary when computing the weighted mean prediction for any spectral quantity.lar (s, V ⋆ , N ) parameter set, the collection of trajectories from all shells is not representative of all possible successful trajectories, since we have sampled more densely at small ϕ i .So, we must weight the means and variances computed from each shell accordingly.
If q l is the simple mean of predictions for a spectral quantity q from trajectories with ϕ i in shell l, then the weighted mean µ q considers all shells as where w l is a shell weight.The weights are the probability that a random successful trajectory from a given parameter set has ϕ i ∈ l.We can compute this probability density function by multiplying the probability P S (ϕ i ) that a trajectory succeeds -converges to the origin with enough inflation -with a volume element to account for more volume in shells farther from the ϕ-origin.P S is an interesting function of ϕ i , plotted in Figure 7, and is computed by taking the simple ratio of successful trajectories to total attempted runs.The probability that a trajectory is successful vanishes in both ϕ i ↓ 0 and ϕ i ↑ ∞ limits (because of e-fold and convergence requirements, respectively) and peaks around ϕ i = 0.7s for each considered parameter set.With this probability evaluated in each shell, we interpolate linearly for a smooth PDF and find the weights to be where SA N (ϕ i ) ∝ R N −1 is the surface area of an Nsphere of radius ϕ i .Similar to µ q , the variance of q is σ 2 q = l w l q 2 l l w l − µ 2 q . (49) 4.2.Power Spectra Predictions We now present the statistical predictions for observable quantities computed from a Gaussian random inflaton potential with parameters sets ( 44)-( 46).We assume any error resulting from compression methods is accurately represented in the variance of the outputted spectral quantities.First, we examine Figure 8, which shows weighted, normalized histograms for each of A S , n S , n T , r, and β iso for each parameter set, colored by initial condition.A few observations to note: • Changing N does not significantly affect the mean or variance of any quantity, but weak trends are observed in n S (positive), A T (negative), n T (negative), and r (negative).Future work into higher N is required to confirm the sensitivity of perturbation spectra to the dimension of V .
• Distributions within individual ϕ i shells are not necessarily the same as the combined distribution for the parameter set.This is especially apparent in the low ϕ i regime (darker blues); simulations with low ϕ i predict values of n S , n T , and r that are significantly offset from the overall mean (which is dominated by simulations with intermediate ϕ i ).
• Predicted mean values are in variable tension with Planck observations, detailed in Table 2.For example, A S is lower than the Planck prediction for all N = 1, 2, 3, while n S , r, and β iso predictions are accommodated by all three parameter sets.Note that this does not indicate that the Gaussian random potential is unable to accommodate all Planck predictions, but that the particular parameter sets (44)-( 46) are not the best fit.
• Distributions roughly normal, but some have slight negative skewness (in particular A S and n S ).Thus, the variance may not be the best metric of statistical spread for each quantity.
The predictions visualized in each histogram are quantified in Table 2.We also present a corner plot of these 5 spectral quantities for our N = 3 simulations.These panels further emphasize a few points: • Most of the quantities appear to be uncorrelated, however weak correlations are observed in n T vs. n S (positive) as well as n T vs. r (negative).Some of the weighted KDEs exhibit curving, non-trivial contours.
• Weighting by Equation ( 48) exaggerates the contributions from simulations with ϕ i ∈ [0.5s, 1.2s] and minimizes contributions from those outside this range, leaving tight kernel density estimations around the means.

CONCLUSIONS
We have presented a method for exploration into Ndimensional inflationary models in which the inflaton potential is a Gaussian random field.We employ several techniques to increase computational efficiency and suggest another for implementation in future work.To be specific, we avoid generating the random potential all at once, but efficiently generate along the path evolution of the inflaton field.The covariance matrix of constraints Γ can grow large and unwieldy, especially with large N , so we propose two compression methods that may ease the computational intensity by exploiting redundant constraints within Γ.We implement "unprincipled forgetting" steps-deleting old constraints to shrink Γ-into our solver and compile thousands of simulations to compute statistical predictions of curvature, tensor, and isocurvature power spectra.We list our main conclusions below: • The Gaussian random potential is a highly versatile inflationary model with a rich volume of promising parameter space.Parameters N and ϕ i each change the statistical predictions for perturbation spectra in unique ways.For example, higher N tends to increase n S while decreasing n T , A T , and r and lower ϕ i predicts lower n S and higher n T .
Parameter sets (44)-( 46) have been demonstrated to accommodate modern constraints on n S , r, and β iso , but not A S (see Table 2 and Figures 8, 9).
• Methods such as the one we describe allow for broad exploration of the space of possible outcomes in any given model.In particular, we do not restrict our attention in advance to trajectories that have certain desired properties; rather, we explore the space of all possible trajectories (subject only to the anthropic constraint that inflation end).This is the right approach if one wants to make probabilistic statements about the probability distribution of any given set of observables in a given model.
• Algorithmic data compression (Section 3.2.2) allows for the possibility of significantly more efficient generation of V in high-N regimes.Unprincipled forgetting is limited in that it relies on the assumption that earlier data contributes weaker constraints on new generation and less rigorous in that it involves the deletion of data rather than compression.Principled forgetting methods, if implemented efficiently, could remedy these issues.
The examples we have presented in this paper are intended as illustrations of the method and are far from exhaustive.For example, the models we consider always involve large excursions (significantly greater than the Planck mass) in field space.It is possible that using a fixed effective field theory over such large excursions is incorrect; in any case, it is of interest to explore other     The spectral quantities n S , n T , r, and β iso are shown to be independent of each other or weakly correlated, at best, for this parameter set.
parameter choices.Similarly, we restrict our attention to the case of a flat field-space metric, but our method could be extended to accommodate any given non-flat metric.Most importantly, it is of course of great interest to explore higher values of N , which will be possible with a combination of improvements in efficiency and simply throwing more CPU cycles at the problem.

Future Work
This project could be further explored in a number of ways, ranked below by priority.
1. Explore parameter space.Carry out a large number of simulations with different parameter sets which vary s and V ⋆ in addition to N .These supplemental runs could quantify the degree to which the Gaussian random potential model is able to accurately reproduce modern constraints on primordial power spectra.
2. Implement principled forgetting.Unprincipled forgetting offers significant speed increases, but at the expense of rigorous, algorithmic data compression.Principled forgetting methods detailed in Section 3.2.2have the capacity to provide enormous compression to Γ, possibly enough to efficiently explore string-theory-esque models with N ∼ 100.
3. Incorporate other observables.The spectral quantities computed from our trajectories are but a subset of the independent quantities constrained by modern observations.Incorporating quantities such as the non-Gaussianity of the CMB could further verify whether the Gaussian random potential is a plausible inflationary model of our universe (Renaux-Petel 2009).Computing non-Gaussianity may require modifying our stopping criterion to include oscillations around the potential minimum.
4. Translate to other languages.The bushwhack code is currently written in Python, but could be optimized by translating to more efficient languages such as C or FORTRAN.

Fig. 4 .
Fig. 4.-A sample trajectory in a Gaussian random potential with parameters s = 30, V 2 ⋆ = 5 × 10 −9 , and N = 2.Left: the inflaton field meanders in the direction of ∇ ϕ V (hatches along the curve) from some initial ϕ i = s to the seeded minimum at the origin.Middle: slow roll parameter ε and inwardness ι, tracked to determine when to stop evolving ϕ.In this case, ϵ triggers when it accelerates to 1, a sign that the trajectory converged to the origin.The inwardness parameter detects that ϕ moved nearly perpendicular to the origin at first, but never triggers.Right: power-law fits to curvature (top) and tensor (bottom) power spectra.Analytic curve for P T (k) is shown in blue.Note that the approximation is only valid near k 0 = 0.05 Mpc −1 ≈ 1.31 × 10 −58 l −1 Pl .
Fig.7.-The probability that a trajectory starting at an initial distance ϕ i from the origin will both converge and accumulate enough inflation.This function is necessary when computing the weighted mean prediction for any spectral quantity.
Fig.8.-Weighted histograms of observable quantities A S , n S , n T , r, and β iso .Different colors on each histogram correspond to different initial positions ϕ i .Each row of plots corresponds to a particular spectral quantity and each column to a parameter set.Horizontal axes are the same along each row.

Fig. 9 .
Fig.9.-Corner plot of histograms and Gaussian kernel density estimations (KDEs) of correlations between spectral quantities for the N = 3 parameter set.On the diagonal, histograms are identical to those plotted in the right column of Figure8.On the off-diagonals, KDEs are computed from weighted scatterplots of ∼ 10 4 simulations: lighter regions indicate the most frequently outputted pairs of parameters while darker regions are uncommon outputs.Two closed black curves are included with each KDE to clarify the subtleties in its shape.The spectral quantities n S , n T , r, and β iso are shown to be independent of each other or weakly correlated, at best, for this parameter set.

TABLE 2
Tabulated results and comparisons.Shown are the weighted means, standard deviations, and tensions with Planck data of various spectral quantities for each parameter set (44)-(