# Non-Gaussian Signatures of a Thermal Big Bang

Maria Mylova  , Marianthi Moschou  , Niayesh Afshordi  and João Magueijo

(today)

###### Abstract

What if Big Bang was hot from its very inception? This is possible in a
bimetric theory where the source of fluctuations is thermal, requiring the
model to live on a critical boundary in the space of parameters and can be
realized when an anti-DBI brane moves within an $EAdS_{2}\times E_{3}$
geometry. This setup renders the model unique, with sharp predictions for the
scalar spectral index and its running. We investigate the non-Gaussian
signatures of this thermal bimetric model, or “bi-thermal” for short. We adapt
the standard calculation of non-Gaussianities for $P(X,\phi)$ models to the
thermal nature of the model, emphasising how the bi-thermal peculiarities
affect the calculation and alter results. This leads to precise predictions
for the shape and amplitude of the three-point function of the bi-thermal
model (at tree-level): $f^{\rm local}_{\rm NL}=-3/2$ and $f^{\rm equil}_{\rm
NL}=-2+4\sqrt{3}\pi/9\simeq 0.4$. We also discover a new shape of flattened
non-gaussianity $\propto(k_{1}+k_{2}-k_{3})^{-3/2}+$ permutations, which is
expected due to the excited thermal initial conditions. These results, along
with our earlier predictions for the scalar power spectrum, provide sharp
targets for the future generation of cosmological surveys.

##  1 Introduction

The inflationary paradigm as proposed in [1] correctly accounts for the Cosmic
Microwave Background (CMB) anisotropies as observed by Planck [2]. The initial
condition is provided by quantum fluctuations which freeze during a rapid
phase of the exponential expansion of the universe. This provides the seeds
responsible for the structure formation of the universe as we know it today.

An alternative approach is to consider that the initial conditions were
thermal in nature, as it was originally proposed in [3]. In this case the
universe starts in a thermal bath of finite temperature $T$ (see also [4, 5]).
Assuming the thermal fluctuations experience a varying speed of sound [6, 7,
8], where the acoustic horizon ($kc_{s}<Ha$) rapidly shrinks, modes that
originally were oscillating find themselves frozen, i.e., outside of the
Hubble horizon ($k<aH$).

Figure 1: Current (forecasted) constraints on spectral index and running from
Planck 2018++ [2] (CMB S4 [9]), compared to predictions from bi-thermal model
(Equation 1.1) and $R^{2}$ (Starobinsky) inflation [10], given the allowed
reheating range.

In this work, we focus on the thermal anti-DBI model, “bi-thermal”for short,
previously discussed in [11]. The bi-thermal model is very interesting in that
all its free parameters and predictions are fundamentally fixed [11]:

$n_{s}=0.96478(64),~{}~{}\frac{dn_{s}}{d\ln k}=-1.8\times 10^{-3},$ (1.1)

by the amplitude of the scalar fluctuations, $A_{s}\simeq 2.1\times 10^{-9}$
(at 0.05 /Mpc) [12] (see Fig. 1). This is partly because the fluctuations in
this scenario have a thermal origin, and thus the model dispensing with a
reheating phase, in contrast to inflationary models. In the latter, the
reheating temperature (or number of e-foldings) introduces a constrained
uncertainty in the predictions, for example, of the scalar spectral index
$n_{s}$ and gravitational waves [13].

Obviously, inflationary models could in principle have all the properties of
reheating fixed by knowledge of the particle physics parameters controlling
the reheating process. However, these parameters are typically extraneous to
the field producing inflation, so that this claim can only be substantiated in
the rare cases where inflation is embedded in a given particle physics model,
such as the Higgs inflation model [14]. None of these problems affects the
predictivity of the bi-thermal model.

By mapping it into an “Einstein frame”, the bi-thermal model can be understood
as an anti-DBI model [6, 15, 16]. But as much as the model differs from
inflation, it is also to be contrasted with these alternative theories in that
the model is “critical”, in the sense that it lives on a discontinuous
boundary. The models in [15, 16] (based on a vacuum quantum fluctuations) do
not have this property, so that deviations from exact scale-invariance can be
continuously dialled [16]. In the critical thermal model not only is exact
scale invariance non-achievable, but the discontinuity imposes a very specific
pattern of deviations from scale-invariance, with the amplitude of the
fluctuations fixing where, on the steep slope of the discontinuity, the
fluctuations we do observe must lie. This fixes the spectral index $n_{s}$ of
the observed fluctuations, now deemed to be not only “close to 1”, or
“generically red”, but also with an $n_{s}$ and its scale-dependence, fixed to
several figures (see equation (1.1)). The combination of the two predictions
renders the model a perfect target for upcoming experiments and observations
[9].

One important matter that was not investigated in [11] is the non-Gaussian
signature of these models. This aspect is particularly important for the non-
critical, vacuum based anti-DBI models, because it is one of the few aspects
in terms of which they are predictive [16]. Working out the 3-point function
for the critical model of [11] is the purpose of this paper.

To start with, the matter looks trivial. Famous last words… The problem is
remarkably complex and this is the reason for the ungainly time delay between
the publication of [11] and the present paper. The trivial aspects of the
problem are presented in Sections (2) and (3), where we review the two-point
function set up and its extension to the vertices of the theory. The 3-point
function is evaluated in Section (4). Unfortunately, the down to earth aspects
of the calculation that follows the laying out of the formalism are far from
trivial, as we find in Sections (5) and (6). This does not happen for
particularly deep reasons, but (at the technical level) boils down to the fact
that we need to evaluate integrals of triple products of Bessel functions away
from the orders usually found in other models. Hence, a number of assumptions
behind the approximations made in the literature are simply not valid here.

Nonetheless, using a combination of analytical and numerical methods we are
able to bring the calculation to good port. The remarkable conclusions are as
follows: the non-Gaussianity of the bi-thermal model is described by three
non-Gaussianity shapes; equilateral, squeezed and a new flattened shape.
Therefore, the bi-thermal model can leave several non-Gaussian fingerprints in
the early universe, with the bispectrum in the equilateral limit having an
amplitude $f_{NL}^{\text{equil}}=-\frac{2}{9}\quantity(9+2\sqrt{3}\pi)$, while
in the squeezed limit $f_{NL}^{\rm local}=-\frac{3}{2}$, both well in
agreement with current observations [17]. The flattened non-Gaussianity is
found to be $\propto(k_{1}+k_{2}-k_{3})^{-3/2}$. Finally, we conclude in
Section (7). A technical Appendix can be found in Section (A).

##  2 The bi-thermal model

###  2.1 The model

The bi-thermal model, as presented in [11], can be nicely written as a scalar-
tensor theory that takes on the following ‘anti-DBI’ form

$S=\int
d^{4}x\sqrt{-g}\left[\frac{1}{B(\phi)}\sqrt{1+2B(\phi)X}-V(\phi)\right],$
(2.1)

where $\phi$ is a scalar field, $X\equiv-\partial_{\mu}\phi\partial^{\mu}\phi$
is its kinetic term (here we are using a mostly + metric convention) and $g$
is the metric determinant. The ‘anti-DBI’ name tag derives from the fact that
this is essentially a DBI action with an opposite sign warp factor $B(\phi)$.
Using Bianchi identities to fix the form of the potential and fixing the
residual ambiguity by choosing a so-called “critical solution” (for details
see [11]) we arrive at a particular form for the warp factor and potential

$\displaystyle B(\phi)$ $\displaystyle=B_{0}\left(\frac{\phi}{M_{\rm
Pl}}\right)^{2}$ $\displaystyle V(\phi)$
$\displaystyle=\frac{3}{4B_{0}}\ln\left(\frac{\phi}{M_{\rm Pl}}\right)^{2}.$
(2.2)

With this as model-specific input, one can then use standard results for
$P(X,\phi)$ theories to investigate the correlation functions this model gives
rise to in the early universe.

When considering the non-Gaussian signatures associated with the bi-thermal
model, three features of this model and its ‘critical’ solution (2.2) are
worth highlighting at this point. Firstly, the model can be understood as
coming from the DBI action of a 3-brane embedded in a $EAdS_{2}\times E_{3}$
geometry. As such the action will be protected by a non-linear realisation of
the associated symmetries and the particular DBI-like form will in fact lead
to one of the vertices in the cubic action vanishing – we will discuss this in
Section (3). Secondly, no scaling solutions, so no exactly scale invariant
solutions can be found, which is a special feature of the potential
(2.2).111Other power-law warp factors lead to power-law potentials with
scaling solutions. This is important, because analytic expressions for the
non-Gaussian amplitude of a $P(X,\phi)$ theory are easiest to derive for exact
scale invariance [18], in which case slow-roll need not be assumed, or up to
leading order in slow-roll parameters [19], in which case scale-invariance
need not be assumed. Finally, the speed of sound of such models will tend to
infinity/diverge in the UV/at early times. As such it will also be varying
rapidly and a slow-roll approximation becomes inappropriate (for closely
related non-Gaussian computations in more standard $P(X,\phi)$ models, see
[16, 20]). We will therefore be interested in a particular setup, where we
have an almost, but not exactly scale-invariant solution and an extremely
large speed of sound at the time CMB correlations are sourced.

###  2.2 The quadratic action and mode functions

With an eye on computing non-Gaussian statistics later on (i.e. three-point
statistics descending from the cubic effective action in perturbations), it is
useful to quickly consider the quadratic regime first and identify the mode
functions that will also enter higher-point statistics. We define the
following two variables

$\displaystyle q$ $\displaystyle\equiv\frac{a\sqrt{2\epsilon}}{\sqrt{c_{s}}},$
$\displaystyle y$ $\displaystyle\equiv\int\frac{c_{s}dt}{a}=\int c_{s}d\eta,$
(2.3)

where $q$ is expressed in terms of the scale factor $a$, the speed of sound of
scalar perturbations $c_{s}$ and the standard slow-roll parameter
$\epsilon\equiv-\dot{H}/H^{2}$, while the ‘sound horizon time’ variable $y$ is
related to physical time $t$ and conformal time $\eta$. In terms of these we
can define the so-called Mukhanov variable (a canonically-normalized scalar
variable) $v=M_{\rm Pl}q\zeta$, where $\zeta$ is the usual coming curvature
perturbation. This way we can write the Mukhanov-Sasaki equation

