Reheating in two-sector cosmology

TL;DR

This paper investigates the dynamics of reheating in two-sector cosmology, focusing on quantum statistical effects and inflaton-mediated energy transfer, providing new power laws and analytic estimates for temperature evolution and asymmetry.

Contribution

It introduces comprehensive numerical solutions including quantum statistics and develops new analytical estimates for temperature asymmetries in two-sector reheating scenarios.

Findings

Quantum statistical effects significantly influence temperature evolution.

Inflaton-mediated scattering is dominant at specific temperature regimes.

Classical reheating approximates asymmetry well at low temperatures.

Abstract

We analyze reheating scenarios where a hidden sector is populated during reheating along with the sector containing the Standard Model. We numerically solve the Boltzmann equations describing perturbative reheating of the two sectors, including the full dependence on quantum statistics, and study how quantum statistical effects during reheating as well as the non-equilibrium inflaton-mediated energy transfer between the two sectors affects the temperature evolution of the two radiation baths. We obtain new power laws describing the temperature evolution of fermions and bosons when quantum statistics are important during reheating. We show that inflaton-mediated scattering is generically most important at radiation temperatures $T\sim M_\phi/4$, and build on this observation to obtain analytic estimates for the temperature asymmetry produced by asymmetric reheating. We find that for reheating temperatures $T_{\textrm{rh}} \ll M_{\phi}/4$, classical perturbative reheating provides an excellent approximation to the final temperature asymmetry, while for $T_{\textrm{rh}}\gg M_{\phi}/4$, inflaton-mediated scattering dominates the population of the colder sector and thus the final temperature asymmetry. We additionally present new techniques to calculate energy transfer rates between two relativistic species at different temperatures.

Equations

$$\frac{\text{d}{\rho }_{\varphi }}{\text{d}t}+3H{\rho }_{\varphi }$$

$$\frac{\text{d}\rho }{\text{d}t}+4H\rho$$

$$H=\frac{1}{\sqrt{3}{M}_{\text{Pl}}}\sqrt{\rho +{\rho }_{\varphi }}.$$

$$4H\rho =\frac{4}{4+\frac{q}{1-p}}\Gamma {\rho }_{\varphi },$$

$$q(a)=\frac{\partial \ln (\Gamma {\rho }_{\varphi }/H)}{\partial \ln a},\text{and}$$

$$\displaystyle\rho=\Big{(}\frac{2}{5}\Gamma_{0}M_{\textrm{Pl}}\sqrt{\rho_{\phi,I}}\Big{)}\left(\frac{a}{a_{I}}\right)^{-3/2},$$

$$\frac{\text{d}}{\text{d}a}\left(\rho a^4\right)$$

$${\rho }_{\text{max}}=0.24M_{\text{Pl}}{\Gamma }_0\sqrt{{\rho }_{\varphi ,I}}.$$

$$\rho =\frac{{\pi }^2{g}_{\ast }}{30}{T}^4\equiv \alpha T^4,$$

$$\Gamma (T)={\Gamma }_0\frac{\exp (\frac{M_{\varphi }}{2T})\pm 1}{\exp (\frac{M_{\varphi }}{2T})\mp 1},$$

$$\displaystyle T=M_{\phi}\Big{(}\frac{2\sqrt{3}}{\alpha}\frac{{M}_{\textrm{Pl}}}{M_{\phi}^{4}}{\Gamma}_{0}\sqrt{{\rho}_{\phi,I}}\Big{)}^{1/3}\left(\frac{a}{a_{I}}\right)^{-1/2}.$$

$$\displaystyle T(a)=M_{\phi}\Bigg{(}\frac{2\sqrt{3}}{\alpha}\frac{{M}_{\textrm{Pl}}}{M_{\phi}^{4}}{\Gamma}_{0}\sqrt{{\rho}_{\phi,I}}\bigg{[}\bigg{(}\frac{a}{a_{I}}\bigg{)}^{-3/2}-\bigg{(}\frac{a}{a_{I}}\bigg{)}^{-3}\bigg{]}\Bigg{)}^{1/3}.$$

$$\displaystyle T_{\textrm{max}}=0.95M_{\phi}\Big{(}\frac{1}{\alpha}\frac{{M}_{\textrm{Pl}}}{M_{\phi}^{4}}{\Gamma}_{0}\sqrt{{\rho}_{\phi,I}}\Big{)}^{1/3},$$

$$\displaystyle T=M_{\phi}\Big{(}\frac{5\sqrt{3}}{56\alpha}\frac{{M}_{\textrm{Pl}}}{M_{\phi}^{4}}{\Gamma}_{0}\sqrt{{\rho}_{\phi,I}}\Big{)}^{1/5}\left(\frac{a}{a_{I}}\right)^{-3/10}.$$

$$\displaystyle T(a)=M_{\phi}\Bigg{(}\frac{5\sqrt{3}}{56\alpha}\frac{{M}_{\textrm{Pl}}}{M_{\phi}^{4}}{\Gamma}_{0}\sqrt{{\rho}_{\phi,I}}\bigg{[}\bigg{(}\frac{a}{a_{I}}\bigg{)}^{-3/2}-\bigg{(}\frac{a}{a_{I}}\bigg{)}^{-5}\bigg{]}\Bigg{)}^{1/5}.$$

$$\displaystyle T_{\textrm{max}}=0.58M_{\phi}\Big{(}\frac{1}{\alpha}\frac{{M}_{\textrm{Pl}}}{M_{\phi}^{4}}{\Gamma}_{0}\sqrt{{\rho}_{\phi,I}}\Big{)}^{1/5},$$

$$\displaystyle T_{\textrm{rh}}=\begin{cases}\dfrac{4\sqrt{3}}{\sqrt{\alpha}}\dfrac{{M}_{\textrm{Pl}}\Gamma_{0}}{M_{\phi}}&\textrm{boson }\\ \bigg{(}\dfrac{\sqrt{3}}{4\sqrt{\alpha}}M_{\textrm{Pl}}\Gamma_{0}M_{\phi}\bigg{)}^{1/3}&\textrm{fermion }\end{cases},$$

$$\displaystyle T_{\textrm{rh}}=\bigg{(}\frac{\sqrt{3}}{\sqrt{\alpha}}{M}_{\textrm{Pl}}\Gamma_{0}\bigg{)}^{1/2},$$

$$\frac{\text{d}{\rho }_{\varphi }}{\text{d}t}+3H{\rho }_{\varphi }$$

$$\frac{\text{d}{\rho }_a}{\text{d}t}+4H{\rho }_a$$

$$\frac{\text{d}{\rho }_b}{\text{d}t}+4H{\rho }_b$$

$$H=\frac{1}{\sqrt{3}{M}_{\text{Pl}}}\sqrt{{\rho }_a+{\rho }_b+{\rho }_{\varphi }},$$

$${\alpha }_a=\frac{{\pi }^2{g}_{\ast ,a}}{30}={\alpha }_b=\frac{{\pi }^2{g}_{\ast ,b}}{30}=30,$$

$${\rho }_{\varphi ,I}=\frac{1}{2}{M}_{\varphi }^2⟨{\varphi }_I^2⟩$$

$$x=T_b/T_a$$

$$x_{\text{rh}}=\frac{T_b(a^{\prime })}{T_a(a^{\prime })},$$

$$w={\Gamma }_{0,b}/{\Gamma }_{0,a}$$

$$\frac{d{\rho }_{\varphi }}{dt}+3H{\rho }_{\varphi }$$

$$\frac{\text{d}{\rho }_a}{\text{d}t}+4H{\rho }_a$$

$$\frac{\text{d}{\rho }_b}{\text{d}t}+4H{\rho }_b$$

$${\mathcal{C}}_E=$$

$$\times [f_1(p_1)f_2(p_2)(1\pm f_3(p_3))(1\pm f_4(p_4))-f_3(p_3)f_4(p_4)(1\pm f_1(p_1))(1\pm f_2(p_2))]$$

$$\equiv$$

$${\mathcal{L}}_{\text{int}}=\frac{1}{2}{\mu }_a\varphi {\chi }_a^2+\frac{1}{2}{\mu }_b\varphi {\chi }_b^2.$$

Full text

1 Introduction

A host of observational data from galactic rotation curves to the cosmic microwave background (CMB) points towards an unknown form of matter comprising more than 85% of the matter quota in our Universe Ade:2015xua ; Aghanim:2018eyx . While this dark matter (DM) could be a single particle, or a family of particles missing from the standard model of particle physics (SM), to date no effects have been detected that require dark matter interact non-gravitationally with the SM. Further, the traditional weakly interacting massive particle (WIMP) paradigm is under ever-increasing pressure from the dearth of observational signatures in collider, direct, and indirect detection experiments Escudero:2016gzx ; Alexander:2016aln .

A straightforward alternative to the WIMP scenario is that the dark matter is in a ‘hidden sector’ containing new particles and forces that may or may not interact non-gravitationally with the standard model and its extensions. Hidden sector DM refers to a class of models where the DM resides in an internally thermalized hidden sector, and has a relic abundance that is determined by physics unrelated to its couplings to the visible sector. When the hidden sector is out of equilibrium with the SM in the early universe, a wide variety of possible cosmic histories for dark matter are newly possible, together with avenues for discovery distinct from standard WIMP search strategies. Exotic thermal evolution in a decoupled dark sector can lead to substantial changes in DM properties relative to a traditional WIMP (e.g., Carlson:1992fn ; Feng:2008mu ; Sigurdson:2009uz ; Cheung:2010gj ; Pappadopulo:2016pkp ; Dror:2016rxc ; Berlin:2016gtr ; Dror:2017gjq ; Georg:2019jld ; Faraggi:2000pv ; Forestell:2016qhc ; Dienes:2016vei ; Forestell:2018dnu ; Kang:2019izi ; Blanco:2019eij ). Decoupled dark sectors can also induce early departures from the standard radiation-dominated evolution of the universe Zhang:2015era ; Berlin:2016gtr , or admit substantial amounts of relic dark radiation without violating the stringent bounds from Planck, $\Delta N_{\rm eff}<0.36$ at 95% confidence Aghanim:2018eyx .

There is now strong evidence from observations of the fluctuations in the CMB that the thermal era was preceded by an epoch of early accelerated expansion—inflation. Inflation exponentially dilutes any pre-existing matter and radiation leaving the Universe cold and empty. The population of otherwise decoupled sectors cannot therefore be put in ‘by hand’ as an initial condition. Instead it must be generated dynamically in the post-inflationary evolution of the Universe. In the simplest scenarios, the accelerated expansion is driven by a single fundamental scalar degree of freedom, whose weak couplings to matter reheat the Universe via perturbative decays. One of the simplest mechanisms for populating hidden sectors is to couple them to the inflaton so that they are populated at reheating along with the visible sector. By arranging the couplings so that the hidden sector couples differently to the inflaton than the SM, reheating can be asymmetric, whereby the SM and the hidden sector are heated to different temperatures Hodges:1993yb ; Berezhiani:1995am ; Adshead:2016xxj ; Halverson:2019kna . However, coupling both the SM and a hidden sector to the inflaton in the UV necessarily results in inflaton-mediated interactions between the two sectors. As demonstrated in reference Adshead:2016xxj , this irreducible inflaton-mediated scattering can thermalize the two sectors under fairly generic conditions.

In this work, we extend the analysis of reference Adshead:2016xxj to explore the effects of out-of-equilibrium inflaton-mediated interactions on asymmetric reheating. Along the way, we develop and implement methods to numerically solve the Boltzmann equations describing the reheating of two otherwise-decoupled sectors from the perturbative decay of the inflaton. In particular, we develop accurate approximations (including the effects of quantum statistics) for the collision terms that describe the inflaton-mediated scattering between thermalized gases of fermionic and bosonic particles. We take an effective field theory approach and consider combinations of trilinear scalar, Yukawa, and pseudo-scalar couplings between fermions, bosons, and the inflaton. When inflaton couplings to matter become sufficiently large, both non-perturbative effects such as preheating and collective effects in the radiation baths such as Landau damping and thermal masses can provide important corrections to the inflaton decay rate and hence the evolution of the temperature asymmetry, particularly at very high radiation temperatures Kolb:2003ke ; Drewes:2013iaa ; Hardy:2017wkr . However, as we show here, both inflaton decays and inflaton-mediated scattering furnish cosmological attractor solutions during the perturbative phase of reheating, making the final temperature asymmetry largely sensitive to the dynamics of the system at and below the perturbative reheat temperature $T_{\rm rh}$. Thus the perturbative reheating process that we analyze in the present paper will often serve as a good guide to the final temperature asymmetry despite the presence of richer dynamics at early times.

This paper is organized as follows. In section 2, we review standard perturbative reheating and extend the well-known single-sector results to include the effects of quantum statistics on the decay width of the inflaton. In section 3, we begin our study of reheating into two sectors by considering the effects of quantum statistics on reheating into combinations of otherwise-decoupled fermionic and bosonic sectors. In section 4, we reintroduce the inflaton-mediated interactions between the two sectors (required by self-consistency) and show how inter-sector scattering dominates over any features from quantum statistics in most of the parameter space. We conclude in section 5. Appendix A elaborates on the construction of cosmological attractor solutions, while appendix B demonstrates that the novel behavior found for bosonic radiation baths in section 3 can indeed be realized during perturbative reheating. In appendices C and D, we lay out the procedure for analytically evaluating energy transfer rates between two relativistic particles at different temperatures mediated by a massive scalar field.

We work in units where $\hbar=c=k_{B}=1$, and denote by $M_{\rm Pl}=2.435\times 10^{18}$ GeV the reduced Planck mass.

2 Quantum statistics in single-sector reheating

