# Opening the reheating box in multifield inflation

Jérôme Martin  Lucas Pinol

(today)

###### Abstract

The robustness of multi-field inflation to the physics of reheating is
investigated. In order to carry out this study, reheating is described in
detail by means of a formalism which tracks the evolution of scalar fields and
perfect fluids in interaction (the inflatons and their decay products). This
framework is then used to establish the general equations of motion of the
background and perturbative quantities controlling the evolution of the system
during reheating. Next, these equations are solved exactly by means of a new
numerical code. Moreover, new analytical techniques, allowing us to interpret
and approximate these solutions, are developed. As an illustration of a
physical prediction that could be affected by the micro-physics of reheating,
the amplitude of non-adiabatic perturbations in double inflation is
considered. It is found that ignoring the fine-structure of reheating, as
usually done in the standard approach, can lead to differences as big as $\sim
50\%$, while our semi-analytic estimates can reduce this error to $\sim 10\%$.
We conclude that, in multi-field inflation, tracking the perturbations through
the details of the reheating process is important and, to achieve good
precision, requires the use of numerical calculations.

##  1 Introduction

One of the most important result recently obtained in the field of Cosmology
is that the physical conditions that prevailed in the very early Universe can
be convincingly described by the theory of cosmic inflation [1, 2, 3, 4, 5, 6,
7, 8, 9, 10]. Moreover, and this is even more impressive, detailed pieces of
information about the mechanism responsible for this inflationary phase have
started to be gathered. So far, all astrophysical data are compatible with
single-field models [11] and the corresponding potential is known to be of the
plateau shape (or close to it) [12, 13, 14], the prototypical scenario
satisfying these physical requirements being the Starobinsky model [1].

Does that mean that the final word has been said about the mechanism that
drives inflation? There are reasons to believe this is not the case. In
particular, since inflation proceeds at very high energies, much higher
energies than the ones probed in accelerators, its description must be based
on extensions of the standard model of particle physics. Generically, these
frameworks predict the presence of several (scalar) fields. So a natural arena
for designing well-motivated models seems to be multi-field inflation (see,
e.g., Refs. [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27] for works
that are particularly relevant for the present article). In this sense, a
fundamental and so far open problem consists in explaining how, from this
multi-field description, an effective single-field scenario emerges, as
revealed by the data.

It is also interesting to notice that the above described situation could,
maybe, change in the near future. Indeed, the next generation of experiments
will probe even further the micro-physics of inflation, for instance by
measuring in refined details the properties of the large scale structures
and/or of the Cosmic Microwave Background (CMB) anisotropies [28, 29, 30].
These experiments will soon come on line. They could reveal that the physical
nature of inflation, so far compatible (and, to some extent, hidden) in an
effective single-field framework, is in fact multi-field at a deeper level.
Given the present efficiency of the single-field description, however, the
corresponding observational signatures, if ever detected, are likely to
manifest themselves as small deviations. The previous considerations indicate
that it is important to know what physical predictions are expected from
multi-field inflation and to calculate those predictions with enough accuracy.

The typical signatures of multi-field inflation, namely those that would allow
us to distinguish multi-field inflation from single-field inflation, have
already been studied in great details in the literature. Among the emblematic
predictions of multi-field inflation that are not compatible with a single-
field description are a violation of consistency relation
$r=-8n_{\scriptscriptstyle{\mathrm{T}}}$ ($r$ is the tensor to scalar ratio
and $n_{\scriptscriptstyle{\mathrm{T}}}$ is the tensor spectral index) [21],
the presence, at a level to be determined, of Non-Gaussianity (NG) [31] and of
non-adiabatic perturbations [32, 33, 34]. In fact, to be more precise, NG is
also compatible with single-field models provided the fields have non-
canonical terms but non-canonical and multi-field NG may be of different type
and therefore can be, at least in principle, distinguished (see, e.g., Refs.
[35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50] for specific
signatures of multi-field inflation in the NG signal). On the other hand, the
presence of non-adiabatic modes in the data would establish, beyond any doubt,
that the dynamics during inflation involves several degrees of freedom.

At this stage, it ought to be mentioned that there is a fundamental difference
between single-field and multi-field inflation: in multi-field inflation, the
details of what happens during the reheating epoch [51, 52, 53, 54, 55, 56,
57, 58, 59] affect the behavior of the perturbations even on super-Hubble
scales [60, 61, 62, 63, 64, 65] (note, however, that even in single-field
inflation, the small scales can be affected by the details of metric
preheating [66, 67, 68, 69]; for interesting references on preheating in
multi-field inflation, see Refs. [70, 71, 72, 73, 74, 75, 76]). Technically,
this is because, in presence of non-adiabatic perturbations, curvature
perturbations can evolve, even on large scales. This implies that calculating
the behavior of the system during inflation is a priori not sufficient to
establish the predictions of a multi-field scenario: we need to model the
reheating and to follow the perturbations during this phase. What we have just
described has important implications for the above discussion. Since the
reheating phase is a priori complicated, so is an accurate calculation of the
corresponding properties of multi-field scenarios. However, as we have already
mentioned, we need such accurate estimates in order to optimize the use of
cosmological data and to constrain inflationary scenarios beyond the single-
field framework.

In the literature, considerable efforts have been made to understanding what
happens during the (multi-field) inflationary phase but much less efforts have
been devoted to investigating the influence of reheating. This will be the
main focus of the present paper. Very often, the dynamics of reheating is
simply ignored and the predictions follow from simple assumptions about the
continuity of the relevant quantities used to describe the inflationary
scenario under scrutiny. In this article, on the contrary, we will model
reheating in detail by carefully following the decay of the inflaton fields
during this epoch. This will require the use of a formalism which explicitly
includes interactions between scalar fields and perfect fluids (since perfects
fluids will be assumed to be a good description of the inflaton fields decay
products) both at the background and perturbative levels [77, 65, 78]. Then,
we will follow the fate of the perturbations through this detailed description
of reheating and study how it depends on the parameters of the models, for
instance the decay rates. This will also allow us to compare approaches where
the details of reheating are ignored to the framework of this article where
the dynamics of reheating is taken into account. As a consequence, we will be
able to assess how much precision on the final predictions of a model is lost
by using a simplified treatment of reheating. As already mentioned, since the
curvature perturbation is not necessarily conserved during reheating, the
corresponding effect is expected to be relevant.

As noticed above, multi-field inflationary scenarios can be quite complicated
and lead to several predictions that differ from those of single-field models.
As a simple illustration of how the predictions in multi-field inflation
depend on reheating, we will therefore not try to be exhaustive nor to
consider only fully realistic scenarios (that is to say, necessarily fully
compatible with the most recent astrophysical data) but, rather, we will
concentrate on one observable, namely the amplitude of non-adiabatic
perturbations. Moreover, when it comes to concrete numbers, we will focus on a
specific scenario, namely double inflation [15, 16, 17, 20]. We will calculate
exactly the evolution of the system by mean of two numerical codes (in order
to check our conclusions). Then, we will compare these exact results with
several analytical approaches, some of them being already present in the
literature [15, 17] and some others being introduced here. The previous
program will allow us to address the main question of the paper, namely
assessing the robustness of the predictions in multi-field inflation to
changes in the physics of reheating.

This article is organized as follows. After this foreword, Sec. 1, we discuss
how exchanges between fluids can be introduced and modeled, at the background
level in Sec. 2, and at the perturbed level in Sec. 3. Then, in Sec. 4, we
consider the previous formalism in the general case of reheating after a model
of two-field inflation and, in Sec. 5, we apply this approach to a concrete
model, namely two-field inflation. In Sec. 6, we compare our results to the
ones already existing in the literature. Finally, we present our conclusions
in Sec. 7. We end the paper with an appendix A where we develop new analytical
techniques to describe the behavior of the background at the time of the heavy
field decay.

##  2 Scalar fields and fluids in presence of energy-momentum exchanges

###  2.1 General equations

In this paper, we consider a situation where there are several fluids living
in a space-time described by the four-dimensional metric tensor
$g_{\mu\nu}(x^{\kappa})$; each fluid is characterized by its own energy-
momentum tensor, $T_{\mu\nu}^{(\alpha)}$. The index “$(\alpha)$” is written
between parenthesis to indicate that it is not a space-time index but a label
identifying the fluid. However, in the following, when no confusion is
possible, we will not write the parenthesis in order to avoid cluttered
notations. Also, the fluid label will indifferently appear either up or down.
Concretely, the different fluids will be either scalar fields or perfect
fluids with constant equations of state (although one could easily accommodate
time-dependent equations of state). A crucial aspect of the physical situation
considered in the present article is that energy-momentum transfers between
the fluids will be possible. This will be described by the following
conservation equations, based on a detailed balance analysis [77, 65, 78]

$\displaystyle\nabla_{\nu}T^{\mu\nu}_{(\alpha)}=\sum_{(\beta)}\left[Q^{\mu}_{(\alpha)\rightarrow(\beta)}-Q^{\mu}_{(\beta)\rightarrow(\alpha)}\right],$
(2.1)

where $\nabla_{\nu}$ denotes the covariant derivative associated with the
metric $g_{\mu\nu}$. The vector $Q^{\mu}_{(\alpha)\rightarrow(\beta)}$
describes the transfer from the fluid “$(\alpha)$” to the fluid “$(\beta)$”
while, evidently, the vector $Q^{\mu}_{(\beta)\rightarrow(\alpha)}$ describes
the transfer from the fluid “$(\beta)$” to the fluid “$(\alpha)$”. The vector
$Q^{\mu}_{(\alpha)\rightarrow(\beta)}$ can always be written in a covariant,
non-perturbative way, as [77, 65, 78]

$\displaystyle
Q^{\mu}_{(\alpha)\rightarrow(\beta)}=Q_{(\alpha)\rightarrow(\beta)}u^{\mu}+f^{\mu}_{(\alpha)\rightarrow(\beta)},$
(2.2)

with $f^{\mu}_{(\alpha)\rightarrow(\beta)}u_{\mu}=0$ where $u^{\mu}$ is the
total velocity of matter. We see that $Q^{\mu}_{(\alpha)\rightarrow(\beta)}$
has been decomposed in terms of a scalar, $Q_{(\alpha)\rightarrow(\beta)}$,
and a vector, $f^{\mu}_{(\alpha)\rightarrow(\beta)}$.

We have just mentioned that individual energy momentum tensors are not
conserved due to the presence of possible transfers between the fluids. But,
clearly, the total energy momentum tensor must be conserved because it sources
Einstein equations

$\displaystyle
G_{\mu\nu}=\frac{1}{M_{\scriptscriptstyle{\mathrm{Pl}}}^{2}}\sum_{(\alpha)}T_{\mu\nu}^{(\alpha)},$
(2.3)

where, of course, $G_{\mu\nu}\equiv R_{\mu\nu}-Rg_{\mu\nu}/2$ is the Einstein
tensor ($R_{\mu\nu}$ is the Ricci tensor and $R$ the scalar curvature) while
$M_{\scriptscriptstyle{\mathrm{Pl}}}$ is the reduced Planck mass. In order to
satisfy the Bianchi identities, one must have
$\nabla_{\nu}[\sum_{(\alpha)}T_{(\alpha)}^{\mu\nu}]=0$. Physically, as already
mentioned, this expresses the fact that, if energy momentum transfers can
occur among fluids, the total energy momentum must be conserved and the net
sum of all transfers must equal zero.

###  2.2 Scalar fields

In this section, we consider the case where the fluids under consideration are
all scalar fields. Therefore, we assume that we have $N_{\mathrm{field}}$
scalars with canonical kinetic terms, that we denote $\phi_{a}$ with
$a=1,\cdots,N_{\mathrm{field}}$ (as it was the case in the last subsection for
indices between parenthesis, $a$ is not a space-time index and, in the
following, for notational convenience, will be displayed either up or down).
As is well-known, the corresponding stress-energy tensor can be written as

$\displaystyle
T_{\mu\nu}=\sum_{a=1}^{N_{\mathrm{field}}}\left(\partial_{\mu}\phi_{a}\partial_{\nu}\phi_{a}-\frac{1}{2}g_{\mu\nu}g^{\alpha\beta}\partial_{\alpha}\phi_{a}\partial_{\beta}\phi_{a}\right)-g_{\mu\nu}V(\phi_{1},\cdots,\phi_{N_{\mathrm{field}}}),$
(2.4)

where $V(\phi_{a})$ is the potential function that we do not need to specify
at this stage. A crucial remark is that, because the potential term
$V(\phi_{a})$ is a priori non-separable, the above stress-energy tensor cannot
be written as the sum of individual stress-energy tensors, namely cannot be
written as $\sum_{a=1}^{N_{\mathrm{field}}}T_{\mu\nu}^{a}$ where
$T_{\mu\nu}^{a}$ would be the stress-energy tensor associated with the field
$\phi_{a}$. However, a collection of $N_{\mathrm{field}}$ scalar fields can
also be viewed as a collection of $N_{\mathrm{field}}+1$ perfect fluids with
constant equations of state [77, 65, 78]. In this approach, we have
$N_{\mathrm{field}}$ “kinetic fluids” with stress-energy tensor111In
accordance with the remark made in Sec. 2.1, we have written the stress-energy
tensors without a parenthesis around the labels “$K_{a}$” and “$V$” since, in
this case, there is no possible confusion with a space-time index.

$\displaystyle
T_{\mu\nu}^{K_{a}}=\partial_{\mu}\phi_{a}\partial_{\nu}\phi_{a}-\frac{1}{2}g_{\mu\nu}g^{\alpha\beta}\partial_{\alpha}\phi_{a}\partial_{\beta}\phi_{a},$
(2.5)

and one “potential fluid” with stress-energy tensor

$\displaystyle
T_{\mu\nu}^{V}=-g_{\mu\nu}V(\phi_{1},\cdots,\phi_{N_{\mathrm{field}}}),$ (2.6)

and the total stress-energy tensor can be expressed as the sum of the stress-
energy tensors of the kinetic and potential fluids,
$T_{\mu\nu}=\sum_{a=1}^{N_{\mathrm{field}}}T_{\mu\nu}^{K_{a}}+T_{\mu\nu}^{V}$.
More precisely, the kinetic fluids have energy density and pressure

$\displaystyle\rho_{K_{a}}=p_{K_{a}}=-\frac{1}{2}g^{\alpha\beta}\partial_{\alpha}\phi_{a}\partial_{\beta}\phi_{a},$
(2.7)

which shows that each of them has a constant, “stiff”, equation of state,
$w_{K_{a}}\equiv p_{K_{a}}/\rho_{K_{a}}=1$. The kinetic fluids have velocity
$u_{\mu}^{K_{a}}$ defined by
$u_{\mu}^{K_{a}}u_{\nu}^{K_{a}}=-\partial_{\mu}\phi_{a}\partial_{\nu}\phi_{a}/(g^{\alpha\beta}\partial_{\alpha}\phi_{a}\partial_{\beta}\phi_{a})$
and, clearly, one verifies that $u_{\mu}^{K_{a}}u^{\mu}_{K_{a}}=-1$ as
expected. On the other hand, the energy density and pressure of the potential
fluid are given by

$\displaystyle\rho_{V}=-p_{V}=V(\phi_{1},\cdots,\phi_{N_{\mathrm{field}}}),$
(2.8)

which means that this fluid has a vacuum equation of state, namely $w_{V}=-1$.
The velocity of the potential fluid is not defined, which is not problematic
since it will be shown that, in fact, this quantity never appears in the
equations of motion and is, therefore, not relevant.

For the previous description to hold, there is however a price to pay: we must
assume that the kinetic and potential fluids interact. This is indeed
necessary in order to recover the correct equations of motion (namely, the
Klein-Gordon equation) of the scalar fields. This can be shown as follows.
Using Eqs. (2.1) and (2.6), the conservation equation for the potential fluid
can be expressed as

$\displaystyle\nabla_{\nu}T_{V}^{\mu\nu}=-g^{\mu\nu}\partial_{\nu}V=-\sum_{a=1}^{N_{\mathrm{field}}}\partial^{\mu}\phi_{a}\frac{\partial
V}{\partial\phi_{a}}=\sum_{a=1}^{N_{\mathrm{field}}}\left(Q^{\mu}_{V\rightarrow
K_{a}}-Q^{\mu}_{K_{a}\rightarrow V}\right),$ (2.9)

which is satisfied by

$\displaystyle Q^{\mu}_{V\rightarrow K_{a}}=0,\quad Q^{\mu}_{K_{a}\rightarrow
V}=\partial^{\mu}\phi_{a}\frac{\partial V}{\partial\phi_{a}}\,,$ (2.10)

where no sum is meant despite the repeated index “$a$” (recall that “$a$” is
not a space-time index). Then, if one writes the conservation equation for the
kinetic fluids, using Eq. (2.5), one arrives at

$\displaystyle\nabla_{\nu}T_{K_{a}}^{\mu\nu}=\partial^{\mu}\phi_{a}\nabla_{\nu}\nabla^{\nu}\phi_{a}=\sum_{b\neq
a}\left(Q^{\mu}_{K_{a}\rightarrow K_{b}}-Q^{\mu}_{K_{b}\rightarrow
K_{a}}\right)+Q^{\mu}_{K_{a}\rightarrow V}-Q^{\mu}_{V\rightarrow K_{a}},$
(2.11)

and, recalling Eqs. (2.10), this reduces to the known Klein-Gordon equation
only if $Q^{\mu}_{K_{a}\rightarrow K_{b}}=0$ for any index $a$ and $b$.

We conclude that a situation with $N_{\mathrm{field}}$ scalar fields is in
fact equivalent to a situation with $N_{\mathrm{field}}+1$ perfect fluids
(with constant equations of state) in interaction [77, 65, 78]. This effective
interaction is such that the interaction among kinetic fluids vanishes and the
energy-momentum transfer only proceeds from the kinetic fluids to the
potential one (and not the opposite).

###  2.3 Scalar fields and perfect fluids in interaction

In this article, since our goal is to study the reheating in multi-field
inflationary models, where, at the end of inflation, the inflaton fields decay
in various channels the physical properties of which can be described by means
of hydro-dynamical considerations, we are interested in a situation where
there are scalar fields and perfect fluids in the Universe. We have just seen
that scalar fields can in fact be viewed as a collection of perfect fluids,
provided those fluids interact in a specific way. Therefore, a situation with
scalar fields and perfect fluids is in fact equivalent to a situation where
there are only perfect fluids, some of them being “fictitious” and some others
being “real”. In the previous section 2.2, we have written the interactions
between the “fictitious” kinetic and potential fluids obtained from the
requirement that the usual equations of motion describing the behavior of
scalar fields are recovered. In this section, we consider the “real”
interaction between a scalar field and a “real” fluid.

Firstly, we notice that, on general grounds, the interaction between two
fluids can be characterized in a covariant, non-perturbative way, by the
following exchange vector [77, 65, 78]

$\displaystyle
Q^{\mu}_{(\alpha)\rightarrow(\beta)}=\Gamma_{(\alpha)\rightarrow(\beta)}T^{\mu\nu}_{(\alpha)}u_{\nu}^{(\alpha)},$
(2.12)