$v^{\prime\prime}_{k}+\left(k^{2}-\frac{q^{\prime\prime}}{q}\right)v_{k}=0,$
(2.4)

which governs the evolution of $\zeta$ at the level of the quadratic action
and where the $q$ parameter, we already encountered above, satisfies

$\frac{q^{\prime\prime}}{q}=\frac{\nu^{2}-1/4}{y^{2}}.$ (2.5)

As such, the mode equation becomes

$v^{\prime\prime}_{k}+\left[k^{2}-\frac{\nu^{2}-1/4}{y^{2}}\right]v_{k}=0.$
(2.6)

In the usual Bunch-Davies solution, one would set $\nu=3/2$ in order to
recover scale invariance. However, for the thermal case we are considering
here, we have $\nu=1$ for a scale-invariant solution. This means we now have a
modified differential equation for $q$, with solutions

$q\sim(-y)^{-1/2}\qquad{\rm or}\qquad q\sim y^{3/2}\,,$ (2.7)

where $y$ runs from $-\infty$ to $0$, so the growing mode solution is now
$q\sim(-y)^{-1/2}$. Correspondingly, we now have the following solution for
the Mukhanov variables

$v_{k}(y)=\frac{\sqrt{-\pi y}}{2}e^{-i\gamma_{\nu}}H^{(2)}_{\nu}(ky),$ (2.8)

where $\gamma_{\nu}=\frac{\pi}{4}(2\nu+1)$ and $\nu\sim 1$ for near scale-
invariant solutions, as discussed above.

The transition from a “v” variable to a “u” variable (directly linked to
$\zeta$) is now given by [18, 16]:

$\displaystyle u_{k}(y)$ $\displaystyle=$
$\displaystyle\frac{1}{a}\left(\frac{c_{s}}{2\epsilon}\right)^{1/2}v_{k}(y)$
(2.9) $\displaystyle=$
$\displaystyle\frac{1}{a}\left(\frac{c_{s}}{2\epsilon}\right)^{1/2}\frac{\sqrt{-\pi
y}}{2}e^{-i\gamma_{\nu}}H^{(2)}_{\nu}(ky),$

where $\gamma_{\nu}=\frac{\pi}{4}(2\nu+1)$ and where the argument $\nu$ of the
Hankel function will eventually be set to (approximately) one, implying a
(near) scale-invariant two-point function. At a constant comoving temperature
$T_{c}=aT/c_{s}$, the power spectrum of scalar perturbations is given by:

${\cal P}_{\zeta}(k)=\lim_{y\rightarrow
0^{-}}\frac{k^{3}}{2{\pi}^{2}}\frac{{\left|v_{k}\right|}^{2}}{q^{2}M^{2}_{P}}[2\langle
n_{k}\rangle_{T_{c}}+1]\simeq\frac{c_{s}T_{c}(-y)^{1-2\nu}}{2\pi^{3}a^{2}\epsilon
M^{2}_{P}}\bigg{|}_{y\sim-k^{-1}},$ (2.10)

where we used the Bose-Einstein distribution for phonon occupation number
$\langle n_{k}\rangle_{T_{c}}\simeq T_{c}/k$ for $k\ll T_{c}$.

##  3 Non-Gaussian vertices

Having computed the mode functions above, we are now in a position to start
working with the cubic effective action and to compute non-Gaussian signals.

###  3.1 Preliminaries

We would like to evaluate the three-point function in the interaction picture,
following the procedure outlined by [21, 19, 22]. In the interaction picture
we are, therefore, after

$\langle\zeta(t,\textbf{k}_{1})\zeta(t,\textbf{k}_{2})\zeta(t,\textbf{k}_{3})\rangle=-i\int_{t_{0}}^{t}{\rm
d}t^{\prime}\langle[\zeta(t,\textbf{k}_{1})\zeta(t,\textbf{k}_{2})\zeta(t,\textbf{k}_{3}),H_{\rm
int}(t^{\prime})]\rangle~{},$ (3.1)

where we will specify the range of integration later. This result is true to
first order in $H_{\rm int}$ (higher order corrections can be computed using
‘nested’ commutators with $H_{\rm int}$) and for the cubic interactions we
will be interested in $H_{\rm int}=-{\cal L}_{\rm int}$, after all terms
proportional to the equations of motion have been removed via field
redefinitions.

In order to explicitly compute (3.1), we need to expand the curvature
perturbation $\zeta$ in terms of raising and lowering operators, as usual

$\zeta(y,{\bf k})=u_{k}(y)a({\bf k})+u_{k}^{*}(y)a^{\dagger}(-{\bf k}).$ (3.2)

Note that the mode functions $u_{k}$ in (3.2) satisfy the corresponding
Mukhanov-Sasaki equation and we will be applying the commutation relations
$[a({\bf k}),a^{\dagger}({\bf k}^{\prime})]=(2\pi)^{3}\delta^{3}({\bf k}-{\bf
k}^{\prime})$, where the factor of $(2\pi)^{3}$ implicitly fixes our Fourier
convention.

###  3.2 Leading order vertices

The cubic effective action, from which we will calculate the non-Gaussian
amplitudes, is the usual [21, 19, 22]

$\displaystyle S_{(3)}$ $\displaystyle=$ $\displaystyle M_{\rm Pl}^{2}\int{\rm
d}t{\rm
d}^{3}x\left\\{-a^{3}\left[\Sigma\left(1-\frac{1}{c_{s}^{2}}\right)+2\lambda\right]\frac{\dot{\zeta}^{3}}{H^{3}}+\frac{a^{3}\epsilon}{c_{s}^{4}}(\epsilon-3+3c_{s}^{2})\zeta\dot{\zeta}^{2}\right.$
$\displaystyle+$
$\displaystyle\frac{a\epsilon}{c_{s}^{2}}(\epsilon-2\epsilon_{s}+1-c_{s}^{2})\zeta(\partial\zeta)^{2}-2a\frac{\epsilon}{c_{s}^{2}}\dot{\zeta}(\partial\zeta)(\partial\chi)$
$\displaystyle+$
$\displaystyle\left.\frac{\epsilon}{2a}(\partial\zeta)(\partial\chi)\partial^{2}\chi+\frac{\epsilon}{4a}(\partial^{2}\zeta)(\partial\chi)^{2}+2f(\zeta)\left.\frac{\delta
L_{(2)}}{\delta\zeta}\right|_{1}\right\\}~{}.$

Now, the first vertex contribution $(\dot{\zeta}^{3})$ vanishes for anti-DBI
(and indeed DBI) models and the final term can be removed via a field
redefinition. Here, dots denote derivatives with respect to proper time $t$,
$\partial$ is a spatial derivative, we have explicitly imposed that
$\eta=\dot{\epsilon}/H\epsilon=0$ and $\chi$ is defined via

$\partial^{2}\chi=\frac{a^{2}\epsilon}{c_{s}^{2}}\dot{\zeta}\,.$ (3.4)

We now perform the field re-definition removing the final term (this only
introduces additional terms proportional to $\eta$, i.e. terms which we will
drop), so that for our model we are left with

$\displaystyle S_{(3)}$ $\displaystyle=$ $\displaystyle M_{\rm Pl}^{2}\int{\rm
d}t{\rm
d}^{3}x\left\\{\frac{a^{3}\epsilon}{c_{s}^{4}}(\epsilon-3+3c_{s}^{2})\zeta\dot{\zeta}^{2}\right.+\frac{a\epsilon}{c_{s}^{2}}(\epsilon-2\epsilon_{s}+1-c_{s}^{2})\zeta(\partial\zeta)^{2}-2a\frac{\epsilon}{c_{s}^{2}}\dot{\zeta}(\partial\zeta)(\partial\chi)$
(3.5) $\displaystyle+$
$\displaystyle\left.\frac{\epsilon}{2a}(\partial\zeta)(\partial\chi)\partial^{2}\chi+\frac{\epsilon}{4a}(\partial^{2}\zeta)(\partial\chi)^{2}\right\\}~{},$

We can now re-write this, explicitly, in terms of $y$-time and find the
following cubic effective action

$\displaystyle S_{(3)}$ $\displaystyle=$ $\displaystyle M_{\rm Pl}^{2}\int{\rm
d}y{\rm
d}^{3}x\frac{a}{c_{s}}\left\\{\frac{a\epsilon}{c_{s}^{2}}(\epsilon-3+3c_{s}^{2})\zeta{\zeta}^{{}^{\prime}}{}^{2}\right.+\frac{a\epsilon}{c_{s}^{2}}(\epsilon-2\epsilon_{s}+1-c_{s}^{2})\zeta(\partial\zeta)^{2}$
(3.6) $\displaystyle-$ $\displaystyle
2\frac{\epsilon}{c_{s}}{\zeta}^{{}^{\prime}}(\partial\zeta)(\partial\hat{\chi})+\left.\frac{\epsilon}{2a}(\partial\zeta)(\partial\hat{\chi})\partial^{2}\hat{\chi}+\frac{\epsilon}{4a}(\partial^{2}\zeta)(\partial\hat{\chi})^{2}\right\\}~{},$

where ${\zeta}^{{}^{\prime}}\equiv(d/dy)\zeta$ and we now have

$\partial^{2}\hat{\chi}=\frac{a\epsilon}{c_{s}}{\zeta}^{{}^{\prime}}\,.$ (3.7)

Ignoring the leading $1/c_{s}$ factor in (3.6), all terms in the second line
go as $1/c_{s}^{2}$, whereas the top line has contributions independent of
$c_{s}^{2}$. In the large $c_{s}$ limit, which is what we are interested in
here, these terms will dominate. So we can write

$\displaystyle S_{(3)}^{\text{leading order}}$ $\displaystyle=$ $\displaystyle
M_{\rm Pl}^{2}\int{\rm d}y{\rm
d}^{3}x\frac{a^{2}\epsilon}{c_{s}}\left\\{3\zeta{\zeta}^{{}^{\prime}}{}^{2}-\zeta(\partial\zeta)^{2}\right\\}.$
(3.8)

