The Thermal Feedback Effects on the Temperature Evolution during Reheating

Many properties of the observing universe today can be understood as the result of the processes happening during the early stage of our universe, e.g. the Big Bang Nucleosynthesis Tanabashi:2018oca and the abundance of the Standard Model matter Davidson:2002qv ; Fukugita:1986hr ; Canetti:2014dka as well as Dark Matter Acharya:2008bk ; Watson:2009hw ; Acharya:2009zt . In the inflationary cosmology, reheating Kofman:1997yn ; Kofman:1994rk ; Shtanov:1994ce is a mechanism to set up the hot big bang initial conditions at the onset of the radiation-dominated era. After inflation Starobinsky:1980te ; Guth:1980zm ; Linde:1981mu , the scalar inflaton field decays and dissipates its energy into the produced particles via the interaction between them, and (re)heats the universe. The history–especially the thermal history–of the universe is strongly affected by the matter production dynamics after inflation, and the abundance of the constituents of the plasma is sensitive to the details of the thermal history Bodeker:2020ghk . Thus the knowledge of the temperature evolution of the early universe is crucial for cosmology. To understand in more detail about the thermal history of the early universe, we would like to investigate the temperature evolution during the reheating process.

The reheating process begins when the universe stops accelerating and ends at the moment when the energy density of the radiation is comparable to that of the expectation value $\varphi=\langle\phi\rangle$ of the scalar inflaton field $\phi$, $\rho_{\gamma}=\rho_{\varphi}$. It is well known that the particle production could be non-perturbative and in highly non-equilibrium with large couplings between the inflaton and other fields, and thus other mechanisms, such as the parametric resonance in preheating Greene:1997fu ; Felder:1998vq ; Podolsky:2005bw , need to be considered. In those models, the effective dissipation rate $\Gamma_{\varphi}$ that leads to the energy transfer from $\varphi$ to radiation may have a complicated time dependence, which makes it hard to study the evolution of the effective temperature $T$ as an effective parameter of the hot plasma.

In this work, we assume that the plasma reaches equilibrium fast between the two zero crossings of the inflaton at the minimum of its potential, so that the occupation numbers can be expressed in terms of an effective temperature $T$ Harigaya:2013vwa ; Davidson:2000er . If the reheating process starts from a cold universe with all the energy stored in the zero modes of the inflaton field, i.e. the initial value of $T$ and $\rho_{\gamma}$ are assumed to be zero, then after inflation the particle production process is dominated by the vacuum decay of the inflaton, which can be expressed by its temperature-independent decay rate $\Gamma_{0}$. As the temperature increases, the thermal feedback effects of the produced particles on $\Gamma_{\varphi}$, which are crucial and can lead to enhanced production of particles, should be considered. We will also discuss the situation where the reheating phase starts with non-zero initial temperature in Section IV. In general, $\Gamma_{\varphi}$ has a complicated dependence on $T$ as shown in detailed studies Drewes:2013iaa ; Drewes:2015eoa ; Mukaida:2012bz , as different processes such as decays and scatterings can contribute to the dissipative rate and the phase space is temperature dependent due to the thermal correction to the effective masses. In general cases, we need to carefully study the particle production, thermalization, and resulting dissipation effect on the inflaton coherent oscillation.

The reheating temperature $T_{R}$ as the initial temperature of the radiation-dominated era is important for cosmology and is usually not the highest temperature during the early universe Giudice:2000ex . During reheating, the inflaton oscillates around the minimum of its effective potential and loses a small fraction of its energy per oscillation due to the relative smallness of $\Gamma\dot{\varphi}$. However, the total amount of energy transferred from $\varphi$ to radiation is the largest at early times because $\rho_{\varphi}$ is huge in the beginning and red-shifted at later times. In common cases, the universe becomes radiation-dominated when $\Gamma_{\varphi}=H$, with $H\equiv\dot{a}/a$ being the Hubble parameter in terms of the scale factor of the universe $a$, and shortly afterwards $\rho_{\gamma}$ exceeds $\rho_{\varphi}$. Since temperature evolution strongly affects the abundance of thermal relics which are sensitive to the thermal history, a quantitative understanding of the reheating process is necessary.