where $\Gamma_{(\alpha)\rightarrow(\beta)}>0$ is a coefficient that controls
the strength of the interaction and which can also be interpreted as a decay
rate. If the fluid has a perfect fluid form (in this article, we do not
consider fluids which have anisotropic stress), namely
$T_{\mu\nu}^{(\alpha)}=\left[\rho_{(\alpha)}+p_{(\alpha)}\right]u_{\mu}^{(\alpha)}u_{\nu}^{(\alpha)}+p_{(\alpha)}g_{\mu\nu}$
then, it is easy to show that

$\displaystyle
Q^{\mu}_{(\alpha)\rightarrow(\beta)}=-\Gamma_{(\alpha)\rightarrow(\beta)}\rho_{(\alpha)}u^{\mu}_{(\alpha)}\,.$
(2.13)

Secondly, we can apply the previous considerations to the questions studied in
this article. The crucial ingredient is to realize that the interaction
between a scalar field and a “real” fluid can in fact be viewed as an
interaction between the corresponding “fictitious” kinetic fluid and the
“real” fluid. In other words, this interaction will be characterized by
$Q^{\mu}_{K_{a}\rightarrow(\beta)}$ given by Eq. (2.12) and
$Q^{\mu}_{(\beta)\rightarrow K_{a}}=0$. The “fictitious” potential fluid will
remain decoupled from all the “real” fluids present in space-time, namely
$Q^{\mu}_{V\rightarrow(\beta)}=Q^{\mu}_{(\beta)\rightarrow V}=0$.

###  2.4 The Homogeneous and Isotropic FLRW Universe

The results discussed in the previous subsections are valid for any metric
tensor, namely for any space-time. We now assume that the Universe is
homogeneous and isotropic on large scales and, therefore, is described by a
Friedmann-Lemaitre-Robertson-Walker (FLRW) metric ${\rm d}s^{2}=-{\rm
d}t^{2}+a^{2}(t)\delta_{ij}{\rm d}x^{i}{\rm d}x^{j}$, where $t$ is the cosmic
time and $a(t)$ is the scale factor. The metric can also be written ${\rm
d}s^{2}=a^{2}(\eta)(-{\rm d}\eta^{2}+\delta_{ij}{\rm d}x^{i}{\rm d}x^{j})$,
where $\eta$ is the conformal time related to cosmic time $t$ by ${\rm
d}t=a(\eta){\rm d}\eta$. For a FLRW Universe (using conformal time), the total
velocity of matter is given by $u^{\mu}=(1/a,{\bm{0}})$,
$u_{\mu}=(-a,{\bm{0}})$. On very general grounds, for any type of fluid, the
relation $f^{\mu}_{(\alpha)\rightarrow(\beta)}u_{\mu}=0$ implies that, in Eq.
(2.2), one has $f^{0}_{(\alpha)\rightarrow(\beta)}=0$ and, since
$f^{i}_{(\alpha)\rightarrow(\beta)}$ must vanish in an homogeneous and
isotropic background222Indeed, this quantity can be decomposed into a scalar
and a transverse 3-vector,
$f^{i}_{(\alpha)\rightarrow(\beta)}=\partial_{i}f_{(\alpha)\rightarrow(\beta)}+\tilde{f}^{i}_{(\alpha)\rightarrow(\beta)}$
with $\partial_{i}\tilde{f}^{i}_{(\alpha)\rightarrow(\beta)}=0$. At the
background level, the scalar part must vanish because the universe is
homogeneous and the transverse 3-vector must also vanish because the universe
is isotropic., we reach the conclusion that
$f^{\mu}_{(\alpha)\rightarrow(\beta)}=0$. Therefore, this quantity simply does
not appear at the background level. As a consequence, for any type of fluid
living in a FLRW Universe, one has (using conformal time)
$Q^{0}_{(\alpha)\rightarrow(\beta)}=Q_{(\alpha)\rightarrow(\beta)}/a$ and
$Q^{i}_{(\alpha)\rightarrow(\beta)}=0$, with, using Eq. (2.2),

$\displaystyle
Q_{(\alpha)\rightarrow(\beta)}=-\Gamma_{(\alpha)\rightarrow(\beta)}\rho_{(\alpha)},$
(2.14)

since at the background level $u^{\mu}_{(\alpha)}=u^{\mu}$ for any $(\alpha)$
and, as we have just seen, $f^{\mu}_{(\alpha)\rightarrow(\beta)}=0$.

####  2.4.1 Scalar fields in the FLRW Universe

Now, let us see how the formalism described in Secs. 2.1, 2.2 and 2.3 can be
applied in an homogeneous and isotropic Universe filled with
$N_{\mathrm{field}}$ scalar fields. The total energy density and pressure can
be written as

$\displaystyle\rho=\sum_{a=1}^{N_{\mathrm{field}}}\frac{\phi_{a}^{\prime}{}^{2}}{2a^{2}}+V(\phi_{1},\cdots,\phi_{N_{\mathrm{field}}}),\quad
p=\sum_{a=1}^{N_{\mathrm{field}}}\frac{\phi_{a}^{\prime}{}^{2}}{2a^{2}}-V(\phi_{1},\cdots,\phi_{N_{\mathrm{field}}}),$
(2.15)

where a prime denotes a derivative with respect to conformal time. As already
discussed, it is not possible to decompose $\rho$ and $p$ as a sum over
$N_{\mathrm{field}}$ individual energy densities and pressures $\rho_{a}$ and
$p_{a}$, one for each scalar field, unless the potential is separable, namely
unless we deal with the particular case where
$V(\phi_{1},\cdots,\phi_{N_{\mathrm{field}}})=\sum_{a=1}^{N_{\mathrm{field}}}V_{a}(\phi_{a})$.
In general, the conservation equation leads to the equations of motion for the
fields $\phi_{a}$, namely

$\displaystyle\sum_{a=1}^{N_{\mathrm{field}}}\frac{\phi_{a}^{\prime}}{a}\left(\phi_{a}^{\prime\prime}+2{\cal
H}\phi_{a}^{\prime}+a^{2}V_{\phi_{a}}\right)=0,$ (2.16)

where the quantity ${\cal H}=a^{\prime}/a$ is the conformal Hubble parameter,
related to the Hubble parameter $H=\dot{a}/a$ (a dot standing for a derivative
with respect to cosmic time) by ${\cal H}=aH$ and where
$V_{\phi_{a}}\equiv\partial V/\partial\phi_{a}$. Eq. (2.16) is satisfied if,
for each field $\phi_{a}$, one has

$\displaystyle\phi_{a}^{\prime\prime}+2{\cal
H}\phi_{a}^{\prime}+a^{2}V_{\phi_{a}}=0,$ (2.17)

which is the standard form of the Klein-Gordon equations in a FLRW Universe.

As discussed in the previous sections 2.2 and 2.3, another way to proceed is
to introduce $N_{\mathrm{field}}$ kinetic fluids and one potential fluid,
which are perfect fluids with energy density and pressure given by

$\displaystyle\rho_{K_{a}}=\frac{\phi_{a}^{\prime}{}^{2}}{2a^{2}}=p_{K_{a}},\quad\rho_{V}=V=-p_{V}.$
(2.18)

Then, one can write the total energy density and pressure as
$\rho=\sum_{a=1}^{N_{\mathrm{field}}}\rho_{K_{a}}+\rho_{V}$ and
$p=\sum_{a=1}^{N_{\mathrm{field}}}p_{K_{a}}+p_{V}$. In order to recover the
Klein-Gordon equations, we must introduce interactions between these
$N_{\mathrm{field}}+1$ fluids which, in an homogeneous and isotropic Universe,
take the following form [we recall that all
$Q^{i}_{(\alpha)\rightarrow(\beta)}=0$ at the background level]

$\displaystyle aQ_{K_{a}\rightarrow V}$
$\displaystyle=-\phi_{a}^{\prime}V_{\phi_{a}},\quad aQ_{V\rightarrow
K_{a}}=0,\quad aQ_{K_{a}\rightarrow K_{b}}=0.$ (2.19)

Indeed, with the above transfer vectors, the most general conservation
equation for the kinetic fluid, namely

$\displaystyle\rho_{K_{a}}^{\prime}+3{\cal
H}\left(\rho_{K_{a}}+p_{K_{a}}\right)$
$\displaystyle=a\left(Q_{K_{a}\rightarrow V}-Q_{V\rightarrow
K_{a}}\right)+a\sum_{b\neq a}\left(Q_{K_{a}\rightarrow
K_{b}}-Q_{K_{b}\rightarrow K_{a}}\right),$ (2.20)

simply reproduces the Klein-Gordon equation (2.17) for $\phi_{a}$. On the
other hand, the conservation equation for the potential fluid

$\displaystyle\rho_{V}^{\prime}=a\sum_{a=1}^{N_{\mathrm{field}}}\left(Q_{V\rightarrow
K_{a}}-Q_{K_{a}\rightarrow
V}\right)=\sum_{a=1}^{N_{\mathrm{field}}}\phi_{a}^{\prime}V_{\phi_{a}},$
(2.21)

is identically satisfied.

####  2.4.2 Scalar fields and perfect fluids in interactions in the FLRW
Universe

We now take into account the interactions between $N_{\mathrm{field}}$ scalar
fields and $N_{\mathrm{fluid}}$ “real” fluids. Since we have these
“fictitious” and “real” fluids in the Universe, the total energy density and
pressure can now be expressed as

$\displaystyle\rho=\sum_{a=1}^{N_{\mathrm{field}}}\rho_{K_{a}}+\rho_{V}+\sum_{(\alpha)=1}^{N_{\mathrm{fluid}}}\rho_{(\alpha)},\quad
p=\sum_{a=1}^{N_{\mathrm{field}}}p_{K_{a}}+p_{V}+\sum_{(\alpha)=1}^{N_{\mathrm{fluid}}}P_{(\alpha)}\,.$
(2.22)

Then, one must rewrite the conservation equation (2.20) for the kinetic fluid
and add the terms describing the exchanges between fields (or “fictitious”
fluids) and “real” fluids. This leads to

$\displaystyle\rho_{K_{a}}^{\prime}+3{\cal
H}\left(\rho_{K_{a}}+p_{K_{a}}\right)$
$\displaystyle=a\left(Q_{K_{a}\rightarrow V}-Q_{V\rightarrow
K_{a}}\right)+a\sum_{b\neq a}\left(Q_{K_{a}\rightarrow
K_{b}}-Q_{K_{b}\rightarrow K_{a}}\right)$
$\displaystyle+a\sum_{(\alpha)=1}^{N_{\mathrm{fluid}}}\left[Q_{K_{a}\rightarrow(\alpha)}-Q_{(\alpha)\rightarrow
K_{a}}\right],$ (2.23)

where $Q_{K_{a}\rightarrow(\alpha)}$ is given by Eq. (2.14),
$Q_{(\alpha)\rightarrow K_{a}}=0$ and other interactions verify Eq. (2.19). As
a consequence, one obtains the following modified Klein-Gordon equation

$\displaystyle\phi_{a}^{\prime\prime}+\left[2{\cal
H}+\frac{a}{2}\sum_{(\alpha)=1}^{N_{\mathrm{fluid}}}\Gamma_{K_{a}\rightarrow(\alpha)}\right]\phi_{a}^{\prime}+a^{2}V_{\phi_{a}}=0\,,$
(2.24)

where it is clear that these new interactions result in extra friction terms
for the scalar fields, parameterized by the decay rates
$\Gamma_{K_{a}\rightarrow(\alpha)}$. On the other hand, the conservation
equation for $\rho_{V}$ is not modified by any new terms since the potential
fluid does not interact with the “real” fluids. As a consequence, it is still
identically satisfied.

If we now consider the “real” fluids, the conservation equation (2.1) leads to
only one non-trivial equation, its time-component. In a homogeneous and
isotropic Universe, it reads

$\displaystyle\rho_{(\alpha)}^{\prime}+3{\cal
H}\left[\rho_{(\alpha)}+p_{(\alpha)}\right]=a\sum_{a=1}^{N_{\mathrm{field}}}\left[Q_{(\alpha)\rightarrow
K_{a}}-Q_{K_{a}\rightarrow(\alpha)}\right]+a\sum_{(\beta)=1}^{N_{\mathrm{fluid}}}\left[Q_{(\alpha)\rightarrow(\beta)}-Q_{(\beta)\rightarrow(\alpha)}\right],$
(2.25)

where we recall that $\rho_{(\alpha)}$ and $p_{(\alpha)}$ are the energy
density and pressure of the “real” fluid $(\alpha)$. Using again Eq. (2.14),
one obtains

$\displaystyle\rho_{(\alpha)}^{\prime}+3{\cal
H}\left[\rho_{(\alpha)}+p_{(\alpha)}\right]-\sum_{a=1}^{N_{\mathrm{field}}}\frac{\Gamma_{K_{a}\rightarrow(\alpha)}}{2a}\phi_{a}^{\prime
2}+a\sum_{(\beta)=1}^{N_{\mathrm{fluid}}}\left[\Gamma_{(\alpha)\rightarrow(\beta)}\rho_{(\alpha)}-\Gamma_{(\beta)\rightarrow(\alpha)}\rho_{(\beta)}\right]=0\,,$
(2.26)

where the decays of the scalar fields result in an enhancement of the energy
densities of the “real” fluids, as should be the case during reheating, while
the self-interactions of the “real” fluids can add extra complexity in the
system with positive and negative contributions. As already mentioned, the
space component of the conservation equation is identically satisfied.

The formalism described above is particularly well suited to describe the
transfers of energy that occur, at the background level, from scalar fields to
cosmological fluids at the end of inflation, namely during the reheating. This
formalism will therefore be very useful to study this epoch of the
inflationary scenario which, we recall, is the main target of this paper.
Notice also that the interactions introduced before need not be turned on by
hand at the end of inflation: they are negligible but present during inflation
and dynamically become relevant only when the decay rates become of the order
of the Hubble parameter. This will be exemplified in the following when we
study the case of double inflation. Finally, as an additional remark, let us
stress that this formalism can also be used for the warm inflation scenario
[79, 80].

We now turn to the description of linear fluctuations around a FLRW background
in the presence of interactions between scalar fields and “real” fluids. This
is indeed crucial in order to establish reliable cosmological predictions in
multi-field models since, in that case, and contrary to single-field
scenarios, the curvature perturbation can evolve during reheating, even on
super-Hubble scales. It is therefore important to track the behavior of the
perturbations when the interactions between the inflaton fields and their
decay products play an important role in the dynamics of the Universe. This is
the goal of the next section.

##  3 Theory of cosmological perturbations in presence of energy-momentum
exchanges

###  3.1 General equations

In this section, we consider a Universe which is no longer homogeneous and
isotropic and which is filled with various fluids that can interact with each
others. We assume that the deviations from homogeneity and isotropy are small
and, therefore, can be treated perturbatively [81]. This leads to the theory
of cosmological perturbations in presence of energy-momentum exchanges [77,
65, 78]. As expected, this only differs from the standard approach to
cosmological perturbations in the fact that the perturbed conservation
equations of the various fluids (“fictitious” ones describing scalar fields,
“real” fluids, etc.) acquire new terms to describe these exchanges, in very
much the same way as described before for the background. As a consequence,
the corresponding perturbed line element is written in the standard way,
namely ${\rm d}s^{2}=a^{2}(\eta)\\{-(1-2\phi){\rm
d}\eta^{2}+2(\partial_{i}B){\rm d}x^{i}{\rm
d}\eta+[(1-2\psi)\delta_{ij}+2\partial_{i}\partial_{j}E]{\rm d}x^{i}{\rm
d}x^{j}\\}$ and is characterized by four functions, $\phi$, $B$, $\psi$ and
$E$ that are time and space dependent. Here, we have taken into account only
scalar perturbations and, in principle, the perturbed line element should also
contain a vector and a tensor parts. In this article, we focus on the scalar
sector only, that is decoupled from vectors and tensors at this order of
perturbation theory. In order to be consistent, the matter sector is also
perturbed and we introduce perturbed scalar fields
$\delta\phi_{a}(\eta,{\bm{x}})$ and/or perturbed energy density
$\delta\rho_{(\alpha)}(\eta,{\bm{x}})$ and pressure, $\delta
p_{(\alpha)}(\eta,{\bm{x}})$, for the fluids. For the velocity, one has
$\delta u^{\mu}_{(\alpha)}=(-\phi/a,v^{i}_{(\alpha)}/a)$ [and $\delta
u_{\mu}^{(\alpha)}=(-a\phi,av_{i}^{(\alpha)}+a\partial_{i}B)$], where the
index of $v_{i}$ is raised by $\delta_{ij}$ at this order of perturbation
theory. One can check that it preserves the normalization of the four-velocity
at the perturbative level. As usual, not all perturbed quantities are
physically meaningful because of the gauge problem. In the following, we will
deal with this issue by making use of the so-called gauge-invariant formalism
for cosmological perturbations [82, 81]. The gravity sector will be described
by the Bardeen potentials, $\Phi=\phi+[a(B-E^{\prime})]^{\prime}/a$ and
$\Psi=\psi-{\mathcal{H}}(B-E^{\prime})$. In the matter sector, as is well-
known [82], there are different ways to define gauge-invariant perturbed
energy densities. Here, we work in terms of
$\delta\rho_{(\alpha)}^{\mathrm{(gi)}}=\delta\rho_{(\alpha)}+\rho_{(\alpha)}^{\prime}(B-E^{\prime})$.
For the pressure, we have $\delta p_{(\alpha)}^{\mathrm{(gi)}}=\delta
p_{(\alpha)}+p^{\prime}_{(\alpha)}(B-E^{\prime})$ and for the velocity,
$v^{(\mathrm{gi})}_{(\alpha)}=v_{(\alpha)}+E^{\prime}$ with
$v_{i}^{(\alpha)}=\partial_{i}v^{(\alpha)}$ (and a similar definition for the
gauge-invariant part) since we consider scalar perturbations only. Finally, if
matter is described by scalar fields, the inhomogeneous field fluctuations can
be described in terms of the quantity
$\delta\phi^{(\mathrm{gi})}_{a}=\delta\phi_{a}+\phi_{a}^{\prime}(B-E^{\prime})$.
It can be checked that all these quantities, $\Phi,\Psi$ and all the
perturbations with a “$(\mathrm{gi})$” symbol are gauge-invariant at linear
order in cosmological perturbation, and coincide with the corresponding
quantities in the longitudinal gauge [82].

Let us now examine the equations of motion controlling the behavior of the
perturbations. We obviously have the perturbed Einstein equations

$\displaystyle\delta
G_{\mu\nu}=\frac{1}{M_{\scriptscriptstyle{\mathrm{Pl}}}^{2}}\sum_{(\alpha)}\delta
T_{\mu\nu}^{(\alpha)},$ (3.1)

and, perturbing the conservation equation (2.1), we also obtain another set of
equations, namely

$\displaystyle\delta\left[\nabla_{\nu}T^{\mu\nu}_{(\alpha)}\right]=\sum_{(\beta)}\left[\delta
Q^{\mu}_{(\alpha)\rightarrow(\beta)}-\delta
Q^{\mu}_{(\beta)\rightarrow(\alpha)}\right].$ (3.2)

In this formula, the perturbed stress-energy tensor can be calculated in the
standard way and expressed in terms of the gauge-invariant quantities
introduced before. For the exchange vector, using the covariant, non-
perturbative Eq. (2.2), one has