Note that this means the non-Gaussian amplitude is ‘slow-roll suppressed’ in
this model, although of course there is no requirement that $\epsilon\ll 1$
here, so this is not necessarily a suppression.

##  4 Direct Evaluation of 3-point function

In this section, we will directly compute the 3-point function, starting from
the third order action. Let us combine the quadratic action for $\zeta$ with
the leading correction (3.8):

$S=\frac{M_{\rm Pl}^{2}}{2}\int dyd^{3}xq^{2}(y)\left[(1+3\zeta)\zeta^{\prime
2}-(1+\zeta)(\partial\zeta)^{2}+{\cal O}(\zeta^{4})\right].$ (4.1)

The Euler-Lagrange equations are then given by:

$q^{-2}[q^{2}(1+3\zeta)\zeta^{\prime}]^{\prime}=\nabla\cdot[(1+\zeta)\nabla\zeta]+3\zeta^{\prime
2}-(\partial\zeta)^{2}+{\cal O}(\zeta^{3}),$ (4.2)

which simplifies to

$q^{-2}(q^{2}\zeta^{\prime})^{\prime}-\partial^{2}\zeta=-3\zeta\left[q^{-2}(q^{2}\zeta^{\prime})^{\prime}-\partial^{2}\zeta\right]-2\zeta\partial^{2}\zeta+{\cal
O}(\zeta^{3}).$ (4.3)

Using the first order equation of motion, the first term on the right hand
side vanishes at second order, which yields:

$q^{-2}(q^{2}\zeta^{\prime})^{\prime}-\partial^{2}\zeta=-2\zeta\partial^{2}\zeta+{\cal
O}(\zeta^{3}),$ (4.4)

and equivalently in Fourier space:

$q^{-2}(q^{2}\zeta_{\bf k}^{\prime})^{\prime}+k^{2}\zeta_{\bf
k}=2\int\frac{d^{3}k^{\prime}}{(2\pi)^{3}}|{\bf k-k^{\prime}}|^{2}\zeta_{\bf
k^{\prime}}\zeta_{\bf k-k^{\prime}}+{\cal O}(\zeta^{3}).$ (4.5)

The next step is to solve this equation perturbatively. The retarded Green’s
function for the linear $\zeta$ equation can be written in terms of the mode
function $u_{k}(y)$

$\begin{split}G_{k}(y,w)&=\frac{q^{2}(w)\left[u_{k}(y)u_{k}^{*}(w)-u_{k}^{*}(y)u_{k}(w)\right]}{q^{2}(y)\left[u^{\prime}_{k}(y)u_{k}^{*}(y)-u^{\prime*}_{k}(y)u_{k}(y)\right]}\Theta(y-w)\\\
&=iq^{2}(w)\left[u_{k}(y)u_{k}^{*}(w)-u_{k}^{*}(y)u_{k}(w)\right]\Theta(y-w),\end{split}$
(4.6)

which can be used to solve the second order equation of motion, at leading
order:

$\begin{split}\zeta^{(2)}_{\bf k}(y)&=\zeta^{(1)}_{\bf
k}(y)+2i\int_{-\infty}^{y}dwq^{2}(w)\int\frac{d^{3}k^{\prime}}{(2\pi)^{3}}|{\bf
k-k^{\prime}}|^{2}\\\
&\times\left[u_{k}(y)u_{k}^{*}(w)-u_{k}^{*}(y)u_{k}(w)\right]\zeta^{(1)}_{\bf
k^{\prime}}(w)\zeta^{(1)}_{\bf k-k^{\prime}}(w).\end{split}$ (4.7)

We can now compute the 3-point function:

$\begin{split}\lim_{y\rightarrow 0^{-}}\langle\zeta_{\bf k_{1}}\zeta_{\bf
k_{2}}\zeta_{\bf k_{3}}\rangle=(2\pi)^{3}\delta^{3}({\bf
k_{1}+k_{2}+k_{3}})(5/3)^{3}B_{\Phi}(k_{1},k_{2},k_{3})\end{split}$ (4.8)

where the bispectrum is given by

$\begin{split}B_{\Phi}(k_{1},k_{2},k_{3})&=-\frac{432\pi^{4}P^{2}_{\zeta}}{125(k_{1}k_{2}k_{3})^{2}}\Bigg{\\{}\int_{-\infty}^{0}dw\,q^{2}(w)k^{2}_{1}(k_{2}^{2}+k_{3}^{2})\\\
&\times\frac{\imaginary\left[u_{k_{1}}(0)u_{k_{1}}^{*}(w)\right]\real\left[u_{k_{2}}(0)u_{k_{2}}^{*}(w)\right]\real\left[u_{k_{3}}(0)u_{k_{3}}^{*}(w)\right]}{k_{2}k_{3}|u_{k_{2}}(0)|^{2}|u_{k_{3}}(0)|^{2}}+{\rm
2~{}perm.}\Bigg{\\}}.\end{split}$ (4.9)

##  5 Evaluating the bispectrum

For now, we focus on the integral

$\begin{split}\mathcal{I}&=\int^{0}_{-\infty}\differential{w}q^{2}(w)\frac{\imaginary[u_{k_{1}}(0)u^{*}_{k_{1}}(w)]\real[u_{k_{2}}(0)u^{*}_{k_{2}}(w)]\real[u_{k_{3}}(0)u^{*}_{k_{3}}(w)]}{\absolutevalue{u_{k_{2}}(0)}^{2}\absolutevalue{u_{k_{3}}(0)}^{2}},\end{split}$
(5.1)

where we have omitted the permutations and $i\epsilon$ prescription. Taking
the real and imaginary parts gives

$\begin{split}\mathcal{I}&=-i\int^{0}_{-\infty}\differential{w}q^{2}(w)\frac{[u_{k_{1}}(0)u_{k_{1}}^{*}(w)-u_{k_{1}}^{*}(0)u_{k_{1}}(w)][u_{k_{2}}(0)u_{k_{2}}^{*}(w)+u_{k_{2}}^{*}(0)u_{k_{2}}(w)]}{8u_{k_{2}}(0)u_{k_{2}}^{*}(0)u_{k_{3}}(0)u_{k_{3}}^{*}(0)}\\\
&\times[u_{k_{3}}(0)u_{k_{3}}^{*}(w)+u_{k_{3}}^{*}(0)u_{k_{3}}(w)].\end{split}$
(5.2)

We change the domain of integration to $(0,\infty)$ by making explicit that in
conformal time $w=-\absolutevalue{w}$ (from here onwards we drop the absolute
sign). We use the analytic continuation formula (10.11.8) in [23], given by

$\begin{split}H_{n}^{(2)}(ze^{m\pi
i})=(-1)^{mn}\quantity(mH_{n}^{(1)}(z)+(m+1)H_{n}^{(2)}(z)),\end{split}$ (5.3)

which allows us to take the complex conjugates

$\begin{split}H_{\nu}^{(1)}(\bar{z})=\overline{H_{\nu}^{(2)}(z)},\mbox{\quad
and\quad}H_{\nu}^{(2)}(\bar{z})=\overline{H_{\nu}^{(1)}(z)},\end{split}$ (5.4)

for positive argument and real order $\nu$. The mode functions become

$\begin{split}&u_{k}(w)=-\frac{1}{q(w)}\frac{\sqrt{\pi
w}}{2}e^{\displaystyle-\frac{i3\pi}{4}}\quantity(H^{(1)}_{\nu}(kw)+2H^{(2)}_{\nu}(kw)),\\\
&u_{k}^{*}(w)=-\frac{1}{q(w)}\frac{\sqrt{\pi
w}}{2}e^{\displaystyle\frac{i3\pi}{4}}\quantity(H^{(2)}_{\nu}(kw)+2H^{(1)}_{\nu}(kw)),\end{split}$
(5.5)

which results to

$\begin{split}\lim_{w\rightarrow
0}u_{k}(0)=\frac{(-1)^{\frac{3}{4}}}{k\sqrt{\pi}},\end{split}$ (5.6)

after the small argument limit has been considered. Then, using expressions
(5.3)–(5.6), one can simplify the integral (5.2) to

$\begin{split}\mathcal{I}&=-\frac{3k_{2}k_{3}\pi^{2}}{64k_{1}}\int_{0}^{\infty}\differential{w}w^{2}\quantity(H^{(1)}_{\nu}(k_{1}w)+H^{(2)}_{\nu}(k_{1}w))\\\
&\times\quantity(H^{(1)}_{\nu}(k_{2}w)-H^{(2)}_{\nu}(k_{2}w))\quantity(H^{(1)}_{\nu}(k_{3}w)-H^{(2)}_{\nu}(k_{3}w))\\\
&=\frac{3k_{2}k_{3}\pi^{2}}{8k_{1}}\int_{0}^{\infty}\differential{w}w^{2}J_{1}(k_{1}w)Y_{1}(k_{2}w)Y_{1}(k_{3}w).\end{split}$
(5.7)

In the last step, we set the order of the Hankel functions to $\nu=1$ and
expanded in terms of Bessel functions and simplified. This resulted to an
integral over a triple product of Bessel functions. Thus, the 3-point function
is now written as222From here onwards we follow the conventions in [24] and
drop the $4\pi^{2}$ factor from the bispectrum expression.

$\begin{split}B_{\Phi}(k_{1},k_{2},k_{3})&=-(2\pi)^{3}\delta^{3}(\mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3})\quantity(\frac{5}{3})^{3}\frac{1}{(k_{1}k_{2}k_{3})^{3}}\mathcal{P}_{\zeta}^{2}\frac{108}{125}k_{1}^{3}(k_{2}^{2}+k_{3}^{2})\\\
&\times\frac{3k_{2}k_{3}\pi^{2}}{8k_{1}}\int_{0}^{\infty}\differential{w}w^{2}J_{1}(k_{1}w)Y_{1}(k_{2}w)Y_{1}(k_{3}w)+\mbox{\quad
2 perms\quad}.\end{split}$ (5.8)

