$S_{\mathrm{m}\gamma,A}$. It can be expressed as

$P_{S_{\mathrm{m}\gamma}}\propto|S_{\mathrm{m}\gamma,\ell}|^{2}+|S_{\mathrm{m}\gamma,\mathrm{h}}|^{2}\simeq|S_{\mathrm{m}\gamma,\ell}|^{2}\left\\{1+\left|\frac{C_{1,\mathrm{h}}}{C_{1,\ell}}\right|^{2}\frac{1}{R^{4}}\left[\frac{\phi_{\ell,\mathrm{ini}}^{2}}{\phi_{\mathrm{h}}(N_{\mathrm{m}})}\right]^{4}\right\\}\,.$
(5.97)

Quantity / Value |  Numeric |  Analytic $\big{(}N_{\mathrm{decay,h}}\big{)}$ |  error |  Analytic $\big{(}N_{\mathrm{decay,h}}^{{}_{\mathrm{WKB}}}\big{)}$ |  error |  Analytic $\big{(}N_{\mathrm{decay,h}}^{\mathrm{ana,sd}}\big{)}$ |  error
---|---|---|---|---|---|---|---
$\Re\mathrm{e}\,S_{\mathrm{m}\gamma,\ell}$ |  $-0.269$ |  $-0.251$ |  $6.5\%$ |  $-0.241$ |  $10.5\%$ |  $-0.195$ |  $27.4\%$
$\Im S_{\mathrm{m}\gamma,\ell}$ |  $-0.245$ |  $-0.229$ |  $6.6\%$ |  $-0.219$ |  $10.6\%$ |  $-0.177$ |  $27.5\%$
$\Re\mathrm{e}\,S_{\mathrm{m}\gamma,\mathrm{h}}\times 10^{3}$ |  $7.35$ |  $7.29$ |  $0.9\%$ |  $6.98$ |  $5.1\%$ |  $5.66$ |  $23.1\%$
$\Im S_{\mathrm{m}\gamma,\mathrm{h}}\times 10^{3}$ |  $6.60$ |  $6.33$ |  $4.2\%$ |  $6.06$ |  $8.3\%$ |  $4.91$ |  $25.7\%$
$P_{S_{\mathrm{m}\gamma}}$ |  $0.132$ |  $0.116$ |  $12.7\%$ |  $0.106$ |  $20.0\%$ |  $0.070$ |  $47.4\%$
Table 2: Numerical results versus analytical approximations following
different treatments to find the time of decay of the heavy field, for the
non-adiabatic fluctuations $S_{\mathrm{m}\gamma,A}$ and its power spectrum
$P_{S_{\mathrm{m}\gamma}}=\left|S_{\mathrm{m}\gamma,\ell}\right|^{2}+\left|S_{\mathrm{m}\gamma,\mathrm{h}}\right|^{2}$
at the end of reheating and beginning of the radiation-domination era. The
analytical approximations crucially depend on the value of the time of the
decay for the heavy field, so three possibilities are shown:
$N_{\mathrm{decay,h}}=-10.81$ for the value from the numerical evolution of
the background, $N_{\mathrm{decay,h}}^{{}_{\mathrm{WKB}}}=-10.35$ for the
analytical treatment with the WKB approximation presented in the App. A (see
Fig. 22) and $N_{\mathrm{decay,h}}^{\mathrm{ana,sd}}=-8.39$ with the simple
analytical treatment that was shown in Sec. 5.2.4 not to be very precise (see
Fig. 3).

From Eqs. (5.95), (5.96) and (5.97), we see that the crucial quantity that
needs to be calculated in order to evaluate the amplitude of the non-adiabatic
perturbations is $N_{\mathrm{decay},\mathrm{h}}$ or
$N_{\mathrm{end}}-N_{\mathrm{decay},\mathrm{h}}$. As discussed in Sec. 5.2.4,
the decay time of the heavy field is the time before which the amount of cold
dark matter should be negligible, and after which the heavy field should be
negligible. Therefore, in Sec. 5.2.4, this time was taken to be
$\rho_{\mathrm{h}}(N_{\mathrm{decay},\mathrm{h}})=\rho_{\mathrm{m}}(N_{\mathrm{decay},\mathrm{h}})$.
In this paper, given the importance of calculating
$N_{\mathrm{decay},\mathrm{h}}$, we have developed different methods to
estimate it. First, there are of course numerical calculations leading to a
semi-analytical determination of $P_{S_{\mathrm{m}\gamma}}$, namely
$P_{S_{\mathrm{m}\gamma}}$ is found analytically according to Eq. (5.97) but,
in this equation, $N_{\mathrm{decay},\mathrm{h}}$ is evaluated numerically. A
second way is to use the simple considerations presented in Sec. 5.2.4, which
leads to Eq. (5.45). Finally, we have also worked out a more precise approach
in Appendix A, based on the WKB approach, which leads to a third estimate of
the non-adiabatic perturbations. These results are summarized in Table. 2. We
see that the first method leads to an error of $\sim 10\%$, the WKB-based one
to an error of $\sim 20\%$ and, finally, the estimate of Eq. (5.45) to a $\sim
50\%$ error. In the case of the first method, since
$N_{\mathrm{decay},\mathrm{h}}$ is calculated exactly, the $10\%$ error
originates from the condition (5.92). Since it is difficult to see how other
analytical approximations could be designed, it represents, to some extent,
the incompressible error of any analytical approach. Note that these numbers
apply to the set of parameters
($R=5,\Gamma_{\mathrm{hm}}=10^{-5}M_{\scriptscriptstyle{\mathrm{Pl}}}$) that
was studied in detail in this section. In Sec. 6, we explore other values of
these parameters for which the analytical approach may be less precise, in
which case the numerical resolution appears to be the only reliable method to
accurately predict the amount of non-adiabatic perturbations at the end of
inflation.

##  6 Comparison with the standard approach

In this section, we compare our (numerical and analytical) calculations of the
non-adiabatic perturbations to the calculations of Refs. [15, 16] and Ref.
[17], which also derive an analytical formula for $P_{S_{\mathrm{m}\gamma}}$
in the double inflation scenario.

Before carrying out this comparison, we first recall, very briefly, how the
results of Refs. [15, 16] and Ref. [17] are obtained. Since the details of
reheating and of the decay process are neglected in those references, each
fluid is always individually conserved and the amplitude of the non-adiabatic
modes can be defined by Eq. (3.30), namely

$\displaystyle
S_{\mathrm{m}\gamma}=\frac{\delta\rho_{\mathrm{m}}}{\rho_{\mathrm{m}}}-\frac{3}{4}\frac{\delta\rho_{\gamma}}{\rho_{\gamma}}\simeq\frac{\delta\rho_{\mathrm{m}}}{\rho_{\mathrm{m}}},$
(6.1)

where the last expression is again the fact, discussed in the previous
section, that the amplitude of the non-adiabatic perturbation is dominated by
what happens at the heavy field decay. The question is then to relate
$\delta\rho_{\mathrm{m}}/\rho_{\mathrm{m}}$ to
$\delta{\rho_{\mathrm{h}}}/\rho_{\mathrm{h}}$. In the literature,
$\delta\phi_{\mathrm{h}}/\phi_{\mathrm{h}}$ is taken to be constant during the
oscillations of the heavy field (because $\delta\phi_{\mathrm{h}}$ and
$\phi_{\mathrm{h}}$ obey the same equation of motion in this regime) which
covers the phase between the heavy field decay and the end of inflation. At
the end of inflation, one can thus write

$\displaystyle\frac{\delta\rho_{\mathrm{m}}}{\rho_{\mathrm{m}}}\sim\frac{\delta\rho_{\mathrm{h}}}{\rho_{\mathrm{h}}}\sim
2\frac{\delta\phi_{\mathrm{h}}}{\phi_{\mathrm{h}}},$ (6.2)

which comes from the continuity of energy density and $\rho_{\mathrm{h}}\sim
m_{\mathrm{h}}^{2}\langle\phi_{\mathrm{h}}\rangle^{2}/2$.

Then, our next goal is to make the connection between the previous
considerations and the super Hubble slow-roll solution given by Eq. (5.58),
which can be written as

$\displaystyle\delta\phi_{\mathrm{h}}=\frac{C_{1}}{H}\dot{\phi}_{\mathrm{h}}+2HC_{3}\frac{V_{\ell}}{V_{\ell}+V_{\mathrm{h}}}\dot{\phi}_{\mathrm{h}}=-2C_{1}\frac{M_{\scriptscriptstyle{\mathrm{Pl}}}^{2}}{\phi_{\mathrm{h}}}\frac{V_{\mathrm{h}}}{V_{\mathrm{h}}+V_{\ell}}-\frac{2}{3}C_{3}\frac{\partial
V_{\mathrm{h}}}{\partial\phi_{\mathrm{h}}}\frac{V_{\ell}}{V_{\mathrm{h}}+V_{\ell}}.$
(6.3)

From this expression, we conclude that, if the light field dominates
($V_{\ell}\gg V_{\mathrm{h}}$), then the first branch of the solution $\propto
C_{1}$ can be ignored and one has

$\displaystyle\frac{\delta\phi_{\mathrm{h}}}{\phi_{\mathrm{h}}}\simeq-\frac{2}{3}C_{3}m_{\mathrm{h}}^{2}.$
(6.4)

As a consequence, the quantity $\delta\phi_{\mathrm{h}}/\phi_{\mathrm{h}}$
also remains constant once the light field has started to dominate, which
occurs slightly before the decay of the heavy field (if one can really define
a time of decay for the heavy field, see the discussion in the previous
sections). In some sense, these considerations allow us to fix the value of
the constant $\delta\phi_{\mathrm{h}}/\phi_{\mathrm{h}}$ at the end of
inflation in Eq. (6.2) and one can write

$\displaystyle
S_{\mathrm{m}\gamma}\simeq\frac{\delta\rho_{\mathrm{m}}}{\rho_{\mathrm{m}}}\simeq-\frac{4}{3}m_{\mathrm{h}}^{2}C_{3}.$
(6.5)

One recovers the fact that non-adiabatic perturbations solely depend on the
constant $C_{3}$.

Then, the constant $C_{3}$ is evaluated with the help of Eq. (5.67). In this
formula, the quantities $\delta\phi_{\ell}$ and $\delta\phi_{\mathrm{h}}$ are
written

$\displaystyle\delta\phi_{\ell}=\frac{H_{*}}{\mathchoice{{\hbox{$\displaystyle\sqrt{2k^{3}\,}$}\lower
0.4pt\hbox{\vrule
height=8.74889pt,depth=-6.99915pt}}}{{\hbox{$\textstyle\sqrt{2k^{3}\,}$}\lower
0.4pt\hbox{\vrule
height=8.74889pt,depth=-6.99915pt}}}{{\hbox{$\scriptstyle\sqrt{2k^{3}\,}$}\lower
0.4pt\hbox{\vrule
height=6.14998pt,depth=-4.92001pt}}}{{\hbox{$\scriptscriptstyle\sqrt{2k^{3}\,}$}\lower
0.4pt\hbox{\vrule
height=4.7611pt,depth=-3.8089pt}}}}e_{\ell},\quad\delta\phi_{\mathrm{h}}=\frac{H_{*}}{\mathchoice{{\hbox{$\displaystyle\sqrt{2k^{3}\,}$}\lower
0.4pt\hbox{\vrule
height=8.74889pt,depth=-6.99915pt}}}{{\hbox{$\textstyle\sqrt{2k^{3}\,}$}\lower
0.4pt\hbox{\vrule
height=8.74889pt,depth=-6.99915pt}}}{{\hbox{$\scriptstyle\sqrt{2k^{3}\,}$}\lower
0.4pt\hbox{\vrule
height=6.14998pt,depth=-4.92001pt}}}{{\hbox{$\scriptscriptstyle\sqrt{2k^{3}\,}$}\lower
0.4pt\hbox{\vrule height=4.7611pt,depth=-3.8089pt}}}}e_{\mathrm{h}},$ (6.6)

where we recall that $H_{*}$ denotes the value of the Hubble parameter at
horizon crossing and $e_{i}$ ($i=\ell,\mathrm{h}$) are classical random
Gaussian processes with $\langle e_{i}({\bm{k}})\rangle=0$ and $\langle
e_{i}({\bm{k}})e_{j}^{*}({\bm{k}}^{\prime})\rangle=\delta_{ij}\delta({\bm{k}}-{\bm{k}}^{\prime})$.
It follows that Eq. (5.67) can be re-written as