$\displaystyle\delta Q^{\mu}_{(\alpha)\rightarrow(\beta)}=\delta
Q_{(\alpha)\rightarrow(\beta)}u^{\mu}+Q_{(\alpha)\rightarrow(\beta)}\delta
u^{\mu}+\delta f^{\mu}_{(\alpha)\rightarrow(\beta)}.$ (3.3)

As it was the case for the background, we must satisfy the constraint that the
vector $f^{\mu}_{(\alpha)\rightarrow(\beta)}$ is perpendicular to the total
velocity of matter $u^{\mu}$. At the perturbative level, this means that
$\delta\left[f^{\mu}_{(\alpha)\rightarrow(\beta)}u_{\mu}\right]=0$. This
implies that

$\displaystyle\delta
f^{\mu}_{(\alpha)\rightarrow(\beta)}u_{\mu}+f^{\mu}_{(\alpha)\rightarrow(\beta)}\delta
u_{\mu}=0.$ (3.4)

The second term is vanishing because we have seen that, at the background
level, $f^{\mu}_{(\alpha)\rightarrow(\beta)}=0$. Moreover, $u_{0}\neq 0$,
$u_{i}=0$ and, therefore, one has $\delta
f^{0}_{(\alpha)\rightarrow(\beta)}=0$. From the above considerations, we
deduce that the time and space components of $\delta
Q^{\mu}_{(\alpha)\rightarrow(\beta)}$ are given by

$\displaystyle\delta Q^{0}_{(\alpha)\rightarrow(\beta)}$
$\displaystyle=\frac{1}{a}\delta
Q_{(\alpha)\rightarrow(\beta)}-\frac{\phi}{a}Q_{(\alpha)\rightarrow(\beta)},$
(3.5) $\displaystyle\delta Q^{i}_{(\alpha)\rightarrow(\beta)}$
$\displaystyle=Q_{(\alpha)\rightarrow(\beta)}\delta u^{i}+\delta
f^{i}_{(\alpha)\rightarrow(\beta)}.$ (3.6)

In the following, since we consider scalar perturbations, we will work in
terms of $v$ and $\delta f_{(\alpha)\rightarrow(\beta)}$ defined by
$v_{i}=\partial_{i}v$ where $\delta u^{i}=v^{i}/a$ and $\delta
f_{i}^{(\alpha)\rightarrow(\beta)}=\partial_{i}\delta
f_{(\alpha)\rightarrow(\beta)}$. Let us emphasize again that, in the above
expressions, $\delta u^{i}$ is the space component of the total velocity and
not the space component of the velocity of some individual fluid. In order to
clarify further this point, let us explain how these quantities are related.
The time-space component of the perturbed Einstein equations can be written as

$\displaystyle\Psi^{\prime}+{\cal
H}\Phi=-\frac{a^{2}}{2M_{\scriptscriptstyle{\mathrm{Pl}}}^{2}}\sum_{(\alpha)}\left[\rho_{(\alpha)}+p_{(\alpha)}\right]v_{(\alpha)}^{\rm(gi)}=-\frac{a^{2}}{2M_{\scriptscriptstyle{\mathrm{Pl}}}^{2}}\left(\rho+p\right)v^{\rm(gi)},$
(3.7)

where $\rho=\sum_{(\alpha)}\rho_{(\alpha)}$ and
$p=\sum_{(\alpha)}p_{(\alpha)}$ are the total energy density and pressure,
respectively. This immediately implies that

$\displaystyle
v^{\rm(gi)}=\frac{1}{\rho+p}\sum_{(\alpha)}\left[\rho_{(\alpha)}+p_{(\alpha)}\right]v_{(\alpha)}^{\rm(gi)},$
(3.8)

can be understood as describing the total velocity.

Finally, we must define gauge-invariant exchange vectors or, rather, define
them in terms of gauge-invariant quantities. For the scalar $\delta
Q_{(\alpha)\rightarrow(\beta)}$, we introduce the gauge-invariant expressions
$\delta Q_{(\alpha)\rightarrow(\beta)}^{(\mathrm{gi})}=\delta
Q_{(\alpha)\rightarrow(\beta)}+Q^{\prime}_{(\alpha)\rightarrow(\beta)}(B-E^{\prime})$.
This expression is similar to the definition of the gauge-invariant scalar
field fluctuations, perturbed energy densities, etc. and this is of course due
to the fact that, in all these cases, we deal with scalar quantities.
Regarding $\delta f_{(\alpha)\rightarrow(\beta)}$, the situation is even
simpler: this quantity is already gauge-invariant since it has no background
counterpart.

###  3.2 Perturbed conservation equations for fluids

We now consider Eq. (3.2) and write it explicitly for a perfect fluid. Let us
first start with the time component. It can be expressed as

$\displaystyle\delta\rho_{(\alpha)}^{\rm(gi)}{}^{\prime}$
$\displaystyle+3{\cal H}\left[\delta\rho_{(\alpha)}^{\rm(gi)}+\delta
p_{(\alpha)}^{\rm(gi)}\right]-3\left[\rho_{(\alpha)}+p_{(\alpha)}\right]\Psi^{\prime}+\left[\rho_{(\alpha)}+p_{(\alpha)}\right]\nabla^{2}\left[v_{(\alpha)}^{\rm(gi)}\right]$
$\displaystyle=\sum_{(\beta)}\biggl{[}a\delta
Q_{(\alpha)\rightarrow(\beta)}^{\rm(gi)}+aQ_{(\alpha)\rightarrow(\beta)}\Phi-a\delta
Q_{(\beta)\rightarrow(\alpha)}^{\rm(gi)}-aQ_{(\beta)\rightarrow(\alpha)}\Phi\biggr{]}.$
(3.9)

It is interesting to notice that, even if the velocity appears in the left
hand side of this equation, the scalar $\delta f_{(\alpha)\rightarrow(\beta)}$
is absent. Contrary to the background case, the space component of the
conservation equation leads, at the perturbative level, to an interesting
equation which controls the evolution of the perturbed velocity. It reads

$\displaystyle\left[\rho_{(\alpha)}+p_{(\alpha)}\right]v_{(\alpha)}^{\rm(gi)}{}^{\prime}$
$\displaystyle+\left[\rho_{(\alpha)}^{\prime}+p_{(\alpha)}^{\prime}\right]v_{(\alpha)}^{\rm(gi)}+4{\cal
H}\left[\rho_{(\alpha)}+p_{(\alpha)}\right]v_{(\alpha)}^{\rm(gi)}+\left[\rho_{(\alpha)}+p_{(\alpha)}\right]\Phi+\delta
p_{(\alpha)}^{\rm(gi)}$
$\displaystyle=\sum_{(\beta)}\biggl{[}aQ_{(\alpha)\rightarrow(\beta)}v^{\rm(gi)}+\delta
f_{(\alpha)\rightarrow(\beta)}-aQ_{(\beta)\rightarrow(\alpha)}v^{\rm(gi)}-\delta
f_{(\beta)\rightarrow(\alpha)}\biggr{]}.$ (3.10)

Eqs. (3.2) and (3.2) are two equations which, when considered together with
the perturbed Einstein equations, form a complete set allowing us to follow
the evolution of cosmological perturbations. Note that, in practice, we will
rather use the re-scaled velocities
$\varsigma_{(\alpha)}\equiv[\rho_{(\alpha)}+p_{(\alpha)}]v_{(\alpha)}^{\rm(gi)}$
for numerical integration.

###  3.3 Perturbed conservation equations for scalar fields

We now consider a collection of $N_{\mathrm{field}}$ scalar fields. As
explained in Sec. 2.2, $N_{\mathrm{field}}$ scalar fields can always be seen
as $N_{\mathrm{field}}+1$ perfect fluids provided those ones interact in a
specific way. Of course, this “technical trick” remains true at the
perturbative level and one of the goals of this section is to determine the
form of the corresponding perturbed exchange vectors. As it was the case for
the background, this is achieved by requiring the equations of motion to
reduce to the equations of motion for perturbed scalar fields, namely the
perturbed Klein-Gordon equation. The great advantage of working with
$N_{\mathrm{field}}+1$ fluids (as opposed to $N_{\mathrm{field}}$ scalar
fields) is that, as it was the case for the background before, this allows us
to consistently implement the interaction between fields and “real” fluids at
the perturbative level.

The perturbed kinetic fluids energy density, pressure and velocity are given
by the following expressions:

$\displaystyle\delta\rho_{K_{a}}^{\rm(gi)}$ $\displaystyle=\delta
p_{K_{a}}^{\rm(gi)}=\frac{\phi_{a}^{\prime}}{a^{2}}\delta\phi_{a}^{\rm(gi)}{}^{\prime}-\frac{\phi_{a}^{\prime}{}^{2}}{a^{2}}\Phi,\quad
v_{K_{a}}^{\rm(gi)}=-\frac{\delta\phi_{a}^{\rm(gi)}}{\phi_{a}^{\prime}},$
(3.11)

while the same quantities for the potential fluid can be expressed as

$\displaystyle\delta\rho_{V}^{\rm(gi)}$ $\displaystyle=-\delta
p_{V}^{\rm(gi)}=\sum_{a=1}^{N_{\mathrm{field}}}V_{\phi_{a}}\delta\phi_{a}^{(\mathrm{gi})}.$
(3.12)

Notice that, given the form of its individual energy-momentum tensor, it is
impossible to define a velocity $v_{V}^{\rm(gi)}$ for the potential fluid.
However, as already mentioned after Eq. (2.8), since this quantity never
appears in the equations of motion, this is actually not an issue. We also
notice that the sound speed, defined by $c^{2}_{{}_{\mathrm{S}}}\equiv
p^{\prime}/\rho^{\prime}$, of the kinetic and potential fluids is equal to the
equation of state parameter, $c_{{}_{\mathrm{S}}}^{2}=w$, and (anticipating a
little bit over the following considerations) that those fluids have no
intrinsic entropy perturbations.

Our next move is to consider Eq. (3.2). We take the expression of the energy
density, pressure and velocity for the kinetic fluids given by Eqs. (3.11) and
insert them in Eq. (3.2). This leads to the standard perturbed Klein-Gordon
equation

$\displaystyle\delta\phi_{a}^{\rm(gi)}{}^{\prime\prime}+2{\cal
H}\delta\phi_{a}^{\rm(gi)}{}^{\prime}-\nabla^{2}\left[\delta\phi_{a}^{\rm(gi)}\right]+a^{2}\sum_{b=1}^{N_{\mathrm{field}}}V_{\phi_{a}\phi_{b}}\delta\phi_{b}^{\rm(gi)}{}=4\phi_{a}^{\prime}\Psi^{\prime}-2a^{2}V_{\phi_{a}}\Psi,$
(3.13)

provided the perturbed exchanges take the following form

$\displaystyle a\delta Q_{K_{a}\rightarrow V}$
$\displaystyle=-V_{\phi_{a}}\delta\phi_{a}^{\rm(gi)}{}^{\prime}+V_{\phi_{a}}\phi_{a}^{\prime}\Psi-\sum_{b=1}^{N_{\mathrm{field}}}V_{\phi_{a}\phi_{b}}\phi_{a}^{\prime}\delta\phi_{b}^{\rm(gi)},$
(3.14) $\displaystyle a\delta Q_{V\rightarrow K_{a}}$ $\displaystyle=0,\quad
a\delta Q_{K_{a}\rightarrow K_{b}}=0.$ (3.15)

If one does the same for the potential fluid, namely insert the perturbed
energy density and perturbed pressure given by Eq. (3.12) in Eq. (3.2) with
the exchanges just established in Eqs. (3.14), (3.15), then one verifies that
the corresponding equation is identically satisfied. The choices (3.14) and
(3.15) are therefore consistent. The next step is to proceed in a similar
fashion for the coefficients $\delta f_{(\alpha)\rightarrow(\beta)}$. Clearly,
this has to be done by considering the momentum conservation equation.
Therefore, we insert the expressions (3.11) in Eq. (3.2) and find that it is
automatically satisfied provided one takes

$\displaystyle\delta f_{K_{a}\rightarrow V}$
$\displaystyle=V_{\phi_{a}}\delta\phi_{a}^{\rm(gi)}-\phi_{a}^{\prime}V_{\phi_{a}}\frac{1}{\sum_{b=1}^{N_{\mathrm{field}}}\phi_{b}^{\prime}{}^{2}}\sum_{c=1}^{N_{\mathrm{field}}}\phi_{c}^{\prime}\delta\phi_{c}^{\rm(gi)},\quad\delta
f_{V\rightarrow K_{a}}=\delta f_{K_{a}\rightarrow K_{b}}=0.$ (3.16)

The expression for the only non-vanishing coefficient, $\delta
f_{K_{a}\rightarrow V}$, is quite complicated but can be simplified if one
notices that, in the case under consideration in this section where one has
only scalar fields (and no “real” fluids), the total velocity reads

$\displaystyle v^{\rm(gi)}$
$\displaystyle=\sum_{(\alpha)}\frac{\rho_{(\alpha)}+p_{(\alpha)}}{\rho+p}v_{(\alpha)}^{\rm(gi)}=\frac{1}{\sum_{b}\phi^{\prime}_{b}{}^{2}/a^{2}}\sum_{a=1}^{N_{\mathrm{field}}}\frac{\phi_{a}^{\prime}{}^{2}}{a^{2}}\left[-\frac{\delta\phi_{a}^{\rm(gi)}}{\phi_{a}^{\prime}}\right].$
(3.17)

We see that this reproduces exactly the expression appearing in Eq. (3.16). As
a consequence, it follows that

$\displaystyle\delta f_{K_{a}\rightarrow V}$
$\displaystyle=-\phi_{a}^{\prime}V_{\phi_{a}}v_{K_{a}}^{\rm(gi)}+\phi_{a}^{\prime}V_{\phi_{a}}v^{\rm(gi)}=-aQ_{K_{a}\rightarrow
V}\left[v^{\rm(gi)}-v_{K_{a}}^{\rm(gi)}\right]$ (3.18)

It is important to keep in mind that, here, $v^{\rm(gi)}$ is the total
velocity in the case where they are only scalar fields. In a situation where
they are scalar fields and fluids, we have the same formula (3.18) but the
expression of $v^{\rm(gi)}$ is modified since the velocities of the fluids
participate to the total velocity. Therefore, it is fair to acknowledge that
another nonequivalent possibility would be to define the momentum exchange
with Eq. (3.16). Moreover, let us notice that the exchanges between the
fictitious kinetic and potential fluids cannot be written in the form used in
Eq. (2.12).

###  3.4 Perturbed conservation equations in presence of scalar fields and
fluids in interaction

We now consider the situation which is probably the most relevant for the
present article, namely the case where there are scalar fields and “real”
fluids in interaction. As already mentioned, this interaction is described
covariantly and non-perturbatively by Eq. (2.13), modeling the exchange
between the “fictitious” fluids representing the scalar fields and the “real”
fluids, as well as the “real” fluids self-interactions. We have already shown
that this implies
$Q_{(\alpha)\rightarrow(\beta)}=-\Gamma_{(\alpha)\rightarrow(\beta)}\rho_{(\alpha)}$
for the background, see Eq. (2.14). Then, at the perturbed level, one has

$\displaystyle\delta
Q_{(\alpha)\rightarrow(\beta)}^{\mu}=-\Gamma_{(\alpha)\rightarrow(\beta)}\delta\rho_{(\alpha)}u^{\mu}-\Gamma_{(\alpha)\rightarrow(\beta)}\rho_{(\alpha)}\delta
u^{\mu}.$ (3.19)

This implies that the time component of the perturbed exchange vector can be
expressed as

$\displaystyle\delta
Q^{0}_{(\alpha)\rightarrow(\beta)}=-\frac{1}{a}\Gamma_{(\alpha)\rightarrow(\beta)}\delta\rho_{(\alpha)}+\Gamma_{(\alpha)\rightarrow(\beta)}\rho_{(\alpha)}\frac{\phi}{a},$
(3.20)

and, comparing to Eq. (3.5), this leads to

$\displaystyle\delta
Q_{(\alpha)\rightarrow(\beta)}^{(\mathrm{gi})}=-\Gamma_{(\alpha)\rightarrow(\beta)}\delta\rho_{(\alpha)}^{(\mathrm{gi})}.$
(3.21)

Then, the scalar $\delta f_{(\alpha)\rightarrow(\beta)}$ which is associated
with the spatial component of the exchange vector remains to be determined.
This last one can be written as

$\displaystyle\delta Q^{i}_{(\alpha)\rightarrow(\beta)}$
$\displaystyle=-\Gamma_{(\alpha)\rightarrow(\beta)}\rho_{(\alpha)}\delta
u^{i}_{(\alpha)}=-\Gamma_{(\alpha)\rightarrow(\beta)}\rho_{(\alpha)}\delta
u^{i}+\Gamma_{(\alpha)\rightarrow(\beta)}\rho_{(\alpha)}\left[\delta
u^{i}-\delta u^{i}_{(\alpha)}\right].$

Comparing to Eq. (3.6), this implies that $\delta
f^{i}_{(\alpha)\rightarrow(\beta)}=\Gamma_{(\alpha)\rightarrow(\beta)}\rho_{(\alpha)}\left[\delta
u^{i}-\delta u^{i}_{(\alpha)}\right]$, that is to say

$\displaystyle\delta
f_{(\alpha)\rightarrow(\beta)}=a\Gamma_{(\alpha)\rightarrow(\beta)}\rho_{(\alpha)}\left[v^{\rm(gi)}-v_{(\alpha)}^{\rm(gi)}\right].$
(3.22)

This completes our derivation of the exchange vector describing the
interaction between scalar field and “real” fluids.

###  3.5 Conserved quantities for scalar fields and fluids

####  3.5.1 General case

In this section, we turn to the definition of the so-called “conserved
quantities” (a name which might not be totally appropriate in the present case
since these quantities will not always be “conserved”, or “constant”, in the
presence of several scalar fields and fluids). These quantities play an
important role in the theory of cosmological perturbations for two reasons.
Firstly, their behavior is, at least in some regimes, quite simple which is of
great help to understand the behavior of the perturbations and to test and
check numerical calculations. Secondly, they appear in the definition of non-
adiabatic or entropy perturbations which correspond to typical signatures of
multi-fluid systems.

We introduce two gauge-invariant “conserved” quantities for a given fluid
“$(\alpha)$”, $\zeta_{(\alpha)}$ and ${\cal R}_{(\alpha)}$, which are related
to the individual energy density and velocity perturbations, respectively.
Their definition reads [82, 81, 77, 65, 78]

$\displaystyle\zeta_{(\alpha)}=-\Psi-{\cal
H}\frac{\delta\rho_{(\alpha)}^{(\mathrm{gi})}}{\rho^{\prime}_{(\alpha)}},\quad{\cal
R}_{(\alpha)}=\Psi-{\cal H}v_{(\alpha)}^{(\mathrm{gi})}.$ (3.23)

One can also introduce the corresponding “total” quantities (as opposed to
individual) by means of the following expressions, $\zeta=-\Psi-{\cal
H}\delta\rho^{(\mathrm{gi})}/\rho^{\prime}$ and ${\cal R}=\Psi-{\cal
H}v^{(\mathrm{gi})}$, where we recall that
$\delta\rho^{(\mathrm{gi})}=\sum_{(\alpha)}\delta\rho_{(\alpha)}^{(\mathrm{gi})}$,
$\rho^{\prime}=\sum_{(\alpha)}\rho_{(\alpha)}^{\prime}$ and
$v^{(\mathrm{gi})}$ is the total velocity defined in Eq. (3.8). It is then
easy to show that