We parameterize the effects as

$$\Gamma_{\varphi}=\sum_{n=0}^{\infty}\Gamma_{n}\left(\frac{T}{m_{\phi}}\right)^{n}~{},$$

(1)

where $m_{\phi}$ is the inflaton mass and $n$ a positive integer. We can always Taylor expand $\Gamma_{\varphi}$ to be $\sum_{n=0}^{\infty}\Gamma_{n}(T/m_{\phi})^{n}$ around a centain value of $T$. We study the set of scenarios where this can be done and the thermal part in $\Gamma_{\varphi}$ is dominated by one or few terms. In this way, it is possible to solve the Boltzmann equations for the energy densities $\rho_{\varphi}$ and $\rho_{\gamma}$ analytically, by considering the phases during which $\Gamma_{0}$ and $\Gamma_{n}$ dominate respectively and matching the solutions at the boundary. This piece-wise method approximates the complicated $\Gamma_{\varphi}$ to be a temperature-dependent monomial locally Drewes:2014pfa ; Drewes:2015coa and has been used to obtain improved analytic estimates of the reheating temperature and maximal temperature in the early universe in the $n=2$ case. In Co:2020xaf , an arbitrary temperature dependence as well as a dependence of the dissipation rate on the scale factor are considered, and the increasing temperature is observed when the rate depends on the scale factor in some physical scenarios.

Our present paper generalize the results of Drewes:2014pfa and extend the research to the cases where the thermal term is a monomial of arbitrary $n$ powers. The $n=1$ case can appear in the corrections from the induced $1\rightarrow 2$ decays. As will be seen, for a Bose enhanced $1\rightarrow 2$ decay one generally has the form

$$\Gamma_{\varphi}=\Gamma_{0}\left(1+2f_{B}\left(\frac{m_{\phi}}{2}\right)\right)~{},$$

(2)

then we get a linear term in $T$ by expanding the Bose-Einstein distribution function $f_{B}$ with $m_{\phi}<T$. The $n=2$ case appears quite generically in the high temperature regime when different processes such as $1\rightarrow 3$ decays or $2\rightarrow 2$ scatterings are included. For example in the $\alpha\phi\chi^{3}$ model one has Drewes:2013iaa

$$\Gamma_{\varphi}\approx\frac{\alpha^{2}m_{\phi}}{3072\pi^{3}}+\frac{\alpha^{2}T^{2}}{768\pi m_{\phi}}~{},$$

(3)

Where $\alpha$ is the coupling constant between the inflaton and another scalar field $\chi$. The higher power cases with $n>2$ can exist in some certain models Garcia:2020wiy ; Mukaida:2012bz . In general, the interactions between inflaton and other fields in those models are higher-order operators in a non-renormalizable theory by dimensional arguments. Yet an analysis with arbitrary positive or even negative $n$ is meaningful especially in intermediate temperature regimes for the realistic models in which $\Gamma_{\varphi}$ is a complicated function of $T$. For example the thermal correction with $\Gamma_{\varphi}\propto T^{n}$ where $n=3,~{}4$ or $5$ may solve the cosmological moduli problem Yokoyama:2006wt ; Bodeker:2006ij . Besides in curved spacetime with a finite temperature background, the effective dissipation rate might even contain negative $n$ powers functions of $T$ Ming:2019qty .

Our work generalizes the attempt to analytically solve the Boltzmann equations for the energy density of inflaton and radiation. The novel points of this paper are

1. 1.

   We describe the entire reheating process.

2. 2.

   We present the dependence of the maximal temperature on model parameters.

3. 3.

   We discuss the impact of thermal effects on expansion history and the CMB.

4. 4.

   We discuss the range of validity for our approach.

We will illustrate these points in the explicit models: the interaction between inflaton and other scalars or the gauge boson production from axion-like coupling. The details are discussed in Section V.