$\displaystyle
S_{\mathrm{m}\gamma}=\mathchoice{{\hbox{$\displaystyle\sqrt{\frac{2}{k^{3}}\,}$}\lower
0.4pt\hbox{\vrule
height=10.81888pt,depth=-8.65514pt}}}{{\hbox{$\textstyle\sqrt{\frac{2}{k^{3}}\,}$}\lower
0.4pt\hbox{\vrule
height=7.58609pt,depth=-6.06891pt}}}{{\hbox{$\scriptstyle\sqrt{\frac{2}{k^{3}}\,}$}\lower
0.4pt\hbox{\vrule
height=5.60275pt,depth=-4.48222pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\frac{2}{k^{3}}\,}$}\lower
0.4pt\hbox{\vrule
height=5.60275pt,depth=-4.48222pt}}}H_{*}\left(\frac{e_{\mathrm{h}}}{\phi_{\mathrm{h}}}-R^{2}\frac{e_{\ell}}{\phi_{\ell}}\right).$
(6.7)

As a consequence, using Eqs. (5.16), (5.18) and (5.21), one arrives at

$\displaystyle
P_{S_{\mathrm{m}\gamma}}=\left|S_{\mathrm{m}\gamma}\right|^{2}=\frac{1}{3k^{3}}\frac{m_{\ell}^{2}}{M_{\scriptscriptstyle{\mathrm{Pl}}}^{2}}\left[1+(R^{2}-1)\sin^{2}\theta_{*}\right]\left(\frac{1}{\sin^{2}\theta_{*}}+\frac{R^{4}}{\cos^{2}\theta_{*}}\right),$
(6.8)

an equation that should be compared to Eq. (5.97). Numerically, this equation
gives (we use $\theta_{*}=0.962$ from $N_{*}=4.67$ determined numerically and
$-N_{*}+N_{\mathrm{end}}\simeq C/\mathrm{cos}^{2}\theta_{*}$)
$P_{S_{\mathrm{m}\gamma}}=0.072$, to be compared with the numerical result
$P_{S_{\mathrm{m}\gamma}}=0.132$ and the analytical approximations presented
in this work: the relative error found with Eq. (6.8) is $45.5\%$, which is
comparable to our worst analytical estimate with
$N_{\mathrm{decay,h}}^{\mathrm{ana,sd}}$, see Table. 2.

Therefore, we conclude that modeling the reheating mechanism is crucial in
order to get an accurate prediction for non-adiabatic perturbations after the
end of multi-field inflation, in the post-inflationary Universe. Ignoring the
details of reheating can indeed lead to an error as big as $\sim 50\%$ (or
worse for some other values of the parameters, see below)

It is also interesting to compare the amplitude of adiabatic and non-adiabatic
perturbations. In the radiation-dominated era, we have
$\zeta_{\gamma}\simeq\zeta_{\ell}^{\mathrm{ad}}=-C_{1}$, see Eq. (5.60). Using
Eq. (5.66), one has

$\displaystyle\zeta_{\gamma}=\frac{1}{2M_{\scriptscriptstyle{\mathrm{Pl}}}^{2}}\frac{H_{*}}{\mathchoice{{\hbox{$\displaystyle\sqrt{2k^{3}\,}$}\lower
0.4pt\hbox{\vrule
height=8.74889pt,depth=-6.99915pt}}}{{\hbox{$\textstyle\sqrt{2k^{3}\,}$}\lower
0.4pt\hbox{\vrule
height=8.74889pt,depth=-6.99915pt}}}{{\hbox{$\scriptstyle\sqrt{2k^{3}\,}$}\lower
0.4pt\hbox{\vrule
height=6.14998pt,depth=-4.92001pt}}}{{\hbox{$\scriptscriptstyle\sqrt{2k^{3}\,}$}\lower
0.4pt\hbox{\vrule
height=4.7611pt,depth=-3.8089pt}}}}\left(\phi_{\mathrm{h}}e_{\mathrm{h}}+\phi_{\ell}e_{\ell}\right),$
(6.9)

from which we deduce

$\displaystyle
P_{\zeta_{\gamma}}=\frac{C^{2}}{3k^{3}}\frac{m_{\ell}^{2}}{M_{\scriptscriptstyle{\mathrm{Pl}}}^{2}}\left[1+(R^{2}-1)\sin^{2}\theta_{*}\right]\frac{(\tan\theta_{*})^{4/(R^{2}-1)}}{\cos^{4}\theta_{*}}.$
(6.10)

It follows, using the expression (5.27) of the constant $C$, that the non-
adiabatic to adiabatic power spectrum ratio can be expressed as

$\displaystyle\frac{P_{S_{\mathrm{m}\gamma}}}{P_{\zeta_{\gamma}}}=\frac{16M_{\scriptscriptstyle{\mathrm{Pl}}}^{4}}{\phi_{\ell}^{4}|_{\mathrm{ini}}}\left(\tan\theta_{*}\right)^{4/(1-R^{2})}\left(\frac{\cos^{4}\theta_{*}}{\sin^{2}\theta_{*}}+R^{4}\cos^{2}\theta_{*}\right).$
(6.11)

For our fiducial parameters used above (leading to $\theta_{*}\simeq 0.962$),
one finds $P_{S_{\mathrm{m}\gamma}}/P_{\zeta_{\gamma}}\simeq 0.75$. This
should be compared to our results,
$P_{S_{\mathrm{m}\gamma}}/P_{\zeta_{\gamma}}\simeq 1.65$ numerically and
$P_{S_{\mathrm{m}\gamma}}/P_{\zeta_{\gamma}}\simeq(1.60,1.46,0.96)$
analytically with the three possible choices of $N_{\mathrm{decay,h}}$ as
explained in Sec. 5.3.4. Notice that the precision on the non-adiabatic power
spectrum to adiabatic power spectrum ratio seems better than the numbers
quoted in Tables. 1 and 2 (for instance, it is $3\%$ with an exact
$N_{\mathrm{decay,h}}$). This is due to the fact that both quantities are
underestimated and, as a consequence, the error in their ratio is smaller than
the error in each of them.

(a) Ratio $P_{S_{\mathrm{m}\gamma}}/P_{\zeta_{\gamma}}$ as a function of $R$

(b) Ratio $P_{S_{\mathrm{m}\gamma}}/P_{\zeta_{\gamma}}$ as a function of
$\Gamma_{\mathrm{hm}}$ (in $M_{\scriptscriptstyle{\mathrm{Pl}}}$ units)

Figure 18: Ratio $P_{S_{\mathrm{m}\gamma}}/P_{\zeta_{\gamma}}$ as a function
of the parameter $R$, panel (a), and as a function of the parameter
$\Gamma_{\mathrm{hm}}$ (in $M_{\scriptscriptstyle{\mathrm{Pl}}}$ units), panel
(b). The blue line represents the exact result, obtained by means of numerical
calculations. The orange line corresponds to the analytical approach developed
in this paper with $N_{\mathrm{decay,h}}$ determined numerically. Finally, the
green line corresponds to the standard formulation, see Eq. (6.11).

In Fig. 18, we have represented the ratio
$P_{S_{\mathrm{m}\gamma}}/P_{\zeta_{\gamma}}$ versus $R$, for a fixed
$\Gamma_{\mathrm{hm}}$ (taken at its fiducial value
$\Gamma_{\mathrm{hm}}=1.0\times 10^{-5}M_{\scriptscriptstyle{\mathrm{Pl}}}$ )
and versus $\Gamma_{\mathrm{hm}}$ for a fixed $R$ (also taken at its fiducial
value $R=5$) for three cases: the blue curve corresponds to the exact
(numerical) result, the orange curve corresponds to the semi-analytical result
obtained in this paper and the green curve is obtained from Eq. (6.11). In
panel (a), we see that our semi-analytical result reproduces well the exact
result while the standard formulation poorly performs, especially for
$R\gtrsim 8$. The two drops that can be seen on the orange curve are due to
the way $N_{\mathrm{decay,h}}$ is calculated. $N_{\mathrm{decay,h}}$ is
obtained by determining when $\rho_{\mathrm{h}}$ and $\rho_{\mathrm{m}}$
intersects. Since these two quantities oscillate in this regime, the
intersection can occur around one specific oscillation and then, suddenly,
when the value of $R$ changes, jumps to another oscillation resulting in the
behavior observed in Fig. 18. In panel (b), we notice that the agreement
between the exact and analytical results (blue and orange lines) is still
quite good but, admittedly, less good than in panel (a). On the other hand,
the standard formula (green line) still poorly performs. Notice that the
latter is almost independent of $\Gamma_{\mathrm{hm}}$, a slight dependence
being nevertheless present through the value of $\theta_{*}$. We also remark
that the agreement between the numerical and analytical results deteriorates
when $\Gamma_{\mathrm{hm}}$ increases compared to our fiducial value
$\Gamma_{\mathrm{hm}}=1.0\times 10^{-5}M_{\scriptscriptstyle{\mathrm{Pl}}}$.
However, when this is the case, one enters a regime where, in the Klein-Gordon
equation, the contribution of the term proportional to $\Gamma_{\mathrm{hm}}$
is no longer negligible compared to $H_{\mathrm{ini}}$ at the onset of
inflation. This means that we are in a situation which is reminiscent of warm
inflation where the considerations presented in this work are no longer valid.
Moreover, although the departure is less visible, it seems also that the
analytical approximation is less accurate when $\Gamma_{\mathrm{hm}}$
decreases too much. But again, this can be understood: if
$\Gamma_{\mathrm{hm}}$ is very small, then the instantaneous decay
approximation for the heavy field is even less justified, and even though one
can still define the time $N_{\mathrm{decay,h}}$ as the time when
$\rho_{\mathrm{h}}=\rho_{\mathrm{m}}$, it becomes less relevant to use it as a
particularly interesting time. This shows once more that, in the general case,
the system should be solved numerically in order to derive reliable
predictions for the cosmological observables.

(a) Full numerical result.

(b) Analytical (approximated) result.

Figure 19: Contour levels of the ratio
$P_{S_{\mathrm{m}\gamma}}/P_{\zeta_{\gamma}}$ as a function of $R$ and
$\Gamma_{\mathrm{hm}}$ (in $M_{\scriptscriptstyle{\mathrm{Pl}}}$ units). The
figure in the left panel is obtained by numerical calculations and is the
results of $10000$ simulations (namely $100$ values for $R$ and $100$ values
for $\Gamma_{\mathrm{hm}}$). The right panel is obtained by means of the
analytical approximations discussed in this work.

In Fig. 19, we present a comparison of the ratios
$P_{S_{\mathrm{m}\gamma}}/P_{\zeta_{\gamma}}$, as a function of the parameters
$R$ and $\Gamma_{\mathrm{hm}}$, obtained numerically and analytically. This
confirms that our analytical approximations correctly reproduce the exact
results, especially for $R\lesssim 7$. For larger values of $R$, the agreement
is less good and we observe a complicated structure in parameter space. As
already mentioned above, this structure is due to difficulties in handling
$N_{\mathrm{decay,h}}$. This structure is the two-dimensional manifestation of
the drops observed in the orange line in Fig. 18, panel (a). Fig. 19 also
confirms that the dependence of $P_{S_{\mathrm{m}\gamma}}/P_{\zeta_{\gamma}}$
is stronger in the $R$-direction than in the $\Gamma_{\mathrm{hm}}$, a
property which can also be noticed in Fig. 18. However, numerical calculations
indicate that this dependence is less flat than suggested by analytical
approximations.

We conclude by stressing again that a detailed modeling of the reheating epoch
can affect the predictions of multi-field scenarios and in a regime in which
non-adiabatic perturbations represent a non-negligible fraction of the
adiabatic ones.

##  7 Conclusions

In multi-field inflation, because curvature perturbations are not necessarily
conserved on large scales, one expects the predictions of a model to depend on
what happens during reheating, and the main goal of this article was to
investigate whether this effect can indeed be relevant. Moreover, beyond being
a question of principle, this issue is also important, and timely, because
future, high-accuracy, experiments could maybe, open up a new window on the
micro-physics of inflation.

In order to carry out the aforementioned task at the practical level, we have
introduced a detailed description of reheating which explicitly takes into
account the interactions between the inflaton fields and their decay products.
As an illustration of a property that can be affected by the physical
processes at play during reheating, we have chosen to consider the amplitude
of non-adiabatic perturbations in the post-inflationary Universe.

We have then developed new numerical codes, and new analytical techniques,
allowing us to track the behavior of the perturbations throughout reheating.
Concretely, these tools were used in the context of a specific model, namely
double inflation.