$\displaystyle\zeta=\sum_{(\alpha)}\frac{\rho^{\prime}_{(\alpha)}}{\rho^{\prime}}\zeta_{(\alpha)},\quad{\cal
R}=\sum_{(\alpha)}\frac{\rho_{(\alpha)}+p_{(\alpha)}}{\rho+p}{\cal
R}_{(\alpha)},$ (3.24)

where $\rho$ and $p$ are the total energy density and pressure, respectively.
The quantities $\zeta$ and ${\cal R}$ are nothing but the weighted sum of the
individual $\zeta_{(\alpha)}$ and ${\cal R}_{(\alpha)}$. We notice, however,
that the weight is not the same for $\zeta$ and ${\cal R}$ unless the
individual fluids are separately conserved in which case
$\rho^{\prime}_{(\alpha)}/\rho^{\prime}=[\rho_{(\alpha)}+p_{(\alpha)}]/(\rho+p)$.
The quantities $\zeta$ and ${\cal R}$ correspond to the well-known curvature
perturbations on constant energy density slices
($\zeta=-\psi|_{\delta\rho=0}$) and on co-moving slices
($\mathcal{R}=\psi|_{v+B=0}$), respectively.

Our next move is to derive the equation of motion for the individual conserved
quantities. Using the equations of motion established before, straightforward
but lengthy calculations lead to

$\displaystyle\zeta_{(\alpha)}^{\prime}$ $\displaystyle=-\frac{1}{3{\cal
H}}\nabla^{2}\left[\Psi-{\cal R}_{(\alpha)}\right]+\frac{3{\cal
H}^{2}}{\rho^{\prime}_{(\alpha)}}\delta
p_{\mathrm{nad}}^{(\alpha)}-\frac{a}{3{\cal
H}\rho_{(\alpha)}^{\prime}}\sum_{(\beta)}\left[Q_{(\alpha)\rightarrow(\beta)}-Q_{(\beta)\rightarrow(\alpha)}\right]\nabla^{2}{\cal
R}_{(\alpha)}$ $\displaystyle-\frac{a}{\rho_{(\alpha)}^{\prime}}\left({\cal
H}-\frac{{\cal H}^{\prime}}{\cal
H}\right)\sum_{(\beta)}\left[Q_{(\alpha)\rightarrow(\beta)}-Q_{(\beta)\rightarrow(\alpha)}\right]\left[\zeta_{(\alpha)}-\zeta\right]$
$\displaystyle-\frac{a{\cal
H}}{\rho_{(\alpha)}^{\prime}}\sum_{(\beta)}\left\\{\delta
Q_{(\alpha)\rightarrow(\beta)}^{(\mathrm{gi})}-\delta
Q_{(\beta)\rightarrow(\alpha)}^{(\mathrm{gi})}-\frac{\delta\rho_{(\alpha)}}{\rho^{\prime}_{(\alpha)}}\left[Q^{\prime}_{(\alpha)\rightarrow(\beta)}-Q^{\prime}_{(\beta)\rightarrow(\alpha)}\right]\right\\}.$
(3.25)

Several comments are in order at this stage. Firstly, the quantity $\delta
p_{\mathrm{nad}}^{(\alpha)}$ is defined by $\delta
p_{\mathrm{nad}}^{(\alpha)}=\delta
p_{(\alpha)}-c_{(\alpha)}^{2}\delta\rho_{(\alpha)}$, where $c_{(\alpha)}$ is
the sound speed for the fluid “$(\alpha)$” [the sound speed has already been
defined in the text, after Eq. (3.12)]. Physically, it represents intrinsic
entropy (or non adiabatic -nad-) perturbations. For a perfect fluid with
constant equation of state, such as the “real” and “fictitious” fluids
considered before, this quantity is vanishing. However, this is not the case
for a scalar field as can be checked by direct inspection. Secondly, the term
between curled brackets in the last line of the above equation can be seen as
a kind of “intrinsic exchange perturbations”. For interactions characterized
by Eqs. (2.12), (2.13) and (2.14), it vanishes. For more complicated
interactions, such as the ones we have introduced between the “fictitious”
fluids, it is not necessarily zero and can contribute to the evolution of
$\zeta_{(\alpha)}$. Thirdly, in absence of exchanges between the fluids and in
absence of intrinsic entropy perturbations, we see that $\zeta_{(\alpha)}$ is
a conserved quantity on large scales. However, if interactions between the
fluids are present, $\zeta_{(\alpha)}$ can evolve even on large scales, thanks
to the term in the second line. Fourthly, it is also possible to establish the
equation satisfied by the “total” $\zeta$. Indeed, $\zeta$ corresponds to a
fluid which is conserved (since its energy momentum tensor has to satisfy the
Einstein equations, see above). As a consequence, all terms related to
exchange vectors have to disappear. Moreover, for the same reason,
$\rho^{\prime}$ can be replaced with $-3{\cal H}(\rho+p)$. It follows that
$\zeta$ obeys the equation

$\displaystyle\zeta^{\prime}$ $\displaystyle=-\frac{1}{3{\cal
H}}\nabla^{2}\left(\Psi-{\cal R}\right)-\frac{{\cal H}}{\rho+p}\delta
p_{\mathrm{nad}},$ (3.26)

with $\delta p_{\mathrm{nad}}=\delta p-c_{{}_{\mathrm{S}}}^{2}\delta\rho$ and
$c_{{}_{\mathrm{S}}}^{2}=p^{\prime}/\rho^{\prime}$, these definitions
involving the total physical quantities (namely total, background and
perturbed, energy density and pressure). We see that $\zeta$ can evolve on
large scales in presence of non-adiabatic pressure.

Let us now present the equation of motion of the individual quantity ${\cal
R}_{(\alpha)}$. Straightforward but quite lengthy manipulations lead to

$\displaystyle{\cal R}_{(\alpha)}^{\prime}$ $\displaystyle=\left({\cal
H}-\frac{{\mathcal{H}}^{\prime}}{{\mathcal{H}}}\right)\left[{\cal R}-{\cal
R}_{(\alpha)}\right]-c_{(\alpha)}^{2}\frac{\rho_{(\alpha)}^{\prime}}{\rho_{(\alpha)}+p_{(\alpha)}}\left[{\cal
R}_{(\alpha)}+\zeta_{\alpha}\right]$
$\displaystyle+\frac{a}{\rho_{(\alpha)}+p_{(\alpha)}}\sum_{(\beta)}\left[Q_{(\alpha)\rightarrow(\beta)}-Q_{(\beta)\rightarrow(\alpha)}\right]\left[{\cal
R}-{\cal R}_{(\alpha)}\right]+\frac{{\cal
H}}{\rho_{(\alpha)}+p_{(\alpha)}}\delta p_{\mathrm{nad}}^{(\alpha)}$
$\displaystyle-\frac{{\cal
H}}{\rho_{(\alpha)}+p_{(\alpha)}}\sum_{(\beta)}\left[\delta
f_{(\alpha)\rightarrow(\beta)}-\delta f_{(\beta)\rightarrow(\alpha)}\right].$
(3.27)

We remark that, contrary to the case of $\zeta_{(\alpha)}$, the scalar $\delta
f_{(\alpha)\rightarrow(\beta)}$ now participates to the equation of motion of
${\cal R}_{(\alpha)}$. The same reasoning as the one used for $\zeta$ allows
us to derive the equation of motion of ${\cal R}$: all terms depending on the
exchange vectors should not be present in that equation since this is an
equation for a fluid (the total fluid) which is conserved and $\rho^{\prime}$
can be replaced by $-3{\cal H}(\rho+p)$. One obtains that

$\displaystyle{\cal R}^{\prime}$ $\displaystyle=3{\cal
H}c_{{}_{\mathrm{S}}}^{2}\left({\cal R}+\zeta\right)+\frac{{\cal
H}}{\rho+p}\delta p_{\mathrm{nad}}.$ (3.28)

As it was already the case before, we notice that the presence of non-
adiabatic pressure can cause the evolution of ${\cal R}$ on large scales.
However, in the present case, one may also wonder about the impact of the
first term in Eq. (3.28). At first sight, it could also be responsible for an
evolution of ${\cal R}$ on large scales. However, one can show that this term
can in fact be expressed in terms of the gradient of the Bardeen potential,
see for instance Eq. (4.40). Therefore, in fact, the presence of non-adiabatic
pressure is the only cause for a potential evolution of ${\cal R}$ on large
scales, as it was the case for $\zeta$.

####  3.5.2 Definitions of non-adiabatic perturbations

We conclude this section 3.5 by introducing the definitions of non-adiabatic
perturbations. These definitions are directly related to the individual
conserved quantities discussed before. The energy density entropy
perturbations $S_{(\alpha)(\beta)}$ and the velocity entropy perturbations
$V_{(\alpha)(\beta)}$ are respectively defined by

$\displaystyle
S_{(\alpha)(\beta)}=3\left[\zeta_{(\alpha)}-\zeta_{(\beta)}\right],\quad
V_{(\alpha)(\beta)}={\cal R}_{(\alpha)}-{\cal R}_{(\beta)}.$ (3.29)

Using Eq. (3.23), one can also write $S_{(\alpha)(\beta)}=-3{\cal
H}[\delta\rho_{(\alpha)}^{(\mathrm{gi})}/\rho_{(\alpha)}^{\prime}-\delta\rho_{(\beta)}^{(\mathrm{gi})}/\rho_{(\beta)}^{\prime}]$
which, if each fluid is separately conserved, reduces to

$\displaystyle
S_{(\alpha)(\beta)}=\frac{\delta\rho_{(\alpha)}^{(\mathrm{gi})}}{\rho_{(\alpha)}+p_{(\alpha)}}-\frac{\delta\rho_{(\beta)}^{(\mathrm{gi})}}{\rho_{(\beta)}+p_{(\beta)}}.$
(3.30)

This matches the standard expression for energy density entropy perturbations.
In the same way, velocity entropy perturbations can be re-expressed as

$\displaystyle V_{(\alpha)(\beta)}=-{\cal
H}\left[v_{(\alpha)}^{(\mathrm{gi})}-v_{(\beta)}^{(\mathrm{gi})}\right],$
(3.31)

and, as expected (and as the name indicates), is related to the difference in
velocities of the two fluids.

##  4 Reheating after multi-field inflation

Having introduced the formalism describing a physical situation where scalar
fields and fluids are present and interacting, we are now in a position where
we can make use of this formalism to study reheating after multi-field
inflation. The equations derived in the previous sections are completely
general and can, in principle, be applied to a situation with an arbitrary
number of fields and fluids in interaction. In this section, however, we
restrict ourselves to two-field inflationary scenarios where each field can
decay into two fluids with constant but otherwise arbitrary equation of state.
This does not restrict the generalities of the results established in this
article since, as already mentioned, from the previous equations, it would be
straightforward to apply the following considerations to a scenario with more
fields and/or fluids. The main advantage of this assumption is that this will
allow to write and work with concrete equations of motion and to discuss
several physical questions in an especially convenient manner. Moreover, this
simple framework covers in fact a very large landscape of models.

###  4.1 Equations of motion

In the following, the two fields will be denoted $\phi_{\ell}$ and
$\phi_{\mathrm{h}}$ and the two perfect fluids will be named “fluid 1” and
“fluid 2”. $\phi_{\ell}$ and $\phi_{\mathrm{h}}$ is the notation used for the
double inflation model where one field is said to be “light” and the other
“heavy”. We will consider this model in great details in what follows, see
Sec. 5, but, at this stage, the potential of the model
$V(\phi_{\ell},\phi_{\mathrm{h}})$ remains arbitrary and the indices “$\ell$”
and “h” must only be viewed as a convenient way to distinguish the two fields
without carrying the meaning it will acquire when, later, we come explicitly
to double inflation. The model we have just introduced is described by the
following Lagrangian

$\displaystyle{\cal L}_{\mathrm{twofields-
inf}}=-\frac{1}{2}\partial^{\mu}\phi_{\ell}\partial_{\mu}\phi_{\ell}-\frac{1}{2}\partial^{\mu}\phi_{\mathrm{h}}\partial_{\mu}\phi_{\mathrm{h}}-V(\phi_{\ell},\phi_{\mathrm{h}})+{\cal
L}_{\mathrm{matter}}+{\cal L}_{\mathrm{int}}.$ (4.1)

As explained before, the two fields responsible for inflation also interact
with other components of matter that, phenomenologically, we represent by
perfect fluids. These fluids are described by ${\cal L}_{\mathrm{matter}}$ and
the interaction between matter and the inflaton fields is supposed to be given
by the term ${\cal L}_{\mathrm{int}}$. Although evidently ever-present, it
becomes relevant only at the end of inflation and will account for the
disintegration of the two inflaton fields explaining the “graceful exit”, i.e.
how the Universe smoothly evolves from inflation to the standard Big Bang
model. Following the formalism used in this article, we will not describe the
interaction between the scalar fields and the fluids by specifying ${\cal
L}_{\mathrm{int}}$ but we will rather proceed as reviewed earlier, at the
level of the non-conservation of the individual energy-momentum tensors.
Indeed, according to the previous considerations, the two scalar fields are in
fact equivalent to three fluids, two kinetic ones and one potential one, with
energy densities and pressures given by

$\displaystyle\rho_{K_{\ell}}$
$\displaystyle=\frac{\phi_{\ell}^{2}{}^{\prime}}{2a^{2}}=p_{K_{\ell}},\quad\rho_{K_{\mathrm{h}}}=\frac{\phi_{\mathrm{h}}^{2}{}^{\prime}}{2a^{2}}=p_{K_{\mathrm{h}}},\quad\rho_{V}=V(\phi_{\ell},\phi_{\mathrm{h}})=-p_{V}.$
(4.2)

We know from Eqs. (3.14) and (3.15) that the exchanges between those three
fluids are given by

$\displaystyle aQ_{K_{\ell}\rightarrow V}$
$\displaystyle=-\phi_{\ell}^{\prime}V_{\phi_{\ell}},\quad aQ_{V\rightarrow
K_{\ell}}=0,\quad aQ_{K_{\mathrm{h}}\rightarrow
V}=-\phi_{\mathrm{h}}^{\prime}V_{\phi_{\mathrm{h}}},\quad aQ_{V\rightarrow
K_{\mathrm{h}}}=0,$ (4.3)

in order to recover the usual equations of motion for the fields. Then, as
already discussed at length before, the crucial ingredient is that the
interactions between scalar fields and fluids are obtained by coupling “real”
cosmological fluids to the “fictitious” kinetic ones only (and not to the
potential fluid). In practice, using Eq. (2.14), one takes

$\displaystyle Q_{K_{\ell}\rightarrow(1)}$ $\displaystyle=-\Gamma_{\ell
1}\rho_{K_{\ell}},\quad Q_{(1)\rightarrow K_{\ell}}=0,\quad
Q_{K_{\ell}\rightarrow(2)}=-\Gamma_{\ell 2}\rho_{K_{\ell}},\quad
Q_{(2)\rightarrow K_{\ell}}=0,$ $\displaystyle
Q_{K_{\mathrm{h}}\rightarrow(1)}$
$\displaystyle=-\Gamma_{\mathrm{h}1}\rho_{K_{\mathrm{h}}},\quad
Q_{(1)\rightarrow K_{\mathrm{h}}}=0,\quad
Q_{K_{\mathrm{h}}\rightarrow(2)}=-\Gamma_{\mathrm{h}2}\rho_{K_{\mathrm{h}}},\quad
Q_{(2)\rightarrow K_{\mathrm{h}}}=0,$ (4.4)

Finally, we make also the hypothesis that the decay products (the “real”
fluids) can interact among themselves. This will be described by

$\displaystyle Q_{(1)\rightarrow(2)}$
$\displaystyle=-\Gamma_{12}\rho_{1},\quad
Q_{(2)\rightarrow(1)}=-\Gamma_{21}\rho_{2}.$ (4.5)

To summarize, the parameters of the model will be the parameters appearing in
the potential $V(\phi_{\ell},\phi_{\mathrm{h}})$ that, as already mentioned,
we do not specify, the equation of state of the two perfect “real” fluids,
$w_{1}$ and $w_{2}$, the parameters describing the decay of the scalar fields
into fluids one and two, $\Gamma_{\ell 1}$, $\Gamma_{\ell 2}$,
$\Gamma_{\mathrm{h}1}$, $\Gamma_{\mathrm{h}2}$ and the possible interaction
between the decay products, $\Gamma_{12}$, $\Gamma_{21}$.

####  4.1.1 Background equations of motion for two-field inflation

We now describe the equations that control the evolution of this system at the
background and perturbative levels. At the background level, we have the
Friedman equation

$\displaystyle
H^{2}=\frac{1}{3M_{\scriptscriptstyle{\mathrm{Pl}}}^{2}}\left(\rho_{K_{\ell}}+\rho_{K_{\mathrm{h}}}+\rho_{V}+\rho_{1}+\rho_{2}\right)=\frac{1}{3M_{\scriptscriptstyle{\mathrm{Pl}}}^{2}}\left[\frac{1}{2}\dot{\phi}^{2}_{\ell}+\frac{1}{2}\dot{\phi}^{2}_{\mathrm{h}}+V\left(\phi_{\mathrm{h}},\phi_{\ell}\right)+\rho_{1}+\rho_{2}\right],$
(4.6)

where $\rho_{1}$ and $\rho_{2}$ are the energy densities associated with the
fluids one and two, respectively. We also have the two Klein-Gordon equations
for the $\phi_{\ell}$ and $\phi_{\mathrm{h}}$ fields, namely

$\displaystyle\ddot{\phi}_{\ell}+\left[3H+\frac{1}{2}\left(\Gamma_{\ell
1}+\Gamma_{\ell 2}\right)\right]\dot{\phi}_{\ell}+\frac{\partial
V}{\partial{\phi_{\ell}}}$
$\displaystyle=0,\quad\ddot{\phi}_{\mathrm{h}}+\left[3H+\frac{1}{2}\left(\Gamma_{\mathrm{h}1}+\Gamma_{\mathrm{h}2}\right)\right]\dot{\phi}_{\mathrm{h}}+\frac{\partial
V}{\partial{\phi_{\mathrm{h}}}}=0.$ (4.7)

Finally, we have the two conservation equations for fluid one and fluid two
which can be written as

$\displaystyle\dot{\rho}_{1}+3H(1+w_{1})\rho_{1}$
$\displaystyle=\frac{1}{2}\Gamma_{\ell
1}\dot{\phi}_{\ell}^{2}+\frac{1}{2}\Gamma_{\mathrm{h}1}\dot{\phi}_{\mathrm{h}}^{2}-\Gamma_{12}\rho_{1}+\Gamma_{21}\rho_{2},$
(4.8) $\displaystyle\dot{\rho}_{2}+3H(1+w_{2})\rho_{2}$
$\displaystyle=\frac{1}{2}\Gamma_{\ell
2}\dot{\phi}_{\ell}^{2}+\frac{1}{2}\Gamma_{\mathrm{h}2}\dot{\phi}_{\mathrm{h}}^{2}+\Gamma_{12}\rho_{1}-\Gamma_{21}\rho_{2}.$
(4.9)

