Chilly Dark Sectors and Asymmetric Reheating

In a broad class of theories, the relic abundance of dark matter is determined by interactions internal to a thermalized dark sector, with no direct involvement of the Standard Model (SM). We point out that these theories raise an immediate cosmological question: how was the dark sector initially populated in the early universe? Motivated in part by the difficulty of accommodating large amounts of entropy carried in dark radiation with cosmic microwave background measurements of the effective number of relativistic species at recombination, $N_{\mathrm{eff}}$, we aim to establish which admissible cosmological histories can populate a thermal dark sector that never reaches thermal equilibrium with the SM. The minimal cosmological origin for such a dark sector is asymmetric reheating, when the same mechanism that populates the SM in the early universe also populates the dark sector at a lower temperature. Here we demonstrate that the resulting inevitable inflaton-mediated scattering between the dark sector and the SM can wash out a would-be temperature asymmetry, and establish the regions of parameter space where temperature asymmetries can be generated in minimal reheating scenarios. Thus obtaining a temperature asymmetry of a given size either restricts possible inflaton masses and couplings or necessitates a non-minimal cosmology for one or both sectors. As a side benefit, we develop techniques for evaluating collision terms in the relativistic Boltzmann equation when the full dependence on Bose-Einstein or Fermi-Dirac phase space distributions must be retained, and present several new results on relativistic thermal averages in an appendix.


