MultiModeCode: An efficient numerical solver for multifield ... - arXiv

8 downloads 216 Views 994KB Size Report
Oct 2, 2014 - we build on ModeCode [19–21], which was developed to test ... Pyflation [25–28] is an object-oriented Python code that uses the ..... where v is given in terms of the momentum density of the stress-energy tensor T ..... the real precision of the computer, which results in a finite difference error in the numerical.
Prepared for submission to JCAP

arXiv:1410.0685v1 [astro-ph.CO] 2 Oct 2014

MultiModeCode: An efficient numerical solver for multifield inflation Layne C. Price,1 Jonathan Frazer,2,3 Jiajun Xu,4 Hiranya V. Peiris,5 and Richard Easther1 1 Department

of Physics, University of Auckland Private Bag 92019, Auckland, New Zealand 2 Department of Theoretical Physics, University of the Basque Country, UPV/EHU 48040 Bilbao, Spain 3 IKERBASQUE, Basque Foundation for Science 48011 Bilbao, Spain 4 Department of Physics, University of Wisconsin-Madison Madison, WI 53706, USA 5 Department of Physics and Astronomy, University College London London WC1E 6BT, UK E-mail: [email protected], [email protected], [email protected], [email protected], [email protected]

Abstract. We present MultiModeCode,1 a Fortran 95/2000 package for the numerical exploration of multifield inflation models. This program facilitates efficient Monte Carlo sampling of prior probabilities for inflationary model parameters and initial conditions and is the first publicly available code that can efficiently generate large sample-sets for inflation models with O(100) fields. The code numerically solves the equations of motion for the background and first-order perturbations of multi-field inflation models with canonical kinetic terms and arbitrary potentials, providing the adiabatic, isocurvature, and tensor power spectra at the end of inflation. For models with sum-separable potentials MultiModeCode also computes the slow-roll prediction via the δN formalism for easy model exploration and validation. We pay particular attention to the isocurvature perturbations as the system approaches the adiabatic limit, showing how to avoid numerical instabilities that affect some other approaches to this problem. We demonstrate the use of MultiModeCode by exploring a few toy models. Finally, we give a concise review of multifield perturbation theory and a user’s manual for the program.

1

Available at www.modecode.org.

Contents 1 Introduction

1

2 Features of MultiModeCode

3

3 A brief review of multifield perturbation theory 3.1 The highlights 3.2 Classical background 3.3 Mode equations 3.4 Power spectra 3.5 δN formalism 3.6 Bundle width

4 4 6 7 8 11 13

4 The method

14

5 Numerical results 5.1 Isocurvature stability 5.2 A case study: Nf -flation with a step

16 16 17

6 Conclusion

19

A Non-canonical kinetic terms

21

B MultiModeCode usage

22

1

Introduction

Many simple models of inflation adeptly reproduce the observed properties of the primordial cosmological perturbations [1–4], predicting a nearly scale-invariant power spectrum and minimal amounts of primordial non-Gaussianity. In the slow-roll, single-field paradigm the predictions of a given model are easily determined as an algebraic function of the field’s potential V and its derivatives in terms of a hierarchy of slow-roll parameters. The resulting observables are simple to compute and easy to interpret. However, relaxing any of the basic assumptions of the slow-roll, single-field models complicates this simple analysis. In particular, for many inflationary scenarios (e.g., multifield inflation, gauge inflation, and non-minimal couplings), the background and mode equations are complex systems of coupled, nonlinear ODEs, making analysis difficult in all but a few cases. Furthermore, while slow-roll, single-field inflation is a simple and easily understood model, it may not necessarily be considered natural in the context of high-energy theories. For example, low energy effective theories derived from string theory generically contain hundreds of scalar fields with complicated interactions, and many theories consider nonminimal couplings to the Ricci scalar (for a recent review, see Ref. [5]). While analytical studies have been able to overcome subsets of these problems, most of the techniques that have been used are situation-specific, which limits their applicability to novel models.

–1–

While significant progress can be made in the slow-roll limit, only numerical techniques can explore the full predictions of more complex inflation models. Even in the purely homogeneous limit, numerically solving the nonlinear Klein–Gordon equation for the homogeneous background fields reveals many interesting features that do not arise in slow-roll analyses, e.g., sensitivity to initial conditions [6–9]. These complications lead naturally to the numerical exploration of inflationary models. In this paper we present and describe MultiModeCode,1 an efficient Fortran 95/2000 package that numerically solves the equations of motion for the background fields and the first-order perturbations for multifield inflation models in which the fields have canonical kinetic terms and are minimally coupled to gravity. MultiModeCode calculates the adiabatic, tensor, and various isocurvature power spectra as a function of scale k, but does not evaluate higher order correlators. If the potential is sum-separable, MultiModeCode uses the solution to the background equations of motion to evaluate the slow-roll δN predictions for the scalar and tensor power spectra and their derivatives near the pivot scale k∗ , also giving the slow-roll results for ns , r, fNL , etc. The code has been extensively tested with various compilers, including the open-source GNU Fortran compiler. Several numerical codes have been developed to study single-field models [10–18]. Here, we build on ModeCode [19–21], which was developed to test single-field inflation models and interfaced with tools such as CAMB [22], CosmoMC [23], and MultiNest [24]. ModeCode was designed for the Bayesian analysis of inflation and used by the Planck collaboration [4] to obtain the posterior probabilities and marginal likelihoods for inflation models. Moving to the multifield case significantly increases the numerical demands on the solver, and puts a premium on efficiency due to the much greater computational resources required by these analyses. A few codes exist to analyze multified models, but the publicly available codes are inadequate for models with many fields and arbitrary potentials. Notably, Pyflation [25–28] is an object-oriented Python code that uses the same method we employ here for solving the perturbation equations, but cannot easily generate large samples due to the speed constraints imposed by a dynamic programming language. This significant extension to ModeCode can be used to study the power spectra of analytically intractable multifield inflationary potentials, and to explore the generic predictions of complex models by marginalizing over large numbers of possible parameters. Complementing currently available codes [25–28], MultiModeCode specializes in obtaining large Monte Carlo samples of initial conditions and parameter prior probabilities. To help users familiarise themselves with MultiModeCode the package includes initial conditions priors used in Refs. [8, 29, 30]. The ability of this code to efficiently generate large Monte Carlo samples has permitted studies of the generic predictions of multifield inflation models with more than 100 fields [30, 31]. In practice, the code can simulate the evolution of the mode equations for O(102 ) fields,2 but will become inefficient for significantly more fields due to the increasing dimensionality of the system, which increases with the number of fields as O(Nf2 ). However, it can efficiently sample the evolution of the background equations of motion for at least O(103 ) fields. While solving just the background equations allows the exploration of background dynamics for such a large number of fields, if the model is sum-separable, then it will also give the slowroll predictions for the adiabatic curvature power spectrum, as well as fNL and τNL , in terms of the δN approximation. This should be valid when the fields are much lighter than H at 1 2

Publicly available at www.modecode.org. Estimates regarding field number are based on Nf -quadratic inflation, which is not numerically intensive.

–2–

horizon crossing and slow-roll holds throughout the duration of inflation. MultiModeCode is released with several example models already implemented and it is straightforward to add to this number. In §5, we demonstrate the features of MultiModeCode with an Nf –flation potential with a sharp step, which we parametrize by a hyperbolic tangent function, following Refs. [10, 32]. We show that, in addition to oscillatory features in the adiabatic curvature power spectrum that are expected from the single-field analysis [10, 33], with more than one field there are also oscillatory features in the isocurvature spectra, which might result in non-trivial evolution of the power spectrum after inflation. We also show that the numerical computation of isocurvature modes results in an inherent numerical instability, since some definitions of isocurvature perturbations involve computing the difference between two quantities that are of the same order of magnitude. This induces a dominant numerical error when these two quantities begin to approach the adiabatic limit. We overcome this problem by implementing a modified definition of isocurvature perturbations [30], which is numerically stable to many more orders of magnitude than some alternative definitions. We also implement a geometrical optics indicator of isocurvature evolution as first presented in Ref. [34]. While this measure only relies on background quantities and also does not suffer from instabilities, as implemented here it does not provide an absolute value of isocurvature, only an indicator of its growth or decay. Finally, in §3 we provide a concise review of multifield perturbation theory with the aim of dispelling misconceptions that exist about this topic, which the enlightened reader can skip.

2

Features of MultiModeCode

We begin by highlighting some of the useful characteristics of MultiModeCode. Speed: The purpose of MultiModeCode is to provide a fast and efficient solver that is well-tested and can be applied to a wide range of possible inflationary scenarios. MultiModeCode is written in Fortran 95/2000, increasing its capabilities relative to existing codes [25–28] and making it tractable to investigate models with many fields or to obtain large Monte Carlo samples from a model’s parameter space. In particular, prototype versions of this program were used in Refs. [30, 31] to analyze large samples of 100-field Nf -monomial inflation. Generality: The code facilitates Bayesian approaches to studying inflation, where the model’s parameters are drawn from prior probabilities from which we can compute a probability distribution for specified observable associated with the model. We consider simple situations, e.g., evolving a model given fixed model parameters and initial conditions, as subcases of the more general Bayesian framework. To facilitate the use of general priors we have implemented the sampling routines in modules which are simple to adapt and restructure for the user’s purposes. Robustness: The program exits gracefully when encountering fatal errors of either a technical or cosmological nature, while also catching specific errors that might only affect one particular configuration of the model. We have extensively checked the program output on various Macintosh and Linux machines with both the gfortran and ifort compilers, and include both a fourth-order Runge-Kutta integrator and an implicit backward-difference formula method, which is suitable for stiff problems.

–3–

Statistics: MultiModeCode provides pivot-scale observables, summarized in Table 1 and can sample the adiabatic and isocurvature power spectra as a function of scale k. We have implemented a variety of numerically stable indicators of the amount of isocurvature present in the system. Slow-roll comparison: If the potential V is sum-separable, MultiModeCode can also calculate observables using the δN approximation, which assumes slow-roll. Since these quantities rely only on solutions of the background equations of motion they are efficient and simple to calculate, scaling with the number of fields as O(Nf ). Consequently, if the model is well-described by the slow-roll approximation between horizon crossing and the end of inflation, computing observables in the δN formalism is efficient and easy.

3

A brief review of multifield perturbation theory

We begin with a short review of first-order, non-interacting multifield perturbation theory before describing MultiModeCode and the dynamics of many-field inflation. There are some substantial differences between single-field and multifield inflation, which we highlight in Section 3.1. Table 1 gives a list of the pivot-scale observables that MultiModeCode computes. There are a few excellent reviews of this topic [26, 35–37] and we particularly recommend Refs. [38, 39] for more information. We first present the nuts-and-bolts of the mode function approach to first-order, multifield perturbations, which is implemented in MultiModeCode. Then we describe the widely-used δN -formalism, which has also been implemented for ease of use and for comparison to the perturbation solutions. 3.1