The integrals involved in our analysis are highly oscillatory and divergent,
causing a great complexity in the analytical and numerical methods of
integration. Therefore, alternative approaches should be used to deal with
this issue. Next, we examine each shape individually, i.e. we independently
address the equilateral, local and isosceles configurations.

###  5.1 The equilateral configuration

To obtain the equilateral shape we let $k_{1}=k_{2}=k_{3}=k$. Then, the
integral in (5.7) simplifies to

$\begin{split}\mathcal{I}_{equil}&=\frac{3k\pi^{2}}{8}\int_{0}^{\infty}\differential{w}w^{2}J_{1}(kw)Y_{1}^{2}(kw),\end{split}$
(5.9)

which can be evaluated exactly, namely
$\mathcal{I}_{equil}=\frac{9+2\sqrt{3}\pi}{18k^{2}}$. Finally, accounting for
all permutations, together with the prefactor involved in (5.8), the
bispectrum in the equilateral configuration is given by the expression

$\begin{split}B_{\Phi}&=-(2\pi)^{3}\delta^{3}(\mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3})\frac{\mathcal{P}_{\zeta}^{2}}{(k_{1}k_{2}k_{3})^{3}}\frac{4}{3}k^{3}\quantity(9+2\sqrt{3}\pi).\end{split}$
(5.10)

We can also verify the above result by solving the integral, numerically, mode
by mode. To this end, we firstly need to regulate the oscillatory nature of
(5.9), at infinity, by multiplying the integrand with
$e^{-\epsilon\absolutevalue{w}^{2}/2}$, where $\epsilon\ll 1$. This has the
effect of suppressing the highly oscillatory behaviour at infinity and forcing
convergence of the numerical integration (see Appendix (A.2)). We find that
for very small values of $\epsilon$, the numerical results are in full
agreement with the analytic results for a large range of comoving scales which
could be observed with current and future CMB experiments.

As we will see later, the equilateral configuration is subdominant since the
bispectrum peaks in the squeezed and folded configurations. Nevertheless,
discussing it in detail allows us to introduce the regulating procedure, due
to the simplicity of the problem, in this case. Next, we will use this method
in place of the $i\epsilon$ prescription, to demonstrate convergence of the
analytic result, for the bispectrum amplitude, in the squeezed limit.

###  5.2 Local configuration

Next, we move on to the local configuration by letting the two long wave modes
to be equal, that is $k_{1}=k_{2}=k$ and take the limit of the remaining mode
to zero, $k_{3}\rightarrow 0$. Similar to what presented in the previous
subsection, there are three permutations of momenta to consider, two of which
coincide owing to the symmetry of the triple product of the Bessel functions
involved in the original integrand (5.7). Therefore, after taking the small
argument approximation of the corresponding Bessel function for $k_{3}\to 0$,
the two integrals to be solved are

$\begin{split}&\mathcal{I}_{sq_{1}}=-\frac{3\pi}{4}\int_{0}^{\infty}\differential{w}wJ_{1}(kw)Y_{1}(kw),\end{split}$
(5.11)

and

$\begin{split}&\mathcal{I}_{sq_{2}}=\frac{3}{16}\pi^{2}k^{2}\int_{0}^{\infty}\differential{w}w^{3}Y_{1}^{2}(kw).\end{split}$
(5.12)

First, let us focus on the integral (5.11). Using the following formula (see
Section 10.22. in [23])

$\begin{split}\int\differential{z}z\mathscr{B}_{\mu}(az)\mathscr{D}_{\mu}(az)=\frac{1}{4}z^{2}\quantity[2\mathscr{B}_{\mu}(az)\mathscr{D}_{\mu}(az)-\mathscr{B}_{\mu+1}(az)\mathscr{D}_{\mu+1}(az)-\mathscr{B}_{\mu+1}(az)\mathscr{D}_{\mu-1}(az)],\end{split}$
(5.13)

allows us to compute the integral (5.11) in terms of its antiderivative. Thus,

$\begin{split}&\mathcal{I}_{sq_{1}}=\frac{3\pi}{4}\left(\lim_{w\to\infty}\mathcal{I}_{sq1}^{\text{anti}}(w)-\mathcal{I}_{sq1}^{\text{anti}}(0)\right),\end{split}$
(5.14)

where

$\begin{split}\mathcal{I}_{sq_{1}}^{\text{anti}}(w)=\frac{1}{4}w^{2}\quantity[J_{2}(kw)Y_{0}(kw)-2J_{1}(kw)Y_{1}(kw)+J_{0}(kw)Y_{2}(kw)].\end{split}$
(5.15)

Finally, the approximations for large and small arguments of Bessel functions
can be used accordingly, such that

$\mathcal{I}_{sq_{1}}=\frac{3}{4k^{2}}-\frac{3}{8k^{2}}\sin(2kw)+\mathcal{O}(w),$
(5.16)

where the second term, involving the sine function, is the contribution from
infinity.

As expected, the last result recovers the oscillations at infinity, which
should be regulated after applying the $i\epsilon$ prescription, leaving
behind the leading order term, that is

$\mathcal{I}_{sq_{1}}=\frac{3}{4k^{2}}+\mathcal{O}(w),\quad\text{ as }w\to 0.$
(5.17)

Alternatively, the oscillatory terms at infinity can be regulated by computing
the result numerically, following the approach discussed in the equilateral
configuration case. That is, to multiply the integrand in (5.11) by a
regulator of the form $e^{-\epsilon\absolutevalue{w}^{2}/2}$, $\epsilon\ll 1$,
which tames the highly oscillatory terms and forces convergence (see Appendix
(A.2)). Indeed, our numerical result is in full agreement with the analytic
result in (5.17).

Finally, including the prefactor and accounting for both permutations, gives

$\begin{split}B_{\phi_{sq_{1}}}&=-(2\pi)^{3}\delta^{3}(\mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3})\frac{1}{(k_{1}k_{2}k_{3})^{3}}\mathcal{P}_{\zeta}^{2}6k^{3}.\end{split}$
(5.18)

Next, we consider the remaining permutation, $\mathcal{I}_{sq_{2}}$. Using the
formula (see section 10.22. in [23])

$\begin{split}&\int\differential{z}z^{\mu+\nu+1}\mathscr{B}_{\mu}(az)\mathscr{D}_{\nu}(az)=\frac{z^{\mu+\nu+2}}{2(\mu+\nu+1)}\quantity[\mathscr{B}_{\mu}(az)\mathscr{D}_{\nu}(az)+\mathscr{B}_{\mu+1}(az)\mathscr{D}_{\nu+1}(az)],\end{split}$
(5.19)

where $\mu+\nu\neq-1$, we find the antiderivative

$\begin{split}\mathcal{I}^{\text{anti}}_{sq_{2}}=\frac{3}{16}\pi^{2}k^{2}\times\frac{1}{6}w^{4}\quantity(Y^{2}_{1}(kw)+Y^{2}_{2}(kw)).\end{split}$
(5.20)

Likewise to our analysis above, the appropriate approximations of the Bessel
functions can be used, along with a regulation process, resulting to

$\begin{split}\mathcal{I}_{sq_{2}}=\frac{1}{2k^{2}}+\mathcal{O}(w),\quad\text{
as }w\to 0.\end{split}$ (5.21)

Putting it together with the prefactor (it is customary to ignore the
$(k_{1}k_{2}k_{3})^{-3}$ term in the bispectrum amplitude calculations when
taking the squeezed limit, $k_{3}\rightarrow 0$) gives a vanishing result in
the limit $k_{3}\rightarrow 0$

$\begin{split}&\lim_{k_{3}\rightarrow
0}B_{\phi_{sq_{2}}}\simeq\frac{108}{125}k_{3}^{3}\bigg{|}_{\lim
k_{3}\rightarrow 0}=0.\end{split}$ (5.22)

In other words, the result that contributes to the bispectrum amplitude, in
the squeezed limit, is the one obtained from the first two permutations in
(5.18).

Note that again we found the main contribution to the bispectrum comes from
approximating the integral near the lower limit. This can be also understood
in another way. The subhorizon modes are oscillating and cancel out, leaving
only the dominant contribution for $0<w\lesssim w_{\Lambda}$ where
$w_{\Lambda}$ is some cut-off. In general, the cutoff $w_{\Lambda}\sim 1/k$
gives results that are consistent with analytic calculations, but as we shall
see in Section (6), this fails in the flattened limit where there is resonance
between the modes (i.e. oscillations don’t cancel) and thus a physical cutoff
will be needed to regulate the integrals.

###  5.3 Isosceles configuration

The last case to be discussed is the isosceles configuration. In principle,
this is the most general case, from which one can derive all the other
standard shapes. Therefore, solving for the isosceles configuration will not
only give a more general result, but also verify all the previous findings,
i.e those for the equilateral and local configurations.

By conservation of momentum (ensured by the argument of the $\delta$-function)
we have three permutations of the momenta of the closed triangle
$k_{1},k_{2},k_{3}$. As before, two of them coincide due to the symmetry of
the triple product of the Bessel functions involved in the original integrand
(5.7). That is, setting $k_{1}=k_{2}=k$, the integrals to solve are

$\begin{split}\mathcal{I}_{iso_{1}}&=\frac{3k_{3}\pi^{2}}{8}\int_{0}^{\infty}\differential{w}w^{2}J_{1}(kw)Y_{1}(kw)Y_{1}(k_{3}w),\end{split}$
(5.23)

and

$\begin{split}\mathcal{I}_{iso_{2}}&=\frac{3k^{2}\pi^{2}}{8k_{3}}\int_{0}^{\infty}\differential{w}w^{2}J_{1}(k_{3}w)Y_{1}^{2}(kw),\end{split}$
(5.24)

where the integral (5.23) is repeated twice when all three permutations are
summed. The integrals (5.23)–(5.24) can be expressed in terms of the very
general Meijer G-functions (see Appendix in (A.1)) and then be evaluated as

