Department of Physics Ben Gurion University of the Negev  Doctor of Philosophy
LaTeX

# Numerical analysis of the Primordial Power Spectrum for inflationary
potentials

Ira Wolfson

###### Abstract

In this work we study small field models of inflation, which, against previous
expectations, yield significant Gravitational Wave (GW) signal, while
reproducing other measured observable quantities in the Cosmic Microwave
Background (CMB). We numerically study these, using previously published
analytic works as general guidelines. We first discuss the framework necessary
to understand the model building procedure and some of its motivations. We
review the slow-roll paradigm, derive the slow-roll parameters and discuss
different formulations thereof. We further present the Lyth bound and its
theoretical descendants and finally, we outline the small/large field taxonomy
and their characterization in the current nomenclature.
We proceed to present our models and the methods used in their building and
examination. We employ MCMC simulations to evaluate model likelihoods and by
process of marginalization extract the most probable coefficients for these
inflationary potentials. An additional method applied is a multinomial fit,
where we create a functional correspondence between coefficients and
observables. This allows us to use the observable values directly to yield the
most likely coefficients. We compare the results of the two methods and
evaluate the level of tuning required for these models.
We discuss an apparent discrepancy between analytical approaches of evaluating
Primordial Power Spectrum (PPS) observables and the precise numerical results
in our models. We identify some of the sources of this discrepancy and remark
on their meaning in the age of precision cosmology.
Finally, we present the results of our study, for the most likely inflationary
models with polynomial potentials of degree 5, and 6. We demonstrate our
ability to produce potentials that yield GW with a tensor-to-scalar ratio
$r=0.03$. This is a realistic expectation of GW detection sensitivity in the
near future. Detecting GW of a primordial source will provide a direct
indication of the energy scale of inflation, and therefore an interesting
probe into physics beyond the standard model.

###### keywords:

LaTeX PhD Thesis Engineering University of Cambridge

For Netta, Tamar, Meir and Daya.

I hereby declare that except where specific reference is made to the work of
others, the contents of this dissertation are original and have not been
submitted in whole or in part for consideration for any other degree or
qualification in this, or any other university. This dissertation is my own
work and contains nothing which is the outcome of work done in collaboration
with others, except as specified in the text and Acknowledgements. This
dissertation contains fewer than 65,000 words including appendices,
bibliography, footnotes, tables and equations and has fewer than 150 figures.

###### Acknowledgements.

I want to acknowledge, first and foremost my advisor, Prof. Ramy Brustein, for
taking a chance on a young student coming hat-in-hand to start working in the
mostly barren field of numerical cosmology in Israel. As is the case
oftentimes, professional disagreements were not hard to come by, but the
personal and professional guidance was always available. Instrumental in my
scietific coding ‘upbringing’ were Dr. Rahul "write your own" Kumar, and Dr.
Daniel "you’re wrong" Ariad. These people were there when I was but a seedling
of the computational physicist I am today, and I am forever in their debt.
In the later stages of my doctoral adventure, several people were instrumental
in guiding me or, at the very least, providing a sounding board for my
craziness. Dr. Ido Ben-Dayan contributed time and effort to bring me up to
speed in all matters of observational and experimental cosmology. Snir
"Erdelyi" Cohen, Tomer "Deadpool" Ygael , and Yossi "Pigal" Naor were the
usual victims of my (temporary, physics induced) insanity.
When PhD has grown terminal, two main people were instrumental in ‘bringing it
all together’. Prof. Eiichiro Komatsu, who, with a 30-minute talk crystallized
many insights, and sharpened them into a point-like focus of what it meant to
be a physicist. And Dr. Yuval Ben-Abu, who never, not even for a minute, let
me rest. He always pushed me towards the next great challenge.

There is still a long list of thanks due, which is in no way exhaustive:
Prof. Antonio Riotto, for his support and invaluable advice.
Members of my advisory committee, professors Eduardo Guendelman and Uri
Keshet.
Sifu and Shermu, Donald and Cheryl-Lynn Rubbo.
And my family and friends, chief among them are my parents - Joshua and Pnina
Wolfson.

In this as in everything I owe my life to my wife, Netta Schramm.

###### Contents

  1. 1 Scientific Introduction - State of the art
    1. 1.1 The Big Bang
    2. 1.2 Problems in heaven
      1. 1.2.1 Horizon
      2. 1.2.2 Flatness
      3. 1.2.3 Relic
    3. 1.3 Inflation
      1. 1.3.1 Features of inflation
    4. 1.4 Solving inflation
      1. 1.4.1 Friedmann equations
      2. 1.4.2 Scalar field inflation
      3. 1.4.3 Equations of motion
    5. 1.5 Slow-roll
      1. 1.5.1 Quantifying slow-roll
      2. 1.5.2 Cosmic perturbations
    6. 1.6 Primordial power spectrum
      1. 1.6.1 Matter power spectrum
      2. 1.6.2 Transfer function
      3. 1.6.3 Scalar power spectrum
    7. 1.7 Stewart-Lyth power spectrum
      1. 1.7.1 Derivation
      2. 1.7.2 Analytical term for index running $\alpha_{s}$
    8. 1.8 Less common formulations for $n_{s}$
      1. 1.8.1 Dodelson-Stewart
      2. 1.8.2 Stewart-Gong
      3. 1.8.3 Schwarz-Escalante-Garcia
    9. 1.9 Lyth bound
    10. 1.10 Small and Large field models taxonomy
      1. 1.10.1 Large field models
      2. 1.10.2 Small field models
  2. 2 Our Models
    1. 2.0.1 Inflaton potentials with $r=0.001$
    2. 2.0.2 Reduced parameter space
    3. 2.0.3 Inflaton potentials with $r=0.01$
    4. 2.0.4 Inflaton potentials with $r=0.03$
  3. 3 Methods
    1. 3.0.1 Coefficient extraction methods
    2. 3.0.2 Monte Carlo analysis of Cosmic Microwave Background with running of running
  4. 4 The INSANE code
    1. 4.1 Code architecture
      1. 4.1.1 Background evolution
      2. 4.1.2 Cosmic Perturbation
      3. 4.1.3 Usage
      4. 4.1.4 Benchmarking
    2. 4.2 Code improvements
      1. 4.2.1 Finding the mode ‘freeze-out’ point
      2. 4.2.2 Forgetting the phase
      3. 4.2.3 Pre-Integration
  5. 5 Fresh eyes on $n_{s}$
    1. 5.1 The $n_{s}$ discrepancy
      1. 5.1.1 Slow-roll hierarchy
      2. 5.1.2 Another degree of separation - slow-roll in potential language
      3. 5.1.3 Does the $\alpha_{s}$ formulation work?
    2. 5.2 A Green’s function approach
      1. 5.2.1 Possible explanations of the source of deviation between precise results and analytical estimates.
  6. 6 Results
    1. 6.0.1 Models with $r=0.001$
    2. 6.1 6th degree models with $r=0.01$
      1. 6.1.1 Results for degree six polynomials that yield $r=0.01$
      2. 6.1.2 Most likely potentials
      3. 6.1.3 Observable dependence
    3. 6.2 6th degree models with $r=0.03$
      1. 6.2.1 Methods
      2. 6.2.2 Results
  7. 7 Conclusion
  8. A Derivation of Mukhanov-Sasaki equation
  9. B The INSANE code - usage examples

###### List of Figures

  1. 1.1 The horizon problem
  2. 1.2 Solving flatness
  3. 2.1 Efold generation regimes
  4. 2.2 Efold generation as $r$ changes
  5. 2.3 Polynomial potentials
  6. 2.4 Efold generation for models with $r\simeq 0.001,0.01,0.03$
  7. 3.1 Covering $n_{s}-\alpha_{s}$ with $r=0.01$ models.
  8. 3.2 Gaussian fitting of coefficients
  9. 3.3 Power spectra in $k$-space
  10. 3.4 Power spectra in $C_{l}$ representation.
  11. 4.1 Difference between cosmic time and efold based integration
  12. 4.2 $n_{s}$ benchmarking in power law inflation
  13. 4.3 Heuristic $\alpha_{s}$ benchmarking in power law inflation
  14. 4.4 Benchmarking of quadratic potentials
  15. 4.5 Freeze-out point
  16. 5.1 $n_{s}$ errors
  17. 5.2 Examples of power spectra discrepancy
  18. 5.3 SEG $n_{s}$ evaluation map
  19. 5.4 $n_{s}$ discrpency with different formulation
  20. 5.5 Within SEG approximation still yield $n_{s}$ discrepancy.
  21. 5.6 Satisfying the condition $\epsilon_{1}|\epsilon_{2}|\times\Delta N<10^{-2\sim 3}$ still yields discrepancy.
  22. 6.1 Covering $n_{s}-\alpha_{s}$ with $r=0.001$ models.
  23. 6.2 Mapping $n_{s}-\alpha_{s}$ phase-space with $\eta_{0}=0$
  24. 6.3 PPS reconstruction for most likely $r=0.01$ model.
  25. 6.4 Dependence of $\beta$ on other observables.
  26. 6.5 Covering $n_{s}-\alpha_{s}$ with $r=0.003$ models.
  27. 6.6 Likelihood analysis for $r=0.03$ models.
  28. 6.7 PPS results of most likely models with $r=0.001,0.01,0.03$.

###### List of Tables

  1. 1.1 Comparison of slow-roll parameters SEG to SG
  2. 3.1 MCMC results
  3. 4.1 Background evolution parameters
  4. 4.2 Background evolution output
  5. 5.1 $n_{s}$ discrepancy for a small sample of models.
  6. 5.2 $n_{s}$ discrepancy of $1\%$, for likely models.
  7. 5.3 Leading slow-roll parameters, numerical vs. analytic.
  8. 5.4 $n_{s}$ derived in different methods
  9. 6.1 Comparison of PPS from different methods
  10. 6.2 Most likely coefficients

##  Chapter 1 Scientific Introduction - State of the art

Even from the infancy stages of humanity, the question of origins was an ever-
present one. Every religion in existence and in fact, every culture,
incorporates an origin story for the universe.

Up until the new era of scientific renaissance, this field was dominated
predominantly by theologists, mystics and story-tellers. Some may attribute
this to proper education being exclusive to members of the upper religious
castes.

###  1.1 The evidence for a hot Big Bang

It is not until Einstein’s theory of relativity that we can construct an
origin story of sufficient precision, to exclude divine intervention. At least
a manifest one.

In the early 20th century, Alexander Friedmann, using Einstein’s own theory of
relativity, was able to derive an equation of motion for the evolution of the
universeFriedmannEqs; 1999GReGr..31.1991F . It took about seven years,
Friedmann’s death, and Georges Lemaître’s independent work
1927ASSB...47...49L; 1931MNRAS..91..483L, for Einstein to accept the idea that
his new physics proposed a mechanism for the evolution of the cosmos itself.
Coupled with Edwin Hubble’s observations and the Hubble law HubbleLaw, this
initiated a flood of scientific theories regarding the origin of our universe,
with three main competing scenarios:

  * •

A static universe - Most famously connected to Einstein’s “biggest blunder",
the cosmological constant Einstein:1917:KBA.

  * •

An ever expanding universe where matter is constantly created, to fill in the
newly created space.

  * •

An evolving universe which is currently expanding, originated by a primordial
explosion-like event.

The last of which was ridiculed by the pre-eminent cosmologist Fred Hoyle, who
coined the then derogatory phrase “Big Bang theory" - in response to the
notion of matter and energy created ex-nihilo at the onset of our universe, to
give rise to rapid expansion.

This last “Big Bang theory", became the lead contender in 1965 when Arno
Penzias and Robert Wilson, of Bell labs, stumbled upon the Cosmic Microwave
Background Radiation (CMBR or CMB) Penzias_Wilson1965. Dicke, Peebles and
Wilkinson, themselves working on a microwave band antenna at the time, were
able to provide the theoretical framework that would explain this as a remnant
of the “Big Bang" DickeWilkinsonPeebles1965. This, together with different
mounting evidence of light-element abundance in the universe was sufficient
confirmation of the “Big Bang" as the most probable physical scenario.

###  1.2 Problems with the old big bang model

Several measurements of the CMB temperature spectrum were made starting with
Penzias and Wilson Penzias_Wilson1965, through Thaddeus Thaddeus:1972an and
others, but the first overwhelming evidence of the black-body spectrum of the
CMB was found by the Far-Infrared Absolute Spectrophotometer (FIRAS)
instrument on the Cosmic Background Explorer (COBE) satellite Mather:1998gm.
From these measurements, with systematic dominated error margins of no more
than 0.3% at peak brightness, and an rms value of 0.01%. This extraordinary
precision in measuring the black-body nature of the CMB put the "Big Bang
theory" on as sound a footing as anyone could hope for. However, as with all
scientific discoveries, one answer gives rise to a multitude of other
problems.

####  1.2.1 The horizon problem

The horizon problem can roughly be put as the following: if we are now at the
light horizon of the ancient "Big Bang" event, the information that is
propagated to us from one end of the CMB, cannot propagate to the other end of
the CMB. Thus these two patches of CMB cannot possibly have exchanged
information during the expansion period of the universe (See Fig. 1.1).
However, they are thermalized and homogeneous to one part in $10^{5}$. If they
can’t exchange information how is that possible?

Fig. 1.1:  Light emitted from the "Big Bang" event is just now reaching our
eyes, after travelling a distance equal to the comoving distance
$\chi(t)=c\int_{t_{BB}}^{t_{now}}\frac{dt^{\prime}}{a(t^{\prime})}$. There are
regions in the CMB surface that are thus causally disconnected but are
nonetheless homogeneous to one part in $10^{5}$. This is suggestive of
information exchange, where no exchange is possible.

Formalising this requires to consider the Friedmann-Robertson-Walker (FRW)
metric:

$\displaystyle ds^{2}=-c^{2}dt^{2}+a(t)^{2}d{\bf x}^{2}.$ (1.1)

When one considers the distance travelled by light, one sets the interval to
zero, thus for light travelling along the $x$ axis yielding:

$\displaystyle x=\int dx=\int c\frac{dt^{\prime}}{a(t^{\prime})}.$ (1.2)

So the comoving distance is defined as:

$\displaystyle\chi(t)=\int_{t_{0}}^{t_{now}}c\frac{dt^{\prime}}{a(t^{\prime})}\Rightarrow\chi(a)=\int_{a}^{1}\frac{da^{\prime}}{a^{\prime
2}H(a^{\prime})},$ (1.3)

where the speed of light $c$ is set to $1$, $a(t)$ is the scale factor of the
FRW metric, which is set such that the current scale is $a(t_{now})=1$ and $H$
is the Hubble parameter, which varies for different historical stages of the
universe. One of the useful quantities to define is conformal time, in which
we rewrite the FRW metric as follows:

$\displaystyle ds^{2}=dt^{2}-a^{2}d{\bf x}^{2}\equiv a^{2}(d\tau^{2}-d{\bf
x}^{2}),$ (1.4)

so $\int cd\tau$ is the comoving distance. So, roughly speaking, unless we are
in a closed universe which is enclosed within a fully connected sphere of
radius no more than the observed $\chi$, there are patches in the CMB that are
causally disconnected. Since these patches are unable to exchange information
since the Big Bang event, how can they be in thermal equilibrium with each
other? Yet we observe them to be in thermal equilibrium up to one part in
$10^{5}$.

####  1.2.2 The flatness problem

Suppose we want to measure distances between different CMB patches. Since we
are limited to a two-dimensional spherical view of the CMB, it is intuitive to
measure distance by the connection between angle and arc length:

$\displaystyle l=d_{A}\theta,$ (1.5)

where $l$ is the arc length, $d_{A}$ is called the angular diameter distance,
and $\theta$ is the angle subtended. In a completely flat universe Eq. (1.3)
is connected to $d_{A}$ by:

$\displaystyle d_{A}^{flat}=a\chi.$ (1.6)

However, if we live in an open or closed universe in which the curvature is
non-zero, the proper generalization is given by:

$\displaystyle
d_{A}=\left\\{\begin{array}[]{lcr}\frac{a}{H_{0}}\sqrt{|\Omega_{k}|}\sinh\left(\sqrt{\Omega_{k}}H_{0}\chi\right)&&\Omega_{K}>0\\\
\frac{a}{H_{0}}\sqrt{|\Omega_{k}|}\sin\left(\sqrt{-\Omega_{k}}H_{0}\chi\right)&&\Omega_{K}<0\end{array}\right.,$
(1.9)

where $\Omega_{k}$ is the curvature density, and $H_{0}$ is the current Hubble
parameter. While the energy density in radiation scales as $a^{-4}$, and
matter density scales as $a^{-3}$, the energy density in curvature scales as
$a^{-2}$. This means that if we can construct a scale ladder spanning several
epochs, in which $a$ sufficiently changes, we can assess the energy density in
curvature. Such a ladder is provided to us by Type Ia Supernovae Riess:1994nx;
Dunkley:2008ie, and so we are able to evaluate the luminosity distance
$d_{L}=\frac{\chi}{a}$ to an object and compare it to the angular distance
$d_{A}$. By observing a large sample of these pairs we can extract the overall
behaviour of $a$ at different times, and which of the three angular distance
functions are the best fit. Using this method the most recent Planck data sets
the curvature density $\Omega_{k}$ at

$\displaystyle\Omega_{k}=0.001\pm 0.002.$ (1.10)

This value is very close to zero, which means a flat universe. Of all possible
values the curvature density can take, which can be of the order $\pm 1$,
corresponding to a universe dominated by the curvature component, why is it
that we find our universe so close to a flat one? This is not the most finely-
tuned quantity in our universe Carroll:2014uoa, but it is fine-tuned
nonetheless.

####  1.2.3 The relic problem

Every Grand-Unified-Theory (GUT) that includes electromagnetism inevitably
produces super-heavy magnetic monopoles Zeldovich:1978wj; PhysRevLett.43.1365;
PhysRevLett.44.631; mukhanov2005physical, typically of magnetic charge
$\propto\frac{1}{e}$, where $e$ is the fundamental electric charge and is of
order $e\simeq 10^{-19}\;C$. Such a large magnetic charge would have
observable effects on the universe and would be easily detectable. However so
far we have not found any such effect and deduce the lack of any monopoles in
the observed universe. During the GUT phase of the universe, the energy
density in magnetic monopoles should have been such that today we should have
measured $\Omega_{\mathrm{Mono}}\sim 10^{13}$. How is it possible that we
detect no monopoles then?

###  1.3 Inflation saves the day

It was Starobinsky Starobinsky:1979ty closely followed by Guth GuthInflation,
that suggested the idea of inflation, an epoch of rapid and accelerated
expansion of the universe. This mechanism at once fixes the above problems.
Let us consider:

$\displaystyle c\tau=\int\frac{da}{a}\frac{c}{aH},$ (1.11)

which is the logarithmic integral of the comoving Hubble radius. When the
distance between two objects is larger than the comoving Hubble radius, they
cannot currently be in causal connection with each other. If particles are
however separated from each other by more than $c\tau$ they could have never
exchanged information. If we can find a way in which two particles are
currently separated more than the comoving Hubble radius, but separated by
less than $c\tau$, we can solve the horizon problem. In other words, we want
today:

$\displaystyle\frac{1}{aH(a)}<1,$ (1.12)

and a scale $a^{\prime}$ at some time in the evolution of the universe such
that

$\displaystyle\frac{1}{a^{\prime}H(a^{\prime})}>1.$ (1.13)

Thus there is a phase in cosmic history where the term $aH$ is increasing:

$\displaystyle\frac{d(aH)}{dt}=\frac{d^{2}a}{dt^{2}}>0,$ (1.14)

which means the scale factor of the universe was growing in an accelerated
fashion. However, this could not have occurred during matter or radiation
dominated eras, so there is something else that gives rise to this phenomenon.
This solution to the horizon problem naturally solves the flatness problem as
well as the relic problem. Rarefying the energy density stored in curvature,
by accelerated growth of the scale of the universe, the currently observed
curvature strongly tends towards zero. This is similar to the process of
deriving a smooth function. By sufficiently ’zooming-in’ ($dx\rightarrow 0$),
the curvature vanishes, and we are left with the linear term alone. Similarly,
a universe initially populated by relics could have been sufficiently diluted
such that we do not observe any relics today.

Fig. 1.2: By continually rescaling the universe, distances which might have
been connected by a curved interval are now connected by a less curved
interval, thus diluting the energy density in curvature. This is analogous to
’zooming-in’ in order to examine the linear approximation of a function at a
point.

####  1.3.1 Features of inflation

There are additional consequences for a period of accelerated growth of the
scale factor, beyond the observations suggested above. Perhaps one of the most
looked at today is the notion of primordial gravitational waves. In a later
section we will quantitatively look at the production of gravitational waves
from inflation. However, at this point we want to look at the general idea.
Gravitational waves are generated whenever the metric at some point in space
undergoes acceleration and is not completely uniform. To see that clearly we
consider small perturbations over the flat Minkowski metric:

$\displaystyle g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu},$ (1.15)

We then define the trace-reversed perturbation:

$\displaystyle\tilde{h}_{\mu\nu}=h_{\mu\nu}-h\eta_{\mu\nu}.$ (1.16)

Moving to the transverse-traceless gauge we get:

$\displaystyle\tilde{h}^{TT}_{\mu\nu}=h^{TT}_{\mu\nu}.$ (1.17)

By then using the Lorentz gauge we set:

$\displaystyle\partial_{\mu}\tilde{h}^{\mu\nu}=0$ (1.18)

We thus get the following a linearized version of the Einstein equation, which
is just a non-homogeneous wave equation:

$\displaystyle\square\tilde{h}_{\mu\nu}=-16\pi GT_{\mu\nu},$ (1.19)

where $T_{\mu\nu}$ is the stress-energy tensor. Solving this equation is
possible by finding the Green’s function for the wave equation, and
incorporating a retarded quantity $t_{r}=t-|{\bf x-y}|$. Finally, the term for
gravitational wave production is proportional to the second temporal
derivative of the Quadrupole moment of the energy density:

$\displaystyle I_{ij}(t)=\int y^{i}y^{j}T^{00}(t,{\bf y})d^{3}y\;;\hskip
50.0pt\tilde{h}_{ij}(t,{\bf x})=\frac{2G}{r}\frac{d^{2}I_{ij}}{dt^{2}}(t_{r})$
(1.20)

This is not the main focus of this manuscript, however we clearly see that
when we have an accelerating quadrupole moment, like in inflation, we have
gravitational wave production. The physics that source primordial GW is slight
spatial perturbation in curvature, sourced by slight quantum perturbation of
the inflaton field. The common notion is that these perturbations are small
such that the overall scale factor is given by $a(t)$ (no spatial dependence).
The way to derive these equations is to consider the FRW metric with slight
spatial perturbation such that:

$\displaystyle
g_{00}=-1\;;\;g_{ij}=a^{2}(t)\left(\begin{array}[]{ccc}1+h_{+}&h_{\times}&0\\\
h_{\times}&1+h_{+}&0\\\ 0&0&1\end{array}\right).$ (1.24)

Going on calculate the Christoffel symbols, Ricci tensor and scalar and
throwing out second amd higher order terms, we get a linearized EFE. The zero-
order contribution are the leading order Friedmann equations. Replacing the
spatial derivative terms with $k^{2}$ by going into spatial Fourier-space and
after some algebra we arrive at:

$\displaystyle\partial_{tt}h_{(+,\times)}+3H\partial_{t}h_{(+,\times)}+\frac{k^{2}h_{(+,\times)}}{a^{2}(t)}=0.$
(1.25)

Rephrased in conformal time, this takes on the follwoing form:

$\displaystyle
h^{\prime\prime}_{\zeta}+2\frac{a^{\prime}}{a}h^{\prime}_{\zeta}+k^{2}h_{\zeta}=0,$
(1.26)

where $\zeta\in(+,\times)$ and the prime repesents $\partial_{\tau}$. It is
beneficial to recast this equation into a harmonic oscillator form. To do this
we redefine:

$\displaystyle\tilde{h}=\frac{ah}{\sqrt{16\pi G}},$ (1.27)

where the numerical coefficient comes from physical Action considerations.
Disregarding this numerical coffeicient we have:

$\displaystyle h=\frac{\tilde{h}}{a},$ (1.28) $\displaystyle
h^{\prime}=\frac{\tilde{h}^{\prime}}{a}-\frac{a^{\prime}}{a^{2}}\tilde{h},$
(1.29) $\displaystyle
h^{\prime\prime}=\frac{\tilde{h}^{\prime\prime}}{a}-2\frac{a^{\prime}}{a^{2}}\tilde{h}^{\prime}-\frac{a^{\prime\prime}}{a^{2}}\tilde{h}+2\frac{(a^{\prime})^{2}}{a^{3}}\tilde{h}\;.$
(1.30)

Inserting these into Eq. (1.26), yields:

$\displaystyle\tilde{h}^{\prime\prime}+\left(k^{2}-\frac{a^{\prime\prime}}{a}\right)\tilde{h}=0.$
(1.31)

This equation is similar to the Mukhanov-Sasaki (MS) equation, to be discussed
later, that governs the production of scalar perturbation. However the pump
field in this equation is the scale factor $a$, whereas the pump field in the
MS equation is $z=a\dot{\phi}/H$.

###  1.4 Solving inflation

Let us first take a look at the simplest form of inflation: A constant energy
density term over an FRW background metric. This is also known as the de-
Sitter case. In this case the metric is given by
$ds^{2}=dt^{2}-a^{2}(t)d\bf{x}^{2}$ where the metric itself is time dependent.
Since the system is explicitly time dependent, using a Hamiltonian formalism
is disfavoured. Thus we use the Einstein Field Equations:

$\displaystyle G_{\mu\nu}=8\pi GT_{\mu\nu},$ (1.32)

along with the Lagrangian formulation for the equations of motion of the
associated fields. When we ’solve’ inflation we intend to say, we calculate
the evolution of all inflationary quantities, which in this case are
$t,H(t),\phi(t)$ from which we construct all associated quantities:
$a(t),\dot{a}(t),N(t),\dot{\phi}(t)$, as well as the slow-roll parameters as
defined
($\varepsilon_{H}=\frac{-\dot{H}}{H^{2}},\delta_{H}=\frac{\ddot{\phi}}{H\dot{\phi}}$,
etc.) versus as potential and derivative terms
($\epsilon_{V}=\frac{1}{2}\left(\frac{V^{\prime}}{V}\right)^{2},\delta_{V}=\frac{V^{\prime\prime}}{V}$
etc.) We also construct the pump field $z$ that is later used for the quantum
perturbations calculation:

$\displaystyle z=\frac{a\dot{\phi}}{H}.$ (1.33)

####  1.4.1 General Relativity, Friedmann equations in conformal and cosmic
time

General Relativity, has shifted our understanding of gravitation from that of
a force, to that of a characteristic of space-time, where the force of
gravitation is nothing but Newton’s first law as applied in curvilinear
coordinates(carroll2003spacetime, Chapter 1).

This is done by coupling a metric (i.e. a manifold that describes space-time)
to the total energy present.

This is encoded mathematically by using (psuedo) Riemannian geometry, most
commonly using the Einstein notation:

$R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}\mathbf{R}=8\pi GT_{\mu\nu},$ (1.34)

where $R^{\mu\nu}$ is the Ricci Tensor (2nd order Tensor form), and
$\mathbf{R}$ is the Ricci scalar which is given by:

$\mathbf{R}=g^{\mu\nu}R_{\mu\nu}.$ (1.35)

These might sometimes be simplified by using the Einstein tensor:

$\displaystyle G_{\mu\nu}\equiv R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}\mathbf{R},$
(1.36)

thus we arrive at the form:

$\displaystyle G_{\mu\nu}=8\pi GT_{\mu\nu}.$ (1.37)

These applied to the universe as a whole yield the Friedmann equations.

##### FRW metric and Friedmann equations - cosmic time

The Friedmann-Robertson-Walker metric (FRW metric) is given by allowing the
spatial coordinates to change scale, as a function of time. In Cartesian
coordinates, this takes on the form of:

$\displaystyle g_{\mu\nu}=\left(\begin{array}[]{cccc}-1&&&\\\ &a^{2}(t)&&\\\
&&a^{2}(t)&\\\ &&&a^{2}(t)\end{array}\right).$ (1.42)

After going through the process of deriving the Ricci Tensor and Scalar, this
metric gives two equations. Taking the temporal-temporal coordinate
($\mu\nu=00$), and any of the spatial-spatial coordinate ($\mu\nu=ii$) we
have:

$\displaystyle\left(\frac{\dot{a}}{a}\right)^{2}=\frac{8\pi G}{3}\rho,$ (1.43)
$\displaystyle-2\left(\frac{\ddot{a}}{a}\right)-\left(\frac{\dot{a}}{a}\right)^{2}=8\pi
GP,$ (1.44)

where we use the diagonalized Stress-Energy tensor:

$\displaystyle
T_{\mu\nu}=g_{\mu\alpha}T^{\alpha}\,_{\nu}=\left(\begin{array}[]{cccc}\rho&&&\\\
&a^{2}P&&\\\ &&a^{2}P&\\\ &&&a^{2}P\end{array}\right),$ (1.49)

with

$\displaystyle T^{\mu}\;_{\nu}=\left(\begin{array}[]{cccc}-\rho&&&\\\ &P&&\\\
&&P&\\\ &&&P\end{array}\right).$ (1.54)

This is where one usually points out the possibility of an accelerated
expansion/contraction solution to the Friedmann equation. This is given by an
energy density term, which is constant, i.e.:

$\displaystyle\left(\frac{\dot{a}}{a}\right)^{2}=\frac{8\pi
G}{3}\rho\Rightarrow\frac{\dot{a}}{a}=\sqrt{\frac{8\pi G\rho}{3}}=Const,$
(1.55) $\displaystyle\Rightarrow a(t)=\left(a_{0}e^{\pm\sqrt{\frac{8\pi
G\rho}{3}}t}\right).$

Ultimately we are dealing with the positive exponential solution since this
fits the observations, and solves the Horizon and Flatness problems as
previously explained.

However a constant $\rho$ means permanent inflation, so $\rho$ has to change
in time, albeit slowly enough to facilitate exponential inflation. This is
called the slow-roll condition. It will be discussed at length in section 1.5.

##### FRW metric and Friedmann equations - conformal time

Moving to conformal time, puts the time coordinate $\tau$ on equal footing
with the other spatial coordinates by defining:

$\displaystyle dt=ad\tau.$ (1.56)

Thus the interval, or line element, takes on the form of:

$\displaystyle
dS^{2}=-dt^{2}+a^{2}(t)\left(dx^{2}+dy^{2}+dz^{2}\right)=a^{2}(\tau)\left(-d\tau^{2}+dx^{2}+dy^{2}+dz^{2}\right),$
(1.57)

so, in essence we are reducing the problem to Minkowski space over a time
dependant overall scale factor.
In this context, the metric is given by:

$\displaystyle g_{\mu\nu}=\left(\begin{array}[]{cccc}-a^{2}(\tau)&&&\\\
&a^{2}(\tau)&&\\\ &&a^{2}(\tau)&\\\
&&&a^{2}(\tau)\end{array}\right)=a^{2}(\tau)\left(\begin{array}[]{cccc}-1&&&\\\
&1&&\\\ &&1&\\\ &&&1\end{array}\right),$ (1.66)

and the derived Friedmann equations are given by:

$\displaystyle\left(\frac{a^{\prime}}{a}\right)^{2}=\frac{8\pi G\rho}{3}a^{2}$
(1.67)
$\displaystyle\frac{\left(-2aa^{\prime\prime}+{a^{\prime}}^{2}\right)}{a^{2}}=8GP\pi
a^{2}$ (1.68)

It is interesting to see that a constant $\rho$ in conformal time produces the
following solution:

$\displaystyle a(\tau)=\frac{1}{1-\sqrt{\frac{8\pi G\rho}{3}}\tau}$ (1.69)

which explodes to infinity at $\tau=\sqrt{\frac{3}{8\pi G\rho}}$, thus
$a(\tau)$ experiences greater than exponential growth in conformal time. While
analytically this stands on equal footing with a cosmic time solution,
numerical integration considerations will later lead us to prefer cosmic over
conformal time.

####  1.4.2 The scalar field inflation

Inflation is driven by an energy density term that permeates all of space and
is slowly varying, even taking into account the rapid expansion of space.

There might be any number of physical drivers that manifest this behavior, but
arguably the simplest form of mechanism which gives rise to inflation is the
scalar field one. It has been shown by Albrecht and Steinhardt, that a Higgs-
like field, can produce such a scalar effective field Albrecht:1982wi. The
simple Lagrangian (density) form of a scalar field is given by:

$\mathcal{L}=\frac{1}{2}\left(\partial_{\mu}\phi\right)^{2}-\frac{m^{2}}{2}\phi^{2}+\mathcal{O}\left(\phi^{3}\right)$
(1.70)

Where the linear term in the potential can be set to zero because we can
always redefine $\tilde{\phi}$ such that the derivatives are unchanged and the
quadratic term is changed accordingly. It is a similar procedure acting on the
interplay between the second and fourth powers of the Lagrangian, that
constitutes the Goldstone gauge Boson procedure NamboGoldstonBoson:1960;
Goldston:1962, and ultimately the Higgs Mechanism EnglertBrout:1964;
Higgs:1964; theRestOfHiggs:1964.

An almost straightforward generalization, using Legendre transform, and going
to covariant form yields (dodelson:2003, Page 152):

$\displaystyle T^{\alpha}\,_{\beta}=g^{\alpha\nu}\frac{\partial\phi}{\partial
x^{\mu}}\frac{\partial\phi}{\partial
x^{\beta}}-g^{\alpha}\,_{\beta}\left[\frac{1}{2}g^{\mu\nu}\frac{\partial\phi}{\partial
x^{\mu}}\frac{\partial\phi}{\partial x^{\nu}}+V\left(\phi\right)\right],$
(1.71)

where $V(\phi)$ is the interaction potential. Considering this field to be
mostly homogeneous, and disregarding the small perturbations to this quantity,
the spatial derivatives vanish, and we are left with:

$\displaystyle
T^{(0)\alpha}\,_{\beta}=-g^{\alpha}\,_{0}g^{\beta}\,_{0}\left(\frac{d\phi^{0}}{dt}\right)^{2}+g^{\alpha}\,_{\beta}\left[\frac{1}{2}\left(\frac{d\phi^{0}}{dt}\right)^{2}-V\left(\phi^{0}\right)\right].$
(1.72)

Identifying this term with Eq. (1.54), we may now denote:

$\displaystyle\rho=\frac{1}{2}\left(\frac{d\phi^{0}}{dt}\right)^{2}$
$\displaystyle+V\left(\phi^{0}\right),$ (1.73) $\displaystyle
P=\frac{1}{2}\left(\frac{d\phi^{0}}{dt}\right)^{2}-V\left(\phi^{0}\right),$
(1.74)

which is in cosmic time FRW. The conformal time version is given in Eq.
(1.81,1.82)
from here on we will forego the superscript (0) for simplicity

####  1.4.3 Equations of motion

When one is tasked with simulating inflation, whether numerically or
analytically, the first step is to evaluate the equations of motion for the
scale factor ,$a$ , and the driving fields, in our case the inflaton $\phi$.
This is done by following these steps:

  * •

Derive Friedmann equations from metric assumptions (FRW/FRW+conformal).

  * •

Replace the quantities $\rho$ and $P$ with the driving fields equivalents.

  * •

Identify the set of differential equations to work with.