The highlights

Multifield inflation differs from the single field case in the following important respects. Isocurvature: Multifield inflation generally permits both adiabatic and isocurvature perturbations. Adiabatic perturbations are related by a gauge transformation to the curvature perturbation on comoving hypersurfaces R, while isocurvature perturbations are entropic perturbations between different matter components on flat hypersurfaces. In single-field inflation there is only one matter component, so there are only adiabatic perturbations. Super-horizon evolution: Isocurvature perturbations source adiabatic perturbations, causing them to evolve even on super-horizon scales. While this can generate novel signatures such as non-Gaussianity, this can also be problematic for comparing the predictions of a model with observation: unless isocurvature modes decay into an adiabatic limit before the end of inflation, the curvature perturbation does not become conserved and is thus sensitive to post-inflationary physics. The two-index mode function: With more than one field, either (a) the direct interaction between fields or (b) the gravity-mediated interaction will mix the particle creation and annihilation operators as a function of time [38]. Instead of a single index mode function, we therefore need to solve for a mode matrix ψIJ , where δφI = ψIJ aJ , for Nf annihilation operators aJ .

–4–

Power spectra (PS)

Type

Reference

PR (k) . . . . . . . . . . . .

Adiabatic scalar spectrum . . . . . . . . . . . . . .

Eq. (3.23)

PS (k) . . . . . . . . . . . .

Isocurvature spectrum . . . . . . . . . . . . . . . . . .

Eq. (3.25)

PδP,nad (k) . . . . . . . . .

Non-adiabatic pressure spectrum . . . . . . .

Eq. (3.33)

Pent (k) . . . . . . . . . . .

Entropic spectrum . . . . . . . . . . . . . . . . . . . . . .

Eq. (3.36)

PRS (k) . . . . . . . . . . .

Adiabatic–non-adiab. cross spectrum . . .

Eq. (3.26)

Ph (k) . . . . . . . . . . . .

Tensor spectrum . . . . . . . . . . . . . . . . . . . . . . . .



Observable at k∗

Name

Description

As . . . . . . . . . . . . . . .

Scalar amplitude . . . . . . . . . . . . . . . . . . . . . . .

PR (k∗ )

Aiso . . . . . . . . . . . . . .

Isocurvature ampl. . . . . . . . . . . . . . . . . . . . . .

PS (k∗ )

APnad . . . . . . . . . . . .

Non-adiab. pressure ampl. . . . . . . . . . . . . .

PδP,nad (k∗ )

Aent . . . . . . . . . . . . . .

Entropy ampl. . . . . . . . . . . . . . . . . . . . . . . . . .

Pent (k∗ )

ACross . . . . . . . . . . . .

Cross spectra ampl. . . . . . . . . . . . . . . . . . . . .

PRS (k∗ )

ns . . . . . . . . . . . . . . . .

Scalar spectral index . . . . . . . . . . . . . . . . . . .

D∗ log PR + 1

nt . . . . . . . . . . . . . . . .

Tensor spectral index . . . . . . . . . . . . . . . . . . .

D∗ log Ph

niso . . . . . . . . . . . . . .

Isocurvature spectral index . . . . . . . . . . . . .

D∗ log PS

nent . . . . . . . . . . . . . .

Entropy spectral index . . . . . . . . . . . . . . . . .

D∗ log Pent

nPnad . . . . . . . . . . . . .

Non-adiab. pressure spectral index . . . . .

D∗ log PδP,nad

αs . . . . . . . . . . . . . . . .

Scalar running . . . . . . . . . . . . . . . . . . . . . . . . . .

D∗2 log PR

r .................

Tensor-to-scalar ampl. . . . . . . . . . . . . . . . . . .

Ph (k∗ )/PR (k∗ )

Θ ................

Bundle width . . . . . . . . . . . . . . . . . . . . . . . . . . .

Eq. (3.58)

cos ∆ . . . . . . . . . . . . .

ω-s correlation angle . . . . . . . . . . . . . . . . . . .

Eq. (3.27)

Table 1. Typical observables at the pivot scale k∗ . The derivative D∗ ≡ d/d log k is evaluated at k = k∗ . MultiModeCode can also generate the full power spectra as a function of scale P(k).

Initial conditions dependence: Multifield inflation models have an infinite number of possible inflationary solutions each of which can, in principle yield a different perturbation spectrum. Consequently, the observable spectra for multifield models can depend on their initial conditions in ways that have no direct analogue in slow-roll, single-field models, which have only one possible trajectory in field-space. Inherently stochastic predictions: Even if the potential V is completely fixed, multifield models will give an inherent spread of predictions due to the allowed variance in the fields’ initial conditions. In general, multifield models will predict a variety of spectra, unless the stochasticity in the initial conditions can be controlled a priori.

–5–

3.2

Classical background

Consider Nf scalar fields φI with the matter sector of the action given by   Z 1 4 √ µ I S = d x −g − ∂µ φI ∂ φ − V (φI ) , 2

(3.1)

where we use the Einstein summation convention over repeated indices. Greek indices describe spacetime, going from 0, . . . , 3, upper-case Latin indices describe the number of fields, going from 1, . . . , Nf , and lower-case Latin indices describe space, going from 1, . . . , 3. The field space indices are raised using the Kronecker delta δ IJ . The determinant of the spatial metric gµν is g. In this paper we only consider inflation models with minimal coupling to Einstein gravity and a matter sector described by scalar fields. The current incarnation of MultiModeCode only solves models with canonical kinetic terms, but we give the equations of motion for models with a non-trivial field-space metric in Appendix A. Implementing these general field-space metrics is straightforward since MultiModeCode has been written modularly, but is left for future work. First-order, non-interacting perturbation theory separates the homogeneous, classical background from the spatially-dependent modes as φI (t, ~x) → φI (t) + δφI (t, ~x), where we assume that these two components can be treated independently. The homogeneous background fields obey the Klein–Gordon equations ∂V φ¨I + 3H φ˙ I + = 0, ∂φI

(3.2)

2 = where an overdot indicates a derivative with respect to cosmic time t and we use MPl (8πG)−1 = 1 throughout this paper. The 0-0 Einstein field equation gives the Friedmann equation 1 3H 2 = φ˙ I φ˙ I + V (φI ), (3.3) 2 which can be differentiated with respect to t to yield

2H˙ = −φ˙ 20 .

(3.4)

In Eq. (3.4) we have used the inflaton trajectory velocity, φ˙ 20 ≡ φ˙ I φ˙ I . We can regard the composite field φ0 as the clock of multifield inflation. It is the classical field defined along the inflaton trajectory, and represents the length of the classical field-space path. In practice, if the dynamics are inflationary, it is numerically convenient to evolve the equation with the number of e-folds Ne ≡ ln a(t) as the independent variable, giving d2 φI dφI 1 ∂V + (3 − ) + 2 I = 0, 2 dNe dNe H ∂φ

(3.5)

where we have defined the first slow-roll parameter as ≡−

H˙ 1 dφI dφI = . H2 2 dNe dNe

(3.6)

The Friedmann equation (3.3) can then also be expressed as H2 =

V . 3−

–6–

(3.7)

If V ≈ 0, Eq. (3.7) requires  ≈ 3, which will result in numerical instability whenever we try to set initial conditions that are dominated by their kinetic energy. We side-step this issue by using the cosmic time Eq. (3.2) and H as defined in Eq. (3.3). Solving Eq. (3.5) therefore only requires the initial conditions φI and dφI /dNe , because the dependence on the scale factor a is explicitly removed by the 0-0 Einstein equation (3.7) as a result of assuming a flat FLRW spacetime. As mentioned in §3.1, the perturbation spectrum depends on these initial conditions, which are specified as a prior probability distribution P (φI , φ0I ). 3.3

Mode equations

To obtain the first-order equation of motion for the perturbations δφI , we need to expand the action (3.1) to second-order and include the first-order scalar perturbations to the flat FLRW metric, given by

where

  ds2 = − (1 + 2Φ) dt2 − 2 a2 B,i dt dxi + a2 (1 − 2Ψ) δij − 2∂hi ∂ji E dxi dxj , 1 ∂hi ∂ji E ≡ ∂i ∂j E − δij ∇2 E 3

(3.8)

(3.9)

is trace-free. We choose the spatially-flat gauge, so that Ψ = E = 0, and vary the expanded action δSφ with respect to the perturbations δφI (t, ~x) to get the first-order equation of motion for the free-field perturbations. After Fourier-transforming the scalar perturbations to δφI (k), the mode equations in this gauge are d2 δφI k2 dδφI + (3 − ) + 2 2 δφI + CIJ δφJ = 0, 2 dNe dNe a H where CIJ

∂I ∂J V 1 ≡ + 2 2 H H



dφI dφJ ∂J V + ∂I V dNe dNe



+ (3 − )

dφI dφJ dNe dNe

(3.10)

(3.11)

and ∂I ≡ ∂/∂φI . The equation of motion for the tensor metric perturbations can be derived similarly; since the non-gauge degrees of freedom are massless and only minimally coupled to the matter sector, the resulting equations of motion are identical to the case of single-field inflation. To solve the perturbation equations, it is usually convenient to work with the Mukhanov– Sasaki variable uI ≡ aδφI . The mode equation for uI is d2 uI duI + (1 − ) + 2 dNe dNe



 k2 − 2 +  uI + CIJ uJ = 0 a2 H 2

(3.12)

with CIJ as in Eq. (3.11). Since the mass matrix, defined as m2IJ ≡ ∂I ∂J V , is not necessarily diagonal, the perturbation equations (3.12) mix the annihilation operators for all of the fields [38]. We therefore need to expand each perturbation mode uI (k) and u†I (k) using Nf harmonic oscillators aJ (k): uI (k, Ne ) = ψI J (k, Ne )aJ (k)

and

–7–

u†I (k, Ne ) = ψI J,∗ (k, Ne )a†J (k),

(3.13)