The previous formulas form a closed set of equations which, when solved,
provides a complete solution for the background behavior.

####  4.1.2 Perturbed equations of motion for two-field inflation

We now turn our attention to the equations of motion for the perturbations. We
have seen that, at the perturbative level, the exchanges between the
inhomogeneous fluids are described by the two scalars $\delta
Q_{(\alpha)\rightarrow(\beta)}^{(\mathrm{gi})}$ and $\delta
f_{(\alpha)\rightarrow(\beta)}$. In the case considered here, upon using Eqs.
(3.21), these coefficients can be expressed as

$\displaystyle\delta Q_{K_{\ell}\rightarrow(1)}^{(\mathrm{gi})}$
$\displaystyle=-\Gamma_{\ell
1}\delta\rho_{K_{\ell}}^{(\mathrm{gi})},\quad\delta
Q_{K_{\ell}\rightarrow(2)}^{(\mathrm{gi})}=-\Gamma_{\ell
2}\delta\rho_{K_{\ell}}^{(\mathrm{gi})},$ $\displaystyle\delta
Q_{K_{\mathrm{h}}\rightarrow(1)}^{(\mathrm{gi})}$
$\displaystyle=-\Gamma_{\mathrm{h}1}\delta\rho_{K_{\mathrm{h}}}^{(\mathrm{gi})},\quad\delta
Q_{K_{\mathrm{h}}\rightarrow(2)}^{(\mathrm{gi})}=-\Gamma_{\mathrm{h}2}\delta\rho_{K_{\mathrm{h}}}^{(\mathrm{gi})},$
$\displaystyle\delta Q_{(1)\rightarrow(2)}^{(\mathrm{gi})}$
$\displaystyle=-\Gamma_{12}\delta\rho_{1}^{(\mathrm{gi})},\quad\delta
Q_{(2)\rightarrow(1)}^{(\mathrm{gi})}=-\Gamma_{21}\delta\rho_{2}^{(\mathrm{gi})},$
(4.10)

all the other quantities $\delta
Q_{(\alpha)\rightarrow(\beta)}^{(\mathrm{gi})}$ being zero. Regarding the
coefficients related to momentum transfers, they are given by Eq. (3.22) and
read

$\displaystyle\delta f_{K_{\ell}\rightarrow(1)}$ $\displaystyle=a\Gamma_{\ell
1}\frac{\phi_{\ell}^{\prime}{}^{2}}{2a^{2}}\left[v^{(\mathrm{gi})}-v_{K_{\ell}}^{(\mathrm{gi})}\right],\quad\delta
f_{K_{\ell}\rightarrow(2)}=a\Gamma_{\ell
2}\frac{\phi_{\ell}^{\prime}{}^{2}}{2a^{2}}\left[v^{(\mathrm{gi})}-v_{K_{\ell}}^{(\mathrm{gi})}\right],$
(4.11) $\displaystyle\delta f_{K_{\mathrm{h}}\rightarrow(1)}$
$\displaystyle=a\Gamma_{\mathrm{h}1}\frac{\phi_{\mathrm{h}}^{\prime}{}^{2}}{2a^{2}}\left[v^{(\mathrm{gi})}-v_{K_{\mathrm{h}}}^{(\mathrm{gi})}\right],\quad\delta
f_{K_{\mathrm{h}}\rightarrow(2)}=a\Gamma_{\mathrm{h}2}\frac{\phi_{\mathrm{h}}^{\prime}{}^{2}}{2a^{2}}\left[v^{(\mathrm{gi})}-v_{K_{\mathrm{h}}}^{(\mathrm{gi})}\right],$
(4.12) $\displaystyle\delta f_{(1)\rightarrow(2)}$
$\displaystyle=a\Gamma_{12}\rho_{1}\left[v^{(\mathrm{gi})}-v_{1}^{(\mathrm{gi})}\right],\quad\delta
f_{(2)\rightarrow(1)}=a\Gamma_{21}\rho_{2}\left[v^{(\mathrm{gi})}-v_{2}^{(\mathrm{gi})}\right],$
(4.13)

the other coefficients $\delta f_{(\alpha)\rightarrow(\beta)}$ vanishing. In
the above expressions, $v^{(\mathrm{gi})}$ is the total velocity defined by
Eq. (3.8) which, in the present context, reads

$\displaystyle v^{(\mathrm{gi})}$
$\displaystyle=\frac{1}{\phi_{\mathrm{h}}^{\prime}{}^{2}/a^{2}+\phi_{\ell}^{\prime}{}^{2}/a^{2}+(1+w_{1})\rho_{1}+(1+w_{2})\rho_{2}}\biggl{[}-\frac{\phi_{\mathrm{h}}^{\prime}{}^{2}}{a^{2}}\frac{\delta\phi_{\mathrm{h}}^{(\mathrm{gi})}}{\phi_{\mathrm{h}}^{\prime}}-\frac{\phi_{\ell}^{\prime}{}^{2}}{a^{2}}\frac{\delta\phi_{\ell}^{(\mathrm{gi})}}{\phi_{\ell}^{\prime}}$
$\displaystyle+(1+w_{1})\rho_{1}v_{1}^{(\mathrm{gi})}+(1+w_{2})\rho_{2}v_{2}^{(\mathrm{gi})}\biggr{]}.$
(4.14)

We see that quantities related to fluids and scalar fields all participate to
the expression of the total velocity.

We are now in a position where the equations of motion can be derived. From
the perturbed Einstein equations (the time-space component), we get

$\displaystyle\frac{\dot{\Phi}}{H}$
$\displaystyle=-\Phi-\frac{1}{2M_{\scriptscriptstyle{\mathrm{Pl}}}^{2}H}\left[-\dot{\phi}_{\ell}\delta\phi_{\ell}^{(\mathrm{gi})}-\dot{\phi}_{\mathrm{h}}\delta\phi_{\mathrm{h}}^{(\mathrm{gi})}+a\rho_{1}(1+w_{1})v_{1}^{(\mathrm{gi})}+a\rho_{2}(1+w_{2})v_{2}^{(\mathrm{gi})}\right].$
(4.15)

Notice that, in the present context, the two Bardeen potentials are equal:
$\Phi=\Psi$. From the time component of the conservation equation (3.2) for
the perturbed “fictitious” kinetic fluids, we obtain the perturbed Klein-
Gordon equations for the two fields in presence of energy-momentum exchanges.
They read

$\displaystyle\ddot{\delta\phi}_{\ell}^{(\mathrm{gi})}$
$\displaystyle+\left[3H+\frac{1}{2}\left(\Gamma_{\ell 1}+\Gamma_{\ell
2}\right)\right]\dot{\delta\phi}_{\ell}^{(\mathrm{gi})}+\frac{\partial^{2}V}{\partial\phi_{\ell}^{2}}\delta\phi_{\ell}^{(\mathrm{gi})}+\frac{\partial^{2}V}{\partial\phi_{\ell}\partial\phi_{\mathrm{h}}}\delta\phi_{\mathrm{h}}^{(\mathrm{gi})}+\frac{k^{2}}{a^{2}}\delta\phi_{\ell}^{(\mathrm{gi})}+2\frac{\partial
V}{\partial\phi_{\ell}}\Phi$ $\displaystyle+\frac{1}{2}\left(\Gamma_{\ell
1}+\Gamma_{\ell
2}\right)\dot{\phi}_{\ell}\Phi-4\dot{\phi}_{\ell}\dot{\Phi}=0,$ (4.16)
$\displaystyle\ddot{\delta\phi}_{\mathrm{h}}^{(\mathrm{gi})}$
$\displaystyle+\left[3H+\frac{1}{2}\left(\Gamma_{\mathrm{h}1}+\Gamma_{\mathrm{h}2}\right)\right]\dot{\delta\phi}_{\mathrm{h}}^{(\mathrm{gi})}+\frac{\partial^{2}V}{\partial\phi_{\mathrm{h}}^{2}}\delta\phi_{\mathrm{h}}^{(\mathrm{gi})}+\frac{\partial^{2}V}{\partial\phi_{\mathrm{h}}\partial\phi_{\ell}}\delta\phi_{\ell}^{(\mathrm{gi})}+\frac{k^{2}}{a^{2}}\delta\phi_{\mathrm{h}}^{(\mathrm{gi})}+2\frac{\partial
V}{\partial\phi_{\mathrm{h}}}\Phi$
$\displaystyle+\frac{1}{2}\left(\Gamma_{\mathrm{h}1}+\Gamma_{\mathrm{h}2}\right)\dot{\phi}_{\mathrm{h}}\Phi-4\dot{\phi}_{\mathrm{h}}\Phi=0.$
(4.17)

In a similar way, Eq. (3.2), written for the fluid one and two lead to the
equations controlling the evolution of the perturbed energy density for those
fluids

$\displaystyle\dot{\delta\rho_{1}}^{(\mathrm{gi})}$
$\displaystyle+3H(1+w_{1})\delta\rho_{1}^{(\mathrm{gi})}-3(1+w_{1})\rho_{1}\dot{\Phi}-a(1+w_{1})\rho_{1}\frac{k^{2}}{a^{2}}v_{1}^{(\mathrm{gi})}-\Gamma_{\ell
1}\left[\delta\rho_{K_{\ell}}^{(\mathrm{gi})}+\rho_{K_{\ell}}\Phi\right]$
$\displaystyle-\Gamma_{\mathrm{h}1}\left[\delta\rho_{K_{\mathrm{h}}}^{(\mathrm{gi})}+\rho_{K_{\mathrm{h}}}\Phi\right]+\Gamma_{12}\left[\delta\rho_{1}^{(\mathrm{gi})}+\rho_{1}\Phi\right]-\Gamma_{21}\left[\delta\rho_{2}^{(\mathrm{gi})}+\rho_{2}\Phi\right]=0,$
(4.18) $\displaystyle\dot{\delta\rho_{2}}^{(\mathrm{gi})}$
$\displaystyle+3H(1+w_{2})\delta\rho_{2}^{(\mathrm{gi})}-3(1+w_{2})\rho_{2}\dot{\Phi}-a(1+w_{2})\rho_{2}\frac{k^{2}}{a^{2}}v_{2}^{(\mathrm{gi})}-\Gamma_{\ell
2}\left[\delta\rho_{K_{\ell}}^{(\mathrm{gi})}+\rho_{K_{\ell}}\Phi\right]$
$\displaystyle-\Gamma_{\mathrm{h}2}\left[\delta\rho_{K_{\mathrm{h}}}^{(\mathrm{gi})}+\rho_{K_{\mathrm{h}}}\Phi\right]+\Gamma_{21}\left[\delta\rho_{2}^{(\mathrm{gi})}+\rho_{2}\Phi\right]-\Gamma_{12}\left[\delta\rho_{1}^{(\mathrm{gi})}+\rho_{1}\Phi\right]=0.$
(4.19)

Notice that the quantity $\delta\rho_{K_{a}}^{(\mathrm{gi})}+\rho_{K_{a}}\Phi$
that often appears in the above formulas can also be re-written as
$\delta\rho_{K_{a}}^{(\mathrm{gi})}+\rho_{K_{a}}\Phi=\dot{\phi}_{a}\dot{\delta\phi_{a}}^{(\mathrm{gi})}-\dot{\phi}_{a}^{2}\Phi/2$.
We also remark that the time component of the conservation equation does not
explicitly depend on the coefficients $\delta f_{(\alpha)\rightarrow(\beta)}$.
As a consequence, if one only tracks perturbed energy densities (and/or energy
density non-adiabatic perturbations), one could be under the impression that
they can be ignored and that the above equations are sufficient. However, this
is not the case because the velocities affect the evolution of the energy
densities and the velocities equation of motion, see for instance Eq. (3.2),
do depend on the coefficients $\delta f_{(\alpha)\rightarrow(\beta)}$.
Therefore, even if at the end one is only interested in the density
perturbations, the momentum exchanges must be specified.

Finally, we now turn to the space component of the conservation equations, see
Eq. (3.2). Let us first write the space component of the conservation equation
for the “fictitious” kinetic fluid associated to, say, $\phi_{\ell}$ (the
conclusion is the same for $\phi_{\mathrm{h}}$). Using Eq. (3.11), this gives

$\displaystyle-\frac{\phi_{\ell}^{\prime\prime}}{a^{2}}\delta\phi_{\ell}^{(\mathrm{gi})}-\frac{2}{a^{2}}{\cal
H}\phi_{\ell}^{\prime}\delta\phi_{\ell}^{(\mathrm{gi})}$
$\displaystyle=a\left[Q_{K_{\ell}\rightarrow
V}+Q_{K_{\ell}\rightarrow(1)}+Q_{K_{\ell}\rightarrow(2)}\right]v^{(\mathrm{gi})}+\delta
f_{K_{\ell}\rightarrow V}$ $\displaystyle+\delta
f_{K_{\ell}\rightarrow(1)}+\delta f_{K_{\ell}\rightarrow(2)}.$ (4.20)

However, the time component of the background equation of motion gives
$\phi_{K_{\ell}}^{\prime}\left(\phi_{K_{\ell}}^{\prime\prime}+2{\cal
H}\phi_{K_{\ell}}^{\prime}\right)/a^{2}=a\left[Q_{K_{\ell}\rightarrow
V}+Q_{K_{\ell}\rightarrow(1)}+Q_{K_{\ell}\rightarrow(2)}\right]$. Plugging
this expression in the space component of the conservation equation leads to

$\displaystyle\delta f_{K_{\ell}\rightarrow V}+\delta
f_{K_{\ell}\rightarrow(1)}+\delta
f_{K_{\ell}\rightarrow(2)}=a\left[Q_{K_{\ell}\rightarrow
V}+Q_{K_{\ell}\rightarrow(1)}+Q_{K_{\ell}\rightarrow(2)}\right]\left[v_{K_{\ell}}^{(\mathrm{gi})}-v^{(\mathrm{gi})}\right],$
(4.21)

which is automatically satisfied, thanks to Eq. (4.11). This result makes
sense because, for scalar fields, there is no other equation of motion than
the Klein-Gordon equation, which is of second order in time and which we have
already derived before.

Let us now study the space component of the conservation equation for the
“real” fluids. Making use of Eq. (3.2), one arrives at

$\displaystyle\dot{v}_{1}^{(\mathrm{gi})}$
$\displaystyle+\left(1-3w_{1}\right)Hv_{1}^{(\mathrm{gi})}+\frac{\Phi}{a}+\frac{w_{1}}{1+w_{1}}\frac{1}{a}\frac{\delta\rho_{1}^{(\mathrm{gi})}}{\rho_{1}}=-\frac{w_{1}}{(1+w_{1})\rho_{1}}Q_{(1)\rightarrow(2)}v_{1}^{(\mathrm{gi})}$
$\displaystyle+\frac{1}{\rho_{1}}Q_{K_{\ell}\rightarrow(1)}\left[v_{1}^{(\mathrm{gi})}-\frac{v_{K_{\ell}}^{(\mathrm{gi})}}{1+w_{1}}\right]+\frac{1}{\rho_{1}}Q_{K_{\mathrm{h}}\rightarrow(1)}\left[v_{1}^{(\mathrm{gi})}-\frac{v_{K_{\mathrm{h}}}^{(\mathrm{gi})}}{1+w_{1}}\right]$
$\displaystyle+\frac{1}{\rho_{1}}Q_{(2)\rightarrow(1)}\left[v_{1}^{(\mathrm{gi})}-\frac{v_{2}^{(\mathrm{gi})}}{1+w_{1}}\right],$
(4.22) $\displaystyle\dot{v}_{2}^{(\mathrm{gi})}$
$\displaystyle+\left(1-3w_{2}\right)Hv_{2}^{(\mathrm{gi})}+\frac{\Phi}{a}+\frac{w_{2}}{1+w_{2}}\frac{1}{a}\frac{\delta\rho_{2}^{(\mathrm{gi})}}{\rho_{2}}=-\frac{w_{2}}{(1+w_{2})\rho_{2}}Q_{(2)\rightarrow(1)}v_{2}^{(\mathrm{gi})}$
$\displaystyle+\frac{1}{\rho_{2}}Q_{K_{\ell}\rightarrow(2)}\left[v_{2}^{(\mathrm{gi})}-\frac{v_{K_{\ell}}^{(\mathrm{gi})}}{1+w_{2}}\right]+\frac{1}{\rho_{2}}Q_{K_{\mathrm{h}}\rightarrow(2)}\left[v_{2}^{(\mathrm{gi})}-\frac{v_{K_{\mathrm{h}}}^{(\mathrm{gi})}}{1+w_{2}}\right]$
$\displaystyle+\frac{1}{\rho_{2}}Q_{(1)\rightarrow(2)}\left[v_{2}^{(\mathrm{gi})}-\frac{v_{1}^{(\mathrm{gi})}}{1+w_{2}}\right].$
(4.23)

Eqs. (4.15), (4.1.2), (4.1.2), (4.1.2), (4.1.2), (4.1.2) and (4.1.2) represent
a complete set allowing us to track the perturbations during inflation and
afterwards. We stress that these equations are valid for any model of two-
field inflation and, in this sense, are quite general. To our knowledge,
although implicitly present in Ref. [83], this is the first time that they are
written explicitly. Similar equations have been studied in Ref. [84] but for
one scalar field and one fluid only. Maybe the closest related work is Ref.
[23] which investigates the generation of entropy fluctuations after multi-
field inflation (this paper studies models where there is a large number of
fields during inflation). But the derivation of the conservation equations
does not follow from a systematic formalism as we have done and, moreover, the
space component of the conservation equations is not presented because Ref.
[23] does not study how velocity non-adiabatic perturbations are produced
(although, as we have discussed before, it is not possible to ignore them
since they influence the evolution of the perturbed density contrasts). We
notice that, when the comparison is possible, the equations presented here are
consistent with those of Ref. [23].

###  4.2 Quantization of cosmological perturbations in multi-field inflation

According to the theory of inflation, all the inhomogeneities (CMB
anisotropies, large scale structures, …) that we observe today originate from
quantum fluctuations that were amplified by gravitational instability and
stretched to astrophysical distances by the cosmic expansion during inflation
[7, 8, 9] (see also Refs. [85, 86, 87]).

In the context of single-field inflation, this means that the perturbed
metric, $\delta\hat{g}_{\mu\nu}(\eta,{\bf x})$ [in practice, the Bardeen
potential $\hat{\Phi}(\eta,{\bm{x}})]$ and the quantity representing matter,
the scalar field fluctuations $\delta\hat{\phi}(\eta,{\bf x})$, must be
quantum operators. In that case, everything can be reduced to the study of a
single variable, the so-called Mukhanov-variable [7, 88, 81]333Here, we denote
this variable by $q(\eta,{\bm{x}})$ and not by the more traditional symbol
$v(\eta,{\bm{x}})$ in order to avoid a possible confusion with the velocities
of the fluids., $\hat{q}(\eta,{\bm{x}})=z(\eta)\hat{\cal R}(\eta,{\bm{x}})$,
where
$z(\eta)=a(\eta)\mathchoice{{\hbox{$\displaystyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=3.22221pt,depth=-2.57779pt}}}M_{\scriptscriptstyle{\mathrm{Pl}}}$.
Then, this quantity is expanded in terms of creation and annihilation
operators,