In this section we revisit the perturbative reheating of a radiation bath. After reviewing the classic treatment, we demonstrate that at temperatures $T>M_{\phi}/4$, where $M_{\phi}$ is the inflaton mass, quantum effects such as Bose enhancement and Pauli blocking can significantly affect the evolution of the temperature of the radiation bath during reheating. We show that the effects of quantum statistics disappear once $T$ drops below $M_{\phi}/4$, and thus alter the outcome of reheating only when $T_{\rm rh}\gtrsim M_{\phi}/4$. While we refer to the decaying particle as an inflaton and have post-inflationary reheating primarily in mind, our results apply also to other “reheatons” such as curvatons or moduli (see also Reece:2015lch ).

2.1 Perturbative reheating

A generic scenario of inflation Guth:1980zm ; Linde:1981mu ; Albrecht:1982wi consists of one or more scalar fields $\phi_{i}$ slowly rolling on a sufficiently flat potential, $V(\phi_{i})$ (see, for example, Martin:2013tda and references within). Inflation ends when the slow-roll conditions are violated, and the fields $\phi_{i}$ roll quickly to the potential minima and start oscillating. For this work, we assume that only one field $\phi$ is relevant during the reheating process, and that its potential is analytic and can be expanded in a Taylor series about its minimum. We further assume that only the leading quadratic term in this Taylor series is needed.111We are explicitly ignoring anharmonic corrections to the inflaton potential that may be relevant during reheating. These anharmonic terms can be important for non-perturbative effects during reheating, such as the formation of oscillons, as recently reviewed in Amin:2014eta . Conversely, the absence of a quadratic minima generically leads to a radiation equation-of-state very quickly following inflation Lozanov:2016hid . The time-averaged equation of state of a field oscillating in a quadratic potential is that of a stationary massive particle, and thus the Universe undergoes a period of matter domination while the inflaton energy density dominates Turner:1983he . During this oscillating phase, the inflaton condensate starts to decay through its couplings to matter, initiating reheating. If these couplings are large enough, the first stage of reheating can proceed through a period of parametric resonance known as preheating Traschen:1990sw ; Kofman:1994rk . In the preheating regime, particle production is non-perturbative and typically requires numerical treatment (however, see Emond:2018ybc ). As the amplitude of inflaton oscillations decreases, due to both Hubble friction and inflaton decay, preheating ceases and particle production can be treated perturbatively. Unless preheating is violent enough to drain an $\mathcal{O}(1)$ fraction of energy out of the inflaton condensate, this final epoch of perturbative reheating typically dominates the properties of the radiation bath produced by inflaton decays.

For this work, we thus consider perturbative reheating in a quadratic potential Abbott:1982hn ; Albrecht:1982mp . We consider the generic case where all particle masses besides the inflaton mass are negligible at the energies we consider, and therefore treat all matter species as radiation. We further neglect inverse decays from radiation into inflaton quanta; this is a good approximation provided the number of species in the radiation bath is large, $g_{*}\gg 1$. With these approximations, the Boltzmann equations describing reheating read (see, for example, Chung:1998rq )

[TABLE]

where the Hubble rate ${H}$ is given by the Friedmann equation

[TABLE]

The inflaton width is denoted by $\Gamma$, and $\rho_{\phi}$ and $\rho$ are the inflaton and radiation energy densities, respectively. These equations are (approximately) valid from the end of inflation at some scale factor $a=a_{I}$, which we take as our initial point. The radiation sector is initially empty,222This is generally a good approximation for models with tri-linear scalar couplings and Yukawa interactions with fermions, as in these cases the daughter fields get a large mass during inflation, shutting off inflaton decays. However, for a pseudo-scalar inflaton coupling to either fermions or gauge bosons, there can be significant energy density already in the radiation sector as inflation ends (see, for example, Adshead:2015pva ; Adshead:2015kza ). However, as we demonstrate below, the specific initial conditions are largely irrelevant for the detailed outcome of reheating. $\rho_{,I}=0$, whereas the initial energy density of the inflaton is given in terms of the mean-square value of the inflaton field just after the end of inflation, $\langle\phi^{2}_{I}\rangle$, as $\rho_{\phi,I}=M^{2}_{\phi}\langle\phi^{2}_{I}\rangle/2$.

In figure 1, we plot the inflaton and radiation energy densities obtained by numerically solving eqns. (1) and (2) with a constant inflaton decay width, $\Gamma=\Gamma_{0}$. Initially, $\Gamma\ll H$ and therefore inflaton decays are inefficient. Thus the inflaton energy density during this phase can be well approximated as diluting only through redshifting, $\rho_{\phi}\approx\rho_{\phi,I}(a/a_{I})^{-3}$. The evolution of the radiation sector, however, is dominated by the energy injection from inflaton decays. Initially, the radiation energy density grows rapidly until the rate at which energy is injected into the radiation bath by inflaton decays, governed by ${\Gamma}{\rho}_{\phi}$, matches the rate at which the radiation bath loses energy due to the expansion of the universe, governed by $4H\rho$. After this point, the evolution of the radiation sector follows an attractor solution, which realizes a quasi-static equilibrium between energy injection and dilution (see appendix A),

[TABLE]

where

[TABLE]

We call the above evolution imposed on the radiation sector by inflaton decays the reheating attractor curve. For a temperature-independent decay width, $\Gamma=\Gamma_{0}$, the factors $q$ and $p$ are readily determined to be constants, $q=-3/2$ and $p=0$, yielding the relation $4H\rho=(8/5)\Gamma_{0}\rho_{\phi}$. On this attractor solution the radiation bath evolves as Chung:1998rq

[TABLE]

when $a\gg a_{I}$. The attractor nature of this solution means that the evolution of the energy density of the radiation bath during reheating is relatively insensitive to its initial conditions. Radiation baths with energy density initially below the attractor solution rise rapidly to meet the attractor. Meanwhile radiation densities initially above the attractor curve redshift as $\rho\propto a^{-4}$ until they meet the attractor, as can occur in (e.g.) scenarios where a modulus comes to dominate the post-reheating universe Reece:2015lch ; Kane:2015jia ; Co:2015pka .

The attractor solution can also be obtained by solving the Boltzmann equations during reheating. Well before reheating, the inflaton condensate dominates the energy budget of the Universe; its comoving energy density is approximately constant and the Boltzmann equation describing the radiation bath can be simplified to yield Chung:1998rq

[TABLE]

Solving eq. (7) with the initial condition $\rho(a_{I})=0$ also allows us to determine the maximum energy density attained by the radiation sector Chung:1998rq ,

[TABLE]

2.2 Quantum statistics during single-sector reheating

The preceding discussion neglected the possible effects of quantum statistics during reheating. Typically, the inflaton decays at rest, producing pairs of particles at a fixed energy $M_{\phi}/2$. To quantify the possible effects of Pauli blocking or Bose enhancement of the inflaton decay, we need to specify the phase space distribution in the radiation sector. For simplicity, we assume the radiation is in thermal equilibrium,

[TABLE]

which amounts to assuming that the thermalization time scale for the radiation sector is much faster than any other time scale in the problem. This is in some sense a conservative assumption for the purpose of analyzing the scattering and reheating processes discussed in this paper: a less equilibrated sector has a greater fraction of particles with energies concentrated near $M_{\phi}/2$, making both inflaton-mediated scattering and quantum statistics more important. However, as we demonstrate below, the post-reheating properties of the radiation baths are typically determined by the late-time behavior of the system, making the detailed approach to thermal equilibrium within each radiation bath largely immaterial for the final outcome of reheating. This separation of timescales generally makes prompt thermalization a robust assumption.

For temperatures $T\lesssim M_{\phi}/4$, a constant (zero-temperature) inflaton decay width is a good approximation. At these temperatures, the phase space where particles are injected by inflaton decays, $E\sim M_{\phi}/2$ , is sparsely populated due to the fast thermalization of the injected particles, and thus the effects of Pauli blocking or Bose enhancement are negligible. However, at higher temperatures, $T\gtrsim M_{\phi}/4$, the equilibrium thermal distributions have significant support at $E\sim M_{\phi}/2$, and the inflaton decay rate can be significantly altered. The partial decay width of a parent scalar to pairs of particles in equilibrium at finite temperature is given by

[TABLE]

where $\Gamma_{0}$ is the zero temperature decay width, and the upper (lower) sign holds for bosons (fermions) in the final state. At high temperatures, the decay width is enhanced (suppressed) for bosons (fermions) due to Bose enhancement (Pauli blocking). We now consider reheating to boson and fermionic radiation separately.

Bosonic reheating:

In the case of decays to bosons, for $T\gg M_{\phi}/4$ the inflaton decay width is approximately given by $\Gamma\approx 4T\Gamma_{0}/M_{\phi}$. Using this decay width in eq. (4) immediately yields a new quasi-static equilibrium solution for the radiation bath ($q=-3/2$ and $p=1/4$), with power law evolution

[TABLE]

This reheating attractor curve can again be found analytically by solving the approximate Boltzmann equation describing the radiation bath during reheating, analogous to the Maxwell-Boltzmann result. With initial condition $\rho(a_{I})=0$, the full temperature evolution is

[TABLE]

The radiation bath attains its maximum temperature,

[TABLE]

at $a=1.6a_{I}$. For bosons, the inflaton decay width decreases with temperature, making energy injection into the radiation sector less efficient as the temperature decreases. This results in the temperature dropping as $T\propto a^{-1/2}$, faster than the classical result $T\propto a^{-3/8}$ (eq. (6)).

Obtaining radiation temperatures high enough to realize this power law requires a relatively large inflaton width. One might be concerned that such large couplings to matter place the inflaton in the regime where preheating dominates over perturbative reheating. However, we demonstrate in appendix B that, for an inflaton with trilinear couplings to scalars, there is indeed a region of parameter space where the novel Bose power law $T\propto a^{-1/2}$ is realized during perturbative reheating.

Fermionic reheating:

For an inflaton decaying to fermions at $T\gg M_{\phi}$, the decay width can be well approximated by $\Gamma\approx\Gamma_{0}M_{\phi}/(4T)$, which gives $q=-3/2$ and $p=-1/4$. In this regime, the radiation sector evolves as

[TABLE]

The full solution to the Boltzmann equations with initial condition $\rho(a_{I})=0$ is

[TABLE]

The maximum temperature attained by the radiation bath is

[TABLE]

attained at $a=1.4a_{I}$.

In the analytic treatment of the Boltzmann equations for reheating in the fermionic and bosonic cases above, we have taken $T(a_{I})=0$ as our initial condition. Strictly this is inconsistent with the high temperature expansion used for the inflaton width. A more complete analytic treatment would use the zero-temperature inflaton width to describe the early evolution of the radiation bath until its temperature rises to $M_{\phi}$ before implementing the high temperature expansion. However, such a procedure only alters the scale factor at which the maximum temperature is attained and not its value. Moreover, since the maximum temperature is attained very quickly compared to other timescales in our problem, the error due to this simplifying assumption is negligible. Perhaps the more consequential assumption in this region is that we have taken the radiation bath to attain internal thermal equilibrium nearly instantaneously. In the very early periods of reheating, the thermalization rate is likely to be smaller than the very rapid rate at which the energy density of the radiation bath grows. The simple solutions presented here for the decay width and the initial evolution of the energy densities are thus probably incorrect for describing these very early regions.

Once the temperature falls below the inflaton mass scale, the temperature dependence of the inflaton decay width in eq. (10) becomes unimportant as inflaton decays now populate sparsely occupied regions of phase space. Subsequently the radiation sector evolves as $T\propto a^{-3/8}$.

Reheating completes when the inflaton decays become efficient, $\Gamma\sim H$, and the inflaton energy density decreases exponentially. During this epoch the Universe transitions from the matter-dominated era of reheating to a radiation-dominated expansion, where the temperature of the radiation sector redshifts adiabatically as $T\propto a^{-1}$. If $\Gamma\sim H$ occurs while $T\gtrsim M_{\phi}/4$, then the temperature of the radiation sector directly transitions to $T\propto a^{-1}$ without going through the classical $T\propto a^{-3/8}$ regime. In this scenario, the resulting reheat temperature depends on the quantum statistics of the inflaton decay products. Estimating the reheat temperature by setting $H=\Gamma(T)$ and taking $H$ to be dominated by the radiation bath, we find for $T_{\textrm{rh}}\gg M_{\phi}$

[TABLE]

in contrast to the classical result

[TABLE]

which holds for $T_{\textrm{rh}}\ll M_{\phi}$. We summarize these different power law behaviors of the radiation temperature in figure 1. In the left panel we show a case where the reheat temperature $T_{\textrm{rh}}$ is below the inflaton mass. In this case, quantum statistics are unimportant for determining $T_{\textrm{rh}}$, as all three scenarios converge onto the attractor solution governing classical perturbative reheating, eq. (4). In the right panel of figure 1 we show a case where $T_{\textrm{rh}}$ is above the inflaton mass. As the inflaton decay width gets significant corrections from quantum statistics at these high temperatures, we observe the different reheat temperatures of eq. (17) expected for different quantum statistics at fixed zero-temperature decay width.

In the above scenario we have assumed that all particles coupled to the inflaton have the same quantum statistics (bosons or fermions). If the inflaton couples to both bosons and fermions then the energy density of the radiation sector as a whole evolves depending on the total inflaton decay width. In this scenario, the inflaton width is dominated by the Bose-enhanced partial widths at very high temperatures, and hence the radiation sector evolves according to the bosonic power law ($T\propto a^{-1/2}$). If the zero-temperature partial-width into fermions is larger than that to bosons, $\Gamma^{\rm fermion}_{0}>\Gamma^{\rm boson}_{0}$, then (assuming $T_{\textrm{rh}}<M_{\phi}/4$) there is a temperature, $T_{*}$, for which $\Gamma^{\rm boson}(T_{*})=\Gamma^{\rm fermion}(T_{*})$, while $\Gamma^{\rm boson}(T)<\Gamma^{\rm fermion}(T)$ for $T<T_{*}$. For $M_{\phi}/4<T<T_{*}$, the radiation bath transitions to the power law $T\propto a^{-3/10}$ (characteristic of high-temperature fermionic reheating) before ultimately transitioning to the classical $T\propto a^{-3/8}$ for $T<M_{\phi}/4$.