In the following section, we will consider the case where $n$ is an arbitrary positive integer, and give analytic solutions to the Boltzmann equations for the energy densities of the inflaton and radiation, $\rho_{\varphi}$ and $\rho_{\gamma}$. The latter can determine the temperature evolution during reheating. With the analytic solution of the effective temperature $T$, we discuss the maximal temperature in different cases. We will also compare the analytic results to the numerical ones, showing that the piece-wise approximation is sufficient to catch the important characteristics of the reheating process. In Section III, we fix the model parameters in realistic models. Section IV discuss the secenario in which there exists a preheating phase after inflation. The conclusions are presented in Section VI.

The Boltzmann equations for the energy densities $\rho_{\varphi}$ and $\rho_{\gamma}$, averaged over momentum and time across a few cycles of oscillations, are

	$\displaystyle\frac{d\rho_{\varphi}}{dt}+3H\rho_{\varphi}+\Gamma_{\varphi}\rho_{\varphi}=0$		(4)
	$\displaystyle\frac{d\rho_{\gamma}}{dt}+4H\rho_{\gamma}-\Gamma_{\varphi}\rho_{\varphi}=0$		(5)

with

$$\rho_{\gamma}=\frac{\pi^{2}g_{*}}{30}T^{4}~{}.$$

(6)

Introducing the variables $\Phi\equiv\rho_{\varphi}a^{3}/m_{\phi}$, $R\equiv\rho_{\gamma}a^{4}$ and $x\equiv am_{\phi}$, the above equations can be rewritten as

	$\displaystyle\frac{d\Phi}{dx}=-\frac{\Gamma_{\varphi}}{Hx}\Phi$		(7)
	$\displaystyle\frac{dR}{dx}=\frac{\Gamma_{\varphi}}{H}\Phi$		(8)

with

	$\displaystyle H=\left(\frac{8\pi}{3}\right)^{1/2}\frac{m_{\phi}^{2}}{M_{P}}\left(\frac{R}{x^{4}}+\frac{\Phi}{x^{3}}\right)^{1/2}$		(9)
	$\displaystyle T=\frac{m_{\phi}}{x}\left(\frac{30}{\pi^{2}g_{*}}R\right)^{1/4}~{}.$		(10)

Here $M_{P}$ and $g_{*}$ are the Planck mass and the number of degrees of freedom in the thermal bath, respectively. With these conventions, $x$ measures the dimensional scale factor $a$ in units of its value at the end of inflation, $x=a/a_{\rm end}$, and the initial value of $R$ is $R=0$ when $x=1$.

The reheating temperature $T_{R}$ is defined to be the temperature when reheating ends and $\rho_{\varphi}=\rho_{\gamma}$, i.e.

$$\frac{R}{x^{4}}=\frac{\Phi}{x^{3}}~{}.$$

(11)

This usually coincides in good approximation with the moment when $\Gamma_{\phi}=H$, prior to which the fractional energy loss of $\varphi$ is negligible and one can in good approximation set $\Phi=\Phi_{I}$ as a constant and view it as an external source in (8). If $\Gamma_{\varphi}$ is independent of $T$, i.e. $\Gamma_{\varphi}=\Gamma_{0}={\rm const}$, the usual estimated expression for $T_{R}$ can be found using $H=\sqrt{8\pi^{3}g_{*}/90}T^{2}/M_{P}$,

$$T_{R}=\sqrt{\Gamma_{0}M_{P}}\left(\frac{90}{8\pi^{3}g_{*}}\right)^{1/4}~{}.$$

(12)

In this work, we consider the situation where one of the temperature dependent terms dominates in the relevant range of temperatures

$$\Gamma_{\varphi}=\Gamma_{0}+\Gamma_{n}\left(\frac{T}{m_{\phi}}\right)^{n}~{}.$$

(13)

We use a local piece-wise approximation of $\Gamma_{\varphi}$ by monomials in $T$ when considering only one term at a time, and the thermal behaviour of the early universe after inflation can be divided into two phases, during which $\Gamma_{0}$ and $\Gamma_{n}$ dominate respectively. A critical value of temperature at which the finite temperature effects start to dominate $\Gamma_{\varphi}$ can be defined as

