\section{Introduction}
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~\cite{Tanabashi:2018oca} and the abundance of the Standard Model matter~\cite{Davidson:2002qv, Fukugita:1986hr, Canetti:2014dka} as well as Dark Matter~\cite{Acharya:2008bk, Watson:2009hw, Acharya:2009zt}. In the  inflationary cosmology, reheating~\cite{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~\cite{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~\cite{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~\cite{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$~\cite{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~\ref{PRE}. In general, $\Gamma_\varphi$ has a complicated dependence on $T$ as shown in detailed studies~\cite{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~\cite{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 necessar
.

We parameterize the effects as
\begin{equation}
    \Gamma_\varphi=\sum_{n=0}^{\infty}\Gamma_n\left(\frac{T}{m_\phi}\right)^n~,
\end{equation}
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~\cite{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~\cite{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~\cite{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\rightarrow2$ decays. As will be seen, for a Bose enhanced $1\rightarrow2$ decay one generally has the form
\begin{equation}
    \Gamma_\varphi=\Gamma_0\left( 1+2f_B\left(\frac{m_\phi}{2}\right)\right)~,
\end{equation}
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\rightarrow3$ decays or $2\rightarrow2$ scatterings are included. For example in the $\alpha\phi\chi^3$ model one has~\cite{Drewes:2013iaa}
\begin{equation}
    \Gamma_\varphi\approx\frac{\alpha^2m_\phi}{3072\pi^3}+\frac{\alpha^2T^2}{768\pi m_\phi}~,
\end{equation}
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~\cite{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~\cite{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$~\cite{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
\begin{enumerate}
	\item We describe the entire reheating process.
	\item We present the dependence of the maximal temperature on model parameters.
	\item We discuss the impact of thermal effects on expansion history and the CMB.
	\item We discuss the range of validity for our approach.
\end{enumerate}
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~\ref{diss}.

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~\ref{par}, we fix the model parameters in realistic models. Section~\ref{PRE} discuss the secenario in which there exists a preheating phase after inflation. The conclusions are presented in Section~\ref{conc}.
\section{The solutions to the Boltzmann equations}\label{solu}
The Boltzmann equations for the energy densities $\rho_\varphi$ and $\rho_\gamma$, averaged over momentum and time across a few cycles of oscillations, are
\begin{eqnarray}
&&\frac{d\rho_\varphi}{dt}+3H\rho_\varphi+\Gamma_\varphi \rho_\varphi =0 \label{e22}\\
&&\frac{d\rho_\gamma}{dt}+4H\rho_\gamma-\Gamma_\varphi \rho_\varphi=0 \label{e23}
\end{eqnarray}
with
\begin{equation}
    \rho_\gamma=\frac{\pi^2 g_*}{30}T^4~. \label{e24}
\end{equation}
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
\begin{eqnarray}
    &&\frac{d\Phi}{dx}=-\frac{\Gamma_\varphi}{H x}\Phi \label{e25}\\
    &&\frac{d R}{dx}=\frac{\Gamma_\varphi}{H }\Phi \label{e26}
\end{eqnarray}
with
\begin{eqnarray}
    &&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} \label{e27}
    \\
    &&T=\frac{m_\phi}{x}\left(\frac{30}{\pi^2g_*}R\right)^{1/4}~. \label{e28}
    \end{eqnarray}
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.
\begin{equation}
    \frac{R}{x^4}=\frac{\Phi}{x^3} \label{R=Phi}~.
\end{equation}
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~(\ref{e26}). 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^3g_*/90}T^2/M_P$,
\begin{equation}
    T_R=\sqrt{\Gamma_0 M_P}\left(\frac{90}{8\pi^3g_*}\right)^{1/4}~.
    \label{TR}
\end{equation}

In this work, we consider the situation where one of the temperature dependent terms dominates in the relevant range of temperatures
\begin{equation}
\Gamma_\varphi=\Gamma_0+\Gamma_n\left(\frac{T}{m_\phi}\right)^n~.
\label{Gp}
\end{equation}
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
\begin{equation}
    T_n=m_\phi\left(\frac{\Gamma_0}{\Gamma_n}\right)^{1/n}~.
\end{equation}
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~\cite{Giudice:2000ex}
\begin{eqnarray}
    &&R=\frac{2}{5}A_0\left(x^{5/2}-1\right) \label{R0}\\
    &&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}\label{T0}\\
    &&x_{\rm max}=\left(\frac{8}{3}\right)^{2/5} \\
    &&T_{\rm max}=T|_{x=x_{\rm max}}\approx0.6\left(\frac{\Gamma_0M_P}{g_*}\right)^{1/4}V_I^{1/8} \label{Tmax}~,
\end{eqnarray}
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
\begin{equation}
    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}~.
    \label{An}
\end{equation}
A useful relation between $A_n$ and $A_0$ can be obtained as
\begin{equation}
    A_n=\frac{A_0\Gamma_n}{\Gamma_0}\left(\frac{30}{\pi^2g_*}\right)^{n/4}~.
\end{equation}

If $T_n<T_{\rm max}$, the temperature reaches $T_n$ when $x=x_n$ and $R=R_n$, where
\begin{eqnarray}
    &&x_n\approx1+\left(\frac{\Gamma_0}{\Gamma_n}\right)^{4/n}\frac{\pi^2 g_*}{30A_0}\label{x_i}~,\\
    &&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} \label{R_i}~.
\end{eqnarray}
Here we have expanded the solution (\ref{T0}) around $x=1$ since $1<x_n<x_{\rm max}\approx1.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~\ref{Tn=0} and~\ref{Hn=0}. In Figure~\ref{Tn=0} we present the numerical solution for $T/m_\phi$ from (\ref{e26}) as the solid red line, and the analytic solution (\ref{T0}) as the dotted blue line, while the black dotted line represents the value of reheating temperature (\ref{TR}). 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~\ref{Hn=0} 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 (\ref{T0}) 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$.
\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.8\linewidth]{Figures/T_n=0.pdf}
    \caption{The time dependence of $T/m_\phi$ with the double logarithmic coordinates when $\Gamma_\varphi=\Gamma_0$. We show
    the comparison of the analytic solution (\ref{T0}) to a numerical solution of (\ref{e25})-(\ref{e28}). The numerical solution for (\ref{e26}) is in the red solid line and the analytic solution (\ref{T0}) in the blue dotted line respectively. The black dotted line is the value of $T_R$ in (\ref{TR}), and the black dashed line is the value of $T_{\rm max}$ in (\ref{Tmax}). The choices for the parameters are: $m_\phi=10^9~{\rm GeV}$, $\Gamma_0=10^{-8}~{\rm GeV}$, $\Phi_I=10^{20}$ and $g_*=100$.}
    \label{Tn=0}
\end{figure}
\begin{figure}[htbp]
\centering
    \includegraphics[width=0.8\linewidth]{Figures/H_n=0.pdf}
    \caption{The time dependence of $\rho_\gamma/\rho_\varphi$ in the red line and $\Gamma_\varphi/H$ in the blue line with the double logarithmic coordinates when $\Gamma_\varphi=\Gamma_0$. After the moment $\rho_\gamma/\rho_\varphi=1$, the reheating ends and the universe becomes radiation dominated. It also coincides with the moment $\Gamma_\varphi=H$ and when (\ref{T0}) starts to fail. The choices for the parameters are: $m_\phi=10^9~{\rm GeV}$, $\Gamma_0=10^{-8}~{\rm GeV}$, $\Phi_I=10^{20}$ and $g_*=100$.}
    \label{Hn=0}
\end{figure}

Below we will present the solutions for the $\Gamma_n$ dominated phase with the initial condition (\ref{R_i}) for arbitrary $n$, assuming that $T_n<T_{\rm max}$. During this phase, the equation (\ref{e26}) is
\begin{equation}
    \frac{dR}{dx}=A_n R^{n/4}x^{3/2-n} \label{Req}~.
\end{equation}
\subsection{Case $n<4$}
In this case, the solution of the above equation is
\begin{equation}
    R=\left(A_n\frac{1-n/4}{5/2-n}\left(x^{5/2-n}-x^{5/2-n}_n\right)+R^{1-n/4}_n\right)^{\frac{4}{4-n}}~,
    \label{me28}
\end{equation}
and therefore the solution for the temperature during the $\Gamma_n$ dominated phase is
\begin{equation}
    T=m_\phi\left(\frac{30}{\pi^2g_*}\right)^{1/4}\frac{1}{x}\left(A_n\frac{1-n/4}{5/2-n}\left(x^{5/2-n}-x^{5/2-n}_n\right)+R^{1-n/4}_n\right)^{\frac{1}{4-n}}~.
    \label{me29}
\end{equation}
If $T_R<T_n<T_{\rm max}$, then the temperature $T_{\rm max}$ given by (\ref{Tmax}) is no longer applicable when $x>x_n$ since the $\Gamma_n$-term cannot be neglected after $T$ reaches $T_n$. The solution for all time will be a piecewise function given by (\ref{T0}) and (\ref{me29}), and the temperature reaches its maximum at the moment
\begin{equation}\begin{aligned}
    \tilde{x}_{\rm max}&=\left(\frac{8-2n}{3}x^{5/2-n}_n+\frac{4(2n-5)}{3A_n}R^{1-n/4}_n\right)^{\frac{1}{5/2-n}}\\
    &\approx\left(\frac{8-2n}{3}\right)^{\frac{2}{5-2n}}x_n^\frac{n}{n-4}~.
    \label{me30}
\end{aligned}
\end{equation}
Considering that
\begin{equation}\begin{aligned}
 &A_n\frac{1-n/4}{5/2-n}\left(x^{5/2-n}_{\rm max}-x^{5/2-n}_n\right)+R^{1-n/4}_n \\
 =&A_n\frac{1-n/4}{5/2-n}\left(\frac{8-2n}{3}x^{5/2-n}_n+\frac{4(2n-5)}{3A_n}R_n^{1-n/4}-x^{5/2-n}_n\right)+R_n^{1-n/4} \\
 =&\frac{2(1-n/4)}{3}A_nx^{5/2-n}_n+\frac{2(n-4)}{3}R^{1-n/4}_n+R^{1-n/4}_n \\
 =&\frac{2(1-n/4)}{3}A_nx_n^{5/2-n}+\frac{2n-5}{3}R^{1-n/4}_n \\
 =&\frac{A_n}{4}\tilde{x}_{\rm max}^{5/2-n}~,
 \end{aligned}
\end{equation}
one gets
\begin{equation}\begin{aligned}
    \tilde{T}_{\rm max}&=m_\phi\left(\frac{30}{\pi^2g_*}\right)^{1/4}\tilde{x}^{-1}_{\rm max}\left(\frac{A_n}{4}\tilde{x}_{\rm max}^{5/2-n}\right)^{\frac{1}{4-n}} \\
    &=m_\phi\left(\frac{30}{\pi^2g_*}\right)^{1/4}\tilde{x}^{-\frac{3}{2(4-n)}}_{\rm max}\left(\frac{A_n}{4}\right)^{\frac{1}{4-n}} \\
    &\approx m_\phi\left(\frac{30}{\pi^2g_*}\right)^{1/4}\left(\frac{A_n}{4}\right)^{\frac{1}{4-n}}\left(\frac{8-2n}{3}\right)^{\frac{3}{(n-4)(5-2n)}}x_n^\frac{3n}{2(n-4)^2}
    ~.
    \label{me32}
    \end{aligned}
\end{equation}

In certain $n$ case, we can express $\tilde{T}_{\rm max}$ in terms of the model parameters \{$m_\phi$, $\Gamma_0$, $\Gamma_n$, $g_*$ and $\Phi_I$\} and $\Gamma_n$ further in terms of the inflaton coupling in realistic models. This will be discussed in more details in Section~\ref{par} and~\ref{diss}.

If $T_n<T_R<T_{\rm max}$, then the universe gets reheated very quickly due to the so called thermal resonance, and the radiation dominated era ($\rho_\gamma>\rho_\varphi$) starts in the $\Gamma_n$ regime
. In this situation, the maximum temperature and the reheating temperature are both roughly $\tilde{T}_{\rm max}$ because the universe reheats almost instantaneously at the moment $x=x_{\rm crit}$ when (\ref{R=Phi}) stands~\cite{Drewes:2014pfa}, i.e. $R=\Phi_I x$, and it corresponds to the solution of
\begin{equation}
    \left(A_n\frac{1-n/4}{5/2-n}\left(x^{5/2-n}-x^{5/2-n}_n\right)+R^{1-n/4}_n\right)^{\frac{4}{4-n}}=\Phi_I x~.
    \label{eqx_cr}
\end{equation}
In the following we will focus on the cases where $T_R<T_n<T_{\rm max}$ and give the analytic solutions for $T$ as well as its maximum $\tilde{T}_{\rm max}$.

More specifically, in the case $n<4$ we have $n=1,~2,~3$, below we present the results respectively. In the $n=1$ and $n=2$ cases, we will pin down the parameters by physical arguments and compare the analytic solution with the numerical one in realistic models.

\begin{itemize}
    \item $n=1$

The solutions (\ref{me28}), (\ref{me29}), (\ref{me30}) and (\ref{me32}) are
    \begin{eqnarray}
    &&R=\left(\frac{A_1}{2}\left(x^{3/2}-x^{3/2}_1\right)+R_1^{3/4}\right)^{4/3}\\
    &&T=m_\phi\left(\frac{30}{\pi^2g_*}\right)^{1/4}\frac{1}{x}\left(\frac{A_1}{2}\left(x^{3/2}-x^{3/2}_1\right)+R_1^{3/4}\right)^{1/3}\label{T1}\\
    &&\tilde{x}_{\rm max}=\left(2x_1^{3/2}-4\frac{R_1^{3/4}}{A_1}\right)^{2/3}\approx\left(\frac{4}{x_1}\right)^{1/3}\\
    &&\tilde{T}_{\rm max}\approx m_\phi\left(\frac{30}{\pi^2g_*}\right)^{1/4}\left(\frac{A_1}{8}\right)^{1/3}x_1^{1/6} \label{TM1}~,
    \end{eqnarray}
while the boundary condition is
\begin{equation}
    T_1=m_\phi\frac{\Gamma_0}{\Gamma_1}~~~~{\rm at}~~~~x_1\approx1+\left(\frac{\Gamma_0}{\Gamma_1}\right)^4\frac{\pi^2g_*}{30A_0}.\label{Ti1}
\end{equation}

Below we present the comparison between the analytic solution and the numerical solution with the above set up. The temperature evolution can be seen in Figure~\ref{Tn=1}, and the evolution for $\rho_\gamma/\rho_\varphi$ and $\Gamma_\varphi/H$ in Figure~\ref{Hn=1}. At the beginning, the temperature is below $T_1$ and $\Gamma_\varphi$ is dominated by $\Gamma_0$, thus the evolution of $T$ is the same as that of the vacuum decays, (\ref{T0}). In the range $T_1<T<\tilde{T}_{\rm max}$, the dissipation can be described by the $\Gamma_1$ term, and the solution is (\ref{T1}) with the boundary conditions (\ref{Ti1}). It is clear that the $\Gamma_0$ term will dominate again after the temperature in the $\Gamma_1$ regime drops below $T_1$, and this extra $\Gamma_0$ phase holds until it reaches $T_R$. When $T_R<T_1<T_{\rm max}$ is fulfilled by the parameters, the piece-wise approximation for $\Gamma_\varphi$ works extremely well compared to the numerical solution.
\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.8\linewidth]{Figures/T_n=1.pdf}
    \caption{The time dependence of $T/m_\phi$ with the double logarithmic coordinates when $\Gamma_\varphi=\Gamma_0+\Gamma_1(T/m_\phi)$. The numerical solution for (\ref{e26}) is in the red solid line, the analytic solution (\ref{T0}) in the blue dotted line and (\ref{T1}) in the blue dashed line respectively. The black dotted line is the value of $T_R$ in (\ref{TR}), the black dashed line is the value of $T_{\rm max}$ in (\ref{Tmax}), the black solid line is the value of $T_1$ in (\ref{Ti1}), and the black dotted-dashed line is the value of $\tilde{T}_{\rm max}$ in (\ref{TM1}). The choices for the parameters are: $m_\phi=10^9~{\rm GeV}$, $\Gamma_0=8\times10^{-5}~{\rm GeV}$, $\Gamma_1=3.2\times10^{-4}~{\rm GeV}$, $\Phi_I=10^{20}$ and $g_*=100$. The choice for these values will be discussed within a realistic model in Section~\ref{par}.}

    \label{Tn=1}
\end{figure}
\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.8\linewidth]{Figures/H_n=1.pdf}
    \caption{The time dependence of $\rho_\gamma/\rho_\varphi$ in the red line and $\Gamma_\varphi/H$ in the blue line when $\Gamma_\varphi=\Gamma_0+\Gamma_1(T/m_\phi)$. The moment $\rho_\gamma/\rho_\varphi=1$ happens when $\Gamma_\varphi=H$ and reheating ends. The choices for the parameters are the same with~\ref{Tn=1}: $m_\phi=10^9~{\rm GeV}$, $\Gamma_0=8\times10^{-5}~{\rm GeV}$, $\Gamma_1=3.2\times10^{-4}~{\rm GeV}$, $\Phi_I=10^{20}$ and $g_*=100$.}

    \label{Hn=1}
\end{figure}

\item $n=2$

The solutions (\ref{me28}), (\ref{me29}), (\ref{me30}) and (\ref{me32}) are
 \begin{eqnarray}
 &&R=\left(A_2\left(\sqrt{x}-\sqrt{x_2}\right)+\sqrt{R_2}\right) \\
 &&T=m_\phi\left(\frac{30}{\pi^2g_*}\right)^{1/4}\frac{1}{x}\sqrt{A_2\left(\sqrt{x}-\sqrt{x_2}\right)+R_2}\label{T2}\\
  &&\tilde{x}_{\rm max}=\frac{\left(4A_2\sqrt{x_2}-4\sqrt{R_2}\right)^2}{9A_2^2}\approx\frac{16}{9x_2}\\
 &&\tilde{T}_{\rm max}=m_\phi\left(\frac{30}{\pi^2g_*}\right)^{1/4}\left(\frac{3}{4}\right)^2\sqrt{\frac{A_2}{3}}x_2^{3/4}\label{TM2}~,
 \end{eqnarray}
 while the boundary condition is
 \begin{equation}
     T_2=m_\phi\sqrt{\frac{\Gamma_0}{\Gamma_2}}~~~~{\rm at}~~~~x_2\approx1+\left(\frac{\Gamma_0}{\Gamma_2}\right)^{2}\frac{\pi^2 g_*}{30A_0}.\label{Ti2}
 \end{equation}

Comparing our results of case $n=2$ with those in the section 2.2.4 of~\cite{Drewes:2014pfa}, we find that there exist some minor differences in this previous research for the expressions of $\tilde{x}_{\rm max}$, $\tilde{T}_{\rm max}$ and $x_2$.

Again, we will extract the parameters in an illustrative model.
and present the comparison between the analytic solution and the numerical solution.

The time evolution of the temperature and the ratio between energy densities are shown in Figure~\ref{Tn=2} and~\ref{Hn=2}.

Figure~\ref{Tn=2} shows the time evolution of the effective temperature in three stages during reheating. At the beginning, the temperature is below $T_2$ and $\Gamma_\varphi$ is dominated by $\Gamma_0$, thus the evolution of $T$ is the same as that of the vacuum decays, (\ref{T0}). In the range $T_2<T<\tilde{T}_{\rm max}$, the dissipation can be described by the $\Gamma_2$ term, and the solution is (\ref{T2}) with the boundary conditions (\ref{Ti2}). The $\Gamma_0$ term will dominate again after the temperature drops below $T_2$ until the reheating ends at $T_R$.
 \begin{figure}[htbp]
    \centering
    \includegraphics[width=0.8\linewidth]{Figures/T_n=2.pdf}
    \caption{The time dependence of $T/m_\phi$ with the double logarithmic coordinates when $\Gamma_\varphi=\Gamma_0+\Gamma_2(T/m_\phi)^2$. The numerical solution for (\ref{e26}) is in the red solid line, the analytic solution (\ref{T0}) in the blue dotted line and (\ref{T2}) in the blue dashed line respectively. The black dotted line is the value of $T_R$ in (\ref{TR}), the black dashed line is the value of $T_{\rm max}$ in (\ref{Tmax}), the black solid line is the value of $T_2$ in (\ref{Ti2}), and the black dotted-dashed line is the value of $\tilde{T}_{\rm max}$ in (\ref{TM2}). The choices for the parameters are: $m_\phi=10^9~{\rm GeV}$, $\Gamma_0=10^{-6}~{\rm GeV}$, $\Gamma_2=4\times10^{-5}~{\rm GeV}$, $\Phi_I=10^{20}$ and $g_*=100$. The choice for these values will be discussed within a realistic model in Section~\ref{par}.}

    \label{Tn=2}
\end{figure}
 \begin{figure}[htbp]
    \centering
    \includegraphics[width=0.8\linewidth]{Figures/H_n=2.pdf}
    \caption{The time dependence of $\rho_\gamma/\rho_\varphi$ in the red line and $\Gamma_\varphi/H$ in the blue line when $\Gamma_\varphi=\Gamma_0+\Gamma_2(T/m_\phi)^2$. The moment $\rho_\gamma/\rho_\varphi=1$ happens when $\Gamma_\varphi=H$ and reheating ends. The choices for the parameters are the same with~\ref{Tn=2}: $m_\phi=10^9~{\rm GeV}$, $\Gamma_0=10^{-6}~{\rm GeV}$, $\Gamma_2=4\times10^{-5}~{\rm GeV}$, $\Phi_I=10^{20}$ and $g_*=100$.}

    \label{Hn=2}
\end{figure}

  \item $n=3$

The solutions (\ref{me28}), (\ref{me29}), (\ref{me30}) and (\ref{me32}) are
  \begin{eqnarray}
  &&R=\left(-\frac{A_3}{2}\left(\frac{1}{\sqrt{x}}-\frac{1}{\sqrt{x_3}}\right)+R_3^{1/4}\right)^4\\
  &&T=m_\phi\left(\frac{30}{\pi^2g_*}\right)^{1/4}\frac{1}{x}\left(\frac{A_3}{2\sqrt{x_3}}+R_3^{1/4}-\frac{A_3}{2\sqrt{x}}\right)\label{T3}\\
  &&\tilde{x}_{\rm max}=\left(\frac{2}{3\sqrt{x_3}}+\frac{4}{3A_3}R_3^{1/4}\right)^{-2}\approx\frac{9}{4x^{3}_3}\\
  &&\tilde{T}_{\rm max}=m_\phi\left(\frac{30}{\pi^2g_*}\right)^{1/4}\frac{A_3}{4}\left(\frac{2}{3}\right)^2x_3^{9/2}~,\label{TM3}
  \end{eqnarray}
  while the boundary condition is
  \begin{equation}
      T_3=m_\phi\left(\frac{\Gamma_0}{\Gamma_3}\right)^{1/3}~~~~{\rm at}~~~~x_3=1+\left(\frac{\Gamma_0}{\Gamma_3}\right)^{4/3}\frac{\pi^2 g_*}{30A_0}~.\label{Ti3}
  \end{equation}

\end{itemize}

\subsection{Case $n=4$}
In this case, the equation (\ref{Req}) becomes
    \begin{equation}
        \frac{dR}{dx}=A_4 R x^{-5/2}~,
    \end{equation}
    and its solution  is
    \begin{equation}
        R=R_4 \exp\left(-\frac{2}{3}A_4\left(x^{-3/2}-x^{-3/2}_4\right)\right)~.\label{R4}
    \end{equation}
Therefore the time dependent function of the temperature is \begin{equation}
     T=m_\phi\left(\frac{30R_4}{\pi^2g_*}\right)^{1/4}\frac{1}{x}\exp\left(-\frac{1}{6}A_4\left(x^{-3/2}-x^{-3/2}_4\right)\right) \label{T4}~.
 \end{equation}
With the boundary condition
\begin{equation}
T_4=m_\phi\left(\frac{\Gamma_0}{\Gamma_4}\right)^{1/4}~~~~{\rm at}~~~~x_4=1+\frac{\Gamma_0}{\Gamma_4}\frac{\pi^2 g_*}{30A_0}~,\label{Ti4}
\end{equation}
the maximal temperature can be found to be at the moment $\tilde{x}_{\rm max}=(A_4/4)^{2/3}$,
\begin{equation}
\tilde{T}_{\rm max}=m_\phi\left(\frac{30R_4}{\pi^2g_*}\right)^{1/4}\left(\frac{4}{A_4}\right)^{2/3}\exp\left(-\frac{2}{3}+\frac{A_4}{6}x_4^{-3/2}\right)~.
\label{TM4}
\end{equation}

\subsection{Case $n>4$}
 In this case, the solutions for $R$ and $T$ become
\begin{eqnarray}
&R=\left|A_n\frac{1-n/4}{5/2-n}\left(x^{5/2-n}-x^{5/2-n}_n\right)+R^{1-n/4}_n\right|^{\frac{4}{4-n}}\label{R<4}  \\
&T=m_\phi\left(\frac{30}{\pi^2g_*}\right)^{1/4}\frac{1}{x}\left|A_n\frac{1-n/4}{5/2-n}\left(x^{5/2-n}-x^{5/2-n}_n\right)+R^{1-n/4}_n\right|^{\frac{1}{4-n}}~.\label{T5}
\end{eqnarray}
Here the reason for the absolute value function is that the variable inside may be negative, and we need to be careful when calculating the roots. The discussion for $n<4$ applies when the argument in the absolute value function is positive. However, one should notice that the above function of $T$ can have a pole in a certain range of parameters because the argument appears in the denominator due to $1/(4-n)<0$ and
\begin{equation}
    \begin{aligned}
    &R_n^{1-n/4}-\frac{1-n/4}{5/2-n}A_nx_n^{5/2-n} \\
    =&A_n\left(\frac{R_n^{1-n/4}}{A_n}-\frac{1-n/4}{5/2-n}x_n^{5/2-n}\right)\\
    \approx&A_n\left(\left(\frac{\Gamma_0}{\Gamma_n}\right)^{4/n-1}\left(\frac{\pi^2g_*}{30}\right)^{1-n/4}\frac{\Gamma_0}{A_0\Gamma_n}\left(\frac{\pi^2g_*}{30}\right)^{n/4}\right.\\
 \hspace{2cm}
    &\left.-\frac{1-n/4}{5/2-n}\left(1+\left(\frac{5}{2}-n\right)\left(\frac{\Gamma_0}{\Gamma_n}\right)^{4/n}\frac{\pi^2g_*}{30A_0}\right)\right)\\
    =&A_n\left(\frac{n/4-1}{5/2-n}+\frac{n}{4}\left(\frac{\Gamma_0}{\Gamma_n}\right)^{4/n}\frac{\pi^2g_*}{30A_0}\right)~.
    \end{aligned}
\end{equation}
The above expression will be negative when the second term in the bracket is smaller than $(n/4-1)/(n-5/2)$, while the latter goes from $0.1$ to $0.25$ in the range $n\geq5$. Thus in this situation, the denominator of (\ref{T5}) would be zero at some finite moment $x_0$. After the temperature reaches $T_n$, it will increase monotonically to infinity as $x\rightarrow x_0$. However, this solution was obtained under the assumption that the change in $\Phi_I$ is neglected, which is a good assumption until $\Gamma_\varphi= H$. This would break down before the temperature reaches the divergence simply because of energy conservation. The temperature must remain lower than the temperature that one would get if the inflaton would instantaneously transfer all its energy into radiation.
This means that reheating would end soon after $x_n$ at some moment $x_{\rm crti}$, the solution of (\ref{eqx_cr}), since the energy of inflaton field is finite and thus the dissipation ends before $x_0$. This corresponds to $T_n<T_R$. In most cases, there does not exist an analytic solution for $x_{\rm crit}$, and the maximal temperature and the reheating temperature are both given by (\ref{T5}) at $x_{\rm crti}$.

\section{Fixing the parameters in realistic models}\label{par}
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
\begin{equation}
	\mathcal{L}_{\rm int}=\frac{\alpha}{\Lambda}\phi F_{\mu\nu}\tilde{F}^{\mu\nu}~,
\end{equation}
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~\cite{Carenza:2019vzg}
\begin{equation}
	\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)~.
\end{equation}
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$,
\begin{equation}
	\Gamma_\varphi\approx\frac{\alpha^2}{4\pi}\frac{m_\phi^3}{\Lambda^2}\left(1+\frac{4T}{m_\phi}\right)~.
\end{equation}
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~\cite{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
\begin{equation}
	\Gamma_0=\frac{1}{4}\Gamma_1=\frac{\tilde{\alpha}^2m_\phi}{4\pi}\approx8\times10^{-5}~{\rm GeV}.\label{G1A}
\end{equation}
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
\begin{equation}
	V_I\approx\frac{3}{4}\rho_\varphi=\frac{3}{4}\rho_\varphi x^3=\frac{3}{4}\rho_\varphi a^3m_\varphi^3=\frac{3}{4}\Phi_I m_\phi^4~.\label{VI}
\end{equation}
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
\begin{equation}
	\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}~.\label{PI}
\end{equation}
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
\begin{equation}
	\mathcal{L}_{\rm int}=\frac{1}{3!}\alpha\phi\chi^3
\end{equation} 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~\cite{Drewes:2013iaa}
\begin{equation}
	\Gamma_\varphi\approx\frac{\alpha^2m_\phi}{3072\pi^3}+\frac{\alpha^2T^2}{768\pi m_\phi}~.
\end{equation}
Thus in our notation we have
\begin{equation}
	\Gamma_2=4\pi^2\Gamma_0=\frac{\alpha^2 m_\phi}{768\pi}~.
\end{equation}
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,
\begin{equation}
	\Gamma_0\approx10^{-6}~{\rm GeV}~~,~~\Gamma_2\approx4\times10^{-5}~{\rm GeV}~.
\end{equation}
\section{The scenarios with the existence of preheating}\label{PRE}
In Section~\ref{solu}, 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
\begin{equation}
    \rho_\gamma=\epsilon\rho_\varphi~,
\end{equation}
where $0\leq\epsilon\leq1$. 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 (\ref{PI}) at the end of inflation)
\begin{equation}
    R_I=\rho_\gamma a^4=\epsilon\rho_\varphi a^4=\epsilon\Phi_I a m_\phi\simeq\epsilon\Phi_I~,
\end{equation}
and
\begin{equation}
    \Phi'_I=(1-\epsilon)\rho_\varphi \frac{a^3}{m_\phi}=(1-\epsilon)\Phi_I~.
\end{equation}
In this way the effective temperature at the beginning of reheating is then
\begin{equation}
    T_I=\left(\frac{30}{\pi^2g_*}\rho_\gamma\right)^{1/4}=\left(\frac{30}{\pi^2g_*}\epsilon\Phi_I m_\phi^4\right)^{1/4}~.
\end{equation}
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
\begin{equation}
    \frac{30}{\pi^2g_*}\epsilon\Phi_I >\left( \frac{\Gamma_0}{\Gamma_n}\right)^{4/n}~.
\end{equation}
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.

\section{Discussions}\label{diss}
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
\begin{enumerate}
	\item We describe the entire reheating process.
	\item We present the dependence of the maximal temperature on model parameters.
	\item We discuss the impact of thermal effects on expansion history and the CMB.
	\item We discuss the range of validity for our approach.
\end{enumerate}
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.

\subsection{The description of the entire reheating process}
To solve the entire reheating process we need to locally use the power law approximation and match the solutions, as done in Section~\ref{solu} by the piece-wise approximation. In general, a $\Gamma_0$ piece is needed otherwise it is impossible to reheat the universe unless one assumes the presence of an intial radiation bath from preheating. Besides, the piece-wise approximation can be extended to the cases where $\Gamma_\varphi$ has a complicated dependence on $T$ other than a power law, and we can expand it to contain more terms and match the solutions piece by piece.

In the situation where the temperature starts from $T=0$, we present the proper matching of the piece-wise solutions. For example, in the case $n=1$,
Figure~\ref{Tn=1} shows the time evolution of the effective temperature in three stages during reheating. At the beginning, the temperature is below $T_1$ and $\Gamma_\varphi$ is dominated by $\Gamma_0$, thus the evolution of $T$ is the same as that of the vacuum decays, (\ref{T0}). In the range $T_1<T<\tilde{T}_{\rm max}$, the dissipation can be described by the $\Gamma_1$ term, and the solution is (\ref{T1}) with the boundary conditions (\ref{Ti1}). The $\Gamma_0$ term will dominate again after the temperature drops below $T_1$ until the reheating ends at $T_R$. Figure~\ref{Hn=1} also indicates the locally power law approximation works well during the entire reheating process.

\subsection{The dependence of the maximal temperature on model parameters}
Our analytic results permits to discuss explicitly the dependence of the maximal temperature on model parameters. Taking (\ref{TM1}) as an example, in the case $n=1$ the dependence of $\tilde{T}_{\rm max}$  on the model parameters \{$m_\phi$, $\Gamma_0$, $\Gamma_1$, $g_*$ and $\Phi_I$\} can be obtained explicitly by (\ref{An}) and (\ref{Ti1}),
\begin{equation}
	\begin{aligned}
	\tilde{T}_{\rm max}&\approx m_\phi\left(\frac{30}{\pi^2g_*}\right)^{1/4}\left(\frac{A_1}{8}\right)^{1/3}x_1^{1/6}\\
	&=\frac{m_\phi}{2}\left(\frac{30}{\pi^2g_*}\right)^{1/3}\frac{\Gamma_1^{1/3}M_P^{1/3}\Phi_I^{1/6}}{m_\phi^{2/3}}\left(\frac{3}{8\pi}\right)^{1/6}\left(1+\left(\frac{\Gamma_0}{\Gamma_1}\right)^4\frac{\pi^2g_*}{30A_0}\right)^{1/6}\\
	&\approx\left(\frac{3}{8\pi}\right)^{1/6}\left(\frac{30}{\pi^2}\right)^{1/3}\frac{m_\phi^{1/3}\Gamma_1^{1/3}M_P^{1/3}\Phi_I^{1/6}}{2g_*^{1/3}} ~,\label{dep}
	\end{aligned}
\end{equation}
the last line is calculated by ignoring the second term in the bracket above,
considering that $1<x_1<x_{\rm max}\approx1.48$ and thus $x_1^{1/6}\approx1$. Together with (\ref{G1A}) and (\ref{VI}), we have in the axion-like model
\begin{equation}
	\tilde{T}_{\rm max}\propto \alpha^{2/3} m_\phi^{2/3}V_I^{1/6}g_*^{-1/3}~.
\end{equation}
It is clear that the maximal temperature increases with the axion-like coupling, the inflaton mass and the value of inflaton potential at the end of inflation, and decreases with the number of degrees of freedom in the bath. If more accurate dependence is required, one can take the second term in the bracket of (\ref{dep}) into account and get another term depending on $\Gamma_0$, which further depends on $\alpha$.

\subsection{The impact of thermal effects on expansion history and the CMB}
The thermal feedback of the produced particles has impact not only on the thermal history of the universe, but also on the expansion history. The expansion history can be indeed affected significantly for $n<2$ and the reheating ends much earlier than expected. To see this, we consider that the Hubble parameter $H$ is larger than the dissipation rate $\Gamma_\varphi$ at the moment $\tilde{x}_{\rm max}$ when the temperature reaches its maximum. In the case $n\ge2$, $H$ approximately scales as $\propto T^2$ near the end of reheating and $\Gamma_\varphi$ scales as $\propto T^n$ during the $\Gamma_n$ dominated stage. Thus $\Gamma_\varphi$ decreases faster than $H$ with a decreasing $T$, and $\Gamma_\varphi=H$ happens during the $\Gamma_0$ dominated stage after $T$ drops below $T_n$, as shown in Section~\ref{solu}. However, in the $n=1$ case, $\Gamma_\varphi$ decreases slower than $H$ and it is probable that they meet each other in the $\Gamma_1$ dominated stage. In this situation, the duration of reheating is significantly shortened and the reheating temperature has to be modified, corresponding to $T_1<T_R<\tilde{T}_{\rm max}$.

We can further find the constraints on $\Gamma_1$ using the solution (\ref{T1}) of $T$ during the $\Gamma_1$ stage. Together with (\ref{e27}) and that $\Gamma_\varphi\simeq\Gamma_1 T/m_\phi$, the reheating temperature $T_R$ is obtained at the moment $\Gamma_\varphi=H$,
\begin{equation}
	T_R=\frac{\delta+\sqrt{\delta^2+4\gamma^3A_1}}{2\gamma^2}~,
\end{equation}
where
\begin{eqnarray}
	&&\delta\equiv\frac{4A_1}{m_\phi^3}\left(\frac{\pi^2g_*}{30}\right)^{3/4}\left(-\frac{1}{2}+\frac{1}{4}\left(\frac{\Gamma_0}{\Gamma_1}\right)^4\frac{\pi^2g_*}{30A_0}\right)\\
	&&\gamma\equiv2\sqrt{\frac{8\pi}{3}}\frac{\sqrt{\Phi_I}}{\Gamma_1M_P}\left(\frac{\pi^2g_*}{30}\right)^{3/4}~.
\end{eqnarray}
Considering that $T_1\equiv m_\phi\Gamma_0/\Gamma_1$, the condition $T_1<T_R$ then gives $\Gamma_1>\Gamma_{\rm 1,crit}$ with
\begin{equation}
	\begin{aligned}
	\Gamma_{1,{\rm crit}}&=\left(\frac{32\pi\Phi_Im_\phi^2\Gamma_0^2}{3M_P^2}-\frac{\pi^2g_*\Gamma_0^4}{30m_\phi^2}\right)^{1/4}\left(\frac{2\Phi_I}{m_\phi^2}-\frac{2\Gamma_0M_P\sqrt{\Phi_I}}{m_\phi^4}\sqrt{\frac{3}{8\pi}}\right)^{-1/4}\\
	&\approx\left(\frac{16\pi}{3}\right)^{1/4}m_\phi\sqrt{\frac{\Gamma_0}{M_P}}~.
	\end{aligned}
\end{equation}
The last step uses $\Gamma_0\ll m_\phi$ and $\Phi_I\sim M^2_P/m_\phi^2$ to drop the second terms in the two brackets.

The reheating phase affects the cosmic microwave background (CMB) via its effects on the expansion history of the universe, since the equation of state during reheating is different from the one in inflationary or radiation dominated era, and the expansion history affect the way how physical scales at present time and during inflation relate to each other.
While we can only observe the time-integrated effect, the quantities that the CMB is directly sensitive to are the duration of the reheating era in terms of the e-folds $N_{\rm re}$ and the averaged equation of state $\bar{\omega}_{\rm re}$ during reheating, with~\cite{Drewes:2019rxn}
\begin{equation}
	\bar{\omega}_{\rm re}=\frac{1}{ N_{\rm re}}\int_0^{ N_{\rm re}}\omega(N)dN~.
\end{equation}
In the case $n=1$, when $\Gamma_1$ is so small that $T_1>T_{\rm max}$, the temperature evolution will just be the one of perturbative decay, a $\Gamma_0$ dominated stage. With $\Gamma_1$ being in a intermediate range, the reheating process contains several stages as decsribed. However, if $\Gamma_1>\Gamma_{1,{\rm crit}}$, the duration of reheating $N_{\rm re}$ will be significantly shortened, the thermal feedback effects has impact on not only the thermal history but also the expansion history, and therefore affects the predication of the observable CMB.

\subsection{The discussion of the range of validity}
We now discuss under what circumstances can our approach be applied. In the $n=1$ case, if we assume the perturbative reheating is dominated by the axion-like coupling, then by (\ref{G1A}) we have
\begin{equation}
	T_1\equiv m_\phi\frac{\Gamma_0}{\Gamma_1}=\frac{m_\phi}{4}~,
\end{equation}
while
\begin{equation}
\Gamma_0=\frac{1}{4}\Gamma_1=\frac{\alpha^2m_\phi^3}{4\pi M_P^2}~.
\end{equation}
To avoid the parametric resonance, the coupling $\alpha$ needs to be relatively small. In the Section 3.4 of~\cite{Drewes:2019rxn}, this upper bound can be found as
\begin{equation}
\tilde{\alpha}\equiv\alpha\frac{m_\phi}{\Lambda}\ll\min\left(\frac{m_\phi}{4\varphi_{\rm end}},~\frac{m_\phi}{5\sqrt{\varphi_{\rm end}M_P}}\right)~.
\end{equation}
If $\Lambda$ and $\varphi_{\rm end}$ are both set to be $M_P$, then one gets
\begin{equation}
	\alpha\ll0.2~.
\end{equation}
On the other hand, with $V_I\sim m_\phi^2 M_P^2/2$ the $T_{\rm max}$ in (\ref{Tmax}) can be expressed as
\begin{equation}
T_{\rm max}\simeq0.6\frac{\alpha^{1/2}m_\phi^{3/4}}{\left(4\pi\right)^{1/4}M_P^{1/2}}\frac{M_P^{1/4}}{g_*^{1/4}}m_\phi^{1/4}M_P^{1/4}=\frac{0.6m_\phi\sqrt{\alpha}}{\left(4\pi\right)^{1/4}g_*^{1/4}}~,
\end{equation}
then the condition $T_1<T_{\rm max}$ translates into
\begin{equation}
\alpha>\left(\frac{(4\pi)^{1/4}g_*^{1/4}}{4*0.6}\right)^{2}\simeq0.58\sqrt{g_*}~.
\end{equation}
These two conditions are contradictory with $g_*\ge1$, concluding that the thermal effects cannot modify the thermal history without triggering a parametric resonance. However, this conclusion relies on the assumptions that $\Lambda=M_P$ and $\varphi_{\rm end}=M_P$. While the former is reasonable, the latter may not be true in small field models in which the inflation ends at a sub-Planckian field value. Besides, if the perturbative reheating is dominated by more than one kind of interactions, the relation between $\Gamma_0$ and $\Gamma_1$ is usually not simply linear. For example if there exist a large Yukawa coupling and a small axion coupling, then the main contribution for $\Gamma_0$ comes from the Yukawa interaction and $\Gamma_1$ from the axion-like coupling since the former is not Bose enhanced. Thus we conclude that the thermal effects can only modify the thermal history if $\varphi_{\rm end}$ is not too large, if one has perturbative reheating in mind. The stepwise evolution of $T$ presented in the above figures still holds schematically, at least within the EFT framework. All in all, these does not inherently invalidate the method used, as long as the non-perturbative dissipation rate of the reheating process can locally be written in the form $\Gamma_n(T/m_\phi)^n$.

\section{Conclusions}\label{conc}
The reheating phase, which populates the universe with hot plasma and sets up the initial conditions for the radiation dominated era, affects the expansion and thermal history of the early universe and may produce relics such as gravitational waves, baryon asymmetry or dark matter. The information of the reheating process, especially the maximum effective temperature during reheating and the reheating temperature at the beginning of the radiation era, is crucial to understand how our universe evolves and how the matters are produced.

In this work, we study the time dependence of the temperature during the reheating process.
The effective dissipation rate $\Gamma_\varphi$ of the inflaton is assumed to be a polynomial containing the vacuum decay of the inflaton represented by $\Gamma_0$, as well as a thermal term $\Gamma_n (T/m_\phi)^n$ which is an arbitrary power monomial of the effective temperature. By the piece-wise approximation of $\Gamma_\varphi$, we then solve the Boltzmann equations for the energy densities of inflaton and radiation, $\rho_\varphi$ and $\rho_\gamma$, and give analytic solutions for the time-dependent temperature. In this way, we describe the thermal history of the reheating process. The maximum temperature in the universe and its dependence on model parameters
are discussed. We also compare the analytic results to the numerical ones, which shows that the piece-wise approximation method is powerful and sufficient for us to study the information of reheating process.

Our work generalizes previous research in which the case $n=2$ was studied. We solve the equations governing the evolution of the energy densities of inflaton and radiation in arbitrary $n$ case, and we give practical examples in which the numerical results match the analytic ones well when $n=1$ or $2$. We find that some expressions in~\cite{Drewes:2014pfa} need to be modified in the $n=2$ case. We also study the situation where there exists a preheating phase after inflation: it changes the initial conditions of reheating for the effective temperature, and therefore the reheating could be dominated by the thermal feedback effects from the beginning. The impact of thermal feedback effects on the expansion history of the universe and therefore on the predication of CMB is discussed.

The result can contribute to the knowledge of the early epochs of cosmic history, and also its thermal relics including dark matter. For example, a large plasma temperature allows the production of heavier particles from inflaton, while this is forbidden via perturbative decays for kinematic reasons. Besides, the abundance of relics that is out of equilibrium before freeze-out is sensitive to the thermal history.

\vspace{6pt}

\acknowledgments
This research was funded in parts by the NSF China (11775110, 11690034), the European Union’s Horizon 2020 research and innovation programme (RISE) under the Marie Sklodowska-Curie grant agreement (644121), and the Priority Academic Program Development for Jiangsu Higher Education Institutions (PAPD). The author would like to thank Marco Drewes, Edna Cheung and Hui Xu for their helpful comments on the manuscript.


\section{Introduction}

Autonomous systems with self-adaptive or learned components have applications in a variety of tasks due to the increasing interest in industrial automation, industry 4.0, or cognitive robotics.
However, for the trustful application of these systems reliability and safety are crucial factors.

Reinforcement Learning (RL) has been successfully applied to control a variety of applications~\cite{Arulkumaran2017}, such as game playing~\cite{Justesen2019}, robotics control~\cite{Kober2013}, or combinatorial optimization~\cite{Dai2017}. %
The RL agent explores its environment by trying random actions and learns through reward signals about its dynamics.
Through sufficient exploration and discovering successful strategies, it converges towards an optimal policy and solves the problem that is defined by the environment and its reward function.
As an RL agent learns from trial-and-error, there is no possibility to employ the agent with hard rules or safety guidelines from the beginning on.
And also after training, there are no given guarantees for the safe and reliable behavior of the agent.
Furthermore, training RL agents is often data-inefficient and takes many iterations of training runs in sometimes computation-intensive simulations or even more costly on actual physical devices.
These properties of RL agents make it challenging to deploy them, either as stand-alone agents or as trainable components within a larger system architecture~\cite{Dulac-Arnold2019}.
Their integration requires handling uncertainties and to ensure they are applied within their designated environment in a more strict way than common software components.
The RL agent is not aware that it is receiving inputs that it is not sufficiently trained for and will be select an action in any case even without having been in the same or similar state before.
An exception to this is Bayesian RL \cite{Ghavamzadeh2015} where uncertainty estimations are considered, but these algorithms have not yet reached the same level as the more common model-free RL algorithms we consider.

Contrary to RL, symbolic AI, such as expert systems or constraint solving, can make decisions based on an explicit constraint model of the environment, but depending on its complexity, it is often infeasible to model and solve the environment for every decision.
Still, constraint models are data-efficient, as there is no training, and always avoid unsafe behaviors as defined by the model.

In this work, we make an argument for the integration of both techniques, such that the RL agent, also the \emph{learner}, is guided by an encapsulating constraint model that describes safety constraints for the agent's task.
We refer to this integration as \emph{constraint-guided reinforcement learning} (CGRL) and it has the goal to design a reliable and safe agent.
Within this hybrid agent, safety models are used to augment the interaction of the RL agent with its environment to avoid unsafe or undesired actions and behaviors.
More specifically, we discuss three interfaces to augment the agent-environment interaction:
1) the interface between environment and agent, where external state information about the environment is observed;
2) the external interface between agent and environment, where selected actions are transferred to the environment;
3) the internal decision-making within the agent, to directly adjust which actions are taken.