The main result of this study is that, indeed, in the context of multi-field
inflation, the predictions can be quite sensitive to the details of what is
going on during reheating. In the case studied here, we have estimated that
different modelings of the reheating phase can result in predictions that
differ by almost $\sim 50\%$. Of course, this number is to be taken with a
grain of salt given that it was obtained for a particular model and for a
particular observable. However, we think it illustrates well the main
conclusion of this article. Another conclusion is that it is possible to
develop analytical techniques which can reasonably approximate the exact
calculations, at the $\sim 10\%$ level. However, for a more realistic and,
therefore, more complicated model (with, for instance, potential and kinetic
interactions between the inflaton fields, more complicated potentials, more
decay channels …), it is very likely that no analytical approach will be able
to reach a good level of accuracy. In this respect, the codes that we have
developed for this paper should be very useful for future investigations in
multi-field inflation.

As a final remark, let us notice that it would be interesting to test the
robustness of the above mentioned lessons by exploring whether other effects
in double inflation, new models and/or new observables are also dependent on
the details of reheating at the level found before. For instance, one
possibility would be that the inflaton fields decay in other scalar fields
rather than in perfect fluids. For example, the light field could decay into
the field $\chi$ with $-{\cal
L}_{\mathrm{matter}}=(\partial\chi)^{2}/2+V(\chi)$ and
$V(\chi)\propto\mathrm{exp}\left[-\mathchoice{{\hbox{$\displaystyle\sqrt{3(1+w)\,}$}\lower
0.4pt\hbox{\vrule
height=7.5pt,depth=-6.00003pt}}}{{\hbox{$\textstyle\sqrt{3(1+w)\,}$}\lower
0.4pt\hbox{\vrule
height=7.5pt,depth=-6.00003pt}}}{{\hbox{$\scriptstyle\sqrt{3(1+w)\,}$}\lower
0.4pt\hbox{\vrule
height=5.25pt,depth=-4.20003pt}}}{{\hbox{$\scriptscriptstyle\sqrt{3(1+w)\,}$}\lower
0.4pt\hbox{\vrule
height=3.75pt,depth=-3.00002pt}}}\chi/M_{\scriptscriptstyle{\mathrm{Pl}}}\right]$,
where $w$ is a free parameter [94]. This potential is used in the model of
power-law inflation. An important feature of this model is that the equation
of state of $\chi$ is constant and given by $w$ (hence the name). In the
context of power-law inflation, $w$ is chosen sufficiently close to minus one
which makes the potential very flat as required for an inflationary model.
However, if $w=1/3$, and, therefore, if the potential is rather steep, then
when the inflaton fields have decayed away and the energy budget of the
Universe is dominated by the field $\chi$, a radiation-dominated epoch starts.
In other words, a field $\chi$ with the above potential provides an example of
a microscopic description of radiation. Then, the interaction between $\chi$
and the inflaton field could typically be taken as $-{\cal
L}_{\mathrm{int}}=\sigma\phi_{\ell}\chi^{2}$, where $\sigma$ is a constant of
dimension one. This Lagrangian describes the decay of one “$\phi_{\ell}$\-
particle” into two “$\chi$-particles”. This represents an interesting example
where the Lagrangian (5.1) is entirely explicit. Then, the calculations of the
present article could be carried out again with this new model in order to see
if similar conclusions are reached. We hope to return to these questions soon.

## Acknowledgments

We would like to thank David Kaiser for insightful discussions after the
publication of the first version of this manuscript.

##  Appendix A Background evolution at heavy field decay

In this appendix, we consider again the problem of modeling the behavior of
the background around the decay time of the heavy field. This question was
already treated in Sec. 5.2.4 and our aim in this appendix is to improve this
treatment. Here, we study in detail the behavior of the system when the heavy
field leaves the slow-roll regime and starts to oscillate and to decay. As
argued before, see Sec. 5.3.4, this is important because tracking what happens
to the background in this regime with good accuracy turns out to be crucial in
order to predict the level of non-adiabatic perturbations that remains after
the end of inflation.

###  A.1 Time-dependent frequency

Figure 20: Squared pulsation $\omega_{\mathrm{h}}^{2}$, defined in Eq. (A.3),
versus number of e-folds: exact numerical result (blue line), quadratic
approximation $\omega_{\mathrm{h,ana,quad}}^{2}$ (orange line) given by Eq.
(A.7) and that becomes roughly valid after the time represented by the orange
vertical dashed line, and linear approximation
$\omega_{\mathrm{h,ana,lin}}^{2}$ (green line) given by Eq. (A.13) and that is
roughly valid until the time represented by the green vertical dashed line.
Both approximations are faithful in the time domain between the two vertical
dashed lines, which also coincides roughly with the time at which the WKB
solution becomes correct (see the next figure).

Since we know that the heavy field decays and oscillates, as it was already
done in Sec. 5.2.4, it is especially convenient to write this quantity as

$\displaystyle\phi_{\mathrm{h}}(t)=\phi_{\mathrm{h},\mathrm{p}}\left(\frac{a_{\mathrm{p}}}{a}\right)^{3/2}e^{-\Gamma_{\mathrm{hm}}(t-t_{\mathrm{p}})/4}\,g_{\mathrm{h}}(t),$
(A.1)

with “$\mathrm{p}$” meaning evaluated at a particular time. The quantity
$\phi_{\mathrm{h,p}}$ is the vacuum expectation value of the heavy field at
this particular time. This particular time will be specified soon. Eq. (A.1)
defines the new function $g_{\mathrm{h}}(t)$. Then, using Eq. (5.7), it is
straightforward to show that $g_{\mathrm{h}}(t)$ obeys the equation [see also
Eqs. (5.40) and (5.41) where these quantities were already discussed in the
main text],

$\displaystyle\ddot{g}_{\mathrm{h}}(t)+\omega^{2}_{\mathrm{h}}(t)g_{\mathrm{h}}(t)=0,$
(A.2)

with

$\displaystyle\omega^{2}_{\mathrm{h}}(t)=m_{\mathrm{h}}^{2}\left(1-\frac{3}{2}\epsilon_{1}\frac{H^{2}}{m_{\mathrm{h}}^{2}}-\frac{9}{4}\frac{H^{2}}{m_{\mathrm{h}}^{2}}-\frac{1}{16}\frac{\Gamma_{\mathrm{hm}}^{2}}{m_{\mathrm{h}}^{2}}-\frac{3}{4}\frac{H}{m_{\mathrm{h}}}\frac{\Gamma_{\mathrm{hm}}}{m_{\mathrm{h}}}\right).$
(A.3)

This effective frequency is represented in Fig. 20 (blue solid line).
Analytically, viewed as a function of cosmic time, the explicit form of
$\omega^{2}_{\mathrm{h}}(t)$ is complicated and, as a consequence, Eq. (A.2)
cannot be integrated exactly. We see that, at the beginning of inflation,
$\omega^{2}_{\mathrm{h}}(t)$ is negative, vanishes at a time that we denote in
the following by $t_{\mathrm{osc}}$, which is a turning point, and then
becomes positive which entails the presence of oscillations in the solution.

The position of the turning point can be determined with the help of the
following considerations. In the neighborhood of the turning point, one can
use an approximation for $\omega_{\mathrm{h}}(t)$ and write
$\omega^{2}_{\mathrm{h}}\simeq m_{\mathrm{h}}^{2}-9H^{2}/4$. The turning point
is therefore located at a time when $H_{\mathrm{osc}}\simeq 2m_{\mathrm{h}}/3$
and using the expression (5.21), it occurs at a value of
$\tan\theta_{\mathrm{osc}}$ which is a solution to the following equation

$\displaystyle\left(\tan\theta_{\mathrm{osc}}\right)^{2/(R^{2}-1)}\left(1+R^{2}\tan^{2}\theta_{\mathrm{osc}}\right)=\frac{2R^{2}}{3C}.$
(A.4)

This equation is a polynomial equation in $\tan^{2}\theta_{\mathrm{osc}}$
which, unfortunately, cannot be solved exactly. However, one can use Newton’s
method to find an approximate solution. Starting from the initial guess
$\tan\theta_{\mathrm{osc}}^{2}\sim R^{-2}$, the first iteration leads to

$\tan^{2}\theta_{\mathrm{osc}}\simeq\frac{3-R^{2}+2R^{2}(R^{2}-1)R^{2/(R^{2}-1)}/(3C)}{R^{2}(1+R^{2})}\,.$
(A.5)

The exact (numerical) solution to Eq. (A.4), for the fiducial parameters, is
$\tan^{2}\theta_{\mathrm{osc}}\simeq 0.0101$ to be compared with the solution
(A.5), $\tan^{2}\theta_{\mathrm{osc}}\simeq 0.0104$, and corresponding to
$N_{\mathrm{osc}}-N_{\mathrm{end}}=-13.3$, to be compared with the exact time
found numerically, $N_{\mathrm{osc}}-N_{\mathrm{end}}=-12.89$.

As noticed above, the form of $\omega_{\mathrm{h}}^{2}(t)$ is complicated.
However, depending on the regimes considered, some useful approximations can
be found. For instance, in the regime dominated by the light field, the Hubble
parameter simplifies and is given by Eq. (5.35),
$H^{2}=2m_{\ell}^{2}(N_{\mathrm{end}}-N)/3$. Moreover, given that
$\mathrm{d}N=H\mathrm{d}t$, the relationship between the number of e-folds and
the cosmic time can be expressed as [this equation was already derived, see
Eq. (5.46), and we reproduce it here for convenience]