$\begin{split}\mathcal{I}_{iso_{1}}&=-\frac{3\pi^{3/2}}{2k_{3}^{2}}\text{MeijerG}\left[\Big{\\{}\\{0\\},\\{1/2,1/2\\}\Big{\\}},\Big{\\{}\\{0,1\\},\\{1/2\\}\Big{\\}},\frac{4k^{2}}{k_{3}^{2}}\right]\\\
&=\frac{3}{2k_{3}}\times\left[\left(\frac{1}{4k^{2}-k_{3}^{2}}\right)\left(k3+\frac{4k\text{
arccsc}\left(\frac{2k}{k_{3}}\right)}{\sqrt{4-\frac{k3^{2}}{k^{2}}}}\right)\right],\end{split}$
(5.25)

and

$\begin{split}\mathcal{I}_{iso_{2}}&=\frac{3k^{2}\pi^{3/2}}{k_{3}^{4}}\text{MeijerG}\left[\Big{\\{}\\{-1\\},\\{-1/2,1/2\\}\Big{\\}},\Big{\\{}\\{-1,1\\},\\{-1/2\\}\Big{\\}},\frac{4k^{2}}{k_{1}^{2}}\right]\\\
&=\frac{6k^{2}k_{3}-3k_{3}^{3}-\frac{24k^{3}\text{arccsc}\left(\frac{2k}{k_{3}}\right)}{\sqrt{4-\frac{k_{3}^{2}}{k^{2}}}}}{8k^{2}k_{3}^{3}-2k_{3}^{5}},\end{split}$
(5.26)

where both (5.25) and (5.26) are valid for $k_{3}<2k$. Finally, the total
result for the bispectrum in the isosceles limit, including the correct
prefactors for each permutation, is

$\begin{split}B_{\Phi_{iso}}&=-(2\pi)^{3}\delta^{(3)}(\mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3})P_{\zeta}^{2}\frac{1}{(k_{1}k_{2}k_{3})^{3}}\frac{1}{k_{3}\quantity(4k^{2}-k_{3}^{2})^{\frac{3}{2}}}\\\
&\times\quantity[12k^{2}\quantity(k_{3}\sqrt{4k^{2}-k_{3}^{2}}\quantity(k^{3}+2k^{2}k_{3}+kk_{3}^{2}-k_{3}^{3})+4k^{3}(k-k_{3})^{2}\arccsc\quantity(\frac{2k}{k_{3}}))]\end{split}$
(5.27)

From the last expression we can derive the squeezed limit by letting $k_{3}\to
0$, as well as any other shapes in between. In particular, we will find the
bi-thermal bispectrum peaks in the flattened limit, that is the limit when
$k_{3}\to 2k$. Indeed, one can verify that by taking $k_{3}\rightarrow 0$ in
(5.27), we obtain the local configuration expression

$\begin{split}B_{\Phi_{sq}}&=-(2\pi)^{3}\delta^{(3)}(\mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3})\frac{1}{(k_{1}k_{2}k_{3})^{3}}P_{\zeta}^{2}6k^{3},\end{split}$
(5.28)

which fully agrees with our result in (5.18). Note that adopting this
approach, which involves Meijer G-functions and their general definition via a
line integral in the complex plane (see Appendix (A.1)) does not necessitate
the use of a regulation process at infinity, since this is dealt with the
complex analysis arguments. Also, we find that the second permutation, given
in (5.26), does not contribute to the bispectrum amplitude in this limit,
which is in agreement with our result in (5.22).

##  6 Bispectrum results

Let us first recall the following formulae for calculating the amplitude and
shape of the bispectrum [24]. The bispectrum is defined, as

$\begin{split}\langle\Phi(\mathbf{k}_{1})\Phi(\mathbf{k}_{2})\Phi(\mathbf{k}_{3})\rangle=(2\pi)^{3}\delta^{(3)}(\mathbf{k}_{1}+\mathbf{k}_{1}+\mathbf{k}_{1})B_{\Phi}(k_{1},k_{2},k_{3}),\end{split}$
(6.1)

where $\Phi=\frac{3}{5}\zeta$. Assuming scale invariant statistics, the shape
function is defined, as

$\begin{split}S(k_{1},k_{2},k_{3})=\frac{(k_{1}k_{2}k_{3})^{2}}{A^{2}}B_{\Phi}(k_{1},k_{2},k_{3}),\end{split}$
(6.2)

where $A$ is the dimensionless power spectrum. We first look at the
equilateral limit. Using our result in (5.10), we see that in the equilateral
limit the bispectrum goes like $B_{\phi}\sim k^{-6}$. Therefore, we can use
the standard template

$\begin{split}B_{\Phi}(k_{1},k_{2},k_{3})&=6A^{2}f_{NL}^{\text{equil}}\bigg{[}-\quantity(\frac{1}{k_{1}^{3}k_{2}^{3}}+\frac{1}{k_{2}^{3}k_{3}^{3}}+\frac{1}{k_{3}^{3}k_{1}^{3}})-\frac{2}{(k_{1}k_{2}k_{3})^{2}}+\bigg{(}\frac{1}{k_{1}k_{2}^{2}k_{3}^{3}}+\frac{1}{k_{1}k_{3}^{2}k_{2}^{3}}\\\
&+\frac{1}{k_{2}k_{1}^{2}k_{3}^{3}}+\frac{1}{k_{2}k_{3}^{2}k_{1}^{3}}+\frac{1}{k_{3}k_{1}^{2}k_{2}^{3}}+\frac{1}{k_{3}k_{2}^{2}k_{1}^{3}}\bigg{)}\bigg{]}.\end{split}$
(6.3)

Setting $k_{1}=k_{2}=k_{3}=k$ in the above, gives
$B_{\Phi}(k_{1},k_{2},k_{3})=6A^{2}f_{NL}^{\text{equil}}k^{-6}$. Identifying
$A$ with $P_{\zeta}$, we find the bispectrum amplitude in the equilateral
limit, is given by 333The negative sign in the bispectrum amplitudes simply
indicates negative correlation of the fluctuations.

$\begin{split}f_{NL}^{\text{equil}}&=-\frac{2}{9}\quantity(9+2\sqrt{3}\pi).\end{split}$
(6.4)

Using (6.2), the normalised shape in the equilateral configuration, takes the
form

$\begin{split}\frac{S^{\text{equil}}(x_{2},x_{3})}{S^{\text{equil}}(1,1,1)}=\bigg{[}-2+\frac{1}{x_{2}}+x_{2}+\frac{1}{x_{3}}-\frac{1}{x_{2}x_{3}}+\frac{x_{2}}{x_{3}}-\frac{x_{2}^{2}}{x_{3}}+x_{3}+\frac{x_{3}}{x_{2}}-\frac{x_{3}^{2}}{x_{2}}\bigg{]},\end{split}$
(6.5)

where $x_{2}=k_{2}/k_{1}$ and $x_{3}=k_{3}/k_{1}$. This is plotted in Fig.
(2).

Next, using (5.18), we find that in the squeezed limit the bispectrum goes
like $B_{\Phi}\sim k^{3}$. Therefore, we can use the standard template

$\begin{split}B_{\Phi}^{\text{local}}(k_{1},k_{2},k_{3})&=2A^{2}f_{NL}^{\text{local}}\quantity(\frac{1}{k_{1}^{3}k_{2}^{3}}+\frac{1}{k_{2}^{3}k_{3}^{3}}+\frac{1}{k_{1}^{3}k_{3}^{3}}).\end{split}$
(6.6)

Taking $k_{3}\rightarrow 0,k_{2}\sim k_{1}$ of the above, gives the bispectrum
in the squeezed limit. Comparing this to our result in (5.18), we find

$\begin{split}f_{NL}^{\rm local}=-\frac{3}{2},\end{split}$ (6.7)

which is well within experimental bounds ([17]). Using (6.2), the local shape
is given by

$\begin{split}S^{\text{local}}(x_{2},x_{3})&=\frac{1}{3}f_{NL}^{\text{local}}\quantity(\frac{1}{x_{2}x_{3}}+\frac{x_{2}^{2}}{x_{3}}+\frac{x_{3}^{2}}{x_{2}}),\end{split}$
(6.8)

where we ensure this is normalised for $S(1,1,1,)=1$. This gives the familiar
shape in Fig. (2).

(a) Equilateral

(b) Local

Figure 2: Plot of the shape function in the equilateral (equation 6.5) and
local (Equation 6.9) configuration (normalised).

So far, we have found that in the equilateral and squeezed limits the
bispectrum takes a standard form. This enabled us to obtain values for the
bispectrum amplitude $f_{NL}$ in those limits. Next, we examine the general
form of the expression in (5.27).

Before we start, a few words are in order. Looking at the denominator in
(5.27), we find that there is a divergence at the limit $k_{3}\to 2k$. Now, if
we naively restore the average of the $k$ momenta this becomes $k_{3}\to
k_{1}+k_{2}$, which is the flattened limit of the bispectrum. This effect is
not entirely surprising. Such behaviour has been encountered before in [25],
when non-Bunch-Davies vacua were considered, as well in [26] they found that
the bispectrum peaks in the flattened limit when starting from excited initial
states (see also [27, 28]). As we discussed above, the reason for this
divergence is that excited modes, e.g., with acoustic frequencies $k_{1}$ and
$k_{2}$ are in resonance with $k_{3}$, if $k_{3}\simeq k_{1}+k_{2}$ and thus
can excite $k_{3}$ particles in long before horizon crossing. If the resonance
is not perfect, we expect the number of $k_{3}$ particles excited to scale as
an inverse power of $k_{1}+k_{2}-k_{3}$. Alternatively, the divergence can be
limited by the finite time that the modes spend in the thermal bath. Indeed,
the early thermal vacuum is understood as a collection of excited states in
equilibrium while the sound horizon is taken to be infinite at $t\to 0$ (or
$y\to-\infty$). If we want to be more realistic we may expect a hard cut-off
or regulator, corresponding to a quantum gravity phase where the bi-thermal
classical description is no longer valid. Solving numerically, this ensures a
finite result in the folded limit for a wide range of modes (in agreement with
the conclusions in [25]).

A possibly more surprising feature of our model is the appearance of local (or
squeezed) non-Gaussianity, usually expected in multi-field models, even though
we only have a single scalar field. However, we note the temperature of the
field, effectively, acts as a second clock in this model, which allows for
coupling of short and long wavelength modes.