$$T_{n}=m_{\phi}\left(\frac{\Gamma_{0}}{\Gamma_{n}}\right)^{1/n}~{}.$$

(14)

Therefore the moment $x_{n}$ when $T=T_{n}$ can be viewed as the end of the $\Gamma_{0}$ phase, and the time-dependent temperature determined by $\Gamma_{0}$ at $x_{n}$ is the initial condition of the $\Gamma_{n}$ phase.

The solutions for the $\Gamma_{0}$ phase, i.e. the case without thermal feedback, can be found as Giudice:2000ex

	$\displaystyle R=\frac{2}{5}A_{0}\left(x^{5/2}-1\right)$		(15)
	$\displaystyle T=m_{\phi}\left(\frac{2}{5}A_{0}\frac{30}{\pi^{2}g_{*}}\right)^{1/4}\left(x^{-3/2}-x^{-4}\right)^{1/4}$		(16)
	$\displaystyle x_{\rm max}=\left(\frac{8}{3}\right)^{2/5}$		(17)
	$\displaystyle T_{\rm max}=T|_{x=x_{\rm max}}\approx 0.6\left(\frac{\Gamma_{0}M_{P}}{g_{*}}\right)^{1/4}V_{I}^{1/8}~{},$		(18)

here $V_{I}$ is the value of the inflaton potential at the end of inflation, $x_{\rm max}$ is the moment when $T$ reaches its maximum $T_{\rm max}$, and $A_{n}$ is defined to be

$$A_{n}=\frac{\Gamma_{n}M_{P}}{m_{\phi}^{2}}\sqrt{\Phi_{I}}\left(\frac{30}{\pi^{2}g_{*}}\right)^{n/4}\left(\frac{3}{8\pi}\right)^{1/2}~{}.$$

(19)

A useful relation between $A_{n}$ and $A_{0}$ can be obtained as

$$A_{n}=\frac{A_{0}\Gamma_{n}}{\Gamma_{0}}\left(\frac{30}{\pi^{2}g_{*}}\right)^{n/4}~{}.$$

(20)

If $T_{n}<T_{\rm max}$, the temperature reaches $T_{n}$ when $x=x_{n}$ and $R=R_{n}$, where

	$\displaystyle x_{n}\approx 1+\left(\frac{\Gamma_{0}}{\Gamma_{n}}\right)^{4/n}\frac{\pi^{2}g_{*}}{30A_{0}}~{},$		(21)
	$\displaystyle R_{n}=\frac{2}{5}A_{0}(x_{n}^{5/2}-1)\approx\left(\frac{\Gamma_{0}}{\Gamma_{n}}\right)^{4/n}\frac{\pi^{2}g_{*}}{30}~{}.$		(22)

Here we have expanded the solution (16) around $x=1$ since $1<x_{n}<x_{\rm max}\approx 1.48$. Otherwise if $T_{n}>T_{\rm max}$ then the temperature can never reach the value where the $\Gamma_{n}$ term starts to dominate. This means the solution with $\Gamma_{\varphi}=\Gamma_{0}={\rm const}$ is valid, and the thermal feedback is negligible. The results for this case can be seen in Figure 1 and 2. In Figure 1 we present the numerical solution for $T/m_{\phi}$ from (8) as the solid red line, and the analytic solution (16) as the dotted blue line, while the black dotted line represents the value of reheating temperature (12). The initial value of $T$ at the moment $x=1$ is set to be zero since we do not consider preheating here. With the numerical solution, Figure 2 shows the time dependence of $\rho_{\gamma}/\rho_{\varphi}$ and $\Gamma_{\phi}/H$ as the red and blue line respectively. It is clear that when $T_{n}>T_{\rm max}$, the reheating process is dominated by $\Gamma_{0}$ and (16) fully describes the time evolution of the effective temperature. The reheating process ends at the moment $\rho_{\varphi}=\rho_{\gamma}$, and this coincides well with $\Gamma_{\phi}=H$.