$\displaystyle\hat{q}(\eta,{\bm{x}})=\frac{1}{(2\pi)^{3/2}}\int\mathrm{d}{\bm{k}}\left[q_{k}(\eta)\hat{c}_{\bm{k}}+q_{k}^{*}(\eta)\hat{c}_{-{\bm{k}}}^{\dagger}\right]e^{i{\bm{k}}\cdot{\bm{x}}},$
(4.24)

with
$\left[\hat{c}_{\bm{k}},\hat{c}^{\dagger}_{{\bm{k}}^{\prime}}\right]=\delta^{(3)}\left({\bm{k}}-{\bm{k}}^{\prime}\right)$.
We have also introduced the mode function $q_{k}(\eta)$ that controls the
evolution of the operator $\hat{q}(\eta,{\bm{x}})$. Notice that, here, we
denote the mode function with the same symbol as the corresponding quantum
operator. No confusion can arise with the Fourier transform of the quantum
operator (which also carries an index ${\bm{k}}$) since the mode function does
not carry the hat symbol. In this framework, the perturbed Einstein equations
discussed in the previous sections must be viewed as operator equations.
However, upon inserting the expansion (4.24) in these operator equations, one
obtains ordinary differential equations for the mode functions $q_{k}(\eta)$.
Of course, the differential equations for the mode functions have in fact the
same form as the operator differential equations. The two-point correlation
function of the Mukhanov-Sasaki variable leads to the definition of the power
spectrum ${\cal P}_{q}$, namely

$\displaystyle\left\langle\hat{q}^{2}(\eta,{\bm{x}})\right\rangle=\int_{0}^{+\infty}\frac{\mathrm{d}k}{k}{\cal
P}_{q}(k),$ (4.25)

with

$\displaystyle{\cal
P}_{q}=\frac{k^{3}}{2\pi^{2}}\left|q_{k}(\eta)\right|^{2}.$ (4.26)

This procedure can be extended to more degrees of freedom, for instance in a
multi-field inflationary setup, which is the case of interest in this paper.
With $N_{\mathrm{field}}$ scalar fields labeled by $a$, there exists a set of
$N_{\mathrm{field}}$ independent annihilation and creation operators that,
however, need not be aligned with the scalar fields themselves, i.e. there
exists a frame labeled by numbers $A$ in which
$\left[\hat{c}_{\bm{k}}^{A},\hat{c}^{\dagger
B}_{{\bm{k}}^{\prime}}\right]=\delta^{AB}\delta^{(3)}\left({\bm{k}}-{\bm{k}}^{\prime}\right)$.
Then, the multi-field Sasaki-Mukhanov variables $q^{a}(\eta,{\bm{x}})$ are
promoted to quantum operators and we have the expansion [19, 89, 24]:

$\displaystyle\hat{q}^{a}(\eta,{\bm{x}})=\frac{1}{(2\pi)^{3/2}}\int\mathrm{d}{\bm{k}}\sum_{A=1}^{N_{\mathrm{field}}}\left[q_{A,k}^{a}(\eta)\hat{c}_{A,{\bm{k}}}+q_{A,k}^{a*}(\eta)\hat{c}_{A,-{\bm{k}}}^{\dagger}\right]e^{i{\bm{k}}\cdot{\bm{x}}},$
(4.27)

where we now have $N_{\mathrm{field}}{}^{2}$ mode functions
$q^{a}_{A,{\bm{k}}}(\eta)$. The calculation of the two-point correlation
function,

$\displaystyle\left\langle\hat{q}^{a}(\eta,{\bm{x}})\hat{q}^{b}(\eta,{\bm{x}})\right\rangle=\int_{0}^{+\infty}\frac{\mathrm{d}k}{k}{\cal
P}_{q}^{ab}(k),$ (4.28)

leads to a definition of the “generalized” power spectrum ${\cal P}_{q}^{ab}$
which reads

$\displaystyle{\cal
P}_{q}^{ab}=\frac{k^{3}}{2\pi^{2}}\sum_{A=1}^{N_{\mathrm{field}}}q_{A,k}^{a}(\eta)\left[q_{A,k}^{b}(\eta)\right]^{*}.$
(4.29)

Interestingly enough, we see that the two-point correlation function receives
contributions from all independent modes labeled by “$A$”.

In this article, we deal with a situation which is slightly different from
what was described before. The reason is of course that we have, at the same
time, the presence of scalar fields and fluids and their associated
fluctuations. Cosmological fluids can be quantized in a way that is very
similar to what was described above, once their Mukhanov-Sasaki variables with
canonical kinetic terms have been identified, see Refs. [88, 81] for details.
As for a scalar field, in the case of a single fluid, the Mukhanov-Sasaki
variable is defined from the conserved quantity, ${\cal R}$, by
$q(\eta,{\bm{x}})=z(\eta){\cal R}(\eta,{\bm{x}})$, with
$z(\eta)=a(\eta)\mathchoice{{\hbox{$\displaystyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=3.22221pt,depth=-2.57779pt}}}M_{\scriptscriptstyle{\mathrm{Pl}}}/c_{\mathrm{S}}$,
see Ref. [88, 81]. Notice that the quantity $z(\eta)$ differs from the case of
a scalar field since there is an additional factor $c_{\mathrm{S}}$
participating to its expression. In the present context, the situation is more
complicated since we deal with several fluids. One should, therefore,
introduce several Mukhanov-Sasaki variables, one for each fluid. To our
knowledge, this question has not been (yet) studied in a comprehensive manner
in the literature but it is clear that each degree of freedom is related to
the individual ${\cal R}_{(\alpha)}$. We will come back to this problem in the
section where we discuss the initial conditions for the perturbations, see
Sec. 4.3. In any case, the individual ${\cal R}_{(\alpha)}$ become quantum
operators $\hat{\cal R}_{(\alpha)}(\eta,{\bm{x}})$ that can be expanded in
terms of creation and annihilation operators and that possess their own mode
functions.

As explained before, once the system has been quantized, the individual
$\hat{\cal R}_{(\alpha)}$ become quantum operators with their associated mode
functions. Clearly, for consistency, this is also the case for any other
quantity participating to the description of matter since they are all related
to the $\hat{\cal R}_{(\alpha)}$. For instance, in the case of scalar fields,
the field fluctuations $\delta\hat{\phi}^{a}(\eta,{\bm{x}})$ can be expanded
as

$\displaystyle\delta\hat{\phi}^{a}(\eta,{\bm{x}})=\frac{1}{(2\pi)^{3/2}}\int\mathrm{d}{\bm{k}}\sum_{A=1}^{N_{\mathrm{field}}+N_{\mathrm{fluid}}}\left[\delta\phi_{A,k}^{a}(\eta)\hat{c}_{A,{\bm{k}}}+\delta\phi_{A,k}^{a*}(\eta)\hat{c}_{A,-{\bm{k}}}^{\dagger}\right]e^{i{\bm{k}}\cdot{\bm{x}}},$
(4.30)

where $\delta\phi^{a}_{A,k}(\eta)$ are the associated mode functions. In a
similar way, for a hydrodynamical fluid, the perturbed energy density
$\delta\hat{\rho}_{(\alpha)}$ and/or the perturbed velocity
$\hat{v}_{(\alpha)}$, must also be viewed as quantum operators. For instance
the perturbed energy density can be expressed as

$\displaystyle\delta\hat{\rho}_{(\alpha)}(\eta,{\bm{x}})=\frac{1}{(2\pi)^{3/2}}\int\mathrm{d}{\bm{k}}\sum_{A=1}^{N_{\mathrm{field}}+N_{\mathrm{fluid}}}\left[\delta\rho_{A,k}^{(\alpha)}(\eta)\hat{c}_{A,{\bm{k}}}+\delta\rho_{A,k}^{(\alpha)*}(\eta)\hat{c}_{A,-{\bm{k}}}^{\dagger}\right]e^{i{\bm{k}}\cdot{\bm{x}}},$
(4.31)

where we have introduced the mode functions
$\delta\rho_{A,k}^{(\alpha)}(\eta)$ that are obviously different from, say,
the mode functions of the perturbed scalar field. However, it is crucial to
remark that the creation and annihilation operators in Eqs. (4.30) and (4.31)
are really the same quantities. In general, with $N_{\mathrm{field}}$ scalar
fields and $N_{\mathrm{fluid}}$ cosmological fluids, the index $A$ runs from
$1$ to $N_{\mathrm{field}}+N_{\mathrm{fluid}}$. In our case of interest with
$2$ scalar fields and $2$ cosmological fluids, the index $A$ runs on the $4$
independent oscillators.

The derivation of the differential equations obeyed by the mode functions of
the system is straightforward. As already mentioned, Eqs. (4.15), (4.1.2),
(4.1.2), (4.1.2), (4.1.2), (4.1.2) and (4.1.2) must now be viewed as
differential equations for quantum operators. Then, we can insert the
canonical expansions of each operators and the same equations for the mode
functions is obtained, however duplicated for each independent oscillator
“$A$”. Concretely, this gives for the perturbed Einstein equations

$\displaystyle\frac{\dot{\Phi}_{A}}{H}$
$\displaystyle=-\Phi_{A}-\frac{1}{2M_{\scriptscriptstyle{\mathrm{Pl}}}^{2}H}\left[-\dot{\phi}_{\ell}\delta\phi_{\ell,A}^{(\mathrm{gi})}-\dot{\phi}_{\mathrm{h}}\delta\phi_{\mathrm{h},A}^{(\mathrm{gi})}+a\rho_{1}(1+w_{1})v_{1,A}^{(\mathrm{gi})}+a\rho_{2}(1+w_{2})v_{2,A}^{(\mathrm{gi})}\right]\,.$
(4.32)

Since, as already mentioned, $A$ runs on four independent oscillators, the
above equation really means four independent equations. For the mode functions
of the field fluctuations operators, one obtains

$\displaystyle\ddot{\delta\phi}_{\ell,A}^{(\mathrm{gi})}$
$\displaystyle+\left[3H+\frac{1}{2}\left(\Gamma_{\ell 1}+\Gamma_{\ell
2}\right)\right]\dot{\delta\phi}_{\ell,A}^{(\mathrm{gi})}+\frac{\partial^{2}V}{\partial\phi_{\ell}^{2}}\delta\phi_{\ell,A}^{(\mathrm{gi})}+\frac{\partial^{2}V}{\partial\phi_{\ell}\partial\phi_{\mathrm{h}}}\delta\phi_{\mathrm{h},A}^{(\mathrm{gi})}+\frac{k^{2}}{a^{2}}\delta\phi_{\ell,A}^{(\mathrm{gi})}+2\frac{\partial
V}{\partial\phi_{\ell}}\Phi_{A}$ $\displaystyle+\frac{1}{2}\left(\Gamma_{\ell
1}+\Gamma_{\ell
2}\right)\dot{\phi}_{\ell}\Phi_{A}-4\dot{\phi}_{\ell}\dot{\Phi}_{A}=0\,,$
(4.33) $\displaystyle\ddot{\delta\phi}_{\mathrm{h},A}^{(\mathrm{gi})}$
$\displaystyle+\left[3H+\frac{1}{2}\left(\Gamma_{\mathrm{h}1}+\Gamma_{\mathrm{h}2}\right)\right]\dot{\delta\phi}_{\mathrm{h},A}^{(\mathrm{gi})}+\frac{\partial^{2}V}{\partial\phi_{\mathrm{h}}^{2}}\delta\phi_{\mathrm{h},A}^{(\mathrm{gi})}+\frac{\partial^{2}V}{\partial\phi_{\mathrm{h}}\partial\phi_{\ell}}\delta\phi_{\ell,A}^{(\mathrm{gi})}+\frac{k^{2}}{a^{2}}\delta\phi_{\mathrm{h},A}^{(\mathrm{gi})}+2\frac{\partial
V}{\partial\phi_{\mathrm{h}}}\Phi_{A}$
$\displaystyle+\frac{1}{2}\left(\Gamma_{\mathrm{h}1}+\Gamma_{\mathrm{h}2}\right)\dot{\phi}_{\mathrm{h}}\Phi_{A}-4\dot{\phi}_{\mathrm{h}}\Phi_{A}=0\,.$
(4.34)

Finally, it remains the equations of motion controlling the behavior of the
mode functions associated to the “real” fluids present in the system. They
read

$\displaystyle\dot{\delta\rho}_{1,A}^{(\mathrm{gi})}+3H(1+w_{1})\delta\rho_{1,A}^{(\mathrm{gi})}-3(1+w_{1})\rho_{1}\dot{\Phi}_{A}-a(1+w_{1})\rho_{1}\frac{k^{2}}{a^{2}}v_{1,A}^{(\mathrm{gi})}$
$\displaystyle-\Gamma_{\ell
1}\left[\delta\rho_{K_{\ell},A}^{(\mathrm{gi})}+\rho_{K_{\ell}}\Phi_{A}\right]-\Gamma_{\mathrm{h}1}\left[\delta\rho_{K_{\mathrm{h}},A}^{(\mathrm{gi})}+\rho_{K_{\mathrm{h}}}\Phi_{A}\right]+\Gamma_{12}\left[\delta\rho_{1,A}^{(\mathrm{gi})}+\rho_{1}\Phi_{A}\right]$
$\displaystyle-\Gamma_{21}\left[\delta\rho_{2,A}^{(\mathrm{gi})}+\rho_{2}\Phi_{A}\right]=0\,,$
(4.35)
$\displaystyle\dot{\delta\rho}_{2,A}^{(\mathrm{gi})}+3H(1+w_{2})\delta\rho_{2,A}^{(\mathrm{gi})}-3(1+w_{2})\rho_{2}\dot{\Phi}_{A}-a(1+w_{2})\rho_{2}\frac{k^{2}}{a^{2}}v_{2,A}^{(\mathrm{gi})}$
$\displaystyle-\Gamma_{\ell
2}\left[\delta\rho_{K_{\ell},A}^{(\mathrm{gi})}+\rho_{K_{\ell}}\Phi_{A}\right]-\Gamma_{\mathrm{h}2}\left[\delta\rho_{K_{\mathrm{h}},A}^{(\mathrm{gi})}+\rho_{K_{\mathrm{h}}}\Phi_{A}\right]+\Gamma_{21}\left[\delta\rho_{2,A}^{(\mathrm{gi})}+\rho_{2}\Phi_{A}\right]$
$\displaystyle-\Gamma_{12}\left[\delta\rho_{1,A}^{(\mathrm{gi})}+\rho_{1}\Phi_{A}\right]=0\,,$
(4.36) $\displaystyle\dot{v}_{1,A}^{(\mathrm{gi})}$
$\displaystyle+\left(1-3w_{1}\right)Hv_{1,A}^{(\mathrm{gi})}+\frac{\Phi_{A}}{a}+\frac{w_{1}}{1+w_{1}}\frac{1}{a}\frac{\delta\rho_{1,A}^{(\mathrm{gi})}}{\rho_{1}}=-\frac{w_{1}}{(1+w_{1})\rho_{1}}Q_{(1)\rightarrow(2)}v_{1,A}^{(\mathrm{gi})}$
$\displaystyle+\frac{1}{\rho_{1}}Q_{K_{\ell}\rightarrow(1)}\left[v_{1,A}^{(\mathrm{gi})}-\frac{v_{K_{\ell},A}^{(\mathrm{gi})}}{1+w_{1}}\right]+\frac{1}{\rho_{1}}Q_{K_{\mathrm{h}}\rightarrow(1)}\left[v_{1,A}^{(\mathrm{gi})}-\frac{v_{K_{\mathrm{h}},A}^{(\mathrm{gi})}}{1+w_{1}}\right]$
$\displaystyle+\frac{1}{\rho_{1}}Q_{(2)\rightarrow(1)}\left[v_{1,A}^{(\mathrm{gi})}-\frac{v_{2,A}^{(\mathrm{gi})}}{1+w_{1}}\right]\,,$
(4.37) $\displaystyle\dot{v}_{2,A}^{(\mathrm{gi})}$
$\displaystyle+\left(1-3w_{2}\right)Hv_{2,A}^{(\mathrm{gi})}+\frac{\Phi_{A}}{a}+\frac{w_{2}}{1+w_{2}}\frac{1}{a}\frac{\delta\rho_{2,A}^{(\mathrm{gi})}}{\rho_{2}}=-\frac{w_{2}}{(1+w_{2})\rho_{2}}Q_{(2)\rightarrow(1)}v_{2,A}^{(\mathrm{gi})}$
$\displaystyle+\frac{1}{\rho_{2}}Q_{K_{\ell}\rightarrow(2)}\left[v_{2,A}^{(\mathrm{gi})}-\frac{v_{K_{\ell},A}^{(\mathrm{gi})}}{1+w_{2}}\right]+\frac{1}{\rho_{2}}Q_{K_{\mathrm{h}}\rightarrow(2)}\left[v_{2,A}^{(\mathrm{gi})}-\frac{v_{K_{\mathrm{h}},A}^{(\mathrm{gi})}}{1+w_{2}}\right]$
$\displaystyle+\frac{1}{\rho_{2}}Q_{(1)\rightarrow(2)}\left[v_{2,A}^{(\mathrm{gi})}-\frac{v_{1,A}^{(\mathrm{gi})}}{1+w_{2}}\right]\,.$
(4.38)

The above differential equations are the equations to be integrated in order
to follow the evolution of the system. In particular, once the mode functions
are known, the two-point correlation functions of any combination of operators
can be evaluated. However, in order to be able to carry out this task, we need
to know the initial conditions for each quantity. We now turn to this
question.

###  4.3 Initial conditions for the perturbations

####  4.3.1 Warm-up: the case of a single fluid

One of the great advantage of the inflationary theory is its ability to
suggest natural and well-defined initial conditions. These initial conditions
can be introduced in several different ways. Here, as a warm-up, we discuss
one method and show how it is related to the formalism introduced before in
the simple (and standard) case where there is only one scalar field or one
perfect fluid (namely, one degree of freedom).

Let us start with the case of one scalar field. The equations for the
conserved quantities $\zeta$ and ${\cal R}$ (here, obviously, the total
$\zeta$ and ${\cal R}$ are the same as the individual conserved quantities
since, in our example, there is only one degree of freedom) are given by Eq.
(3.26) and Eq. (3.28). In the present context, the non-adiabatic pressure that
appears in those equations is the intrinsic non-adiabatic pressure of a scalar
field. It is given by

$\displaystyle\delta
p_{\mathrm{nad}}=-2M_{\scriptscriptstyle{\mathrm{Pl}}}^{2}(1-c_{{}_{\mathrm{S}}}^{2})\frac{k^{2}}{a^{2}}\Phi_{\bm{k}}.$
(4.39)

We see that it is non-vanishing unless $c_{{}_{\mathrm{S}}}=1$, namely a
scalar field is a fluid with non-vanishing intrinsic non-adiabatic pressure,
but it is proportional to $k^{2}/a^{2}$ and, therefore, becomes irrelevant on
large scales. Then, using the above expression of $\delta p_{\mathrm{nad}}$ in
Eq. (3.28) and the fact that the Bardeen potential can be written as

$\displaystyle\Phi_{\bm{k}}=-\frac{3}{2}\frac{a^{2}(\rho+p)}{k^{2}M_{\scriptscriptstyle{\mathrm{Pl}}}^{2}}(\zeta_{\bm{k}}+{\cal
R}_{\bm{k}}),$ (4.40)