We proceed to examine the shape of (5.27). We can naively restore the momenta
as $k=1/2(k_{1}+k_{2})$. Using this and substituting for $x_{2}=k_{2}/k_{1}$
and $x_{3}=k_{3}/k_{1}$ in the standard expression for the bispectrum shape
(6.2), we obtain

$\begin{split}S_{\text{bi-
thermal}}(x_{2},x_{3})&=-\frac{1}{8x_{2}x_{3}^{2}\quantity(1+x_{2}-x_{3})^{\frac{3}{2}}\quantity(1+x_{2}+x_{3})^{\frac{3}{2}}}\Bigg{[}3\quantity(1+x_{2})^{2}\\\
&\times\Bigg{(}x_{3}\sqrt{\quantity(1+x_{2}-x_{3})\quantity(1+x_{2}+x_{3})}\bigg{(}\quantity(1+x_{2})^{3}+4\quantity(1+x_{2})^{2}x_{3}\\\
&+4\quantity(1+x_{2})x_{3}^{2}-8x_{3}^{3}\bigg{)}+\quantity(1+x_{2})^{3}\quantity(1+x_{2}-2x_{3})^{2}\arccsc\quantity(\frac{1+x_{2}}{x_{3}})\Bigg{)}\Bigg{]}.\end{split}$
(6.9)

The normalised shape is plotted in Fig. (3). We see that the bispectrum peaks
in the folded and squeezed limit, as expected. Please note, this is not a
template for general $x_{2},x_{3}$ but only valid in the limit
$x_{3}\rightarrow 0$ and on the line $x_{3}=2x_{2}$.

Finally, setting $x_{2}=1$ and expanding around $x_{3}\sim 2$, we find that in
the flattened limit the bispectrum shape goes as
$\propto(k_{1}+k_{2}-k_{3})^{-3/2}$.

Figure 3: Our final bispectrum shape function for the bi-thermal model
(equation 6.9), extrapolated from exact isosceles result. The bispectrum peaks
at the folded and squeezed limit.

##  7 Conclusions

The bi-thermal model of [11] is the pinnacle of a proud lineage of models
playing with the speed of light [29], most notably rendering the speed of
massless matter particles different from the speed of gravity. Starting from
[30, 31, 32] these models became progressively better defined [33] and more
concrete and solid in their predictions [34, 6, 15], eventually matching
inflation. But it was not until [11] that this general class of models met
with the prospect of potentially supplanting inflation in predictivity [13].
This claim rests on past work on the power spectrum of the fluctuations. To
this we have now added predictions involving the 3-point function.

Specifically we found the following results:

  * •

Squeezed non-Gaussianity with $f^{\rm local}_{NL}=-\frac{3}{2}$, which is
unique for single field scenarios.

  * •

Equilateral non-Gaussianity with $f^{\rm
equil}_{NL}=-\frac{2}{9}\quantity(9+2\sqrt{3}\pi)$.

  * •

We obtained a new shape that diverges as $(k_{1}+k_{2}-k_{3})^{-\frac{3}{2}}$
in the flattened limit.

These results are quite spectacular, in that they leave no leeway for fiddling
with parameters, as they are distinctive and unique to this model. Hence, not
only is the bi-thermal model more predictive than inflation, but it is also
capable of making a prediction which cannot be manufactured within inflation,
without significant fine-tuning.

Although, it is true that current data can only, at most, disprove different
classes of models, we hope that the unique qualities of the bi-thermal model
combined with its rich non-Gaussian footprint will help, in the future, to
differentiate it from other models of inflation.

Finally, we would like to add that the presence of a divergence in the
flattened limit necessitates the application of a physical cutoff. We argue
that this divergence is simply an artefact of the effective field theory.
Indeed, there should be extra corrections coming from the effective field
theory of the DBI action (or the geometric picture). It would be interesting
to see how the flattened limit is modified in that case. Ideally, the UV
completion of the bi-thermal model should, almost certainly, be able to
regulate this divergence.

There are currently weak observational constraints for non-Gaussianities in
the flattened limit. The observability of this signal is something interesting
we should study in the future.

##  8 Acknowledgments

We would like to thank Johannes Noller for initial collaboration and Jinn-Ouk
Gong for discussions and advice in relation to this paper. This work was
supported by the STFC Consolidated Grant ST/L00044X/1 (JM), and the National
Research Foundation of Korea Grant 2019R1A2C2085023 (Maria M), the University
of Waterloo, Natural Sciences and Engineering Research Council of Canada
(NSERC) and the Perimeter Institute for Theoretical Physics (NA). Research at
Perimeter Institute is supported in part by the Government of Canada through
the Department of Innovation, Science and Economic Development Canada and by
the Province of Ontario through the Ministry of Colleges and Universities.

##  Appendix A Appendix

###  A.1 Analytic continuation of Meijer G-functions

A Meijer G-function is defined as follows

$\displaystyle G_{p,q}^{m,n}\left(z\left|\,\begin{matrix}a_{1}&\cdots&a_{p}\\\
b_{1}&\cdots&b_{q}\end{matrix}\right.\right)=$
$\displaystyle\text{MeijerG}\left[\Big{\\{}\\{a_{1},\dots
a_{n}\\},\\{a_{n+1},\dots a_{p}\\}\Big{\\}},\Big{\\{}\\{b_{1},\dots
b_{m}\\},\\{b_{m+1},\dots b_{q}\\}\Big{\\}},z\right]$
$\displaystyle=\frac{1}{2\pi i}$
$\displaystyle\int_{\gamma}\frac{\Gamma(1-a_{1}-s)\cdots\Gamma(1-a_{n}-s)\Gamma(b_{1}+s)\cdots\Gamma(b_{m}+s)}{\Gamma(a_{n+1}+s)\cdots\Gamma(a_{p}+s)\Gamma(1-b_{m+1}-s)\cdots\Gamma(1-b_{q}-s)}\,z^{-s}\text{d}s,$
(A.1)

where the contour $\gamma$ is illustrated in Fig. (4).

Figure 4:  The contour of a Meijer G-function lies between the poles of
$\Gamma(1-a_{i}-s)$ and the poles of $\Gamma(b_{i}+s)$.

In many special cases, Meijer G-function is converted to other known
functions, as it can be seen from the evaluation of the Bessel integrals
below. The contour integral can be viewed as an _inverse Mellin transform_ and
can be treated as such. One can find the poles of the integrand and then by a
_residue analysis_ can evaluate the integral as an infinite series. There
should be paid attention to the domain of convergence, since depending on the
values of $z$, only some of the poles will be considered in the integration
process.

The integral in (5.23) evaluates to

$\displaystyle\mathcal{I}_{iso_{1}}(k,k_{3})=$
$\displaystyle\frac{3k_{3}\pi^{2}}{8}\int_{0}^{\infty}\text{d}w\,\,w^{2}J_{1}(kw)Y_{1}(kw)Y_{1}(k_{3}w)$
$\displaystyle=$
$\displaystyle-\frac{3\pi^{3/2}}{2k_{3}^{2}}\text{MeijerG}\left[\Big{\\{}\\{0\\},\\{1/2,1/2\\}\Big{\\}},\Big{\\{}\\{0,1\\},\\{1/2\\}\Big{\\}},\frac{4k^{2}}{k_{3}^{2}}\right],$
(A.2)

which utilizes the Meijer G-function, as defined in (A.1), for $n=1$, $p=3$,
$m=2$, $q=3$ and $z=\frac{4k^{2}}{k_{3}^{2}}$, with $a_{1}=a_{n}=0$,
$a_{2}=a_{n+1}=1/2$, $a_{3}=a_{p}=1/2$ and $b_{1}=0$, $b_{2}=b_{m}=1$,
$b_{3}=b_{m+1}=b_{q}=1/2$. For the simplification from the integral form to
the Meijer G representation we used [35]. Writing expression (A.2) in its line
integral form, the original integral for $\mathcal{I}_{iso_{1}}$ is now
written as

$\displaystyle\mathcal{I}_{iso_{1}}(k,k_{3})=$
$\displaystyle-\frac{3\pi^{3/2}}{2k_{3}^{2}}\frac{1}{2\pi
i}\int_{\gamma}\frac{\Gamma(1-s)\Gamma(s)\Gamma(1+s)}{\Gamma(\frac{1}{2}+s)\Gamma(\frac{1}{2}+s)\Gamma(\frac{1}{2}-s))}\,\left(\frac{4k^{2}}{k_{3}^{2}}\right)^{-s}\text{d}s,$
(A.3) $\displaystyle=$
$\displaystyle-\frac{3\pi^{3/2}}{2k_{3}^{2}}\frac{1}{2\pi
i}\int_{\gamma}\frac{\cos(\pi s)}{\sin(\pi
s)}\frac{\Gamma(1+s)}{\Gamma(\frac{1}{2}+s)}\left(\frac{4k^{2}}{k_{3}^{2}}\right)^{-s}\text{d}s,$
(A.4)

where the contour $\gamma$ can be seen in the left of Fig. (5). In the same
plot, the poles of the nominator in (A.3) can be seen. Note that for
$2k=k_{3}$ the above integral does not converge. Therefore, we should proceed
under the condition that either $2k<k_{3}$ or $2k>k_{3}$, such that
convergence is achieved. If $2k<k_{3}$, then the integral converges for
$\real(s)<0$ and therefore we should deform the contour $\gamma$ to the left
half $s$-plane. On the other hand, if $2k>k_{3}$, then the integral converges
for $\real(s)>0$ and therefore we should deform the contour $\gamma$ to the
right half $s$-plane. Since we are interested in the bispectrum shapes for
which $2k>k_{3}$, we will deform our original contour $\gamma$ to the right,
where $\real(s)>0$ (see Fig. (5)).

(a)

(b)