Despite model-dependent uncertainties associated with the initial evolution of the radiation baths, the attractor nature of these perturbative reheating solutions renders the later temperature evolution, and the resulting reheat temperatures, insensitive to variations to the initial conditions and early evolution provided that the attractor solution is obtained. Reaching the attractor solution requires that 1) the energy density of the oscillating inflaton dominates the Hubble rate for some time, during which inflaton decays become a cosmologically important source for the radiation bath, and 2) that the thermalization timescale is short compared to the duration of inflaton domination. As we demonstrate in the remainder of this paper, similarly general results can be obtained for the more complicated scenarios that arise in two-sector reheating as well.

3 Two-sector reheating and quantum statistics

In section 2, we demonstrated how quantum statistics alter the temperature evolution of the radiation sector prior to reheating. In this section we explore the extent to which these quantum statistical effects can impact the final temperature asymmetry in two-sector reheating. To isolate the effects of quantum statistics in the inflaton decay width on the resulting temperature asymmetries, in this section we temporarily ignore inflaton-mediated scattering. We begin by writing down the Boltzmann equations for two-sector reheating, and then consider the case where the inflaton couples to fields with the same statistics in both sectors and the case where the inflaton couples to fields with different statistics in each sector.

3.1 Two-sector reheating

We restrict ourselves to the regime of perturbative reheating and neglect inflaton-mediated scattering. For simplicity we continue to assume that the inflaton is only coupled to one species of particle (boson or fermion) in each sector. Ignoring any inflaton quanta, the Boltzmann equations in this limit read

[TABLE]

Here ${H}$ is the Hubble rate, which is given by the Friedmann equation

[TABLE]

and $\Gamma_{a,b}$ are the (temperature-dependent) decay rates of $\phi$ to the respective sectors. The subscript ‘$a$’ denotes the sector that attains the larger temperature at the end of reheating. Generically, this corresponds to the sector with the largest zero-temperature decay width. However, as we demonstrate below, this is not the case when the inflaton couples to bosons in one sector and fermions in the other and the resulting reheat temperature is large, $T_{\rm rh}\gtrsim M_{\phi}/4$. ‘Reheat temperature’ in this context refers to the temperature of the hotter sector when the universe transitions from matter to radiation domination. We define the transition from matter domination to radiation domination at the point where energy density in the radiation becomes equal to the energy density in the inflaton, $\rho_{a}(a_{\textrm{rh}})+$ $\rho_{b}(a_{\textrm{rh}})=\rho_{\phi}(a_{\textrm{rh}})$.

Although no scattering terms appear in eq. (3.1), the radiation baths are still coupled gravitationally through the Friedmann equation, eq. (20). Prior to reheating, when the Hubble rate is dominated by the inflaton, the two sectors are effectively decoupled and evolve independently of each other. During this phase, the temperature evolution in both radiation sectors is determined by their respective reheating attractor curves, as discussed in section 2. The reheat temperature is eventually determined by the hotter sector, and following reheating both sectors evolve adiabatically.

For the numerical results in the rest of the paper we adopt a common reference set of numerical values for the inflaton mass and initial energy density as well as the number of degrees of freedom in each radiation bath,

[TABLE]

We assume for simplicity that $\alpha_{a}$ and $\alpha_{b}$ are constant over the range of temperatures we consider. While in what follows we have fixed the value of $M_{\phi}$, our results are broadly independent of its precise value. As we demonstrate below, our results for the final temperature asymmetry depend on $M_{\phi}$ only through $T_{\mathrm{rh}}$ and the ratio $T_{\mathrm{rh}}/M_{\phi}$. The specific value of $M_{\phi}$ is generally only important insofar as smaller values of $M_{\phi}$ make it easier to obtain larger $T_{\mathrm{rh}}/M_{\phi}$.

When solving the Boltzmann equations numerically, we start the computation immediately after the end of inflation with $\rho_{\phi}(a_{I})=\rho_{\phi,I}$. However, instead of assuming the radiation sectors to be initially empty, we take them to be at their maximum temperatures, $\left.T_{a,b}(a_{I})=T_{a,b}\right|_{\textrm{max}}$ (given by eqs. (13), (16), and (8)), to improve computational speed. Since the maximum temperatures are obtained very quickly around $a\approx 1.5a_{I}$, the above approximation deviates only marginally from the exact result. Moreover, the attractor nature of the reheating process ensures that such small deviations in the early evolution of the radiation baths leave their late time evolution entirely unaffected.

3.2 Temperature asymmetries for sectors with the same quantum statistics

We first consider the case where the inflaton couples to fields with the same statistics in both sectors. Figure 2 shows the evolution of the temperatures of the two sectors as a function of the scale factor when the inflaton is coupled to either bosons (left panels) or fermions (right panels) in both sectors. The temperature evolution changes as the temperatures of the sectors drop below $M_{\phi}/4$ (indicated in the figure by asterisks), and quantum statistics become unimportant. As the two sectors reach this critical temperature $M_{\phi}/4$ at different times, the temperature ratio

[TABLE]

changes during the period where quantum statistics are important in one sector but not the other. This is displayed in the right panels of figure 2. Note that $x$ increases in the case of bosons (lower left panel) while it decreases for fermions (lower right panel). This is because Bose enhancement causes the hotter sector to redshift more quickly, $T\propto a^{-1/2}$, compared to the classical behavior $T\propto a^{-3/8}$ at $T\lesssim M_{\phi}/4$. On the other hand, fermions redshift more slowly due to Pauli blocking, $T\propto a^{-3/10}$. After reheating, both sectors evolve adiabatically with $T\propto 1/a$, and the temperature asymmetry between the two sectors is frozen in.

We define the final post-reheating temperature ratio as

[TABLE]

where $a^{\prime}>a_{\textrm{rh}}$ is some value of the scale factor at which both radiation sectors are evolving adiabatically, i.e. $\rho\propto a^{-4}$. In figure 3, we display the contours of $x_{\textrm{rh}}$ as a function of the ratio of the zero-temperature decay widths into each sector,

[TABLE]

for inflaton decays to bosons (left) or fermions (right) in both sectors. When $T_{\textrm{rh}}\ll M_{\phi}$ and $T_{\textrm{rh}}\gg M_{\phi}$, the final temperature ratio depends only on the ratio of the zero-temperature decay widths. In the intermediate region, $T_{\textrm{rh}}\sim M_{\phi}$, the contours tilt to the right (left) for the boson (fermion) case, since here reheating occurs when $T_{a}\gtrsim M_{\phi}/4$ but $T_{b}\lesssim M_{\phi}/4$.

In the next section, we show that these results, obtained in the absence of scattering, are a good guide to the temperature asymmetries obtained in the full theory in the regime when $T_{\textrm{rh}}\lesssim M_{\phi}$ (see figure 8 and 9, below). Above $T_{\textrm{rh}}\gtrsim M_{\phi}$, however, inflaton-mediated scattering between the two sectors is significant in determining the final temperature asymmetries.

3.3 Temperature asymmetries for sectors with differing quantum statistics

Finally, we consider the case where the inflaton is coupled to fermions in one sector and bosons in the other. For this scenario, we label the sector with bosons ‘$a$’, and the sector with fermions ‘$b$’.

As in the cases discussed above, prior to reheating the two radiation baths evolve independently. When $T\gg M_{\phi}$ the different quantum statistics results in different temperature evolution in the two sectors. As the temperature drops below the inflaton mass scale, $T\lesssim M_{\phi}/4$, the inflaton decay width becomes, to an excellent approximation, temperature-independent, and the temperature of both sectors evolves as $T\propto a^{-3/8}$. For reheating at temperatures $T\lesssim M_{\phi}/4$, the attractor nature of the temperature evolution in both sectors ensures that, regardless of the earlier history of both sectors, the final temperature asymmetry is set by the ratio of the zero-temperature decay widths. However, when $T_{\rm rh}\gtrsim M_{\phi}/4$, the differing temperature evolution in the two sectors is frozen into the resulting temperature asymmetries, as we now discuss.

In the case where the zero-temperature decay width is the same to both sectors, ${\Gamma}_{0b}={\Gamma}_{0a}$, if reheating occurs after quantum statistics become unimportant, the two sectors subsequently attain the same temperature. This is evident in the left panel of figure 1. However, for reheating at $T_{\textrm{rh}}>M_{\phi}/4$, the details of the quantum statistics become important. This is demonstrated in the left panel of figure 4, which shows how early reheating can freeze in a temperature asymmetry due to quantum statistics in an otherwise symmetric reheating scenario.

If the zero-temperature partial width to fermions is larger than that to bosons, ${\Gamma}_{0b}\geq{\Gamma}_{0a}$, the bosonic sector could still initially be hotter (due to the Bose-enhanced larger initial $T_{\textrm{max}}$) than the fermionic sector, and then evolve to become cooler. An example of this behavior is illustrated in the right panel of figure 4. In this case, the identity of the ‘hotter’ sector depends on the temperature at which reheating happens. In particular, for $T_{\textrm{rh}}\gtrsim M_{\phi}/4$, the fermionic sector can have a smaller temperature than the bosonic sector even when the inflaton has a larger zero-temperature partial width into fermions.

In figure 5, we plot the post-reheating temperature ratio, $x_{\textrm{rh}}=$ $\textrm{min}(T_{b}/T_{a},T_{a}/T_{b})$, as a function of ${\Gamma}_{a}$ and $w={\Gamma}_{b}/{\Gamma}_{a}$. When reheating happens below the scale of the inflaton mass, $T_{\textrm{rh}}\ll M_{\phi}$, the final temperature asymmetry is insensitive to the effects of quantum statistics and the final contours depend only on the ratio of the widths, $w$. For $T_{\textrm{rh}}\gtrsim M_{\phi}/4$, on the other hand, $x_{\rm rh}$ depends on both ${\Gamma}_{a}$ and $w$ due to the asymmetric effects from quantum statistics. In the blue region on the left, bosons are the hotter sector, while in the blue region on the right, fermions are the hotter sector. The transition between these situations occurs in the yellow region. As in the cases above with identical statistics, these no-scattering results provide an accurate guide to the final temperature asymmetries for $T_{\textrm{rh}}\lesssim M_{\phi}$ (see figure 12, below).

Unless there are additional non-adiabatic processes at late times, the final temperature asymmetry between the sectors is determined by their temperatures at reheating. Since the radiation temperatures are governed by the reheating attractor curve $4H\rho\propto\Gamma\rho_{\phi}$, the temperature asymmetry depends only on the ratio of decay widths at reheating. Hence, if reheating happens at temperatures $T_{\textrm{rh}}\lesssim M_{\phi}/4$ the final temperature asymmetry is determined by the ratio of the zero-temperature decay widths, $x=(\alpha_{a}{\Gamma}_{0,b}/\alpha_{b}{\Gamma}_{0,a})^{1/4}$. However, larger reheat temperatures probe the quantum statistics of the inflaton decay products, which can leave imprints in enhanced or suppressed temperature asymmetries.

4 Two-sector reheating with inflaton-mediated interactions

We now incorporate inflaton-mediated scattering between the two sectors and study its effect on the temperature evolution in both sectors and the final temperature asymmetry. As we demonstrate in this section, inflaton-mediated energy transfer between sectors also yields an attractor solution for the temperature of the colder radiation bath, which allows us to make analytic predictions for the final temperature asymmetries in the regime where inflaton-mediated scattering is important.

We begin by establishing our notation. Introducing the scattering terms in the Boltzmann equations, we have Adshead:2016xxj

[TABLE]

where $\mathcal{C}_{E}$ is the collision term describing the energy transfer from the hotter radiation sector to the colder radiation sector via two-to-two scattering processes of the form $1+2\rightarrow 3+4$,

[TABLE]

Here $\mathcal{C}_{E}^{f}$ and $\mathcal{C}_{E}^{b}$ are the collision terms for forward and backward reactions respectively, $|{\overline{\mathcal{M}}}|^{2}$ is the spin-summed scattering amplitude as determined by the particular inflaton-radiation interaction, and $S$ is a symmetry factor accounting for possible identical particles in the initial and/or final state. We retain the full dependence on quantum statistics to accurately describe the energy transfer between two relativistic radiation baths Adshead:2016xxj , which makes the evaluation of the collision term more challenging. In appendix C, we analytically evaluate eq. (4) for $s$-channel processes between two relativistic species at different temperatures for the specific interactions we analyze below. In appendix D, we describe a novel procedure to reduce the integral in eq. (4) for a $t$-channel scattering amplitude and verify explicitly that the energy transfer rate through $t$-channel processes is always orders of magnitude smaller than the energy transfer via $s$-channel processes.

We next demonstrate that the inflaton-mediated energy transfer yields a cosmological attractor solution for the colder radiation bath, using the model where the inflaton has trilinear couplings to scalar fields in both sectors as an illustrative example. We then analyze in detail how the interplay between this scattering attractor solution and the reheating attractor curve of the previous section determines the final temperature asymmetry. We then extend this analysis to other forms of the inflaton couplings to matter. In particular, we consider theories where the inflaton has: Yukawa couplings to fermions in both sectors; axion-like couplings to gauge bosons in both sectors; and a mixed scenario with a trilinear coupling to scalars in one sector and Yukawa coupling to fermions in the other.

The numerical analysis in the following sections uses the same parameter values specified in eq. (21). For the collision term, $\mathcal{C}_{E}$, we use the analytic approximations derived in appendix C. For simplicity we continue to assume that the inflaton only couples to a single species in each sector.

4.1 Scalar trilinear couplings

We begin by considering a theory where the inflaton is coupled to scalar fields in both sectors, $\chi_{a}$ and $\chi_{b}$, via trilinear couplings

[TABLE]

This interaction results in zero-temperature decay widths given by

[TABLE]