where (†) and (∗) represent Hermitian and complex conjugation, respectively.3 We can then define canonical commutation relations [aJ (k), a†I (k0 )] = (2π)3 δIJ δ (3) (k − k0 ). The mode matrix ψIJ evolves according to  2  k d2 ψIJ dψIJ + − 2 +  ψIJ + CIL ψ LJ = 0. (3.14) + (1 − ) dNe2 dNe a2 H 2 Finding the perturbation spectrum requires setting initial conditions in Eq. (3.14) and using the background equations (3.5) to find the time Ne,k when the mode k leaves the horizon, which also depends on the moment at which the pivot scale k∗ leaves the horizon, N∗ e-folds before the end of inflation. The usual initial condition is the Bunch-Davies state [43], which assumes the field basis has been chosen such that the ψIJ are originally diagonal and sets the initial condition for Eq. (3.14) as if the mode matrix were freely oscillating in Minkowski space. This is wellmotivated, since for modes deep in the horizon k  aH, the mode matrix ψIJ obeys the free wave equation in conformal time d2 ψIJ + k 2 ψIJ = 0, (3.15) dτ 2 where dτ ≡ a dt. If we assume that the mode matrix is initially diagonal at τ = −∞, then Eq. (3.15) yields two solutions  1  (3.16) ψIJ = √ C1 eikτ + C2 e−ikτ δIJ . 2k Translating to e-fold time, the initial conditions can be set by ψIJ

Ne =0

1 = √ (C1 + C2 ) δIJ 2k

and

r dψIJ i k = (C1 − C2 ) δIJ . dNe Ne =0 aH 2

(3.17)

The Bunch-Davies initial condition is equivalent to choosing C1 = 0 and C2 = 1. While only the Bunch-Davies initial condition is implemented in MultiModeCode, non–Bunch-Davies modes could be easily accommodated.4 Although the uI ’s are convenient for short wavelength modes, they grow exponentially after the modes exit the horizon. So once the mode is outside the horizon, MultiModeCode switches from uI to δφI by matching boundary conditions at a time Ne∗ just after horizon crossing with   duI dδφI Ne∗ Ne∗ uI ∗ = e δφI ∗ and δφI + (3.18) =e ∗. dNe Ne∗ dNe Ne Ne Ne

3.4

Power spectra

Unlike single-field inflation, the multifield power spectrum involves contractions of the mode matrix. Using the canonical commutation relations above, the two-point VEV of the field perturbations yields the power spectrum   k3 1 IJ Pδφ (k) = 2 2 ψ IL ψ JL,∗ . (3.19) 2π a 3

An alternative approach is to simply bypass this issue by solving for the field correlation functions directly rather than the individual modes, as in the transport method [34, 40–42]. 4 One would do this by changing the modes’ initial conditions in the set background and mode ic() subroutine in modpk.f90.

–8–

When the field trajectories are not turning, on super-horizon scales the fields φI and their momenta πI commute, indicating that they have transitioned to a regime where Eq. (3.19) can be interpreted as an expectation value over realizations of classical, random fields. To relate this field-space power spectrum to gauge-invariant perturbation variables [44– 46], we first define the curvature perturbation on comoving hypersurfaces R by 1 R ≡ Ψ + ∇2 E + aH (B + v) , 3

(3.20)

where v is given in terms of the momentum density of the stress-energy tensor T µν as  ∂v T i0 ≡ ρ¯ + P¯ δ ij j , ∂x

(3.21)

where ρ¯ and P¯ are the background energy and pressure densities, respectively. If we evaluate Eq. (3.20) on spatially-flat hypersurfaces during inflation, R reduces to R=−

H ωI δφI , ˙ φ0

(3.22)

where ωI ≡ φ˙ I /φ˙ 0 is a basis vector that projects δφI along the direction of the classical background trajectory, given by the solutions to Eq. (3.5). The vector ω ~ and a complementary set of (Nf −1) mutually orthonormal basis vectors ~sK form the kinematic basis [47, 48], where the separation between the adiabatic perturbations in Eq. (3.22) and transverse, isocurvature perturbations is made explicit. Since ω ~ depends on the nonlinear background evolution, in MultiModeCode we find the ~sK numerically by Gram–Schmidt orthogonalization. IJ along the field The adiabatic curvature power spectrum PR is then the projection of Pδφ vector ωI , scaled by the pre-factor in Eq. (3.22), leaving PR (k) =

1 IJ (k). ωI ωJ Pδφ 2

(3.23)

The gauge-invariant scalar density spectrum in Eq. (3.23) is the final result for the adiabatic two-point function to first-order in perturbation theory. Since Eqs. (3.22) and (3.23) are projected along ω ~ , a simple definition for the isocurvature perturbations SK is the orthogonal projection along the ~sK directions SK ≡ −

H J sK δφJ . φ˙ 0

(3.24)

IJ onto all the directions s that are orthogonal to ω and scaling the result By projecting Pδφ K I as in Eq. (3.23), leads to the isocurvature power spectrum: Nf −1 Nf 1 X X K L IJ PS (k) = sI sJ Pδφ (k), 2 KL

(3.25)

IJ

where we have left the summations explicit to indicate that the isocurvature basis vectors are (Nf − 1)–dimensional. We include this definition of isocurvature because it is numerically stable, as we discuss in Sect. 5.1.

–9–

Similarly, we define the adiabatic-isocurvature cross-spectra PRS , which is the crosscorrelation between the comoving curvature perturbation and the total isocurvature perIJ with both ω and the isocurvature basis vectors turbation, given by the contraction of Pδφ sK Nf −1 Nf  1 X X IJ JI PRS (k) = ωI sJK Pδφ + Pδφ . (3.26) 2 K

IJ

Cross-correlations are generically expected if the background trajectory is curved as modes of interest leave the horizon. By parametrizing Eq. (3.26) with the scalar value cos ∆ ≡ √

PRS , PR PS

(3.27)

it was shown in Ref. [49] that, for the case of Nf = 2, the value of r is suppressed relative to the single-field, slow-roll expectation by r ≈ 16 sin2 ∆, to first-order in slow-roll. In principle, ∆ may be detectable from CMB observations [50, 51]. However, by differentiating Eq. (3.20) with respect to time t, the comoving curvature perturbation will not necessarily be constant even for k  aH. Instead, H R˙ = − 2 δPnad , φ˙ 0

(3.28)

where δPnad is the non-adiabatic pressure perturbation [52–54]. This quantity is the difference between the total pressure perturbation δP = φ˙ I δφ˙ I − φ˙ I φ˙ I Φ − V,I δφI ,

(3.29)

and the adiabatic pressure perturbation δPad = c2s δρ, where the speed of sound is c2s = P˙ /ρ˙ and the lapse function is 1 ˙ φI δφI , (3.30) Φ= 2H defined in the spatially-flat gauge [55]. Given the total density perturbation δρ = φ˙ I δφ˙ I − φ˙ I φ˙ I Φ + V,I δφI ,

(3.31)

the non-adiabatic pressure power spectrum PδP,nad reduces to PδP,nad (k) =

k3 2π 2 a2

 I J L ∗ 0 A A ψI ψLJ + AI B J ψI L ψLJ

(3.32) 

0 + B I B J ψ 0 L ψ ∗0 , +B I AJ ψJ∗ L ψLI I LJ

where (0) indicates a derivative with respect to e-foldings Ne and we have defined the vectors     1 1 0 0 0,L 2 0 2 0,M φ −3H φL − ∂L V ∂I V + H ∂M V φ δLI + φL φI (3.33) AI = 3aH 2  2 and

 BI = 1 − c2s H 2 φ0I .

(3.34)

By analogy to Eq. (3.24), we can build an entropy perturbation from the non-adiabatic pressure perturbation [26, 47, 56], with δS =

H δPnad . P˙

– 10 –

(3.35)

From this we obtain our final definition of isocurvature, the comoving entropy spectrum, given by  2 H Pent (k) = PδP,nad . (3.36) P˙ 3.5

δN formalism

The separate-universe assumption [53, 57–61], often referred to as δN , states that when smoothed on some physical scale much larger than the horizon, the evolution of each smoothed patch can be computed using only background quantities. By identifying that ζ = δN , where ζ is the curvature perturbation on constant density hypersurfaces and δN measures the variation in the number of e-folds between an initial flat hypersurface and a subsequent constant density hypersurface, Lyth and Rodriguez demonstrated that this assumption can be taken advantage of when computing correlation functions by performing a Taylor expansion in terms of the initial conditions [62]. 1 ζ = N,I δφI∗ + N,IJ δφI∗ δφJ∗ + . . . . 2

(3.37)

The main difficulty in this approach lies in computing the derivatives of the number of efolds (N,I ≡ ∂Ne /∂φI,∗ , N,IJ etc.). However for sum-separable models these expressions can be computed analytically [63, 64]. For models with fields much lighter than H at horizon crossing, the numerically intensive calculation of solving for the modes may therefore be unnecessary. MultiModeCode implements this δN slow-roll formalism where we assume that t∗ is the moment when the pivot-scale k∗ leaves the horizon and that the field perturbations at this time are uncorrelated, with a power spectrum IJ Pδφ

=



H 2π

2

δ IJ .

(3.38)

We also assume that the tensor modes, which are massless and uncoupled to the matter sector, have a power spectrum Ph = 8 (H/2π)2 . At least to first order, on super-horizon scales ζ = R [65], which allows us to compare the predicted power spectrum for ζ using the δN formalism to the adiabatic power spectrum in Eq. (3.23). If the potential V is sum-separable so that V =

X

VI (φI ),

(3.39)

I

then we can use the Klein–Gordon equations (3.5) for the scalar fields to obtain a sumseparable expression for the amount of expansion between the two surfaces X Z c VI Ne = − (3.40) 0 dφI , ∗ VI I

where VI0 ≡ dVI /dφI . If V were not sum-separable, the derivatives of Ne would in general have to be obtained numerically by evolving the background equations of motion (3.5) on a stencil in field-space and taking the finite difference. We have not implemented this feature in MultiModeCode as it is more computationally intensive than solving the mode equations.

– 11 –

When the potential is sum-separable, the derivatives of Ne can be simplified into the expressions [63, 64] 1 V ∗ + Zc N,I = p ∗ I ∗ I (3.41) V 2I and

N,IJ = δIJ



η∗ 1 − I∗ 2I



VI∗ + ZIc V∗

where ZIc = V c c ZIJ =

V2 − c∗ V

r



2  J

Nf X

K

K=1



I





cI − VIc , c

− δIK

X

and η≡

I =

I

(3.42)

(3.43)

   ηK  J , − δJK 1 −  

 

and the slow-roll parameters are

≡

∂ZJc 1 +p ∗ ∗ , 2J V ∂φ∗I

(3.44)

c

1 X VI02 2 V2

(3.45)

X V 00

(3.46)

I

X

ηI =

I

I

V

I

.

The contribution from the EOI surface is therefore completely encoded in the functions ZI and its derivatives. The relationship (3.37) and the expansion equation (3.40) allow us to define pivot-scale observables for the scalar perturbations ζ. We will focus on the observables obtainable only through the first and second derivatives of Ne , and express our results only to the lowest order in slow-roll. We start with the ζ power spectrum Pζ = N,I N ,I and the tensor-to-scalar ratio r=