Figure 5: (a) An example of a phase-portrait. Here, we illustrate the phase-
portrait of the complex function f(z) = 1/z (simple pole at z). The colour
variation shows how the phase continuously varies while encircling a simple
pole, i.e. red indicates negative real values ($\pm\pi$ phase), cyan the
positive real values (zero phase), etc (b) A phase-portrait of the product of
$\Gamma(1-s)\Gamma(s)\Gamma(1+s)\left(4k^{2}/k_{3}^{2}\right)^{-s}$ can be
seen (nominator in (A.3)). Here, $\frac{2k}{k_{3}}=1.5$, i.e. $2k>k_{3}$. The
solid black line represents the original $\gamma$ contour and the dashed one
the deformed contour. Note that in the negative real $z$-axis the non-
convergence effect is clearly visible. Here, we used the complex variable $z$
instead of $s$.

Now, let’s apply a standard residue analysis to the above contour integral,
under the assumption that $2k>k_{3}$.

Let

$f(s)=\frac{\cos(\pi s)}{\sin(\pi
s)}\frac{\Gamma(1+s)}{\Gamma(\frac{1}{2}+s)}\left(\frac{4k^{2}}{k_{3}^{2}}\right)^{-s}.$

Then, according to the _Residue Theorem_ from complex analysis, we have that

$\frac{1}{2\pi i}\int_{\gamma}f(s)\text{d}s=\frac{1}{2\pi
i}\int_{\gamma^{+}}f(s)\text{d}s=-\sum_{j=1}^{\infty}\text{Res}(f,s_{j}),$

where $s_{j}$, $j=1,2,3,\dots$ are the poles which are enclosed from the
contour $\gamma^{+}$ and

$\text{Res}(f,s_{j})=\lim_{s\to s_{j}}(s-s_{j})f(s).$ (A.5)

In our case the poles are known exactly, namely $s_{j}=j$, where
$j=1,2,3,\dots$. Applying L’Hospital’s rule in order to compute the limit
(A.5), one finds that the contour integral (A.4) is given from an infinite
series which can be evaluated exactly. That is,

$\displaystyle\mathcal{I}_{iso_{1}}(k,k_{3})$
$\displaystyle=-\frac{3\pi^{3/2}}{2k_{3}^{2}}\times\frac{1}{2\pi
i}\int_{\gamma}f(s)\text{d}s$
$\displaystyle=-\frac{3\pi^{3/2}}{2k_{3}^{2}}\times\left[-\sum_{j=1}^{\infty}\frac{\Gamma(1+j)}{\pi\Gamma(\frac{1}{2}+j)}\left(\frac{2k}{k_{3}}\right)^{-2j}\right]$
$\displaystyle=-\frac{3\pi^{3/2}}{2k_{3}^{2}}\times\left[-\left(\frac{k_{3}}{\pi^{3/2}(4k^{2}-k_{3}^{2})}\right)\left(k3+\frac{4k\text{
arccsc}\left(\frac{2k}{k_{3}}\right)}{\sqrt{4-\frac{k3^{2}}{k^{2}}}}\right)\right]$
$\displaystyle=\frac{3}{2k_{3}}\times\left[\left(\frac{1}{4k^{2}-k_{3}^{2}}\right)\left(k3+\frac{4k\text{
arccsc}\left(\frac{2k}{k_{3}}\right)}{\sqrt{4-\frac{k3^{2}}{k^{2}}}}\right)\right].$
(A.6)

Finally, in the _sqeezed limit_ , which corresponds to the local
configuration, equation (A.6) results to

$\mathcal{I}_{iso_{1}}=\frac{3}{4k^{2}}+\mathcal{O}(k_{3}),\quad\text{ as
}k_{3}\to 0,$ (A.7)

while for the _equilateral case_ we find

$\mathcal{I}_{iso_{1}}=\frac{9+2\sqrt{3}\pi}{18k^{2}}+\mathcal{O}(k_{3}-k),\quad\text{
as }k_{3}\to k.$ (A.8)

The next integral to solve is

$\displaystyle\mathcal{I}_{iso_{2}}(k,k_{3})=$
$\displaystyle\frac{3k^{2}\pi^{2}}{8k_{3}}\int_{0}^{\infty}\text{d}w\,\,w^{2}J_{1}(k_{3}w)(Y_{1}(kw))^{2}$
$\displaystyle=$
$\displaystyle\frac{3k^{2}\pi^{3/2}}{k_{3}^{4}}\text{MeijerG}\left[\Big{\\{}\\{-1\\},\\{-1/2,1/2\\}\Big{\\}},\Big{\\{}\\{-1,1\\},\\{-1/2\\}\Big{\\}},\frac{4k^{2}}{k_{3}^{2}}\right],$

which again utilizes the Meijer G-function, as defined in (A.1), for $n=1$,
$p=3$, $m=2$, $q=3$ and $z=\frac{4k^{2}}{k_{3}^{2}}$, with $a_{1}=a_{n}=-1$,
$a_{2}=a_{n+1}=-1/2$, $a_{3}=a_{p}=1/2$ and $b_{1}=-1$, $b_{2}=b_{m}=1$,
$b_{3}=b_{m+1}=b_{q}=-1/2$. Then,

$\displaystyle\mathcal{I}_{iso_{2}}(k,k_{3})=$
$\displaystyle\frac{3k^{2}\pi^{3/2}}{k_{3}^{4}}\frac{1}{2\pi
i}\int_{\gamma}\frac{\Gamma(2-s)\Gamma(s-1)\Gamma(1+s)}{\Gamma(s-\frac{1}{2})\Gamma(\frac{1}{2}+s)\Gamma(\frac{3}{2}-s))}\,\left(\frac{4k^{2}}{k_{3}^{2}}\right)^{-s}\text{d}s,$
(A.9) $\displaystyle=$
$\displaystyle\frac{3k^{2}\pi^{3/2}}{k_{3}^{4}}\frac{1}{2\pi
i}\int_{\gamma}\frac{\cos(\pi s)}{\sin(\pi
s)}\frac{\Gamma(1+s)}{\Gamma(\frac{1}{2}+s)}\left(\frac{4k^{2}}{k_{3}^{2}}\right)^{-s}\text{d}s,$
(A.10)

where the contour $\gamma$ is similar to the one illustrated in Fig. (5), but
shifted to the right of the pole $s_{j}=1$. As before, for $2k=k_{3}$ the
above integral does not converge, which suggests that either $2k<k_{3}$ or
$2k>k_{3}$.

Now, let’s apply a standard residue analysis to the above contour integral,
under the assumption that $2k>k_{3}$. We will proceed similarly with the
analysis of the integral (A.4). Let

$f(s)=\frac{\cos(\pi s)}{\sin(\pi
s)}\frac{\Gamma(1+s)}{\Gamma(\frac{1}{2}+s)}\left(\frac{4k^{2}}{k_{3}^{2}}\right)^{-s}$

and apply the _Residue Theorem_

$\frac{1}{2\pi i}\int_{\gamma}=-\sum_{j=1}^{\infty}\text{Res}(f,s_{j}),$

where $s_{j}$, ($j=1,2,3,\dots$) are the poles which are enclosed from the
contour $\gamma^{+}$ and

$\text{Res}(f,s_{j})=\lim_{s\to s_{j}}(s-s_{j})f(s).$ (A.11)

Likewise with what presented before, the location of the poles is known
exactly, namely $s_{j}=j$, where $j=2,3,4\dots$ 444Note that the first pole to
consider is $s_{2}=2$. This is because of the position of the contour $\gamma$
which is slightly shifted to the right, such that it doesn’t include the pole
$s_{1}=1$.. Applying L’Hospital’s rule in order to compute the limit (A.11),
one finds that the contour integral (A.10) is given from an infinite series
which then can be evaluated exactly. That is,

$\displaystyle\mathcal{I}_{iso_{2}}(k,k_{3})$
$\displaystyle=\frac{3k^{2}\pi^{3/2}}{k_{3}^{4}}\times\frac{1}{2\pi
i}\int_{\gamma}f(s)\text{d}s$
$\displaystyle=\frac{3k^{2}\pi^{3/2}}{k_{3}^{4}}\times\left[-\sum_{j=2}^{\infty}\frac{\Gamma(1+j)}{\pi\Gamma(\frac{1}{2}+j)}\left(\frac{2k}{k_{3}}\right)^{-2j}\right]$
$\displaystyle=\frac{6k^{2}k_{3}-3k_{3}^{3}-\frac{24k^{3}\text{arccsc}\left(\frac{2k}{k_{3}}\right)}{\sqrt{4-\frac{k_{3}^{2}}{k^{2}}}}}{8k^{2}k_{3}^{3}-2k_{3}^{5}}.$
(A.12)

As we mentioned earlier, all other shapes can be derived from the isosceles
configuration by considering the appropriate limits. To retrieve the local
configuration we just need to introduce the limiting procedure for $k_{3}$,
while for the equilateral case one should consider $k_{3}\to k$. One should be
careful, since during the limiting process for the equilateral case in the
second permutation $I_{iso_{2}}$, an extra term should be added; a term
correcting the complex argument of our result and ensuring that the function
maintains its analytical properties. This becomes more clear in Fig. (6). Let
us expand on this argument.

(a)

(b)

(c)

Figure 6: Phase-portraits of the results (a) (A.6), (b) (A.12) and (c) (A.13).
Here, the complex variable is $z=k$.

Recall that in the equilateral limit both permutations $I_{iso_{1}}$ and
$I_{iso_{2}}$ should retrieve the same result. Yet, this is not the case here.
Comparing the phase-portraits of results (A.6) and (A.12) when $k_{3}\to k$
((a) and (b) in Fig. (6)) one observes that the two expressions are not
identical, but rather differ in their complex argument. This can be corrected
by adding the extra term

$\mathcal{I}_{extra}(k,k_{3})=\frac{3\pi
k^{4}}{k_{3}^{3}(4k^{2}-k_{3}^{2})^{3/2}}.$

in (A.12). Then,

$\displaystyle\hat{\mathcal{I}}_{iso_{2}}(k,k_{3})$
$\displaystyle=\frac{6k^{2}k_{3}-3k_{3}^{3}-\frac{24k^{3}\text{arccsc}\left(\frac{2k}{k_{3}}\right)}{\sqrt{4-\frac{k_{3}^{2}}{k^{2}}}}}{8k^{2}k_{3}^{3}-2k_{3}^{5}}+\mathcal{I}_{extra}(k,k_{3}),$
(A.13)