where $m_{a,b}$ denotes the mass of the fields $\chi_{a,b}$, which we have assumed to be much smaller than the inflaton mass, $m_{a,b}\ll M_{\phi}$. Our convention is that sector $a$ is the hotter sector, and accordingly we take $\mu_{a}\geq\mu_{b}$ in what follows.

4.1.1 The collision term

The $s$-channel amplitude for $\chi_{a}\chi_{a}\leftrightarrow\chi_{b}\chi_{b}$ scattering mediated by inflaton exchange is given by

[TABLE]

In appendix C.1, we compute the collision term, $\mathcal{C}_{E}$, following from this amplitude.

The collision term in general is a function of both $T_{a}$ and $T_{b}$. However, for large asymmetries $T_{b}\ll T_{a}$, the forward energy transfer term governing energy injection into the colder sector dwarfs the backward energy transfer term. Moreover, in this regime we can also ignore the final state Bose enhancement of $C_{E}^{f}$: while $\chi_{b}$ particles produced in the forward reaction typically have energies of order $\sim T_{a}$, for $T_{b}\ll T_{a}$ those energy levels are mostly unpopulated. In this simplified regime, the collision term thus depends only on $T_{a}$ as

[TABLE]

as derived in appendix C.1; see figure 14. For $T_{a}\gtrsim M_{\phi}$, the collision term is substantially enhanced by the resonant exchange of inflaton particles. The divergence in the Bose-Einstein distributions at $E\to 0$ combined with the resonant peak in the scattering amplitude results in a logarithmic dependence on $T_{a}/M_{\phi}$ for $T_{a}\gg M_{\phi}$. As $T_{a}$ drops below the inflaton mass, the scattering goes off resonance and $\mathcal{C}_{E}$ drops rapidly. Because the scattering is dominated by the energetic tail of the phase space distribution, this fall-off of the energy transfer rate can be accurately described assuming Maxwell-Boltzmann statistics, thus yielding a Bessel function $K_{2}(M_{\phi}/T_{a})\sim(M_{\phi}/T_{a})^{3/2}\exp(-M_{\phi}/T_{a})$. Note that, in the resonant regime, the energy transfer rate depends more strongly on the smaller coupling $\mu_{b}$ than the larger coupling $\mu_{a}$, and in particular, when $\mu_{b}\ll\mu_{a}$, the rate is almost independent of $\mu_{a}$. Below the resonance, the energy transfer rate drops rapidly until it reaches the low-temperature regime $T_{a}\ll M_{\phi}$. In this regime, the inflaton can be integrated out of the theory, leaving a constant scattering amplitude, $|\mathcal{M}(s)|^{2}\approx\mu_{a}^{2}\mu_{b}^{2}/M_{\phi}^{4}$. Thus we obtain the $\mathcal{C}_{E}\propto T_{a}^{5}$ behavior in the last line of eq. (32). Finally, at temperatures low enough that one or both of the scattering species becomes non-relativistic, the energy transfer rate becomes Boltzmann-suppressed; we do not include this effect, as we find that generically the behavior of $\mathcal{C}_{E}$ below $T_{a}<M_{\phi}/4$ is inconsequential to determining the final temperature asymmetry.

Finally, we stress that the expression for $\mathcal{C}_{E}$ given in eq. (32) is a limiting version that neglects its dependence on $T_{b}$. Dependence on $T_{b}$ can enter in two ways: first, via the backward energy transfer term, and second, from Bose enhancement of $\mathcal{C}_{E}^{f}$. The backward energy transfer term becomes important when $T_{b}\gtrsim 0.9T_{a}$, and as the two sectors approach equilibration the net energy transfer rate rapidly drops. The Bose enhancement of the forward energy transfer term is more involved to model. This Bose enhancement largely serves to increase $\mathcal{C}_{E}^{f}$ in the high and low temperature regimes in eq. (32) with increasing $T_{b}$. The middle regime in eq. (32), however, is insensitive to the possible Bose enhancement terms, as that regime is effectively described by Maxwell-Boltzmann statistics. As we show below, this last property enables us to obtain analytic predictions of the final temperature asymmetry without needing to keep track of the full behavior of the Bose enhancements.

4.1.2 The scattering attractor solution

Now we discuss the impact of the collision term on the temperature evolution of both sectors, and derive the corresponding scattering attractor curve for the temperature of the colder sector. We begin by considering scenarios where $T_{\textrm{rh}}\gg M_{\phi}$, where, as we demonstrate, scattering becomes important post-reheating. One such parameter point is shown in figure 6, which plots numerical solutions for the radiation temperatures obtained using the collision term as given in eq. (C.1). To highlight the importance of scattering, we also show the temperature evolution when the scattering has been turned off. In the right panel, we show the evolution of the collision term, $\mathcal{C}_{E}$, in comparison to the combinations $4H\rho_{a,b}$ and $\Gamma_{a,b}\rho_{\phi}$. We mark the point $T_{a}=M_{\phi}/4$ around which $\mathcal{C}_{E}$ begins to exhibit substantial Boltzmann suppression. As the scattering process affects the temperature evolution substantially post-reheating in this example, we can cleanly separate the effects of scattering from the contributions of reheating; in this discussion, reheating itself is only important insofar as it provides initial conditions for the subsequent post-reheating evolution of $T_{a}$ and $T_{b}$.

As figure 6 shows, $T_{b}$ begins to deviate from the no-scattering solution as soon as the fractional energy transfer rate into the colder sector becomes comparable to the Hubble rate, $\Gamma_{E,b}=\mathcal{C}_{E}/\rho_{b}\sim H$. When this happens we say that inflaton-mediated scattering becomes effective. In contrast, when the fractional energy transfer rate out of the hotter sector becomes comparable to the Hubble rate, $\Gamma_{E,a}=\mathcal{C}_{E}/\rho_{a}\sim H$, inflaton-mediated scattering becomes efficient and the two sectors attain thermal equilibrium. In the scenario shown in figure 6, inflaton-mediated scattering becomes effective but never efficient. The solution to the Boltzmann equation for $T_{b}$ when scattering becomes effective is approximated by the quasi-static attractor solution (see appendix A)

[TABLE]

where

[TABLE]

We call this evolution of $\rho_{b}$ the scattering attractor curve. In evaluating $q(a)$, the scale factor dependence in $\mathcal{C}_{E}/H$ comes through $T_{a}$, which in the present scenario is evolving adiabatically. For a given value of $T_{a}$, there is a single corresponding value of $T_{b}$ that satisfies eq. (33). In general, solving eq. (33) for $T_{b}$ is non-trivial given the dependence of $\mathcal{C}_{E}$ on $T_{b}$ through Bose enhancement. The attractor curve exists as long as $4+q/(1-p)>0$, which translates to the condition that $\mathcal{C}_{E}$ falls off more slowly with scale factor than $Ha^{-4}$. At temperatures below $T\sim M_{\phi}/4$ the collision term falls off exponentially (eq. (32)), marking the end of the attractor evolution. Beyond that point, $\rho_{b}$ evolves adiabatically as seen in figure 6. Thus the scattering attractor curve yields a final temperature asymmetry simply given by the asymmetry at $T_{a}\approx M_{\phi}/4$.

To further highlight the attractor nature of the collision term, figure 7 shows the post-reheating evolution of the temperature ratio, $x=(\alpha_{a}\rho_{b}/(\alpha_{b}\rho_{a}))^{1/4}$, for the parameter point of figure 6, but now considering a range of (post-reheating) initial conditions for $\rho_{b}$ (or equivalently $x$). In the left panel, the solid blue line tracks the evolution of $T_{b}/T_{a}$ following from figure 6, where the initial conditions are determined self-consistently from inflaton decays, $x_{i}=x_{\textrm{rh}}=0.02$. The purple dot-dashed line shows the evolution when the initial temperature ratio is instead zero, $x_{1,i}=0$; again, initial densities below the attractor solution rise rapidly to attain the attractor. The yellow dot-dashed line shows the evolution with an initial temperature ratio $x_{2,i}=0.1>x_{\textrm{rh}}$; this solution still attains the scattering attractor curve, eq. (33). The green dot-dashed line denotes evolution with an initial temperature ratio, $x_{3,i}=0.5$, much above the final temperature ratio determined by the scattering attractor solution. In this case $T_{b}$ remains mostly unaffected by inflaton-mediated interactions. Thus we see that the inflaton-mediated interactions impose a minimum value for the final temperature ratio: any initial temperature ratio below this minimum value is increased to that value by the scattering attractor solution while initial temperature ratios above this minimum remain largely unaffected. This minimum final temperature ratio is simply determined by the behavior of $\mathcal{C}_{E}$ near $T_{a}\sim M_{\phi}/4$, as we elaborate below.

The right panel of figure 7 shows the ratio $\mathcal{C}_{E}(T_{a},x)/(Ha^{-4})$ for all scenarios. The different behavior of the curves at early times shows the differing importance of Bose enhancement on the final state in the different scenarios. All curves converge onto a single common solution for $T_{a}\sim M_{\phi}$, when scattering is well described by Maxwell-Boltzmann statistics. As $H\rho_{a}$ evolves with scale factor as $Ha^{-4}$, the value of $a$ where $\mathcal{C}_{E}/Ha^{-4}$ is maximized333The visible hiccup near the peak in $\mathcal{C}_{E}(T_{a},0)$ is a feature of the imperfect fitting used in our analytic approximation of $\mathcal{C}_{E}$, eq. (C.1). coincides with the value of $a$ where the scattering attractor solution for $\rho_{b}$ ends.

To accurately determine the final temperature asymmetry predicted by the scattering attractor-curve, we need to integrate the Boltzmann equations around the region where $T_{a}\lesssim M_{\phi}$. Assuming no $T_{b}$ dependence in $\mathcal{C}_{E}$ or $H$, the Boltzmann equation for $\rho_{b}$ can be directly integrated,

[TABLE]

where we have defined $a_{1}$ as the scale factor at which $T_{a}(a_{1})=M_{\phi}$, and $z=a^{\prime}/a_{1}$. We evaluate $\mathcal{C}_{E}$ by taking the limit $x\rightarrow 0$ around $T_{a}\sim M_{\phi}$ in eq. (C.1) to obtain

[TABLE]

We determine the initial energy density of the colder sector, $({\rho}_{b}a^{4})_{a=a_{1}}$, by assuming that the colder sector is already on the scattering attractor curve defined by evaluating eq. (33) using $\mathcal{C}_{E}=\mathcal{C}_{E}(T_{a},0)$. However, the Maxwell-Boltzmann behavior of the collision term in this temperature range helps to ensure that the final result is insensitive to the specific choice of $x=0$. We find

[TABLE]

The final energy density of the colder sector is then given by

[TABLE]

The final temperature ratio between the two sectors predicted by inflaton-mediated scattering is then

[TABLE]

where we have used $\rho_{a}(a)=\alpha_{a}M_{\phi}^{4}(a_{1}/a)^{4}$. This is the value that the temperature ratios $x,x_{1}$ and $x_{2}$ asymptote to as shown by the horizontal black dashed line in figure 7. Eq. (39) holds as long as $x_{sc}\lesssim 0.9$. Once the temperature ratio approaches unity, backward energy transfer and the contribution of $\rho_{b}$ to the Hubble parameter become important, and the attractor solution no longer captures the full behavior of the system. In these cases, where the two sectors approach thermalization, a more detailed numerical study is required.

Finally, it is worth emphasizing that the scattering attractor curve discussed here is dominated by the resonant behavior of the energy transfer rate, and depends on the properties of the radiation baths at $T_{a}\sim M_{\phi}$. In the trilinear scalar model, a second attractor phase appears at temperatures well below the resonance ($T_{a}\ll M_{\phi}$). This is evident from the late-time increase in $\mathcal{C}_{E}/(Ha^{-4})$ in the right panel of figure 7, after the resonant enhancement ends. This possibility of IR thermalization is a special feature of the trilinear scalar model, where integrating out the inflaton introduces a renormalizeable quartic interaction between the two sectors. In all other cases $\mathcal{C}_{E}$ falls off much faster at lower temperatures due to the higher ($\geq 4$) mass dimensions of the operators that couple the inflaton to the radiation baths. Once $T_{a,b}\sim m_{a,b}$, $\mathcal{C}_{E}$ becomes exponentially suppressed and scattering is cut off. Thus, thermalization in the IR depends on the mass scales in the matter sectors coupled to the inflaton, as well as the inflaton mass and $T_{\rm rh}$. Late-time equilibration through scalar portal interactions is studied in detail in Krnjaic:2015mbs ; Adshead:2016xxj ; Evans:2017kti and we do not discuss it further here.

4.1.3 Final temperature ratios

Both reheating and scattering, considered independently, produce attractor solutions. In most of parameter space, one attractor solution dominates over the other and thus is primarily responsible for determining the final temperature asymmetry. As demonstrated above in section 4.1.2, when $T_{\rm rh}>M_{\phi}$, a good semi-analytic approximation to the final temperature asymmetry is therefore

[TABLE]

Both $x_{\textrm{rh}}$ and $x_{sc}$ can be straightforwardly computed from the Lagrangian parameters without any need to solve the full Boltzmann equations.

In the case $T_{\textrm{rh}}\ll M_{\phi}/4$, the reheating attractor solution dominates, as we now show. As $\mathcal{C}_{E}$ redshifts more slowly than $\Gamma_{b}\rho_{\phi}$ during reheating, it suffices to show that $\mathcal{C}_{E}/\Gamma_{b}\rho_{\phi}<1$ at $T_{a}\approx M_{\phi}/4$ when $\mathcal{C}_{E}$ is maximized. Using eq. (4) for $\rho_{\phi}$ at $T_{a}=M_{\phi}/4$ and eq. (100) for $\mathcal{C}_{E}(T_{a}=M_{\phi}/4)$ (at $x=0$) we find

[TABLE]