Below we will present the solutions for the $\Gamma_{n}$ dominated phase with the initial condition (22) for arbitrary $n$, assuming that $T_{n}<T_{\rm max}$. During this phase, the equation (8) is

$$\frac{dR}{dx}=A_{n}R^{n/4}x^{3/2-n}~{}.$$

(23)

In the previous section, we have given the analytic solutions as well as the comparison with the numerical results. Below we will first discuss how to determine the parameters in realistic models for the presented figures.

For the case $n=1$, let us now consider in a model in which the perturbative reheating is dominated by the axion-like coupling between the scalar inflaton field and the gauge boson. The interaction term in the Lagrangian has the form

$$\mathcal{L}_{\rm int}=\frac{\alpha}{\Lambda}\phi F_{\mu\nu}\tilde{F}^{\mu\nu}~{},$$

(53)

in which $F_{\mu\nu}$ is the field strength tensor of the vector bosons, $\alpha$ a dimensionless number and $\Lambda$ a mass scale. In the range where the inflaton mass is larger than the mass of the bosons and both are smaller than the temperature, the effective dissipation rate can be approximated as Carenza:2019vzg

$$\Gamma_{\varphi}\approx\frac{\alpha^{2}}{4\pi}\frac{m_{\phi}^{3}}{\Lambda^{2}}\left(1+2f_{B}\left(\frac{m_{\phi}}{2}\right)\right)~{}.$$

(54)

The Bose-Einstein distribution $f_{B}(y)$ can be expanded with $y<T$ to be $T/y$, and we get a linear term in temperature for $\Gamma_{\varphi}$,

$$\Gamma_{\varphi}\approx\frac{\alpha^{2}}{4\pi}\frac{m_{\phi}^{3}}{\Lambda^{2}}\left(1+\frac{4T}{m_{\phi}}\right)~{}.$$

(55)

We set the mass scale $\Lambda$ to be Planck mass $M_{p}$, the largest scale in our investigation. Introducing the dimensionless coupling $\tilde{\alpha}\equiv\alpha m_{\phi}/M_{P}$, we would like $\tilde{\alpha}$ to be neither too large to trigger the parametric resonance Drewes:2019rxn , nor too tiny to be unrealistic. In the numerical simulation we set $\tilde{\alpha}=10^{-6}$, $m_{\phi}=10^{9}~{}{\rm GeV}$ and $g_{*}=100$. Then in the above notation we have

$$\Gamma_{0}=\frac{1}{4}\Gamma_{1}=\frac{\tilde{\alpha}^{2}m_{\phi}}{4\pi}\approx 8\times 10^{-5}~{}{\rm GeV}.$$

(56)

For the estimation of the initial value of $\Phi_{I}$, considering that at the beginning of reheating $x=1$, the value of the inflaton potential can be expressed as

$$V_{I}\approx\frac{3}{4}\rho_{\varphi}=\frac{3}{4}\rho_{\varphi}x^{3}=\frac{3}{4}\rho_{\varphi}a^{3}m_{\varphi}^{3}=\frac{3}{4}\Phi_{I}m_{\phi}^{4}~{}.$$

(57)

If the oscillation happens in the quadratic regime, $V_{I}\sim m_{\phi}^{2}\varphi_{\rm end}^{2}/2$, together with the fact that the field value at the end of inflation $\varphi_{\rm end}$ is usually near the Planck scale, one has

$$\Phi_{I}=\frac{4V_{I}}{3m_{\phi}^{4}}\approx\frac{2m_{\phi}^{2}\varphi_{\rm end}^{2}}{3m_{\phi}^{4}}\sim\frac{M_{P}^{2}}{m_{\phi}^{2}}~{}.$$

(58)

Thus we can estimate that $\Phi_{I}=10^{20}$.

For the $n=2$ case, we consider that the perturbative reheating is dominated by the interaction

$$\mathcal{L}_{\rm int}=\frac{1}{3!}\alpha\phi\chi^{3}$$

(59)