it is easy to establish that the quantity ${\cal R}_{\bm{k}}^{\prime}$ can be
re-expressed as ${\cal R}_{\bm{k}}^{\prime}=3{\cal H}({\cal
R}_{\bm{k}}+\zeta_{\bm{k}})$. Then, deriving this expression once, and using
the $\delta p_{\mathrm{nad}}$-independent expression of
$\zeta_{\bm{k}}^{\prime}+{\cal R}_{\bm{k}}^{\prime}$ obtained by summing up
Eq. (3.26) and Eq. (3.28), one arrives at

$\displaystyle{\cal R}_{\bm{k}}^{\prime\prime}+2\frac{z^{\prime}}{z}{\cal
R}_{\bm{k}}^{\prime}+k^{2}{\cal R}_{\bm{k}}=0,$ (4.41)

where we recall that $z\equiv
aM_{\scriptscriptstyle{\mathrm{Pl}}}\mathchoice{{\hbox{$\displaystyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule height=3.22221pt,depth=-2.57779pt}}}$. Finally, introducing
the Mukhanov-Sasaki variable, see above, $q=z{\cal R}$, we obtain the
following equation of a parametric oscillator

$\displaystyle
q_{\bm{k}}^{\prime\prime}+\left(k^{2}-\frac{z^{\prime\prime}}{z}\right)q_{\bm{k}}=0.$
(4.42)

At the beginning of inflation, the physical wavelengths of Fourier modes of
astrophysical interest today were much smaller than the Hubble radius, meaning
that $k^{2}\gg z^{\prime\prime}/z$. As a consequence, in this regime, the
solution of Eq. (4.42) reads

$\displaystyle
q_{\bm{k}}(\eta)=A_{\bm{k}}e^{ik(\eta-\eta_{\mathrm{ini}})}+B_{\bm{k}}e^{-ik(\eta-\eta_{\mathrm{ini}})},$
(4.43)

where $A_{\bm{k}}$ and $B_{\bm{k}}$ are integration constants. Then,
quantizing the fluctuations and assuming that the initial state is the vacuum,
adiabatic or Wentzel Kramer Brillouin (WKB), Bunch-Davies state, it follows
that $A_{\bm{k}}=(2k)^{-1/2}$ and $B_{\bm{k}}=0$.

The above reasoning can easily be repeated if one now assumes that there is
only one perfect fluid in the Universe. The main difference with the scalar
field calculation presented above is that we now have $\delta
p_{\mathrm{nad}}=0$. Let us restart from Eq. (3.28) and take the derivative of
this equation and, then, replace $\zeta^{\prime}_{\bm{k}}$ using Eq. (3.26).
One arrives at

$\displaystyle{\cal R}_{\bm{k}}^{\prime\prime}=\frac{{\cal H}^{\prime}}{{\cal
H}}{\cal
R}_{\bm{k}}^{\prime}+\frac{(c_{{}_{\mathrm{S}}}^{2})^{\prime}}{c_{{}_{\mathrm{S}}}^{2}}{\cal
R}_{\bm{k}}^{\prime}+3{\cal H}c_{{}_{\mathrm{S}}}^{2}{\cal
R}_{\bm{k}}^{\prime}-c_{{}_{\mathrm{S}}}^{2}k^{2}{\cal
R}_{\bm{k}}+c_{{}_{\mathrm{S}}}^{2}k^{2}\Phi_{\bm{k}}.$ (4.44)

The last term in the above equation can be expressed in terms of
$\zeta_{\bm{k}}+{\cal R}_{\bm{k}}$ using Eq. (4.40) and, thanks to Eq. (3.28),
$\zeta_{\bm{k}}+{\cal R}_{\bm{k}}$ is proportional to ${\cal
R}_{\bm{k}}^{\prime}$. As a consequence, one obtains a closed equation for
${\cal R}_{\bm{k}}$ which reads

$\displaystyle{\cal R}_{\bm{k}}^{\prime\prime}=\left[\frac{{\cal
H}^{\prime}}{\cal
H}+\frac{(c_{{}_{\mathrm{S}}}^{2})^{\prime}}{c_{{}_{\mathrm{S}}}^{2}}+3{\cal
H}c_{{}_{\mathrm{S}}}^{2}-{\cal H}\epsilon_{1}\right]{\cal
R}_{\bm{k}}^{\prime}-c_{{}_{\mathrm{S}}}^{2}k^{2}{\cal
R}_{\bm{k}}=-2\frac{z^{\prime}}{z}{\cal
R}_{\bm{k}}^{\prime}-c_{{}_{\mathrm{S}}}^{2}k^{2}{\cal R}_{\bm{k}},$ (4.45)

where
$z=aM_{\scriptscriptstyle{\mathrm{Pl}}}\mathchoice{{\hbox{$\displaystyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule height=3.22221pt,depth=-2.57779pt}}}/c_{{}_{\mathrm{S}}}$.
Notice that $z$ is not well-defined for a pressure-less fluid. Then, as we did
for the case of a scalar field, we define the Mukhanov-Sasaki variable of the
fluid, $q$, by $q=z{\cal R}$ and it follows that

$\displaystyle
q_{\bm{k}}^{\prime\prime}+\left(c_{{}_{\mathrm{S}}}^{2}k^{2}-\frac{z^{\prime\prime}}{z}\right)q_{\bm{k}}=0.$
(4.46)

This equation is very similar to Eq. (4.43), the only difference being the
appearance of the sound velocity in the gradient term and in the definition of
$z(\eta)$ (which, if $c_{{}_{\mathrm{S}}}$ is constant, cancels out in the
term $z^{\prime\prime}/z$). If the Fourier mode under consideration is such
that its wavelength is much smaller than the Hubble horizon initially, then
the solution of the above equation reads
$q_{\bm{k}}(\eta)=A_{\bm{k}}e^{ic_{{}_{\mathrm{S}}}k(\eta-\eta_{\mathrm{ini}})}+B_{\bm{k}}e^{-ic_{{}_{\mathrm{S}}}k(\eta-\eta_{\mathrm{ini}})}$.
Quantizing the hydrodynamical fluctuations and assuming that they are
initially placed in the vacuum state leads to the following initial
conditions: $A_{\bm{k}}=(2c_{{}_{\mathrm{S}}}k)^{-1/2}$ and $B_{\bm{k}}=0$.

Endowed with the initial conditions for the quantity $q_{\bm{k}}(\eta)$
determined above, one can infer the initial conditions of any other variable.
However, this is not straightforward and, as a preparation to the multi-fluid
case, we explain how it can be done in the single fluid case. Let us first
assume that there is only one scalar field. We have two methods to derive the
initial conditions of the Bardeen potential. Using Eq. (4.40) and the fact
that ${\cal R}_{\bm{k}}+\zeta_{\bm{k}}={\cal R}_{\bm{k}}^{\prime}/(3{\cal
H})$, one can easily establish that (notice that this equation is exact)

$\displaystyle\Phi_{\bm{k}}=-\frac{{\cal
H}\epsilon_{1}}{k^{2}}\left(\frac{q_{\bm{k}}}{z}\right)^{\prime},$ (4.47)

or

$\displaystyle a\Phi_{\bm{k}}=-\frac{1}{2k^{2}}\frac{\cal
H}{M_{\scriptscriptstyle{\mathrm{Pl}}}}\mathchoice{{\hbox{$\displaystyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=3.22221pt,depth=-2.57779pt}}}\left(q_{\bm{k}}^{\prime}-\frac{z^{\prime}}{z}q_{\bm{k}}\right),$
(4.48)

from which one can write

$\displaystyle\Phi_{\bm{k}}=-\frac{i}{(2k)^{3/2}}\frac{\cal
H}{M_{\scriptscriptstyle{\mathrm{Pl}}}}\mathchoice{{\hbox{$\displaystyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=3.22221pt,depth=-2.57779pt}}}\left(1-\frac{z^{\prime}}{ikz}\right)e^{ik(\eta-\eta_{\mathrm{ini}})}\simeq-\frac{i}{(2k)^{3/2}}\frac{\cal
H}{M_{\scriptscriptstyle{\mathrm{Pl}}}}\mathchoice{{\hbox{$\displaystyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=3.22221pt,depth=-2.57779pt}}}e^{ik(\eta-\eta_{\mathrm{ini}})},$ (4.49)

the second term between the parenthesis in the above equation being negligible
in the small-scale limit. Therefore, we have that
$\Phi_{\bm{k}}\sim(2k)^{-3/2}$. Then, since
$a\delta\phi_{\bm{k}}^{(\mathrm{gi})}=q_{\bm{k}}-z\Phi_{\bm{k}}$, one has
$\delta\phi_{\bm{k}}^{(\mathrm{gi})}\sim(2k)^{-1/2}$, the second term
proportional to the Bardeen potential giving a sub-dominant contribution.

A second way to derive those results, which was also the method used in Ref.
[84], consists in the following. The single-field version of the Einstein
equation, see Eq. (4.32), can be written as

$\displaystyle\Phi_{\bm{k}}^{\prime}=-{\cal
H}\Phi_{\bm{k}}+\frac{\phi^{\prime}}{2M_{\scriptscriptstyle{\mathrm{Pl}}}^{2}}\delta\phi_{\bm{k}}^{(\mathrm{gi})}.$
(4.50)

We are interested in the sub-Hubble behavior of $\Phi_{\bm{k}}$ and
$\delta\phi_{\bm{k}}^{(\mathrm{gi})}$. In this regime, given the solution
found above for $q_{\bm{k}}$, one can write $\Phi_{\bm{k}}={\cal
A}_{\Phi}e^{ik(\eta-\eta_{\mathrm{ini}})}$ and
$\delta\phi_{\bm{k}}^{(\mathrm{gi})}={\cal
A}_{\delta\phi}e^{ik(\eta-\eta_{\mathrm{ini}})}$, where ${\cal A}_{\Phi}$ and
${\cal A}_{\delta\phi}$ are slowly varying overall amplitudes. Inserting these
expressions into Eq. (4.50) and using the definition of $q_{\bm{k}}$,
$q_{\bm{k}}=a\delta\phi_{\bm{k}}^{(\mathrm{gi})}+z\Phi_{\bm{k}}$, in order to
express ${\cal A}_{\delta\phi}$ in terms of $q_{\bm{k}}$ and ${\cal
A}_{\Phi}$, one obtains

$\displaystyle a{\cal A}^{\prime}_{\Phi}+ik\left[1+\frac{{\cal
H}(1+\epsilon_{1})}{ik}\right]a{\cal A}_{\Phi}=\frac{1}{2}\frac{\cal
H}{M_{\scriptscriptstyle{\mathrm{Pl}}}}\mathchoice{{\hbox{$\displaystyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule height=3.22221pt,depth=-2.57779pt}}}q_{\bm{k}}.$ (4.51)

In this equation, the derivative ${\cal A}_{\Phi}^{\prime}$ can be neglected
since ${\cal A}_{\Phi}$ is changing very slowly; moreover, given that $k/{\cal
H}\gg 1$ on sub-Hubble scales, one arrives at exactly Eq. (4.49). Therefore,
this also completely determines the initial conditions of $\Phi_{\bm{k}}$ in
terms of those of $q_{\bm{k}}$ and gives results similar to the first method.
As a consequence, as announced before, once the initial conditions for the
Mukhanov-Sasaki variable are known, the initial conditions for any other
relevant quantities can be automatically inferred.

The same reasoning can also be applied to the situation where there is only
one perfect fluid in the Universe. In that case, an exact result is

$\displaystyle\Phi_{\bm{k}}=-\frac{{\cal
H}\epsilon_{1}}{c_{{}_{\mathrm{S}}}^{2}k^{2}}\left(\frac{q_{\bm{k}}}{z}\right)^{\prime},$
(4.52)

or

$\displaystyle
a\Phi_{\bm{k}}=-\frac{c_{{}_{\mathrm{S}}}}{2c_{{}_{\mathrm{S}}}^{2}k^{2}}\frac{\cal
H}{M_{\scriptscriptstyle{\mathrm{Pl}}}}\mathchoice{{\hbox{$\displaystyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=3.22221pt,depth=-2.57779pt}}}\left(q_{\bm{k}}^{\prime}-\frac{z^{\prime}}{z}q_{\bm{k}}\right).$
(4.53)

Then, we introduce the same WKB ansatz already discussed before with one
important difference though. Instead of oscillatory terms $\propto
e^{ik\eta}$, one needs to take into account the fact that, in a fluid,
fluctuations propagate with a speed which not the speed of light. As a
consequence, the relevant quantities characterizing the fluid will be written
as a slow-varying amplitude times
$e^{ic_{{}_{\mathrm{S}}}k(\eta-\eta_{\mathrm{ini}})}$. Then, repeating the
previous considerations leads to the following equation

$\displaystyle\Phi_{\bm{k}}=-\frac{ic_{{}_{\mathrm{S}}}}{(2c_{{}_{\mathrm{S}}}k)^{3/2}}\frac{\cal
H}{M_{\scriptscriptstyle{\mathrm{Pl}}}}\mathchoice{{\hbox{$\displaystyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=3.22221pt,depth=-2.57779pt}}}\left(1-\frac{z^{\prime}}{ic_{{}_{\mathrm{S}}}kz}\right)e^{ic_{{}_{\mathrm{S}}}k(\eta-\eta_{\mathrm{ini}})}\simeq-\frac{ic_{{}_{\mathrm{S}}}}{(2c_{{}_{\mathrm{S}}}k)^{3/2}}\frac{\cal
H}{M_{\scriptscriptstyle{\mathrm{Pl}}}}\mathchoice{{\hbox{$\displaystyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=3.22221pt,depth=-2.57779pt}}}e^{ic_{{}_{\mathrm{S}}}k(\eta-\eta_{\mathrm{ini}})},$
(4.54)

and, therefore, that $\Phi_{\bm{k}}\propto(2c_{{}_{\mathrm{S}}}k)^{-3/2}$.
Notice, as expected, that the above equation is exactly Eq. (4.49) if one
takes $c_{{}_{\mathrm{S}}}=1$.

The second method can also be used to check the validity of the result. In
case of one perfect fluid, Eq. (4.32) reduces to

$\displaystyle\Phi_{\bm{k}}^{\prime}=-{\cal
H}\Phi_{\bm{k}}-\frac{a^{2}}{2M_{\scriptscriptstyle{\mathrm{Pl}}}^{2}}\rho(1+w)v^{(\mathrm{gi})},$
(4.55)

and using similar considerations as the ones presented before in the case of a
single scalar field, we obtain

$\displaystyle a{\cal
A}^{\prime}_{\Phi}+ic_{{}_{\mathrm{S}}}k\left[1+\frac{{\cal
H}(1+\epsilon_{1})}{ic_{{}_{\mathrm{S}}}k}\right]a{\cal
A}_{\Phi}=\frac{c_{{}_{\mathrm{S}}}}{2}\frac{\cal
H}{M_{\scriptscriptstyle{\mathrm{Pl}}}}\mathchoice{{\hbox{$\displaystyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule height=3.22221pt,depth=-2.57779pt}}}q_{\bm{k}}.$ (4.56)

that is to say the same formula as for a scalar field except for the presence
of the sound velocity. In the large scale limit, one obtains exactly Eq.
(4.54) for the Bardeen potential and we conclude that
$\Phi_{\bm{k}}\propto(2c_{{}_{\mathrm{S}}}k)^{-3/2}$.

For the velocity, we can use the definition of $q_{\bm{k}}$, namely ${\cal
H}{\cal A}_{v}={\cal A}_{\Phi}-q_{\bm{k}}/z$ and, therefore, ${\cal
A}_{v}\propto(2c_{{}_{\mathrm{S}}}k)^{-1/2}$. In particular, the contribution
coming from the Bardeen potential is sub-dominant as expected.

Finally, the initial conditions for the density contrast remains to be
discussed. In the single fluid case, the perturbed energy density conservation
equation, see Eq. (4.2), reads

$\delta\rho^{\prime}+3{\cal
H}(1+w)\delta\rho-3(1+w)\rho\Phi^{\prime}-(1+w)\rho k^{2}v=0,$ (4.57)

yielding

$\displaystyle{\cal
A}_{\delta\rho}^{\prime}+ikc_{{}_{\mathrm{S}}}\left[1+\frac{3{\cal
H}(1+w)}{ikc_{{}_{\mathrm{S}}}}\right]{\cal
A}_{\delta\rho}=3(1+w)\rho\left({\cal
A}_{\Phi}^{\prime}+ikc_{{}_{\mathrm{S}}}{\cal A}_{\Phi}\right)+(1+w)\rho
k^{2}{\cal A}_{v}.$ (4.58)

In the right hand side, the dominant term is the term proportional to ${\cal
A}_{v}$ which scales as $\propto k^{3/2}$ since the term proportional to
${\cal A}_{\Phi}$ “only” scales $\propto k^{-1/2}$. As a consequence, one can
write

$\displaystyle{\cal
A}_{\delta\rho}=\frac{(1+w)\rho}{ic_{{}_{\mathrm{S}}}^{2}}c_{{}_{\mathrm{S}}}k{\cal
A}_{v},$ (4.59)

and deduce that ${\cal A}_{\delta\rho}\propto k^{1/2}$. It is interesting to
test the consistency of this result with the other conservation equation
which, in the single fluid case, reads [see Eq. (4.2)]

$\displaystyle v^{\prime}+3{\cal
H}(1-3w)v+\Phi+\frac{w}{1+w}\frac{\delta\rho}{\rho}=0.$ (4.60)

Inserting the WKB ansatz in this equation, one obtains

$\displaystyle{\cal A}_{v}^{\prime}+ic_{{}_{\mathrm{S}}}k\left[1+\frac{3{\cal
H}(1-3w)}{ic_{{}_{\mathrm{S}}}k}\right]{\cal A}_{v}=-{\cal
A}_{\Phi}-\frac{w}{(1+w)\rho}{\cal A}_{\delta\rho}.$ (4.61)

In the right hand side the term proportional to ${\cal A}_{\Phi}$ is
subdominant since it scales $\propto k^{-3/2}$ while the term proportional to
${\cal A}_{\delta\rho}$ is proportional to $k^{-1/2}$. As a consequence, one
obtains

$\displaystyle{\cal
A}_{\delta\rho}=-\frac{(1+w)\rho}{w}ic_{{}_{\mathrm{S}}}k{\cal A}_{v},$ (4.62)

which is exactly Eq. (4.59) since $w=c_{{}_{\mathrm{S}}}^{2}$ for a perfect
fluid with constant equation of state. We conclude that this is entirely
consistent with the results obtained from the energy density perturbation
conservation equation and that the relevant initial conditions have now been
completely specified.

####  4.3.2 The many-fluid case

In the previous section, we have shown how to connect the formalism presented
in Sec. 3 to the approach utilizing the Mukhanov-Sasaki variable. In
principle, the generalization to the many-fluid case is straightforward.