H 2π

2

,

(3.47)

8 , N,I N ,I

(3.48)

which have simple expressions only in terms of N,I . The adiabatic and tensor spectral indices ns and nt also have easily evaluated expressions 2 ns − 1 = −2∗ − + N,I N ,I and nt =



−2∗ . 1 − ∗

2 V



V,IJ N ,I N ,J N,K N ,K

(3.49)

(3.50)

The expression for the scalar running αs is more complicated, but straightforward to compute (e.g., Eq. 6.14 in Ref. [66]).

– 12 –

To obtain the amplitude of the predicted non-Gaussianity we further assume that the field perturbations at horizon crossing are purely Gaussian, since the non-Gaussianity generated by sub-horizon evolution of the modes is typically slow-roll suppressed [63, 67], assuming that slow-roll is not violated. Following Refs. [63, 68], we use the non-linearity parameter Q 3 Bζ N,I N,J N ,IJ k 6 − fNL ≡ Pi i3 ≈ , 5 4π 4 Pζ2 (N,K N ,K )2 i ki

(3.51)

where Bζ is the bispectrum. Given Gaussian field perturbations at horizon crossing, the trispectrum amplitude is then parametrized by the non-linearity parameters [69, 70] τNL = and gNL =



N,IJ N ,IK N ,J N,K

25 54

(3.52)

(N,L N ,L )3 

N,IJK N ,I N ,J N ,K (N,L N ,L )3

.

(3.53)

Since gNL ∼ N,IJK we do not compute it here, although it could be implemented by taking the third derivative of Ne as in Ref. [64]. 3.6

Bundle width

An alternative method of monitoring isocurvature is to acknowledge that under slow-roll, the separate universe assumption is precisely analogous to geometrical optics in field space [34]. The smoothed spatial patches described in Sect. 3.5 each correspond to a distinct noninteracting trajectory in field space with perturbed initial conditions with respect to some arbitrary fiducial trajectory. These perturbed trajectories can√then be thought of as forming a bundle moving through a medium with refractive index 2. One can therefore track isocurvature evolution using only background quantities, by associating isocurvature growth and decay with dilation and contraction of the bundle. While the precise analogy with geometrical optics does not remain when slow-roll is violated, one still has a useful set of geometrical quantities for understanding the evolution of field perturbations. Under slow-roll, the Klein–Gordon equation reduces to dφI = −∂I log V, dNe

(3.54)

which is Huygen’s equation and an infinitesimal vector propagated along the beam is called a Jacobi field. If we take δφI to be such a field, we can obtain from Eq. (3.54) how it will propagate: dδφI = − [∂I ∂J log V ] δφJ , (3.55) dNe which is the slow-roll analogue of Eq. (3.12) [34, 42]. Indeed, we could have recast the whole of Sect. 3.3 in this language [42]. The term in square brackets is usually referred to as the expansion tensor and it encodes all information required for tracking field perturbations;5 5

This point is heavily emphasised in the context of the transport method of computing inflationary correlation functions [34, 40–42, 71].

– 13 –

under slow-roll it has a particularly simple geometric interpretation. We can decompose the expansion tensor as θ ∂I ∂J log V = + σIJ + ωIJ , (3.56) Nf where σIJ is the symmetric shear, ωIJ is the antisymmetric twist, and the key quantity for our purposes is the dilation, given by the trace θ = Tr ∂I ∂J log V.

(3.57)

If θ > 0, then isocurvature is growing and if θ < 0, then isocurvature is decaying. We can then find a measure Θ of the bundle width by integrating this along the inflationary trajectory Θ ≡ exp

Z

N

0

θ(N )dN

N0

0



,

(3.58)

which normalizes the bundle width so that Θ(N0 ) ≡ 1. In situations where we only want to solve the background equations of motion, the bundle width is informative for understanding whether or not ζ becomes conserved on superhorizon scales, which is a crucial requirement when comparing the predictions of a model with observation. For two fields Θ → 0 is a necessary and sufficient condition for the approach to an adiabatic limit. However when there are more fields the situation is more complicated, e.g., the bundle may contract to a sheet rather than a caustic. We refer the reader to Ref. [34, 66, 71] for more details.

4

The method

We outline the procedure used to obtain the power spectrum predictions, with the algorithmic structure of MultiModeCode in Algorithm 1. While this largely follows previous implementations, such as Pyflation [25–28], we give the method the sake of clarity and reproducibility. We start by defining the functional form of the potential V and prior probability distribution functions (PDFs) for the parameters that define V , which we call Lagrangian parameters or model parameters, and the background initial conditions φI,0 and φ0I,0 . We treat the simple situation of exactly specifying a set of Lagrangian parameters and initial conditions as a special case, where the prior probability for these parameters is trivial. Given these priors, the program will build a numerical sample by iteration until a pre-defined number of samples is reached. MultiModeCode first solves the background equations of motion (3.5) until the endof-inflation. While we have included the natural condition of  = 1 as the default ending criterion for inflation, there is complete functionality to end inflation by another method, in particular a waterfall transition via the hybrid mechanism [72, 73] at some reference phasespace point. Given a value for the number of e-folds N∗ between when the pivot scale k∗ leaves the horizon and the end-of-inflation, which is either fixed by the user or set in each iteration of the code through the sampling of a prior probability P (N∗ ), we obtain the value of H at horizon crossing by interpolating the numerical background solution. The pivot scale k∗ must be pre-defined by the user and defaults to 0.002 Mpc−1 . From this, we normalize the size of the universe so that k∗ = aH∗ at Ne = Ntot − N∗ .

– 14 –

Algorithm 1 MultiModeCode method define sample size, V , k∗ for all elements in sample do procedure Background Solver: get Lagrangian parameters for V and ICs for Eq. (3.5) from prior PDF with the end-of-inflation (EOI) criterion set by user, solve Eq. (3.5) until EOI check inflation (¨ a > 0) started and ended procedure Scale-factor Normalizer: get N∗ from user or by prior PDF check total inflationary e-folds Ntot ≥ N∗ define a such that k∗ = a∗ H∗ at Ne = Ntot − N∗ before inflation ends procedure δN Calculator: if V is sum-separable, then calculate δN observables near k∗ for all modes k do procedure Mode Initializer: define initial time Ne,0 with k  a0 H0 while the corrections to Eq. (3.15) are above some tolerance: set earlier Ne,0 and check Ne,0 > 0 set Bunch-Davies IC for mode matrix ψIJ (k) at Ne,0 procedure Mode Solver: solve Eq. (3.12) until k ≈ aH change variable as in Eq. (3.18) and solve until EOI calculate power spectra for k procedure k∗ -observable Calculator: calculate amplitudes, spectral indices, etc. at k∗ by finite difference in k-space

For each scale of interest k, we set the modes’ initial conditions at a time Ne,0 when that mode is significantly sub-horizon, k  a0 H0 . For the Bunch-Davis initial state, this point is chosen iteratively by making sure that the relative corrections to Eq. (3.15) that are sub-dominant for k  aH are smaller than a pre-defined tolerance. This tolerance is set to 1 × 10−5 ; from observing the sub-horizon evolution of the modes, using a tolerance at least this tight gives no change to the value of the modes at horizon crossing. We then solve the mode equations (3.12) for the variable ψIJ for the period of time when the modes are smaller than the causal horizon, k & aH, and then switch to a two-index matrix built from the uI in Eq. (3.18) for super-horizon evolution. The power spectra are calculated for each k and various pivot-scale statistics are evaluated by finite-difference between a few scales ki near k∗ . If the potential V is sum-separable, the program also calculates the δN values for the observables described in Section 3.5. Numerous checks are performed on the background and mode equation evolution so that MultiModeCode will either fail gracefully if a fatal exception is raised or declare a particular initial parameter set invalid and iteratively generate a new set of parameters in order to explore cosmologically relevant parameter sets. We have extensively tested the numerical stability of the code and have included a number of easily controllable options allowing the user to control the numerical accuracy, as well as the type of ODE solver.

– 15 –

10−7

10−28

10−13

10−32

10−19

P

10−36

10−25

10−40 10−44

10−31

|δPnad |2 |δPad |2

10−48 0

10

Pent PR PS

10−37

20

30

40

50

0

Ne

10

20

30

40

50

Ne

Figure 1. The evolution of the power spectra during the last 55 e-folds of inflation for a twofield Nf -quadratic model. (Left) The power spectrum for adiabatic (green) and non-adiabatic (blue) pressure perturbations δP . The total pressure spectrum and the adiabatic pressure spectrum are nearly coincident on this scale, so the total pressure spectrum PδP has not been plotted. The gray area is an estimate of the region dominated by double-precision error due to round-off in PδP,nad . (Right) The power spectra for perturbations in the adiabatic curvature PR , the isocurvature PS , and the comoving entropy Pent . Pent is a rescaling of PδP,nad and is numerically unstable for Ne & 30 in this model. PS is numerically stable until the end of inflation.

5 5.1

Numerical results Isocurvature stability

Fig. 1 illustrates a problem that arises when computing the isocurvature spectra PδP,nad and Pent . We have plotted the super-horizon evolution of the power spectra for the adiabatic and non-adiabatic pressure perturbations, as well as the adiabatic curvature, entropic, and isocurvature spectra, with N∗ = 55, for a two-field inflation model with the potential 1 1 V = m21 φ21 + m22 φ22 . 2 2

(5.1)

To match the analysis performed in Refs. [26, 28, 74, 75] we choose m21 = 10−11.7 , m22 = 10−10.0 , and initial conditions φ1,0 = φ2,0 = 12.0 MPl . In particular, Fig. 1 can be compared directly to Figs 1 and 3 in Ref. [26]. With this choice of parameters, the background trajectory evolves primarily along the direction of the heavier field φ2 for Ne . 25, then turns sharply toward the φ1 direction for the remainder of inflation. The effect of this turn on the superhorizon perturbations can be seen clearly in the power spectra in Fig. 1. In general, the calculation of PδP,nad and Pent becomes dominated by numerical error as the isocurvature perturbations decay. From Fig. 1, regardless of the amplitude of the isocurvature modes, the adiabatic pressure perturbations δPad = c2s δρ do not exponentially decay between horizon exit and the end of inflation. For the example model (5.1), the power spectrum for δPad is approximately constant after the turn at Ne ≈ 25. However, the total pressure perturbation δP is approximately equal to δPad during this time and the difference between the two reduces exponentially as the isocurvature modes decay.

– 16 –