Encapsulating RL agents within an outside model for safety or training purposes has been explored before \cite{Garcia2015}, e.g. in the context of shielding reinforcement learning agents~\cite{Alshiekh2018} or imitation learning, where the agent is shown the behavior of experts and tries to mimic them~\cite{Price2003,Ho2016}.
Our focus for the integration within the scope of this work lies on the different interfaces for augmentation and their effects on the system's design and reliability.
We further discuss constraint solving techniques to model acceptable behavior or to define how to act in case of bad decisions from the RL agent.

Developing dedicated encapsulation techniques, represented by the different constraint models on the interfaces, allows setting strict rules and guidelines, while not prohibiting the free exploration and decision making of the agent within the set boundaries.
To underline the relevance and motivate further research in this direction, we discuss the general working mechanisms of reinforcement learning as well as the potential integration points for an external constraint model.
The implementation of CGRL and its effectiveness is shown in two case studies.
Our results confirm that constraint-guidance is beneficial for both faster training and better final performance than pure model-free reinforcement learning.

The main contributions of this paper are:
1) We integrate safety models in the RL process at the intersection between agent and environment to guide the decision-making.
2) We present an architecture to implement these models by using both satisfaction and optimization techniques.
3) We perform an experimental evaluation of the presented method on two case studies: a card game and a grid world with a combinatorial puzzle goal.