The first step consists in defining a Mukhanov-Sasaki variable for each fluid
present in the system. For scalar fields, this is known to be
$q_{a}=a\delta\phi_{a}^{(\mathrm{gi})}+z_{a}\Phi$ with
$z_{a}=a\phi^{\prime}_{a}/{\cal H}$, see Ref. [18]. In the case of a single
scalar field, we recall that $q=a\delta\phi^{(\mathrm{gi})}+z\Phi$, where, as
already introduced before,
$z=a\mathchoice{{\hbox{$\displaystyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=3.22221pt,depth=-2.57779pt}}}M_{\scriptscriptstyle{\mathrm{Pl}}}=a\phi^{\prime}/{\cal
H}$. In some sense, $z_{a}$ is a generalization of the second manner of
writing the $z$ variable (namely a generalization of $a\phi^{\prime}/{\cal H}$
and “not” of
$aM_{\scriptscriptstyle{\mathrm{Pl}}}\mathchoice{{\hbox{$\displaystyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule height=3.22221pt,depth=-2.57779pt}}}$).

For perfect fluids, to our knowledge, the question has not been studied as
thoroughly as in the case of scalar fields. We recall that, for a single fluid
$q=z{\cal R}=z\left[\Psi-{\cal H}v_{(\mathrm{gi})}\right]$, with
$z=a\mathchoice{{\hbox{$\displaystyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=3.22221pt,depth=-2.57779pt}}}M_{\scriptscriptstyle{\mathrm{Pl}}}/c_{{{}_{\mathrm{S}}}}$.
The question is then to define a variable $q_{(\alpha)}$ and a quantity
$z_{(\alpha)}$ for each fluid. One possibility for $z_{(\alpha)}$ would be
$z_{(\alpha)}=a\mathchoice{{\hbox{$\displaystyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{2\epsilon_{1}\,}$}\lower
0.4pt\hbox{\vrule
height=3.22221pt,depth=-2.57779pt}}}M_{\scriptscriptstyle{\mathrm{Pl}}}/c_{(\alpha)}$.
However, there is also another possibility which seems closer to what is done
in the case of scalar fields. Indeed the one-fluid definition of $z$ can also
be written as $z=a^{2}/\cal
H\mathchoice{{\hbox{$\displaystyle\sqrt{(\rho+p)/(p^{\prime}/\rho^{\prime})\,}$}\lower
0.4pt\hbox{\vrule
height=7.5pt,depth=-6.00003pt}}}{{\hbox{$\textstyle\sqrt{(\rho+p)/(p^{\prime}/\rho^{\prime})\,}$}\lower
0.4pt\hbox{\vrule
height=7.5pt,depth=-6.00003pt}}}{{\hbox{$\scriptstyle\sqrt{(\rho+p)/(p^{\prime}/\rho^{\prime})\,}$}\lower
0.4pt\hbox{\vrule
height=5.25pt,depth=-4.20003pt}}}{{\hbox{$\scriptscriptstyle\sqrt{(\rho+p)/(p^{\prime}/\rho^{\prime})\,}$}\lower
0.4pt\hbox{\vrule height=3.75pt,depth=-3.00002pt}}}$, where, here, $\rho$ and
$p$ are the total energy density and pressure, respectively, which
(obviously!) are also the energy density and the pressure of the fluid under
consideration since we assumed there is only one degree of freedom. In the
multi-fluid situation, this suggests the introduction of the quantity
$z_{(\alpha)}$ defined by

$\displaystyle z_{(\alpha)}=\frac{a^{2}}{\cal
H}\mathchoice{{\hbox{$\displaystyle\sqrt{\frac{\rho_{(\alpha)}+p_{(\alpha)}}{p^{\prime}_{(\alpha)}/\rho^{\prime}_{(\alpha)}}\,}$}\lower
0.4pt\hbox{\vrule
height=14.83333pt,depth=-11.86671pt}}}{{\hbox{$\textstyle\sqrt{\frac{\rho_{(\alpha)}+p_{(\alpha)}}{p^{\prime}_{(\alpha)}/\rho^{\prime}_{(\alpha)}}\,}$}\lower
0.4pt\hbox{\vrule
height=10.45833pt,depth=-8.3667pt}}}{{\hbox{$\scriptstyle\sqrt{\frac{\rho_{(\alpha)}+p_{(\alpha)}}{p^{\prime}_{(\alpha)}/\rho^{\prime}_{(\alpha)}}\,}$}\lower
0.4pt\hbox{\vrule
height=8.54166pt,depth=-6.83336pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\frac{\rho_{(\alpha)}+p_{(\alpha)}}{p^{\prime}_{(\alpha)}/\rho^{\prime}_{(\alpha)}}\,}$}\lower
0.4pt\hbox{\vrule height=8.54166pt,depth=-6.83336pt}}}.$ (4.63)

It is important to notice that $z_{(\alpha)}$ can no longer be expressed in
terms of $\epsilon_{1}$ because $\epsilon_{1}$ is now determined by the total
energy density and pressure. Then, one can define $q_{(\alpha)}$ by
$q_{(\alpha)}=z_{(\alpha)}{\cal R}_{(\alpha)}=z_{(\alpha)}\left[\Psi-{\cal
H}v_{(\alpha)}^{(\mathrm{gi})}\right]$.

Having defined the generalized Mukhanov-Sasaki variables for each fluid of the
system, the next step consists in establishing from the general equations of
Sec. 3, the equations satisfied by the $q_{a}$ and the $q_{(\alpha)}$. These
equations will obviously be coupled. Finally, one needs to take the large
scale limit, $k/{\cal H}\rightarrow\infty$, in order to guess the initial
conditions for the $q_{a}$ and the $q_{(\alpha)}$. Technically, this is
clearly a complicated task.

However, there exists a route which is equivalent and much easier. First of
all, we can remark that, initially, the couplings between the scalar fields
and the fluids can be neglected. Technically, this is due to the fact that the
physical momenta go to infinity in this regime and, therefore, become the
leading contribution in the equation of motion for the Mukhanov-Sasaki
variable. As a consequence, the $N_{\mathrm{field}}$ scalar fields and the
$N_{\mathrm{fluid}}$ fluids can be treated separately. The initial conditions
for a collection of $N_{\mathrm{field}}$ perturbed fields have been studied in
great details in the literature and it is standard to show that one has [19,
89, 24]

$\displaystyle
q_{a,k}^{A}(\eta)\rightarrow\frac{1}{\mathchoice{{\hbox{$\displaystyle\sqrt{2k\,}$}\lower
0.4pt\hbox{\vrule
height=6.94444pt,depth=-5.55559pt}}}{{\hbox{$\textstyle\sqrt{2k\,}$}\lower
0.4pt\hbox{\vrule
height=6.94444pt,depth=-5.55559pt}}}{{\hbox{$\scriptstyle\sqrt{2k\,}$}\lower
0.4pt\hbox{\vrule
height=4.8611pt,depth=-3.8889pt}}}{{\hbox{$\scriptscriptstyle\sqrt{2k\,}$}\lower
0.4pt\hbox{\vrule
height=3.47221pt,depth=-2.77779pt}}}}e^{ik(\eta-\eta_{\mathrm{ini}})}\delta_{a}^{A}.$
(4.64)

In this expression, it is worth noticing the presence of $\delta_{a}^{A}$
which means that, initially, the canonical variables are not “mixed”. Only the
time evolution of those variables, in presence of an interaction between them,
will be able to mix them.

The treatment of a collection of $N_{\mathrm{fluid}}$ fluids is less standard.
As already mentioned, in principle, one should establish the equations of
motion of the variables $q_{(\alpha)}$. However, in the large scale limit, we
expect the Bardeen potential not to play an important role. As a consequence,
one can introduce a simplified version of $q_{(\alpha)}$,
$\mathfrak{q}_{(\alpha)}$, defined by
$\mathfrak{q}_{(\alpha)}=-z_{(\alpha)}\mathcal{H}v_{(\alpha)}^{(\mathrm{gi})}$.
Assuming that each fluid has a constant equation of state, which is the case
of interest in this paper, we find that $\mathfrak{q}_{(\alpha)}$ obeys the
equation

$\mathfrak{q}_{(\alpha)}^{\prime\prime}+\mathfrak{q}_{(\alpha)}\left[w_{\alpha}k^{2}-\frac{z_{(\alpha)}^{\prime\prime}}{z_{(\alpha)}}+\mathcal{H}^{2}\epsilon_{1}\epsilon_{2}\right]=0.$
(4.65)

This equation is sufficient to fix the initial conditions in the sub-Hubble
regime and we take

$\displaystyle\mathfrak{q}_{(\alpha),k}^{A}\rightarrow\frac{1}{\mathchoice{{\hbox{$\displaystyle\sqrt{2c_{(\alpha)}k\,}$}\lower
0.4pt\hbox{\vrule
height=6.94444pt,depth=-5.55559pt}}}{{\hbox{$\textstyle\sqrt{2c_{(\alpha)}k\,}$}\lower
0.4pt\hbox{\vrule
height=6.94444pt,depth=-5.55559pt}}}{{\hbox{$\scriptstyle\sqrt{2c_{(\alpha)}k\,}$}\lower
0.4pt\hbox{\vrule
height=4.8611pt,depth=-3.8889pt}}}{{\hbox{$\scriptscriptstyle\sqrt{2c_{(\alpha)}k\,}$}\lower
0.4pt\hbox{\vrule
height=3.47221pt,depth=-2.77779pt}}}}e^{ic_{(\alpha)}k(\eta-\eta_{\mathrm{ini}})}\delta_{(\alpha)}^{A},$
(4.66)

where, as it was the case for the scalar fields, we have introduced “non-
mixing” initial conditions, see the presence of $\delta_{(\alpha)}^{A}$. When
the Bardeen potential is sub-dominant on sub-Hubble scales, one expects
$\mathfrak{q}_{(\alpha),k}^{A}\simeq q_{(\alpha),k}^{A}$.

Here, a comment is in order. We see that Eq. (4.66) is ill-defined when
$c_{(\alpha)}=0$, namely in the case where the fluid under consideration is
pressure-less, which may be relevant in the present context. In that case,
however, there is no oscillatory modes anymore and, moreover, the physical
wavelengths of the Fourier amplitudes are never inside the sonic horizon
(which simply vanishes), a necessary criterion to be able to single out well-
defined initial conditions. Carrying out the quantization of such a system
seems therefore difficult. On the other hand, we need to choose some initial
conditions in order to perform numerical calculations. Here, we will simply
take a small but non-vanishing pressure for the matter fluid:
$c_{\mathrm{m}}=0.01$. Actually, having $c_{\mathrm{m}}\neq 0$ is physically
well-justified since, in a realistic microscopic description, we expect non-
relativistic matter to be made of particles moving at small, but, crucially,
non-vanishing, velocities. Although it may be interesting to discuss this
question more deeply, at the practical level, we have noticed that the
numerical behavior of the system does not depend on those initial conditions.
Therefore the physical conclusions obtained in this paper will not be
sensitive to this issue.

Having chosen the initial conditions for the “canonical variables” $q_{a}$ and
$q_{(\alpha)}$, we now discuss how the initial conditions for the other
relevant quantities can be determined. This is done with the method presented
in Sec. 4.3.1. Notice, however, that the scalar field fluctuations and perfect
fluid fluctuations propagate with different speeds. The different perfect
fluid perturbations also propagate with different speeds since, a priori, the
fluids do not have the same equation of state. To deal with this issue, we
write $\Phi^{A}={\cal A}_{\Phi}^{A}e^{ic_{A}k(\eta-\eta_{\mathrm{ini}})}$,
with $c_{A}=(1,1,c_{1},c_{2})$. In the same manner, we write the two perturbed
scalar fields as $\delta\phi_{\ell,\mathrm{h}}^{A}={\cal
A}_{\ell,\mathrm{h}}^{A}e^{ic_{A}k(\eta-\eta_{\mathrm{ini}})}$, the two
perturbed energy densities as $\delta\rho_{1,2}^{A}={\cal
A}_{1,2}^{A}e^{ic_{A}k(\eta-\eta_{\mathrm{ini}})}$ and, finally, the two
velocities as $v_{1,2,A}^{(\mathrm{gi})}={\cal
A}_{v_{1},v_{2}}^{A}e^{ic_{A}k(\eta-\eta_{\mathrm{ini}})}$. It is interesting
to see that the scalar field fluctuations which, a priori, propagate $\propto
e^{ik(\eta-\eta_{\mathrm{ini}})}$ can also acquire modes $\propto
e^{ic_{1}k(\eta-\eta_{\mathrm{ini}})}$ and $\propto
e^{ic_{2}k(\eta-\eta_{\mathrm{ini}})}$ (the same remark could be done for
perfect fluids perturbations). Whether we initially populate the mode, say,
$\Phi^{1}={\cal A}_{\Phi}^{1}e^{ic_{1}k(\eta-\eta_{\mathrm{ini}})}$ is
precisely the choice of the initial conditions. In a standard situation, one
could also decide to populate this type of mode but, without interactions
between scalar fields and perfect fluids, this is not physically very relevant
and, in any case, does not correspond to the usual choice of initial
conditions. In the present context, however, this discussion is much more
relevant since the interactions are able to sustain such a mode.

The next step consists in inserting the above expressions of $\Phi^{A}$,
$\delta\phi_{\ell,\mathrm{h}}^{A}$, $\delta\rho_{1,2}^{A}$ and
$v_{1,2,A}^{\mathrm{(gi)}}$ in Eqs. (4.32), (4.2), (4.2), (4.2), (4.2), (4.2)
and (4.2), and neglect the interactions terms. As in Sec. 4.3.1, we consider
that the amplitudes ${\cal A}_{\Phi}^{A}$, ${\cal A}_{\ell,\mathrm{h}}^{A}$,
${\cal A}_{1,2}^{A}$ and ${\cal A}_{v_{1},v_{2}}^{A}$ are slowly varying
functions of time. Moreover, we use the initial conditions obtained from the
quantum-mechanical considerations presented earlier in this section. In
particular, the requirement (4.64) implies that

$\displaystyle{\cal
A}_{\ell}^{A}=\frac{K_{\ell}}{a}\delta_{\ell}^{A}-\frac{z_{\ell}}{a}{\cal
A}_{\Phi}^{A},\quad{\cal
A}_{\mathrm{h}}^{A}=\frac{K_{\mathrm{h}}}{a}\delta_{\mathrm{h}}^{A}-\frac{z_{\mathrm{h}}}{a}{\cal
A}_{\Phi}^{A},$ (4.67)

and we also have

$\displaystyle{\cal A}_{v_{1}}^{A}=-\frac{K_{1}}{{\cal
H}z_{1}}\delta_{1}^{A}+\frac{1}{\cal H}{\cal A}_{\Phi,1}^{A},\quad{\cal
A}_{v_{2}}^{A}=-\frac{K_{2}}{{\cal H}z_{2}}\delta_{2}^{A}+\frac{1}{\cal
H}{\cal A}_{\Phi,2}^{A},$ (4.68)

where the precise definitions of the background functions $z_{\ell}$,
$z_{\mathrm{h}}$, $z_{1}$ and $z_{2}$ have been introduced before. In fact, at
this stage, all we need to know is that these quantities are background
quantities and we will see that they do not play a crucial role in determining
the initial conditions. The quantities $K_{\ell}$, $K_{\mathrm{h}}$, $K_{1}$
and $K_{2}$, on the contrary, determine the scaling of the initial conditions
and are fixed by the quantum-mechanical considerations described above, namely
$K_{\ell}=K_{\mathrm{h}}=(2k)^{-1/2}$ and $K_{1,2}=(2c_{1,2}k)^{-1/2}$.

Let us start with the equation for $\Phi^{A}$. Substituting the previous WKB
ansatz in Eq. (4.32), one arrives at

$\displaystyle{\cal A}_{\Phi}^{A}{}^{\prime}+ic_{A}k\left\\{1+\frac{\cal
H}{ic_{A}k}+\frac{1}{2M_{\scriptscriptstyle{\mathrm{Pl}}}^{2}ic_{A}k}\left[\frac{\Phi_{\ell}^{\prime}z_{\ell}}{a}+\frac{\Phi_{\mathrm{h}}^{\prime}z_{\mathrm{h}}}{a}+\frac{a^{2}\rho_{1}(1+w_{1})}{\cal
H}+\frac{a^{2}\rho_{2}(1+w_{2})}{\cal H}\right]\right\\}{\cal A}_{\Phi}^{A}$
$\displaystyle=\frac{1}{2M_{\scriptscriptstyle{\mathrm{Pl}}}^{2}}\left[\frac{\phi_{\ell}^{\prime}}{a}K_{\ell}\delta_{\ell}^{A}+\frac{\phi^{\prime}_{\mathrm{h}}}{a}K_{\mathrm{h}}\delta_{\mathrm{h}}^{A}+\frac{a^{2}\rho_{1}(1+w_{1})}{{\cal
H}z_{1}}K_{1}\delta_{1}^{A}+\frac{a^{2}\rho_{2}(1+w_{2})}{{\cal
H}z_{2}}K_{2}\delta_{2}^{A}\right],$ (4.69)

from which one obtains

$\displaystyle{\cal
A}_{\Phi}^{A}\simeq\frac{1}{2M_{\scriptscriptstyle{\mathrm{Pl}}}^{2}ic_{A}k}\left[\frac{\phi_{\ell}^{\prime}}{a}K_{\ell}\delta_{\ell}^{A}+\frac{\phi^{\prime}_{\mathrm{h}}}{a}K_{\mathrm{h}}\delta_{\mathrm{h}}^{A}+\frac{a^{2}\rho_{1}(1+w_{1})}{{\cal
H}z_{1}}K_{1}\delta_{1}^{A}+\frac{a^{2}\rho_{2}(1+w_{2})}{{\cal
H}z_{2}}K_{2}\delta_{2}^{A}\right].$ (4.70)

We see that each component of ${\cal A}_{\Phi}^{A}$ is determined by the
corresponding fluid, that is to say the component ${\cal A}_{\Phi}^{\ell}$ is
determined by $K_{\ell}$, ${\cal A}_{\Phi}^{\mathrm{h}}$ is determined by
$K_{\mathrm{h}}$ and so on. As for the single fluid case, one has ${\cal
A}_{\Phi}^{A}\propto k^{-3/2}$. Then, using Eqs. (4.67) and (4.68), one
deduces that, at leading order, ${\cal A}_{\ell}^{A}\propto
k^{-1/2}\delta_{\ell}^{A}$, ${\cal A}_{\mathrm{h}}^{A}\propto
k^{-1/2}\delta_{\mathrm{h}}^{A}$, ${\cal A}_{v_{1}}^{A}\propto
k^{-1/2}\delta_{1}^{A}$ and ${\cal A}_{v_{2}}^{A}\propto
k^{-1/2}\delta_{2}^{A}$. At this stage, a comment is in order about the off-
diagonal terms of the components ${\cal A}_{\ell}^{A}$, ${\cal
A}_{\mathrm{h}}^{A}$, ${\cal A}_{v_{1}}^{A}$ and ${\cal A}_{v_{2}}^{A}$.
According to the previous considerations, these quantities are next-to-leading
order. For instance, according to Eq. (4.67), if $A\neq\ell$, one has ${\cal
A}_{\ell}^{A}\propto-z_{\ell}{\cal A}_{\Phi}^{A}/a$, which means a scaling
$\propto k^{-3/2}$ instead of $\propto k^{-1/2}$ for the diagonal component
($A=\ell$). However, one has to remember that these conclusions are based on
the quantization of the variable $\mathfrak{q}_{(\alpha)}$ introduced before,
which is equivalent to $q_{(\alpha)}$ only if the Bardeen contribution is