Since δPnad ≡ δP − δPad and δPad → δP , the numerical accuracy for δPnad is limited by the real precision of the computer, which results in a finite difference error in the numerical calculation of δPnad and a loss of significance. Using standard double precision accuracy, the expected error in δPnad should then be ∆err PδP,nad ∼ O(10−15 ) PδP ∼ O(10−15 ) PδP,ad ,

(5.2)

which is confirmed in Fig. 1. Without correcting for this dominant error term, the value of PδP,nad will oscillate arbitrarily between zero and the limit in Eq. (5.2), which is an upper bound on the amplitude of the non-adiabatic pressure perturbations. Since entropic perturbations are usually defined as [56] SIJ ≡ ζI − ζJ ,

(5.3)

where ζI is the curvature perturbation resulting from the I th fluid, this problem will arise naturally for all calculations of Pent . In contrast, the calculation of PS in Eq. (3.25) is directly proportional to the value of the decaying isocurvature modes in the kinematic basis. Using this isocurvature spectrum largely alleviates the numerical problems with δPnad , yielding a more faithful measure with a higher degree of accuracy. Figure 1 shows the exponential decay of PS after the super-horizon turn at Ne ∼ 25. We compare this to Pent , which becomes numerically unstable at Ne ≈ 30, showing that the two measures PS and Pent are separated by 27 orders of magnitude at the end of inflation, despite being of the same magnitude at horizon crossing.6 5.2

A case study: Nf -flation with a step

We have shown in Refs. [30, 31] that MultiModeCode is able to produce large volume Monte Carlo samples for Nf –monomial inflation with the potential V =

1X λI |φI |p , p

(5.4)

I

for real exponents p [76–85]. In Ref. [30] we focused on the Nf –quadratic case with p = 2 and demonstrated that the predictions for the power spectrum do not sensitively depend on the prior probability chosen for the initial conditions of the fields. In Ref. [31] we further demonstrated this for the general case in Eq. (5.4), while focusing on the gravitational wave consistency relation. We were able to straightforwardly compare the analytical δN results to the numerics, greatly simplifying the procedure for comparing analytical results to the full numerical calculation. We include all of the IC priors used in these papers in MultiModeCode. Since we have already demonstrated the power of MultiModeCode in Monte Carlo sampling, in this paper we will instead focus on a few case studies that are interesting 6

As the adiabatic limit is approached, PS can also receive a dominant contribution from roundoff error in the Gram-Schmidt orthogonalization procedure. If some components of the isocurvature vectors sIK are much smaller than others, this can result in a spurious projection of PR onto the isocurvature directions. In MultiModeCode we have included an optional subroutine renormalize remove smallest in modpk potential.f90, where the components of sIK are set to zero if they do not affect the normalization of sK , i.e., if the value of sIK is indistinguishable from roundoff error. In practice, we have never seen this problem arise, so this option needs to be uncommented in the source code before compilation.

– 17 –

due to their analytic intractability. We present results for a multifield generalization of the inflationary step-potential first used in Ref. [10]. This potential has the form   ¯ I  1X 2 2 φI − Φ V = (5.5) mI φI 1 + cI tanh 2 dI I

¯ I specifying the slope, amplitude, and with masses mI and real constants dI , cI , and Φ th position, respectively, for a step feature in the I field. Phase transitions in sectors coupled only gravitationally to the inflaton sector may generate these hyperbolic-tangent features in V and leave an observable imprint in the primordial density spectra if these symmetry breaking transitions occur during the last O(60) e-folds of inflation [10, 32]. In the sharp-step limit, these features introduce oscillations as a function of k into the adiabatic curvature power spectrum and a scale-dependent, oscillatory bispectrum [10, 33, 86, 87]. To keep V > 0 we require cI < 1 and to satisfy the latest constraints on oscillations in the scalar power spectrum amplitude requires cI . 10−3 , assuming that the step occurs as the scales relevant for the CMB leave the horizon [88–90]. With cI → 0, Eq. (5.5) is an uncoupled assisted inflation model [76, 91], first proposed in Ref. [81]. Models with a step feature are additionally interesting, because they can fit a wider range of data and have been well-studied in the single-field case. In particular, Ref. [33] contains an elegant analytical calculation for the single-field case of Eq. (5.5). However, replicating the same calculation for the general potential would be difficult — if not impossible — with the same techniques, since the possible existence of isocurvature perturbations significantly complicates the analysis. Consequently, a numerical exploration of this model is well-motivated. Fixing the number of fields to Nf = 10, we set the initial conditions to φI,0 = 10, with the initial velocities set in slow-roll. The size and slope of the step are set to cI = 10−3 and dI = 10−2 respectively, and the masses mI relative to the fiducial mass to m ¯ 2 = 4.30 × 10−11 , which in the single-field limit yields As at the best-fit value from the Planck T T data. Following Ref. [80], we choose the masses mI according to the Marˇcenko-Pastur distribution P (m2I ) where

=

1 2πm2I m ¯ 2β

q

β+ − m2I

 p 2 β± = m ¯2 1± β



 m2I − β− ,

(5.6)

(5.7)

with β = 1/2. This distribution of masses is derived in Ref. [80], and has also been used in Refs. [30, 84, 92, 93]. ¯ I for each field at the field-space point where the pivot scale We set the step positions Φ −1 k∗ = 0.002 Mpc leaves the horizon at N∗ = 55 e-folds before the end of inflation in the ¯ I are functions no-step limit, cI → 0. Since the fields have identical initial conditions, the Φ only of the masses, so we plot the step positions versus the mI in Fig. 2. We also present the field-space trajectories according to Eq. (3.5) for the last 75 e-folds of inflation with these parameters. The heavier fields relax more quickly toward their minimum at φI = 0 and the lighter fields have a larger value at horizon crossing. Since cI = 10−3 , the step is not obviously visible at the level of the background trajectory without zooming in significantly. However, Fig. 3 shows the substantial effect on the power spectra due to the steps. We see oscillatory behavior in the adiabatic, isocurvature, and entropic power spectra, but

– 18 –

10

8

8

6

6

φI

¯I Φ

10

4

4

2

2

0

0 0.0

0.5

1.0

1.5

2.0

2.5

3.0

m2I /m ¯2

−20 −10

0

10

20

30

40

50

60

Ne

Figure 2. (Left) The masses mI for each of the 10 fields in Eq. (5.5), drawn from the distribu¯ I for that field, tion (5.6) with m ¯ 2 = 4.3 × 10−11 , compared to the corresponding step positions Φ ¯ I , given the initial which is positioned so that the pivot scale k∗ = 0.002 Mpc−1 leaves the horizon at Φ conditions φI,0 = 10. (Right) The field trajectories (colored lines), with the same initial condition, as a function of e-folding Ne , with k∗ (vertical line) leaving the horizon 55 e-folds before the end of ¯ I are marked in blue and Ne has been renormalized so that k∗ = aH inflation. The step positions Φ at Ne = 0.

almost no change in the tensor spectrum. Furthermore, we can see clearly that PS and Pent exhibit a nearly identical feature, simply scaled by a factor of roughly 65. These features in the isocurvature spectrum may lead to interesting effects during reheating or the subsequent evolution of the post-inflation universe.

6

Conclusion

We present the Fortran 95/2000 code MultiModeCode, designed to maximize computational efficiency when numerically exploring a broad range of multifield inflation models. The code also provides Monte Carlo sampling of prior probabilities for inflationary model parameters and initial conditions, enabling automated model exploration and the computation of probability distributions for observables. The mode equation method has a broad range of applicability, but the computational cost scales with the number of fields as O(Nf2 ). For models with sum-separable potentials, we have also implemented a slow-roll δN calculation, which only requires solving the background equations of motion once in order to obtain the full power spectrum as well as higher order statistics. This drastically improves computation time, since the background equations of motion are only O(Nf ). This code was used to explore the predictions of models with O(100) fields in Refs. [30, 31]; here, we demonstrated its use with an Nf -flation model with a step. We find that a feature in the inflationary potential not only results in a feature in both the adiabatic power spectrum as a function of scale, PR (k), as well as the isocurvature spectra PS , Pent , and PδP,nad , with possible implications for the dynamics of many-field preheating scenarios. Further, we see numerical evidence that the isocurvature spectrum PS is a simple rescaling of the entropic spectrum Pent , indicating that the projection of the mode power spectrum

– 19 –

×10−10

1.4

×10−17

8 1.2 7

PS

PR

1.0

6

0.8

5

cI → 0

1.2

0.6

cI = 10−3

4

0.4

×10−10

×10−19

1.8 1.1

Ph

Pent

1.4

1.0

1.0

0.9 0.8 10−4

10−3



k Mpc

10−2  −1

10−1

0.6 10−4

10−3



k Mpc

10−2  −1

10−1

Figure 3. Features in the power spectra due to the step (5.5), which is positioned so that it affects the power spectra around the pivot scale k∗ = 0.002 Mpc−1 (gray). We compare (dashed, blue) the no-step case with cI = 0, to (solid, green) the case with cI = 10−3 . While there are oscillations in the adiabatic PR , isocurvature PS , and entropic Pent spectra, there is little variation in the tensor spectrum Ph .

onto the isocurvature directions is related to a quantity that sources a change in R on superhorizon scales. MultiModeCode complements codes that currently exist to numerically compute the inflationary power spectra [10, 12, 19–21, 25–28, 71, 94]. We provide a basic usage manual for MultiModeCode in Appendix B to help users to adapt this program to their own problems. The theoretical basis of the method is outlined in Section 4. The ability of MultiModeCode to solve numerically challenging problems, such as the step-potential in §5.2, and to provide large samples of many-field inflationary models adds significantly to the early universe cosmologist’s toolkit for exploring and understanding realistic inflation models.

Acknowledgments We thank Grigor Aslanyan, Adam Christopherson, Mafalda Dias, Ian Huston, Karim Malik, David Mulryne, David Seery, and Jonathan White for many helpful discussions throughout the duration of this work. We acknowledge the use of the dvode f90 m.f90 numerical integration code,7 developed by G. Byrne, S. Thompson, and LLNL, which we redistribute here under the terms of their BSD-like license, and the CSV writing capabilities of the FLIBS 7

www.radford.edu/~thompson/vodef90web/

– 20 –