Introduction
Although we have ample gravitational evidence for the existence of some form of dark matter (DM) constituting 26% of the energy budget of our universe [12], its detailed nature and properties remain one of the greatest outstanding mysteries in particle physics and cosmology. The lack of observational evidence to date of traditional weakly-interacting massive particle (WIMP) dark matter candidates in direct detection, indirect detection, and collider searches has helped to motivate a recent explosion of interest in a much broader range of dark matter theories that exhibit a wide variety of interesting and nontraditional signals. In many of these models, the DM relic abundance is chiefly determined by selfinteractions among a set of fields that live in a thermalized dark sector, with little to no involvement of the Standard Model (SM). Such self-interacting thermal dark sectors can provide novel solutions to long-standing puzzles in particle physics or astrophysics, can yield novel signals, and in more generality represent a generic possibility for the physics of the invisible universe.
Models that invoke a thermalized dark sector immediately raise a cosmological question: how was this dark sector populated in the early universe? One minimal answer is to produce the dark sector through a (small) interaction with the SM, e.g., the frequently-considered kinetic mixing between a dark gauge group and SM hypercharge [1][2][3]. Such an interaction between the hidden sector (HS) and the SM is often required for other reasons, for instance, to enable the deposition of the dark sector's entropy into the SM plasma prior to Big Bang Nucleosynthesis (BBN). A thermal dark sector can be produced when the coupling between the HS and the SM is sufficiently large to bring the two sectors into thermal equilibrium at some temperature in the early universe. This mechanism for populating a dark sector has several attractive features, and in particular is relatively insensitive to the as-yet-unknown evolution of the early universe.
On the other hand, once we specify that a dark sector was once in thermal equilibrium with the SM, we limit the degree to which its temperature T HS can subsequently differ from the temperature of the SM plasma, T SM . This causes difficulties for many models of hidden sector dark matter that rely critically on having a large asymmetry between T HS and T SM in order to prevent hot relics, such as mirror neutrinos or dark gauge bosons, from contributing at unacceptable levels to the expansion of the universe (see, for example, [4][5][6][7][8][9][10][11]). From Planck's current constraints on the effective number of neutrinos [12], we can limit the total number of degrees of freedom in a hidden sector that was once thermalized with the SM to be g * S,HS < O(10), if it contains hot relics that contribute as neutrinos to the expansion of the universe during the formation of the CMB. This constraint already cannot be met by many of the above models, and will only become more stringent with time: the projected capabilities of the next generation of CMB experiments [13][14][15] represent an order of magnitude improvement over Planck's current sensitivity to free-streaming hot relics. While this is one of the most frequently-invoked motivations to consider a decoupled hidden sector, such hidden sectors are a generic possibility for the origin of dark matter, and can lead to qualitatively novel signals [7,[16][17][18].
We can then immediately identify two equally minimal mechanisms to populate a thermal hidden sector that was never in kinetic equilibrium with the SM. First, the HS can have a small interaction with the SM that never enters equilibrium. In this case, the SM plasma will continually inject energy into the hidden sector through out-of-equilibrium processes over a range of temperatures from some initial value T max of order the reheating temperature, until some final temperature T end when the interaction becomes negligible [19][20][21][22]. Second, the physical process that populated the SM itself in the early universe at T max can also simultaneously produce the HS at a lower temperature; we will refer to this mechanism as "asymmetric reheating" [5,6].
In the first case, the dark sector will undergo non-adiabatic evolution during the period of continual energy injection. Depending on the size of the temperature asymmetry and the relation of T end to various scales in the hidden sector, the entropy injection can alter the predictions of hidden sector dark matter theories in interesting ways, e.g., by suppressing indirect detection signals (see also [23]).
The second option, asymmetric reheating, requires that some physics beyond the SM couple to both the SM and the HS. This implies the inevitable existence of reactions that transfer energy between the two sectors, and thus can, if strong enough, erase a would-be temperature asymmetry.
In this paper, we study in detail the energy transfer between two sectors in minimal single-field models of reheating, and identify the regions where potentially resonant inflaton or reheaton exchange efficiently equilibrates the two sectors. Along the way we obtain new results on relativistic thermal averages for particles obeying Bose-Einstein and Fermi-Dirac distributions. Thus we will demonstrate that a sufficiently large temperature asymmetry between two sectors either (1) restricts the possible combinations of reheating temperatures, inflaton masses, and inflaton coupling structures; (2) requires a non-minimal mechanism for reheating, e.g. with multiple fields [6]; or (3) requires an alternative cosmological history for one or both sectors, e.g. with late-time asymmetric entropy release into the SM alone [24].
Because reheating occurs on subhorizon scales that subsequently undergo non-linear evolution, in simple inflationary scenarios the reheating epoch itself generically leaves no direct imprint on scales relevant for cosmology, such as could be observed in the CMB or the distribution of galaxies. 1 Thus information about the reheating phase of the universe is generically very difficult to obtain directly. Aside from model-dependent effects such as the production of gravitational waves, magnetic fields, and primordial black holes (see, for example [27][28][29][30] and references therein), the primary impact of the reheating phase on cosmology enters through the unknown evolution of the scale factor and Hubble rate during the entire period that connects the end of inflation to the period of radiation domination. This epoch of expansion changes how physical length scales in the universe today are related to length scales during inflation, and can thereby alter the precise location on the inflationary potential where the observable fluctuations in the CMB were produced. Therefore, given a specific model for inflation which makes a specific prediction for the spectrum of fluctuations, indirect constraints on reheating can be placed by constraining the expansion history during the reheating period [31][32][33][34][35][36]. However, caution is necessary when translating these into constraints on microphysical particle physics theories, as the detailed properties of the processes responsible for reheating can significantly alter its duration [37][38][39][40][41][42][43][44][45]. By contrast, the approach we take here is to remain as independent of specific inflationary models as possible, and our results are largely insensitive to the details of the potential probed during inflation, or the matter content of the SM and hidden sectors. Connecting the physics of the reheating epoch to admissible cosmologies of dark matter theories helps expand future avenues to further pin down the properties of this phase of our universe's evolution.
We begin by remarking on the currently forecast sensitivity to ∆N eff from future CMB experiments, translating the potential sensitivity into projected constraints on the field content in oncethermalized hidden sectors that contain a free-streaming hot relic in section 2. In section 3 we discuss single-field reheating and describe the setup for our calculations in section 4, which maps out the regions where asymmetric reheating can or cannot yield a temperature asymmetry for a variety of different minimal models. Section 5 contains our conclusions. New analytic and numerical results on relativistic thermal averages are collected in the appendices A and B. In appendix C we present a brief overview of preheating. We work in units where = c = k B = 1, but we explicitly retain the reduced Planck mass M 2 Pl = 1/(8πG).

Future CMB sensitivity to dark radiation in thermal hidden sectors
A long-standing motivation for asymmetric reheating has been the need to prevent stable dark radiation in a dark sector from contributing at unacceptable levels to observables sensitive to the expansion of the universe at early times. The most important such observables are the light element abundances from BBN, which test the number of relativistic species at temperatures T ∼ O(MeV), and the CMB, which is sensitive to the number of relativistic species at temperatures T ∼ O(eV). While the need to reconcile the existence of dark radiation with existing BBN and CMB measurements is hardly the only motivation to consider asymmetric reheating, it arises in a wide variety of models, and is connected to some of the most direct experimental signatures in several models of hidden sector dark matter. For this reason it is worth mentioning that the anticipated sensitivities of future CMB experiments, in particular those currently being discussed for CMB Stage-IV, improve on Planck's current sensitivity to dark radiation by an order of magnitude. This improvement has profound implications for dark sectors with stable relativistic relics that were ever in thermal equilibrium with the SM, as we briefly review here. The sensitivity of CMB measurements to the number of effective relativistic degrees of freedom at recombination is conventionally quoted in terms of the effective number of neutrinos, N eff . Additional relativistic degrees of freedom contribute a shift in the effective number of neutrinos by an amount where g * S,HS , g * ,HS are defined (in analogy to the usual definitions of g * S , g * in the SM) within the dark sector with respect to the temperature of the dark plasma T HS , For simplicity we assume a single dark radiation species, so, denoting by IR the value pertaining during the formation of the CMB, g IR * ,HS = g IR * S,HS , while g IR * S,SM = 3.9. Current forecasts for the capabilities of next-generation CMB experiments, in particular CMB Stage-IV, anticipate an exciting level of sensitivity to ∆N eff , with a projected 68% CL uncertainty in the range σ(N eff ) ∼ 0.015 − 0.03 [13][14][15]. To highlight the potential impact on models of hidden sector dark matter, we briefly illustrate the implications of a SM-like measurement of N eff at this level of precision for the total field content of simple hidden sectors.
Specifically, we are interested in internally-thermalized hidden sectors that contain stable dark state(s) that may include dark radiation as well as dark matter. When the branching fractions of hidden sector states to SM states are negligible, all (or almost all) of the entropy contained in a hidden sector is ultimately carried by the lightest stable dark state. In the generic situation where this lightest stable dark state is a hot relic that is still relativistic during the formation of the CMB, the total entropy in the hidden sector governs the magnitude of the dark radiation's contribution to N eff . If such a dark sector was ever in thermal equilibrium with the SM at some point in its cosmic history, then the temperature at thermal decoupling T D sets a lower bound on the amount of entropy in the hidden sector post-decoupling, and therefore on the dark radiation's ensuing contribution to N eff .
For simplicity, we will take the dark radiation to be free-streaming, and work in the limit of instantaneous decoupling. Under these assumptions, a constraint on N eff immediately translates into a constraint on the total effective number of degrees of freedom in a hidden sector that was last in equilibrium with the SM at temperature T D .
In figure 1, we show the projected values of g * S,HS (T D ) that would be consistent at 2σ with a future measurement of N eff at its SM value in CMB Stage-IV, for σ(N eff ) = 0.015 and σ(N eff ) = 0.03, Figure 1. Projected values of gS, * HS(TD) allowed at 2σ by future CMB experiments such as CMB Stage-IV in a hidden sector that contains free-streaming dark radiation, as a function of the HS-SM kinetic decoupling temperature TD. Allowed regions are shaded, using σ(N eff ) = 0.03 (solid boundary) and σ(N eff ) = 0.015 (dotted boundary), for real scalar (red), Majorana fermion (blue), or massless vector boson (cyan) dark radiation.
and for three different species of dark radiation. With σ(N eff ) = 0.015, CMB Stage-IV will be able to exclude a single massless vector or fermionic dark radiation species at 2σ, but there will remain a small window for real scalars that decouple from the SM at T D > 51 GeV. With σ(N eff ) = 0.03, there is some allowed parameter space remaining for single massless vector or fermionic dark radiation species at sufficiently large T D , while real scalars decoupling from the SM at any temperature above the chiral phase transition are allowed. Note that none of these scenarios allow more than one hidden sector species to be relativistic at T D , which would place stringent constraints on the allowed masses of possible non-relativistic relics-i.e., dark matter-in the hidden sector. These results assume freestreaming dark radiation and instantaneous decoupling, and are thus are only illustrative. In specific models of once-thermalized dark sectors, model-dependent effects such as the distortion of the phase space distribution function of the dark radiation during decoupling [46] or the interactions of dark radiation with itself and with dark matter [47][48][49] can often contribute corrections that will become increasingly important given the unprecedented sensitivity being discussed for CMB-IV and related experiments.
One possible way to relax this constraint on the field content of such hidden sectors is to add new degrees of freedom coupling to the SM in the ultraviolet. Adding the field content of the Minimal Supersymmetric Standard Model (MSSM) below T D relaxes the forecast constraint sufficiently to allow, e.g., a tiny hidden sector with g D * S,HS < 3.004 and g IR * S,HS = 1, taking σ(N eff ) = 0.015. While this hidden sector can now in principle allow for multiple species, it cannot accommodate an entire superfield, suggesting that the scale of supersymmetry breaking must be higher in the dark sector than in the SM.
A far more flexible way to relax the constraints on dark radiation is simply not to let the hidden sector ever thermalize with the SM. If T HS is sufficiently small compared to T SM , models with dark radiation (such as, for example, [7][8][9][10]) can be made compatible with even next-generation measurements of N eff and related observables. Given this strong motivation to consider chilly dark sectors, we now turn to discussing how they can be obtained during the reheating epoch, and the non-trivial demands that are therefore placed on the possible reheating histories of our universe.
In the simplest picture of reheating, dubbed the elementary theory of reheating [39,[50][51][52], reheating occurs through the perturbative decays of the inflaton as it oscillates about the minimum of its potential. This process occurs on a time scale set by the inflaton width, t ∼ 1/Γ φ . During this time, the energy stored in the inflaton field dominates the energy density of the universe, resulting in a matter-dominated expansion. The decays of the inflaton drain the energy from the condensate and damp the oscillations, leaving a relativistic bath of daughter particles. Reheating thus ignites the subsequent radiation-dominated era as relativistic particles are produced and subsequently thermalize. This process can be described by the following Boltzmann equations for the energy densities of the inflaton and the radiation bath, The energy density stored in radiation rapidly increases from zero to a maximum value and then declines as T ∼ a −3/8 , where T can always be formally defined through even before the radiation bath has attained internal thermal equilibrium. Here g * is the number of (effective) degrees of freedom in the radiation bath. The reheat temperature T RH is defined as the temperature where the radiation comes to dominate over the inflaton, at which point the universe enters a standard radiation-dominated phase. A simple estimate of T RH can be obtained by comparing the inflaton decay width Γ φ to the Hubble rate. When the decay width is comparable to the Hubble rate, Γ φ ∼ H, the decay process proceeds efficiently. In the radiation-dominated era starting from T RH , the Hubble parameter is given by which gives the classic expression for T RH as a function of the inflaton width, While the maximum value attained by the radiation energy density depends on the energy scale at the end of inflation, the reheat temperature does not, and neither quantity is directly dependent on the inflaton mass.
This classic picture of reheating has made some reasonable but simplifying assumptions. In particular, • Reheating is driven by the oscillations and decay of a single field φ. This is the minimal possibility, and an assumption we will maintain throughout this paper. For simplicity, we will refer to this field as the inflaton, although it could be a separate "reheaton" field. 2 • The region of the potential seen by the inflaton during its oscillations is purely quadratic. Departures from pure quadratic behavior will alter the effective equation of state of the oscillating inflaton field, which can be parameterized by a field-dependent w(φ) [54]. A non-zero w(φ) alters the expansion rate of the universe during reheating, which will in turn affect the resulting T RH for fixed inflaton mass and width.
• The dissipation of energy from the inflaton condensate occurs through perturbative two-body decays. Depending on the form of the interaction between the inflaton and its daughter fields as well as the details of the interaction between the daughter fields themselves, this assumption may fail during a large fraction of the reheating process. In particular, there may be significant particle production through a non-perturbative "preheating" process involving parametric resonance if the couplings between the inflaton and matter are sufficiently strong [38,39,41]. Preheating can result in a far more efficient transfer of energy from the inflaton condensate to the radiation bath than the tree-level width implies, and will generically increase the ultimate T RH for fixed inflaton mass and width. Other possible effects such as thermal blocking and Landau damping can arise when we relax the additional implicit assumption that collective effects in the radiation bath can be largely neglected [44,[55][56][57].
For our studies in this paper, which serves as the first general study of multi-sector reheating, we will adopt the following reasonable but simplifying assumptions as a starting point to describe the radiation bath: • All states of interest in the radiation sector have effective masses m M φ , T RH . In other words, we consider radiation baths composed of particles with masses small compared to other mass scales in the problem. This covers the majority of the parameter space and allows for simpler analytic approximations to various quantities we compute. While the production of heavy particles (m T RH ) during the reheating process can be of interest for theories of dark matter and baryogenesis (for example, [42,43,[58][59][60][61]), this possibility is tangential to our study. We leave the careful study of the effects of thermal masses to future work.
• Radiation sectors have attained internal thermal equilibrium by T RH . Attaining an equilibrium phase space distribution requires number-changing as well as momentum-changing interactions to be efficient [62,63]. The timescale for the thermalization of a radiation sector is necessarily dependent on the strength and structure of the interactions within that sector, and in some cases can be longer than Γ −1 φ . For simplicity of presentation, we will assume near-instantaneous thermalization in writing reaction rates, but all of our quantitative conclusions use these reaction rates evaluated at T RH and later, where it will be often be a good assumption that local thermal equilibrium has been attained within each sector. In some bosonic cases that we will see below, reheating can be very quick, which can render this assumption questionable. However, in these cases, using a thermal distribution to describe the radiation is conservative for evaluating the degree of inflaton-mediated energy transfer between two radiation sectors: a less-equilibrated sector has a greater fraction of its particles concentrated in the region of phase space where resonant (and Bose-enhanced) inflaton-mediated scattering can occur.
• The states in the radiation bath that dominate the energy dissipation from the inflaton are welldescribed as single (quasi-)particles. This is another manifestation of the implicit assumption that collective effects in the radiation bath can be largely neglected, and ensures that the description of reheating in terms of Boltzmann equations is valid.
In contrast to many previous studies, however, we use full Bose-Einstein or Fermi-Dirac distributions f (p) to describe particles in thermal equilibrium. This is critical, as the thermally averaged reaction rates for particles with different statistics can have differences of multiple orders of magnitude. We give many new results on relativistic thermal averages in the Appendices.
These simplifying assumptions allow us to expose the essential physics without being overwhelmed by model-dependent details. While departures from these assumptions will change the details of our conclusions, we expect the general picture established in the following sections to be preserved.

Thermal history of two-sector reheating
We now develop the Boltzmann equations for two-sector reheating. Our primary aim is to investigate the possibility of thermalization between the SM and a hidden sector in the limit where their only non-gravitational interaction is the inevitable scattering via inflaton exchange that follows from the assumption that both sectors are reheated through couplings to the same field. We take the inflaton to have partial decay widths into two otherwise decoupled sectors. In general the two sectors can have different temperatures, which we label as T a and T b ; where there is a distinction, we will take a to represent the SM and b to represent the HS.
Two-sector reheating can be described by the coupled Boltzmann equations where Γ a + Γ b = Γ φ , and the collision term C describes the energy transfer between sectors that results from inflaton-mediated 2 ↔ 2 scatterings. We consider the simple and illustrative case where the inflaton couples to a single species of particle in each sector. In this case, the leading-order scattering matrix elements are completely determined by the zero-temperature partial widths Γ 0 a,b and M φ , and have contributions in both s and t channels, where the two SM particles are labeled 1, 2, the two HS particles are labeled 3, 4, and n-particle phase space is defined to include the internal degrees of freedom, The inflaton decay rate Γ φ also depends on the phase space distributions f (p) within each radiation sector: e.g., where Γ 0 a is the zero-temperature partial width, and in the second line we have taken particles 1 and 2 to be described by a common distribution f (E). Eqs. (3.6)-(3.8) are thus a full description of the system only when both sectors have attained internal thermal equilibrium, and the f i can be taken to be the equilibrium distributions in each sector; otherwise one must also take into account the differential evolution of the f i .
This brings us to another point: in eqs. (3.6)-(3.8), we have also neglected the production of finite-momentum inflaton quanta through their interactions with the radiation baths. In general these quanta will be present, and should be included in a complete description of the system. We have omitted these interactions for simplicity, as even in the cases where inflaton quanta are efficiently populated, they will not affect our results, as we argue below.
We are interested in establishing the regime where the energy transfer between sectors described by eqs. (3.9)-(3.10) reaches equilibrium. To do so, we consider the energy transfer rate Γ E at T RH , by which time we assume that the two sectors have indeed attained internal thermal equilibrium, and compare it to H(T RH ). If the energy transfer rate is in equilibrium for temperatures T < T RH , when inflaton decays cease to be important, then any initial temperature asymmetry will have been erased, and no subsequent asymmetry will be generated.
For this purpose, it suffices to consider s-channel scattering. The energy transfer in an s-channel process is maximal, the full incoming energy that participates in the reaction: ∆E s = E 1 + E 2 . The energy transfer in a t-channel exchange is ∆E t = E 1 − E 2 , and is comparable to the s-channel transfer only when there is a substantial difference in temperatures between sectors. While t-channel scattering rates can be enhanced in the forward region when M φ T , the energy transfer in this region is small. The rates for s-channel scattering are also resonantly enhanced when T ∼ M φ , where M φ is the mass of the inflaton. Therefore, s-channel processes will give the dominant contribution to the energy transfer, which simplifies our analysis. For the remainder of our discussion we specialize to s-channel inflaton exchange, as the inclusion of t-channel processes will not affect our results.

Determining T RH
As in the single sector case of eq. (3.5), we estimate the end of reheating as occurring when Γ φ (T a , T b ) = H(T a , T b ), which defines a reheat temperature in each sector. Defining this determines T RH,a in terms of the temperature asymmetry ξ through where we have definedg For the inflaton width, we use eq. (3.13) with equilibrium phase space distributions in each sector, giving where Γ 0 a,b is the zero-temperature partial width, and the upper (lower) sign holds for bosons (fermions). In figure 2 we show how quantum statistics affect the determination of T RH for a single sector. Pauli blocking and especially Bose enhancement have a major impact for T RH M φ /2, where T RH can differ by orders of magnitude from the zero-temperature estimate. The maximum achievable T RH corresponds to a near-instantaneous transfer of energy from the oscillating scalar field to the radiation bath, The exponential enhancement of Γ φ (T ) in the case of Bose statistics can result in extremely efficient inflaton decays for temperatures sufficiently bigger than M φ , saturating this upper bound.
To test whether inflaton scattering can thermalize the two sectors, we assume that the two sectors have thermalized, i.e., we take ξ = 1, and then check to see whether the resulting energy transfer rate can exceed the Hubble rate at T RH . If this criterion is met, then the scenario is self-consistent, and thermalization can be achieved. To be quantitative, we will need to assume a value for g * ,b . We take for definiteness g * ,a = g * ,b , but our main results are not sensitive to the details of this choice. It is largely thanks to the large SM value of g * ,a that our results are insensitive to the possible production of inflaton quanta: in the regions where thermalization occurs and we expect inflaton quanta to be copiously produced, the total amount of entropy carried by these inflaton quanta will be an unimportant fraction of the entropy in the radiation bath. Studies of low-scale reheating, where g φ g * ,a (T RH ) may no longer hold, would need to more carefully revisit this point.

Energy transfer rate
We now turn to the collision term in eq. (3.9). Consider the forward process in the s-channel collision term of eq. (3.9), the backward process can be written analogously as n 2 In thermal equilibrium, the forward and backward processes are equal.
A convenient expression for the fractional energy transfer rate is thus whereρ is the (average) energy density carried by the incoming particle. When the particles that interact with the inflaton in both sectors follow the same type of quantum statistics, we have simplȳ (3.20) When sectors a and b follow different quantum statistics, we definē in order to maintain manifest the invariance of the energy transfer rate under exchanging sectors a and b when both sectors are in equilibrium. The evaluation of the integrals involved in the thermal average in eq. (3.18) with Fermi-Dirac or Bose-Einstein distributions is non-trivial, yet can be crucial for capturing the physical behavior of relativistic scatterings. The majority of existing studies on thermally averaged scattering crosssections are in the context of non-relativistic WIMP dark matter near its thermal freeze-out. In this non-relativistic context, the phase space distribution is well approximated by the Maxwell-Boltzmann distribution and phase space blocking or enhancement factors can be neglected, both of which greatly simplify the integration. In the present case, we have no reason to expect either of these approximations to hold. Here we adopt and extend the integration strategies originally developed in [64,65] to study neutrino decoupling in the early universe. We conduct full computations with both Fermi-Dirac and Bose-Einstein distributions for various scattering cross-sections of interest. The s-channel resonance enhancement is found to be dramatic as long as there is sizeable phase space where the center-of-mass energy of the particle collision is around M φ (i.e. T is not too much lower than M φ ). This has a significant impact on thermalization between the two sectors. To the best of our knowledge, many of the results presented here based on full treatment of relativistic scatterings are new to the literature, and highlight the critical importance of treating quantum statistics accurately at high temperature, particularly in the case of scalars. In order to not distract from the physics, we present the somewhat gory details of our calculations in appendix A.

Inflaton-mediated thermalization in simple models of two-sector reheating
We are now ready to use the set-up of section 3.1 to conduct quantitative studies of several simple representative models of two-sector reheating that cover various assignments of quantum numbers to the inflaton and its decay products. For each example model we consider, we derive an approximate analytic expression for the energy transfer rate Γ E , then check to see in what region of parameter space the self-consistent thermalization condition can be satisfied for T ≤ T RH . In most cases, Γ E (T ) decreases faster with T than does H(T ). This means that if thermalization does not happen by T RH , it will not occur at lower temperatures. However, there are two exceptions, as we will discuss in more detail in the following subsections. First, the resonant structure of the energy transfer rate becomes important for T ∼ 2M φ /5, which can allow thermalization even if Γ E (T RH ) < H(T RH ). Second, in the case where the inflaton and its decay products in both sectors are all scalars, Γ E (T ) ∝ T for T M φ , making thermalization possible at late times.
We consider fixed example values of M φ , and for simplicity, we assume all external masses m M φ , T RH to be equal. In practice, we take a small but non-zero value of 10 GeV for these masses to avoid numerical artifacts, unless otherwise indicated. In our scan of parameter space we parametrize the models using the variables where we generally assume k ≤ 1. The parameter w represents the (zero-temperature) inflaton decay width relative to its mass and is related to T RH via eq. (3.15), while the parameter k characterizes the ratio of zero-temperature widths to SM and HS states. To better illustrate the physics in variables more closely related to observables, we will present our results in the plane of (k, T RH ), where T RH has a one-to-one correspondence with w for a given M φ . To be as model-independent as possible, we consider a general range of possible values for w. We restrict w 0.01 for perturbativity, which can limit the achievable T RH in specific models. Among the models we consider is the case where the inflaton has axion-like couplings to gauge bosons and fermions via dimension-five operators, in which case the validity of the effective field theory (EFT) imposes limits on the achievable values of w for a given M φ .
We derive analytic approximate formulae for the thermally averaged energy transfer rate and apply them in our numerical scan. The details of these approximations are presented in appendix B. Our analytic formulae agree well with the results from exact numerical evaluation from T M φ down to T ∼ m, where m is the external particle mass. For temperatures below T ∼ m, we expect that the Boltzmann suppression of the external states will render inflaton-mediated scattering unimportant on cosmological timescales.

Scalar inflaton coupling to scalar pairs in both sectors
We first consider the model where the inflaton is a CP-even scalar which couples to scalars in both sectors via dimensionful trilinear couplings. The relevant Lagrangian in this case is where we have not explicitly written the mass terms for the fields, taking m Sa = m S b ≡ m S M φ . In this minimal model we take the S i to be real singlets; extensions to complex scalars (such as the SM Higgs) are trivial. The zero-temperature partial widths are thus We can express µ a in terms of k and w as (4.5) The amplitude for S a S a ↔ S b S b scattering is given by (4.6) The derivation of the energy transfer rate for this case can be found in appendix B.1. More than any other case, the potentially resonant scattering of scalars to scalars demonstrates the importance of retaining full dependence on the quantum statistics in evaluating the thermal average. The energy transfer rate is plotted as a function of temperature in the left panel of figure 3 (see also figure 10) for two different values of the width. These values correspond to the boundaries of the middle panel of figure 4 at k = 0.5. Strikingly, the rate asymptotes to a constant at high temperature, thanks to the non-vanishing overlap of the zero-momentum singularity in the Bose-Einstein distribution with the pole in the scattering amplitude.
Scattering in the scalar case has the unusual feature that, as it is controlled by super-renormalizable couplings µ a,b , it falls off more slowly than does the Hubble rate as the temperature decreases. Below the scale of the inflaton mass, we can estimate while the Hubble rate decreases as H ∝ T 2 during the radiation-dominated era. We compare the evolution of the two rates as a function of temperature for one choice of parameters in the left panel of figure 3. Thus in this particular case, thermalization can be controlled by mass scales in the infrared: in our simple models, the lowest temperature where thermalization is possible is set by the mass scale of the external states, T ∼ m S . In this model, when we check for thermalization, we check to see whether Γ E (T ) > H(T ) for any T in the range (m S , T RH ). We continue to take ξ = 1 in determining T RH , even though thermalization may not occur until much lower temperatures. This is conservative, as keeping ξ < 1 earlier in the radiation-dominated era would decrease the Hubble rate relative to Γ E , thereby making thermalization easier. The region of parameter space where thermalization can be achieved in this model is shown in figure 4 for three reference choices of the inflaton mass. The maximum T RH shown in these plots is the maximum allowed by energy conservation, i.e., is the value realized by instantaneous reheating; the tree-level couplings are still perturbative.
As is evident, thermalization preferentially occurs for high values of T RH , despite the IR sensitivity in this model. This can be understood simply from dimensional analysis. Using the estimate of eq.  (4.7) and T RH ∼ Γ φ M pl and dropping numerical constants, the criterion for thermalization becomes (4.8) Thus we can make two important observations: first, that late-time thermalization will be more important in theories that yield high T RH , and second, that thermalization is easiest when T RH > M φ . This estimate neglects the possible Bose enhancement in determining T RH (and in scattering at T RH ); this enhancement increases T RH for fixed M φ , Γ 0 φ , further reinforcing these two observations. In figure 4, we observe that for M φ = 10 3 GeV, thermalization occurs near the resonance, and consequently the thermalizing parameter space does not depend on the value of the external m S . The energy transfer rate drops sharply below the resonance, and there are not enough decades of temperature in the remaining range m S < T < M φ to allow the energy transfer rate to come back into equilibrium. Meanwhile for M φ = 10 13 GeV, m S T RH , and if thermalization occurs, it does so in the IR. The value M φ = 10 8 GeV realizes an intermediate case, with lower values of T RH yielding thermalization at T ∼ M φ , while higher values of T RH enable late-time thermalization.
The regions where thermalization can occur are characterized by relatively strong couplings of the inflaton to matter. This is precisely the region of parameter space where we expect non-perturbative particle production -i.e., preheating -to potentially be important. For sufficiently large couplings of the radiation bath to the inflaton, the time-dependence of the oscillating inflaton background will lead to particle production. We would typically expect non-perturbative particle production to become important when the background-induced time variation of the frequency for the radiation modes is no longer small relative to the frequency:ω ω 2 , (4.9) which, for our trilinear model, is parametrically the condition that µΦ M 2 φ . Note that this condition depends not only on the inflaton mass and couplings to radiation, but also the amplitude of the collective oscillations. Thus specifying Γ φ and M φ is not a priori enough information to determine the importance of preheating. However, the Hubble parameter during reheating is bounded: on one hand, reheating occurs for H Γ φ , and on the other hand, we must have H M φ if φ is not to inflate (we assume, as always, that the scalar field dominates the energy density of the universe and that the region of the potential explored by φ is quadratic). This translates parametrically into a finite possible range for Φ, As we review in appendix C.2, the regions in parameter space where preheating is relevant for scalar trilinear couplings are largely determined by the value of the parameter q ≡ 2µΦ/M 2 φ . The violation of adiabaticity of eq. (4.9) corresponds to the condition that q > 1, which defines the "broad resonance" regime. Here preheating can yield a very efficient transfer of energy from the inflaton condensate to the radiation bath, resulting in the universe spending little to no time in the perturbative oscillatory phase. Thus, when the broad resonance regime pertains, we obtain a higher formal value for T RH given fixed inflaton mass and couplings than the perturbative estimate would indicate. A higher value of T RH in turn translates into more efficient thermalization between the two sectors. Of course, we expect that sectors preheated via a broad resonance may not have had time to attain internal thermal equilibrium by the end of reheating. This will also tend to enhance inflaton-mediated scattering, as the regions of phase space with momentum p ∼ M φ /2 will be preferentially occupied relative to the equilibrium distribution.
Given Γ φ and M φ , for any condition on q, we can thus identify three regions: (i) where the condition cannot be met for any Φ in the allowed range of eq. (4.10); (ii) where the condition can be met for some but not all allowed Φ; and (iii) where the condition is met for all allowed values of Φ. In figure 4 we have indicated how these regions for the broad resonance condition q > 1 intersect the parameter space that realizes thermalization. Above the dashed red line, broad resonance is operative for all Φ, while below this line, broad resonance is realized for some Φ.

Scalar inflaton coupling to Dirac fermion pairs in both sectors
We next consider the case where a scalar inflaton has renormalizable couplings to a pair of fermions in each sector. For simplicity, we take these fermions to be Dirac, but results for Majorana fermions (such as, e.g., right-handed neutrinos in the SM) are very similar. The relevant Lagrangian in this case is yielding the zero-temperature partial width The summed, averaged amplitude forψ a ψ a ↔ψ b ψ b scattering is (4.13) The full expression for the energy transfer rate for this case is derived in appendix B.2, and plotted in figure 11. In this case, Maxwell-Boltzmann distributions are not an unreasonable guide to the behavior of the thermal average, giving results that vary by only factors of order unity from the exact result. We use the approximate analytic expression to the Fermi-Dirac result of eq. (B.7) in our numerical scans. In the middle panel of figure 3 we plot the temperature dependence of the energy transfer rate for this model compared to the Hubble rate. In figure 5 we show the parameter space where thermalization can be achieved in this model. We cut the parameter space off when the width exceeds w > 0.01. Note that the maximum T RH obtained in these models is smaller than in the previous case with final state scalars: this highlights the impact of Pauli blocking at T RH > M φ . It is immediately obvious from the figure that in this case thermalization requires T RH > M φ . Here, we can estimate (dropping numerical factors) that, neglecting resonant effects, thermalization at T RH requires (4.14) which cannot be satisfied for T RH < M φ . For T RH > M φ , the availability of the resonance helps boost the energy transfer rate even above the estimate from dimensional analysis, and can enable the two sectors to thermalize. For fermions, Pauli blocking prevents non-perturbative particle production from being an efficient dissipation mechanism for the inflaton condensate, regardless of the strength of the coupling y or the amplitude of oscillation Φ, as we review in appendix C.4. Figure 5. Shaded regions indicate the parameter space where inflaton-mediated scattering can thermalize two otherwise decoupled sectors, in the case where a scalar inflaton couples to fermions in both sectors, for M φ = 10 3 GeV (left), 10 8 GeV (center), and 10 13 GeV (right). The upper boundary on these plots is determined by the maximum reheat temperature that can be obtained with w < 0.01.

Pseudo-scalar inflaton coupling to gauge boson pairs in both sectors
The third case we consider is a pseudo-scalar inflaton with axion-like couplings to gauge bosons in both sectors. We consider Abelian gauge fields for simplicity; the generalization to non-abelian cases is straightforward. The relevant Lagrangian is The zero-temperature decay width to each sector is then where we have included possible finite photon masses from e.g. thermal effects. The summed, averaged amplitude for γ a γ a ↔ γ b γ b scattering is then (4.17) In order for this EFT description of the axion interactions to make sense, we require M φ < Λ i . We will also require T RH < Λ i in our numerical scans over parameter space, in order to maintain the validity of the EFT in describing scattering at high temperatures. Both conditions place nontrivial restrictions on the zero-temperature mass-to-width ratio w.
We discuss the derivation of the resulting energy transfer rate in appendix B.3. In this case, the full use of Bose-Einstein statistics does not lead to dramatic enhancements in the energy transfer rate at high temperatures. This occurs because the non-renormalizability of the inflaton couplings Figure 6. suppresses contributions to the scattering amplitude from low momenta, with the consequence that Maxwell-Boltzmann is a fairly good approximation to the full result, as shown in figure 12 (see also the right panel of figure 3). The analytical approximation we use in our scans is given in eq. (B.11).
In figure 6 we show the region of parameter space where thermalization can occur in this model. We cut the parameter space off somewhat above the EFT bound to show context; for instance, in the case of M φ = 10 13 GeV, thermalization is not possible for T RH < Λ i . Thermalization in models with lower values of T RH is dominated by the resonant enhancement at T M φ , while in models with higher values of T RH , thermalization occurs for T ∼ T RH . This is evident from the right panel of figure  3. In the lower curve in this figure (corresponding to the lowest value of the width that equilibrates at k = 0.5), thermalization occurs at T < T RH as scattering becomes resonant following reheating. As the width, and thus reheating temperature, is increased, the energy transfer rate can exceed the Hubble rate at reheating itself. However, thermalization cannot occur unless T RH M φ , as again can be supported by dimensional analysis.
In this model, non-perturbative particle production can again become the dominant process governing the decay of the inflaton condensate for sufficiently large couplings to matter Λ and sufficiently large field oscillations Φ, as reviewed in appendix C.3. The red dashed line in figure 6 shows the boundary between the regime where broad resonance preheating is important for any physical Φ, and the regime where broad resonance preheating is important for some but not all physical Φ. All thermalizing regions within the range of validity of the EFT lie in the region where preheating can be, but is not necessarily, realized. As for scalars, we expect that when preheating is important, it will tend to enhance thermalization, thanks to both the higher values of T RH and the imperfect equilibration of the radiation sectors.

Hybrid cases
Finally, we consider cases where the particles that couple to the inflaton in the hidden and visible sectors follow different quantum statistics. We consider both the case of a scalar inflaton coupling to both fermions and scalars, and a pseudo-scalar inflaton coupling to fermions and gauge bosons.

Fermion-scalar scattering via scalar exchange
The relevant Lagrangian in this case is where again for simplicity we take the scalar, S, to be real and the fermion, ψ, to be Dirac. Considering S to be in sector a and ψ in sector b, this translates into our scan parameters as where k can now take on any real value (here we have dropped the small external masses). The summed, averaged amplitude for SS →ψψ scattering is The full expression for the energy transfer rate for this case is derived in appendix B.4, and plotted in figure 13. We use the approximate analytic expression of eq. B.14 in our numerical scans. As in our previous examples, dimensional analysis indicates that thermalization of the two sectors through inflaton exchange requires T RH > M φ . The region where thermalization occurs is shown in figure 8. As indicated by the red dashed line, the parameter space that realizes thermalization will also generically undergo broad resonance preheating, with the exception of a limited parameter space at lighter inflaton masses where preheating can but does not necessarily occur.

Fermion-gauge boson scattering via pseudo-scalar exchange
Here we consider a pseudo-scalar inflaton that couples to Dirac fermions in one sector, and (Abelian) gauge bosons in another sector. The relevant Lagrangian is then given by We have parameterized the fermion-pseudoscalar coupling in terms of a chiral-symmetry-breaking mass m ψ , as generally in realistic models this interaction is effective dimension five; we thus typically expect the effective Yukawa coupling m ψ /Λ b 1. 3 The inflaton's zero-temperature partial width to fermions is then In each case, plotted is the energy transfer rate at the lowest value of the width that thermalizes (lower curve), and the maximum allowed value of the width (blue, dashed). For the fermion-scalar case, we show the rate for k = 2, while for the fermion-gauge case we show the rate for k = 0.5. We also show the Hubble rate (red, solid). The black asterisk * indicates the reheat temperature corresponding to the lowest width that thermalizes, while the red # indicates the maximum value of the reheat temperature.
giving Λ a = 1 16 (we again here neglect subleading dependence on the external masses). In principle, k can take any value, but as we typically expect gauge couplings to dominate reheating when they are available, we will consider only k < 1 in our scans. We fix m ψ = 10 −2 M φ and allow Λ b to vary, subject to Λ b > T RH , M φ . This necessitates using an equal external mass for the gauge bosons, as our analytic approximation to the energy transfer rate holds in the limit that all external masses are equal; however, as m γ = m ψ M φ , and thermalization once again becomes most important in the regime T RH > M φ , the values taken for the photon external mass are not numerically important. The summed, averaged amplitude forψψ ↔ γγ scattering is (4.24) The full expression for the energy transfer rate for this case is derived in appendix B.5, and plotted in figure 14 (see also the right panel of figure 7). We use the approximate analytic expression of eq. Figure 8. Parameter region where thermalization can be achieved for the scalar-fermion model for M φ = 10 3 GeV (left), 10 8 GeV (center), and 10 13 GeV (right). The upper boundaries of these regions are set by the maximum attainable reheating temperature in each case. For fermion-dominated decays, the upper limit is set by perturbativity, which limits w < 0.01. For boson-dominated decays, conservation of energy in instantaneous reheating fixes the upper limit. The red dashed line indicates the lower boundary of the region where preheating is always important. Figure 9. Parameter region where thermalization can be achieved for the fermion-gauge model, for M φ = 10 3 GeV (left), 10 8 GeV (center), and 10 13 GeV (right). Above the black dashed line, the EFT description of the inflaton couplings breaks down, while above the red dashed line preheating is always important.
In the right panel of figure 7 we show the evolution of the energy transfer rate and the Hubble rate for this model. Figure 9 shows the regions of parameter space where the energy transfer rate exceeds the Hubble rate for a given inflaton mass. As in the gauge-gauge case, the portion of thermalizing parameter space where the EFT is reliable lies in the regime where broad resonance preheating can be but is not necessarily realized.

Discussion and Conclusions
Asymmetric reheating is one of the minimal cosmological origins for a dark sector that was never in thermal equilibrium with the SM, and as such is well-motivated by a broad and growing class of DM theories that contain multiple stable dark states. Such theories already frequently have trouble accommodating CMB measurements of the effective number of neutrinos, N eff , unless the dark sector is at lower temperatures relative to the SM than are possible if the two sectors were ever in thermal equilibrium at any stage. As future CMB experiments envision an order-of-magnitude improvement over Planck's current sensitivity to N eff , continued SM-like measurements of N eff will go from restrictive to prohibitive for once-equilibrated dark sectors containing dark radiation, as we have sketched in section 2. More generally, the assumption that a dark sector has thermalized with the SM is actually a fairly strong condition on the leading interaction(s) between the two sectors, and from this standpoint dark sectors that never thermalize with the SM represent a generic scenario for the origin of dark matter.
Here we have systematically studied minimal models of asymmetric reheating, demonstrating that the inevitable resulting inflaton (or more generally 'reheaton') exchange can itself thermalize the two sectors, wiping out any would-be temperature asymmetry. Thus realizing a sufficiently large temperature asymmetry between two sectors either imposes new restrictions on possible inflaton masses, coupling structures, and reheat temperatures, or requires a non-minimal cosmology.
We typically find that thermalization becomes important when T RH > M φ . The case where reheating is dominated by trilinear inflaton couplings to scalars in both sectors is an interesting exception, as in this case thermalization can occur in the infrared. In most cases, however, thermalization is driven by potentially resonant inflaton-mediated scattering at temperatures ∼ T RH . This region of parameter space, with T M φ , affords a rich variety of interesting and model-dependent phenomena, including nonperturbative particle production [38,39,41] and thermal effects such as thermal blocking and Landau damping [44,[55][56][57], and as such has been a focus of exploration in detailed studies of reheating. We have argued that we generically expect both incomplete internal thermal equilibration and nonperturbative particle creation in the form of preheating to make inflation-mediated energy transfer between sectors more, rather than less, important.
Our study generally serves as the first systematic study of multi-sector reheating, and as such offers many areas for future study. Chief amongst these is to extend our calculations to quantitatively establish the temperature asymmetry that follows from a given set of zero-temperature inflaton partial widths. A full exploration of two-sector thermalization when collective effects in one or both radiation baths can no longer be largely neglected also represents an obvious avenue of investigation for future investigation, as does the question of scattering in sectors that have not attained internal kinetic or thermal equilibrium.
of Ontario through the Ministry of Research and Innovation. YC is also supported in part by the Maryland Center for Fundamental Physics. We thank the Aspen Center for Physics for its hospitality and support through National Science Foundation Grant No. PHY-1066293. JS additionally thanks the Perimeter Institute for its hospitality during the completion of this work.

A Collision integrals
In this appendix, we provide a detailed derivation of our evaluation of the collision term in the Boltzmann equation. Our derivation draws on the derivation presented in the appendix of ref. [65], which in turn is based on the earlier work of ref. [64]. In this work, we have extended the previous results by giving general expressions for the thermally averaged cross-sections in terms of the Mandelstam variables. Further, we derive some (to our knowledge) new analytic results in the limit that all particles have the same mass.

A.1 The 2 → 2 collision operator
In this work, we are interested in processes in which two standard model particles scatter with two hidden sector particles via inflaton exchange. Ultimately, we are interested in the rate at which these scatterings transfer energy between sectors; however, this rate is closely related to the collision operator describing scattering.
We want to perform the integral 4 |M(12 → 34)| 2 is the matrix element of the process in question and S is a possible symmetrization factor to account for identical particles in the initial and/or final state. For notational simplicity, in our derivations here we will absorb the g i into an effective symmetrization factorŜ. The plus sign here is for bosons, while the minus sign is for fermions, accounting for the effects of Bose enhancement and Pauli blocking, respectively. We assume local thermal equilibrium and neglect chemical potentials, so that the distribution functions take their equilibrium forms where Υ = ∓1 for bosons and fermions, respectively, and Υ = 0 for classical particles. In the limit of equilibrium, the full collision integral vanishes. We are interested in the rate at which collisions occur relative to the Hubble rate. We thus consider only the forward direction, which is equivalent to the product of the thermally averaged cross-section and the number densities of species 1 and 2, We now need to perform the integrations.
We begin by writing the integral of the Lorentz-invariant phase-space of particle 4 as The integral over d 4 p 4 is then done using the overall energy-conserving delta function δ 4 (p 1 + p 2 − p 3 − p 4 ). This then fixes where we have introduced s, t and u, the standard Mandelstam variables. The integration measures for the 3-momenta, p 2 and p 3 , are written d 3 p 2 = p 2 2 dp 2 d cos α dβ, d 3 p 3 = p 2 3 dp 3 d cos θdµ, where the angles are defined as We work with s, t and u instead of the angles α, θ and β, in terms of which The integration measure for t and s is straightforward to obtain. To integrate over u instead of β, we need to solve for sin β. To this end we note 12) and observe that the integration is restricted to regions where the argument of the square root is positive (this is equivalent to cos 2 β ≤ 1). The integration measure for p 2 and p 3 is written d 3 p 2 d 3 p 3 = 1 8 dp 2 dp 3 p 2 1 dµ ds dt du g(s, t, u) Θ(g(s, t, u)). (A.13) Explicitly, where from eq. (A.10), cos α and cos θ are considered functions of s and t, respectively. There is no dependence on µ, and the integral over this variable can be trivially performed to give a factor of 2π. We can combine this with the rest of the expression, finding For a completely general matrix element, this is as far as we can go in full generality. We can make some more progress if we assume, as appropriate for the interactions important in this work, that the matrix element is a function only of either s, t, or u.
For each variable, note that g(s, t, u) is a quadratic function of s, t, and u. After the integration over the delta function is performed, it is a quadratic function of only two of s, t or u, and we can make use of the fundamental identity [64,65] dx where x is any of s, t or u, and a, b and c are (quadratic) functions of the remaining Mandelstam variable. In this work we are interested primarily in energy transfer between sectors. As discussed above in section 3.3, energy transfer is dominated by s-channel processes in the regimes of interest. We thus focus on these. Similar results for u and t channel amplitudes can be readily obtained following the same steps.

A.2 s-channel scattering
Assuming that the matrix element depends only on s, we integrate over u using the delta function, and make use of eq. (A.16) to integrate over t. In this case, the parameter a on the right hand side of eq. (A.16) is The Heaviside function, Θ(b 2 − 4ac), restricts the integration to regions where In this expressions = s − (E 1 + E 2 ) 2 and After switching integration variables to E 1 , E 2 and E 3 , eq. (A.15) can be written To make further progress, we need to find the region of integration where eq. (A.18) holds. To this end, recall that by definition (see eq. (A.10)) so in order to have a non-zero integration region, A > 0, we need to find the regions where both hold.
There are then four cases which can lead to non-zero integration regions with A > 0 The limits in eq. (A.18) can be used to restrict the integration over E 3 . In general, the radicals that appear when the three momentum, p, is replaced by the corresponding energy E in these conditions make them difficult, if not intractible to reduce analytically. However, if we assume all particles have identical masses, some progress can be made. Assuming all four particles have identical masses, m, we find that only the first two regions 1. & 2. above contribute and lead to the following restrictions.
• From the first case above, the integration limits are • For the second case, we find for E 2 < E 1 and for E 1 < E 2 and If we further assume that particles 1 and 2 obey the same statistics, as do particles 3 and 4, we find where Next, we can integrate over the distributions. The integrals over E 3 and E 2 are easily done, yielding where Υ a , Υ b = 1 for Fermi-Dirac statistics, Υ a , Υ b = −1 for Bose-Einstein and Υ a , Υ b = 0 for Maxwell-Boltzmann. In the case where particles 3 and 4 obey different statistics to particles 1 and 2, the remaining three integrations (after inserting Eqns. (A.28) and (A.29) into eq. (A.27)) need to be performed numerically. If all particles obey the same statistics, Υ a = Υ b , we can analytically integrate over E − andẼ − as well, in which case we arrive at In general, these final two integrations must be performed numerically. Similarly, the energy transfer rate is given by the integral

A.2.1 Maxwell-Boltzmann limit
In the limit that all particles obey Maxwell-Boltzmann statistics, Υ a = 0, Υ b = 0, we can evaluate the integral over E + analytically, and recover the known results of Gelmini and Gondolo [66]. In the Maxwell-Boltzmann limit, eq. (A.30) becomes where here and below the K n (x) are modified Bessel functions of the second kind of order n. In the center of mass frame, so that if M depends only on s, Making use of the Maxwell-Boltzmann limit result we find the thermally averaged cross-section in the Maxwell-Boltzmann limit σv = 1 in agreement with the result of Gelmini and Gondolo [66]. We write this as A similar integral form can be found for the energy transfer rate

B Energy transfer rates
In this section we consider specific simple models as discussed in section 4, and develop some approximations that allow for analytic integration of the thermally averaged cross-section, as well as the energy transfer rate, in the Maxwell-Boltzmann limit. Away from the resonance, only small corrections are required to obtain good agreement with the full numerical evaluation including full quantum statistics.
For numerical efficiency, we make use of the approximations derived below in our scans over parameter space in section 4. In all cases, we do not expect that the poorer agreement of our approximation with the exact result near T ∼ M φ will significantly affect the results of our scans. This can be seen from figure 3 and figure 7. Notice that in each of the limiting cases illustrated in these figures, H(T ) intersects Γ E (T ) at the temperature where Γ E (T ) is at a local maximum, which occurs in the region T < M φ 5 where our approximations are excellent. Further, note that whenever our approximations mis-estimate Γ E , the combination of the shape of resonant feature in Γ E and our scan strategy that looks at all T < T RH means that we will not misattribute thermalized regions as unthermalized or vise-versa. However, if one were only interested in the condition Γ E (T RH ) > H(T RH ) at a single fixed T RH , then these approximations would be inadequate.
All numerical result in this section were obtained using the publicly available Cubature, 6 and Cuba 7 [67] libraries for C.

B.1 Scalar boson scattering
We first consider the scattering of real scalars, S a , S b , via a trilinear interaction with a third scalar φ, as described in eqs. (4.3)-(4.6). To a very good approximation, the amplitude for this process can be approximated by treating the resonance as infinitely sharp and taking the matrix element to be simply a constant below the resonance and vanishing above, Somewhat surprisingly, comparison of a direct numerical evaluation of the thermally averaged crosssection in the Maxwell-Boltzmann limit shows that this is an excellent approximation even when Γ φ M φ . 5 This behavior is straightforward to understand. Near the resonance, Γ E takes the form x −3 K 2 (1/x), with x = T /M φ . The function K 2 (x) can be approximated near x ∼ 0 as x 1/2 exp(1/x). The function x −5/2 exp(1/x) peaks at x peak ∼ 2/5, and thus the local maximum due to the resonance occurs at T = 2M φ /5, somewhat below T ∼ M φ . 6 http://ab-initio.mit.edu/wiki/index.php/Cubature 7 http://www.feynarts.de/cuba/ In the case where all external particles obey Maxwell-Boltzmann statistics, it is straightforward to obtain for the thermally averaged scattering cross-section from eq. (A.38) where the functions K i (x) are the modified Bessel functions and m = m S is the mass of the scalar boson. Similarly, we can obtain the thermally averaged energy-transfer rate using eq. (A.39) Here G 3,0 1,3 (x) is the Meijer G-function. It is the subleading part of the expression at high T /m and can be dropped in this regime, while it changes the result by order unity at low T /m ∼ 1.
For Bose-Einstein statistics, the energy transfer rate at eq. (A.31) is well-approximated by where A = 1.06 and B = 0.615. As can be seen in figure 10, this approximation is excellent away from the resonance region. It fails completely for T ∼ m. Unfortunately, we have been unable to derive this result in a consistent fashion by approximating the integrals of eq. (A.30).

B.2 Fermion scattering
We next consider the scattering of Dirac fermions ψ a , ψ b , via Yukawa interactions with a scalar φ, as described in eqs. (4.11)-(4.13). Analogous to the scalar trilinear interaction considered above in section B.1, to a good approximation, we can approximate this matrix element using an infinitely sharp resonance with power law (in s) behaviour away from the resonance region where m = m ψ is the fermion mass.  If we neglect quantum statistics and treat the fermions as classical particles, we can evaluate the energy transfer rate The cross-section can also be found using the same technique. In obtaining this result, as described above, we have expanded in the large T /m and T /M φ limits. Away from the Maxwell-Boltzmann limit, analytic solution of the energy transfer rate is difficult. We can, however, numerically evaluate eq. (A.31), finding that away from the resonance region it is well-approximated by where A ≈ 0.85 and B = 2 √ 2 5 . The results of these approximations are shown in figure 11.  We now consider the scattering of gauge bosons via a dimension-five interaction with a pseudoscalar φ, as described in eqs. (4.15)- (4.17). Analogous to the above cases, this matrix element is well-approximated by

B.3 Gauge boson scattering
where m = m γ is the gauge boson mass.
In this case, if we ignore quantum statistics, all integrations are analytically tractable and the result for the cross-section can be obtained in terms of Meijer G-functions. However, the Meijer Gfunctions are not easily implemented numerically, and do not provide much insight. For these reasons we expand our results as power series in T /M φ and T /m. For m < T and and M φ < T , to an excellent approximation we can expand and keep only the leading order behavior in powers of T where here and below the Heaviside functions are put in by hand, and chosen to cut off the power-laws in the appropriate places. For T ∼ m, this result breaks down. However, further expansions in the low T /m regime are possible to model the result in this region. Similarly, we can calculate the energy transfer rate. With the above approximation for the amplitude, this result can also be found exactly in terms of Meijer G-functions; however, the result is not particularly illuminating and we omit it here. The Meijer G-functions can be expanded in the limit of large T /m and T /M φ to give to a good approximation This is an excellent approximation for T > m. An expansion is possible in the low T /m limit which can capture this limit extremely accurately. We do not present it here. Away from the Maxwell-Boltzmann limit, with Bose-Einstein statistics, eq. (A.31) is difficult to evaluate. However, numerically we find that away from the resonance, the full result using Bose-Einstein statistics is almost indistinguishable from Maxwell-Boltzmann statistics. eq. (B.9) and (B.10) are reasonable approximations for the full energy transfer rates found by evaluating eq. (A.31). These approximations can be improved by correcting the asymptotic power law behavior. By comparing comparing the result of eq. (B.10) to the numerically evaluated eq. (A.31), we find that the T > M φ asymptote should be rescaled by a factor of approximately 4/π, giving In figure 12 we show the results of these approximations.  Our model for fermion-scalar scattering is described in eqs. 4.18-4.20. The amplitude is wellapproximated by where m = m S = m ψ is the mass of the scalar boson and the fermions, which we take to be equal.  Away from the Maxwell-Boltzmann limit, analytic evaluation of the full energy transfer integral with full quantum statistics is difficult. However, as described above, we can compare the full result found by numerically evaluating eq. (A.31) with our Maxwell-Boltzmann approximation in eq. (B.16) and correct the asymptotics to find The results of these approximations are shown in figure 14. We note that away from the T = M φ , this is an excellent approximation.

C Preheating
In this appendix we briefly review the physics of preheating and sketch out the regions of parameter space in each of our simple reheating models where we expect that preheating may be important. The theory of preheating is a well-developed field; see, for example, refs. [30,68,69] for recent reviews.

C.1 Preheating in brief
The details of the reheating phase that immediately follows inflation are model-dependent. In particular, the overall description of particle production is dependent on both details of the inflationary model, such as the scale at the end of inflation V (φ i ) and the structure of the potential probed by the field oscillations, as well as on the details of the particle model describing the couplings of the inflaton to other degrees of freedom. In this section we briefly present the general conditions for the universe to undergo a period of strong, non-perturbative particle production, dubbed preheating, sourced by the time-dependent inflaton background immediately following inflation, and discuss how to understand these general conditions in the parameter space presented in section 3.2.
Preheating can lead to rapid particle production through the parametric resonance of momentum modes in the daughter fields which couple to the inflaton. The oscillating inflaton background causes non-adiabatic changes in the effective frequency of these daughter fields which can result in their copious production. This non-perturbative particle production can be, but is not always, highly efficient at transferring energy from the inflaton condensate into radiation, thereby resulting in a much higher reheating temperature than that obtained from estimates based on perturbative decay. By contrast, when preheating is incomplete (in the sense that it cannot dissipate an O(1) fraction of the initial energy density from the condensate), it will be followed by a phase of perturbative reheating, and the ultimate impact of preheating on T RH will typically be small.
In our study, we have assumed that at the beginning of the reheating phase, the energy density of the universe is overwhelmingly dominated by a single scalar field with a quadratic potential. When the Hubble rate of the universe drops below the mass of the scalar field, M φ , it begins to oscillate approximately sinusoidally with a decaying envelope. While the expansion of the universe is extremely important for determining whether preheating is an effective means of transferring energy out of the inflaton condensate, the essential physics of the preheating process can be understood in flat space. In this limit, daughter fields χ coupled to the oscillating inflaton can typically have their equations of motion recast as the Mathieu equation, Here the derivative acts over z = M φ t + δ, where δ is a phase, and A k and q depend on the parameters of the theory, such as the initial background field amplitudes and the couplings and masses of the external field(s) χ. The Mathieu equation can be solved exactly. The standard arguments used in the analysis of this equation are reminiscent of those used to solve the Schrodinger wave equation in a periodic potential, with solutions given by Bloch waves, where P is periodic in z with period π. The quantity m k is called the Mathieu exponent, and its properties determine the stability of the solutions. For growing solutions, the real part of the Mathieu exponent [m k ] = µ k is always non-negative. If µ k (q) = 0 then |χ k | is stable while if µ k (q) > 0, then the field amplitude |χ k | grows exponentially.
Contour plots of µ k are usually presented as a function of the parameters A and q, which together control the width of the unstable bands and their strength µ k . In general it is useful to distinguish two regimes: • The narrow resonance regime, where q < 1. In this regime, µ k 1, and the growth of particle occupation numbers can be moderate compared to the expansion of spacetime.
• The broad resonance regime, where q > 1. In this regime µ k O(1), and particle production is explosive.
The occupation numbers in the unstable bands grow exponentially as n k ∼ exp(2µ k z) = exp(2µ k M φ t).

(C.3)
In an expanding spacetime, the redshifting of momentum modes will modify this picture, most importantly by shutting off the parametric resonance. Two necessary conditions for preheating to be successful, i.e., for an O(1) fraction of the inflaton energy density to be dissipated through preheating, are: 1. The rate at which modes are amplified must proceed faster than the decay rate of the inflaton, and faster than the rate at which the inflaton oscillations are damped The rate at which preheating proceeds must be faster than the rate at which the modes are redshifted out of the resonance band. The time that the mode remains in the resonance band depends on the equation of state of the matter; however, in the narrow resonance regime, it can be estimated as [39] ∆t ∼ qH −1 .

(C.5)
Requiring that the parametric resonance rate exceeds this, one finds the condition Typically in the narrow resonance region (q < 1) eq. (C.6) is a stronger condition than is eq. (C.4). Both conditions are satisfied in the broad resonance region q > 1. Unless preheating begins with q 1, and ends near q ∼ 1, it is unlikely to be completely effective at reheating the universe, and it is expected that the final reheating temperature will be determined by perturbative decays [69]. However, preheating will alter the details of the distribution functions of fields that are coupled to the inflaton and can affect the timescale for the radiation bath to attain internal thermal equilibrium.
The parameter regions relevant for preheating depend on the specific model in question, so we will consider each case in turn.
(C. 8) In this case, modes with k 2 + m 2 S < µΦ are exponentially enhanced, where here and throughout Φ is the initial inflaton amplitude. For this model, we can identify the Mathieu parameters Since for A k < 2q the frequency of modes with momentum k can become negative during part of the oscillation period, there is a tachyonic instability in this model. The regime in which this instability is important is known as 'tachyonic resonance.' When q 1, many dynamical modes are tachyonic for part of the oscillation, and the resulting period of preheating is extremely efficient and reheating can occur within a few oscillations of the inflaton. However, provided q < 1, the width of the stability and instability bands are comparable, and tachyonic resonance becomes indistinguishable from narrow parametric resonance 8 [70].
As the criteria for efficient preheating depend on both µ and Φ for a fixed value of M φ , specifying M φ and the perturbative value of T RH is not sufficient to determine whether or not preheating is important in a given model: the cosmological history, i.e., the field amplitude Φ at the onset of reheating, must also be specified. As discussed in section 4.1, only a bounded range for Φ is consistent with the minimal model of reheating that we adopt here, which allows us to broadly distinguish regions in the parameter space (M φ , T RH ) where preheating is, may be, or is not relevant.

C.3 Axion couplings to gauge fields
In this work we have considered the coupling of an axion to gauge fields via the dimension-5 operator where Λ is a mass scale associated with the axion. This case is unique, because it involves a derivative coupling to an irrelevant operator.

Preheating constraints
This case is treated in refs. [71,72]. The equation of motion for the gauge modes, B ± k , in flat space is given byB where ± denotes the two possible helicity states of the gauge field, and the derivatives act over time t. Defining z = M φ t, and writingφ = ∓M φ Φ cos(2z), we can again identify the Mathieu parameters: The different signs for the two different helicities have been absorbed into the choice of phase. Note that it may appear that one can get large q here, and thus efficient preheating, by looking at larger wavenumbers. However, as one increases k, A also increases and pushes one into the stable regions of the Mathieu chart. Parametric resonance in these models is most effective for the lowest momentum modes possible, and thus we will use the value of q for the lowest dynamical mode, k ∼ aH ∼ M φ Φ/M pl . For this mode, we have A k < 2q, which is similar to the tachyonic resonance case considered for scalars above. The narrow resonance condition becomes while the condition that narrow resonance is inefficient is Pl . (C.14) For large field inflation, Φ ∼ M Pl , notice that these conditions roughly coincide.

Inflationary constraints
While in the present work we have been agnostic about whether the field that drives reheating is the same field that is responsible for inflation, it is worth pointing out that in the case where the axion does drive inflation, the coupling to gauge fields can lead to the exponential production of gauge bosons during the inflationary phase itself [73,74]. These gauge bosons can scatter off the inflaton condensate and generate fluctuations in the axion field. The axion fluctuations generated this way are strongly non-Gaussian, and can lead to unacceptable levels of non-Gaussianity in the CMB [75,76] and primordial black holes at the end of inflation [77] if the couplings are too large. The conditions on the inflationary production of gauge quanta are quoted in terms of a parameter ξ, defined as The bounds from non-Gaussianity in the CMB [78] are ξ < 3.3 during the observable e-folds of evolution. We can use the constraint on the tensor-to-scalar ratio r < 0.07 − 0.09 [79], and the slow-roll condition r = 16 H to place a bound on the scale Λ, Again, this bound is only applicable if the axion in question is the field that drives the period of slow-roll inflation constrained by the CMB.

C.4 Fermion Preheating
The preheating of fermions via a Yukawa coupling of the form L = y i φψ i ψ i (C. 17) was first studied in detail by Kofman and Greene [80]. The efficiency of fermion preheating is limited by the fact that fermions obey the Pauli exclusion principle, which prevents them piling up in one state. However, fermion preheating can be somewhat effective at populating fermion states. For the Yukawa theory, the analogue of the q parameter is given by (C. 18) Preheating fills up fermion states up to a maximum wave number of [80] k max ∼ a 1/4 q 1/4 M φ , (C. 19) which grows like a 1/4 . While in the bosonic case, taking q > 1 leads to efficient 'broad resonance' preheating, for fermions the effect of a large q parameter is to excite modes up to a higher wavenumber. In order to see if fermion preheating can significantly alter our conclusions, we examine the amount of energy that is transferred non-perturbatively from the inflaton condensate to fermions. The energy density in the preheated fermions can be estimated as Taking E k ∼ k/a, and n k = 1/2 up to k max , we have (C. 21) In the absence of their perturbative decay, the preheated fermions remain a fixed (small) fraction of the total energy density. For m ψ M φ , the zero-temperature inflaton width into fermions is given by so we can express the above as Thus for perturbative couplings the impact of preheated fermions on our calculations is negligible.