$\displaystyle
t-t_{\mathrm{p}}=-\frac{\mathchoice{{\hbox{$\displaystyle\sqrt{6\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{6\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{6\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{6\,}$}\lower
0.4pt\hbox{\vrule
height=3.22221pt,depth=-2.57779pt}}}}{m_{\ell}}\left[\left(N_{\mathrm{end}}-N\right)^{1/2}-\left(N_{\mathrm{end}}-N_{\mathrm{p}}\right)^{1/2}\right].$
(A.6)

As a consequence, inserting those equations in Eq. (A.3), one finds that
$\omega^{2}_{\mathrm{h}}(t)$ reduces to a quadratic polynomial in time, namely

$\displaystyle\omega^{2}_{\mathrm{h}}(t)\simeq\bar{c}+\bar{b}\,m_{\ell}\left(t-t_{\mathrm{p}}\right)+\bar{a}\,m_{\ell}^{2}\left(t-t_{\mathrm{p}}\right)^{2},$
(A.7)

with

$\displaystyle\bar{a}$ $\displaystyle=-\frac{1}{4}m_{\ell}^{2}\,,$ (A.8)
$\displaystyle\bar{b}$
$\displaystyle=\frac{1}{4}m_{\ell}\Gamma_{\mathrm{hm}}+\mathchoice{{\hbox{$\displaystyle\sqrt{\frac{3}{2}\,}$}\lower
0.4pt\hbox{\vrule
height=9.66666pt,depth=-7.73337pt}}}{{\hbox{$\textstyle\sqrt{\frac{3}{2}\,}$}\lower
0.4pt\hbox{\vrule
height=6.76666pt,depth=-5.41336pt}}}{{\hbox{$\scriptstyle\sqrt{\frac{3}{2}\,}$}\lower
0.4pt\hbox{\vrule
height=4.83331pt,depth=-3.86667pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\frac{3}{2}\,}$}\lower
0.4pt\hbox{\vrule
height=4.83331pt,depth=-3.86667pt}}}m_{\ell}^{2}\mathchoice{{\hbox{$\displaystyle\sqrt{N_{\mathrm{end}}-N_{\mathrm{p}}\,}$}\lower
0.4pt\hbox{\vrule
height=6.83331pt,depth=-5.46667pt}}}{{\hbox{$\textstyle\sqrt{N_{\mathrm{end}}-N_{\mathrm{p}}\,}$}\lower
0.4pt\hbox{\vrule
height=6.83331pt,depth=-5.46667pt}}}{{\hbox{$\scriptstyle\sqrt{N_{\mathrm{end}}-N_{\mathrm{p}}\,}$}\lower
0.4pt\hbox{\vrule
height=4.78333pt,depth=-3.82668pt}}}{{\hbox{$\scriptscriptstyle\sqrt{N_{\mathrm{end}}-N_{\mathrm{p}}\,}$}\lower
0.4pt\hbox{\vrule height=3.41666pt,depth=-2.73334pt}}}\,,$ (A.9)
$\displaystyle\bar{c}$
$\displaystyle=m_{\mathrm{h}}^{2}-\frac{1}{2}m_{\ell}^{2}-\frac{\Gamma^{2}_{\mathrm{hm}}}{16}-\mathchoice{{\hbox{$\displaystyle\sqrt{\frac{3}{8}\,}$}\lower
0.4pt\hbox{\vrule
height=9.66666pt,depth=-7.73337pt}}}{{\hbox{$\textstyle\sqrt{\frac{3}{8}\,}$}\lower
0.4pt\hbox{\vrule
height=6.76666pt,depth=-5.41336pt}}}{{\hbox{$\scriptstyle\sqrt{\frac{3}{8}\,}$}\lower
0.4pt\hbox{\vrule
height=4.83331pt,depth=-3.86667pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\frac{3}{8}\,}$}\lower
0.4pt\hbox{\vrule
height=4.83331pt,depth=-3.86667pt}}}\Gamma_{\mathrm{hm}}m_{\ell}\mathchoice{{\hbox{$\displaystyle\sqrt{N_{\mathrm{end}}-N_{\mathrm{p}}\,}$}\lower
0.4pt\hbox{\vrule
height=6.83331pt,depth=-5.46667pt}}}{{\hbox{$\textstyle\sqrt{N_{\mathrm{end}}-N_{\mathrm{p}}\,}$}\lower
0.4pt\hbox{\vrule
height=6.83331pt,depth=-5.46667pt}}}{{\hbox{$\scriptstyle\sqrt{N_{\mathrm{end}}-N_{\mathrm{p}}\,}$}\lower
0.4pt\hbox{\vrule
height=4.78333pt,depth=-3.82668pt}}}{{\hbox{$\scriptscriptstyle\sqrt{N_{\mathrm{end}}-N_{\mathrm{p}}\,}$}\lower
0.4pt\hbox{\vrule
height=3.41666pt,depth=-2.73334pt}}}-\frac{3}{2}m_{\ell}^{2}\left(N_{\mathrm{end}}-N_{\mathrm{p}}\right).$
(A.10)

This approximation is represented and compared to the exact
$\omega^{2}_{\mathrm{h}}(t)$ in Fig. 20 (solid orange line).

###  A.2 The WKB approximation for the heavy field

Figure 21: Time evolution of the quantity $Q/\omega^{2}_{\mathrm{h}}$ (blue
line) with $Q$ defined in Eq. (A.11), before, around and after the time
“$\mathrm{osc}$”. The smallness of this quantity represents a criterion for
the regime of validity of the WKB approximation, whose strictness is
parameterized by the number $\mathcal{C}$ which is either $1$ (horizontal
black line) or $0.1$ (horizontal gray line) in this work. The corresponding
time from which the WKB approximation becomes valid are represented by
vertical dashed lines (respectively black and gray ones). The time from which
the quadratic approximation for the pulsation $\omega_{\mathrm{h}}^{2}$
becomes roughly valid is also represented by the orange vertical dashed lines.

As discussed earlier, Eq. (A.2) cannot be solved exactly due to the complex
analytical time dependence of $\omega_{\mathrm{h}}^{2}(t)$. However, the WKB
approximation can be used to find an approximate solution. Indeed, if one
represents the WKB criterion, $Q/\omega^{2}_{\mathrm{h}}$, where

$\displaystyle
Q=\frac{3}{4}\frac{\dot{\omega}^{2}_{\mathrm{h}}}{\omega^{2}_{\mathrm{h}}}-\frac{\ddot{\omega_{\mathrm{h}}}}{2\omega_{\mathrm{h}}},$
(A.11)

versus time, see Fig. 21, we notice that, quite quickly, on the left and right
hand sides of the turning point (where, by definition, $Q$ blows up), this
quantity tends to zero, signaling that WKB should be a very good
approximation. Usually, the WKB approximation is used in the following way. It
is first used on the left side of the turning point (region I) where the
initial conditions are specified and, as consequence, where the solution is
completely specified. The WKB approximation is also used on the right side
(region III) where the solution, being a solution of a second order
differential equation, depends on two arbitrary integration constants. These
two constants can be fixed if the solution in region III is matched to the
solution in region I. However, as noticed above, the solution in between
(region II) is a priori not known since the WKB approximation fails in this
region. What is usually done to deal with this problem is to linearize the
effective frequency in region II which allows us to find an approximation of
the solution in this regime which, then, permits to connect the solution in
region I to the solution in region III. In this way, one can specified the two
unknown integration constants in region III and have a complete knowledge of
the solution.

Here, we closely follow this WKB procedure with, however, some differences in
order to adapt our method to the problem at hand. In region I, the solution is
known and is given by the slow-roll approximation. So we do not need to employ
the WKB approximation in this regime (although it might be interesting to
study how WKB is connected to slow-roll in this regime, see also Ref. [95]).
As a consequence, we will “start the evolution” not deep in region I (on the
left) but at the onset of region II. In fact, the time of transition,
$t_{\mathrm{trans}}$ turns out to be located precisely at the beginning of
this region and, therefore, we fix the “initial conditions” (those relevant
for our problem here; they have obviously nothing to do with the physical
initial conditions, discussed before, and chosen at the onset of inflation) at
this time. In addition, they are completely known from the previously studied
slow-roll solution and read

$g_{\mathrm{h}}(t_{\mathrm{trans}})=1\,,\quad\dot{g}_{\mathrm{h}}(t_{\mathrm{trans}})=-\frac{m_{\mathrm{h}}^{2}}{3H(t_{\mathrm{trans}})}+\frac{3H(t_{\mathrm{trans}})}{2}+\frac{\Gamma_{\mathrm{hm}}}{4}\,,$
(A.12)

where we used the slow-roll expression of $\dot{\phi}_{\mathrm{h}}$.
Numerically, with our fiducial parameters, given that
$H_{\mathrm{trans}}=4.3\times 10^{-5}$ (obtained analytically within the slow-
roll approximation), we find that $\dot{g}_{\mathrm{trans}}=4.8\times
10^{-5}$. Moreover, one can compute the difference in cosmic time between the
transition [$\theta_{\mathrm{trans}}=\mathrm{arctan}(1/R)$] and the onset of
the oscillations
[$\theta_{\mathrm{osc}}=\mathrm{arctan}(\mathchoice{{\hbox{$\displaystyle\sqrt{0.0104\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{0.0104\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{0.0104\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{0.0104\,}$}\lower
0.4pt\hbox{\vrule height=3.22221pt,depth=-2.57779pt}}})$], with use of Eq.
(5.2.2) which is indeed relevant around the transition time, and one finds
$m_{\ell}(t_{\mathrm{osc}}-t_{\mathrm{trans}})=0.32$.

Let us now turn to the question of calculating $g_{\mathrm{h}}(t)$ in region
II. As explained before, in this region, the WKB approximation breaks down and
we need to use another technique which consists in linearizing the effective
frequency around the turning point. Concretely, we write

$\displaystyle\omega^{2}_{\mathrm{h}}=-\alpha(t-t_{\mathrm{osc}})+\cdots$
(A.13)

with
$\alpha\equiv-\mathrm{d}\omega^{2}_{\mathrm{h}}(t=t_{\mathrm{osc}})/\mathrm{d}t$.
In our case, $\alpha$ is negative. The linear approximation (A.13) is
represented in Fig. 20, see the solid green line. We have seen before that, in
the vicinity of the turning point, $\omega^{2}_{\mathrm{h}}\simeq
m_{\mathrm{h}}^{2}-9H^{2}/4$, and this implies that
$\alpha=-9/2\times(\epsilon_{1}H^{3})_{\mathrm{osc}}$. With the fiducial
parameters, we find analytically $\alpha=-3.0\times 10^{-14}$ in
$M_{\scriptscriptstyle{\mathrm{Pl}}}^{3}$ units, where one has used the slow-
roll formula for $\epsilon_{1}$, namely

$\epsilon_{1}=\frac{1}{2s}\frac{1+R^{4}\mathrm{tan}^{2}\theta}{\cos^{2}\theta\left(1+R^{2}\tan^{2}\theta\right)^{2}}\,.$
(A.14)

Then, plugging the expression (A.13) into the equation (A.2), and introducing
a new variable $y$ defined by the expression
$y\equiv\alpha^{1/3}(t-t_{\mathrm{osc}})$, the equation for
$g_{\mathrm{h}}(t)$ takes the form

$\displaystyle\frac{\mathrm{d}^{2}g_{\mathrm{h}}}{\mathrm{d}y^{2}}-yg_{\mathrm{h}}=0.$
(A.15)

The solution to this equation is well-known and can be expressed in terms of
the Airy functions of first and second kind, $\mathrm{Ai}$ and $\mathrm{Bi}$
[92, 93]. As a consequence, one can write

$\displaystyle g_{\mathrm{h}}(t)=B_{1}\mathrm{Ai}(y)+B_{2}\mathrm{Bi}(y),$
(A.16)

where $B_{1}$ and $B_{2}$ are two integration constants. These constants can
be determined using the values of $g_{\mathrm{h}}$ and its derivative at the
transition point, see Eqs. (A.12). This leads to

$\displaystyle B_{1}$
$\displaystyle=\pi\left[\frac{\mathrm{d}}{\mathrm{d}s}\mathrm{Bi}(y_{\mathrm{trans}})-\frac{1}{\alpha^{1/3}}\dot{g}_{\mathrm{trans}}\mathrm{Bi}(y_{\mathrm{trans}})\right],$
(A.17) $\displaystyle B_{2}$
$\displaystyle=-\pi\left[\frac{\mathrm{d}}{\mathrm{d}s}\mathrm{Ai}(y_{\mathrm{trans}})-\frac{1}{\alpha^{1/3}}\dot{g}_{\mathrm{trans}}\mathrm{Ai}(y_{\mathrm{trans}})\right].$
(A.18)

Numerically, with our fiducial parameters, one finds, $B_{1}\simeq 1.98+3.90i$
and $B_{2}\simeq 3.75-2.26i$. An intermediate step has been to use
$y_{\mathrm{trans}}=e^{i\pi/3}|\alpha|^{1/3}(t_{\mathrm{trans}}-t_{\mathrm{osc}})$,
with the value of $(t_{\mathrm{trans}}-t_{\mathrm{osc}})$ that was found in
the previous section.555These results and the following ones (until the
numerical determination of the constants $C_{\pm}$) are derived with the
choice of a complex cubic root for the negative number
$\alpha=|\alpha|e^{i\pi}$ which is $\alpha^{1/3}=e^{i\pi/3}|\alpha|^{1/3}$.
Note that in this case the time coordinate $y$ has a constant phase $\pi/3$,
as indeed $y=e^{i\pi/3}|\alpha|^{1/3}(t-t_{\mathrm{osc}})$, a fact that is
crucial in order to derive consistently the asymptotic limit shown in Eq.
(A.2). This choice should not matter as long as it is consistent during the
whole calculation, and indeed we have checked that using the different cubic
root $\alpha^{1/3}=-|\alpha|^{1/3}$, leads to the same final results (but
intermediate ones like the values for $y_{\mathrm{trans}}$ or $B_{1,2}$ are
therefore obviously different). This completes our calculation of
$g_{\mathrm{h}}$ (and, therefore, of $\phi_{\mathrm{h}}$) in region II which
is now completely known.

The final step consists in writing the solution in region III. This region is
the most interesting one for us since it corresponds to the phase dominated by
the light field. Obtaining the solution for $g_{\mathrm{h}}(t)$ in this regime
allows us to track the behavior of $\phi_{\mathrm{h}}(t)$ when the light field
drives inflation and the heavy field oscillates and decays. As explained
before, in this region, $Q/\omega_{\mathrm{h}}^{2}\ll 1$ and the WKB
approximation is satisfied. As a consequence, the solution can be written as

$\displaystyle
g_{\mathrm{h}}(t)=\frac{C_{+}}{\omega_{\mathrm{h}}^{1/2}}e^{i\int_{t_{\mathrm{osc}}}^{t}\omega_{\mathrm{h}}(\tau)\mathrm{d}\tau}+\frac{C_{-}}{\omega_{\mathrm{h}}^{1/2}}e^{-i\int_{t_{\mathrm{osc}}}^{t}\omega_{\mathrm{h}}(\tau)\mathrm{d}\tau},$
(A.19)

where $C_{+}$ and $C_{-}$ are two integration constants. This solution is
valid only in region III, namely only for times $t$ where
$Q/\omega_{\mathrm{h}}^{2}\ll 1$ (on the right hand side of the turning
point). Notice, however, that the lower bounds in the integrals have been
taken at the turning point. We first remark that changing the lower bound is
in fact equivalent to redefining the constants $C_{+}$ and $C_{-}$. So
choosing the lower bound to be the turning point is not in contradiction with
the WKB approximation and just amounts to a particular (convenient) definition
of $C_{+}$ and $C_{-}$. At this stage, our first goal is to determine these
two quantities. In order to carry out this calculation, we must check that
there is an overlap between region II and region III, that is to say a
patching region where the linear approximation of the effective frequency is
not too bad and where, at the same time, the WKB regime is satisfied. The
linear approximation is acceptable if the second order term remains smaller
than the first order one which means $t-t_{\mathrm{osc}}\lesssim
2(\mathrm{d}\omega_{\mathrm{h}}^{2}/\mathrm{d}t)/(\mathrm{d}^{2}\omega^{2}_{\mathrm{h}}/\mathrm{d}t^{2})$
or

$\displaystyle
m_{\ell}(t-t_{\mathrm{osc}})\lesssim\frac{3}{R(\epsilon_{1}-3\epsilon_{2})|_{\mathrm{osc}}}.$
(A.20)

On the other hand, if the linear approximation is accurate, the WKB
approximation is satisfied provided
$|\dot{\omega}_{\mathrm{h}}/\omega_{\mathrm{h}}^{2}|<1$, namely

$\displaystyle
m_{\ell}(t-t_{\mathrm{osc}})\gtrsim\left(\frac{m_{\ell}^{3}}{4|\alpha|}\right)^{1/3}.$
(A.21)

This is possible if

$\displaystyle|\alpha|>\frac{m_{\mathrm{h}}^{3}}{108}(\epsilon_{1}-3\epsilon_{2})^{3}|_{\mathrm{osc}},$
(A.22)

which is satisfied since, numerically, one finds $|\alpha|>2.5\times 10^{-14}$
to be compared with the value found here, $|\alpha|=3.0\times 10^{-14}$. In
this patching region, one can use at the same time the solution (A.16) and the
WKB solution (A.19). Since both formula must lead to the same solution, this
allows us to identify the coefficients $C_{+}$ and $C_{-}$ in terms of $B_{1}$
and $B_{2}$. Concretely, if one goes away from the turning point, Eq. (A.16)
leads to

$\displaystyle g_{\mathrm{h}}(t)$
$\displaystyle\simeq\frac{1}{2\mathchoice{{\hbox{$\displaystyle\sqrt{\pi\,}$}\lower
0.4pt\hbox{\vrule
height=4.30554pt,depth=-3.44446pt}}}{{\hbox{$\textstyle\sqrt{\pi\,}$}\lower
0.4pt\hbox{\vrule
height=4.30554pt,depth=-3.44446pt}}}{{\hbox{$\scriptstyle\sqrt{\pi\,}$}\lower
0.4pt\hbox{\vrule
height=3.01389pt,depth=-2.41113pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\pi\,}$}\lower
0.4pt\hbox{\vrule
height=2.15277pt,depth=-1.72223pt}}}}e^{-i\pi/12}|\alpha|^{-1/12}\left(t-t_{\mathrm{osc}}\right)^{-1/4}\biggl{[}(B_{1}+iB_{2})\,e^{-\frac{2i}{3}|\alpha|^{1/2}(t-t_{\mathrm{osc}})^{3/2}}$
$\displaystyle+2B_{2}\,e^{\frac{2i}{3}|\alpha|^{1/2}(t-t_{\mathrm{osc}})^{3/2}}\biggr{]},$
(A.23)

while the WKB solution (A.19) in the linear part of the time-dependent
frequency can be written as

$g_{\mathrm{h}}(t)\simeq|\alpha|^{-1/4}(t-t_{\mathrm{osc}})^{-1/4}\left[C_{+}e^{\frac{2i}{3}|\alpha|^{1/2}(t-t_{\mathrm{osc}})^{3/2}}+e^{-\frac{2i}{3}|\alpha|^{1/2}(t-t_{\mathrm{osc}})^{3/2}}\right]\,.$
(A.24)

Requiring these two expressions to represent the same solution, this
immediately leads to the relationship

$\displaystyle C_{+}$
$\displaystyle=\frac{B_{2}}{\mathchoice{{\hbox{$\displaystyle\sqrt{\pi\,}$}\lower
0.4pt\hbox{\vrule
height=4.30554pt,depth=-3.44446pt}}}{{\hbox{$\textstyle\sqrt{\pi\,}$}\lower
0.4pt\hbox{\vrule
height=4.30554pt,depth=-3.44446pt}}}{{\hbox{$\scriptstyle\sqrt{\pi\,}$}\lower
0.4pt\hbox{\vrule
height=3.01389pt,depth=-2.41113pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\pi\,}$}\lower
0.4pt\hbox{\vrule
height=2.15277pt,depth=-1.72223pt}}}}|\alpha|^{1/6}e^{-i\pi/12}\,,\quad
C_{-}=\frac{B_{1}+iB_{2}}{2\mathchoice{{\hbox{$\displaystyle\sqrt{\pi\,}$}\lower
0.4pt\hbox{\vrule
height=4.30554pt,depth=-3.44446pt}}}{{\hbox{$\textstyle\sqrt{\pi\,}$}\lower
0.4pt\hbox{\vrule
height=4.30554pt,depth=-3.44446pt}}}{{\hbox{$\scriptstyle\sqrt{\pi\,}$}\lower
0.4pt\hbox{\vrule
height=3.01389pt,depth=-2.41113pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\pi\,}$}\lower
0.4pt\hbox{\vrule
height=2.15277pt,depth=-1.72223pt}}}}|\alpha|^{1/6}e^{-i\pi/12}\,,$ (A.25)

Numerically, one finds $C_{+}\simeq 0.0096-0.01i$ and $C_{-}=C_{+}^{*}$, as
required since $g_{\mathrm{h}}(t)$ is a real function. Therefore, we have
reached our goal, namely expressing $g_{\mathrm{h}}(t)$ in the region where
inflation is driven by the light field. The solution is a priori entirely
specified since the constants $C_{+}$ and $C_{-}$ are now related to what
happens at the transition time. However, we must also calculate the WKB phase
which requires the integration of $\omega_{\mathrm{h}}(t)$. It can be
expressed as

$\displaystyle\int_{t_{\mathrm{osc}}}^{t}\omega(\tau)\mathrm{d}\tau=\int_{t_{\mathrm{osc}}}^{t_{\mathrm{match}}}\omega(\tau)\mathrm{d}\tau+\int_{t_{\mathrm{match}}}^{t}\omega(\tau)\mathrm{d}\tau=\Phi+\int_{t_{\mathrm{match}}}^{t}\omega(\tau)\mathrm{d}\tau,$
(A.26)

where we have introduced the time $t_{\mathrm{match}}$. This time is a time at
which the quadratic approximation of $\omega_{\mathrm{h}}^{2}(t)$, see Eq.
(A.7), becomes valid. In principle, the phase $\Phi$ can be calculated but, in
our case, this would require the integration of the exact time-dependent
frequency which is complicated given its form (A.3) (and which would undermine
the interest of having an analytical approximation). A rough approximation is
to use the linear expression in the vicinity of the turning point which
results in

$\displaystyle\Phi\simeq\frac{2}{3}|\alpha|^{1/2}(t_{\mathrm{match}}-t_{\mathrm{osc}})^{3/2}.$
(A.27)

However, this result does not give a very accurate estimation of the phase
while the accuracy of the WKB approximation quite sensitively depends on
$\Phi$. A more careful calculation, using $\omega_{\mathrm{h}}^{2}\simeq
m_{\mathrm{h}}^{2}-9H^{2}/4$ gives, after the change of variable defined in
Eq. (5.22)

$\Phi=-\mathchoice{{\hbox{$\displaystyle\sqrt{6C\,}$}\lower 0.4pt\hbox{\vrule
height=6.83331pt,depth=-5.46667pt}}}{{\hbox{$\textstyle\sqrt{6C\,}$}\lower
0.4pt\hbox{\vrule
height=6.83331pt,depth=-5.46667pt}}}{{\hbox{$\scriptstyle\sqrt{6C\,}$}\lower
0.4pt\hbox{\vrule
height=4.78333pt,depth=-3.82668pt}}}{{\hbox{$\scriptscriptstyle\sqrt{6C\,}$}\lower
0.4pt\hbox{\vrule
height=3.41666pt,depth=-2.73334pt}}}\frac{R}{R^{2}-1}\int_{\theta_{\mathrm{osc}}}^{\theta_{\mathrm{match}}}\frac{\mathchoice{{\hbox{$\displaystyle\sqrt{1+R^{2}\tan^{2}\theta\,}$}\lower
0.4pt\hbox{\vrule
height=8.63776pt,depth=-6.91023pt}}}{{\hbox{$\textstyle\sqrt{1+R^{2}\tan^{2}\theta\,}$}\lower
0.4pt\hbox{\vrule
height=8.63776pt,depth=-6.91023pt}}}{{\hbox{$\scriptstyle\sqrt{1+R^{2}\tan^{2}\theta\,}$}\lower
0.4pt\hbox{\vrule
height=6.0722pt,depth=-4.85779pt}}}{{\hbox{$\scriptscriptstyle\sqrt{1+R^{2}\tan^{2}\theta\,}$}\lower
0.4pt\hbox{\vrule
height=4.70554pt,depth=-3.76445pt}}}}{\sin\theta\times\cos\theta}(\tan\theta)^{\frac{1}{R^{2}-1}}\left[1-\frac{3C}{2R^{2}}\frac{\left(\mathrm{tan}\theta\right)^{\frac{2}{R^{2}-1}}}{\mathrm{cos}^{2}\theta}\right]^{1/2}\mathrm{d}\theta.$
(A.28)

Unfortunately, this integral cannot be computed analytically, but it is
straightforward to compute it numerically. It gives $\Phi=1.093$ for the
fiducial parameters of the model, a value that should be compared with the
numerical result $\Phi=2.039$. Admittedly, the result is still not very
accurate and, therefore, the determination of $\Phi$ constitutes an additional
source of error. Nevertheless, taking into account the phase determined
according to the considerations presented above significantly improves the
agreement between the exact and WKB solutions.

Let us now turn to the calculation of the second term in Eq. (A.26). As
already discussed above, this can be carried out with the help of the
quadratic approximation for $\omega_{\mathrm{h}}^{2}(t)$, see Eq. (A.7). It is
especially convenient to redefine the coefficients $\bar{a}$, $\bar{b}$ and
$\bar{c}$ according to $a=-\bar{a}/m_{\ell}=1/4$, $b=\bar{b}/m_{\ell}^{2}$ and
$c=\bar{c}/m_{\ell}^{2}$. In fact, since $c=R^{2}-1/2-b^{2}$, there remains
only one free coefficient which can be chosen to be $b$, a positive definite
quantity. Then, defining

$\displaystyle
u=\frac{m_{\ell}(t-t_{\mathrm{trans}})-2b}{2\mathchoice{{\hbox{$\displaystyle\sqrt{R^{2}-1/2\,}$}\lower
0.4pt\hbox{\vrule
height=8.63776pt,depth=-6.91023pt}}}{{\hbox{$\textstyle\sqrt{R^{2}-1/2\,}$}\lower
0.4pt\hbox{\vrule
height=8.63776pt,depth=-6.91023pt}}}{{\hbox{$\scriptstyle\sqrt{R^{2}-1/2\,}$}\lower
0.4pt\hbox{\vrule
height=6.0722pt,depth=-4.85779pt}}}{{\hbox{$\scriptscriptstyle\sqrt{R^{2}-1/2\,}$}\lower
0.4pt\hbox{\vrule height=4.70554pt,depth=-3.76445pt}}}}\,,$ (A.29)

the frequency $\omega_{\mathrm{h}}^{2}(t)$ takes the form

$\displaystyle\omega_{\mathrm{h}}(t)=m_{\ell}\left(R^{2}-\frac{1}{2}\right)^{1/2}\left(1-u^{2}\right)^{1/2}.$
(A.30)

Notice that, in the regime we are interested in, $u^{2}$ is always smaller
than one since $\omega_{\mathrm{h}}^{2}(t)$ is always real and positive. Then,
it follows that

$\displaystyle\int_{t_{\mathrm{match}}}^{t}\omega(\tau)\mathrm{d}\tau=2\left(R^{2}-\frac{1}{2}\right)\int_{u_{\mathrm{match}}}^{u}\mathrm{d}u\mathchoice{{\hbox{$\displaystyle\sqrt{1-u^{2}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{1-u^{2}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{1-u^{2}\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{1-u^{2}\,}$}\lower
0.4pt\hbox{\vrule
height=3.44165pt,depth=-2.75334pt}}}=\left(R^{2}-\frac{1}{2}\right)G(u),$
(A.31)

with

$\displaystyle
G(u)=\left[u\mathchoice{{\hbox{$\displaystyle\sqrt{1-u^{2}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{1-u^{2}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{1-u^{2}\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{1-u^{2}\,}$}\lower
0.4pt\hbox{\vrule
height=3.44165pt,depth=-2.75334pt}}}+\mathrm{arcsin}(u)\right]^{u}_{u_{\mathrm{match}}}.$
(A.32)

This completes the calculation of $g_{\mathrm{h}}(t)$ since all the terms of
Eq. (A.19) are now explicitly known.666It is also interesting to notice that,
in region III, if the quadratic approximation for the effective frequency is
used, then an exact solution is available. Let us indeed define
$z=2^{1/2}a^{1/4}[x-b/(2a)]$ and
$A=-1/(2\mathchoice{{\hbox{$\displaystyle\sqrt{a\,}$}\lower 0.4pt\hbox{\vrule
height=4.30554pt,depth=-3.44446pt}}}{{\hbox{$\textstyle\sqrt{a\,}$}\lower
0.4pt\hbox{\vrule
height=4.30554pt,depth=-3.44446pt}}}{{\hbox{$\scriptstyle\sqrt{a\,}$}\lower
0.4pt\hbox{\vrule
height=3.01389pt,depth=-2.41113pt}}}{{\hbox{$\scriptscriptstyle\sqrt{a\,}$}\lower
0.4pt\hbox{\vrule
height=2.15277pt,depth=-1.72223pt}}})[b^{2}/(4a)+c]=1/2-R^{2}$. Concretely,
$z$ can be expressed as  $\displaystyle z$
$\displaystyle=-\frac{1}{2}\frac{\Gamma_{\mathrm{hm}}}{m_{\ell}}-\mathchoice{{\hbox{$\displaystyle\sqrt{6\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{6\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{6\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{6\,}$}\lower
0.4pt\hbox{\vrule
height=3.22221pt,depth=-2.57779pt}}}(N_{\mathrm{end}}-N)^{1/2}.$ (A.33) Then,
it follows that Eq. (A.2) can be written as
$\displaystyle\frac{\mathrm{d}^{2}g}{\mathrm{d}z^{2}}-\left(\frac{z^{2}}{4}+A\right)g=0,$
(A.34) which is the canonical form of the equation defining the Parabolic
cylinder functions [92, 93]. As a consequence, the solution in region III
takes the following form [92, 93] $\displaystyle
g(z)=C_{1}\mathrm{U}(A,z)+C_{2}\mathrm{V}(A,z),$ (A.35) with
$\displaystyle\mathrm{U}(A,z)$ $\displaystyle=\mathrm{D}_{-A-1/2}(z),$ (A.36)
$\displaystyle\mathrm{V}(A,z)$
$\displaystyle=\frac{1}{\pi}\Gamma\left(\frac{1}{2}+A\right)\left[\mathrm{D}_{-A-1/2}(-z)+\sin(\pi
A)\mathrm{D}_{-A-1/2}(z)\right]$ (A.37) where $C_{1}$ and $C_{2}$ are two
integration constant and $\mathrm{D}_{\nu}(z)$ are parabolic functions. Notice
that another canonical form which is also convenient is [92, 93]
$\displaystyle\frac{\mathrm{d}^{2}g}{\mathrm{d}z^{2}}+\left(\nu+\frac{1}{2}-\frac{z^{2}}{4}\right)g=0,$
(A.38) where $A=-\nu-1/2$ or $\nu=R^{2}-1$. In that case, the solution can be
written as  $\displaystyle
g(z)=C_{1}\mathrm{D}_{\nu}(z)+C_{2}\mathrm{D}_{\nu}(-z),$ (A.39) since
$\mathrm{D}_{\nu}(\pm z)$ are linearly independent provided $\nu\neq 0,\pm
1,\cdots$ (we have used the same notations for the integration constants but
obviously they differ from the ones already introduced before).

Figure 22: Comparison of $\rho_{\mathrm{h}}$ and $\rho_{\mathrm{m}}$ (in
$M_{\scriptscriptstyle{\mathrm{Pl}}}^{4}$ units) during the oscillations and
the decay of the heavy field, between their exact numerical results
(respectively the orange and green solid lines), and their analytical
approximations inferred from the WKB approach (respectively the orange and
green dashed lines), see Eq. (A.40) and Eq. (A.3). The time of decay of the
heavy field, defined as the time when $\rho_{\mathrm{h}}=\rho_{\mathrm{m}}$ is
shown with the vertical gray lines both from the numerical approach (solid
line, $N_{\mathrm{decay,h}}=-10.81$) and the analytical one (dashed,
$N_{\mathrm{decay,h}}^{{}_{\mathrm{WKB}}}=-10.35$). This approach clearly
gives much better results than the one presented in body of the paper,
compared to Fig. 3 and the previous estimate
$N_{\mathrm{decay,h}}^{\mathrm{sd}}=-8.39$ , and enables one to derive
reliably and analytically the time of decay $N_{\mathrm{decay,h}}$.

The result of the WKB approximation is represented in Fig. 22 where the
evolution of the heavy field energy density

$\displaystyle\rho_{\mathrm{h}}(t)=\frac{1}{2}\dot{\phi}_{\mathrm{h}}^{2}+\frac{1}{2}m_{\mathrm{h}}^{2}\phi_{\mathrm{h}}^{2}$
(A.40)

is represented and compared to the numerical solution. Evidently, the result
is excellent and the WKB solution follows quite accurately the exact result.

###  A.3 The evolution of the matter energy density

Having calculated the evolution of the heavy field, we now turn to the
behavior of the matter energy density $\rho_{\mathrm{m}}$ whose equation of
motion is Eq. (5.8). As it was already noticed before, its solution can be
expressed as

$\displaystyle
a^{3}(t)\rho_{\mathrm{m}}(t)=a^{3}(t_{\mathrm{ini}})\rho_{\mathrm{m}}(t_{\mathrm{ini}})+\frac{\Gamma_{\mathrm{hm}}}{2}\int_{t_{\mathrm{ini}}}^{t}a^{3}(\tau)\dot{\phi}_{\mathrm{h}}^{2}(\tau)\mathrm{d}\tau.$
(A.41)

We see that this requires the calculation of the integral of the kinetic
energy density of the field, ${\cal
I}=\int_{t_{\mathrm{ini}}}^{t}a^{3}(\tau)\dot{\phi}_{\mathrm{h}}^{2}(\tau)\mathrm{d}\tau$.
In the slow-roll regime, $\dot{\phi}_{\mathrm{h}}^{2}(t)$ evolves slowly and
could be put outside the integral, immediately leading to an explicit
solution, see Sec. 5.2.3. Here, however, it is no longer possible to use this
trick and we must insert the WKB expression of $\phi_{\mathrm{h}}(t)$ in the
above integral and calculate it explicitly. This leads to

$\displaystyle{\cal I}=$ $\displaystyle
a_{\mathrm{p}}^{3}\phi_{\mathrm{p}}^{2}\int_{t_{\mathrm{match}}}^{t}\mathrm{d}\tau
e^{-\Gamma_{\mathrm{hm}}(\tau-
t_{\mathrm{p}})/2}\biggl{(}\frac{9}{4}H^{2}g^{2}_{\mathrm{h}}+\frac{3}{4}H\Gamma_{\mathrm{hm}}g^{2}_{\mathrm{h}}-3Hg_{\mathrm{h}}\dot{g}_{\mathrm{h}}-\frac{\Gamma_{\mathrm{hm}}}{2}g_{\mathrm{h}}\dot{g}_{\mathrm{h}}$
$\displaystyle+\frac{\Gamma_{\mathrm{hm}}^{2}}{16}g^{2}_{\mathrm{h}}+\dot{g}^{2}_{\mathrm{h}}\biggr{)}.$
(A.42)

In the present case, the initial time has been taken to be
$t_{\mathrm{match}}$, the time from which the quadratic assumption becomes
satisfied. Using the WKB solution (A.19) for $g_{\mathrm{h}}(t)$, the above
integral takes the form

$\displaystyle{\cal I}$
$\displaystyle=\phi_{\mathrm{p}}a_{\mathrm{p}}^{3}\int_{t_{\mathrm{match}}}^{t}\mathrm{d}\tau
e^{-\Gamma_{\mathrm{hm}}(\tau-
t_{\mathrm{p}})/2}\Biggl{[}C_{+}^{2}\left(\frac{{\cal
C}}{\omega_{\mathrm{h}}}-i\right)^{2}\omega_{\mathrm{h}}\,e^{2i\int_{t_{\mathrm{osc}}}^{\tau}\omega_{\mathrm{h}}\mathrm{d}\tau}$
$\displaystyle+C_{-}^{2}\left(\frac{{\cal
C}}{\omega_{\mathrm{h}}}+i\right)^{2}\omega_{\mathrm{h}}\,e^{-2i\int_{t_{\mathrm{osc}}}^{\tau}\omega_{\mathrm{h}}\mathrm{d}\tau}+2C_{+}C_{-}\omega_{\mathrm{h}}\left(\frac{{\cal
C}}{\omega_{\mathrm{h}}}-i\right)\left(\frac{{\cal
C}}{\omega_{\mathrm{h}}}+i\right)\Biggr{]},$ (A.43)

with

$\displaystyle\frac{{\cal
C}}{\omega_{\mathrm{h}}}=\frac{3}{2}\frac{H}{\omega_{\mathrm{h}}}+\frac{1}{4}\frac{\Gamma_{\mathrm{hm}}}{\omega_{\mathrm{h}}}+\frac{\dot{\omega}_{\mathrm{h}}}{2\omega^{2}_{\mathrm{h}}}\simeq\frac{3}{2}\frac{H}{\omega_{\mathrm{h}}}+\frac{1}{4}\frac{\Gamma_{\mathrm{hm}}}{\omega_{\mathrm{h}}}.$
(A.44)

We see that we have three terms to calculate and we write them as ${\cal
I}={\cal I}_{1}+{\cal I}_{2}+{\cal I}_{3}$. The next step consists in
inserting the form of $\omega_{\mathrm{h}}(t)$ in the above expression. Since
we are interested in the behavior of $\rho_{\mathrm{m}}$ in region III, one
can use the quadratic approximation of $\omega^{2}_{\mathrm{h}}(t)$, see Eq.
(A.30) expressed in terms of the variable $u$, see Eq. (A.29). Then, the first
term reads

$\displaystyle{\cal
I}_{1}=a_{\mathrm{p}}^{3}\phi_{\mathrm{p}}C_{+}^{2}e^{i2\Phi}(2R^{2}-1)e^{-b\Gamma_{\mathrm{hm}}/m_{\ell}}\int_{u_{\mathrm{match}}}^{u}\mathrm{d}u\,F(u)\,e^{i(2R^{2}-1)G(u)},$
(A.45)

with

$\displaystyle
F(u)=e^{-\Gamma_{\mathrm{hm}}(R^{2}-1/2)^{1/2}u/m_{\ell}}(1-u^{2})^{-1/2}\left[u-i\left(1-u^{2}\right)^{1/2}\right]^{2},$
(A.46)

and where the calculation of the WKB phase (A.31) has been used, see Eq.
(A.32) for the definition of the function $G(u)$. The second term is then
trivially found since this is just the complex conjugate of the first one,
${\cal I}_{2}={\cal I}_{1}^{*}$. Unfortunately, the integral (A.45) cannot be
performed analytically. However, we notice that the integrand has the form of
a smooth function $F(u)$ multiplied by an a priori rapid oscillatory phase
(since $2R^{2}-1$ is a priori quite a large number). As is well-known, several
techniques (Filon method, Levin collocation method etc., see Ref. [96]) have
been developed to deal with this type of integrals. The simplest one just
consists in integrating by part leading to

$\displaystyle\int_{u_{\mathrm{match}}}^{u}\mathrm{d}u\,F(u)\,e^{i(2R^{2}-1)G(u)}$
$\displaystyle=\frac{1}{i(2R^{2}-1)}\biggl{[}\frac{F(u)}{G^{\prime}(u)}e^{i(2R^{2}-1)G(u)}$
$\displaystyle-\frac{F(u_{\mathrm{match}})}{G^{\prime}(u_{\mathrm{match}})}e^{i(2R^{2}-1)G(u_{\mathrm{match}})}\biggr{]}$
$\displaystyle-\frac{1}{i(2R^{2}-1)}\int_{u_{\mathrm{match}}}^{u}\mathrm{d}u\frac{\mathrm{d}}{\mathrm{d}u}\left(\frac{F}{G^{\prime}}\right)e^{i(2R^{2}-1)G(u)}.$
(A.47)

In terms of the (inverse of the) large parameter $2R^{2}-1$, the second term
is of higher order and, therefore, can be neglected thus leading to an
approximation of the integral. If one uses this approximation to estimate
${\cal I}_{1}$, this gives

$\displaystyle{\cal I}_{1}$ $\displaystyle\simeq-
a_{\mathrm{p}}^{3}\phi_{\mathrm{p}}iC_{+}^{2}e^{i2\Phi}e^{-b\Gamma_{\mathrm{hm}}/m_{\ell}}\left[\frac{F(u)}{G^{\prime}(u)}e^{i(2R^{2}-1)G(u)}-\frac{F(u_{\mathrm{match}})}{G^{\prime}(u_{\mathrm{match}})}e^{i(2R^{2}-1)G(u_{\mathrm{match}})}\right].$
(A.48)

In principle, additional (successive) integrations by parts could lead to an
even more accurate approximation but, here, we will restrict ourselves to
first order.

Finally, the third term remains to be calculated. Straightforward
manipulations show that it can be written as

$\displaystyle{\cal I}_{3}$
$\displaystyle=2a_{\mathrm{p}}^{3}\phi_{\mathrm{p}}C_{+}C_{-}(2R^{2}-1)e^{-b\Gamma_{\mathrm{hm}}/m_{\ell}}\int_{u_{\mathrm{match}}}^{u}\mathrm{d}u\,e^{-\aleph
u}\,(1-u^{2})^{-1/2}.$ (A.49)

with

$\displaystyle\aleph=\frac{\Gamma_{\mathrm{hm}}}{m_{\ell}}\left(R^{2}-\frac{1}{2}\right)^{1/2}.$
(A.50)

The integral (A.49) ${\cal
J}=\int_{u_{\mathrm{match}}}^{u}\mathrm{d}u\,e^{-\aleph u}\,(1-u^{2})^{-1/2}$
cannot be performed exactly. But it is however relatively easy to derive an
accurate analytical approximation of ${\cal J}$. One can proceed as follows.
Let us first introduce the new variable $\psi$ such that $u=\cos\psi$. Then
the integral to be calculated reduces to

$\displaystyle{\cal
J}=-\int^{\psi}_{\psi_{\mathrm{match}}}e^{-\aleph\cos\Omega}\,\mathrm{d}\Omega.$
(A.51)

The integrand can easily be written as Fourier series, namely

$\displaystyle
e^{-\aleph\cos\Omega}=I_{0}(-\aleph)+2\sum_{n=1}^{+\infty}I_{n}(-\aleph)\cos(n\Omega),$
(A.52)

where $I_{n}$ is a modified Bessel function [92, 93]. Each term of the series
(A.52) can easily be integrated exactly and it follows that the integral
${\cal J}$ can be expressed as

$\displaystyle{\cal
J}=\left[-I_{0}(-\aleph)\psi-2\sum_{n=1}^{+\infty}\frac{1}{n}I_{n}(-\aleph)\sin(n\psi)\right]^{\psi}_{\psi_{\mathrm{match}}}$
(A.53)

In practice, in order to reach a good accuracy, it is sufficient to keep only
the first two or three terms in the above sum. We conclude that

$\displaystyle{\cal I}_{3}$
$\displaystyle=-2a_{\mathrm{p}}^{3}\phi_{\mathrm{p}}C_{+}C_{-}(2R^{2}-1)e^{-b\Gamma_{\mathrm{hm}}/m_{\ell}}\Biggl{[}I_{0}(-\aleph)\arccos
u$ $\displaystyle+2\sum_{n=1}^{+\infty}\frac{1}{n}I_{n}(-\aleph)\sin(n\arccos
u)\Biggr{]}^{u}_{u_{\mathrm{match}}}.$ (A.54)

The final step consists in collecting the three terms calculated before in the
final expression of the matter energy density. One obtains the rather lengthy,
but explicit, following equation [using the fact that
$G(u_{\mathrm{match}})=0$ by definition]

$\displaystyle\rho_{\mathrm{m}}(t)$
$\displaystyle\simeq\rho_{\mathrm{m}}^{\mathrm{match}}\left[\frac{a_{\mathrm{match}}}{a(t)}\right]^{3}+\phi_{\mathrm{p}}\left[\frac{a_{\mathrm{p}}}{a(t)}\right]^{3}\frac{\Gamma_{\mathrm{hm}}}{2}e^{-b\Gamma_{\mathrm{hm}}/m_{\ell}}$
$\displaystyle\times\Biggl{(}2\mathrm{Re}\biggl{\\{}-iC_{+}^{2}e^{i2\Phi}\biggl{[}\frac{F(u)}{G^{\prime}(u)}e^{i(2R^{2}-1)G(u)}-\frac{F(u_{\mathrm{match}})}{G^{\prime}(u_{\mathrm{match}})}\biggr{]}\biggr{\\}}$
$\displaystyle-2\left|C_{+}\right|^{2}(2R^{2}-1)\Biggl{[}I_{0}(-\aleph)\arccos
u+2\sum_{n=1}^{+\infty}\frac{1}{n}I_{n}(-\aleph)\sin(n\arccos
u)\Biggr{]}^{u}_{u_{\mathrm{match}}}\Biggr{)}.$ (A.55)

This expression is represented in Fig. 22 together with the exact (numerical)
solution of $\rho_{\mathrm{m}}$. Obviously, it matches quite well the
numerical result.

###  A.4 Determination of the time of decay

Having determined the WKB form of $\rho_{\mathrm{h}}(t)$, see Eq. (A.40) and
$\rho_{\mathrm{m}}(t)$, see Eq. (A.3), we can now calculate the time of decay.
Recall that it is defined by the condition
$\rho_{\mathrm{h}}(t_{\mathrm{decay,h}})=\rho_{\mathrm{m}}(t_{\mathrm{decay,h}})$.
Given the complexity of Eqs. (A.40) and (A.3), it is clear that it is hopeless
to solve
$\rho_{\mathrm{h}}(t_{\mathrm{decay,h}})=\rho_{\mathrm{m}}(t_{\mathrm{decay,h}})$
analytically. Here, we solve this transcendental equation numerically and the
solution is denoted $N_{\mathrm{decay,h}}^{{}_{\mathrm{WKB}}}$. For the
fiducial parameters, one finds
$N_{\mathrm{decay,h}}^{{}_{\mathrm{WKB}}}=-10.35$ to be compared with the
exact value $N_{\mathrm{decay,h}}=-10.81$. The quantity
$N_{\mathrm{decay,h}}^{{}_{\mathrm{WKB}}}$ is used in Table. 2 for the
analytical estimation of the amplitude of non-adiabatic perturbations after
double inflation.

## References

  * [1] A. A. Starobinsky, _A New Type of Isotropic Cosmological Models Without Singularity_ , _Phys. Lett._ B91 (1980) 99.
  * [2] A. H. Guth, _The Inflationary Universe: A Possible Solution to the Horizon and Flatness Problems_ , _Phys.Rev._ D23 (1981) 347.
  * [3] A. D. Linde, _A New Inflationary Universe Scenario: A Possible Solution of the Horizon, Flatness, Homogeneity, Isotropy and Primordial Monopole Problems_ , _Phys.Lett._ B108 (1982) 389.
  * [4] A. Albrecht and P. J. Steinhardt, _Cosmology for Grand Unified Theories with Radiatively Induced Symmetry Breaking_ , _Phys.Rev.Lett._ 48 (1982) 1220.
  * [5] A. D. Linde, _Chaotic Inflation_ , _Phys.Lett._ B129 (1983) 177.
  * [6] A. A. Starobinsky, _Spectrum of relict gravitational radiation and the early state of the universe_ , _JETP Lett._ 30 (1979) 682.
  * [7] V. F. Mukhanov and G. Chibisov, _Quantum Fluctuation and Nonsingular Universe._ , _JETP Lett._ 33 (1981) 532.
  * [8] S. Hawking, _The Development of Irregularities in a Single Bubble Inflationary Universe_ , _Phys.Lett._ B115 (1982) 295.
  * [9] A. H. Guth and S. Pi, _Fluctuations in the New Inflationary Universe_ , _Phys.Rev.Lett._ 49 (1982) 1110.
  * [10] J. M. Bardeen, P. J. Steinhardt and M. S. Turner, _Spontaneous Creation of Almost Scale - Free Density Perturbations in an Inflationary Universe_ , _Phys.Rev._ D28 (1983) 679.
  * [11] J. Martin, _The Observational Status of Cosmic Inflation after Planck_ , _Astrophys. Space Sci. Proc._ 45 (2016) 41 [1502.05733].
  * [12] J. Martin, C. Ringeval and V. Vennin, _Encyclopædia Inflationaris_ , _Phys. Dark Univ._ 5-6 (2014) 75 [1303.3787].
  * [13] J. Martin, C. Ringeval, R. Trotta and V. Vennin, _The Best Inflationary Models After Planck_ , _JCAP_ 03 (2014) 039 [1312.3529].
  * [14] D. Chowdhury, J. Martin, C. Ringeval and V. Vennin, _Assessing the scientific status of inflation after Planck_ , _Phys. Rev. D_ 100 (2019) 083537 [1902.03951].
  * [15] D. Polarski and A. A. Starobinsky, _Spectra of perturbations produced by double inflation with an intermediate matter dominated stage_ , _Nucl. Phys. B_ 385 (1992) 623.
  * [16] P. Peter, D. Polarski and A. A. Starobinsky, _Confrontation of double inflationary models with observations_ , _Phys. Rev. D_ 50 (1994) 4827 [astro-ph/9403037].
  * [17] D. Langlois, _Correlated adiabatic and isocurvature perturbations from double inflation_ , _Phys. Rev. D_ 59 (1999) 123512 [astro-ph/9906080].
  * [18] C. Gordon, D. Wands, B. A. Bassett and R. Maartens, _Adiabatic and entropy perturbations from inflation_ , _Phys. Rev. D_ 63 (2000) 023506 [astro-ph/0009131].
  * [19] S. Groot Nibbelink and B. J. W. van Tent, _Scalar perturbations during multiple field slow-roll inflation_ , _Class. Quant. Grav._ 19 (2002) 613 [hep-ph/0107272].
  * [20] S. Tsujikawa, D. Parkinson and B. A. Bassett, _Correlation - consistency cartography of the double inflation landscape_ , _Phys. Rev. D_ 67 (2003) 083516 [astro-ph/0210322].
  * [21] D. Wands, N. Bartolo, S. Matarrese and A. Riotto, _An Observational test of two-field inflation_ , _Phys.Rev._ D66 (2002) 043520 [astro-ph/0205253].
  * [22] D. Wands, _Multiple field inflation_ , _Lect. Notes Phys._ 738 (2008) 275 [astro-ph/0702187].
  * [23] K.-Y. Choi, J.-O. Gong and D. Jeong, _Evolution of the curvature perturbation during and after multi-field inflation_ , _JCAP_ 0902 (2009) 032 [0810.2299].
  * [24] A. Achucarro, J.-O. Gong, S. Hardeman, G. A. Palma and S. P. Patil, _Mass hierarchies and non-decoupling in multi-scalar field dynamics_ , _Phys. Rev. D_ 84 (2011) 043502 [1005.3848].
  * [25] D. I. Kaiser, _Nonminimal Couplings in the Early Universe: Multifield Models of Inflation and the Latest Observations_ , _Fundam. Theor. Phys._ 183 (2016) 41 [1511.09148].
  * [26] R. Schimmrigk, _Multifield Reheating after Modular $j$-Inflation_, _Phys. Lett. B_ 782 (2018) 193 [1712.09961].
  * [27] M. Braglia, D. K. Hazra, L. Sriramkumar and F. Finelli, _Generating primordial features at large scales in two field models of inflation_ , _JCAP_ 08 (2020) 025 [2004.00672].
  * [28] L. Amendola et al., _Cosmology and fundamental physics with the Euclid satellite_ , _Living Rev. Rel._ 21 (2018) 2 [1606.00180].
  * [29] O. Doré, J. Bock, P. Capak, R. de Putter, T. Eifler et al., _SPHEREx: An All-Sky Spectral Survey_ , 1412.4872.
  * [30] T. Matsumura et al., _Mission design of LiteBIRD_ , _J. Low Temp. Phys._ 176 (2014) 733 [1311.2847].
  * [31] L. E. Allen, S. Gupta and D. Wands, _Non-gaussian perturbations from multi-field inflation_ , _JCAP_ 01 (2006) 006 [astro-ph/0509719].
  * [32] I. Huston and A. J. Christopherson, _Isocurvature Perturbations and Reheating in Multi-Field Inflation_ , 1302.4298.
  * [33] C. T. Byrnes and D. Wands, _Curvature and isocurvature perturbations from two-field inflation in a slow-roll expansion_ , _Phys. Rev. D_ 74 (2006) 043529 [astro-ph/0605679].
  * [34] K. Schutz, E. I. Sfakianakis and D. I. Kaiser, _Multifield Inflation after Planck: Isocurvature Modes from Nonminimal Couplings_ , _Phys. Rev. D_ 89 (2014) 064044 [1310.8285].
  * [35] X. Chen and Y. Wang, _Large non-Gaussianities with Intermediate Shapes from Quasi-Single Field Inflation_ , _Phys. Rev. D_ 81 (2010) 063511 [0909.0496].
  * [36] X. Chen and Y. Wang, _Quasi-Single Field Inflation and Non-Gaussianities_ , _JCAP_ 04 (2010) 027 [0911.3380].
  * [37] D. Baumann and D. Green, _Signatures of Supersymmetry from the Early Universe_ , _Phys. Rev. D_ 85 (2012) 103520 [1109.0292].
  * [38] V. Assassi, D. Baumann and D. Green, _On Soft Limits of Inflationary Correlation Functions_ , _JCAP_ 11 (2012) 047 [1204.4207].
  * [39] T. Noumi, M. Yamaguchi and D. Yokoyama, _Effective field theory approach to quasi-single field inflation and effects of heavy fields_ , _JHEP_ 06 (2013) 051 [1211.1624].
  * [40] J.-O. Gong, S. Pi and M. Sasaki, _Equilateral non-Gaussianity from heavy fields_ , _JCAP_ 11 (2013) 043 [1306.3691].
  * [41] N. Arkani-Hamed and J. Maldacena, _Cosmological Collider Physics_ , 1503.08043.
  * [42] H. Lee, D. Baumann and G. L. Pimentel, _Non-Gaussianity as a Particle Detector_ , _JHEP_ 12 (2016) 040 [1607.03735].
  * [43] S. Garcia-Saenz, S. Renaux-Petel and J. Ronayne, _Primordial fluctuations and non-Gaussianities in sidetracked inflation_ , _JCAP_ 1807 (2018) 057 [1804.11279].
  * [44] S. Garcia-Saenz and S. Renaux-Petel, _Flattened non-Gaussianities from the effective field theory of inflation with imaginary speed of sound_ , _JCAP_ 1811 (2018) 005 [1805.12563].
  * [45] J. Fumagalli, S. Garcia-Saenz, L. Pinol, S. Renaux-Petel and J. Ronayne, _Hyper non-Gaussianities in inflation with strongly non-geodesic motion_ , 1902.03221.
  * [46] T. Bjorkmo, R. Z. Ferreira and M. D. Marsh, _Mild non-gaussianities under perturbative control from rapid-turn inflation models_ , _Journal of Cosmology and Astroparticle Physics_ 2019 (2019) 036–036.
  * [47] D.-G. Wang, _On the inflationary massive field with a curved field manifold_ , _JCAP_ 01 (2020) 046 [1911.04459].
  * [48] S. Garcia-Saenz, L. Pinol and S. Renaux-Petel, _Revisiting non-Gaussianity in multifield inflation with curved field space_ , 1907.10403.
  * [49] R. Z. Ferreira, _Non-Gaussianities in models of inflation with large and negative entropic masses_ , _JCAP_ 08 (2020) 034 [2003.13410].
  * [50] L. Pinol, _Multifield inflation beyond $N_{\mathrm{field}}=2$: non-Gaussianities and single-field effective theory_, 2011.05930.
  * [51] A. Albrecht, P. J. Steinhardt, M. S. Turner and F. Wilczek, _Reheating an Inflationary Universe_ , _Phys. Rev. Lett._ 48 (1982) 1437.
  * [52] M. S. Turner, _Coherent Scalar Field Oscillations in an Expanding Universe_ , _Phys. Rev. D_ 28 (1983) 1243.
  * [53] Y. Shtanov, J. H. Traschen and R. H. Brandenberger, _Universe reheating after inflation_ , _Phys. Rev._ D51 (1995) 5438 [hep-ph/9407247].
  * [54] B. A. Bassett, S. Tsujikawa and D. Wands, _Inflation dynamics and reheating_ , _Rev. Mod. Phys._ 78 (2006) 537 [astro-ph/0507632].
  * [55] R. Allahverdi, R. Brandenberger, F.-Y. Cyr-Racine and A. Mazumdar, _Reheating in Inflationary Cosmology: Theory and Applications_ , _Ann. Rev. Nucl. Part. Sci._ 60 (2010) 27 [1001.2600].
  * [56] M. A. Amin, M. P. Hertzberg, D. I. Kaiser and J. Karouby, _Nonperturbative Dynamics Of Reheating After Inflation: A Review_ , _Int. J. Mod. Phys._ D24 (2014) 1530003 [1410.3808].
  * [57] J. Martin and C. Ringeval, _First CMB Constraints on the Inflationary Reheating Temperature_ , _Phys. Rev. D_ 82 (2010) 023511 [1004.5525].
  * [58] J. Martin, C. Ringeval and V. Vennin, _Observing Inflationary Reheating_ , _Phys. Rev. Lett._ 114 (2015) 081303 [1410.7958].
  * [59] J. Martin, C. Ringeval and V. Vennin, _Information Gain on Reheating: the One Bit Milestone_ , _Phys. Rev. D_ 93 (2016) 103532 [1603.02606].
  * [60] B. A. Bassett, D. I. Kaiser and R. Maartens, _General relativistic preheating after inflation_ , _Phys. Lett. B_ 455 (1999) 84 [hep-ph/9808404].
  * [61] B. A. Bassett, F. Tamburini, D. I. Kaiser and R. Maartens, _Metric preheating and limitations of linearized gravity. 2._ , _Nucl. Phys. B_ 561 (1999) 188 [hep-ph/9901319].
  * [62] B. A. Bassett, C. Gordon, R. Maartens and D. I. Kaiser, _Restoring the sting to metric preheating_ , _Phys. Rev. D_ 61 (2000) 061302 [hep-ph/9909482].
  * [63] S. Tsujikawa and B. A. Bassett, _A New twist to preheating_ , _Phys. Rev. D_ 62 (2000) 043510 [hep-ph/0003068].
  * [64] S. Tsujikawa and B. A. Bassett, _When can preheating affect the CMB?_ , _Phys. Lett. B_ 536 (2002) 9 [astro-ph/0204031].
  * [65] D. H. Lyth and D. Wands, _Conserved cosmological perturbations_ , _Phys. Rev. D_ 68 (2003) 103515 [astro-ph/0306498].
  * [66] F. Finelli and R. H. Brandenberger, _Parametric amplification of gravitational fluctuations during reheating_ , _Phys. Rev. Lett._ 82 (1999) 1362 [hep-ph/9809490].
  * [67] K. Jedamzik, M. Lemoine and J. Martin, _Collapse of Small-Scale Density Perturbations during Preheating in Single Field Inflation_ , _JCAP_ 1009 (2010) 034 [1002.3039].
  * [68] J. Martin, T. Papanikolaou and V. Vennin, _Primordial black holes from the preheating instability_ , _JCAP_ 2001 (2020) 024 [1907.04236].
  * [69] J. Martin, T. Papanikolaou, L. Pinol and V. Vennin, _Metric preheating and radiative decay in single-field inflation_ , _JCAP_ 05 (2020) 003 [2002.01820].
  * [70] F. Finelli and R. H. Brandenberger, _Parametric amplification of metric fluctuations during reheating in two field models_ , _Phys. Rev. D_ 62 (2000) 083502 [hep-ph/0003172].
  * [71] M. P. DeCross, D. I. Kaiser, A. Prabhu, C. Prescod-Weinstein and E. I. Sfakianakis, _Preheating after Multifield Inflation with Nonminimal Couplings, I: Covariant Formalism and Attractor Behavior_ , _Phys. Rev. D_ 97 (2018) 023526 [1510.08553].
  * [72] M. P. DeCross, D. I. Kaiser, A. Prabhu, C. Prescod-Weinstein and E. I. Sfakianakis, _Preheating after multifield inflation with nonminimal couplings, II: Resonance Structure_ , _Phys. Rev. D_ 97 (2018) 023527 [1610.08868].
  * [73] M. P. DeCross, D. I. Kaiser, A. Prabhu, C. Prescod-Weinstein and E. I. Sfakianakis, _Preheating after multifield inflation with nonminimal couplings, III: Dynamical spacetime results_ , _Phys. Rev. D_ 97 (2018) 023528 [1610.08916].
  * [74] J. Jiang, Q. Liang, Y.-F. Cai, D. A. Easson and Y. Zhang, _Numerical study of inflationary preheating with arbitrary power-law potential and a realization of curvaton mechanism_ , _Astrophys. J._ 876 (2019) 136 [1812.08220].
  * [75] R. Nguyen, J. van de Vis, E. I. Sfakianakis, J. T. Giblin and D. I. Kaiser, _Nonlinear Dynamics of Preheating after Multifield Inflation with Nonminimal Couplings_ , _Phys. Rev. Lett._ 123 (2019) 171301 [1905.12562].
  * [76] J. van de Vis, R. Nguyen, E. I. Sfakianakis, J. T. Giblin and D. I. Kaiser, _Time scales for nonlinear processes in preheating after multifield inflation with nonminimal couplings_ , _Phys. Rev. D_ 102 (2020) 043528 [2005.00433].
  * [77] K. A. Malik, D. Wands and C. Ungarelli, _Large scale curvature and entropy perturbations for multiple interacting fluids_ , _Phys. Rev. D_ 67 (2003) 063516 [astro-ph/0211602].
  * [78] K. A. Malik and D. Wands, _Adiabatic and entropy perturbations with interacting fluids and fields_ , _JCAP_ 0502 (2005) 007 [astro-ph/0411703].
  * [79] A. Berera, _Warm inflation_ , _Phys. Rev. Lett._ 75 (1995) 3218 [astro-ph/9509049].
  * [80] J. Yokoyama and A. D. Linde, _Is warm inflation possible?_ , _Phys. Rev. D_ 60 (1999) 083509 [hep-ph/9809409].
  * [81] V. F. Mukhanov, H. Feldman and R. H. Brandenberger, _Theory of cosmological perturbations. Part 1. Classical perturbations. Part 2. Quantum theory of perturbations. Part 3. Extensions_ , _Phys. Rept._ 215 (1992) 203.
  * [82] J. M. Bardeen, _Gauge Invariant Cosmological Perturbations_ , _Phys. Rev._ D22 (1980) 1882.
  * [83] K. A. Malik and D. Wands, _Cosmological perturbations_ , _Phys. Rept._ 475 (2009) 1 [0809.4944].
  * [84] L. Visinelli, _Cosmological perturbations for an inflaton field coupled to radiation_ , _JCAP_ 1501 (2015) 005 [1410.1187].
  * [85] J. Martin, _Inflationary cosmological perturbations of quantum-mechanical origin_ , _Lect. Notes Phys._ 669 (2005) 199 [hep-th/0406011].
  * [86] J. Martin, _Inflationary perturbations: The Cosmological Schwinger effect_ , _Lect. Notes Phys._ 738 (2008) 193 [0704.3540].
  * [87] J. Martin and V. Vennin, _Quantum Discord of Cosmic Inflation: Can we Show that CMB Anisotropies are of Quantum-Mechanical Origin?_ , _Phys. Rev. D_ 93 (2016) 023505 [1510.04038].
  * [88] H. Kodama and M. Sasaki, _Cosmological Perturbation Theory_ , _Prog. Theor. Phys. Suppl._ 78 (1984) 1.
  * [89] S. Weinberg, _Cosmology_. Oxford Univ. Press, 2008.
  * [90] W. H. Press, S. a. Teukolsky, W. T. Vetterling and B. P. Flannery, _Numerical Recipes in Fortran 77: the Art of Scientific Computing. Second Edition_ , vol. 1. 1996.
  * [91] S. M. Leach, A. R. Liddle, J. Martin and D. J. Schwarz, _Cosmological parameter estimation and the inflationary cosmology_ , _Phys. Rev. D_ 66 (2002) 023515 [astro-ph/0202094].
  * [92] I. S. Gradshteyn and I. M. Ryzhik, _Table of Integrals, Series, and Products_. Academic Press, New York and London, 1965.
  * [93] M. Abramowitz and I. A. Stegun, _Handbook of mathematical functions with formulas, graphs, and mathematical tables_. National Bureau of Standards, Washington, US, ninth ed., 1970.
  * [94] P. González, G. A. Palma and N. Videla, _Covariant evolution of perturbations during reheating in two-field inflation_ , _JCAP_ 12 (2018) 001 [1805.10360].
  * [95] J. Martin and D. J. Schwarz, _WKB approximation for inflationary cosmological perturbations_ , _Phys. Rev. D_ 67 (2003) 083512 [astro-ph/0210090].
  * [96] S. Olver, _Numerical approximation of highly oscillatory integrals_ , Ph.D. thesis, University of Cambridge, 2008.