This ratio is small by assumption for SM-scale values of $\alpha_{a}$. Hence for low reheat temperatures inflaton-mediated scattering is unimportant during reheating. After reheating, scattering cannot thermalize the two sectors as the resonant enhancement has already ended. Thus, for $T_{\textrm{rh}}\ll M_{\phi}$ the final temperature asymmetry is simply given by $x_{\textrm{rh}}=(\alpha_{a}\Gamma_{0b}/\alpha_{b}\Gamma_{0a})^{1/4}$. Although $x_{sc}$, as defined in eq. (39) as the temperature ratio obtained by the post-reheating scattering attractor curve, does not strictly pertain in this case, one can check that its value is always less than $x_{\textrm{rh}}$ when $T_{\textrm{rh}}\ll M_{\phi}/4$. Thereby we can extend eq. (40) to hold for $T_{\textrm{rh}}\ll M_{\phi}/4$ as well.

We are now ready to consider the full numerical solution to the Boltzmann equations, eq. (25). The resulting numerical temperature asymmetry, $x_{f,n}$, is shown in the left panel of figure 8 as a function of the ratio of the zero-temperature widths $w=\Gamma_{0,b}/\Gamma_{0,a}$.

As expected, for low reheat temperatures the reheating attractor curve dominates the evolution. The contours below $T_{\textrm{rh}}\lesssim M_{\phi}$ show the same behavior as in the absence of inflaton-mediated scattering; compare the left panel of figure 3. Inflaton-mediated scattering becomes important roughly for $T_{\textrm{rh}}\gtrsim M_{\phi}$, and dominates for $T_{\textrm{rh}}\gtrsim 10M_{\phi}$. In this high-temperature regime the contours are almost diagonal, reflecting the fact that $x_{sc}$ is dominantly governed by the smaller decay width $\Gamma_{0,b}$. In the right panel of figure 8 we compare our analytic estimate of eq. (40)444The temperature ratio $x_{\textrm{rh}}$ resulting from the reheating attractor curve used in eq. (40) is evaluated by solving eq. (3.1) numerically. with the result obtained from numerically solving the Boltzmann equations. The analytic estimate agrees with the numerical results within 20%. The discrepancies are greatest exactly where the scattering and reheating attractor curves are no longer individually sufficient to capture the full behavior of the system: when both scattering and inflaton decays are important for determining the final asymmetry, around $T_{\textrm{rh}}\sim\,\mathrm{few}\,M_{\phi}$, and when the sectors are approaching (but not obtaining) thermalization, $x_{f}\sim 0.7-0.8$.

Finally, let us note the important point that our numerical results for $x_{f,n}$ are themselves based on analytic approximations to the collision term. Our analytic fits to the collision term deviate from the exact numerical values by almost 50% near $T\sim M_{\phi}$ (see figure 14). As the final temperature ratio is predominantly determined by the behavior of the collision term near $T\sim M_{\phi}$, this error is unfortunately not negligible for our final results. However, this error is made less numerically consequential once we take the fourth root to find the temperature (eq. (39)), inducing uncertainties of up to $\sim 15\%$ in the numerical temperature ratio plots at high $T_{\textrm{rh}}$, figure 8.

4.2 Final temperature asymmetry for other theories

The two key properties of $\mathcal{C}_{E}$—the exponential suppression at $T_{a}\lesssim M_{\phi}/4$ and the weak dependence on $T_{b}$ in this range—that allowed us to analytically determine the final temperature asymmetry for the scalar trilinear case are generic features of resonant $s$-channel interactions. Much of our analysis in the previous section can thus be applied directly to other interaction structures. As we demonstrate, in models where the inflaton has renormalizeable couplings to matter, scattering is only important for determining the final temperature asymmetry when the endpoint of the scattering attractor curve occurs post-reheating. However, scattering during reheating can also be important when the inflaton is a pseudoscalar with dimension-five couplings to gauge fields in both sectors.

4.2.1 Yukawa couplings

We begin with a model where the inflaton has Yukawa couplings to fermions in both sectors,

[TABLE]

This interaction results in zero-temperature inflaton decay widths given by

[TABLE]

where $m_{a,b}\ll M_{\phi}$ denotes the mass of fields $\psi_{a,b}$. The $s$-channel spin-summed scattering amplitude between the two species is

[TABLE]

The total energy transfer collision term, $\mathcal{C}_{E}$, following from this amplitude is discussed in appendix C.2 and shown in figure 15. Unlike the scalar case discussed in section 4.1, the collision term is almost insensitive to the temperature of the colder sector unless the two temperatures are very close and $\mathcal{C}_{E}^{b}$ becomes important. In the limit that the temperature ratio between the two sectors is very small, $x\ll 1$, $\mathcal{C}_{E}$ is approximately given by

[TABLE]

At temperatures much larger than the inflaton mass, the inflaton mass can be neglected and the scattering amplitude is approximately constant, $|\overline{\mathcal{M}}(s)|^{2}\approx y_{a}^{2}y_{b}^{2}$, yielding the $\mathcal{C}_{E}\propto T^{5}$ behavior required from dimensional analysis. At temperatures closer to the inflaton mass, the energy transfer rate is resonantly enhanced, yielding $\mathcal{C}_{E}\propto T^{3}$ behavior. As the temperature drops below the inflaton mass, the energy transfer rate is dominated by resonant scattering in the Boltzmann-suppressed tails. Analogously to the scalar case, $\mathcal{C}_{E}$ can be well modeled in this region using Maxwell-Boltzmann statistics. In the low temperature regime the scattering amplitude can be approximated as $|\overline{\mathcal{M}}(s)|^{2}\approx y_{a}^{2}y_{b}^{2}s^{2}/M_{\phi}^{4}$, yielding the steep $\mathcal{C}_{E}\propto T^{9}$ behavior. Note that, like the scalar trilinear case, the energy transfer rate depends most strongly on the smaller coupling in the resonant regime.

We can again obtain an analytic expression for the final temperature asymmetry due to inflaton-mediated scattering, as we did for scalars in section 4.1.2. Using $\mathcal{C}_{E}$ from eq. (C.2) and taking $x\rightarrow 0$, we obtain

[TABLE]

The final temperature asymmetry can then be estimated using eq. (40), i.e., by comparing the lower bounds from the scattering and reheating attractor solutions. In figure 9 we show numerical final temperature ratios in the left panel and in the right panel compare our analytic estimate to the numerical results. We again observe a transitional region around $T_{\rm rh}\sim\,\mathrm{few}\times M_{\phi}$ where both reheating and scattering are important for determining the final value of $x_{f}$. Note that the analytic estimate from the scattering attractor curve has better agreement with the numerical results in the region near thermalization, $x_{f}\to 1$, than we saw for the scalar case; this is because the Fermi blocking of $C_{E}^{f}$ that occurs here is nowhere near as large an effect as the Bose enhancement we discussed in the previous subsection.

4.2.2 Axionic couplings to gauge bosons

We next consider a theory where a pseudo-scalar inflaton couples to gauge bosons in both sectors,

[TABLE]

This interaction results in zero-temperature decay widths given by

[TABLE]

where $m_{a,b}\ll M_{\phi}$ denotes the mass of the gauge fields, $A^{\nu}_{a,b}$. The $s$-channel spin-summed amplitude for $A_{a}A_{a}\leftrightarrow A_{b}A_{b}$ scattering mediated by inflaton exchange is

[TABLE]

In appendix C, we derive the total energy transfer rate, $\mathcal{C}_{E}$ for this amplitude; see figure 16. When the temperature ratio between the two sectors is very small, $x\ll 1$, the temperature dependence of $\mathcal{C}_{E}$ is approximately given by

[TABLE]

The steep rise in the collision term ($\mathcal{C}_{E}\propto T_{a}^{9}$ ) at high temperatures is a consequence of the high mass-dimension of the operators mediating the interaction. This behavior will be modified when $T_{a}\gtrsim\Lambda_{a}$ and the effective field theory breaks down.

Repeating the calculation from section 4.1.2, using $\mathcal{C}_{E}$ from eq. (124) with $x\rightarrow 0$, we obtain the asymptotic temperature asymmetry resulting from the scattering attractor curve,

[TABLE]

The final temperature asymmetry can then be estimated using eq. (40). In the left panel of figure 10 we show the final temperature ratio determined by numerically solving the Boltzmann equations. In this section, we take $\rho_{\phi,I}=10^{-10}M_{\phi}^{2}M_{\rm Pl}^{2}$ in order to keep $T_{\rm max}<\Lambda_{a}$ in all of our parameter space, thus ensuring that the effective field theory is valid throughout the entire evolution of the system. Due to the attractor nature of the Boltzmann equations describing reheating, larger values of $\rho_{\phi,I}$ do not change the final value of $x$ that one would compute for a given set of Lagrangian parameters. However, changing $\rho_{\phi,I}$ does alter the maximum temperature attained (see eq. (13)), and therefore if we require $T_{\rm max}<\Lambda_{a},\Lambda_{b}$ then we are restricted to parameters that satisfy

[TABLE]

In the left panel of figure 10 the red dot-dashed lines indicate where $T_{\rm max}=\Lambda_{a}$ for different values of $\rho_{\phi,I}$. Above those lines $T_{\rm max}>\Lambda_{a}$, and thus the early evolution of the system probes the theory above the cutoff. In the right panel of figure 10 we compare our analytic estimate to the numerical result.

In the top left corner of the right panel of figure 10, large discrepancies between the analytic estimate and the numerical computation are becoming evident. In the same region in the left panel, the contours of fixed temperature asymmetry are beginning to extend more deeply into the region of small $w$ than the previous examples. Both these features are the consequence of early (i.e. pre-reheating) thermalization, enabled by the UV-dominated energy transfer process ($\mathcal{C}_{E,UV}\propto T_{a}^{9}$) whose effects are not incorporated into the analytic estimate in eq. (40). At sufficiently high temperatures, $T_{a}$, $\mathcal{C}_{E,UV}$ dominates over the energy dumped from the inflaton. This UV behavior can be seen in figure 11, which shows the various contributions to the evolution of the energy density in eq. (25) as a function of scale factor. Because $\mathcal{C}_{E,UV}$ redshifts faster than $\Gamma_{0a,b}\rho_{\phi}$, the energy injection due to inflaton decays can exceed $\mathcal{C}_{E,UV}$ before reheating terminates, $\mathcal{C}_{E,UV}(T_{\textrm{rh}})<\Gamma_{0,b}(T_{\textrm{rh}})\rho_{\phi}$. When this occurs, (see, e.g. the left plot in figure 11), the temperature ratio at the end of reheating is the same as the one obtained due to the reheating attractor. Thus, asymmetric reheating overwhelms the collision term. However, when $\mathcal{C}_{E,UV}(T_{\textrm{rh}})>\Gamma_{0,b}(T_{\textrm{rh}})\rho_{\phi}$, the temperature ratio at the end of reheating, $x_{\textrm{rh}}$, is larger than the case without scattering, i.e. the result obtained from the reheating attractor. This deviation would not be reflected in the final temperature asymmetry if $x_{sc}$ is larger than this modified $x_{{}_{\textrm{rh}}}$ (center plot in figure 11). It is only when the modified $x_{\rm rh}$ due to $\mathcal{C}_{E,UV}(T_{\textrm{rh}})$ is larger than $x_{sc}$ (right plot in figure 11), that the effects from $\mathcal{C}_{E,UV}$ impact the final temperature ratio as we see in the top left corner of figure 10. It is worth recalling that thermal effects beyond the scope of this paper, in particular Landau damping and thermal blocking, can be important for determining the duration and dynamics of reheating in the high-$T_{\textrm{rh}}$ regime where the effects from $\mathcal{C}_{E,UV}$ show up.

4.2.3 Mixed Yukawa and trilinear couplings

Finally, we consider a theory in which inflaton has trilinear couplings to scalars in sector $a$ and Yukawa couplings to fermions in sector $b$,

[TABLE]

This interaction results in zero-temperature partial widths given by

[TABLE]

The spin-summed $s$-channel scattering amplitude between the two sectors is

[TABLE]

Using this scattering amplitude we derive the total energy transfer rate, $\mathcal{C}_{E}$, given in eq. (C.4); see figure 17. The collision term is almost insensitive to $T_{b}$ except when $T_{b}\approx T_{a}$. However, since the two sectors have different quantum statistics, the behavior of the collision term changes depending on which sector is hotter. When there is a large temperature asymmetry between the two sectors ($T_{a}\gg T_{b}$ or $T_{b}\gg T_{a}$), analogous to the cases considered above, the temperature dependence of $\mathcal{C}_{E}$ is approximately

[TABLE]

where the minus signs appear when $T_{b}>T_{a}$, as consistent with our definition of the energy transfer term in eq. (25).

Determining the final temperature asymmetry due to inflaton-mediated scattering as in section 4.1.2, making use of $\mathcal{C}_{E}$ from eq. (C.4) with $x\rightarrow 0$, we find

[TABLE]

where $\alpha_{\textrm{hot}}$ ($\alpha_{\textrm{cold}}$) denotes the value of $\alpha=\pi^{2}g_{*}/30$ corresponding to the hotter (colder) sector. The final temperature asymmetry can then be estimated using eq. (40). In the left panel of figure 12 we show numerical results for the final temperature ratio. In the right panel of figure 12 we show the disagreement between our analytic estimate and the numerical result as a percentage of the numerical result.

5 Summary and conclusion

Asymmetric reheating is a minimal way to populate dark sectors that are otherwise completely decoupled from the SM following inflation. In this work, we have performed the first detailed analysis of perturbative asymmetric reheating. Specifically, by solving the Boltzmann equations describing the perturbative decay of the inflaton into two otherwise decoupled radiation sectors, we have studied in detail the resulting temperature asymmetries attained by the sectors. Scattering processes mediated by inflaton exchange couple the two sectors in the UV, and our self-consistent treatment takes into account the associated collision terms that transfer energy between the radiation sectors. Furthermore, we have carefully accounted for the effects of quantum statistics. At high temperatures (compared to the relevant mass scale in the problem, the inflaton mass) these quantum-statistical effects lead to important corrections in both the inflaton decay rate, as well as the inflaton-mediated scattering processes that transfer energy between the sectors.