We now proceed to do so, for a conformal background geometry, driven by a
scalar field.

The metric is:

$\displaystyle g_{\mu\nu}=\left(\begin{array}[]{cccc}-a(\tau)^{2}&&&\\\
&a(\tau)^{2}&&\\\ &&a(\tau)^{2}&\\\ &&&a(\tau)^{2}\end{array}\right)$ (1.79)

and the Stress-Energy tensor for a scalar field is given by:

$\displaystyle T^{\alpha}\,_{\beta}=g^{\alpha\nu}\frac{\partial\phi}{\partial
x^{\nu}}\frac{\partial\phi}{\partial
x^{\beta}}-\delta^{\alpha}\;_{\beta}\left[\frac{1}{2}g^{\mu\nu}\frac{\partial\phi}{\partial
x^{\mu}}\frac{\partial\phi}{\partial x^{\nu}}+V\left(\phi\right)\right].$
(1.80)

These imply - in conformal time:

$\displaystyle
T^{(0)0}\;_{0}\equiv\rho=\frac{1}{2a^{2}}\left(\frac{d\phi^{0}}{d\tau}\right)^{2}+V(\phi^{(0)}),$
(1.81) $\displaystyle T^{(0)i}\;_{i}\equiv
P=\frac{1}{2a^{2}}\left(\frac{d\phi^{0}}{d\tau}\right)^{2}-V(\phi^{(0)}).$
(1.82)

With the two equations which we get from the Einstein field equations (EFE):

$\displaystyle\left(\frac{1}{a}\frac{da}{d\tau}\right)^{2}=\frac{8\pi
G\rho}{3}a^{2},$ (1.83)
$\displaystyle\frac{\left(-2a\cdot\frac{d^{2}a}{d\tau^{2}}+(\frac{da}{d\tau})^{2}\right)}{a^{2}}=8GP\pi
a^{2},$ (1.84)

one can work through the algebra, to get:

$\displaystyle\left\\{\begin{array}[]{l}\frac{d^{2}\phi}{d\tau^{2}}+2\left(\frac{1}{a}\frac{da}{d\tau}\right)\frac{d\phi}{d\tau}+a^{2}V^{\prime}=0\\\
\\\ \frac{d^{2}a}{d\tau^{2}}=\frac{4\pi
Ga}{3}\left[4a^{2}V-\left(\frac{d\phi}{d\tau}\right)^{2}\right]\\\ \\\
\frac{da}{d\tau}=\sqrt{\frac{8\pi
G}{3}\left[\frac{a^{2}}{2}\left(\frac{d\phi}{d\tau}\right)^{2}+a^{4}V(\phi)\right]}.\end{array}\right.$
(1.90)

Since the first two are second order differential equations, and since $a$ is
not observable, it might be useful to integrate over $A$ as defined by
$A\equiv\frac{d\ln(a)}{d\tau}$. For this aim the second equation is given by:

$\displaystyle\frac{dA}{d\tau}=A^{2}\left(\frac{V-\left(\frac{d\phi}{d\tau}\right)^{2}}{V+\frac{1}{2}\left(\frac{d\phi}{d\tau}\right)^{2}}\right).$
(1.91)

To conclude, in conformal time, a "good" integration scheme might use this set
of equations:

$\displaystyle\left\\{\begin{array}[]{ccc}\phi^{\prime\prime}&=&-2A\phi^{\prime}-\frac{dV}{d\phi}\\\
a^{\prime}&=&A\cdot a\\\ A^{\prime}&=&A^{2}\left(\frac{V-\phi^{\prime
2}}{V+\tfrac{1}{2}\phi^{\prime 2}}\right)\end{array}\right.$ (1.95)

Deriving the analogue in cosmic time gives:

$\displaystyle\left\\{\begin{array}[]{l}\ddot{\phi}+3H\dot{\phi}+V^{\prime}=0\\\
\\\ \dot{a}=\sqrt{\frac{\left(\dot{\phi}^{2}+2V\right)}{6M_{pl}^{2}}}a\\\ \\\
\ddot{a}=-\frac{a}{3M_{pl}^{2}}\left[\dot{\phi}^{2}-V\right]\end{array}\right.\;,$
(1.101)

with an equation for $\dot{H}$ given by:

$\displaystyle\dot{H}=-\frac{\dot{\phi}^{2}}{2M^{2}_{pl}}\;,$ (1.102)

where the quantity $\sqrt{\frac{1}{8\pi G}}$ is neatly wrapped in the term
$\sqrt{\frac{1}{8\pi G}}=m_{pl}$, the reduced Planck mass.
Thus we are left to solve this set of coupled ordinary differential equations
(ODE), to resolve the evolution of the background geometry:

$\displaystyle\left\\{\begin{array}[]{l}\ddot{\phi}=-3H\dot{\phi}-V^{\prime}\\\
\dot{H}=-\frac{\dot{\phi}^{2}}{2m^{2}_{pl}}\end{array}\right.$ (1.105)

###  1.5 The slow-roll paradigm

As was alluded to before when we have a constant energy density term $\rho$
the solution for the temporal-temporal Friedmann equation over an FRW metric
is given by:

$\displaystyle H^{2}=\frac{8\pi G}{3}\rho,$ (1.106)

Since $\rho$ is just some constant we can now unpack $H$ and solve for $a(t)$:

$\displaystyle H=\frac{\dot{a}}{a}=\sqrt{\frac{8\pi G}{3}\rho}$ (1.108)
$\displaystyle\Rightarrow
a=\sum_{j=1}^{2}A_{j}\exp\left((-1)^{j}\sqrt{\frac{8\pi G}{3}\rho}t\right).$
(1.109)

It is customary to disregard the exponentially suppressed solution since they
would be undetected in 2-3 efolds of inflation. Thus we arrive at the pure de-
Sitter (dS) solution:

$\displaystyle a(t)=a_{0}\exp\left(\sqrt{\frac{8\pi G}{3}\rho}\cdot t\right).$
(1.110)

This gives rise to the idea of inflation as a phase of an exponent-like
evolution of the universe, with an energy density that changes slow enough, to
facilitate a perturbative solution, over the baseline of a pure de-Sitter
evolution.

####  1.5.1 Quantifying slow-roll

One way to quantify slow-roll would be to assume the basic solution for $a(t)$
is an exponent function, and then use perturbations over that baseline such
that $a_{0}=\exp(C\cdot t)$, $a_{1}=\lambda f(t)$, and go on in a perturbative
expansion. However, the usual formalism that is used to quantify slow-roll is
based on the idea that the Hubble parameter $H$ should be very slowly changing
with respect to cosmic time. This means that the perturbative degrees are
applied to $H$ rather than $a$. It ,therefore, makes sense to write $H$ as
some power series in $t$.

##### Taylor expansion

As it turns out since $H$ has units of $1/T$, i.e. $H$ is the reciprocal of
cosmic time, it makes sense to develop the function $1/H$ in a Taylor series
around $H_{0}$, where $H_{0}$ is a constant. Thus we have:

$\displaystyle\frac{1}{H}\simeq\frac{1}{H_{0}}-\left.\frac{\dot{H}}{H^{2}}\right|_{H=H_{0}}t-\left.\frac{\ddot{H}H^{2}-2H\dot{H^{2}}}{H^{4}}\right|_{H=H_{0}}\frac{t^{2}}{2}+\ldots$
(1.111)

The first slow-roll parameter can be read directly from the above expansion:

$\displaystyle\varepsilon_{H}\equiv-\frac{\dot{H}}{H^{2}},$ (1.112)

and it is usually, within the framework of slow-roll inflation, small and
positive since it is understood that, naturally, the average energy density
$\rho$ should decline as the universe inflates. This leads $H$ to slightly
decrease, thus slightly increasing $1/H$. We are left, then, with a positive
first-order contribution to the Hubble horizon. Formally slow-roll takes place
when

$\displaystyle 0<\varepsilon_{H}\ll 1.$ (1.113)

The second term in the Taylor expansion needs a bit of handling to find a more
aesthetic term.

##### The slow-roll parameters $\varepsilon_{H}$ and $\delta_{H}$

While the first slow-roll parameter $\varepsilon_{H}$ is a small, non-negative
quantity, given simply by (1.112), in order to simplify the second slow-roll
parameter we need to refer to the Klein-Gordon equation for a time-dependent
scalar field (1.101, top equation). Reason states that where the field $\phi$
is slowly rolling, the drag term that is given by $3H\dot{\phi}$, roughly
counteracts the force term $\frac{dV}{d\phi}$, leading to a terminal velocity-
like state. In that scenario the fields ’acceleration’ is roughly zero, so we
get:

$\displaystyle\ddot{\phi}=-3H\dot{\phi}-\frac{dV}{d\phi}\simeq 0,$ (1.114)

differently stated as

$\displaystyle|3H\dot{\phi}|\simeq\left|\frac{dV}{d\phi}\right|.$ (1.115)

It makes sense, then, to look at the dimensionless quantity

$\displaystyle\delta_{H}\equiv\frac{\ddot{\phi}}{H\dot{\phi}},$ (1.116)

and set an additional slow-roll condition as

$\displaystyle\left|\delta_{H}\right|\ll 1.$ (1.117)

The connection with the Taylor expansion is not immediately obvious, so we
will explicitly derive it. The second order term in the Taylor expansion for
$1/H$ is given by:

$\displaystyle\left.\frac{\ddot{H}H^{2}-2H\dot{H^{2}}}{H^{4}}\right|_{H=H_{0}}\frac{t^{2}}{2}.$
(1.118)

Let us examine this term:

$\displaystyle\frac{\ddot{H}}{H^{2}}=\frac{1}{H^{2}}\frac{d}{dt}\left(-\frac{\dot{\phi}^{2}}{2}\right)$
(1.119)
$\displaystyle=\frac{-1}{H^{2}}\dot{\phi}\ddot{\phi}=-\frac{\dot{\phi}^{2}}{H}\frac{\ddot{\phi}}{H\dot{\phi}}$
$\displaystyle=2H\frac{\dot{H}}{H^{2}}\delta_{H}=-2H\varepsilon_{H}\delta_{H}.$

Thus the aforementioned term (1.118) is given by:

$\displaystyle\left.\frac{\ddot{H}H^{2}-2H\dot{H^{2}}}{H^{4}}\right|_{H=H_{0}}\frac{t^{2}}{2}=-H\left.\varepsilon_{H}\left(\delta_{H}+\varepsilon_{H}\right)\right|_{H=H_{0}}t^{2},$
(1.120)

and the Taylor expansion up to second order is given by:

$\displaystyle\frac{1}{H_{0}}+\varepsilon_{H}t+\varepsilon_{H}\left(\delta_{H}+\varepsilon_{H}\right)H_{0}t^{2},$
(1.121)

where the expansion is around $H=H_{0}$.

##### Switching gears into $\epsilon_{V}$ and $\eta_{V}$

It is beneficial for theoretical physicists to be able to connect the slow-
roll parameters directly to the inflationary potential. This is done because
as we shall see the primordial power spectrum can be quantified up to several
percent precision, directly by the slow-roll parameters. The first slow-roll
parameter in potential and potential derivative language is usually given by

$\displaystyle\varepsilon_{V}\equiv\frac{1}{2}\left(\frac{V^{\prime}}{V}\right)^{2}.$
(1.122)

The two definitions, $\varepsilon_{H}$ and $\varepsilon_{V}$ coincide in both
extreme limits, i.e. when slow-roll is extremely slow and $\ddot{\phi}$
strictly vanishes and $\dot{a}\propto\ddot{a}$, and at the other end of
inflation, where $H$ is no longer slowly varying, and the accelerated
evolution of the scale factor vanishes ($\ddot{a}=0$). The derivation of the
slow-roll parameter is as follows: In slow-roll inflation, the potential
overpowers the kinetic energy immensely thus:

$\displaystyle V\gg\dot{\phi}^{2}$

This simplifies greatly the two major players of inflation, the first of which
is:

$\displaystyle\begin{array}[]{lcr}H^{2}=\frac{\rho}{3}&;(\rho=\frac{\phi^{2}}{2}+V);\Rightarrow&H^{2}\simeq\frac{V}{3}\;.\\\
\end{array}$

The second equation, the equation of motion for the scalar field is usually
given by:

$\displaystyle\ddot{\phi}+3H\dot{\phi}+V^{\prime}=0,$

which is the result of derivation with respect to time of the more primitive
form of:

$\displaystyle\left(\frac{\dot{\phi}^{2}}{2}+V\right)-3H^{2}=0.$

So now, seeing as the first term on the left is dominated by $V$ we can set:

$\displaystyle 0=\left(\frac{\dot{\phi}^{2}}{2}+V\right)-3H^{2}\simeq
V-3H^{2}$

And deriving with respect do time we have:

$\displaystyle\dot{V}-6H\dot{H}\simeq 0$

Now we remind ourselves that $\dot{V}=V^{\prime}\dot{\phi}$ and that
$\dot{H}=-\frac{\dot{\phi}^{2}}{2}$ to get:

$\displaystyle\left(V^{\prime}\dot{\phi}+3H\dot{\phi}^{2}\right)\simeq 0,$

or in its elegant form:

$\displaystyle-V^{\prime}\simeq 3H\dot{\phi}.$

So, to sum it up we have:

$\displaystyle\left\\{\begin{array}[]{c}H^{2}\simeq\frac{V}{3}\\\
3H\dot{\phi}\simeq-V^{\prime}\end{array}\right.$

From the second we get $\dot{\phi}=-\frac{V^{\prime}}{3H}$, and since
$\dot{H}=-\frac{\dot{\phi}^{2}}{2}$ we have:

$\displaystyle\dot{H}=-\frac{V^{\prime 2}}{18H^{2}},$

and so we have:

$\displaystyle\varepsilon_{H}\equiv-\frac{\dot{H}}{H^{2}}=\frac{V^{\prime
2}}{18H^{4}}=\frac{V^{\prime
2}}{18}\cdot\frac{9}{V^{2}}=\frac{1}{2}\left(\frac{V^{\prime}}{V}\right)^{2}\equiv\varepsilon_{V},$

as long as we are in slow-roll. We want to pay attention though to the
assumption in which $V\gg\dot{\phi}^{2}$, when this assumption is weaker, we
need additional terms to mitigate between $\varepsilon_{H}$ and
$\varepsilon_{V}$. In fact, later we will make the adjustments up to second
order and supply a ’recipe’ to extend this to higher orders as well. By the
same kind of process we can look at the term:

$\displaystyle\ddot{\phi}=\frac{d}{dt}\dot{\phi}=\frac{d\phi}{dt}\frac{d}{d\phi}\left(\dot{\phi}\right).$
(1.123)

Since we know that
$\dot{\phi}\simeq\frac{V^{\prime}}{3H}\simeq\frac{V^{\prime}}{\sqrt{3V}}$ we
can now derive with respect to $\phi$ and get:

$\displaystyle\ddot{\phi}\simeq\frac{\dot{\phi}}{\sqrt{3}}\left(\frac{V^{\prime\prime}}{V^{1/2}}-\frac{1}{2}\frac{V^{\prime
2}}{V^{3/2}}\right).$ (1.124)

Thus the quantity

$\displaystyle\delta_{H}\equiv\frac{\ddot{\phi}}{H\dot{\phi}}\simeq\frac{\sqrt{3}}{\sqrt{V}}\frac{1}{\sqrt{3}}\left(\frac{V^{\prime\prime}}{V^{1/2}}-\frac{1}{2}\frac{V^{\prime
2}}{V^{3/2}}\right)=\frac{V^{\prime\prime}}{V}-\frac{1}{2}\left(\frac{V^{\prime}}{V}\right)^{2},$
(1.125)

which is to say, if we define $\eta_{V}$ as

$\displaystyle\eta_{V}\equiv\frac{V^{\prime\prime}}{V},$ (1.126)

we have

$\displaystyle\delta_{H}\simeq\eta_{V}-\varepsilon_{V}\equiv\delta_{V}.$
(1.127)

####  1.5.2 First-order scalar perturbation theory

In the previous sections, we have used the Einstein Field Equations, as well
as a Stress-Energy tensor for a scalar field. The procedure for deriving the
Mukhanov-Sasaki equation uses a first-order perturbation approach to the
EFE’s. We shall go through the main parts, but only in general, whereas the
detailed process is demoted to an appendix status. It should be noted that, in
general, we follow Mukhanov’s derivation as outlined in Mukhanov1992203.

We use the conformal metric -
$g_{\mu\nu}=a^{2}(\tau)\left(\eta_{\mu\nu}+h_{\mu\nu}\right)$, where
$\eta_{\mu\nu}$ is the Minkowski metric, and $h_{\mu\nu}$ is a small
perturbation of the form:

$\displaystyle\begin{array}[]{ccc}h_{00}&=&2\Phi\\\ h_{0i}&=&0\\\
h_{ij}&=&2\Psi\cdot\delta^{i}\;_{j}\end{array}.$ (1.131)

The Action for the scalar field is given by:

$\displaystyle\begin{array}[]{ccc}I&=&\int
d^{4}x\sqrt{-g}\left[\tfrac{1}{2}g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi-V(\phi)\right]\end{array}.$
(1.133)

Varying the Action with respect to $\phi$ yields the equations of motion for
$\phi$ which are nothing but the Klein-Gordon equation for $\phi$ in conformal
time:

$\displaystyle\left(\phi^{\prime\prime}-\nabla^{2}\phi\right)+2H\phi^{\prime}+a^{2}\frac{dV}{d\phi}=0$
(1.134)

Where, using a perturbed conformal metric, as well as a first-order perturbed
$\phi=\phi^{(0)}(\tau)+\delta\phi(\tau,\vec{x})$ gives additional information:

$\displaystyle\begin{array}[]{ccc}2H\delta\phi^{\prime}+\left(\delta\phi^{\prime\prime}-\nabla^{2}\delta\phi\right)-\left(\phi^{(0)\prime}(\Phi^{\prime}+3\Psi^{\prime})\right)&+a^{2}\frac{d^{2}V}{d\phi^{2}}+2a^{2}\Phi\frac{dV}{d\phi}&=0\\\
\end{array}.$ (1.136)

The equations $G^{\mu}\;_{\nu}=8\pi GT^{\mu}\;_{\nu}$, can be broken into
perturbative degrees:

$\displaystyle\begin{array}[]{ccc}^{(0)}G^{\mu}\;_{\nu}&=&8\pi
G^{(0)}T^{\mu}\;_{\nu},\\\ \delta G^{\mu}\;_{\nu}&=&8\pi G\delta
T^{\mu}\;_{\nu},\\\ .&&\\\ \vdots&&\\\ \end{array}$ (1.141)

where the leading-order perturbations in the Stress-Energy tensor are:

$\displaystyle\left\\{\begin{array}[]{cc}\delta
T^{0}\,_{0}=&\frac{1}{a^{2}}\left[\phi^{\prime
0}\delta\phi^{\prime}-\Phi\phi^{\prime(0)2}+a^{2}\frac{dV}{d\phi}\delta\phi\right]\\\
\delta T^{0}\,_{i}=&\frac{\phi^{\prime(0)}\delta\phi_{;i}}{a^{2}}\\\ \delta
T^{i}\,_{j}=&-\frac{1}{a^{2}}\left[\phi^{\prime(0)}\delta\phi^{\prime}-\Phi\phi^{\prime(0)2}-a^{2}\frac{dV}{d\phi}\delta\phi\right]\delta^{i}\,_{j},\end{array}\right.$
(1.145)

and the corresponding perturbed Einstein tensor, after proving $\Psi=\Phi$ in
our case, is given by:

$\displaystyle\begin{array}[]{lcl}\delta
G^{0}\,_{0}&=&\frac{2}{a^{2}}\left[-3H\left(H\Phi+\Phi^{\prime}\right)+\nabla^{2}\Phi\right]\\\
\\\ \delta
G^{0}\,_{i}&=&\frac{2}{a^{2}}\partial_{i}\left[H\Phi+\Phi^{\prime}\right]\\\
\\\ \delta
G^{i}\,_{j}&=&\frac{-2}{a^{2}}\left[\left(2H^{\prime}+H^{2}\right)\Phi+3H\Phi^{\prime}+\Phi^{\prime\prime}\right]\delta^{i}\,_{j}.\end{array}$
(1.151)

Using these to compile 3 equations, and using the equations of motion for
$\phi$ we are left with:

$\displaystyle\Phi^{\prime\prime}-\nabla^{2}\Phi+2\left(H-\tfrac{2\phi^{\prime\prime}}{\phi^{\prime}}\right)\Phi^{\prime}+2\left(H\\-\tfrac{2\phi^{\prime\prime}}{\phi^{\prime}}H\right)\Phi=0,$
(1.152)

Or in an equivalent form:

$\displaystyle\Phi^{\prime\prime}-\nabla^{2}\Phi+2\left(\tfrac{a}{\phi^{\prime}}\right)^{\prime}\left(\tfrac{a}{\phi^{\prime}}\right)^{-1}\Phi^{\prime}+2\phi^{\prime}\left(\tfrac{H}{\phi^{\prime}}\right)^{\prime}\Phi=0.$
(1.153)

Switching to a new variable
$u\equiv\left(\tfrac{a}{\phi^{\prime}}\right)\Phi$, with
$\theta\equiv\frac{H}{a\phi^{\prime}}$ this equation now takes the form:

$\displaystyle
u^{\prime\prime}-\nabla^{2}u-\left(\frac{\theta^{\prime\prime}}{\theta}\right)u=0,$
(1.154)

This is the Mukhanov-Sasaki equation.

To simplify the solution of this equation we decompose $u$ unto Fourier modes
$u_{k}$, going into Fourier space we get:

$\displaystyle
u^{\prime\prime}_{k}+\left(k^{2}-\frac{z^{\prime\prime}}{z}\right)u_{k}=0\;,$
(1.155)

and when we go to Fourier space the pump field in the original space,
$\theta$, becomes the pump field in Fourier space

$\displaystyle z=\frac{a\dot{\phi}}{H}.$ (1.156)

We can wrap up the quantity $\left(k^{2}-\frac{z^{\prime\prime}}{z}\right)$ as
follows:

$\displaystyle\left(k^{2}-\frac{z^{\prime\prime}}{z}\right)\equiv\omega_{k}^{2}(t)\;,$
(1.157)

so the MS equation can be expressed as:

$\displaystyle u^{\prime\prime}_{k}+\omega_{k}^{2}(t)u_{k}=0\;.$ (1.158)

In this form it is easily recognized as a Time-Dependent Harmonic Oscillator
(TDHO).

###  1.6 The primordial power spectrum - a signature of our history

The patterns we currently see in the CMB, are remnants of the time of last
scattering. The CMB photons that are just now detected in our instruments were
last scattered when the universe was approximately $1000\sim 3000$ times
smaller and hotter. At that time the different constituents of the universe
were in thermal equilibrium. Photons were scattered off electrons, and in
turn, electrons were constantly bombarded with photons. In such a state matter
is in a plasma phase, the optical depth is very small and whatever information
a photon carries is erased after a few interactions. The current working
assumption is that the optical depth increased suddenly, due to a sudden
temperature drop that enabled electrons to be captured into previously ionized
hydrogen atoms. At that point, the photons were no longer tightly coupled to
electrons and were free to escape. It is these photons that first broke free
that we now collect in our CMB instruments.

####  1.6.1 Matter power spectrum

However, since the time of last scattering several mechanisms have affected
photons over time. It is customary to decompose the CMB radiation to modes,
each corresponding to a different physical scale. The smallest $k$-mode we can
hope to detect now corresponds to the scale of the current universe. The mode
that is just now becoming causally connected corresponds to the current light
horizon size. It can be shown that any region separated by more than 2 degrees
in the sky today would have been causally disconnected at the time of
decoupling Liddle:1999mq. Modes that are not yet causally connected are
’frozen’ and are roughly measured "as-is". Conversely, causally connected
modes can now discharge energy, dissipate, or otherwise interact with other
energy constituents of the universe. As such the CMB we currently detect
differs from the radiation at the time of last scattering.

####  1.6.2 Transfer function

In order to reflect the different physical processes sub-horizon and super-
horizon modes undergo, as well as the overall expansion of the universe, which
affect all modes, the Transfer function is called for. Given an initial power
spectrum at the time of last scattering:

$\displaystyle P(k)_{\text{matter}}\propto
P(k)_{prim}T^{2}(k)\left(\frac{D(a)}{D(a=1)}\right)^{2},$ (1.159)

where here $P(k)_{\text{matter}}$ is the matter power spectrum, $P(k)_{prim}$
is the PPS, $T(k)$ is the transfer function, and $D(a)$ is the growth
function, which is dependent only on $a$, and reflects the change in modes due
to the overall increase in the scale factor of the universe.

####  1.6.3 Scalar power spectrum

The primordial power spectrum (PPS) is traditionally characterized by its
spectral index $n_{s}$ and the index running $\alpha_{s}$, which are given by
the first and second logarithmic derivatives of the logarithm of the PPS:

$\displaystyle
n_{s}=1+\left.\frac{\partial\ln\left(P_{s}\right)}{\partial\ln\left(k\right)}\right|_{aH=k},$
(1.160)
$\displaystyle\alpha_{s}=\left.\frac{\partial^{2}\ln\left(P_{s}\right)}{\partial\ln\left(k\right)^{2}}\right|_{aH=k}=\left.\frac{\partial
n_{s}}{\partial\ln\left(k\right)}\right|_{aH=k},$ (1.161)

where $aH=k$ denotes the CMB scale. Extending this Taylor expansion, we can
write the PPS as:

$\displaystyle
P_{s}=A_{s}\left(\frac{k}{k_{0}}\right)^{n_{s}-1+\frac{\alpha_{s}}{2}\ln{\left(\frac{k}{k_{0}}\right)}+\frac{\beta_{s}}{6}\ln^{2}{\left(\frac{k}{k_{0}}\right)}+...},$
(1.162)

where $k_{0}$ is the scale around which we Taylor expand, and is called the
pivot scale. $\beta_{s}$ is the next term in the Taylor expansion and is
called the running of running. It is given by:

$\displaystyle\beta_{s}=\left.\frac{\partial\alpha_{s}}{\partial\ln{k}}\right|_{aH=k}\;.$
(1.163)

In most places the subscript denoting the assessment of the power spectrum at
$aH=k$ is suppressed.

###  1.7 Stewart-Lyth power spectrum

####  1.7.1 Derivation

Here we retrace the procedure of deriving the most widely used analytical
expression for $n_{s}$. Recalling the definition for the pump field $z$, and
the MS equation (Eqs. (1.156),(1.155)). The parameter $\nu$ is properly
defined as:

$\displaystyle\nu=+\sqrt{\frac{z^{\prime\prime}}{z}\tau^{2}+\frac{1}{4}}.$
(1.164)

However, in the Stewart-Lyth (SL) formulation, the approximations made lead to
the definition of $\nu$ as the somewhat different:

$\displaystyle\nu_{\text{\tiny{SL}}}=\frac{3+2\delta_{H}+\epsilon_{H}}{2\left(1-\epsilon_{H}\right)}\;.$
(1.165)

Then the MS equation becomes:

$\displaystyle
u_{k}^{\prime\prime}+\left(k^{2}-\frac{\left(\nu^{2}-\frac{1}{4}\right)}{\tau^{2}}\right)u_{k}=0.$
(1.166)

For a constant $\nu$ this becomes the Bessel equation, with known solutions.
The value of $\tfrac{z^{\prime\prime}}{z}$ is formally given by:

$\displaystyle\frac{z^{\prime\prime}}{z}=2a^{2}H^{2}\left(1+\frac{3\delta_{H}}{2}+\epsilon_{H}+\frac{\delta_{H}^{2}}{2}+\frac{\epsilon_{H}\delta_{H}}{2}+\frac{1}{2H}\left(\dot{\epsilon_{H}}+\dot{\delta_{H}}\right)\right).$
(1.167)

In many cases, one assumes that the time derivatives are small and can be
neglected. However, these derivatives yield second-order terms that can
significantly affect the value of $\tfrac{z^{\prime\prime}}{z}$. The full
expression is given by:

$\displaystyle\frac{z^{\prime\prime}}{z}=2a^{2}H^{2}\left(1+\frac{3\delta_{H}}{2}+\epsilon_{H}+\epsilon_{H}^{2}+2\epsilon_{H}\delta_{H}+\frac{1}{2}\frac{\delta_{H}\dddot{\phi}}{H\ddot{\phi}}\right),$
(1.168)

which may differ from Eq. (1.167) after removing the time derivative terms,
when $\delta^{2}_{H}$ and/or $\frac{\delta_{H}\ddot{\phi}}{H\dot{\phi}}$ are
non-negligible. $\epsilon^{2}_{H}$ is usually of the order of $10^{-5}$ or
less, even for models with high $r$.

Applying boundary conditions and taking the small arguments limit we are left
with a power spectrum of:

$\displaystyle\ln\left(P_{R}\right)$
$\displaystyle=-\ln(32\pi^{2}\Gamma^{2}(\tfrac{3}{2}))$ (1.169)
$\displaystyle+2\nu\ln(2)+2\ln(k)+2\ln\left(\Gamma(\nu)\right)$
$\displaystyle+(1-2\nu)\ln(-k\tau),$

which yields the scalar index of:

$\displaystyle
n_{s}=4-2\nu+2\left(\ln(2)+\psi(\nu)\right)\frac{\partial\nu}{\partial\ln(k)},$
(1.170)

with the digamma function $\psi(x)\equiv\frac{\Gamma^{\prime}(x)}{\Gamma(x)}$.
The final expression is heavily dependent on the value and time derivative of
$\nu$. This is a possible source of inaccuracy. In the original text the
author now define:

$\displaystyle\begin{array}[]{ccc}\alpha=\left(\frac{V_{,\phi}}{V}\right)^{2}&\beta=\frac{V_{,\phi\phi}}{V}&\gamma=\frac{V_{,\phi^{3}}}{V_{,\phi}},\\\
\end{array}.$ (1.172)

Having defined these, usually one connects the original slow-roll parameters
with the above quantities by Stewart:1993bc

$\displaystyle\begin{array}[]{c}\epsilon_{H}\simeq\frac{\alpha}{2}-\frac{\alpha^{2}}{3}+\frac{\alpha\beta}{3}\\\
\\\
\delta_{H}\simeq\frac{\alpha}{2}-\beta-\frac{2\alpha^{2}}{3}+\frac{4\alpha\beta}{3}-\frac{\beta^{2}}{3}-\frac{\alpha\gamma}{3}\\\
\\\
\frac{\delta_{H}\dddot{\phi}}{H\ddot{\phi}}\simeq\alpha^{2}-\frac{5\alpha\beta}{2}+\beta^{2}+\alpha\gamma.\end{array}$
(1.178)

With these relations, one can substitute the slow-roll parameters in Eq.
(1.170), for the quantities in Eq. (1.178), to get the most commonly used
analytical expression for the scalar index Lyth:1998xn:

$\displaystyle n_{s}$ $\displaystyle\simeq 1-6\varepsilon_{V}+2\eta_{V}$
$\displaystyle+2\times\bigg{[}\frac{\eta_{V}^{2}}{3}-(8b+1)\varepsilon_{V}\eta_{V}$
$\displaystyle-\left(\frac{5}{3}-12b\right)\varepsilon_{V}^{2}+\left(b+\frac{1}{3}\right)\xi_{V}^{2}\bigg{]},$
(1.179)

where
$\varepsilon_{V}=\tfrac{\alpha}{2}\;;\;\eta_{V}=\beta\;;\;\xi^{2}_{V}=\alpha\gamma$,
and $b=2-\ln{2}-\gamma_{\mathrm{\scalebox{0.4}{Euler}}}$,
$\gamma_{\mathrm{\scalebox{0.4}{Euler}}}$ being the Euler number.

####  1.7.2 Analytical term for index running $\alpha_{s}$

Following the same procedure for the running of the scalar index, gives, to
second-order:

$\displaystyle\alpha_{s}\simeq$
$\displaystyle-16\varepsilon_{V}\eta_{V}+24\varepsilon_{V}^{2}+2\xi_{V}^{2},$
(1.180)

where both (1.179), and (1.180), are evaluated at the CMB point if one wishes
to probe the CMB point. Equivalently if one wishes to compare to observables
as derived from MCMC analyses at some pivot scale $k_{\mathrm{pivot}}$, one
needs to find the correct $\phi_{\mathrm{pivot}}$, and evaluate these
expressions at the pivot scale.

###  1.8 Other formulations for the power spectrum observables

In this section, attention should be paid to different formulations and
notations. Every subsection is self-contained to encapsulate the derivations
fully. Wherever possible we use the notation in the source material. Over the
years, our observations have become more accurate and spanned a longer
interval of physical scales Bennett:2003ba; Leitch:2001yx; Halverson:2001yy;
Pryke:2001yz; Netterfield:2001yq; deBernardis:2001mjh; Lee:2001yp;
Stompor:2001xf; Tegmark:2000qy; Hannestad:2001nu. Even well before the Planck
mission Planck:2006aa was proposed and eventually launched, the prospects of
probing not just the slope of the PPS but the running of the slope was
compelling. As a result, several forays into calculating the running of the
spectral index were made. It was understood early on, that Eq. (1.180), was
put forth as a second order derivation of an expression that by definition
vanishes. As a result different approaches were taken, at the level of
analytic expression for the PPS itself, rather than calculating derivatives of
the $n_{s}$ it predicts. All of these attempts were made around 2011.

####  1.8.1 Dodelson-Stewart

Twin papers written by Ewan Stewart and Scott Dodelson were published in
September and October of 2011. Stewart published a detailed calculation of the
PPS Stewart:2001cd, dropping a previous assumption, that
$|\alpha_{s}|\ll|n_{s}-1|$. Dodelson and Stewart published a joint paper,
based on the calculations made in Stewart:2001cd, that studies the
ramifications of these formulae on inflationary potentials Dodelson:2001sh.
This treatment uses the slow-roll parameters in the efold-flow formulation,
such that:

$\displaystyle\varepsilon=\frac{1}{2}\left(\frac{\dot{\phi}}{H}\right)\hskip
20.0pt\textrm{and}\hskip
20.0pt\delta_{p}\equiv\frac{1}{H^{p}\dot{\phi}}\left(\frac{d}{dt}\right)^{p}\dot{\phi},$
(1.181)

and their evolution is governed by:

$\displaystyle\frac{d\varepsilon}{d\ln{a}}=\frac{d\varepsilon}{dN}=2\left(\varepsilon+\delta_{1}\right)\varepsilon$
(1.182)
$\displaystyle\frac{d\delta_{p}}{d\ln{a}}=\frac{d\delta_{p}}{dN}=\delta_{p+1}+\left(p\varepsilon-\delta_{1}\right)\delta_{p}.$
(1.183)

In this approximation the scalar index is given by:

$\displaystyle
n_{s,\mathrm{\scalebox{0.5}{DS}}}=1-4\varepsilon-2\delta_{1}-2\sum_{p=1}^{\infty}d_{p}\delta_{p+1}+\mathcal{O}(\xi^{2}),$
(1.184)

where $d_{p}$ are numerical coefficients of order unity, and the remainder
includes the terms $\varepsilon^{2},\delta_{1}\varepsilon$ and
$\delta_{1}^{2}$. These $d_{p}$ terms can be calculated by using the following
generating function:

$\displaystyle\sum_{p=0}^{\infty}d_{p}x^{p}=2^{x}\cos{\left(\frac{\pi
x}{2}\right)}\frac{\Gamma(2-x)}{1+x}.$ (1.185)

In this formulation the scalar index running is given by:

$\displaystyle\alpha_{s,\mathrm{\scalebox{0.5}{DS}}}=-2\sum_{p=0}^{\infty}d_{p}\delta_{p+2}+\mathcal{O}(\xi^{2}).$
(1.186)

In terms of inflationary potential and derivatives thereof, the scalar index
and the index running are given by:

$\displaystyle
n_{s,\mathrm{\scalebox{0.5}{DS}}}=-3\left(\frac{V^{\prime}}{V}\right)^{2}+2\sum_{p=0}^{\infty}q_{p}\left(\frac{V^{\prime}}{V}\right)^{p}\frac{V^{(p+2)}}{V}$
(1.187)
$\displaystyle\alpha_{s,\mathrm{\scalebox{0.5}{DS}}}=-2\sum_{p=0}^{\infty}q_{p}\left(\frac{V^{\prime}}{V}\right)^{p+1}\frac{V^{(p+3)}}{V},$
(1.188)

where here $q_{p}$ are generated by the following generating function:

$\displaystyle\sum_{p=0}^{\infty}q_{p}x^{p}\equiv 2^{-x}\cos{\left(\frac{\pi
x}{2}\right)}\frac{3\Gamma(2+x)}{(1-x)(3-x)}.$ (1.189)

As we shall see this formulation is close, but only in the ’horseshoes and
hand grenades’ sense.

####  1.8.2 Stewart-Gong

The Stewart-Gong formulation Gong:2001he is, in essence, similar. They first
set up the Green’s function formalism to tackle the problem of inflation. They
rewrite the MS equation Eq. (1.155) as:

$\displaystyle\frac{d^{2}y}{dx^{2}}+\left(1-\frac{1}{z}\frac{d^{2}z}{dx^{2}}\right)y=0,$
(1.190)

where $x$ and $y$ are defined as111Note that in this derivation $\eta$ denotes
conformal time:

$\displaystyle x=-k\eta$ (1.191) $\displaystyle y=\sqrt{2k}u_{k}.$ (1.192)

They then choose the ansatz $z=\frac{1}{x}f(\ln{x})$, where $f$ is some
function. Defining a new function $g$ such that:

$\displaystyle g=\frac{-3f^{\prime}+f^{\prime\prime}}{f},$ (1.193)

the equation of motion now transforms to:

$\displaystyle\frac{d^{2}y}{dx^{2}}+\left(1-\frac{2}{x^{2}}\right)y=\frac{1}{x^{2}}g(\ln{x})y.$
(1.194)

Now we can take the Green’s function approach to get the eigenfunction
$y_{0}=\left(1+\frac{i}{x}\right)e^{ix}$, and yield the full solution:

$\displaystyle
y(x)=y_{0}(x)+\frac{i}{2}\int_{x}^{\infty}du\frac{1}{u^{2}}g(\ln{u})y(u)\left[y_{0}^{*}(u)y_{0}(x)-y_{0}^{*}(x)y_{0}(u)\right].$
(1.195)

The authors continue to write the full expression for $n_{s}$:

$\displaystyle n_{s,\mathrm{\scalebox{0.5}{SG}}}$
$\displaystyle=1-4\varepsilon_{1}-2\delta_{1}+(8b-8)\varepsilon_{1}^{2}+(10b-6)\varepsilon_{1}\delta_{1}-2b\delta_{1}^{2}+2b\delta_{2}$
(1.196)
$\displaystyle+\left(-16b^{2}+40b-108+\frac{28\pi^{2}}{3}\right)\varepsilon_{1}^{3}+\left(-31b^{2}+60b-172+\frac{199\pi^{2}}{12}\right)\varepsilon_{1}^{2}\delta_{1}$
$\displaystyle+\left(-3b^{2}+4b-30+\frac{13\pi^{2}}{4}\right)\varepsilon_{1}\delta_{1}^{2}+\left(-2b^{2}+8-\frac{5\pi^{2}}{6}\right)\delta_{1}^{3}$
$\displaystyle+\left(-7b^{2}+8b-22+\frac{31\pi^{2}}{12}\right)\varepsilon_{1}\delta_{2}+\left(3b^{2}-8+\frac{3\pi^{2}}{4}\right)\delta_{1}\delta_{2}+\left(-b^{2}+\frac{\pi^{2}}{12}\right)\delta_{3},$

where the hierarchy of slow-roll parameters is defined as:

$\displaystyle\varepsilon_{1}=\frac{-\dot{H}}{H^{2}}\;;\;\delta_{n}=\frac{1}{H^{n}\dot{\phi}}\frac{d^{n+1}\phi}{dt^{n+1}}.$
(1.197)

At the end of this formulation, as we shall see, lay no cigars as well.
It is, however, worth mentioning, that Dvorkin & Hu Dvorkin:2009ne, used this
method iteratively to calculate the PPS to $\sim 0.3\%$ accuracy.

####  1.8.3 Schwarz-Escalante-Garcia

The formulation in Schwarz:2001vv claims not to rely on the slow-roll
expansion, but as we shall see this is not accurate. We first define the
horizon flow parameters as:

$\displaystyle\varepsilon_{0}=\frac{d_{H}(N)}{d_{H_{i}}},$ (1.198)

where $d_{H}\equiv\frac{1}{H}$ is the Hubble distance, and $d_{H_{i}}$ is the
Hubble distance at the start of inflation (at $t_{i}$, where $N(t_{i})=0$).
The hierarchy of functions is defined as:

$\displaystyle\varepsilon_{m+1}\equiv\frac{d\ln{|\varepsilon_{m}|}}{dN},\hskip
20.0pt\;m>0.$ (1.199)

According to this definition, $\varepsilon_{1}=\frac{d\ln{d_{H}}}{dN}$, and it
can be shown that

$\displaystyle\varepsilon_{1}\equiv\varepsilon_{H}=-\frac{\dot{H}}{H^{2}}.$
(1.200)

Upon closer inspection, one finds this is nothing but a rewriting of the slow-
roll parameters of Eq. (1.197), as shown in table .1.1.

Schwarz-Escalante-Garcia | Stewart-Gong
---|---
$\varepsilon_{1}$ | $\varepsilon_{1}$
$\varepsilon_{2}$ | $2\varepsilon_{1}+2\delta_{1}$
$\varepsilon_{2}\varepsilon_{3}$ | $4\varepsilon_{1}^{2}+6\varepsilon_{1}\delta_{1}-2\delta_{1}^{2}+2\delta_{2}$
Table 1.1: Slow-roll parameters in Schwarz-Escalante-Garcia formulation,
translated to Stewart-Gong formalism.

As such their analysis is not so different from previous ones made. However,
they provide a taxonomy of models which divide models into two categories:
Constant-Horizon (CH) models, and Growing-Horizon (GH) models. The CH category
is reserved for models that have a very small time derivative of the Hubble
distance. This implies that $\varepsilon_{1}\ll 1$, but it doesn’t necessarily
mean that $\varepsilon_{m}$ is small for arbitrary $m$. Thus the authors
define $\mathrm{CH}n$ as models in which
$|\varepsilon^{n}_{2}|<\mathrm{max}\left(|\varepsilon_{1}\varepsilon_{2}|,|\varepsilon_{2}\varepsilon_{3}|\right)$.
In this case, the following terms are allowed in the approximation:
$1,\varepsilon_{1},\varepsilon_{2},...,\varepsilon_{2}^{n}$. In the GH
treatment, the approximation scheme is different. This scheme is inspired by
the power-law inflation case in which $\varepsilon_{1}=\frac{1}{p}$, for
inflation with $a\propto t^{p}$, and the other slow-roll parameters strictly
vanish. Thus this scheme defines $n$ as the maximal integer in which
$\varepsilon_{1}^{n}>\mathrm{max}\left(|\varepsilon_{1}\varepsilon_{2}|,|\varepsilon_{2}\varepsilon_{3}|\right)$
holds true. So the slow-roll terms allowed in this approximation scheme are
$1,\varepsilon_{1},...,\varepsilon_{1}^{n},\varepsilon_{2}$.
The power spectra given in the CH case are given by:

$\displaystyle k^{3}P_{s}$ $\displaystyle=$
$\displaystyle\frac{H^{2}}{\pi\varepsilon_{1}}\big{[}a_{0}+a_{1}\ln{(k/k_{0})}+a_{2}\ln^{2}{(k/k_{0})}+a_{3}\ln^{3}{(k/k_{0})}$
$\displaystyle+...\big{]},$ (1.201) $\displaystyle a_{0}$ $\displaystyle=$
$\displaystyle
1-2(1-b)\varepsilon_{1}+b\varepsilon_{2}+\frac{1}{8}\left(4b^{2}+\pi^{2}-8\right)\varepsilon_{2}^{2}$
(1.202)
$\displaystyle-\frac{1}{24}\left(3b(8-\pi^{2})-4b^{3}+14\zeta(3)-16\right)\varepsilon_{2}^{3},$
$\displaystyle a_{1}$ $\displaystyle=$
$\displaystyle-2\varepsilon_{1}-\varepsilon_{2}-b\varepsilon_{2}^{2}-\frac{1}{8}\left(4b^{2}+\pi^{2}-8\right)\varepsilon_{2}^{3},$
(1.203) $\displaystyle a_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\varepsilon_{2}^{2}+\frac{b}{2}\varepsilon_{2}^{3},$
(1.204) $\displaystyle a_{3}$ $\displaystyle=$
$\displaystyle-\frac{1}{6}\varepsilon_{2}^{3},$ (1.205)

with the same $b$ as before:
$b=2-\ln{2}-\gamma_{\mathrm{\scalebox{0.4}{Euler}}}$, and where
$\zeta(3)\approx 1.2021$. The power spectra of GH models are given by:

$\displaystyle k^{3}P_{s}$ $\displaystyle=$
$\displaystyle\frac{H^{2}}{\pi\varepsilon_{1}}\big{[}c_{0}+c_{1}\ln{(k/k_{0})}+c_{2}\ln^{2}{(k/k_{0})}+c_{3}\ln^{3}{(k/k_{0})}$
$\displaystyle+...\big{]},$ (1.206) $\displaystyle c_{0}$ $\displaystyle=$
$\displaystyle
1-2(1-b)\varepsilon_{1}+\frac{1}{2}\big{[}4b(b-1)+\pi^{2}-10\big{]}\varepsilon_{1}^{2}$
(1.207)
$\displaystyle+\frac{1}{3}\big{[}4b^{3}+3b(\pi^{2}-12)-14\zeta(3)+19\big{]}\varepsilon_{1}^{3}+b\varepsilon_{2},$
$\displaystyle c_{1}$ $\displaystyle=$
$\displaystyle-2\varepsilon_{1}+2(1-2b)\varepsilon_{1}^{2}-(4b^{2}+\pi^{2}-12)\varepsilon_{1}^{3}-\varepsilon_{2},$
(1.208) $\displaystyle c_{2}$ $\displaystyle=$ $\displaystyle
2\varepsilon_{1}^{2}+4b\varepsilon_{1}^{3},$ (1.209) $\displaystyle c_{3}$
$\displaystyle=$ $\displaystyle-\frac{4}{3}\varepsilon_{1}^{3}$ (1.210)

with the same values of $b$ and $\zeta(3)$.
While this approach at least gives us consistent conditions for the voracity
of different approximation schemes, it is no closer to the mark than the
previous ones.

###  1.9 The Lyth bound, original and extended versions

One of the most useful heuristics out there, to build and evaluate
inflationary models, is the so-called Lyth bound. In the original paper
Lyth:1996im, there are several conclusions. However we concentrate on the
following argument: Given that the PPS can be written as

$\displaystyle P_{s}\simeq\left.\frac{8\pi
GH^{2}}{9\varepsilon_{H}k^{3}}\right|_{k=aH},$ (1.211)

and the tensor power spectrum is given by:

$\displaystyle P_{T}\simeq\left.\frac{8\pi GH^{2}}{k^{3}}\right|_{k=aH},$
(1.212)

we can define the tensor-to-scalar ratio, $r$, as the power spectra ratio at
some scale $k_{0}$ of our choosing. Thus $r$ is set as:

$\displaystyle r=9\varepsilon_{H}.$ (1.213)

The numerical coefficient is subject to some finer detail, but the order of
magnitute is always the same. In Lyth:1996im we have $r=2\cdot
6.9\varepsilon_{h}$, whereas other sources state $r=12\sim 16\varepsilon_{H}$.
This statement, along with the (extreme) slow-roll assumption in which
$\varepsilon_{H}=\varepsilon_{V}$ yields the following:

$\displaystyle\left(\frac{r}{6.9}\right)^{\frac{1}{2}}=\left|\frac{V^{\prime}}{V}\right|.$
(1.214)

In order to complete the argument, we must look at the quantity

$\displaystyle\frac{d\phi}{dN}=\frac{d\phi}{Hdt}=\frac{\dot{\phi}}{H}=\sqrt{2\varepsilon_{H}}.$
(1.215)

So, plugging in the former equation we get:

$\displaystyle\frac{d\phi}{dN}=\left(\frac{r}{6.9}\right)^{\frac{1}{2}},$
(1.216)

where we adhere to the numerical coefficient in Lyth:1996im.
The scales on which the tensor power spectrum is most prominent, and therefore
is more likely to be discovered at is $1<l\lesssim 100$. With
$k\sim\frac{H_{0}l}{2}$, the interval in $k$-space over which we expect to
find tensor modes is $\Delta\ln{k}\simeq 4.6$. Taking $H$ to be roughly
constant, when we evaluate $k=aH$ we have

$\displaystyle\ln{k}=\ln{aH}=\ln{H}+\ln{a}=\ln{H}+N,$ (1.217)

thus

$\displaystyle d\ln{k}=dN.$

If we additionally take $d\phi/dN$ to be approximately constant along the
aforementioned $k$-interval we have:

$\displaystyle\frac{d\phi}{dN}\simeq\frac{\Delta\phi}{\Delta N}.$ (1.218)

Putting these together yields:

$\displaystyle\Delta\phi\simeq
4.6\left(r/6.9\right)^{1/2}=0.46\left(r/0.07\right)^{1/2}.$ (1.219)

This means that given some $r$ we can approximate the field excursion, in
Planck units, during the first $\sim 4.6$ efolds of inflation. For example,
suppose we find $r\lesssim 0.03$, which is possible in the foreseeable future.
In that scenario we have:

$\displaystyle\Delta\phi\sim 0.3\;.$ (1.220)

Applying Eq. (1.214), we can also get a one-to-one relation between $r(k)$ and
the slope of the potential at a chosen $k_{0}$, as long as $k_{0}$ is inside
the $1<l\lesssim 100$ window. This relation is formalized in BenDayan:2009kv
as

$\displaystyle\frac{V^{\prime}}{V}=-\sqrt{\frac{r}{8}},$ (1.221)

where the numerical coefficient is subject to slight changes, and the minus
sign comes from the notion of the field rolling towards the larger positive
$\phi$ values.

This original argument was extended in Easther:2006qu; Efstathiou:2005tq. In
Easther:2006qu, the argument was amended to include higher slow-roll
parameters, to allow for models that change by a significant amount during the
tensor mode window. This extended Lyth bound is given by:

$\displaystyle\Delta\phi\gtrsim\sqrt{\frac{r}{4\pi}}\left[2-n_{s}-\frac{r}{8}\right],$
(1.222)

and in the case of $r=0.03$ we have $\Delta\phi\gtrsim 0.05$. In
Efstathiou:2005tq however, some assumptions were added, in order to extend the
bound to the full efold extent of inflation, i.e. $N=55\sim 60$. In this case
the authors extend the Lyth bound such that:

$\displaystyle\Delta\phi\approx 6r^{1/4}.$ (1.223)

This would imply that models that predict $r=0.03$, should incur a field
excursion of $\Delta\phi\sim 2.5$.

###  1.10 Small and Large field models taxonomy

Current nomenclature distinguishes between ’Large field models’ and ’Small
field models’. This distinction is based on Dodelson:1997hr, where several
models were studied and a taxonomy of models was presented. In the interest of
brevity and accessibility, we will first reiterate the key concepts needed to
understand this distinction. First, the slow-roll condition on the potential
is given by:

$\displaystyle\varepsilon_{V}\equiv\frac{1}{2}\left(\frac{V^{\prime}}{V}\right)^{2}<1,$
(1.224)

where $V^{\prime}=\frac{dV(\phi)}{d\phi}$, and usually it is understood as a
strong requirement i.e. we prefer $\varepsilon_{V}\ll 1$, in order to be in
the slow-roll regime. The number of efolds generated within a certain
inflation is given by:

$\displaystyle
N=\int_{t_{\mathrm{init}}}^{t_{\mathrm{fin}}}Hdt\simeq\int_{\phi_{\mathrm{init}}}^{\phi_{\mathrm{fin}}}-\frac{V}{V^{\prime}}d\phi,$
(1.225)

where the minus sign is due to our interpretation of the field $\phi$ rolling
down the potential. Finding $\phi_{\mathrm{fin}}$ given the potential is
usually straightforward, by setting $\varepsilon_{V}=1$ we identify the end of
slow-roll inflation, thus:

$\displaystyle\left(\frac{V^{\prime}}{V}\right)=\pm\sqrt{2},$ (1.226)

from which we can extract $\phi_{\mathrm{fin}}$.
In Dodelson:1997hr, the authors divide the models studied into 3 categories:

  * 1\.

Large field models: in which the inflaton field $\phi$ is displaced far from
its minimum, and rolls down the potential towards a minimum at the origin. The
models $V(\phi)=\Lambda^{4}\left(\phi^{p}\right)\;;\;p>1$ and
$V(\phi)=\Lambda^{4}\exp\left(-\sqrt{2\phi^{2}/p}\right)\;;\;p>0$ are classic
examples of large field models.

  * 2\.

Small field models: in which the field is initially near the origin, and rolls
down towards a minima that is removed from the origin ($<\phi>\neq 0$), these
types of models would be expected as a result of spontaneous symmetry
breaking.

  * 3\.

Hybrid models of inflation: in which the field evolves towards a minima with a
non-zero vacuum energy. Usually, hybrid models are realized with a number of
dynamic fields. However during inflation, in most cases one can identify a
dominant field, and thus treat this type of models as single field inflation
models, at least a posteriori. We do not discuss hybrid models of inflation
within the context of this work.

####  1.10.1 Large field models

Since these type of models, as defined in Dodelson:1997hr, usually have
initial $\phi$ values of $\sim 10\;m_{pl}$ or more, the taxonomy had evolved
to the following simplified meaning. Large field models are models in which
the field excursion during inflation is more than a few $m_{pl}$. Let us
examine the class of Large field polynomial potentials to understand why that
is:

$\displaystyle V(\phi)=\Lambda^{4}\left(\phi^{p}\right)\;;\;p>1.$ (1.227)

By virtue of Eq. (1.226) we have $\phi_{\mathrm{fin}}=\pm\frac{p}{\sqrt{2}}$,
and by applying Eq. (1.225), we recover the following condition on the initial
$\phi$ value:

$\displaystyle\sqrt{p\left(2N+\frac{1}{\sqrt{2}}\right)}=\phi_{\mathrm{init}}.$
(1.228)

With $N\gtrsim 50$, and even when $p=1$, this produces
$\phi_{\mathrm{init}}\sim 10\;m_{pl}$ and $\Delta\phi$ approximately the same.
Thus in these models $\Delta\phi\gtrsim 10\;m_{pl}$. With the exponential
models the end of inflation occurs near $\phi_{\mathrm{fin}}=\sqrt{p}$. The
relation between the starting and end positions of the field is given by:

$\displaystyle\sqrt{p}e^{\frac{2N}{p}}=\phi_{\mathrm{init}},$ (1.229)

and a quick inspection of this as a function of $p$ reveals a minimum at
$p=4N$, which means the minimum field excursion is given by:

$\displaystyle\Delta\phi\gtrsim 2\sqrt{N}\left(e^{\frac{1}{2}}-1\right)\approx
9\;m_{pl}$ (1.230)

when $N=50$. This is enough to justify the notion of large field models, as
models in which the field excursion during inflation is of several Planck
masses.

####  1.10.2 Small field models

The methodology of evaluating small field models is the same as the above,
with the only difference of $\phi_{\mathrm{init}}$ starting near the origin:

$\displaystyle<\phi_{\mathrm{init}}>\simeq 0.$ (1.231)

The usual representative of this class of models is the "small-field
polynomial model"

$\displaystyle V(\phi)=\Lambda^{4}\left(1-(\phi/\mu)^{p}\right)\;\phi\ll\mu\ll
m_{pl}\;\text{and}\;p>2.$ (1.232)

This model is not to be confused with the models we study, as it is actually a
monomial, meaning there is only one term in which $\phi$ appears. Note that
since this is the most usual representative of the small field models, there
is a tendency to study this class and apply these conclusions to the entire
class. This is evident in Martin:2013tda for example. In these types of models
it is easily discerned that the field excursion can be of the order of
$\Delta\phi\lesssim 1\sim 2\;m_{pl}$, and usually much smaller.

A different more observable-oriented classification of models can be found in
Schwarz:2004tz.

##  Chapter 2 Our Models

Small field models of inflation in which inflation occurs near a flat feature,
a maxima, or a saddle point are studied (see Boubekeur:2005zm for a review).
This class of models is interesting because they appear in many fundamental
physics frameworks, effective field theory, supergravity Yamaguchi:2011kg and
string theory Baumann:2014nda in successive order of complexity. Our focus on
such models is also motivated by the expected properties of the moduli
potentials in string theory. More generally speaking these type of models can
be viewed as a Taylor expansion approach to other models Dodelson:1997hr. A
different more observable-oriented classification of models can be found in
Schwarz:2004tz, in which analysis our models fall into the toward-exit class.

In general, inflation will occur in a multi-dimensional space. However, the
results for multifield inflation cannot usually be obtained simply. In many
known cases, it is possible to identify a-posteriori a single degree of
freedom along which inflation takes place. To gain some insight about the
expected typical results effective single field potentials can be used.

Generic small field models predict a red spectrum of scalar perturbations,
negligible spectral index running and non-gaussianity. They also predict a
characteristic suppression of tensor perturbations BenDayan:2008dv. Hence,
they were not viewed as candidate models for high-$r$ inflation. Large field
models of inflation are thus the standard candidates for high-$r$ inflation.
For more detailed model building considerations one can review Makarov:2005uh
and Lesgourgues:1998mq.

In BenDayan:2009kv, a new class of more complicated single small field models
of inflation was considered (see also Hotchkiss:2011gz) that can predict,
contrary to popular wisdom Lyth:1996im; Martin:2013tda, an observable GW
signal in the CMB (see also Cicoli:2008gp.) The notion that observable signal
GW precludes small field models partly stems from Martin:2013tda and similar
analyses that study monomial potential models as small field models. The
spectral index, it’s running, the tensor to scalar ratio and the number of
e-folds were claimed to cover all the parameter space currently allowed by
cosmological observations. The main feature of these models is that the high
value of $r$ is accompanied by a relatively strong scale dependence of the
resulting power spectrum. Another unique feature of models in this class is
their ability to predict, again contrary to popular wisdom Easther:2006tv, a
negative spectral index running. The single observable consequence that seems
common to all single field models is the negligible amount of non-gaussianity.
In Lesgourgues:2007gp the inflationary potential was Taylor-expanded up to
order $4$. The approach applied in Lesgourgues:2007gp is similar to ours,
however only potentials that are monotonic in the entire CMB window were
considered. The family of potentials we study can easily contain members that
have a shallow minimum point followed by an equally shallow maximum. As long
as there is enough kinetic energy to clear that interval in $\phi$, we will
clear the potential ‘valley’ and not be trapped in a scenario where inflation
doesn’t end. A concrete example would be a third degree small field monomial
of the form $V(\phi)=V_{0}(1-M\phi^{3})$, after a shift such that
$\phi\rightarrow\phi-\tilde{\phi}$. In this case $\phi$ and $\phi^{2}$ terms
will appear in the shifted potential. Following this shift it is an easy
matter to add some small terms of degree four, five and six to the potential,
such that we get a six degree polynomial member with a shallow minimum. Though
we haven’t targeted potentials with minima especially, we have not ruled these
out a-priori, thus some potentials as the one in the example were manufactured
by the computerized model building package, and were tested for compliance.

The current work yields corrected predictions of this class of models by a
systematic high-precision analysis, thus providing a viable alternative to the
large field-high $r$ option. The analysis of BenDayan:2009kv is extended, in
preparation for a subsequent detailed comparison of the models to data. This
is done in order to simplify the parametrization of the potential and
facilitate a comprehensive numerical study.

####  2.0.1 Inflaton potentials with $r=0.001$

The following class of polynomial inflationary potentials proposed in
BenDayan:2009kv is:

$\displaystyle V(\phi)=V_{0}\left(1+\sum_{p=1}^{5}a_{p}\phi^{p}\right)\;.$
(2.1)

The virtue of these models from a phenomenological point-of-view is the
ability to separate the CMB region from the region of large e-fold production.
Hence, these potentials can produce a very different spectrum early on, than
in the later stages of inflation. Fig. 2.1 illustrates this point, with
separate CMB region and e-fold generation region. In the context of both
classification systems mentioned, current observational data weakly support
these Martin:2014lra; Vennin:2015eaa. However the small field model studied in
Martin:2014lra are monomial potential models of the form $V\propto
1-a_{p}\phi^{p}$, which are different from many of our models.

Fig. 2.1: A graph depicting $-1/\sqrt{2\epsilon}$ as a function of the
inflaton $\phi$ for a model for which $r_{0}=0.001$. The CMB interval is
covered by $\sim 8$ e-folds generated while the field changes by about
$\Delta\phi\sim 0.1$. Most of the e-folds are generated when $\phi$ reaches
$\sim 0.4$.

In many models, $\varepsilon_{V}\sim 1/N^{2}$, $\eta_{V}\sim 1/N^{2}$, and the
time derivative $\frac{d}{Hdt}$ can approximately be replaced with a factor of
$\frac{1}{N^{2}}$ Kosowsky:1995aa. In the above models, this standard
hierarchal dependence is broken, they have a more complicated dependence while
obeying the slow-roll conditions $\epsilon_{H}$, $\delta_{H}\ll 1$. In
BenDayan:2009kv it was shown that these models can be written as:

$\displaystyle
V(\phi)=V_{0}\left(1-\sqrt{\frac{r_{0}}{8}}\phi+\frac{\eta_{0}}{2}\phi^{2}+\frac{\alpha_{0}}{3\sqrt{2r_{0}}}\phi^{3}+a_{4}\phi^{4}+a_{5}\phi^{5}\right).$
(2.2)

Here $r_{0}$, $\eta_{0},\alpha_{0}$ are defined as
$r=8\left(\frac{V^{\prime}}{V}\right)^{2}$, $\eta=\frac{V^{\prime\prime}}{V}$,
$\alpha=-2\xi^{2}$, respectively. The subscript $0$ means that these are the
values at the CMB point.

Specifically for a potential of the form $V\propto
1+\sum_{p=1}^{5}a_{p}\phi^{p}$, the SL analytic expressions for the scalar
index and its running (Eqs. (1.179,1.180)) are given by

$\displaystyle n_{s}\simeq$ $\displaystyle 1-3a_{1}^{2}+4a_{2}$ (2.3)
$\displaystyle+2\Bigg{[}\frac{4a_{2}^{2}}{3}-\left(\frac{5}{3}-12b\right)\frac{a_{1}^{4}}{4}$
$\displaystyle-\left(8b+1\right)a_{1}^{2}a_{2}+\left(6b+2\right)a_{1}a_{3}\Bigg{]},$
$\displaystyle\alpha_{s}\simeq$ $\displaystyle
16a_{1}^{2}a_{2}-6a_{1}^{4}-2a_{1}a_{3}.$ (2.4)

####  2.0.2 Reduced parameter space

The potential in Eq. (2.2) is a small field candidate, which after some
scaling and normalization, depends on four free parameters. One parameter is
used for setting $r_{0}$ at the CMB point, and thus the predicted amplitude of
the GW signal produced, while the other two parameters are used to parametrize
the $n_{s}-\alpha_{s}$-plane. The fourth parameter determines the number of
e-folds from the CMB point to the end of inflation. $\phi_{\text{\tiny{end}}}$
is set to $1$ to simplify the analysis. It follows that

$\displaystyle\frac{1}{2}\left(\frac{V^{\prime}}{V}\right)^{2}{|_{\phi=1}}=1.$
(2.5)

Suppose we want inflation to end at $\phi=\zeta$, we can rescale $\phi$:

$\displaystyle\phi\rightarrow\tilde{\phi}=\frac{\phi}{\zeta}.$ (2.6)

In this formulation,

$\displaystyle
V=V_{0}\left(1+\sum_{p}a_{p}\zeta^{p}\tilde{\phi}^{p}\right)=V_{0}\left(1+\sum_{p}\tilde{a_{p}}\tilde{\phi}^{p}\right),$
(2.7)

where $\tilde{a_{p}}=a_{p}\zeta^{p}$. Since this is the same potential, it
follows the same CMB observables are produced. Thus, applying condition (2.5)
can be viewed as a scaling scheme for the different terms in the potential
which does not limit the generality of our results.

Substituting the expression for the potential and its derivative at $\phi=1$
we get:

$\displaystyle-\sqrt{2}=\frac{\sum_{p=1}^{5}p\cdot
a_{p}}{1+\sum_{p=1}^{5}a_{p}}.$ (2.8)

$a_{4}$ is now given in terms of the other coefficients:

$\displaystyle
a_{4}=\frac{-1}{4+\sqrt{2}}\left(\sqrt{2}+\sum_{p\in(1,2,3,5)}\left(p+\sqrt{2}\right)a_{p}\right)$
(2.9)

Using the standard definition for the number of e-folds
$N=\int_{t_{\text{\tiny{CMB}}}}^{t_{\text{\tiny{end}}}}Hdt\simeq-\int_{\phi_{\text{\tiny{CMB}}}}^{\phi_{\text{\tiny{end}}}}\frac{V}{V^{\prime}}d\phi$,
and the approximation $V(\phi)=1+\sum_{p=1}^{5}a_{p}\phi^{p}\simeq 1$ yields a
rough estimate for $a_{5}$ as a function of $N$,

$\displaystyle
N\simeq-\int_{0}^{1}\frac{V(a_{1},a_{2},a_{3},a_{5})}{V^{\prime}(a_{1},a_{2},a_{3},a_{5})}d\phi\simeq-\int_{0}^{1}\frac{d\phi}{V^{\prime}(a_{1},a_{2},a_{3},a_{5})}.$
(2.10)

This estimate is then used as a starting point to refine $a_{5}$ by solving
the background equations iteratively thereby obtaining the accurate
coefficient $a_{5}$ that yields the correct $N$. Thus a 4-dimensional
parameter space $r_{0}$, $a_{2}$, $a_{3}$, $N$ is defined. The parameters
$a_{2},a_{3}$ are constrained by the requirement $|a_{2}|,|a_{3}|\ll 1$,
$a_{1}$ is constrained by the observable value of $r$ and $a_{5}$ is
determined by the other parameters and by the number of e-folds (taken to be
in between $50\sim 60$). The PPS considered is in the range of the first
$\log(2500)\sim 8$ e-folds of inflation.

####  2.0.3 Inflaton potentials with $r=0.01$

In Wolfson:2016vyx, a class of small field inflationary models which can
reproduce the currently measured CMB observables, while also generating an
appreciable primordial GW signal was studied. The existence of such small
field models provides a viable alternative to the large field models that
generate a high Tensor-to-Scalar ratio. Our exact analysis was shown to give
accurate results Wolfson:2016vyx. Models which yield Tensor to Scalar ratio,
of less than $r\lesssim 0.003$ were previously studied in Wolfson:2016vyx. The
initial study additionally demonstrated a significant difference between
analytical Stewart-Lyth Stewart:1993bc; Lyth:1998xn estimates and the exact
results. This result should be confronted with analyses such as in
Martin:2013tda where the Stewart-Lyth expression is relied upon, and
Dodelson:2001sh in which the authors use a Green’s function approach and
perturbation theory, but assume the log of the input is well behaved. Our
method extends and improves the method of the model building technique
employed in BenDayan:2009kv; Hotchkiss:2011gz. Previous analytical work
Choudhury:2014kma; Chatterjee:2014hna; Choudhury:2015pqa has shown that a
fourth-order polynomial potential is sufficient to generate a high tensor-to-
scalar ratio, even up to $r\gtrsim 0.1$. However, it was hard to realize this
numerically. It was discovered in Wolfson:2016vyx, that a fifth-order
polynomial potential was required for generating $0.001\lesssim r\lesssim
0.003$. Furthermore a sixth-order polynomial seems to be required for a
tensor-to-scalar ratio greater than $r\gtrsim 0.003$. A simple explanation is
offered by observing (see Fig. 2.2) that increasing $r$ by factor $\sim 10$,
causes the e-folds per field excursion generated at the CMB window to decrease
by a factor of $\sim 3$. This means widening the CMB window and losing the
decoupling between the CMB window and the e-fold generating peak. Adding the
6th coefficient pushes the peak from $\phi\in[0.4,0.5]$, to higher values of
$\phi$ and decouples these regions.

Fig. 2.2: A graph depicting $-1/\sqrt{2\epsilon}\sim\frac{V}{V_{\phi}}$ as a
function of the inflaton $\phi$ for two fifth-order polynomial models, and a
sixth-order polynomial model. For a fifth-order polynomial model with
$r=0.001$ (Blue line) the CMB window width is $\sim 8$ e-folds, while the
field changes by about $\Delta\phi\sim 0.1$. Most of the e-folds are generated
when $\phi$ reaches $\sim 0.4$. When $r$ is increased the CMB window widens
and approaches the e-fold generating peak (Red dots). While marginally
affecting the CMB window width, the introduction of an additional coefficient,
$a_{6}$, allows shifting the peak to higher values of $\phi$, thereby
decoupling the CMB window and the e-fold generating peak (Green dash).

####  2.0.4 Inflaton potentials with $r=0.03$

We continue our investigations Wolfson:2016vyx; Wolfson:2018lel of a class of
inflationary models that were proposed by Ben-Dayan and Brustein
BenDayan:2009kv and were followed by Hotchkiss:2011gz; Antusch:2014cpa;
Garcia-Bellido:2014wfa. This class of models is compatible with several
fundamental physics considerations. Recently, interest in this class of models
was revived by the discussion about the “swampland conjecture",
Lehners:2018vgi; Garg:2018reu; Kehagias:2018uem; Ben-Dayan:2018mhe which
suggests that small field models are favoured by various string-theoretical
considerations (see Palti:2019pca for a recent review).

In addition, for this class of inflationary models, high values of $r$ result
in a scale dependence of the scalar power spectrum. Future experiments such as
Euclid Amendola:2012ys, and SPHEREx Dore:2014cca aim to probe the running of
the scalar spectral index $\alpha_{s}$ at the level of $1\sigma\simeq 1\times
10^{-3}$. This is a major improvement in comparison to the Planck bounds on
$\alpha_{s}$ which are currently at the level of $1\sigma\simeq 17\times
10^{-3}$. Such future measurements could provide additional constraints on our
models.

Fig. 2.3: Polynomial potentials. The blue (solid) line depicts a potential of
a model that predicts $r\simeq 0.03$, while the red (dash-dotted) line depicts
a potential that predicts $r\simeq 0.01$. The purple (dotted) line depicts a
degree five polynomial potential of a model that predicts $r\simeq 0.001$. All
models are variants of the hilltop model, with a flatter region in which most
e-folds are generated.

The small field models that we study are single-field models. The action of
such models is given by

$\displaystyle S=\int
d^{4}x\sqrt{-g}\left[\frac{R}{2}-\frac{1}{2}\partial^{\mu}\phi\partial_{\mu}\phi-V(\phi)\right].$
(2.11)

The metric is of the FRW form and the potential given by

$\displaystyle V(\phi)=V_{0}\left[1+\sum_{p=1}^{6}a_{p}\phi^{p}\right].$
(2.12)

Previously, in BenDayan:2009kv; Wolfson:2016vyx this class of models was
discussed from a phenomenological and theoretical points of view. In
Wolfson:2016vyx, the technical details of model building and simulation
methods were discussed, while in Wolfson:2018lel, the analysis and the
extraction of the most probable model were discussed. Additionally, in
Wolfson:2018lel, the most likely model which yields $r=0.01$ was identified.

Fig. 2.4: Field excursions in (reduced) Planck units. Different predicted
values of $r$ require different field excursions to generate the $\sim 8$
e-folds probed by the CMB. The model predicting $r\simeq 0.03$ (blue line)
requires an excursion of $\left(\Delta\phi\right)_{CMB}\simeq 0.28$ to
generate the same amount of e-folds which the model predicting $r\simeq 0.01$
model (red dash-dot) generates in $\left(\Delta\phi\right)_{CMB}\simeq 0.2$.
This means more tuning is required for models that predict $r\simeq 0.03$ .
The model predicting $r\simeq 0.001$ (purple dots) requires only
$\left(\Delta\phi\right)_{CMB}\simeq 0.1$ to generate the CMB window.

The small field models previously studied in Wolfson:2016vyx yielded results
that are consistent with observable data up to values of $r\simeq 0.003$.
While these values agree with the current limits on $r$ set by Planck
Ade:2015tva; Ade:2015xua, we are interested in studying models with higher
$r$. For models with $r\gtrsim 0.003$, significant running of running is
found. This means that while three free parameters (corresponding to
$n_{s},\alpha_{s},N$) were previously needed, we now need an additional free
parameter. Therefore we turn to a model of a degree six polynomial potential.
Obviously considering higher degree models complicates the analysis by adding
other tunable parameters. The potential is given by the following polynomial:

$\displaystyle V=V_{0}\left(1+\sum_{p=1}^{6}a_{p}\phi^{p}\right)\;.$ (2.13)

It has been shown BenDayan:2009kv that the potential can be written as:

$\displaystyle
V=V_{0}\left(1-\sqrt{\frac{r_{0}}{8}}\phi+\frac{\eta_{0}}{2}\phi^{2}+\frac{\alpha_{0}}{3\sqrt{2r_{0}}}\phi^{3}+a_{4}\phi^{4}+a_{5}\phi^{5}+a_{6}\phi^{6}\right).$
(2.14)

However, for simplicity, we express the potential as follows:

$\displaystyle
V=V_{0}\left(1-\sqrt{\frac{r_{0}}{8}}\phi+\sum_{p=2}^{6}a_{p}\phi^{p}\right)\;,$
(2.15)

with the subscript $0$ denoting the value at the CMB point. By setting
$\phi_{0}=0\;;\;\phi_{end}=1$ we limit ourselves to small field models in
which $\Delta\phi=1$ in Planck units, with little effect on CMB observables.
According to the Lyth bound Lyth:1996im; Easther:2006qu, given a Tensor-to-
Scalar ratio of $r\simeq 0.01$, the lower bound on the field excursion is
approximately given by $\Delta\phi_{4}\gtrsim 0.03\;m_{pl}$. Here
$\Delta\phi_{4}$ is the field excursion while the first $\sim 4$ efolds are
generated. Our models satisfy this strict bound, as the first 4 efolds or so
typically result in $\Delta\phi_{4}\sim 0.15$ which is well above $0.03$. The
Lyth bound was further extrapolated Efstathiou:2005tq to cover the entire
inflationary period. Applying this approach to models with $r\sim 0.01$ yields
$\Delta\phi\simeq 2\;m_{pl}$. However, in Hotchkiss:2011gz, it was shown that
in models such as the ones we study, the value of $\Delta\phi$ can be smaller
because $\epsilon_{H}$ is non-monotonic. In this case, $\Delta\phi=1\;m_{pl}$
from the CMB point to the end of inflation is consistent with the Lyth bound.

When the coefficients $\\{r_{0},a_{2},a_{3},a_{4}\\}$ are fixed, the remaining
coefficients are related by:

$\displaystyle a_{5}=f_{1}(r_{0},a_{2},a_{3},a_{4},a_{6}),$ (2.16)
$\displaystyle a_{6}=f_{2}(r_{0},a_{2},a_{3},a_{4},N).$ (2.17)

The procedure of finding $f_{1}$ and $f_{2}$ was explained in detail for the
degree 5 polynomial models in Wolfson:2016vyx, and here we follow a similar
procedure for the degree 6 models. So, ultimately, the model is parametrized
by 5 parameters: the two physical parameters $r_{0}$ and $N$ and the three
other parameters ($a_{2},a_{3},a_{4}$) that are used to parametrize the
$n_{s},\alpha_{s},\beta_{s}$ parameter space. It should be pointed out that
$N$ is not an observable, rather $N\sim 50-60$ is a ’soft’ constraint.
Strictly speaking, $N$ depends on the reheating temperature and only its
maximum value can be determined. However, for simplicity, we treat $N$ as an
observable, in order to facilitate the study of a large sample of models.

##  Chapter 3 Methods

####  3.0.1 Coefficient extraction methods

In this section, we explain the two methods for calculating the most likely
coefficients $\\{a_{2},a_{3},a_{4}\\}$, given a large number of simulated
models and the likelihood data for the CMB observables. This data is available
through CosmoMC Lewis:2002ah analysis of CMB data, such as the Planck data
Ade:2015tva.

Fig. 3.1: Small field inflationary potentials which yield $r=0.01$, as well as
PPS observables within 68% and 99% confidence levels. Every
$(n_{s},\alpha_{s})$ pair is accessible using these models. The likelihood
curves are results of a CosmoMC Lewis:2002ah run with the latest BICEP2+Planck
Ade:2015tva data.

##### Likelihood assignment method - Gaussian extraction

To each potential, after calculating the observables
$n_{s},\alpha_{s},\beta_{s}$, we assign a likelihood. For each observable, we
calculate the likelihood according to the MCMC likelihood analysis of the data
sets used. We then assign the product of the likelihoods $L_{(n_{s})}\times
L_{(\alpha_{s})}\times L_{(\beta_{s})}$ to the potential. A concrete example
is the following: suppose we extract the trio
$(n_{s},\alpha_{s},\beta_{s})=(0.96,0.011,0.024)$, we look up the likelihoods:
$(L_{\left(n_{s}=0.96\right)},L_{\left(\alpha_{s}=0.011\right)},L_{\left(\beta_{s}=0.024\right)})$.
We now multiply them, and so the likelihood attached to that specific model
which yielded these observables is given by
$L_{potential}=L_{\left(n_{s}=0.96\right)}\times
L_{\left(\alpha_{s}=0.011\right)}\times L_{\left(\beta_{s}=0.024\right)}\;.$
We proceed to extract likelihoods for the different coefficients by process of
marginalization. The expectation is that this method will yield a (roughly)
Gaussian distribution for each of the values of $a_{2},a_{3},a_{4}$. The
advantage of this method is in yielding not only the most likely value but
also the width of the Gaussian. This width can then be used as an indication
for the level of tuning that is needed in these models.

Fig. 3.2: The calculated likelihoods for the coefficients
$\\{a_{2}$,$a_{3}$,$a_{4}\\}$, in models with $r=0.01$. The most likely
coefficients are given by: $a_{2}=0.04$, $a_{3}=-0.15$, $a_{4}=0.02$. The
tuning level for each coefficient is given by Barbieri-Giudice measure
Barbieri:1987fn and is $(0.375,0.27,5.5)$.

##### Possible pitfalls

This method of likelihood assignment is vulnerable in two ways:

  * a)

To be valid, this method requires a uniform cover of the relevant parameter
space, by the potential parameters. If the cover significantly deviates from
uniform, the results might be skewed by overweighting areas of negligible
likelihood, or underweighting areas of significant likelihood. Fig. 3.1 shows
a mostly uniform cover.

  * b)

Since $L_{n_{s},\alpha_{s},\beta_{s}}\simeq L_{(n_{s})}\times
L_{(\alpha_{s})}\times L_{(\beta_{s})}$ only if the paired covariance is
small, we must make sure that this is the case. In our underlying MCMC
analysis this is indeed the case. The covariance terms are, in general, one to
two orders of magnitude smaller than the likelihoods at the tails of the
Gaussian.

  * c)

We also run the risk of false results if the fit we apply to the data points
produced by the numerical analysis yields a large fitting error. However, the
fitting error of the polynomial function to the $\log PPS-\log k$ data is
usually of the order of $10^{-6}$. This fitting is done over 30 data points
generated by the MS equation numerical evaluation, for each potential. The
error is calculated as $\Delta=\sqrt{\sum_{1}^{30}\left((\log
PPS)_{i}-\mathrm{fit}((\log k)_{i})\right)^{2}}$, thus the error per data
point is of the order of $10^{-6\sim 7}$. We conclude that the $\log PPS$
function is well fitted.

##### Multinomial fit

Another method for calculating the most likely coefficients is by fitting the
simulated data with a multinomial function of the CMB observables. We aim to
find a set of functions $F_{i}$ such that, for example,
$a_{2}=F_{2}(n_{s},\alpha_{s},\beta_{s})$. We assume that this function is
smooth and thus can be expanded in the vicinity of the most likely CMB
observables. Hence, we can find a set of multinomials $(F_{2},F_{3},F_{4})$,
such that:

$\displaystyle\begin{array}[]{ccc}a_{2}&=&F_{2}(n_{s},\alpha_{s},\beta_{s})\\\
a_{3}&=&F_{3}(n_{s},\alpha_{s},\beta_{s})\\\
a_{4}&=&F_{4}(n_{s},\alpha_{s},\beta_{s}).\end{array}$ (3.4)

We have found that a quadratic multinomial is sufficiently accurate and that
using a higher degree multinomial does not improve the accuracy significantly.
Thus we may represent these by a symmetric bilinear form plus a linear term,
as follows:

$\displaystyle F_{i}=OB_{i}O^{{\dagger}}+A_{i}O^{{\dagger}}+p_{0,i}\;,$ (3.5)

where $O=(n_{s},\alpha_{s},\beta_{s})$, $B_{i}$ is the bilinear matrix, and
the linear coefficient vector is $A_{i}$.

Fig. 3.3: Power spectra as recovered using CosmoMC Lewis:2002ah analysis with
latest BICEP2+Planck data Ade:2015tva. Allowed area (68% CL) for a fixed
$n_{s}$ analysis is shown (blue). Similarly, a fixed $\alpha_{s}$ analysis
(red), and a fixed $\beta_{s}$ (green) are shown. The other colors are
intersection areas. The pivot scale in this graph is at
$\log\left(\frac{k}{k_{0}}\right)=0$, where $k_{0}=0.05\;hMpc^{-1}$. The
apparent divergence in high k’s is due to the inability of Planck to constrain
these k’s. This is also shown in Fig. 3.4. With more data, it will be possible
to differentiate between the three possibilities.

##### Pivot scale

Fig. 3.4: Power spectra in the $C_{l}$’s decomposition (upper panel), with a
free $n_{s}$ (thick blue line), free $\alpha_{s}$ (thin cyan dots), and free
$\beta_{s}$ (medium red dash). The lower panel shows the relative difference
between the different cases. The relative difference (lower panel) is bound
from above by $\sim 1\%$. Additionally, the Planck observation error bars are
shown.

So far, we discussed matching potentials and their resulting PPS around the
CMB point. However, in order to correctly compare the results of the PPS to
observables, one has to take into account the pivot scale at which the CMB
observables are defined. Since, in this case, the pivot scale is given by
$k_{0}=0.05\;hMpc^{-1}$, and the CMB point is at $k\sim 10^{-4}\;hMpc^{-1}$ ,
the observables in the CMB point and $k_{0}$ should be related in a simple way
only if the spectrum varies slowly with $k$. This is not true for the case at
hand. Two potentials can yield very different power spectra near the CMB
point, and nevertheless yield the same observables at the pivot scale. These
degeneracies, stem from our limited knowledge of the power spectra on small
scales, and at the CMB point. For concreteness take two PPS functions, one
that is well approximated by a cubic fit near the pivot scale, and the other
that is well approximated only when we consider a quartic fit. Suppose,
additionally, that these two PPS functions have the exact same first three
coefficients, it follows that they yield the exact same observables
$\\{n_{s},\alpha_{s},\beta_{s}\\}$. However, if we go to sufficiently small
scales or high enough $k$ values, these functions will diverge. This is also
true at the large scale end, where the CMB point is set. Hence the degeneracy.
A possible solution to this problem is classifying the resulting power spectra
by the level of minimal good fit. We define a good fit as one in which the
cumulative relative error
$\Delta=\sqrt{\sum_{k}\left(\log(PPS(\log(k)))-fit(\log(k))\right)^{2}}$, is
less than $10^{-7}$. Given a single power spectrum, we fit our result with a
polynomial fit, increasing in order until the accumulated relative error is
sufficiently small. The minimal degree polynomial fit that approximates the
$\log(PPS)$ function to the aforementioned accuracy is called the minimal good
fit. We then study separately power spectra that are well fitted by cubic
polynomials, quartic polynomials etc. In this way we make sure that we compare
non-degenerate cases.

####  3.0.2 Monte Carlo analysis of Cosmic Microwave Background with running
of running

In Cabass:2016ldu it was shown that the inclusion of additional parameters,
i.e., the running of the spectral index ($\alpha_{s}$), and the running of the
running ($\beta_{s}$) resolves much of the tension between different data
sets. In this section, we briefly discuss the effect of considering non-
vanishing $\alpha_{s}$ and $\beta_{s}$ on the most likely shape of the PPS.
First, we find $n_{s}$ when it is the only free parameter. We then use
$n_{s}$, and $\alpha_{s}$ as the free parameters, and finally we conduct an
analysis with $n_{s},\alpha_{s}$ and $\beta_{s}$ as the free parameters. The
shape of the power spectrum changes significantly when running of running is
considered.
The data sets that were used are the latest BICEP2+Planck baseline
Ade:2015tva, along with the low $l$’s Bennett:2012zja, low TEB and lensing
likelihoods. The results of these analyses are given in Table 3.1, as well as
in Fig. 3.3. As expected the resulting power spectra converge at the pivot
scale $k_{0}=0.05\;hMpc^{-1}$. However, for lower $k$’s, the resulting spectra
diverge considerably, consistent with cosmic variance. Notably, the spectra
also diverge at higher $k$’s. This indicates the inability of current
observational data to constrain the models in this range of $k$’s. This
inability is also demonstrated in Fig. 3.4 where, for $l>1500$, the most
restrictive data cannot rule out models with significant running, or running
of running. Figure 3.4 also shows that the three models are virtually
indistinguishable in terms of the observed $C_{l}$’s.
The conclusion is that we will need additional accurate data from smaller
cosmic scales to be able to differentiate between the three scenarios. These
extra e-folds might come from future missions such as Euclid Amendola:2016saw,
or $\mu$-type distortion data Diacoumis:2017hhq; Abitbol:2017vwa.

Parameter (68%) | free $n_{s}$ | free $\alpha_{s}$ | free $\beta_{s}$
---|---|---|---
$\log(10^{10}A_{s})$ | $3.1047\pm 0.0057$ | $3.1073\pm 0.006$ | $3.1061\pm 0.0065$
$n_{s}$ | $0.9751\pm 0.0045$ | $0.973\pm 0.0057$ | $0.9687^{+0.0051}_{-0.006}$
$\alpha_{s}$ | N/A | $-0.009\pm 0.0067$ | $0.008\pm 0.013$
$\beta_{s}$ | N/A | N/A | $0.020\pm 0.013$
Table 3.1: Results from 3 analyses of the latest BICEP2+Planck dataset, each
adding a free parameter in the power spectrum. The results shown are best
fits, within the 68% confidence level for each analysis.

##  Chapter 4 The INSANE code

In order to assess the primordial power spectrum given an inflationary
potential, a stand-alone simulator was built. The code was given the name
INflationary potential Simulator and ANalysis Engine (INSANE) and is a fully
numerical simulator that solves the background and MS equations fully and
precisely, for a wide variety of symbolic potentials. Of course this is by no
means the first foray into numerical cosmology. Salopek, Bond and Bardeen
Salopek:1988qh have calculated power spectra resulting from different
potentials as early as 1989. Adams, Cresswell and Easther Adams:2001vc have
studied PPS responses to features in the potential, while Peiris et. al. have
utilized such codes to analyse the first-year results from WMAP mission
Peiris:2003ff. Mortonson, Dvorkin, Peiris and Hu Mortonson:2009qv have also
utilized such codes to study features of inflation. There are of course many
others who have solved such inflationary systems numerically, some have
published their codes (e.g. for example Price:2014xpa and Ringeval:2005yn;
Martin:2006rs; Ringeval:2007am). However the code presented here differs in
two main features:

  * (a)

Most if not all currently existing codes use the parameter flow equations,
first introduced by Kinney Kinney:2002qn, that replace time as the underlying
quantity over which integration is performed, with the number of efolds. Since
$N=\int Hdt$, it follows that $\frac{d}{dt}=H\frac{d}{dN}$ thus, the equations
of motion are now given by:

$\displaystyle\left\\{\begin{array}[]{rcl}\frac{d^{2}\phi}{dN^{2}}&=&\left(\frac{1}{2}\left(\frac{d\phi}{dN}\right)^{2}-3\right)\frac{d\phi}{N}-\frac{1}{H^{2}}\frac{dV}{d\phi}\\\
\frac{dH}{dN}&=&\frac{-H}{2}\left(\frac{d\phi}{dN}\right)^{2}.\end{array}\right.$
(4.3)

While this formulation does wonders to hasten the integration process, it is
our experience that this flow formulation, while analytically sound, results
in information loss in the computational implementation. For instance in a
tilted natural inflation scheme, near the limit of tilting the model such that
$\epsilon=1$ is accessible, there is a marked disagreement between the results
of integration with the two formulations. This can be seen in figure 4.1.

Fig. 4.1: While analytic integration over efolds, cosmic or conformal time are
equivalent, integration over conformal time is nigh impossible due to tendancy
toward an ’exploding’ solution. Numerical integration over efolds or cosmic
time usually give near identical results, but in some cases, as shown above
they diverge considerably.

To be perfectly clear, we do not state with absolute confidence that one is
better than the other, simply that they may differ. Since neither cosmic time
nor efold number are strictly observables, and we are unable to run an
experiment to validate one version over the other it is hard to prefer one
over the other. However we contend that this does merit some additional
discussion within the computational cosmology community. Further differences,
especially in potentials which sport a plateau, were discussed by Coone,
Roest, and Vennin Coone:2015fha.

  * (b)

At the time of writing the INSANE package, no computing language had both an
embedded symbolic math capability as well as good numerical solvers (the one
exception to this is Mathematica which was, at the time, opaque to ‘under-the-
hood’ scrutiny). Thus by writing this package using Python, we were able to
include fully symbolic parsers.

###  4.1 Code architecture

The code is organized into two main components each with its routines and sub-
routines. The chosen computing language was Python, since at the time it was
the only language to support both numerical integration and symbolic math,
along with an object-oriented-programming (OOP) paradigm. The OOP
compatibility was paramount for the construction of code that is reusable and
can accept any inflationary potential. Initially, the plan was to have the
code accept any physical action along with any underlying metric, construct
the appropriate Friedmann equations and then develop and solve the background
and MS equations to yield the full PPS. However, since we decided on studying
small field scalar potential models, this was deemed a wild overshoot and was
relegated to possible future projects. Additionally, the inclusion of a
Boltzman code, pyCAMB Lewis:1999bs, and a sky-map realization, HEALpy
Gorski:2004by, is a recent implementation and as such have yet to be fully
debugged.

####  4.1.1 Background evolution

This part of the code is given several parameters, for the evaluation of the
background geometry given an inflationary potential. The equations of
background evolution are given by the Friedmann equation, along with the
Klein-Gordon for the inflaton field:

$\displaystyle\left\\{\begin{array}[]{ccc}\dot{H}&=&-\frac{\dot{\phi}^{2}}{2}\\\
&&\\\ \ddot{\phi}&=&-3H\dot{\phi}-\frac{dV(\phi)}{d\phi}\end{array}\right.$
(4.7)

It is implicitly assumed that we are dealing with a FRW metric such that the
scale factor $a(t)$ can be gleaned directly from knowing
$H(t)\equiv\frac{\dot{a}}{a}$. The parameters supplied are the following:

Param. | Type | Default value | Description
---|---|---|---
V0 | double | $1$ |  $V_{0}$ the scalar potential at $\phi=0$
H0 | double | $\sqrt{\frac{V_{0}}{3}}$ |  $H_{0}$, Initial Hubble parameter
phi0 | double | $0$ | Initial $\phi$ value
phidot0 | double | $0$ | Initial $\dot{\phi}$ value
efolds | double | $\log{2500}$ | Number of efolds to simulate
|  |  | in the CMB window
efoldsAgo | double | 60 | Number of efolds in the simulated
|  |  | inflation, minimum is 50
Tprecision | double | 0.01 | Precision of time integration
|  |  | for the Background geometry solver
Table 4.1: Possible parameters to input into the code, for background
evolution simulation.

##### Potential

The inflationary potential that is supplied to the simulator can be in the
form of a single variable potential, for instance,
$V(\phi)=V_{0}\cos\left(c\phi\right)$ with some set of predefined $V_{0},c$,
or a set of numerical coefficients $\\{V_{0},a_{0},a_{1}...a_{n}\\}$ such that
the potential is a finite polynomial representation of some function:

$\displaystyle V(\phi)=V_{0}\sum_{p=0}^{n}a_{p}\phi^{p}.$ (4.8)

Currently, the code accepts polynomials up to degree 5, degree 6 and more
polynomials should be inserted as symbolic functions. The code then parses
this expression, and decides whether to manufacture a lookup table as a stand-
in for the potential at different values of $\phi$, or use the potential as
is, in the case of a polynomial representation.

##### Background Geometry

The Background geometry solver class receives the potential, parses it, and
integrates the Friedmann equations of the zero-order scalar field
$\phi^{(0)}(t)$. The design choice was to use the cosmic time equations, since
the conformal solution explodes to infinity faster than exponentially, close
to the end of inflation, which makes the integration step go down to zero.
Thus using an integration scheme that can handle exponential integration over
a large stretch of time was called for. The equations we integrate are given
by:

    H_rhs=(-1/float(2*Mpl**2))*phidot**2
    phi_rhs=phidot
    phidot_rhs=-3*H*phidot -vd  ,

where the rhs suffix means this is the coordinate derivative such that phi_rhs
means $\dot{\phi}$, and vd is the numerical equivalent of $dV(\phi)/d\phi$.
In order to solve the background equations, initial values must be supplied.
The BackgroundGeometry initiator method __init__ takes the following
arguments:

    Ψ
    (l1=0,l2=0,l3=0,l4=0,l5=0,Hubble0=None,Phi0=None,PhiEnd=1,
    Phidot0=None,Tprecision=0.01,endTime=200,Planck_Mass=1,
    physical=True,selfCheck=False,poly=True,
    Pot=(sym.Symbol(’x’))**0,mode=’silent’) ,

where l1 through l5 are a degree 5 polynomial coefficients such that
$V(\phi)=V_{0}\left(1+\sum_{p=0}^{5}l_{p}\phi^{p}\right)$.
Hubble0,Phi0,Phidot0 are initial values for $H_{0},\phi_{0},\dot{\phi}_{0}$
correspondingly. We can set the integration precision as well as several other
quantities. In general in the Python syntax, where an equality sign appears,
this is the default value, such that Tprecision=0.01 means the integration
step is $0.1$ ‘seconds’, unless otherwise specified.

####  4.1.2 Cosmic Perturbation

As outlined in 1.5.2, the equations for modes of scalar perturbation are given
by the set of MS equations (1.155), one for each $k$-mode. In order to solve
these, we first need to construct the pump field (1.33). The Background solver
contains, among others, the fields in Table 4.2, which are used to construct
the pump field $z$. One important note - the field a in the BackgrounGeometry
class denotes $\log{a(t)}$ instead of $a(t)$, and needs to be exponentiated
before construction of $z$.

Class field | type | Description
---|---|---
BG.a | ndarray(float128) | Instead of $a(t)$, this is $\log{a(t)}$
|  | or equivalently the number of efolds $N$
|  | from the start of inflation.
BG.phi | ndarray(float128) | $\phi(t)$
BG.phidot | ndarray(float128) | $\dot{\phi}(t)$
BG.H | ndarray(float128) | $H(t)$
Table 4.2: Some of the outputs of the BackgrondGeometry class, which are used
to construct the pump field $z=\frac{a\dot{\phi}}{H}$. The instance of
BackgrondGeometry in this table is BG.

After constructing the pump field it is possible now to solve the MS equation
for each mode.

##### The Mukhanov-Sasaki solver

The class responsible for solving the MS equations for each mode, and
construct the PPS is MsSolver. In general, the class is fed the quantities
$a,H,\phi,\dot{\phi},t$ at initiation as well as other parameters. MsSolver
then splines everything, finds the analytic expressions for the observables in
several approximations (Eqs. (1.179,1.180,1.184,1.196, etc.)), and finds the
CMB point by calculating a user-specified number of efolds back from the end
of inflation. It then finds the $k$ corresponding to the CMB point. This is
done by finding the pump field value at $t=t_{CMB}$, by setting the time-
dependent frequency in Eq. (1.157) $\omega_{k}^{2}(t)=0$. we thus find the $k$
number which is just now entering the horizon and call it $k_{0}$. This
corresponds to the scale of the universe and is, therefore, the physical
$k=1.2\cdot 10^{-4}\;\mathrm{h\;Mpc^{-1}}$. If we are interested in physical
inflation, this scale will be the default pivot scale, unless some other scale
was set as the pivot scale. A preparatory phase is immediately initiated after
the class is built, in which an interpolation table is created for the pump
field and all quantities needed for a solution of the MS equation to arbitrary
precision. The precision is predefined by the user, usually due to memory and
running time considerations vs. precision demands. The solver then proceeds to
integrate the MS equation for each $k$-mode, build a two variable function
$U(k,t)$ which is a stack of eigenfunctions $u_{k}(t)$ along the $k$ axis. The
solver finally builds the quantity

$\displaystyle
P_{s}(k,t)=\frac{k^{3}}{2\pi^{2}}\left|\frac{U(k,t)}{z(t)}\right|^{2},$ (4.9)

and evaluates this function at a time $t_{ps}$ much later than the freeze-out
time of the last mode of interest, to yield:

$\displaystyle
P_{s}(k)=\frac{k^{3}}{2\pi^{2}}\left|\frac{U(k,t_{ps})}{z(t_{ps})}\right|^{2}\;.$
(4.10)

This concludes the calculation phase, after which we have an array of
$k$-numbers, corresponding to the $k$-modes, and the PPS as a function of
these $k$-numbers. The analysis phase now begins, in which we assign the
pivot-scale specified. We first reassign $k_{0}$ to the correct physical
scale:

    Ψk=(k/k0)*1.2*10**(-4) .

If the user specified a pivot scale $k_{pivot}$, we further assign:

    Ψk=pivotScale*k/(1.2*10**(-4)) .

This effectively sets the point around which we will Taylor expand as the
pivot scale. The analysis phase continues by fitting the power spectrum with a
polynomial fit around the pivot scale. The polynomial fit degree starts at
$1$, i.e. linear fitting, and if the SSE value is over a set threshold, the
degree is increased. This is done until the SSE is lower than the set
threshold or the polynomial degree reaches 20. We assume that if polynomial of
degree 20 does not accurately fit the power spectrum, we either have some
error along the way, or we are describing a manifestly non-physical scenario.
In order to correctly initialize the MsSolver class the following input needs
to be supplied:

    Ψ
    (a,H,phi,phidot,t,epsilon,eta,xisq,Vd4V,Vd5V,V,
    eh,delh,pp,Tprecision=0.0064,Kprecision=15,
    log=True,pivotK=0.05,efoldsAgo=60,efoldsNum=8,
    check=False,physical=True,mode=’aprox’) ,

where a,H,phi,phidot,t,epsilon,eta,xisq,Vd4V,Vd5V,V,eh,delh,pp are all fields
manufactured in the BackgroundGeometry class, and held as internal fields
within it. Tprecision and Kprecision are precision parameters, for integration
time-step, and for the number of $k$-modes to solve and include in the
construction of the PPS. If log=True the code draws $k$-numbers such that the
distribution of $k$’s will be uniform in a log-log scale, to have equal
fidelity across the CMB window. pivotK,efoldsAgo,efoldsNum are the pivot
scale, number of efolds to consider in the inflation as ‘visible’ inflation,
and the number of efolds in the window to analyse as the CMB window,
correspondingly. We allow for non-physical inflationary scenarios in which
slow-roll is not held along the entire CMB window, in which case we will set
physical=False. Finally we consider the exact number of inflationary efolds by
setting mode=’exact’, otherwise, setting mode=’aprox’ means the CMB point will
be set at $\phi=0$. The difference between both cases is very slight, but may
be of import in some cases.

####  4.1.3 Usage

In order to correctly use INSANE, one must specify several arguments. These
arguments are conveniently edited in the params.in file. The full file along
with example inputs is given in Appendix B.

####  4.1.4 Benchmarking

Several tests are now in order to asses the voracity of this code. We first
test the code against the only known relevant analytical solution - the power
law inflation111There are other analytical solutions (c.f. Martin:2000ei) but
they are not relevant since they produce a very blue tilt, not in the vicinity
of the observable regime..

Fig. 4.2: Benchmarking our code against the analytical case of power-law
inflation reveals relative errors of order $10^{-3}\%$, which we take as a
sign of agreement between analytical and numerical cases.

##### The power law case

The power law inflation is the case where the Hubble parameter $H$ scales as
$H\propto t^{p}$ during inflation. The scalar potential that results in such a
case is given by:

$\displaystyle V(\phi)=A\exp{\left(-\sqrt{\frac{2}{p}}\phi\right)},$ (4.11)

in which $p$ is the same as the power law index. The resultant scalar index
$n_{s}$ of such an inflationary scenario is given by:

$\displaystyle n_{s}=1-\frac{2}{p},$ (4.12)

and is scale independent, such that the index running $\alpha_{s}$ vanishes by
definition. We have run the INSANE code against these scenarios to compare
analytic and numeric solutions. The results are shown in Fig. 4.2, and reveal
an agreement between analytical and numerical analysis to a relative error of
order $10^{-3}\%$. In evaluating the fidelity of $\alpha_{s}$ there is an
intrinsic problem since in the analytic case it vanishes by definition. Thus
no intrinsic scale is available for us to relate to. We, therefore, took the
heuristic approach of assessing $\Delta n_{s}/\Delta\log{k}$. If $\alpha_{s}$,
as recovered by the code, is lower than $\Delta n_{s}/\Delta\log{k}$, we are
satisfied that the code yields a correct $\alpha_{s}$.

Fig. 4.3: A heuristic study of $n_{run}$ (or $\alpha_{s}$), reveals the values
recovered numerically are consistently below the error threshold for $n_{s}$.
We take this as a sign of $n_{run}$ being at the level of numerical error,
hence it is equivalent to a vanishing running.

We find that the $\alpha_{s}$ values recovered by the numerical calculation
are always below the numerical error in $n_{s}$ divided by the overall
$\log{k}$. We take this as a sign that the recovered $\alpha_{s}$ is at the
level of numerical error, thus equivalent to a vanishing $\alpha_{s}$. Fig.
4.3 shows these findings. Additionally, we show the current observational
threshold of $\alpha_{s}$.

##### The quadratic case

Fig. 4.4: Relative error (in percent) between numerical results and the SL
analytical expression (black X’s). The errors converge to $0$ for large values
of $N$. Shifting the number of the efolds by $N\rightarrow N-0.8$ yields a
relative error of the order of $10^{-3}\sim 10^{-4}\%$ (green pluses).

As we aim to study models that produce slow-roll parameters which are time-
dependent, we need to check the precision of the numerical code against such
models. Consequently, we tested the accuracy of our calculations for quadratic
potentials of the type

$\displaystyle V=\frac{1}{2}m^{2}\phi^{2}.$ (4.13)

In these cases, the analytic expression for the scalar index is given by,

$\displaystyle n_{s}=1-\frac{8}{4N+2}+\frac{32b}{\left(4N+2\right)^{2}},$
(4.14)

Here $N$ is the number of efolds and $b$ is the same as in (1.179). Fig. 4.4
presents the results of this study, as relative errors between precise
calculations and the SL analytic expressions. These results are accurate to
$\sim 0.1\%$. However, there is a systematic error that is traced back to the
inaccuracy of the approximation:

$\displaystyle
N=\int_{t_{\text{\tiny{CMB}}}}^{t_{\text{\tiny{end}}}}Hdt\simeq-\int_{\phi_{\text{\tiny{CMB}}}}^{\phi_{\text{\tiny{end}}}}\frac{V}{V^{\prime}}d\phi.$
(4.15)

A shift $N\rightarrow N-0.8$ is sufficient to reduce the systematic error such
that the relative error is of the order of $10^{-3}\sim 10^{-4}\%$. Additional
types of simple potentials, which yield time-dependent slow-roll parameters
were also studied. In all cases, the relative error between calculated results
and the traditional SL expression Eq. (5.20) is bounded from above by $\sim
0.1\%$. Furthermore, a more careful analytical treatment leads to better
accuracy, bounded from above by about 0.02% relative error. Additionally, we
were able to recover the “Cosmic ring” phenomenon, that is the PPS response to
a step function in the potential. This response feature in the PPS was first
studied in Adams:2001vc.
We take all these results as a strong indication of sufficient accuracy of our
calculations.

###  4.2 Improvements - Teaching a PC some physics

In order to speed up the calculation, and enable the study of a class of
inflationary potentials rather than a solitary example, some adaptations of
the code were called for. The most crucial and time-consuming element of
calculating the power spectrum is the numerical integration of the different
$k$-modes that make up the PPS. We, therefore, identified several mechanisms
that, while technically might be interesting, produce little to no effect in
the context of yielding the PPS.

####  4.2.1 Finding the mode ‘freeze-out’ point

The first and most influential improvement was achieved by using the physical
understanding that each mode of quantum oscillation ‘freezes out’ of the
horizon at some point in its evolution. That is, at some point of the
inflationary evolution the light horizon becomes smaller than the mode’s
wavelength. If we take the viewpoint of the classical harmonic oscillator such
that

$\displaystyle U_{k}=A_{k}e^{i\omega t},$ (4.16)

with $\omega^{2}=k^{2}-\frac{z^{\prime\prime}}{z}$, which at some point
becomes negative, it is immediately apparent that $\omega$ becomes purely
imaginary around

$\displaystyle k^{2}=\frac{z^{\prime\prime}}{z}.$ (4.17)

And so, since we care only about the modes after the freeze-out time, we can
in theory start the integration approximately around the freeze-out time of