retrieves the correct limit in the right half complex plane, which is made
clear by comparing the phase-portraits (a) and (c) in Fig. (6). Note that we
are only interested in the right half–plane, which is the domain where
$\real(k)>0$.

Finally, in the limiting cases of interest one obtains

$\mathcal{I}_{iso_{2}}=-\frac{1}{2k^{2}}+\mathcal{O}(k_{3}),\quad\text{ as
}k_{3}\to 0,$ (A.14)

and

$\mathcal{I}_{iso_{2}}=\frac{9+2\sqrt{3}\pi}{18k^{2}}+\mathcal{O}(k_{3}-k),\quad\text{
as }k_{3}\to k,$ (A.15)

where in the last result we have dropped the hat notation.

###  A.2 Regulating oscillatory integrals

The integrals in (5.7) are highly oscillatory and suffer from convergence
issues. A regulator, or damping factor, can force decaying behaviour which
allows us to compute the integral numerically. Take, for example, the integral
in the equilateral configuration, given in (5.9). We multiply the integrand by
a term of the form $e^{-\epsilon\absolutevalue{w}^{2}/2}$, where $\epsilon\ll
1$

$\begin{split}\mathcal{I}_{equil}&=\frac{3k\pi^{2}}{8}\int_{0}^{\infty}\differential{w}w^{2}J_{1}(kw)Y_{1}^{2}(kw)e^{-\frac{\epsilon\absolutevalue{w}^{2}}{2}}.\end{split}$
(A.16)

As seen in Fig. (7(a)), the integrand is highly oscillatory with increasingly
large amplitudes as $w$ approaches infinity. Using the regulator, forces the
integrand to decay at infinity as seen in Fig. (7(b)). This helps the integral
to converge at infinity and therefore, to numerically compute it, giving full
agreement with the analytic result.

Similarly, the oscillatory behaviour of the integral in the squeezed limit,
given in (5.11), can be regulated at infinity, using the same method, as seen
in Figure (8). This ensures convergence for $w\rightarrow\infty$ in the
numerical computation.

The numerical result becomes more accurate as $\epsilon$ takes smaller values
and it is valid for a large range of comoving scales seen by CMB experiments.

Unfortunately, the regulator is not effective on all cases considered here.
Instead, we can use a hard cut-off $w_{\Lambda}\sim 1/k$. (see Fig. (9)). For
suitable choice of $w_{\Lambda}$, this successfully approximates the analytic
result at $w\rightarrow 0$ in (5.21), for a wide range of modes.

(a) Without regulator

(b) With regulator

Figure 7: The integrand in the equilateral configuration can be regulated at
infinity. In Fig. (a) we plot the integrand without a regulator, while in Fig.
(b) we plot it with a regulator.

(a) Without regulator

(b) With regulator

Figure 8: The integrand in the squeezed configuration, for the first two
permutations, can be regulated at infinity. In Fig. (a) we plot the integrand
without a regulator, while in Fig. (b) we plot it with a regulator.

(a) Without cut-off

(b) With cut-off

Figure 9: The integrand in the squeezed configuration, for the last
permutation, can be regulated at infinity by using a hard cutoff. In Fig. (a)
we plot the integrand without a hard cutoff, while in Fig. (b) we plot it with
a hard cutoff.

## References

  * [1] A. H. Guth, _Inflationary universe: A possible solution to the horizon and flatness problems_ , _Phys. Rev. D_ 23 (Jan, 1981) 347–356.
  * [2] N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini et al., _Planck 2018 results_ , _Astronomy & Astrophysics_ 641 (Sep, 2020) A6.
  * [3] J. Magueijo and L. Pogosian, _Could thermal fluctuations seed cosmic structure?_ , _Physical Review D_ 67 (Feb, 2003) .
  * [4] P. Ferreira and J. Magueijo, _Observing the temperature of the Big Bang through large scale structure_ , _Phys. Rev. D_ 78 (2008) 061301, [0708.0429].
  * [5] T. Biswas, R. Brandenberger, T. Koivisto and A. Mazumdar, _Cosmological perturbations from statistical thermal fluctuations_ , _Phys. Rev. D_ 88 023517, [1302.6463].
  * [6] J. Magueijo, _Speedy sound and cosmic structure_ , _Physical Review Letters_ 100 (Jun, 2008) .
  * [7] D. Bessada, W. H. Kinney, D. Stojkovic and J. Wang, _Tachyacoustic cosmology: An alternative to inflation_ , _Physical Review D_ 81 (Feb, 2010) .
  * [8] A. Agarwal and N. Afshordi, _Thermal tachyacoustic cosmology_ , _Physical Review D_ 90 (Aug, 2014) .
  * [9] CMB-S4 collaboration, K. N. Abazajian et al., _CMB-S4 Science Book, First Edition_ , 1610.02743.
  * [10] A. A. Starobinsky, _A new type of isotropic cosmological models without singularity_ , _Physics Letters B_ 91 (Mar., 1980) 99–102.
  * [11] N. Afshordi and J. Magueijo, _The Critical Geometry of a Thermal Big Bang_ , 1603.03312.
  * [12] Planck collaboration, Y. Akrami et al., _Planck 2018 results. X. Constraints on inflation_ , _Astron. Astrophys._ 641 (2020) A10, [1807.06211].
  * [13] G. Gubitosi, M. Lagos, J. Magueijo and R. Allison, _Bayesian evidence and predictivity of the inflationary paradigm_ , _JCAP_ 06 (2016) 002, [1506.09143].
  * [14] F. Bezrukov, A. Magnin, M. Shaposhnikov and S. Sibiryakov, _Higgs inflation: consistency and generalisations_ , _Journal of High Energy Physics_ 2011 (Jan., 2011) 16, [1008.5157].
  * [15] J. Magueijo, _Bimetric varying speed of light theories and primordial fluctuations_ , _Phys. Rev._ D79 (2009) 043525, [0807.1689].
  * [16] J. Magueijo, J. Noller and F. Piazza, _Bimetric structure formation: non-Gaussian predictions_ , _Phys. Rev._ D82 (2010) 043521, [1006.3216].
  * [17] Planck collaboration, Y. Akrami et al., _Planck 2018 results. IX. Constraints on primordial non-Gaussianity_ , _Astron. Astrophys._ 641 (2020) A9, [1905.05697].
  * [18] J. Khoury and F. Piazza, _Rapidly-Varying Speed of Sound, Scale Invariance and Non-Gaussian Signatures_ , _JCAP_ 0907 (2009) 026, [0811.3633].
  * [19] X. Chen, M.-x. Huang, S. Kachru and G. Shiu, _Observational signatures and non-Gaussianities of general single field inflation_ , _JCAP_ 0701 (2007) 002, [hep-th/0605045].
  * [20] J. Noller and J. Magueijo, _Non-Gaussianity in single field models without slow-roll_ , _Phys. Rev._ D83 (2011) 103511, [1102.0275].
  * [21] J. M. Maldacena, _Non-Gaussian features of primordial fluctuations in single field inflationary models_ , _JHEP_ 05 (2003) 013, [astro-ph/0210603].
  * [22] D. Seery and J. E. Lidsey, _Primordial non-Gaussianities in single field inflation_ , _JCAP_ 0506 (2005) 003, [astro-ph/0503692].
  * [23] M. Abramowitz and I. A. Stegun, _Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables_.  Dover, New York, ninth dover printing, tenth gpo printing ed., 1964.
  * [24] P. A. R. Ade, N. Aghanim, M. Arnaud, F. Arroja, M. Ashdown, J. Aumont et al., _Planck2015 results_ , _Astronomy & Astrophysics_ 594 (Sep, 2016) A17.
  * [25] X. Chen, M.-x. Huang, S. Kachru and G. Shiu, _Observational signatures and non-gaussianities of general single-field inflation_ , _Journal of Cosmology and Astroparticle Physics_ 2007 (Jan, 2007) 002–002.
  * [26] R. Holman and A. J. Tolley, _Enhanced non-gaussianity from excited initial states_ , _Journal of Cosmology and Astroparticle Physics_ 2008 (May, 2008) 001.
  * [27] M. Liguori, E. Sefusatti, J. R. Fergusson and E. P. S. Shellard, _Primordial non-gaussianity and bispectrum measurements in the cosmic microwave background and large-scale structure_ , _Advances in Astronomy_ 2010 (2010) 1–64.
  * [28] A. Ashoorioon, G. Geshnizjani and H. J. Kim, _Non-Gaussianities in the Extended EFT of Inflation_ , _arXiv e-prints_ (Oct., 2021) arXiv:2110.06038, [2110.06038].
  * [29] J. Magueijo, _New varying speed of light theories_ , _Rept. Prog. Phys._ 66 (2003) 2025, [astro-ph/0305457].
  * [30] J. W. Moffat, _Superluminary universe: A Possible solution to the initial value problem in cosmology_ , _Int. J. Mod. Phys. D_ 2 (1993) 351–366, [gr-qc/9211020].
  * [31] A. Albrecht and J. Magueijo, _A Time varying speed of light as a solution to cosmological puzzles_ , _Phys. Rev. D_ 59 (1999) 043516, [astro-ph/9811018].
  * [32] J. D. Barrow and J. Magueijo, _Varying alpha theories and solutions to the cosmological problems_ , _Phys. Lett. B_ 443 (1998) 104–110, [astro-ph/9811072].
  * [33] J. Magueijo, _Covariant and locally Lorentz invariant varying speed of light theories_ , _Phys. Rev. D_ 62 (2000) 103521, [gr-qc/0007036].
  * [34] J. W. Moffat, _Scalar-tensor-vector gravity theory_ , _JCAP_ 03 (2006) 004, [gr-qc/0506021].
  * [35] V. Adamchik and O. Marichev, _The algorithm for calculating integrals of hypergeometric type functions and its realization in reduce system_ , in _Proceedings of the international symposium on Symbolic and algebraic computation_ , pp. 212–224, 1990.