The system of Boltzmann equations describing the evolution of the energy densities in the various sectors is a coupled set of three first-order non-linear differential equations, and a general analytic solution is not available. However, in this work we have demonstrated that the system can be accurately analyzed by making use of the attractor nature of its solutions. Broadly, we have identified two classes of quasi-static attractor solutions to which the energy density of the radiation bath evolves depending on the physical process that is dominating the evolution. In a broad range of parameter space and to a good approximation, at any given time the evolution is dominated by either 1) the energy injection from the decay of the inflaton, 2) the transfer of energy between the sectors through inflaton-mediated scattering, or 3) the adiabatic expansion of the Universe. Case 1) leads to a reheating attractor curve, case 2) yields a scattering attractor curve, while in case 3) the radiation density simply redshifts as $\rho\propto a^{-4}$. As we have demonstrated, the utility of these attractor solutions is that they allow for a very accurate semi-analytic determination of the resulting temperature asymmetry between the sectors; the asymmetry is simply determined by the process which dominates the evolution at the latest time.

Our results for the temperature asymmetries generated by asymmetric reheating are surprisingly universal across various coupling structures and particle types. The key property that determines the outcome of asymmetric reheating is the reheating temperature, $T_{\textrm{rh}}$ relative to the inflaton mass-scale, $M_{\phi}$, as follows:

• •

When $T_{\textrm{rh}}\ll M_{\phi}/4$, the temperature asymmetry is solely determined by perturbative reheating process. More specifically, when $T_{\textrm{rh}}<M_{\phi}/10$, the final temperature ratio is simply given by the ratio of the zero-temperature decay widths, $x=T_{b}/T_{a}=(\alpha_{a}\Gamma_{0b}/\alpha_{b}\Gamma_{0a})^{1/4}$. As the reheat temperature is increased (but still $<M_{\phi}$) quantum-statistical corrections to the inflaton decay width begin to significantly affect the final temperature asymmetry. In this region asymmetric reheating can be achieved by quantum statistical effects alone, with otherwise identical couplings.

• •

When $T_{\textrm{rh}}\gg M_{\phi}/4$, the final temperature asymmetry is determined solely by inflaton-mediated scattering. Inflaton-mediated energy transfer between the sectors falls off exponentially when the temperature of the hotter sector falls below $T_{a}<M_{\phi}/4$ due to the $s$-channel scattering process going off-resonance. If the radiation sectors have not thermalized by this time, the colder sector is populated by a freeze-in like process where its final density (or equivalently the temperature ratio $x=T_{b}/T_{a}$) is primarily determined by the collision term and the Hubble rate at $T_{a}=M_{\phi}/4$,

[TABLE]

Because the collision term $\mathcal{C}_{E}(T_{a}=M_{\phi}/4)$ is largely insensitive to the inflaton coupling to the hotter sector as well as (at $T_{a}=M_{\phi}/4$) the quantum statistics of the interacting particles, the final temperature ratio is determined solely by the coupling strength of the colder sector irrespective of its particle identity.

In the region $T_{\textrm{rh}}\sim M_{\phi}/4$, both reheating and scattering are important in determining the final temperature asymmetry. We find that the final temperature asymmetry, as a function of $T_{\textrm{rh}}$ and the ratio of zero temperature partial widths $w$ depends on the inflaton mass only through $T_{\textrm{rh}}/M_{\phi}$. However, lower inflaton masses allow for the consistent realization of higher values of $T/M_{\phi}$ prior to reheating, which can be particularly important for models where the inflaton couples to the radiation baths through non-renormalizeable interactions (as in the axionic coupling to gauge bosons considered here).

The primary goal of this paper was to analyze, in detail, the temperature evolution of two otherwise-decoupled radiation sectors during and after asymmetric reheating, but along the way we obtained a number of other novel results. We found novel power laws describing the evolution of radiation baths during reheating at temperatures larger than the inflaton mass scale, when quantum statistics are important. We developed methods to derive closed form (approximate) analytic expressions for energy transfer rates between two relativistic particles at different temperatures via $s$-channel interactions mediated by a massive scalar field. Finally, we derived reduced integral-expressions for energy-transfer rate between two relativistic sectors at different temperatures via $t$-channel interactions.

The analytic estimates of the final temperature ratio developed here for two-sector reheating can be straightforwardly extended to $N$-sector reheating scenarios Arkani-Hamed:2016rle . In such cosmologies, for each of the sub-dominant sectors, the dominant energy injection from scattering is the collision term determined by the hottest sector. Provided the expansion rate is dominated by a single component (either the inflaton, or a single dominant radiation bath), to a very good approximation, the subdominant sectors are insensitive to each others presence.

In this work, we limited our analysis to perturbative reheating, ignoring the effects from 1) incomplete internal thermalization in the sectors during early reheating, 2) thermal modifications to the inflaton decays from collective effects, such as thermal blocking or Landau damping, 3) back-scatterings into inflaton quanta and 4) preheating. As long as these effects do not significantly alter the final reheating temperature obtained from perturbative reheating, our results for the final temperature asymmetry remain robust. Even in scenarios where such effects do significantly alter the reheat temperatures, the scattering attractor curve provides a strict upper bound to the temperature asymmetry between the sectors, $x\geq x_{sc}$ (see eq. (40)), as long as reheating occurs before inflaton-mediated scattering drops off resonance, $T_{\textrm{rh}}>M_{\phi}/4$. We leave the further study of temperature asymmetries under these potentially disruptive effects to future work. Another possible extension of this work is to study scenarios that include large asymmetries in the number of degrees of freedom in the two sectors. In such a scenario, the sector with the higher temperature could have a sub-dominant energy density, a possibility we explicitly ignored in this work.

Acknowledgements.

We thank Yanou Cui and Scott Watson for useful discussions. The work of PA and PR was supported in part by NASA Astrophysics Theory Grant NNX17AG48G. The work of JS and PR was supported in part by DOE Early Career grant DE-SC0017840.

Appendix A Cosmological attractor solutions

In this section we provide a general discussion of an important class of solutions to the cosmic Boltzmann equations that exhibit extended periods of quasi-static evolution. Such quasi-static equilibria can occur in any scenario where a quantity is fed by an external source at a faster rate than it is diluted by cosmic expansion. These quasi-static equilibria are attractor solutions in a sense that we make precise in this appendix. Our primary interest here is in the energy density contained in a thermal radiation bath, where notable examples of such attractor solutions include the $T\propto a^{-3/8}$ evolution of a radiation bath during (classical, perturbative) reheating Chung:1998rq and the $T\propto a^{-3/4}$ behavior of a radiation bath fed by out-of-equilibrium renormalizeable scattering processes Chu:2011be ; egs . Another important class of examples is realized by various models of freeze-in dark matter McDonald:2001vt ; Hall:2009bx , where the relevant quantity is the number density of DM.

We begin with a representative Boltzmann equation for the energy density of a radiation sector, given by

[TABLE]

where $\mathcal{C}_{E}$ is the net energy input into the sector from an external source. We rewrite this equation in the form

[TABLE]

By defining new variables as

[TABLE]

we can further modify eq. (59) to yield

[TABLE]

This equation dictates how the ratio of $H\rho/\mathcal{C}_{E}$ evolves depending on the functional behavior of $F(\rho,a)=\mathcal{C}_{E}/H$ encoded in $p,q$. Now note that for $p<1$ and $q>-4(1-p)$,

[TABLE]

is a stable attractor solution for this equation, provided that $p$ and $q$ slowly vary with $a$ ($\frac{\partial\ln(p,q)}{\partial\ln a}\ll 1$). This solution is an attractor: radiation baths initially below this steady-state solution rise up very rapidly to meet it, while radiation baths initially above it redshift as $\rho\propto a^{-4}$ until the attractor solution is attained. Thus, the quasi-static behavior of $\rho$ can be found by simply solving the equation

[TABLE]

In cases of cosmological interest $F$ very frequently has power law dependence on $\rho$ and $a$, thus making $\lambda$ a fixed and readily computable constant (usually of $O(1)$). In such cases, the relevant power law describing the temperature evolution can then be quickly obtained by solving $\rho\ \propto\mathcal{C}/H$.

When $q<-4(1-p)$, there is no attractor solution (as $\lambda$ is always positive) and $\lambda$ increases uncontrollably. This corresponds to the cases when $\mathcal{C}_{E}$ is falling faster than the redshifting of the energy density, and the evolution of the radiation bath is approximately adiabatic. On the other hand, when $p>1$, the attractor solution (when it exists) is not stable. If $\mathcal{C}_{E}$ ever came to dominate in this scenario then it would lead to an indefinite explosive rise in $\rho$ due to the positive feedback from $\mathcal{C}_{E}$. The subsequent solution can be obtained by simply solving $\dot{\rho}=\mathcal{C}_{E}$.

We can perform an analogous exercise for number density. The Boltzmann equation we start with here is

[TABLE]

and, defining

[TABLE]

we may rewrite this equation as

[TABLE]

Then the attractor solution is given by

[TABLE]

or

[TABLE]

Example: perturbative reheating of bosons at high temperatures

As an example, we apply the above formalism to reheating a thermal bath of bosons at temperatures $T\gg M_{\phi}$. In this regime,

[TABLE]

The corresponding quantity $F(\rho,a)$ is thus given by

[TABLE]

where again $\alpha=\pi^{2}g_{*}/30$, giving $p=1/4$ and $q=-3/2$ (eq. (60)). Eq. (63) then gives

[TABLE]

Solving for the temperature evolution then yields {T}(a)=\Big{(}2\sqrt{3}{M}_{\textrm{Pl}}{\Gamma}_{0}\sqrt{\rho_{\phi,I}}/(\alpha M_{\phi})\Big{)}^{1/3}(a/a_{I})^{-1/2}, exactly the expected asymptotic behavior.

Appendix B Preheating and the Bose power law

In this appendix we demonstrate that it is possible to realize a radiation bath following the Bose power law $T\propto a^{-1/2}$ of eq. (11) during perturbative reheating. We focus on a theory with an inflaton, $\phi$, coupled to a scalar field, $\chi$, via the trilinear coupling

[TABLE]

This model can experience broad resonance preheating for sufficiently large values of $\mu$ and sufficiently large inflaton oscillations $\Phi=\sqrt{\langle\phi^{2}\rangle}$, in which case energy is effectively drained from the inflaton condensate before perturbative reheating can occur. The condition that no broad resonances are present in the theory is

[TABLE]

When this condition is satisfied, preheating is inefficient and perturbative reheating dominates the properties of the radiation bath Traschen:1990sw ; Kofman:1994rk ; Dufaux:2006ee .

During perturbative reheating, for a given value of $\mu$, the inflaton amplitude $\Phi$ uniquely specifies the temperature of the matter sector. Using the reheating attractor solution eq. (4) along with the quantum statistics correction to the inflaton decay width eq. (10), we obtain for the radiation bath

[TABLE]

where

[TABLE]

Eq. (74) can be solved to yield $T$ as a function of $\Phi$. At high and low temperatures, the above relation simplifies to

[TABLE]

Eq. (74) is valid as long as $\Gamma(T)<H$, or

[TABLE]

In figure 13 we show the resulting parameter space for perturbative reheating for three different values of $M_{\phi}$. We show the equalities corresponding to broad resonance preheating (eq. (73)) and the end of perturbative reheating (eq. (77)) in red and blue respectively; the yellow shaded region represents the region where perturbative reheating dominates the evolution of the radiation bath. We further show contours of $T$, from eq. (74). Above $T=M_{\phi}$, the matter sector realizes the $T\propto a^{-1/2}$ power law. We thus observe that for all three mass points, there is some region of parameter space where reheating is dominated by perturbative processes and the radiation bath realizes the bosonic power law. Lower inflaton masses enable the radiation bath to reach higher temperatures during perturbative reheating.

For a given value of $\mu$, the theory may avoid preheating if the inflaton amplitude at the end of inflation is below the red line in figure 13. As perturbative reheating occurs, the inflaton amplitude decreases due to redshifting in an expanding universe. This redshifting corresponds to traversing downward in the $\mu-\Phi$ parameter space. The temperature of the radiation bath decreases correspondingly along this trajectory. This downward trajectory continues until we reach the blue line and reheating occurs.

Appendix C Collision terms for $s$-channel processes

In this section we derive the total energy transfer rate between two relativistic species at different temperatures, mediated by a massive scalar field. We extend the analysis of references Birrell:2014uka and Adshead:2016xxj to determine the non-equilibrium energy transfer between two sectors at different temperatures. Further, we develop a new procedure to evaluate the phase space integral for a $t$-channel process in appendix D, and confirm explicitly that the contribution from $t$-channel processes is always negligible compared to the $s$-channel process. Finally, for some specific theories where the interactions between the two relativistic species is mediated by a massive scalar field, we find closed-form analytic approximations for the total energy transfer rate.

We start by considering the forward energy transfer for the process

[TABLE]

where 1, 2, 3, and 4 represent the particles participating in the scattering process. The forward collision term for this process is given by

[TABLE]

Here $f_{i}$ is the distribution function for the $i^{th}$ particle, $U$ is the four-velocity of the frame in which we are calculating the collision term, $|\overline{\mathcal{M}}|^{2}$ is the spin-summed scattering amplitude of the process, and $S$ includes any identical particle factors.

This phase space integral simplifies in the center-of-mass (CM) frame, where $U$ is non-trivial. For simplicity we shift to the variables

[TABLE]

In terms of these variables, the Mandelstam variables are $s=p^{2}$, $t=(q-q^{\prime})^{2}/4$ and $u=(q+q^{\prime})^{2}/4$. In the CM frame, $p=(\sqrt{s},0,0,0)$ and consequently $U=\frac{1}{\sqrt{s}}(\sqrt{|\vec{p}|^{2}+s},0,0,|\vec{p}|)$, where $|\vec{p}|$ is the spatial component of $p$ in the frame $U=(1,0,0,0)$. Reference Birrell:2014uka shows that the 12-dimensional phase-space integral of eq. (C) can be reduced to