The remainder of the paper is structured as follows:
After reviewing the necessary background on RL and constraint solving as well as existing works on their intersection in Sections~\ref{sec:background} and \ref{sec:relwork}, we discuss their integration for constraint-guided RL in-depth in Section \ref{sec:integration_cp_rl}.
In Section~\ref{sec:experiments}, we implement CGRL in two environments and evaluate its effects.
Before concluding the paper in Section~\ref{sec:conclusion}, we discuss the impact and current challenges for CGRL as well as opportunities for future research in Section~\ref{sec:discussion}.

\section{Background}
\label{sec:background}

\subsection{Reinforcement Learning}
\label{sec:rl}

In reinforcement learning (RL) an agent interacts with an environment by observing the environment's state $s_t$ at time step $t$ and selecting an action $a_t$.
From the effect of the action, a reward $r_{t+1}$ is received and a new state $s_{t+1}$ is reached.
The goal of the RL agent is to maximize the expected return, described as the cumulative discounted reward at time $t$, that is, $G_t = \sum_{k=0}^{\inf} \gamma^k r_{t+k+1}$, where $\gamma \in [0, 1)$ is the discount factor.
More formally, the environment is characterized by the tuple $\langle S, A, T, R \rangle$, with $S$ being the set of states in the environment, $A$ being the set of actions.
$T : S \times A \rightarrow S$ is the transition function from the execution of an action in one state leading to another state and $R$ is the reward function, which assigns a numerical feedback to an action in a state $r : S \times A \rightarrow \mathbb{R}$.