open-source library,8 developed by Arjen Markus. JF is supported by IKERBASQUE, the Basque Foundation for Science. HVP is supported by STFC and the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007-2013) / ERC grant agreement no 306478-CosmicDawn. The authors acknowledge the contribution of the NeSI high-performance computing facilities and the staff at the Centre for eResearch at the University of Auckland. New Zealand’s national facilities are provided by the New Zealand eScience Infrastructure (NeSI) and funded jointly by NeSI’s collaborator institutions and through the Ministry of Business, Innovation & Employment’s Research Infrastructure programme [http://www.nesi.org.nz]. This work has been facilitated by the Royal Society under their International Exchanges Scheme. This work was supported in part by National Science Foundation Grant No. PHYS-1066293 and the hospitality of the Aspen Center for Physics.

A

Non-canonical kinetic terms

This appendix largely follows Ref. [59] and describes the background and first-order mode equations for inflation models with multiple scalar fields, but a general field space metric GIJ (φK ). These equations are not coded into MultiModeCode, but are an important reference, since multifield models with non-canonical kinetic energy terms fit the Planck data extremely well [95]. These have also been implemented into the single-field version of ModeCode by Ref. [96].9 We recommend Ref. [59] for a thorough derivation of these quantities. We start with the action for Nf scalar fields φI , given by   Z 1 4 √ I µ J I S = d x −g − GIJ ∂µ φ ∂ φ − V (φ ) . (A.1) 2 Again, we assume that Greek indices α ∈ [0, . . . , 3] are for spacetime, upper-case Latin letters index the number of fields, I ∈ [1, . . . , Nf ], and lower-case Latin letters index three-space, i ∈ [1, . . . , 3]. The action (3.1) is a special case of Eq. (A.1), with GIJ = δIJ . To change the equations in §3 to reflect the general field-space metric, we follow the typical procedure of replacing partial derivatives with respect to the fields with covariant derivatives. Varying Eq. (A.1) with respect to the fields gives the background equation of motion DφI + 3H φ˙ I + GIJ V;J = 0, dt

(A.2)

where we have assumed that the background fields are homogeneous, ∂i φJ = 0. In Eq. (A.2) we have used the covariant differential [59] DφI = dφI + ΓIJK φ˙J dφK ,

(A.3)

with the field-space Christoffel symbols 1 ΓIJK = GIL (GLJ,K + GLK,J − GJK,L ) , 2 8

(A.4)

http://flibs.sourceforge.net/ We will update MultiModeCode to give this functionality in the near future and a simple implementation in Mathematica using the Transport method will be available soon [71]. 9

– 21 –

which gives the covariant derivative V;I ≡ DV /dφI = ∂V /∂φI , because V is a scalar. Following the treatment in Section 3, it is convenient to express the equations of motion for the field perturbations δφIk in spatially-flat gauge in terms of the generalized trajectory velocity φ˙ 0 , which can be extended from Eq. (3.4) to φ˙ 20 ≡ GIJ φ˙ I φ˙ J .

(A.5)

From φ˙ 0 , we define the adiabatic unit vector in the kinematic basis [47, 48] as ω I ≡ φ˙ I /φ˙ 0 , which projects vectors along the adiabatic direction in the generalized field-space. Finally, the mode equations in spatially-flat gauge read [59] D2 I D k2 δφk + 3H δφIk + 2 δφIk + C I J δφJk = 0, 2 dt dt a

(A.6)

which is a simple generalization of Eq. (3.10), with the mass matrix  H ω I V;J + ωJ V ;I + 2(3 − )H 2 ω I ωJ , C I J = GIK V;K;J − φ˙ 20 RI JKL ω K ω L + 2 ˙ φ0

(A.7)

where RI JKL is the field-space Riemann tensor, built from GIJ . Again, we note that index contraction implies summation with respect to the field-space metric, XI Y I = GIJ X I Y J . Quantizing the modes using the mode-matrix δφI = ψI J aJ proceeds as in Section 3.3 and the adiabatic power spectrum PR is identical to Eq. (3.23), except contracting with respect to GIJ . The isocurvature directions should again be found by Gram-Schmidt orthogonalization, except now implemented in the curved field-space, and the spectra PδP,nad and Pent can be found as in Eqs. (3.33) and (3.36). The projected isocurvature spectrum PS in Eq. (3.25) can instead be built by replacing the summation over the isocurvature directions ˆ IJ , after by contraction with respect to the isocurvature directions of the field-space metric, G transforming to the kinematic basis.

B

MultiModeCode usage

MultiModeCode is publicly available at www.modecode.org and is released with a Modified BSD License. MultiModeCode has been implemented in a mix of Fortran 95/2000 and has been thoroughly checked on Mac OS and Linux systems with the freely-available GNU Fortran compiler (version 4.6.3+) and with Intel Fortran (version 14.0.2), with the Intel compiler yielding significant improvements in speed. There are no external dependencies necessary to use MultiModeCode. We include a driver file multimodecode driver.f90 that contains the basic structure needed to explore the predictions of a model. The driver has many important routines for calculating the power spectrum and outputting the results. The file parameters multimodecode.txt contains runtime parameters that are often changed between subsequent runs. The parameters are listed in Fortran namelists, so a change here does not require the user to recompile the whole code. There are a mix of basic and advanced parameters available and we describe them here. The &init namelist holds important parameters related to initializing the program, the choice of inflationary potential V , and program output:

– 22 –

&init num_inflaton = 10 potential_choice = 1 vparam_rows = 4 slowroll_infl_end = .true. instreheat = .false. /

The number of fields is set with the variable num inflaton and an indentifying number is chosen for the potential with potential choice. The currently available potentials are listed in the routine pot(phi) in modpk potential.f90, including the multifield Nf -quadratic [76, 81, 97], Nf -monomial inflation with V ∼ λI |φI |n [31], two-field hybrid inflation [72, 73], a product of exponentials [76], and the multifield generalization of the hyperbolic-tangent step potential in Ref. [10], which was used in Section 5.2. Adding a new potential is as easy as providing the potential V and its derivatives in modpk potential.f90 with a new value for potential choice. The array vparams contains information passed to the potential function (e.g., masses and couplings) and has dimensions (vparam rows)×(num inflaton). The values for vparams are set in parameters multimode.txt in the namelist ¶ms; variables related to the pivot scale N pivot and k pivot are also set here. The user can change the conditions for inflation to end by varying slowroll infl end, which when set to .true. evolves the background fields until  = 1. If you do not want this to be the ending criterion, then set slowroll infl end=.false. and adapt the subroutine alternate infl end in modpk odeint.f90 to change the ending condition. Furthermore, with instreheat=.true. N∗ becomes a derived parameter by requiring inflation to thermalize immediately after the end of inflation with w = 1/3, as in Ref. [20]. We use the &analytical namelist to set the options for calculating the power spectrum, either using the δN calculations of §3.5 or evaluating the full mode equations of §3.3 or both. &analytical use_deltaN_SR = .true. use_horiz_cross_approx = .false. evaluate_modes = .true. get_runningofrunning = .false. /

With use deltaN SR=.true. MultiModeCode will calculate the δN observables at the pivot scale (as given in namelist ¶ms) and if use horiz cross approx=.true., it will ignore the contribution to N,I and N,IJ from the end-of-inflation surface via the horizon crossing approximation (HCA) [63, 83]. Setting evaluate modes=.false. will make the program only solve the background equations of motion, relying on the δN calculations to obtain the spectra. If get runningofrunning=.true., then the derivative of αs with respect to ln k will be calculated by a five-point stencil finite difference method, which requires two additional calls to the code that solves the mode equations, significantly slowing down the speed of the program. The way by which the initial conditions are chosen for a given simulation depends on the namelist &ic sampling nml namelist. &ic_sampling_nml ic_sampling = 1

– 23 –

numb_samples = 1 energy_scale = .1 save_iso_N = .false. N_iso_ref = 55 /

The variable ic sampling controls the main behavior of the initial conditions’ prior probability and is set to an identifying number defined in the file modpk icsampling.f90 with the ic samp flags type. The currently available values are below: type :: ic_samp_flags integer :: reg_samp = 1 integer :: eqen_samp = 2 integer :: slowroll_samp = 3 integer :: iso_N = 6 end type

Invoke each case by setting ic sampling equal to the desired number. The functionality of each of these cases is: • reg samp (regular sampling): the simple case of setting a multifield initial condition and calculating the power spectrum. The initial conditions for the fields are set as the variable phi init0 in the ¶ms namelist and the fields’ velocities are assumed to be initially in slow-roll. • eqen samp (equal-energy sampling): quasi–equal-area sampling of a phase-space surface with same initial energy, as in Refs [8, 9, 30]. Set the initial energy with the 2 = (8πG)−1 = 1. To also record energy scale variable in units of MPl , where MPl and output the background field values as the background reaches Ne = N iso ref set save iso N=.true.. The prior ranges for the fields and velocities are set in the namelist &priors. • slowroll samp (slow-roll sampling): choose initial conditions in field space and set velocities by the slow-roll condition. The prior ranges for the fields are again chosen in &priors. P 2 • iso N (sampling e-fold surface): uniformly samples the surface Ne = I φI /2p for Nf -monomial inflation, as in Refs [29, 30]. Set the value for Ne with N iso ref in the &ic sampling namelist. To implement a different initial conditions measure or sampling behavior, add a new identifier for the ic samp flags type. The initial conditions are set in the routine get ic in modpk icsampling.f90 and you will need to implement your sampling technique here, following the examples in the code. Similarly, the prior probabilities on the vparams array that defines the potential are set through the ¶m sampling nml namelist. ¶m_sampling_nml param_sampling = 1 use_first_priorval = .true. vp_prior_min(1,:) = -14 vp_prior_max(1,:) = -12

– 24 –

varying_N_pivot = .false. /

As with the initial conditions sampler, different behaviors are chosen by setting the variable param sampling to a unique integer as specified in the param samp flags type in modpk icsampling.f90. type :: param_samp_flags integer :: reg_constant = 1 integer :: unif_prior = 2 integer :: log_prior = 3 end type

Again, invoke each case by setting param sampling to the desired number. To vary the number of e-folds after the pivot scale leaves the horizon, set varying N pivot=.true. and set the limits on the prior on N∗ in the &priors namelist, where we have assumed a uniform prior. Note that this is overridden if instreheat=.true.. To use a different prior on the vparams array, add a new integer to the param samp flags type and change the routine get vparams in modpk icsampling.f90. The default behavior is: 1. reg constant (regular, constant parameters): the vparams array is kept constant, as specified in the ¶ms namelist. 2. unif prior (uniform prior probability): each column in vparams is chosen with a uniform prior between vp prior min and vp prior max. If use first priorval=.true., then the first entry in the priors are used for all the columns. 3. log prior (logarithmic prior probability): the columns of a dummy array αIJ are chosen with a uniform prior according to the priors in the namelist. The columns in the vparams are then set by vparams =10αIJ . The ¶ms namelist is used to set the vparams array, the fields’ initial conditions phi init0, the pivot scale k pivot, and the number of e-folds between horizon exit and the end of inflation for the pivot scale by N pivot. ¶ms N_pivot = 55.0 k_pivot = 0.002 dlnk = 0.4 phi_init0 = 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 vparams(1,:) = -11.4 -11.0 -10.8 -10.6 -10.5 -10.4 -10.3 -10.2 -10.1 -9.9 vparams(2,:) = 8.8 7.3 6.0 4.8 3.6 2.6 1.7 1.0 5.6 0.16 vparams(3,:) = 1e-3 1e-3 1e-3 1e-3 1e-3 1e-3 1e-3 1e-3 1e-3 1e-3 vparams(4,:) = 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 /

These are all default values and may be overridden with different choices of sampling, as mentioned above. The variable dlnk is the difference in k-space that is used when we evaluate the pivot-scale observables from the mode equations via finite difference. All scales are in units of Mpc−1 and fields are in MPl .

– 25 –

The prior probability ranges for the sampling of the fields’ initial values and N∗ are in the &priors namelist, which is relatively self-explanatory. &priors phi0_priors_min = 2.0 2.0 phi0_priors_max = 30.0 30.0 dphi0_priors_min = -1.262e0 -1.262e0 dphi0_priors_max = 1.262e0 1.262e0 N_pivot_prior_min = 30 N_pivot_prior_max = 70 /

MultiModeCode defaults to obtaining only the pivot-scale observables by taking the numerical derivative of the power spectra at k∗ . However, the full power spectra can also be solved for and provide the full description of the model’s predictions over the scales of interest, as in §5.2. We control this through the &full pk namelist. &full_pk calc_full_pk = .false. steps = 300 kmin = 1.0e-4 kmax = 1.0e0 /

The variable calc full pk options the calculation of the full P(k). The program will interpolate between a number of points, given by the variable steps, between the scales kmin and kmax. MultiModeCode is able to save and output a significant amount of data for later analysis. However, since this has an obvious affect on the speed of the code, the amount and verbosity of this output can be specified by the attributes of the out opt instance of the print options type in the print out namelist. &print_out out_opt%modpkoutput = .true. out_opt%output_reduced = .true. out_opt%output_badic =.false. out_opt%save_traj = .true. out_opt%fields_horiz = .false. out_opt%fields_end_infl = .false. out_opt%spectra = .false. out_opt%modes = .false. /

If you want nothing to write to stdout (the terminal, usually), then set modpkoutput=.false.; if less output is requested, then set output reduce=.true.. If output badic=.false., then any set of parameters that do not lead to a successfully inflating universe are discarded and ignored; otherwise, they will be saved and output into the data files with dummy values for their spectra. The remaining attributes of out opt controls what cosmological quantities are printed to file in addition to the pivot scale observables: save traj records the background trajectory as a function of Ne ; fields horiz saves the field values as the pivot scale crossing the horizon; fields end infl saves the field values at the end of inflation; spectra will

– 26 –

record the superhorizon power spectra as a function of Ne ; and modes prints all the mode functions during the entire evolution. Setting modes=.true. consequently results in a lot of output. The output of MultiModeCode is saved in comma-delimited CSV files, with the first row corresponding to a header that names each column. Subsequent rows correspond to different samples of the same model with different parameters. To change the output, simply find the point where the header is written in the code, add a new column(s), and print out the desired quantity in the correct order. Finally, some technical options are controllable via the &technical namelist, by changing the attributes of the tech opt instance of the tech options type. In particular, the choice of numerical integration scheme can be changed by the tech opt%use dvode integrator flag. Setting this to .true. will invoke a backwards-difference formula method, which is suitable for stiff problems, while .false. will use a fourth-order Runge-Kutta integrator. Various accuracy settings for the integrators can also be controlled in this namelist, with the global behavior set by the variable accuracy setting =−1, 1, 2. Using 1 sets a minimal amount of accuracy, which we find is suitable for calculating the adiabatic power spectrum for simple models. Increasing this to 2, increases the accuracy and, in particular, increases the accuracy as the evolution moves out of slow-roll, which we find is necessary to obtain resolved isocurvature spectra, as in Fig. 3. Lastly, you can override our choices for the absolute and relative error tolerances by setting accuracy setting=-1 and manually changing the remaining attributes of tech opt in the namelist, which is self-explanatory.

References [1] WMAP Collaboration, G. Hinshaw et al., Nine-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Cosmological Parameter Results, Astrophys.J.Suppl. 208 (2013) 19, [arXiv:1212.5226]. [2] Planck Collaboration, P. Ade et al., Planck 2013 results. I. Overview of products and scientific results, arXiv:1303.5062. [3] Planck Collaboration, P. Ade et al., Planck 2013 results. XVI. Cosmological parameters, arXiv:1303.5076. [4] Planck Collaboration, P. Ade et al., Planck 2013 results. XXII. Constraints on inflation, arXiv:1303.5082. [5] D. Baumann and L. McAllister, Inflation and String Theory, arXiv:1404.2601. [6] R. Easther and K.-i. Maeda, Chaotic dynamics and two field inflation, Class.Quant.Grav. 16 (1999) 1637–1652, [gr-qc/9711035]. [7] S. Clesse, C. Ringeval, and J. Rocher, Fractal initial conditions and natural parameter values in hybrid inflation, Phys.Rev. D80 (2009) 123534, [arXiv:0909.0402]. [8] R. Easther and L. C. Price, Initial conditions and sampling for multifield inflation, JCAP 1307 (2013) 027, [arXiv:1304.4244]. [9] R. Easther, L. C. Price, and J. Rasero, Inflating an Inhomogeneous Universe, arXiv:1406.2869. [10] J. A. Adams, B. Cresswell, and R. Easther, Inflationary perturbations from a potential with a step, Phys.Rev. D64 (2001) 123514, [astro-ph/0102236].

– 27 –

[11] WMAP Collaboration, H. Peiris et al., First year Wilkinson Microwave Anisotropy Probe (WMAP) observations: Implications for inflation, Astrophys.J.Suppl. 148 (2003) 213, [astro-ph/0302225]. [12] J. Martin and C. Ringeval, Inflation after WMAP3: Confronting the Slow-Roll and Exact Power Spectra to CMB Data, JCAP 0608 (2006) 009, [astro-ph/0605367]. [13] L. M. Hall and H. V. Peiris, Cosmological Constraints on Dissipative Models of Inflation, JCAP 0801 (2008) 027, [arXiv:0709.2912]. [14] R. Bean, X. Chen, H. Peiris, and J. Xu, Comparing Infrared Dirac-Born-Infeld Brane Inflation to Observations, Phys.Rev. D77 (2008) 023527, [arXiv:0710.1812]. [15] L. Lorenz, J. Martin, and C. Ringeval, Brane inflation and the WMAP data: A Bayesian analysis, JCAP 0804 (2008) 001, [arXiv:0709.3758]. [16] J. Martin and C. Ringeval, First CMB Constraints on the Inflationary Reheating Temperature, Phys.Rev. D82 (2010) 023511, [arXiv:1004.5525]. [17] J. Martin, C. Ringeval, and R. Trotta, Hunting Down the Best Model of Inflation with Bayesian Evidence, Phys.Rev. D83 (2011) 063524, [arXiv:1009.4157]. [18] J. Martin, C. Ringeval, R. Trotta, and V. Vennin, The Best Inflationary Models After Planck, JCAP 1403 (2014) 039, [arXiv:1312.3529]. [19] M. J. Mortonson, H. V. Peiris, and R. Easther, Bayesian Analysis of Inflation: Parameter Estimation for Single Field Models, Phys.Rev. D83 (2011) 043505, [arXiv:1007.4205]. [20] R. Easther and H. V. Peiris, Bayesian Analysis of Inflation II: Model Selection and Constraints on Reheating, Phys.Rev. D85 (2012) 103533, [arXiv:1112.0326]. [21] J. Norena, C. Wagner, L. Verde, H. V. Peiris, and R. Easther, Bayesian Analysis of Inflation III: Slow Roll Reconstruction Using Model Selection, Phys.Rev. D86 (2012) 023505, [arXiv:1202.0304]. [22] A. Lewis, A. Challinor, and A. Lasenby, Efficient computation of CMB anisotropies in closed FRW models, Astrophys.J. 538 (2000) 473–476, [astro-ph/9911177]. [23] A. Lewis and S. Bridle, Cosmological parameters from CMB and other data: A Monte Carlo approach, Phys.Rev. D66 (2002) 103511, [astro-ph/0205436]. [24] F. Feroz, M. Hobson, and M. Bridges, MultiNest: an efficient and robust Bayesian inference tool for cosmology and particle physics, MNRAS 398 (2009) 1601–1614. [25] I. Huston and K. A. Malik, Numerical calculation of second order perturbations, JCAP 0909 (2009) 019, [arXiv:0907.2917]. [26] I. Huston and A. J. Christopherson, Calculating Non-adiabatic Pressure Perturbations during Multi-field Inflation, Phys.Rev. D85 (2012) 063507, [arXiv:1111.6919]. [27] I. Huston and K. A. Malik, Second Order Perturbations During Inflation Beyond Slow-roll, JCAP 1110 (2011) 029, [arXiv:1103.0912]. [28] I. Huston and A. J. Christopherson, Isocurvature Perturbations and Reheating in Multi-Field Inflation, arXiv:1302.4298. [29] J. Frazer, Predictions in multifield models of inflation, arXiv:1303.3611. [30] R. Easther, J. Frazer, H. V. Peiris, and L. C. Price, Simple predictions from multifield inflationary models, Phys.Rev.Lett. 112 (2014) 161302, [arXiv:1312.4035]. [31] L. C. Price, H. V. Peiris, J. Frazer, and R. Easther, Gravitational wave consistency relations for multifield inflation, arXiv:1409.2498. [32] J. A. Adams, G. G. Ross, and S. Sarkar, Multiple inflation, Nucl.Phys. B503 (1997) 405–425, [hep-ph/9704286].

– 28 –

[33] P. Adshead, C. Dvorkin, W. Hu, and E. A. Lim, Non-Gaussianity from Step Features in the Inflationary Potential, Phys.Rev. D85 (2012) 023531, [arXiv:1110.3050]. [34] D. Seery, D. J. Mulryne, J. Frazer, and R. H. Ribeiro, Inflationary perturbation theory is geometrical optics in phase space, JCAP 1209 (2012) 010, [arXiv:1203.2635]. [35] D. H. Lyth and A. Riotto, Particle physics models of inflation and the cosmological density perturbation, Phys.Rept. 314 (1999) 1–146, [hep-ph/9807278]. [36] D. Langlois and S. Renaux-Petel, Perturbations in generalized multi-field inflation, JCAP 0804 (2008) 017, [arXiv:0801.1085]. [37] D. H. Lyth and A. R. Liddle, The primordial density perturbation: Cosmology, inflation and the origin of structure, . [38] D. Salopek, J. Bond, and J. M. Bardeen, Designing Density Fluctuation Spectra in Inflation, Phys.Rev. D40 (1989) 1753. [39] B. A. Bassett, S. Tsujikawa, and D. Wands, Inflation dynamics and reheating, Rev.Mod.Phys. 78 (2006) 537–589, [astro-ph/0507632]. [40] D. J. Mulryne, D. Seery, and D. Wesley, Moment transport equations for non-Gaussianity, JCAP 1001 (2010) 024, [arXiv:0909.2256]. [41] D. J. Mulryne, D. Seery, and D. Wesley, Moment transport equations for the primordial curvature perturbation, JCAP 1104 (2011) 030, [arXiv:1008.3159]. [42] D. J. Mulryne, Transporting non-Gaussianity from sub to super-horizon scales, JCAP 1309 (2013) 010, [arXiv:1302.3842]. [43] T. Bunch and P. Davies, Quantum Field Theory in de Sitter Space: Renormalization by Point Splitting, Proc.Roy.Soc.Lond. A360 (1978) 117–134. [44] J. M. Bardeen, Gauge Invariant Cosmological Perturbations, Phys.Rev. D22 (1980) 1882–1905. [45] J. M. Bardeen, P. J. Steinhardt, and M. S. Turner, Spontaneous creation of almost scale-free density perturbations in an inflationary universe, Phys. Rev. D 28 (Aug, 1983) 679–693. [46] S. Groot Nibbelink and B. van Tent, Density perturbations arising from multiple field slow roll inflation, hep-ph/0011325. [47] C. Gordon, D. Wands, B. A. Bassett, and R. Maartens, Adiabatic and entropy perturbations from inflation, Phys.Rev. D63 (2001) 023506, [astro-ph/0009131]. [48] S. Groot Nibbelink and B. van Tent, Scalar perturbations during multiple field slow-roll inflation, Class.Quant.Grav. 19 (2002) 613–640, [hep-ph/0107272]. [49] C. T. Byrnes and D. Wands, Curvature and isocurvature perturbations from two-field inflation in a slow-roll expansion, Phys.Rev. D74 (2006) 043529, [astro-ph/0605679]. [50] N. Bartolo, S. Matarrese, and A. Riotto, Adiabatic and isocurvature perturbations from inflation: Power spectra and consistency relations, Phys.Rev. D64 (2001) 123504, [astro-ph/0107502]. [51] D. Wands, N. Bartolo, S. Matarrese, and A. Riotto, An Observational test of two-field inflation, Phys.Rev. D66 (2002) 043520, [astro-ph/0205253]. [52] J. Garc´ıa-Bellido and D. Wands, Metric perturbations in two field inflation, Phys.Rev. D53 (1996) 5437–5445, [astro-ph/9511029]. [53] D. Wands, K. A. Malik, D. H. Lyth, and A. R. Liddle, A new approach to the evolution of cosmological perturbations on large scales, Phys.Rev. D62 (2000) 043527, [astro-ph/0003278]. [54] K. A. Malik, D. Wands, and C. Ungarelli, Large scale curvature and entropy perturbations for multiple interacting fluids, Phys.Rev. D67 (2003) 063516, [astro-ph/0211602].

– 29 –

[55] K. A. Malik, Cosmological perturbations in an inflationary universe, astro-ph/0101563. [56] K. A. Malik and D. Wands, Adiabatic and entropy perturbations with interacting fluids and fields, JCAP 0502 (2005) 007, [astro-ph/0411703]. [57] A. A. Starobinsky, Multicomponent de Sitter (Inflationary) Stages and the Generation of Perturbations, JETP Lett. 42 (1985) 152–155. [58] D. Lyth, Large Scale Energy Density Perturbations and Inflation, Phys.Rev. D31 (1985) 1792–1798. [59] M. Sasaki and E. D. Stewart, A General analytic formula for the spectral index of the density perturbations produced during inflation, Prog.Theor.Phys. 95 (1996) 71–78, [astro-ph/9507001]. [60] D. Salopek and J. Bond, Nonlinear evolution of long wavelength metric fluctuations in inflationary models, Phys.Rev. D42 (1990) 3936–3962. [61] M. Sasaki and T. Tanaka, Superhorizon scale dynamics of multiscalar inflation, Prog.Theor.Phys. 99 (1998) 763–782, [gr-qc/9801017]. [62] D. H. Lyth and Y. Rodriguez, The Inflationary prediction for primordial non-Gaussianity, Phys.Rev.Lett. 95 (2005) 121302, [astro-ph/0504045]. [63] F. Vernizzi and D. Wands, Non-gaussianities in two-field inflation, JCAP 0605 (2006) 019, [astro-ph/0603799]. [64] T. Battefeld and R. Easther, Non-Gaussianities in Multi-field Inflation, JCAP 0703 (2007) 020, [astro-ph/0610296]. [65] K. A. Malik and D. Wands, Cosmological perturbations, Phys.Rept. 475 (2009) 1–51, [arXiv:0809.4944]. [66] M. Dias, J. Frazer, and A. R. Liddle, Multifield consequences for D-brane inflation, JCAP 1206 (2012) 020, [arXiv:1203.3792]. [67] D. Seery and J. E. Lidsey, Primordial non-Gaussianities from multiple-field inflation, JCAP 0509 (2005) 011, [astro-ph/0506056]. [68] J. M. Maldacena, Non-Gaussian features of primordial fluctuations in single field inflationary models, JHEP 0305 (2003) 013, [astro-ph/0210603]. [69] L. Alabidi and D. H. Lyth, Inflation models and observation, JCAP 0605 (2006) 016, [astro-ph/0510441]. [70] C. T. Byrnes, M. Sasaki, and D. Wands, The primordial trispectrum from inflation, Phys.Rev. D74 (2006) 123519, [astro-ph/0611075]. [71] M. Dias, J. Frazer, and D. Seery, A code to compute observables in curved multifield models of inflation, In preparation. [72] A. D. Linde, Hybrid inflation, Phys.Rev. D49 (1994) 748–754, [astro-ph/9307002]. [73] E. J. Copeland, A. R. Liddle, D. H. Lyth, E. D. Stewart, and D. Wands, False vacuum inflation with Einstein gravity, Phys.Rev. D49 (1994) 6410–6433, [astro-ph/9401011]. [74] Z. Lalak, D. Langlois, S. Pokorski, and K. Turzynski, Curvature and isocurvature perturbations in two-field inflation, JCAP 0707 (2007) 014, [arXiv:0704.0212]. [75] A. Avgoustidis, S. Cremonini, A.-C. Davis, R. H. Ribeiro, K. Turzynski, et al., The Importance of Slow-roll Corrections During Multi-field Inflation, JCAP 1202 (2012) 038, [arXiv:1110.4081]. [76] A. R. Liddle, A. Mazumdar, and F. E. Schunck, Assisted inflation, Phys.Rev. D58 (1998) 061301, [astro-ph/9804177].

– 30 –

[77] P. Kanti and K. A. Olive, On the realization of assisted inflation, Phys.Rev. D60 (1999) 043502, [hep-ph/9903524]. [78] P. Kanti and K. A. Olive, Assisted chaotic inflation in higher dimensional theories, Phys.Lett. B464 (1999) 192–198, [hep-ph/9906331]. [79] N. Kaloper and A. R. Liddle, Dynamics and perturbations in assisted chaotic inflation, Phys.Rev. D61 (2000) 123513, [hep-ph/9910499]. [80] R. Easther and L. McAllister, Random matrices and the spectrum of N-flation, JCAP 0605 (2006) 018, [hep-th/0512102]. [81] S. Dimopoulos, S. Kachru, J. McGreevy, and J. G. Wacker, N-flation, JCAP 0808 (2008) 003, [hep-th/0507205]. [82] S. A. Kim and A. R. Liddle, Nflation: multi-field inflationary dynamics and perturbations, Phys.Rev. D74 (2006) 023513, [astro-ph/0605604]. [83] S. A. Kim and A. R. Liddle, Nflation: Non-Gaussianity in the horizon-crossing approximation, Phys.Rev. D74 (2006) 063522, [astro-ph/0608186]. [84] S. A. Kim and A. R. Liddle, Nflation: observable predictions from the random matrix mass spectrum, Phys.Rev. D76 (2007) 063515, [arXiv:0707.1982]. [85] D. Wenren, Tilt and Tensor-to-Scalar Ratio in Multifield Monodromy Inflation, arXiv:1405.1411. [86] X. Chen, R. Easther, and E. A. Lim, Large Non-Gaussianities in Single Field Inflation, JCAP 0706 (2007) 023, [astro-ph/0611645]. [87] X. Chen, R. Easther, and E. A. Lim, Generation and Characterization of Large Non-Gaussianities in Single Field Inflation, JCAP 0804 (2008) 010, [arXiv:0801.3295]. [88] R. Easther and R. Flauger, Planck Constraints on Monodromy Inflation, JCAP 1402 (2014) 037, [arXiv:1308.3736]. [89] P. D. Meerburg, D. N. Spergel, and B. D. Wandelt, Searching for Oscillations in the Primordial Power Spectrum: Perturbative Approach (Paper I), Phys.Rev. D89 (2014) 063536, [arXiv:1308.3704]. [90] P. D. Meerburg and D. N. Spergel, Searching for Oscillations in the Primordial Power Spectrum: Constraints from Planck (Paper II), Phys.Rev. D89 (2014) 063537, [arXiv:1308.3705]. [91] E. J. Copeland, A. Mazumdar, and N. Nunes, Generalized assisted inflation, Phys.Rev. D60 (1999) 083506, [astro-ph/9904309]. [92] D. Battefeld and S. Kawai, Preheating after N-flation, Phys.Rev. D77 (2008) 123507, [arXiv:0803.0321]. [93] T. C. Bachlechner, M. Dias, J. Frazer, and L. McAllister, A New Angle on Chaotic Inflation, arXiv:1404.7496. [94] C. Ringeval, P. Brax, C. van de Bruck, and A.-C. Davis, Boundary inflation and the wmap data, Phys.Rev. D73 (2006) 064035, [astro-ph/0509727]. [95] D. I. Kaiser and E. I. Sfakianakis, Multifield Inflation after Planck: The Case for Nonminimal Couplings, Phys.Rev.Lett. 112 (2014), no. 1 011302, [arXiv:1304.0363]. [96] S. Li and A. R. Liddle, Observational constraints on K-inflation models, JCAP 1210 (2012) 011, [arXiv:1204.6214]. [97] A. D. Linde, Chaotic Inflation, Phys.Lett. B129 (1983) 177–181.

– 31 –