[TABLE]

where $r$ is the magnitude of $\vec{q}$ and $y=\cos\phi,\theta$ give the direction of $\vec{q}$ with respect to $\vec{p}$, while $r^{\prime}$, $y^{\prime}$, and $\theta^{\prime}$ denote the corresponding quantities for $q^{\prime}$. The spatial and temporal magnitudes of $q$ and $q^{\prime}$ are given by

[TABLE]

For scenarios where the scattering amplitude is a function only of $s$, eq. (C) further simplifies to

[TABLE]

where $s_{0}=\textrm{max}((m_{1}+m_{2})^{2},(m_{3}+m_{4})^{2})$. To evaluate these final integrals, we need to specify the distribution functions. We take particles 1 and 2 to be of species $a$ at a temperature $T_{a}$ and particles 3 and 4 of species $b$ at temperature $T_{b}$. Inserting the corresponding equilibrium distribution functions, we obtain

[TABLE]

Here $\zeta_{a,b}=\pm 1$ depending on whether the respective particles are fermions/bosons and

[TABLE]

Next, we scale out the temperature of the hotter sector, $T_{a}$, by defining

[TABLE]

This isolates the temperature dependence in the integral, which becomes

[TABLE]

The temperature $T_{a}$ enters the integrand only through ${\mathcal{M}}(T_{a}^{2}\tilde{s})$ and $\tilde{m}_{a,b}$.

The total collision term describing net energy transfer is

[TABLE]

where $S^{\prime}=S/(4(2\pi)^{5})$. Given ${\mathcal{M}},\zeta_{a}$ and $\zeta_{b}$, we can now numerically evaluate this integral to obtain the total energy transfer.

Even though we are interested in the regime where all external particles are relativistic, $\tilde{m}_{a,b}\ll 1$, we have explicitly retained $\tilde{m}_{a,b}\neq 0$ in the integral. Retaining finite masses can be important for regulating bosonic scattering amplitudes: for bosons ($\zeta=-1$), the integrand diverges as $\tilde{s},\tilde{p}\rightarrow 0$ when $\tilde{m}=0$, reflecting the zero-momentum singularity in the Bose-Einstein distribution. However, if $|\overline{\mathcal{M}}(\tilde{s})|^{2}$ is finite at $\tilde{s}\rightarrow 0$ then the divergence vanishes after integration over $\tilde{s},\tilde{p}$ and the collision term remains finite as $\tilde{m}_{a,b}\to 0$.555In the scattering collision term, the cancellation of the divergence depends on summing both forward and backward processes; the forward scattering collision term alone retains a logarithmic dependence on $\tilde{m}$. In general one expects thermal self-energies to regulate this behavior when $T_{a,b}\gg m_{a,b}$; see also Evans:2017kti . In our calculations below we can thus freely work in the limit $\tilde{m}_{a,b}\to 0$; however, in numerical work we retain small finite external particle masses for ease of computation.

In the subsequent sections we present analytic estimates for the collision term given in eq. (89) in various theories. We consider the full collision term, i.e., both forward and backward contributions.

C.1 Trilinear scalar couplings

We consider two scalar species, $\chi_{a}$ and $\chi_{b}$, interacting via

[TABLE]

Here $\phi$ is a massive scalar (inflaton) mediator with mass $M_{\phi}$. As both the coupled fields are scalars (and hence bosons), we take $\zeta_{a,b}=-1$ and $S^{\prime}=1/(16(2\pi)^{5})$.

The scattering amplitude for the s-channel process in this theory, for $m_{a,b}\ll M_{\phi}$, is given by

[TABLE]

For $\mu_{a,b}\ll M_{\phi}$ we can approximate the scattering amplitude as Adshead:2016xxj

[TABLE]

where $w=\Gamma_{0b}/\Gamma_{0a}=\mu_{b}^{2}/\mu_{a}^{2}$. To analytically estimate the behavior of $\mathcal{C}_{E}$, we combine the simplified form of the scattering amplitude given in eq. (92) along with approximations $\tilde{M}_{\phi}\ll 1$ and $\tilde{M}_{\phi}\gg 1$ at high and low temperatures respectively.

High temperature limit, $T_{a}\gg M_{\phi}$.

In the high-temperature limit $\tilde{M}_{\phi}\rightarrow 0$, the contribution to the integral in eq. (89) from the $\Theta$ function term in eq. 92 is dwarfed by the contribution from the Dirac delta term. Subsequently, in the high temperature limit we can to good approximation retain only the Dirac delta portion, giving

[TABLE]

To evaluate the above integral we separate it into two domains: $\tilde{p}<0.1$ and $\tilde{p}>0.1$. In the latter region we approximate $\tilde{p}\gg\tilde{M}_{\phi}$ to give

[TABLE]

where

[TABLE]

To evaluate the integral in the region $\tilde{p}<0.1$ we first consider the case where $T_{a},T_{b}\gg M_{\phi}$, allowing us to approximate $\tilde{M}_{\phi}\ll x\ll 1$. Next, note that the integrand $D(\tilde{M}_{\phi}^{2},\tilde{p},x,0,0)$ is peaked near $\tilde{p}\sim\tilde{M}_{\phi}$. Near this peak we can use the approximation $\tilde{p}\ll x$. Assuming the contribution from the peak dominates the integral, we extend the approximation $\tilde{p}\ll x$ to the entire integration range $\tilde{p}\in(0,0.1)$, yielding

[TABLE]

In the case $T_{b}\ll M_{\phi}\ll T_{a}$ the assumptions we used above no longer hold. One can instead use the approximations $\tilde{p},\tilde{M}_{\phi}\ll 1$ along with $e^{\tilde{M}_{\phi}/x}\gg 1$ to simplify the integral and show that its contribution is always dwarfed by the contribution from $\tilde{p}>0.1$. For brevity we do not show the calculations here. Thus we can neglect contributions from $\tilde{p}<0.1$ in eq. (93) when $T_{b}<M_{\phi}$666In fact, even when $\tilde{M}_{\phi}\ll x$, the contributions from the $p<0.1$ integral remains sub-dominant until extremely large temperatures, $T_{a}\geq 10^{5}M_{\phi}$.. We find empirically that using eq. (C.1) for all $T_{b}$ helps improve the agreement between the analytic estimate and the full numerical calculation for $T_{a}$ as low as $T_{a}\sim M_{\phi}$. Thus we approximate the full collision term at high temperatures as

[TABLE]

From the asymptotic behavior at small $x$ we see that ${\mathcal{C}}_{\textrm{high-T}}$ is largely insensitive to $T_{b}$. At extremely small $x$, the logarithmic term is dominant. However, at large temperatures $T\gtrsim 10^{2}M_{\phi}$, as $x$ increases to $x\sim 0.5$, the higher powers of the logarithm take over and enhance the collision term by roughly two orders of magnitude. This enhancement is due to the Bose enhancement of the forward energy transfer. As $x$ further increases towards unity, the backward collision term starts catching up to forward collision term, eventually completely cancelling it at $x=1$.

In the left panel of figure 14, we compare our high temperature estimate with the numerically evaluated collision term. The Bose-Einstein enhancement over the classical Maxwell-Boltzmann result is clearly visible at high temperatures.

Intermediate temperatures, $\bm{T_{a}\leq M_{\phi}}$.

As $\tilde{M}_{\phi}$ begins to exceed unity, the Dirac delta contribution to the matrix element ensures that the integral of eq. (89) has support dominantly at $\tilde{s}=\tilde{M}_{\phi}>1$. However, here the phase space distribution functions, contributing through the factor $D$, are exponentially suppressed. This Boltzmann suppression causes the collision term to fall sharply. In other words, in the intermediate temperature regime the integral receives its dominant contribution from an energy scale much larger than either temperature, which means that to excellent approximation the scattering here can be described using classical statistics.

Using classical statistics, the overall integral over $\tilde{p}$ (eq. (93)) can be performed exactly,

[TABLE]

where $K_{2}$ is the modified Bessel function of second kind. Again, we can see that at small $x$ the collision term becomes insensitive to variations in the colder sector. As the temperatures fall further below the inflaton mass, the collision term becomes Boltzmann-suppressed.

Low temperature limit, $\bm{m_{a,b}\ll T\ll M_{\phi}}$.

In the low-temperature regime, the integral is dominated by off-shell inflaton scattering, described by the Heaviside term in eq. (92). Thus at low temperatures we need to evaluate

[TABLE]

As $D$ is exponentially suppressed at large values of $\tilde{s}$, we can take the upper limit of the $\tilde{s}$ integral to infinity with negligible errors. Both the integrand and the limits of integration thus become independent of temperature, giving

[TABLE]

Again, as required, we find that the energy transfer function becomes insensitive to the colder sector as $x\rightarrow 0$.

Total collision term.

To get a complete analytic estimate of $\mathcal{C}_{E}$ over all temperature ranges we combine the analytic estimates as

[TABLE]

where $\mathcal{C}_{\textrm{high-T}},\mathcal{C}_{\textrm{MB}}$ and $\mathcal{C}_{\textrm{low-T}}$ are given in eq. (C.1), (100) and (C.1). The Heaviside functions ensure that each function contributes only in its region of validity. This function describes the collision term for all temperature ranges as long as the scattering particles remain relativistic ($T_{a,b}\gg m_{a,b}$).

In the left panel of figure 14, we compare the analytic approximations to the energy transfer collision term derived in high-, low-, and intermediate temperature regions with the exact numerical value. For illustrative purposes, each approximation is shown over a range larger than that taken in eq. (C.1). Together these approximations accuratelly describe the behavior of $\mathcal{C}$ in their respective regions. In the right panel of figure 14, we compare the ratio of our analytic approximation, eq. (C.1), to the numerically evaluated collision term using the full matrix element of eq. (91). The largest deviation occurs during the transition from $\mathcal{C}_{\textrm{high-T}}$ to $\mathcal{C}_{\textrm{MB}}$ between $M_{\phi}/4<T_{a}<M_{\phi}$ and is of the order $\sim 50\%$.

C.2 Yukawa couplings

We next consider two Dirac fermions, $\psi_{a}$ and $\psi_{b}$, interacting with a scalar inflaton $\phi$ via

[TABLE]

In this case we have $\zeta_{a,b}=1$ and $S^{\prime}=1/(4(2\pi)^{5})$.

The $s$-channel scattering amplitude in this theory, for $m_{a,b}\ll M_{\phi}$, is given by

[TABLE]

where

[TABLE]

For small $y_{a,b}$ the scattering amplitude can be approximated as Adshead:2016xxj

[TABLE]

where $w=\Gamma_{0,b}/\Gamma_{0,a}=y_{b}^{2}/y_{a}^{2}$.

To estimate $\mathcal{C}_{E}$ analytically, we combine the simplified form of the scattering amplitude given in eq. (C.2) along with the approximations $\tilde{M}_{\phi}\ll 1$ and $\tilde{M}_{\phi}\gg 1$ at high and low temperatures respectively. Moreover, since for fermions the contribution of the distribution functions to the integrand, $D$, is regular at $\tilde{s},\tilde{p}=0$, the limit ${m}_{a,b}=0$ does not need any special attention.

High temperature limit, $\bm{T_{a}\gg M_{\phi}}$.

In the high temperature limit only the Dirac delta term and the second Heaviside theta term in eq. (C.2) contribute to the integral. The collision term then becomes

[TABLE]

Note that as $\tilde{s}\rightarrow 0$ the integrand $D(\tilde{s},\tilde{p},x,0)$ asymptotes to a finite value over all $\tilde{p}$. Thus, we can safely approximate $\tilde{M}_{\phi}=0$ in the integrand, making the integrals independent of $T_{a}$,

[TABLE]

where,

[TABLE]

We can check that at small $x$ we see that $\mathcal{C}_{\textrm{high-T}}$ is in this limit insensitive to $T_{b}$, and at $x=1$ all these functions go to zero as backward energy transfer exactly balances the forward energy transfer. The collision term at very high temperatures in this case is not sensitive to the inflaton mass.

Intermediate regime, $\bm{T_{a}\lesssim M_{\phi}}$.

For $T_{a}\sim M_{\phi}$ the Dirac delta part of the scattering amplitude will dominate the collision term. As discussed in section C.1 above for scalars, the collision term can be well approximated using Maxwell-Boltzmann statistics as the temperature drops below the inflaton mass scale, $T_{a}<M_{\phi}$. Thus, the collision term can be simply written as

[TABLE]

where $K_{2}$ is the modified Bessel function of the second kind.

Low temperature regime, $\bm{m_{a,b}\ll T\ll M_{\phi}}$.

In the low temperature regime, the integral is dominated by the $\Theta(\tilde{M}_{\phi}^{2}-\tilde{s})$ term. Just as for scalars, we can to a good approximation replace $\tilde{M}_{\phi}\rightarrow\infty$ in the limit of integration. This yields

[TABLE]

Total collision term.

We combine the analytic estimates derived above to approximate the collision term for all temperatures as

[TABLE]

where $\mathcal{C}_{\textrm{high-T}},\mathcal{C}_{\textrm{MB}}$, and $\mathcal{C}_{\textrm{low-T}}$ are as described in eq. (109),(112) and (C.2). Eq. (C.2) predicts the collision term for all temperature ranges as long as the particles remain relativistic ($T_{a,b}\gg m_{a,b}$).

In left panel of figure 15, we compare our analytic estimate of the energy transfer collision term with the exact numerical value. Together these approximations accurately model the behavior of $\mathcal{C}_{E}$ in their respective regions. In the right panel of figure 15, we compare the ratio of our analytic approximation, eq. (C.2), to the collision term numerically evaluated (eq. (89)) with the scattering amplitude in eq. (105). The largest deviation occurs during the transition from $\mathcal{C}_{\textrm{high-T}}$ to $\mathcal{C}_{\textrm{MB}}$ between $M_{\phi}/4<T_{a}<M_{\phi}$ and is of the order $\sim 50\%$.