The agent follows a trainable policy $\pi$ from which it chooses the action in each state $\pi : S \times A \rightarrow [0, 1]$.
The policy assigns each action $a$ a probability to be taken and from this probability distribution either the action with highest probability is picked greedily or it is randomly sampled to encourage exploration.
This general agent-environment interaction process is shown in Figure \ref{fig:integration}, including extensions for CGRL which we will discuss at a later point.
The agent interacts with the environment in finite episodes, which means that some states are terminal states and once these states are reached, the episode is over and a new episode starts. Terminal states can represent both desired goal states, e.g. in a maze where a certain location has to be reached, but also undesired states, e.g. the agent performs illegal or unsafe behavior.
For a detailed introduction to RL, we refer the interested reader to the book by Sutton and Barto (2018).

Over the recent years, deep RL, where the action selection policy is approximated by deep neural networks, has received large attention and many RL algorithms have been proposed in the literature \cite{Arulkumaran2017}, for example, Proximal Policy Optimization (PPO) \cite{Schulman2017} which we use throughout the experiments in this paper.
These algorithms share, that they represent the state as input to the neural network and receive a vector of action scores as output.
These scores are interpreted as probabilities and the action is randomly sampled or greedily chosen by the highest score.
The state representation is usually problem-specific.
For problems where the state resembles a form of spatial information, e.g. images, convolutional neural networks (CNN) are commonly used, whereas otherwise fully-connected feed-forward neural networks are common.