between the inflaton and another scalar field $\chi$, where $\alpha$ is the dimensionless coupling constant. In the limit of vanishing $\chi$-quasiparticle width and with $m_{\chi}\ll m_{\phi}$, the effective dissipation rate can be approximated as Drewes:2013iaa

$$\Gamma_{\varphi}\approx\frac{\alpha^{2}m_{\phi}}{3072\pi^{3}}+\frac{\alpha^{2}T^{2}}{768\pi m_{\phi}}~{}.$$

(60)

Thus in our notation we have

$$\Gamma_{2}=4\pi^{2}\Gamma_{0}=\frac{\alpha^{2}m_{\phi}}{768\pi}~{}.$$

(61)

We set the parameters to be $\alpha=10^{-5}$, $m_{\phi}=10^{9}~{}{\rm GeV}$, $\Phi_{I}=10^{20}$ and $g_{*}=100$. In this set up,

$$\Gamma_{0}\approx 10^{-6}~{}{\rm GeV}~{}~{},~{}~{}\Gamma_{2}\approx 4\times 10^{-5}~{}{\rm GeV}~{}.$$

(62)

In Section II, we solve the Boltzmann equations for $\rho_{\varphi}$ and $\rho_{\gamma}$ in cases where $n$ is arbitrary. The initial value of $R$ (and thus $T$) at the beginning of reheating is set to be zero since we are only concerned with the perturbative reheating and the thermal feedback effects of the produced particles. However, the parametric resonance generally exists with relatively large couplings between the inflaton and other fields. With the resonant matter production in the short preheating phase before the reheating process, the initial conditions need to be modified. Assuming that the inflaton dissipates parts of its energy into other degrees in preheating, quantified as

$$\rho_{\gamma}=\epsilon\rho_{\varphi}~{},$$

(63)

where $0\leq\epsilon\leq 1$. If the parametric resonance is extremely efficient then $\epsilon=1$. On the other hand, $\epsilon=0$ corresponding to our previous assumption. Considering that the duration of preheating is usually much shorter than reheating, the initial conditions of the reheating process can be written as (the estimation of $\Phi_{I}$ still follows (58) at the end of inflation)

$$R_{I}=\rho_{\gamma}a^{4}=\epsilon\rho_{\varphi}a^{4}=\epsilon\Phi_{I}am_{\phi}\simeq\epsilon\Phi_{I}~{},$$

(64)

and

$$\Phi^{\prime}_{I}=(1-\epsilon)\rho_{\varphi}\frac{a^{3}}{m_{\phi}}=(1-\epsilon)\Phi_{I}~{}.$$

(65)

In this way the effective temperature at the beginning of reheating is then

$$T_{I}=\left(\frac{30}{\pi^{2}g_{*}}\rho_{\gamma}\right)^{1/4}=\left(\frac{30}{\pi^{2}g_{*}}\epsilon\Phi_{I}m_{\phi}^{4}\right)^{1/4}~{}.$$

(66)

Adopting the polynomial approximation of $\Gamma_{\varphi}$, one would find that the thermal term $\Gamma_{n}$ could dominate at the beginning of the reheating process, if $T_{I}$ is larger than $T_{n}$ in a certain $n$ case. And this corresponds to

$$\frac{30}{\pi^{2}g_{*}}\epsilon\Phi_{I}>\left(\frac{\Gamma_{0}}{\Gamma_{n}}\right)^{4/n}~{}.$$

(67)

Therefore the evolution of the effective temperature will be modified with the initial conditions determined by $\epsilon$ and the values of $\Gamma_{n}$ in practical models and the method used in the last section still applies.

Our work generalizes the attempt to analytically solve the Boltzmann equations for the energy density of inflaton and radiation. The novel points in this paper are

1. 1.

   We describe the entire reheating process.

2. 2.

   We present the dependence of the maximal temperature on model parameters.

3. 3.

   We discuss the impact of thermal effects on expansion history and the CMB.

4. 4.

   We discuss the range of validity for our approach.

We have illustrated these points in the following explicit models: the interaction between inflaton and other scalars or the gauge boson production from axion-like coupling. Below we will discuss them in details.