C.3 Axionic couplings to gauge bosons

Next, we consider two (Abelian) gauge bosons interacting with a (pseudoscalar) inflaton $\phi$ via

[TABLE]

In this case we have $\zeta_{a,b}=-1$ and $S^{\prime}=1/(16(2\pi)^{5})$.

The $s$-channel scattering amplitude in this theory, for $m_{a,b}\ll M_{\phi}$, is given by

[TABLE]

where

[TABLE]

For $\Lambda_{a,b}\gg 1$, we approximate the scattering amplitude as Adshead:2016xxj

[TABLE]

where $w=\Gamma_{0b}/\Gamma_{0a}=\Lambda_{a}^{2}/\Lambda_{b}^{2}$.

To estimate $\mathcal{C}_{E}$ analytically, we combine the simplified scattering amplitude given in eq. (118) along with high and low-temperature approximations in the limits $\tilde{M}_{\phi}\ll 1$ and $\tilde{M}_{\phi}\gg 1$ respectively.

High temperature limit, $\bm{T_{a}\gg M_{\phi}}$.

At high temperatures both the Dirac delta contribution and the Heaviside theta term $\propto{\tilde{s}}^{2}$ contribute importantly to the integral. Subsequently we approximate the scattering amplitude as

[TABLE]

In the above equation we have already assumed $\tilde{m}_{a,b}=0$. The first integral on the RHS is exactly the same as that evaluated for scalars at high temperatures. In the second integral, the integrand vanishes as $s\rightarrow 0$, allowing us to freely take $\tilde{M}_{\phi}\approx 0$. This yields

[TABLE]

where

[TABLE]

Here the $Y_{i}$ are defined in eq. (95), (96) and (97).

In the left panel of figure 16 we compare this high temperature estimate with the numerically evaluated collision term.

Intermediate temperatures, $\bm{T_{a}\leq M_{\phi}}$.

In the intermediate regime near $T_{a}\sim M_{\phi}$, integral can again be well approximated with Maxwell-Boltzmann distribution. Thus the collision term can be simply written as

[TABLE]

where $K_{2}$ is the modified Bessel function of the second kind.

Low temperature limit $\bm{m_{a,b}\ll T\ll M_{\phi}}$

In this regime, the $\Theta(\tilde{M}_{\phi}^{2}-\tilde{s})$ term dominates in eq. (89). We can again take $\tilde{M}_{\phi}\rightarrow\infty$ in the limit of integration. This yields

[TABLE]

Total collision term.

We combine the analytic estimates in the following manner,

[TABLE]

where $\mathcal{C}_{\textrm{high-T}},\mathcal{C}_{\textrm{MB}}$, and $\mathcal{C}_{\textrm{low-T}}$ are as described in eq. (C.3),(122) and (C.3). This function describes the collision term for all temperature ranges as long as the scattering particles remain relativistic ($T_{a,b}\gg m_{a,b}$).

In the left panel of figure 16, we compare our analytic estimate derived in the three regions with the numerical value. In the right panel of figure 16, we compare the ratio of our analytic estimate, eq. (124), to the numerically evaluated collision term (eq. (89)) with scattering amplitude given in eq. (116). The largest deviation occurs during transition from $\mathcal{C}_{\textrm{high-T}}$ to $\mathcal{C}_{\textrm{MB}}$ between $M_{\phi}/4<T_{a}<M_{\phi}$ and is of the order $\sim 50\%$.

C.4 Mixed Yukawa and scalar trilinear couplings

In this case we consider a Dirac fermion $\psi$ and a scalar field $\chi$ that interact with the inflaton $\phi$ via

[TABLE]

Note that sector $a$ is not necessarily hotter in this scenario. For this theory $\zeta_{a}=-1$, $\zeta_{b}=+1$, and $S^{\prime}=1/(8(2\pi)^{5})$.

The spin-summed $s$-channel scattering amplitude in this theory, for $m_{a,b}\ll M_{\phi}$, is given by

[TABLE]

where

[TABLE]

For $y,\mu\ll 1$ we approximate the scattering amplitude as Adshead:2016xxj

[TABLE]

where $w=\Gamma_{0b}/\Gamma_{0a}=4y_{b}^{2}M_{\phi}^{2}/{\mu}_{a}^{2}$.

To estimate the behavior of $\mathcal{C}_{E}$ analytically, we combine the simplified form of the scattering amplitude given in eq. (128) along with the high- and low-temperature approximations $\tilde{M}_{\phi}\ll 1$ and $\tilde{M}_{\phi}\gg 1$.

High temperature limit, $\bm{T\gg M_{\phi}}$.

When the temperature of the hotter sector, $T=\textrm{max}(T_{a},T_{b})$, is larger than the inflaton mass, the Dirac-delta term of the scattering amplitude dominates in the integral in eq. (89), giving

[TABLE]

As the scalars $\chi_{a}$ follow BE distribution, $D(\tilde{M}_{\phi}^{2},\tilde{p},x,0)$ has a pole as $\tilde{M}_{\phi}^{2}\rightarrow 0$. We cannot simply approximate $\tilde{M}_{\phi}=0$, and we proceed analogously to the case of section C.1 and split the integral as,

[TABLE]

The integrand in the first term on RHS vanishes for small $\tilde{p}$. Hence, for $\tilde{M}_{\phi}\ll 1$ the first integral can be ignored. Expanding in $\tilde{p}\gg\tilde{M}_{\phi}$ limit, the second integral yields

[TABLE]

As the integrand above diverges only logarithmically as $\tilde{p}\rightarrow 0$, the integral is insensitive to its lower limit, which can be replaced with [math] with negligible error. The collision term at high temperatures can then be simply written as

[TABLE]

where

[TABLE]

Note that the high temperature collision term is approximately insensitive to the colder sector in this case. The collision term’s logarithmic sensitivity on the inflaton mass depends on whether or not the scalars are hotter than the fermions. However, in both the cases the collision term is IR-sensitive due to the dependence on inflaton mass.

In left panels of figure 17 we compare our high temperature estimate with the numerically evaluated collision term.

Intermediate temperatures, $\bm{T\leq M_{\phi}}$.

For temperatures near $T\sim M_{\phi}$ the Dirac delta part of the scattering amplitude dominates the behavior of collision term. As discussed above, in this region, the distribution functions are well approximated by Maxwell-Boltzmann distributions as temperature drops below the inflaton mass scale, $T<M_{\phi}$. The collision term can therefore be simply computed as,

[TABLE]

where $K_{2}$ is the modified Bessel function of second kind. Again, we can see that at small/large $x$ the collision term becomes insensitive to variations in the colder sector.

Low temperature limit, $\bm{m_{\chi,\psi}\ll T\ll M_{\phi}}$.

In the low-temperature regime, the contribution from the Dirac delta part of the matrix element falls below the one from $\Theta(\tilde{M}^{2}_{\phi}-\tilde{s})$ term in the integral in eq. (89). Just like in scalar case, to a good approximation we can replace $\tilde{M}_{\phi}\rightarrow\infty$ in the integral limit. This yields

[TABLE]

Similar to our high temperature estimate, we find that the energy transfer function becomes insensitive to the colder sector as $x\rightarrow 0,\infty$.

Total collision term.

To get an analytic estimate of $\mathcal{C}_{E}$ over all temperature ranges we combine the analytic estimates as

[TABLE]

where $T=\textrm{max}(T_{a},T_{b})$ and $\mathcal{C}_{\textrm{high-T}},\mathcal{C}_{\textrm{MB}}$ and $\mathcal{C}_{\textrm{low-T}}$ are as described in eq. (132), (134) and (135).

In the left panels of figure 17, we compare our analytic estimate with the exact numerical value. In the right panels of figure 17, we compare the ratio of our analytic fit, eq. (C.4), with the numerical evaluation of the full expression (eq. (89)) with scattering amplitude given in eq. (128). The largest deviation occurs during the transition from $\mathcal{C}_{\textrm{high-T}}$ to $\mathcal{C}_{\textrm{MB}}$ between $M_{\phi}/4<T<M_{\phi}$ and is of the order $\sim 50\%$.

Appendix D Collision terms for $t$-channel processes

In this section we calculate the energy transfer rate between two relativistic thermal bath via $t$-channel scatterings and explicitly demonstrate that its contribution will always be dwarfed by the $s$-channel.

For the $2\to 2$ scattering process $1\ +\ 2\ \rightarrow\ 3\ +\ 4$, where 1, 2, 3, and 4 represent the particles participating in the scattering process, the forward collision term for a $t$-channel process is given by

[TABLE]

Here $f_{i}$ is the momentum space distribution function for the $i^{th}$ particle, $U$ is the 4-velocity of the frame in which we are calculating the collision term, $|\overline{\mathcal{M}}|^{2}$ is the spin summed scattering amplitude of the process, $S$ includes the identical particle factors of the process. Unlike the $s$-channel calculation, we do not transform to the center of mass frame, keeping $U=(1,0,0,0)$. To simplify the integral, we perform the following change of variables:

[TABLE]

Correspondingly the Mandelstam variables will be given by $s=(q+q^{\prime})^{2}/4$, $t=p^{2}$ and $u=(q-q^{\prime})^{2}/4$. These substitutions also simplify the phase space element in eq. (D) to,

[TABLE]

Integrating out $p^{\prime}$ and combining the conditions in the step functions we find

[TABLE]

Here $dI_{13}$ and $dI_{24}$ are phase space elements identical up to the different masses of the scattering particles. We first focus on simplifying the $dI_{13}$ phase space element and then swap the masses to get the result for $dI_{24}$.

The delta functions in $dI_{13}$ impose the conditions

[TABLE]

Taking the $z$-axis to be along the spatial component of $p$, the above equations can be re-expressed as

[TABLE]

where $\vec{q}_{xy}$ is the component of $\vec{q}$ along $x-y$ plane. Hence, the delta function in $dI_{13}$ can be transformed to impose conditions on $q^{0}$ and $|\vec{q}_{xy}|$. The corresponding Jacobian for transforming the Dirac delta conditions is

[TABLE]

Subsequently, using the Dirac delta conditions, we integrate over the $|\vec{q}_{xy}|$ and $q^{0}$ to yield,

[TABLE]

where $\theta_{xy}$ denotes the angle between $\vec{q}_{xy}$ and $x$ axis.

The measure $dI_{24}$ has the same form as $dI_{13}$, except with the substitutions $m_{1,3}\to m_{2,4}$ and $q\to q^{\prime}$. Inserting these simplified phase-space elements back in eq. (D), using eq. (D) and making use of the fact that the scattering amplitude only depends on $t=p^{2}$, we get

[TABLE]

In the first step above we have integrated over the phase space angles $\theta_{xy}$ and $\theta_{xy}^{\prime}$. The result in eq. (148) holds for any process whose scattering amplitude is only a function of $t$.

Eq. (148) can be further simplified if particles 1 and 3 are identical (henceforth, particle $a$), and and similarly particle 2 and 4 (particle $b$). With these assumptions, $F_{a,b}=F_{13,24}$ simplifies to

[TABLE]

where $\beta_{a,b}=\sqrt{1-\frac{4m_{a,b}^{2}}{p^{2}}}$. Furthermore, if particles $a$ and $b$ are thermally distributed with temperatures $T_{a}$ and $T_{b}$, respectively, the integral in $F_{a,b}$ can be evaluated exactly to yield

[TABLE]

where $\zeta_{a,b}=\pm 1$, depending on whether particle $a,b$ is a fermion or boson respectively. Substituting this result into the forward collision term and converting the integral over $p^{0}$ to an integral over $t$ using

[TABLE]

we find

[TABLE]

In the case where sectors $a$ and $b$ have thermalized, i.e. $T_{a}=T_{b}$, the forward collision term is identically zero, as expected from a $t$-channel process. The backward collision term has the same form, except with an overall minus sign. The total collision term, $\mathcal{C}_{E}=\mathcal{C}_{E}^{f}-\mathcal{C}_{E}^{b}$, is therefore twice the forward rate. Finally, we make the integrand dimensionless by pulling overall factors of $T_{a}^{5}$ to give

[TABLE]

where $S^{\prime}=S/4/(2\pi)^{5}$, $x=T_{b}/T_{a}$ and a tilde over a variable indicates that it has been made dimensionless by dividing by the appropriate power of $T_{a}$.

D.1 Trilinear scalar couplings

In this section we evaluate the $t$-channel contribution to the collision term and compare it with the $s$-channel rate for the theory in which the two sectors interact via

[TABLE]

The $t$-channel scattering amplitude for this theory is given by

[TABLE]

In figure 18 we compare the resulting $t$-channel contribution to the collision term with the $s$-channel contribution for this theory at fixed $x=T_{b}/T_{a}$. The $s$-channel contribution always dominates over the $t$-channel contribution, both below the resonance as well as above.

At low temperatures, $T_{a}\ll M_{\phi}$, we can approximate the scattering amplitude as a constant just like we did for $s$-channel,

[TABLE]

The temperature dependence in this region can then be isolated,

[TABLE]

The low temperature behavior derived above is similar to the one we had in $s$-channel case, eq. (C.1). However, unlike $f(x)$ in eq. (C.1), $f_{t}(x)$ decreases as $x\rightarrow 0$, rather than approaching a constant; this reflects the suppressed probability of finding an initial scattering particle in the colder sector as the thermal abundance drops. We find the $t$-channel rate to be suppressed by a factor of six relative to the $s$-channel rate even when $x\rightarrow 1$. In a similar manner one can also show that the $t$-channel contribution remains subdominant for other theories considered in appendix C. Thus our neglect of the $t$-channel contribution to $\mathcal{C}_{E}$ is justified.