We use PPO in an actor-critic way~\cite{Mnih2016}, where the agent is modeled by two neural networks.
The first network, the actor, outputs the scores for each action as described above.
The second network has only a single output that predicts the expected reward for the current state.
Both networks usually share the same architecture and also the weights of the initial layers with the only difference being the final layers, where the actor has as many outputs as actions, the critic only has a single real-valued output.
Using a critic network helps to stabilize and accelerate the learning process since it allows to compare the current policy against the earlier made experiences and supports the improvement of the policy.

\subsection{Constraint Solving}

Throughout this paper, we refer to constraint solving in the context of modeling the acceptable behavior of a learner, like the RL agent, which is either a satisfaction problem, when asking "Does the learner behave according to the specification?", or an optimization problem when asking "What is the state representation or action that complies with the specification with only minimal changes?".
Modeling a satisfaction problem allows to identify undesired states and actions, but does not guide the learner, except when designing each individual action as separate problems.
The optimization model instead allows augmenting the state or action towards specification compliance, according to its constraints.

Formally, we define a \emph{constraint satisfaction problem} (CSP) as a triple $\langle \mathcal{X},\,\mathcal{C},\, \mathcal{D} \rangle$ where $\mathcal{X} = \{\mathbf{x}_1, \dots, \mathbf{x}_n\}$ is a set of variables, $\mathcal{C}$ is a set of constraints $\mathcal{C}=\{c_1,\dots,c_m\}$, where each constraint $c_i$ involves a subset of the variables.
Each variable $\mathbf{x}_i$ is also associated with a finite integer domain $\mathcal{D}(\mathbf{x}_i)$ of all its possible values.
A variable assignment $\theta$ satisfies the CSP $\langle \mathcal{X},\,\mathcal{C},\, \mathcal{D} \rangle$, if each variable $x$ has a value assigned from its domain $D(\mathbf{x}_i)$, such that none of the constraints $\mathcal{C}$ are violated.
A \emph{constraint optimization problem} (COP) extends the CSP formulation by an objective function $f_z$, which is calculated over a subset of the variables, with an objective value $z$, which, in our context, is to be minimized, while the variable assignment $\theta$ still satisfies the constraints $\mathcal{C}$.
For a detailed overview on constraint solving, we refer the interested reader to \cite{Rossi2006}.

The constraint model can be modelled and solved by dedicated constraint solving libraries, e.g. MiniZinc or Essence. %
This has the advantage of having solver-provided guarantees for the correctness of results - under the assumption of a correct model - and tools optimized for the task, but with the disadvantage of increasing system complexity due to additional libraries and related dependencies.
The alternative would therefore be to implement the constraint model within the programming paradigm of the autonomous system.
In this case, there are no additional external dependencies, but the implementation is potentially more error-prone since constraint modeling languages are usually more expressive for the tasks than general-purpose programming languages.
We will come back to the discussion between these two alternatives in the discussion part of the paper in Section~\ref{sec:discussion}.

\section{Related Work}
\label{sec:relwork}

The main areas of related work are safe, ethical, or social reinforcement learning and the design of autonomous systems.
The work closest to ours is the method of \emph{Shielded Reinforcement Learning} by Alshiekh~\emph{et al.}\xspace~\cite{Alshiekh2018} .
The authors propose a safety framework where a shielding component is deployed either between the environment and the agent or after the agent has made its action selection.
The shield represents a safety specification, based on temporal logic, that assures the correctness of the system.
This is related to our method, which also considers integrating a constraint-guidance component at the same interfaces.

Similar works have been based on using barrier functions~\cite{Cheng2019} or provable formal theories of the acceptable behavior specifications of the RL agent~\cite{Fulton2018}, or focus on specifics of the RL environment, e.g. stochastic dynamics in robotics~\cite{Li2020a}.
These approaches focus mostly on the provable aspects of the formalization, and strongly motivate the necessity and feasibility of hybrid approaches, but focus less on the aspects of how they are implemented and expressed.
A further avenue is the design of RL methods that have safety inherently in their functional methodology, e.g. \cite{Achiam2017,Chow2018}.
Unlike these methods, we provide safety extensions to any form of RL agent, even without inherent safety notion.
We see an advantage in a modular approach, however, we do not argue that one or the other is necessarily better.
Finally, for an in-depth overview on safe reinforcement learning, we refer to the survey by García and Fernández~\cite{Garcia2015}.

Understanding and steering how agents behave and learn in unknown environments also affects ethical concerns and is an aspect of active research to include ethical standards as boundaries for the agents' behavior~\cite{Rossi2019}.
For example, Noothigattu~\emph{et al.}\xspace explore the effects of imposing ethical constraints on a Pacman agent, such that it is prohibited to eat ghosts~\cite{Noothigattu2018}.
The consideration of ethics for learning agents is of relevance and imposing constraint-guidance, based on models that capture ethical constraints, proposes an approach towards this challenge.
Another area of interest is to let the agent learn which actions are either not safe, not desired, or just unfeasible in a given state.
Zahavy~\emph{et al.}\xspace~\cite{Zahavy2018} propose a method where the agent learns which actions not to choose in a given state.
Learning when an action is infeasible and should be avoided or when it is just not available has been part of teaching an agent's policy, where the agent learns which actions to avoid, but without an explicit specification of wrong and forbidden actions.

\section{Integrating Constraints and RL}
\label{sec:integration_cp_rl}

To achieve constraint-guided reinforcement learning, the constraint model needs to be integrated within the agent-environment process.
In this section, we identify and discuss the three interfaces for constraint-guidance (Figure~\ref{fig:integration}) and their impacts on the functionality and the autonomy of the encapsulated RL agent (Figure~\ref{fig:space_reduction}).

\begin{figure}[t]
    \centering
    \resizebox{0.9\columnwidth}{!}{\input{figures/agent_environment_process.tex}}
    \caption{The Agent-Environment Process with constraint-guidance integration: (1) Environment-Agent, (2) External Agent-Environment, (3) Internal Agent-Environment.}
    \label{fig:integration}
\end{figure}

\begin{figure*}[t]
    \centering
    \subfloat[Environment-Agent]{%
        \includegraphics[width=0.38\textwidth]{figures/observation_space_venn.pdf}}\hfill
    \subfloat[External Agent-Environment]{%
        \includegraphics[width=0.3\textwidth]{figures/action_replacement.pdf}}\hfill
    \subfloat[Internal Agent-Environment]{%
        \includegraphics[width=0.3\textwidth]{figures/action_mask.pdf}}
    \caption{Effects of the augmentation interfaces on the observation and action space. The Environment-Agent interface augments the observation space to map unseen or unsafe observations to a familiar environment in which the agent can safely act. The Agent-Environment interfaces replace unsafe actions by a safe default (External) or limit the availability of actions for the agent before selection (Internal). The safety of an action depends on the observation and can vary from step to step.}
    \label{fig:space_reduction}
\end{figure*}

One integration point lies at the interface between the environment and the agent, where observations and rewards are received.
The second integration point addresses the interaction from the agent to the environment, i.e., when an action is issued by the agent.
These two interfaces wrap the RL agent as a black-box and manipulate the inputs and outputs of the agent.
The third integration point, on the other hand, lies within the agent and allows to apply constraint-guidance in the internal mechanisms of the RL agent.
Relating to the existing literature, the first two interfaces are similar to the work by Alshiekh \emph{et al.}\xspace \cite{Alshiekh2018} when their proposed shielding is introduced either before or after the action selection part.
The constraint models for each of the interfaces are modeled based on a risk- and safety-analysis of the system.
The approach is to formalize the most essential boundaries for the agent's behavior without restricting the agent's freedom within these boundaries more than necessary.
We discuss each of the three integration types in the following and consider the interfaces that have to be implemented by a constraint-guidance mechanism.

\subsection{Environment-Agent Interface}

The \textit{Environment-Agent Interface} augments the transfer of state information from the environment to the agent in a form of state preprocessing.
By augmenting the state information, either additional information can be added to the state, or parts of the state can be masked or modified to avoid unsafe decisions.
One application of this is to reduce the decision space for the RL agent.
Aspects of the state, which would lead the agent to act in a specific, but potentially dangerous way, are replaced by safe dummy values or a representation that does not lead to a certain action, but also does not misrepresent the actual state space.
For example, when having an agent with an inventory of items, those items which are dangerous to be used in the current state can be omitted from the augmented state that is given to the agent.

\begin{definition}[Environment-Agent Interface]
The interface function $f : S \rightarrow S$ augments the observed state, such that the difference diff($s_t, s'_t$) is minimized while following side constraints $\mathcal{C}_{Env}$, and returns a modified state $s'_t$.
\end{definition}

Technically, the augmentation component is based on a constraint model $c(s_t) = s'_t$, that specifies potentially dangerous states, and can make modifications
$c(s)$ receives a state $s_t$ and returns a state $s'_t$, which can be identical to the original state if $s_t$ is not considered dangerous, or modified to exclude or mask the undesired parts of the state.

To implement this interface in an \emph{Observation Mask} model, that augments the observed state for the agent, a constraint optimization model can be used.
The instance parameters are the observed state, the decision variables are the augmented state, and the objective function measures the changes between the observed and augmented state.
During optimization, the augmented state is adjusted such that the state complies with the modeled specification, e.g. unusable items from an agent's inventory are masked while making as few modifications as possible.
Making only minimal changes is relevant to keep the augmented state close to the actual observation because large changes can lead the agent to make an unsuitable decision for the actual state and thereby hinder the success within the environment.

\subsection{External Agent-Environment Interface}

Once the agent selected an action, this action is passed back to the environment to be executed.
At this point, the \emph{Agent-Environment Interface} augments the selection and can interfere by observing both the most recent state and the selected action in order to decide whether the chosen action is acceptable or violates some guidance criterion, e.g. being unsafe or potentially harmful.
In case the action is not accepted, the interface can replace the selected action following a heuristic procedure.
This replacement action might not be optimal as it does not follow the learned strategy but can avoid unsafe outcomes instead.

\begin{definition}[External Agent-Environment Interface]
The interface function $f : S \times A \rightarrow A$ receives the last observed state $s_t$ and the taken action $a_t$ as input and returns a, probably modified, action $a'_t$.
Action $a'_t$ must obey the constraints $\mathcal{C}_{Env}$ for the environment and the last observed state.
\end{definition}

An \emph{Action Replacement} model can be implemented either as a satisfaction model, that validates whether an action follows the specification and replaces invalid actions by a default action or as an optimization model that has a notion of distance between actions and can optimize the action change to be minimal while following the specification.
The Agent-Environment Interface forms thereby a fallback strategy that monitors the RL agent's behavior and only if it does not behave according to the specified acceptable behavior, its decisions are overruled.

\subsection{Internal Agent-Environment Interface}

The first two augmentation interfaces are decoupled from the RL agent and see it as a black-box.
The third interface, \emph{Internal Agent-Environment Interface}, is integrated within the action selection mechanism of the RL agent.

While augmenting the observation and checking the selected action against a constraint model allow to prohibit undesired or unsafe behavior, the RL agent is not aware of this constraint-guidance and does not learn that its action was overruled by an external mechanism.
That means the agent does not learn to avoid this action the next time the same state is encountered and the observed reward does not match the selected action.
The internal integration into action selection overcomes this, but at the cost of requiring access inside the agent, unlike the previous interfaces.

\begin{definition}[Internal Agent-Environment Interface]
The interface function $f : S \rightarrow \mathbb{R}^{n}$ where $n$ is the number of available actions: $f(s_t) = m$.
Each action is assigned a masking score, which is inserted into the policy $\pi$ and multiplied with the network outputs before selecting an action.
From the masked output, the action selection is continued with a bias towards actions favored by $f$.
\end{definition}

Using constraint-guidance, the scores of the individual actions can be adjusted to either discourage unwanted actions by decreasing their score or to encourage the selection of desired actions by increasing their score.
In the simplest case, $m$ is a boolean vector where each action is masked as allowed/available or unavailable for the current state.
Any of these adjustments steer the RL agent towards selecting desired actions for which then a reward is received and the policy is trained with the normal training procedure.

An \emph{Action Mask} model, that implements the internal agent-environment interface, can be interpreted as a constraint satisfaction model with one decision variable per action, where a solution enables all allowed actions for the given state.

\subsection{Deriving Constraint Models}
\label{sec:deriving_constraints}

All interfaces require the existence of a constraint model for augmentation.
There are different ways to derive these constraint models and how to formulate them.
The first option, which we consider mostly in this paper and for the experimental evaluation, is to have the constraint model handcrafted from expert knowledge.
This knowledge can either stem from precisely analyzing error causes and avoiding these situations or setup generic rules for unsafe states and actions that are known to cause unrecoverable error situations.

The second option targets to iteratively refine the constraint model through a dedicated learning process, for example using constraint acquisition or learning~\cite{Bessiere2015}.
This process would require access to a safe simulation environment as the initial set of rules cannot be expected to be safe.
Over subsequent iterations, the RL agent can act using the initial constraint model, and, in case of unsafe situations and error cases, the provoking states and actions are used to refine the constraint model via a constraint learning step.
We do not further consider this option in this work, but leave it as a promising direction for future work on CGRL.

\section{Experiments}
\label{sec:experiments}
In the following section, we apply constraint-guided reinforcement learning in two environments.
The structure for each case study environment follows the same approach.
At first, we introduce the problem setting, the available state information for the RL agent, and the reward structure.
Afterwards, we discuss the constraint-guidance models to integrate with the RL agent.
In each environment, we consider three constraint-guidance models, one for each of the integrations previously discussed, and a vanilla RL agent without constraint-guidance as a baseline for comparison.

\begin{figure}[t]
    \centering
    \includegraphics[width=0.9\columnwidth]{figures/architecture.pdf}
    \caption{Architecture: The constraint-guidance model acts as a wrapper and has first access to observations and final control over the selected actions.}
    \label{fig:architecture}
\end{figure}

The RL algorithm is PPO (Proximal Policy Optimization) \cite{Schulman2017}\footnote{The PPO implementation is available here: \url{github.com/lcswillems/torch-ac}}.
All runs use the same hyperparameters, which closely follow the defaults proposed by the authors, except for the network architecture, which is adjusted for each environment.

\subsection{Environment I: Card Game}
\label{sec:thegame}

\subsubsection{Scenario \& Rules}
In the first experiment, we focus on a card game, modeled after ``The
Game''\footnote{\textit{The Game: Are you ready to play The Game?} (2015). Created by Steffen Benndorf. Published by N{\"u}rnberger-Spielkarten-Verlag.}.
The goal of the game is to distribute 98 unique cards, ranging from 2 -- 99, onto four card piles.
Two of these card piles are to be filled in increasing order from 1 -- 100, the other two in decreasing order from 100 -- 1.
It is only allowed to place cards with a higher respectively lower card value on each of the piles, otherwise the action is invalid.
There is one exception to this rule.
A lower respectively higher card value can be played if the difference in card value is exactly 10, i.e. if the increasing pile has card 42 on top, card 32 is allowed on top of it.

The game ends when all cards have been distributed or there is no further move possible.
The original card game can be played with 1--5 players, in this case study we focus on the single player version, where the player has 8 cards on their hand and refills their hand from the stack after having played two cards.

\subsubsection{Environment \& Agent}

The state is represented by a 12-di\-men\-sio\-nal vector, consisting of the four card piles and the eight hand cards, sorted in increasing order.
If a slot for a hand card is empty, because a card was played and no new cards were drawn yet, it is replaced by a placeholder value of $-1$.
The agent can place any hand card on any of the card piles, leading to a total of $8 \times 4 = 32$ possible actions.
The agent's policy is represented by a feed-forward neural network with three fully connected layers of each 32 hidden nodes with tanh activations.
Both the actor and the critic network share the weights of the initial layers.

The agent receives a positive reward of $r = 0.1$ for each valid move it plays and a negative reward of $r = -0.1$ for each invalid move.
At the end of the round the final reward equals the number of played cards, i.e. at most 98.

\subsubsection{Constraint-Guidance Models}

We define the CGRL models such that invalid actions are avoided, but do not impose are restrictions on the inner RL agent.
If an action has to be selected through constraint-guidance, it greedily picks that action which has the smallest difference between the card piles and the hand card value.

\textbf{Observation Mask}
The Environment-Agent Interface is implemented such that cards with only negative effects, i.e. they cannot be played on any of the card piles, are hidden from the RL agent by showing that these card slots are empty.
If all cards can be placed on at least one of the piles, no changes are made.

\textbf{Action Replacement}
The External Agent-Environment Interface monitors the selected action for invalid card placements, which would result in a negative reward penalty. These actions are heuristically replaced by that action which incurs the smallest step on any of the card piles, i.e. that card which the minimum difference between card and pile value.

\textbf{Action Mask}
The internal environment-agent interface behaves similar to action replacement, but instead of replacing the action after it has been selected, it masks those actions which are invalid.
This enforces the agent to always select a valid action and focuses the exploration towards understanding good game strategies without having to learn about the validity of actions.
Further variants can be considered, for example, it is possible to enforce playing cards where the special rule of decreasing/increasing the card pile in the "wrong" direction can be used.
But, this goes beyond the usage of constraint-guidance for safety aspects and enforces a particular playing strategy upon the agent, which would be less comparable to the other models.

\subsubsection{Results}

The reward curves of the four different agents over a training time of 50 million frames are shown in Figure~\ref{fig:thegame_training}.
We train the RL agent such that it runs multiple environments in parallel to collect experiences for training.
This rollout both allows to run the training in parallel on multiple CPU cores as well as it stabilizes it, because more diverse states and rewards are available.
We show both the best and the average rewards of 64 training environments.

\begin{figure}[t]
    \centering
    \resizebox{0.8\columnwidth}{!}{\input{figures/rewards_TheGame-v0.pgf}}
    \caption{Card Game: Maximum (solid lines) and average (dotted) rewards with and without masking. We show the mean and std. dev. over 10 training runs.}
    \label{fig:thegame_training}
\end{figure}

The results show a strong impact of constraint-guidance over the RL-only agent, but they also show differences between the different integration techniques.
The best performance is reached with Action Mask model, which implements the internal agent-environment interface, and manipulates the available actions by removing invalid actions.
Its performance surpasses that of the other agents from the beginning and quickly reaches a high level.

The other models follow a similar learning curve but are shifted in their base performance, which we attribute to the effects of the different augmentation models.
Augmenting the observations, i.e. removing hand cards that are not possible to be played in the current state, is the least beneficial mask in this scenario, and also causes instability in the training process visible by the high fluctuation at around 49 million frames.
This is caused by the game mechanics, where only in few cases a card is actually invalid to be played on all four card piles. It is much more common that a card is only invalid for a few piles, but this cannot be modeled in our setup by augmenting only the observation, but requires adjusting the actions, like the Action Mask does it.
Still, it can be helpful and performs better than RL only.

Regarding the training duration until the agent's performance converges and does no longer improve by a large margin, we observe faster convergence when using the action mask model compared to the external augmentation.
The RL only and the Observation Mask models steadily improve over a longer time, but their main performance level is reached in a similar time than the others.

\subsection{Environment II: Grid World Puzzle}
\subsubsection{Scenario \& Rules}

The second environment is a puzzle variant of the grid world problem.
The agent is placed in a random 2D grid world environment, consisting of four rooms, and has the goal to collect a number of distinct items before finding the exit door of the grid world (see Figure~\ref{fig:grid world_environment}).
It can move forward, turn left or right, or pickup an item that is in front of the agent.
Every item is described by its shape (key or circle) and its color.
The goal is to pick one item of each type and reach the exit door in as few steps as possible.

\begin{figure}[t]
    \centering
    \includegraphics[width=0.22\textwidth]{figures/gridworld_environment.pdf}
    \caption{Grid World Puzzle: The agent (red arrow with observable space highlighted) has to collect 16 items (circles and keys in 8 colors) before reaching the goal (green square).}
    \label{fig:grid world_environment}
\end{figure}

This environment poses two challenges for the agent:
1) The visible state of the grid world is limited to the eight surrounding tiles next to the agent's position, which requires exploration of the environment;
2) the agent needs to learn to avoid items that it has already collected.
The inventory of already collected items is part of the agent's state and given to the agent at each step.

\subsubsection{Environment \& Agent}

The total grid world environment is $32 \times 32$ tiles, divided into four equally sized rooms connected with open doors.
Both the placement of doors and the items is randomly generated for each new environment.
The state of the environment is represented by the agent's field of view, $7 \times 7$ fields, and its current inventory, listing all collected items.
Grid world tiles can be either a plain field, a field with an item, or the exit door, all of which can be accessed by the agent at any time, or a wall, which blocks the agent and cannot be accessed.
The set of actions is $a \in \{\text{forward}, \text{left}, \text{right}, \text{pickup}\}$.
The agent's policy is represented by a small convolutional neural network with three 2D convolutional layers followed by two fully connected layers.
Each layer is connected with ReLU activation functions.
Again, both the actor and the critic network share the weights of the initial layers.

The reward function is modeled such that the agent receives a positive reward of $1/16$ for each distinct object it collects.
Collecting a duplicate item, that is already in the inventory, causes a negative reward of $-1$ and ends the episode.
The maximum length of an episode are $8192$ steps, once this number of steps is reached, the episode is also terminated.
When the exit door is reached, the agent receives an additional reward of $1 - 0.8 \times \text{Number of steps}/8192$.
This reward function does only return sparse rewards at specific events instead of continuously giving feedback to the agent, which is contrary to the card game environment, where every action of the agent is rewarded by the environment.
The implementation of this environment is based on the \emph{Mini\-Grid} project \cite{gym_minigrid} with extensions for the additional puzzle objective.

\subsubsection{Constraint-Guidance Models}

The goal of the constraint-guidance models for the grid world puzzle environment is to avoid collecting duplicate items, i.e. the undesired behavior.%

\textbf{Observation Mask}
The model for augmenting the observed states from the environment hides any duplicate item from the observed state and masks it as a wall.
This signals to the agent, that the masked tile is not of interest and encourages the agent to further explore or focus on another, unmasked item.

\textbf{Action Replacement}
When the agent decides to move to a tile with a duplicate item, this model observes this specification violation and replaces the move by turning right.

\textbf{Action Mask}
To avoid picking duplicate items, the action mask disables the action to pick an item as long there is no distinct new item in front of the agent.
Thereby the agent can either move, as long as there is no new item in front of it, or it is set to pick a new item, which disables the movement.
This guidance model removes the logic for picking items from the learned agent, which is then only concerned with navigation towards distinct items.

\subsubsection{Results}

The results, in form of training curves, are shown in Figure~\ref{fig:gridworld_training}.
As in the previous experiment, we trained the agents with 64 environments and parallel and report both the best and average performance over time, repeated with 10 different random seeds.

\begin{figure}[t]
    \centering
    \resizebox{0.8\columnwidth}{!}{\input{figures/rewards_MiniGrid-Combination-Picker-32x32-v0.pgf}}
    \caption{Grid World Puzzle: Maximum (solid lines) and average (dotted) rewards with and without masking. Shown is the mean and standard deviation over 10 training runs.}
    \label{fig:gridworld_training}
\end{figure}

Models using constraint-guidance perform better than the RL-only agent, with the average CGRL performance being competitive to the best rewards of the RL-only.
This results show the general difficulty of the environment for a standard RL algorithm.
The average reward is close to the minimum reward that is received for successfully reaching the exit door, but does not show signs for increased attention towards successfully collecting unique items.
The CGRL models on the other hand achieve higher rewards, with the Observation Mask model, that excludes already collected items from the observed state, performing best before the Action Mask and Action Replace models.
We further observe a faster performance increase in the Action Mask model than in the eventually better Observation Mask model.
While the observation mask excludes already collected items from the observed state by presenting them as walls, the agent still can walk towards them and collecting them, as the selected action is not augmented. As the agent needs to learn that walls cannot be accessed, the observation mask is not as effective as the action mask in the beginning.
The action mask, which was most effective in the card game, disallows the agent to access fields with already collected items, but only when the agent stands directly in front of the item. Before that state is reached, the agent can still plan to collect the item, if it has not learned to avoid duplicates. This can cause the lower performance of this mask.

\section{Discussion}
\label{sec:discussion}

Both performed experiments showed the effectiveness of CGRL for guiding a RL agent in complex environments.
This result confirms our expectations, since we introduce external domain knowledge into the training process and the agent's behavior.
However, we also observe a better data-efficiency of the training and the avoidance of unsafe behavior, as confirmed by higher rewards.
Better data-efficiency is important when a single step of the RL agent is more costly than in the small and computationally cheap environments used here, e.g. when running a complex physical simulation or even using physical components, there it is especially useful to model the acceptable behavior via a constraint model.
In addition to the physical environment, where unsafe behavior might be harmful for the machine or the physical environment, the avoidance of unsafe behavior is crucial for the application and deployment to either safety-critical or ethically challenging environments.

A challenge in the application of CGRL can occur from the additional runtime overhead added from the constraint solving for augmentation.
As the model has to be run for every step of the RL agent a small overhead can accumulate to a longer runtime.
However, this limitation mostly applies only in cases where the overhead is significant in relation to the general computational cost of one step in the environment.
In an environment where every step is already costly, the benefits of faster training and less undesired behavior will exceed the disadvantage of a runtime overhead.

As future work on CGRL, besides the challenges raised, we plan to extend the method to environments with more complex state representations.
Here it can be necessary to additionally preprocess the observation before constraint reasoning can be applied.
Another direction is to adopt principles from imitation learning, where expert behavior can be observed, to learn augmentation rules, e.g. via constraint acquisition \cite{Bessiere2015}, as addressed in Section~\ref{sec:deriving_constraints}.
This way, CGRL would first learn a static constraint model for guidance and then allow an RL agent to learn a strong policy within its specification.

\section{Conclusion}
\label{sec:conclusion}

We presented Constraint-Guided Reinforcement Learning (CGRL), a method to combine constraint solving and reinforcement learning into a hybrid agent that behaves according to a specification, but is otherwise free to explore and learn the best policy for the problem.
CGRL augments the interfaces between the RL agent and the environment through constraint models, which specify acceptable behavior that the agent should comply to.
Constraint solving acts here as a safety fallback for the encapsulated RL agent.

Augmentation can happen by changing the state observation from the environment to exclude certain aspects in order to encourage the agent towards accepted behavior, but it can also observer the actions made by the agent and adjust to avoid behavior that violates the specification.
Finally, a third potential interface for constraint-guidance is to augment the action selection process within the agent -- unlike the first two approaches, which consider the RL agent as black-box -- and mask actions that would violate the specification.

\printbibliography
\end{document}