Symmetric and Asymmetric Reheating

We study models in which the inflaton is coupled to two otherwise decoupled sectors, and the effect of preheating and related processes on their energy densities during the evolution of the universe. Over most of parameter space, preheating is not disrupted by the presence of extra sectors, and even comparatively weakly coupled sectors can get an order 1 fraction of the total energy at this time. If two sectors are both preheated, the high number densities could also lead to inflaton mediated thermalisation. If only one sector is preheated, Bose enhancement of the late time inflaton decays may cause significant deviations from the perturbative prediction for their relative reheat temperatures. Meanwhile, in Non-Oscillatory inflation models resonant effects can result in exponentially large final temperature differences between sectors that have similar couplings to the inflaton. Asymmetric reheating is potentially relevant for a range of beyond the Standard Model physics scenarios. We show that in dark matter freeze-in models, hidden sector temperatures a factor of 10 below that of the visible sector are typically needed for the relic abundance to be set solely by freeze-in dynamics.


Introduction
Dark sectors decoupled from the Standard Model are well motivated both from UV and IR considerations. String compactifications that contain the Standard Model (SM) often also include dark sectors. Meanwhile, dark sectors can include viable dark matter (DM) candidates, with the possibility of self interactions [1,2], or distinct indirect detection signals [3][4][5]. However if a hidden sector contains light states, observable cosmology could be disrupted. For example, energy in light dark sector degrees of freedom increases the effective number of neutrino species N eff , and the current bound on contributions beyond the SM is ∆N eff 0.25 [6,7]. New sub-keV states are further constrained by measurements of the Lyman-α spectrum, which gives the leading lower mass limit on thermal DM, m DM 3 keV [8]. A straightforward way to relax tension with these observations is through the dark sector having a temperature T DM substantially lower than that of the visible sector T SM [9,10].
Further, some DM models, as well as other motivated beyond the SM scenarios, need significant temperature differences between the visible sector and a dark sector. In freezein DM the dark and visible sectors are almost decoupled [11], and the DM relic density is produced by a small interaction between the dark sector and the SM thermal bath. As we show, for the relic abundance to be determined only by the form of this interaction a dark sector that is initially significantly colder than the visible sector is required, even if the DM has annihilations to light hidden sector states. In particular, we impose that any annihilations are frozen out before the freeze-in starts. Models in which the DM is a glueball of a strongly coupled hidden sector gauge theory also need the visible sector to be hot compared to the other sectors [1]. Stable glueballs from such sectors are expected in some compactifications of string theory [12][13][14][15]. Additionally, differences between the energy in separate sectors at early times could have important consequences, independent of their final relative temperatures.
For example, a sector having a high temperature in the early universe can lead to processes such as phase transitions or production of heavy or weakly coupled states [16][17][18][19].
The simplest reheating mechanism is perturbative inflaton decays, the speed of which is fixed by the single inflaton particle decay rate [20][21][22]. If the inflaton dominantly decays to a particular sector, this will be reheated to a higher temperature than others [23], with a resulting temperature ratio that has the approximate parametric dependence T 1 /T 2 ∼ g 1 /g 2 , where g i is the inflaton coupling to kinematically accessible states in the two sectors. 1 However, in models with relatively large couplings between the inflaton and other states, nonperturbative processes can lead to resonantly enhanced energy transfer during a period prior to perturbative decays, called preheating. The dynamics of this has been studied extensively in single sector models [18,[25][26][27][28].
In this paper, we study the impact of preheating, and other associated processes, on the relative energy in two sectors that are decoupled, apart from both being coupled to the inflaton. The effects are complex, involving non-perturbative dynamics and particle distributions that are far from thermal equilibrium. To make progress we have to resort to a range of approximations and for preheating itself numerical simulations are useful. Consequently we do not obtain precise results, aiming instead to identify interesting physical processes that can be important. We also do not attempt to study complete realistic models of the visible and hidden sectors. Typically we just take two copies of a simple toy model sector, identical other than their coupling to the inflaton, and try to find their relative energy densities as a function of their couplings to the inflaton as the universe evolves. Already this leads to several possible effects, although further study of more complex models would be interesting.
One issue is whether the dynamics of preheating to a sector are modified by the presence of other sectors coupled to the inflaton also being preheated. This can happen since backreaction or rescattering from a second sector might cause preheating to end sooner than it otherwise would. Immediately after preheating, thermalisation between sectors by off-shell inflaton scattering could also change the distribution of energy, potentially enhanced by the high number and energy densities produced by preheating. In many models there will subsequently be a long period of matter domination, due to the energy remaining in inflaton states. In this case, the late time relative temperature of two sectors will be close to the perturbative prediction unless inflaton mediated thermalisation is efficient, although as noted the early time dynamics might still lead to physically important effects. 2 However, if the inflaton has large trilinear couplings it might decay perturbatively before matter domination is reached, and the dynamics of this could be affected by preheating. In such models, several complicated effects might occur, and final temperatures either closer together or further apart than perturbative inflaton decays alone are possible. As an example we study the effect of Bose enhancement of the inflaton decay rate.
We also note that exponentially large differences in the final temperatures of sectors with order 1 differences in their couplings to the inflaton are possible in Non-Oscillatory models of inflation [30,31]. In such theories, there are no perturbative inflaton decays. Meanwhile the energy transferred non-perturbatively is exponentially suppressed below a threshold coupling, and large temperature differences can occur without tuning parameters to particular values.
The structure of the paper is as follows: In Section 2 we motivate studying temperature difference between sectors by calculating the relative initial temperature ratio needed for the DM relic abundance to be set solely by freeze-in. In Section 3 we analytically and numerically study preheating in models with the inflaton coupled to two sectors. In Section 4 we consider effects such as thermalisation that can take place once preheating is finished. In Section 5 we describe how Non-Oscillatory inflation models with instant preheating can lead to large temperature ratios between sectors. Finally we discuss phenomenological implications of our work in Section 6.

Initial conditions for freeze-in dark matter
In models of freeze-in dark matter [11,[32][33][34][35][36][37][38], prior to freeze-in the abundance of DM is assumed negligible and the DM is generated by a portal operator connecting the hidden and visible sectors. The final relic abundance is controlled only by the DM mass and the portal operator if this is renormalisable (otherwise production is UV sensitive and may depend on the SM reheat temperature).
Unless the inflaton has no decays to the dark sector, reheating leads to a population of DM, which we call the primordial component. For freeze-in to set the DM relic density, the primordial component of DM (PDM) should be substantially less than the DM abundance observed today Ω PDM Ω B today Ω DM Ω B observed 5 . (2.1) The conditions after reheating needed to satisfy this depend on the details of the hidden sector. We consider two example scenarios, depending on if the dark sector has number changing interactions.

Dark sectors with no number changing interactions
One possibility is that reheating occurs via perturbative decays and the DM has no number changing interactions other than those involving the inflaton or freeze-in. If each inflaton decay produces order one DM state, then after reheating where Br φ→DM = Γ φ→DM /Γ is the inflaton (φ) branching fraction to DM, m φ is the inflaton mass, and ρ φ is the inflaton energy density at H Γ, with H the Hubble parameter.
Assuming that the visible sector dominates the entropy of the universe (and that there are no further entropy injections, for example through a late decaying scalar) the initial condition for freeze-in DM, Eq. (2.1), is then is the visible sector reheat temperature. For simplicity, we have assumed that the visible sector consists of just the SM at all energy scales up to T (RH) SM . Since the baryon asymmetry is small ((n B − nB)/s 0.88 × 10 −10 at the present time) a large temperature asymmetry is needed unless the reheat temperature is very low. Immediately after being produced by perturbative inflaton decays, the DM states will typically have extremely high kinetic energy. If the dark sector has non-number changing interactions, this could be redistributed amongst the DM including the freeze-in component. The warming generated by this is negligible, provided the branching ratio is 1 as is required by the relic density constraint almost all of parameter space. Otherwise if there are no interactions at all, the primordial DM will either red-shift to be non-relativistic or form a small hot DM component, depending on the parameters of the model.

Dark sectors with annihilations
Alternatively the DM might have number changing interactions, and we focus on the case in which it thermalises with a bath of light dark sector states. The contribution of the relic energy density of the new light states to the effective number of neutrinos can be sufficiently small, provided the hidden sector is slightly colder than the visible sector (which is needed for freeze-in anyway) [39,40]. Alternatively, it might be possible to construct models such that they decay to the SM. If the dark sector states are relatively heavy 100 MeV, constraints from collider observations easily allow decays before BBN [41]. Meanwhile, for lighter masses collider and astrophysics bounds require much longer decay times, and there is greater danger of observable energy injection [42], although again this may be allowed if the hidden sector temperature is relatively low. However, further study is required to ensure the coupling between the sectors does not disrupt freeze-in, and we leave this for future work.
While the dark sector is at sufficiently high temperatures, number changing interactions keep the DM number density close to its equilibrium value. However, as the hidden sector temperature drops, the DM abundance will freeze-out, similarly to the usual WIMP scenario. In order that the DM relic abundance is set purely by the freeze-in operator, there are two constraints that must be satisfied. First, the component of DM left from thermal freeze-out must be substantially smaller than the observed quantity. Second, DM annihilations must be negligible by the time freeze-in of the DM relic density happens (and remain small at later times). 3 Assuming no entropy injections or large changes in the effective number of relativistic degrees of freedom before freeze-in, the ratio of hidden sector and visible sector 3 A slightly weaker possible constraint is that annihilations of the thermal population must finish before temperatures is a constant ξ = T (RH) DM /T (RH) SM , determined by dynamics at an early time such as reheating. The fraction of the visible sector's energy transferred to the hidden sector per Hubble time through the freeze-in operator is largest when freeze-in is taking place, and at earlier times interactions through this coupling do not increase the relative temperature of the hidden sector significantly.
The calculation of the hidden sector freeze-out from the Boltzmann equations is straightforward, and closely follows [43]. The underlying physics can be understood from an approximate analytic solution accurate to the % level. Freeze-out finishes when the hidden sector temperature is m DM The parameter α is where M Pl is the (unreduced) Planck mass, and g v is the visible sector effective number of relativistic degrees of freedom, which are assumed to dominate the energy of the universe, and g is the effective DM number of degrees of freedom (for bosons g = g DM , and for fermions g = 7 8 g DM ). The annihilation cross section has been parameterised as σv ≡ (σv) 0 (m DM /T DM ) −n . The primordial relic density gives a yield y = n/s after freezeout of .
(2.6) From Eq. (2.4), the freeze-out temperature is related to the DM mass by a parameter that is typically less than 1 and is a slowly varying function of the properties of the model, similarly to the usual WIMP case for which T FO m DM /25. The relic density and freeze-out temperature constraints together lead to a maximum ratio of hidden sector and visible sector temperatures ξ, apart from for very large DM masses, for which the unitarity bound on annihilations is stronger [44]. The largest ξ are allowed if freeze-in takes place as late as possible, corresponding to models in which freeze-in happens at a temperature T SM m DM . In this case, the temperature constraint and Eq. (2.4) requires that ξ < T (FO) DM /m DM (which is an implicit equation since the freeze-out temperature depends on ξ). The relic density constraint is imposed by demanding that Ω PDM < 0.1 Ω DM in Eq. (2.6). In Fig. 1 (left) we illustrate the interaction of these constraints, as a function of the annihilation cross section, assuming n = 0 and a DM mass of 10 GeV.
The maximum allowed ξ as a function of the DM mass is plotted in Fig. 1 (right), for swave annihilation. At small masses this is approximately constant because of the logarithmic dependence of Eq. (2.4) on the parameters of the model, and for DM masses below a PeV, the hidden sector must be colder by a factor ∼ 10 relative to the visible sector. For masses freeze-in, allowing annihilations after the freeze-in component is generated. Such models are interesting [32,33], and the constraints obtained on the hidden sector initial temperature are close to those we obtain. 10 -10 10 -9 10 -8 10 -7 10 -3 In the blue shaded region, the primordial relic density exceeds 0.1Ω DM , while in the orange shaded region the DM annihilations freeze-out after freeze-in has occurred. In both cases the DM relic abundance is not set purely by the freeze-in mechanism. Right: The maximum allowed value of T RH DM /T RH SM as a function of the DM mass, if the relic abundance is to be set by freeze-in alone, assuming s-wave annihilation. At high DM masses 100 TeV the unitarity bound on the annihilation rate results in strong constraints on the hidden sector temperature.
above ∼ PeV the unitarity bound, σv r < 4π/ m 2 DM v r for s-wave annihilation and relative velocity v r [44], constrains the maximum possible cross-section, and much smaller hidden sector temperatures are needed. The unitarity bound becomes important at slightly higher masses than the usual WIMP case because, for a given DM mass, the constraint on the annihilation freeze-out requires a smaller annihilation cross section than that which would lead to the correct DM relic abundance for ξ = 1. Light DM masses 100 MeV, and hidden sector states for them to annihilate into, may require colder hidden sector temperatures to avoid constraints in complete models.
In other models, freeze-in will dominantly occur at temperatures T (FI) SM > m DM , leading to stronger constraints on ξ. For example freeze-in through a portal generated by a heavy state, with mass above the DM mass, that leads to a non-renormalisable operator linking the two sectors when integrated out. In this case, the bound on the freeze-out temperature is to order of magnitude given by since the hidden sector freeze-out temperature is typically T DM ∼ m DM /20. Precise bounds can be straightforwardly calculated. Meanwhile, if freeze-in is via non-renormalisable operators that are not UV completed below the reheating temperatures, the dominant freeze-in production occurs immediately after reheating. Consequently, for such models number changing interactions must be absent at all relevant temperatures and the constraints are as for Section 2.1.
If a hidden sector thermalises through 2 → n and n → 2 self interactions, rather than with a bath of light states, the dynamics are changed. In such models the DM number density decreases over time, resulting in the hidden sector heating up. For the motivated case of the glueball of a hidden sector being the DM, the relic density constraint is where Λ is the hidden sector confinement scale [13]. To avoid observational constraints on DM self-interactions Λ 100 MeV is needed [45], requiring a large initial temperature asymmetry.

Preheating in models with two sectors
In many models with relatively large couplings between the inflaton and other states, the first stage of reheating happens through a period of parametric resonance [27]. 4 In this section we apply existing results from single sector models to preheating in theories with multiple inflaton decay sectors. An especially interesting aspect is whether there are regions of parameter space where efficient preheating would be expected in the case of only one sector, but there is no efficient energy transfer because of the dynamics from the other sector. We focus on models in the broad resonance regime (defined shortly), since in an expanding universe the alternative of narrow resonance does not transfer a significant fraction of the inflaton energy. As is well known, preheating terminates early when an order 1 fraction of the total energy in still in inflaton modes, and this remains true in models with multiple sectors. The system then evolves in a complicated way, with the expectation that each sector will move towards a thermal distribution. There is also the possibility of extra energy transfer from the inflaton during this time, and narrow resonances could be important for the dynamics. However, due to the high number density of decay products, which have energy distributions far from equilibrium, at this stage the system must be studied using lattice simulations. Ultimately for reheating to complete perturbative inflaton decays, and therefore trilinear inflaton couplings, are needed [22,46]. These decays will have an important effect on the final relative temperatures of the two sectors and we study this, and other late time processes, in Section 4.
In broad resonance, the coherent oscillations of the inflaton act as a changing background mass m χ for the inflaton's decay products χ (the non-expectation value induced mass of χ must be smaller for efficient resonance). Over most of an oscillation, the evolution of the decay products' mass is slow compared to 1/m χ , and close to adiabatic, so that particle production is weak. However, when the inflaton is near the minimum of its potential, m χ is small and evolving fast. Consequently the system can be highly non-adiabatic, and the number density of the decay products may grow extremely fast. Assuming the decay products are bosons also results in a Bose enhancement and exponential growth in number density. In particular, during preheating, one oscillation of the inflaton leads to particle production where n i k is the number density of the momentum mode k of χ after i inflaton oscillations, and µ i k is a numerical coefficient. µ i k itself is a fast fluctuating function of the momentum k, the velocity of the inflaton in the non-adiabatic part of its potential, and the coupling between the inflaton and χ [28]. Preheating to fermionic decay products is also possible, due to similar non-adiabatic processes. However in this case, though the rate of particle production is parametrically faster than perturbative decays, the exponential growth in number density is absent [18,47].
If the rate of resonant energy transfer exceeds that of perturbative reheating, and happens sufficiently fast that the produced modes can reach large occupation number before being redshifted away from the resonant momentum band, preheating will have a significant effect. The details of resonance in an expanding universe are complex because the decrease in the inflaton oscillation amplitude from the expansion of the universe can lead to the parameter µ k in Eq. (3.1) changing significantly even during the period of a single inflaton oscillation. However, since particle production only take place during a small fraction of each oscillation, analytic results are possible, and these have been studied carefully in [28,48] (for subsequent literature see the reviews [49][50][51], and references therein).
As an example that can be analysed analytically, and is also compatible with current observations, we study a simple model with inflaton φ and potential where χ 1 is a real scalar in one sector, and χ 2 a real scalar in a second otherwise decoupled sector. A quadratic inflaton potential is commonly found in, for example, Chaotic Inflation models [52], and preheating in the single sector version of this model has been carefully studied in [28] where full details may be found. Rather than attempt to review all details, we simply quote the relevant physical results. The inflaton mass is fixed to m φ0 = 10 −6 M Pl , motivated by the COBE normalisation for a pure quadratic inflaton potential [53], and the masses of χ 1,2 are assumed to both be much less than m φ0 . For simplicity, we take the same initial conditions as [28] t 0 = π/ (2m φ0 ) and Φ (t 0 ) = 2M Pl / √ 3π 3 . While the number density of the produced states is sufficiently small that they do not significantly affect the inflaton equation of motion, Φ (t) = Φ (t 0 ) t 0 /t due to the expansion of the universe, where Φ (t) is the magnitude of the inflaton oscillations.

Single sector preheating
Properties of the resonances are determined by a parameter where Φ (t), and m φ (t) includes contributions from expectation values, which decreases with time. Since we are considering models with two different decay products, these will have separate q i determining the resonances to that sector. If q i 1, resonances are wide, with width ∆k ∼ m φ , and in an expanding universe µ k fluctuates randomly within one inflaton oscillation due to Φ (t) changing. This regime is known as stochastic resonance, and the number density of a mode is determined by the effective average value of a random parameter µ k (t). Meanwhile for smaller values of q i , but still 1, and at later times when the universe is expanding slower, the values of the µ k fluctuate less. The resonance bands remain wide, but the dynamics are closer to broad resonance in a non-expanding universe with the µ k constant and calculable over timescales longer than 1/m φ . In both the stochastic and broad resonance regimes particle creation dominantly occurs at times when φ (t) 0 and the evolution is non-adiabatic. For later use we also note that the typical momentum of modes produced is In contrast, if q i 1 resonance bands are narrow, with width parametrically ∆k ∼ q i m, and µ k ∼ q i /2 for the resonance with smallest k. However, for small q (t) energy is not transferred efficiently, since momentum modes are redshifted out of the resonance band faster than they are generated, preventing an exponential growth. In particular, this happens if q 2 i m φ > H [28], and for the model we consider H ∼ m φ /10, so fast energy transfer does not occur for q 1/4. 5 Since q (t) decreases, a theory starting with q (t 0 ) 1 will begin in the stochastic resonance regime, and then probably pass to broad resonances (unless q 1 is quickly reached, while the universe is still expanding fast). This transition is not important for our purposes since its effects can be absorbed in the averaged µ k values. Subsequently, the system will evolve into the narrow resonance regime, and shortly after exponentially efficient energy transfer will end, so evolution to a narrow resonance regime is a crucial event.
If the dominant effect reducing q i (t) is Hubble expansion decreasing Φ (t), the end of preheating happens at a time Alternatively, the decay products of the inflaton can themselves modify the properties of the resonance and potentially stop it. One important effect, called backreaction, is the generation of an expectation value χ 2 i , which affects the inflaton equation of motion. A careful analysis of the impact on the inflaton oscillations is needed since χ 2 i is itself time varying, and leads to χ 2 i n χi / (g i Φ) (where n χi is the total number density of χ i ) and Importantly for our work, backreaction can potentially lead to two distinct stages of preheating (both of which take place in the stochastic/broad regime).
Considering a model with only a single inflaton decay product, and q (t 0 ) > 1, first there will be a period when the contribution to the inflaton mass Eq. (3.6) is negligible. If q (t 0 ) is relatively small but still 1, this will be the only stage and q ∼ 1/4 will be reached before the number density n χ is large enough for the expectation value induced inflaton mass to be comparable to m φ0 .
For larger q (t 0 ), backreaction will be important and at a time t b leads to ∆m 2 φ m 2 φ0 . After this, the effective inflaton mass squared will increase exponentially fast, due to the exponential increase in n χ . Therefore the inflaton will oscillate much faster, and q (t f ) 1/4 will be reached in 10 oscillations of the inflaton for typical values of q 1 (t b ), which take place over a very short time , so at this time the energy is automatically shared equally between the inflaton potential, kinetic energy of the χ (and also the interaction energy ∼ g 2 Φ 2 χ 2 ) [28]. This is unlike if q = 1/4 is reached before backreaction, in which case m φ is independent of χ 2 , and the energy in χ states is much smaller.
Preheating could also be stopped by scattering effects where a large φ 2 develops due to the production of finite momentum φ states from interactions with the χ i . This would increase the mass of χ, and for models with large couplings could stop resonance before t b (we consider this possibility later). Depletion of the inflaton zero mode by scattering will also reduce the energy transfer rate, however typically the resonance is stopped by backreaction before this becomes important.

Two sector preheating and backreaction
Having reviewed this known physics, we now consider models with two sectors, and in particular if backreaction leads to effects between sectors. For simplicity we assume that g 1 g 2 and the equivalent averaged parameters for the fastest growing modes (that is, the ones that determine the eventual χ i number density) for each sector satisfy µ 1 µ 2 , where a subscript labels the sector, not momentum. The extension to other cases is straightforward. For couplings g i in the range 10 −4 ÷ 10 −3 typical values of µ are distributed randomly in the range 0.1 ÷ 0.14, with fluctuations up to µ ∼ 0.2 for some narrow windows of couplings. We also assume throughout that g i 10 −6 so that q i (t 0 ) 1 and at least some preheating can happen, otherwise the corresponding sector simply has no energy transferred during preheating.
The simplest scenario is if both the first and the second sector reach q i = 1/4 and finish preheating before backreaction is reached. This approximately happens if for both sectors. In this case, the energy transferred from the inflaton to the χ i is small, and the presence of another sector has no effect. During this stage of preheating, the number densities of χ i are given by [28] n χi (t) 8) and the time for q i = 1/4 to be reached for each sector is given by Eq. (3.5). Therefore the post-preheating relative number densities of the χ 1 and χ 2 is where the number density of χ 2 has been evolved from t H,2 to the later time t H,1 . Since this scenario assumes 3 × 10 5 µ i g i 12, the exponent can have a large, but not enormous, effect. The total kinetic energy in the state χ i , which is k 2 * ,i χ 2 i , can be compared to the remaining energy in the inflaton oscillations where parameters are normalised relative to g 0 = 3 × 10 −4 and µ 0 = 0.13. For comparison between models, both the kinetic energy density and the inflaton energy in Eq. (3.10) are evaluated at a time t b,0 , although the qualitative features of the result are not sensitive to this choice. 6 Due to the early termination of the resonance, only a small fraction of the inflaton energy is transferred during preheating in this scenario, unless a sector is very close to this boundary of reaching backreaction. The more interesting case is if a sector gets blocked by backreaction, which by assumption happens fastest for sector 1. We do not explicitly study the possibility that the dynamics in the two sectors are close enough to both contribute significantly to the blocking, but this is straightforward to include by taking Eq. (3.6) to have two sources.
If the second sector has a sufficiently small q 2 (t 0 ), resonance will stop before backreaction happens, in particular if t H,2 < t b,1 . The time for sector 1 to reach backreaction t b,1 can be well approximated by [28] and t H,2 is given by Eq. (3.5). t b,1 is only weakly dependent on g 1 , and µ 1 typically does not vary by more than a factor ∼ 2 at most. Therefore imposing t H,2 < t b,1 is almost equivalent to requiring t H,2 < t b,2 , that is the second sector would not reach backreaction if it was the only sector. In models in which this is satisfied, sector 1 gets an order 1 fraction of the energy density since it reaches q 1 = 1/4 in the backreaction regime. Meanwhile, the energy in sector 2 is given by Eq. (3.10), and is exponentially suppressed for small couplings. 7 If g 1 and g 2 are both relatively large 3 × 10 −4 , then q 2 (t b,1 ) 1/4 and the second stage of preheating will transfer energy to both sectors. Energy transfer to sector 1 continues until m φ increases enough that q 1 = 1/4. Meanwhile the energy transfer to sector 2 will stop earlier once q 2 = 1/4. During this process Φ (t) decreases due to the energy transfer to the χ i , and can be calculated using conservation of energy (this stage lasts much less than a Hubble time). However, during this time Φ (t) will only drop by a factor of 1/q 1 (t b,1 ) 1/4 which is ∼ few, and most of the decrease happens during the last one or two inflaton oscillations. Therefore, to get approximate analytic results we assume Φ has a constant value = Φ (t b,1 ) during this period.
The expectation value χ 2 1 increases proportional to e 4πµ 1 N , where N is the number of inflaton oscillations at this stage, and the q i decrease at the same rate. Consequently the total number of inflaton oscillations that each sector is preheated for during this stage N 1,2 are related by (3.12) Therefore, using Eqs. (3.8) and (3.11), the final relative number density of χ 1 and Provided µ 1 is not too different to µ 2 the exponentials are not enormous (since N 1 10). Therefore, if an analytic analysis is accurate, in this part of parameter space there can be significant differences between the number density of states in each sector, but not far beyond those that would come from perturbative decays. The total kinetic energy density in each sector is approximately k * ,i n χi . Since k * ,1 /k * ,2 (g 1 /g 2 ) 1/2 , and an order 1 fraction of the total energy finishes in χ 1 kinetic energy, the approximate energy fraction in the χ 2 sector is straightforwardly obtained from Eq. (3.13).
To summarise, if preheating is assumed to be stopped only by backreaction, then a sector that would be preheated efficiently if it was the only one in a theory will typically still be partially preheated. However, if another sector is more strongly coupled to the inflaton, the energy transferred to the less strongly coupled sector will be suppressed by powers of coupling constants (and potentially exponentials of O(1) numbers). This is in contrast to if the model only contained the less strongly coupled sector, in which case it would get an order 1 fraction of the total energy.

Rescattering and other effects
However, as mentioned backreaction is not the only effect that can end preheating. Another possibility, called rescattering, is that scattering of χ states off the inflaton zero mode produces a significant density of finite momentum inflaton states. These can induce an effective mass for χ i , which if it is larger than the typical resonant momentum k * ,i , prevents efficient preheating to the χ i sector. The computation of this process is subtle because only inflaton modes with momentum k k * ,i /4 (called hard modes) are distinguishable from the inflaton zero mode for the purposes of preventing the resonance. Again careful discussion may be found in [28]. Briefly reviewing the relevant physics, preheating to χ i is disrupted if the induced χ i mass is sufficiently large (3.14) The hard inflaton modes have kinetic energy ∼ 1 2 k 2 * ,i δφ 2 , and therefore the efficient initial stage of preheating stops when approximately 1/256 of the total energy in the system is in this form. After this time energy transfer will continue through more complex dynamics at a slower rate [55][56][57][58][59]. Following [28], since the energy of the hard modes comes from that of the χ i , a lower bound on the time at which this happens is the time t r at which 1/256 of the systems energy is in χ i kinetic energy. It can be straightforwardly shown that this corresponds to χ 2 1/256 Φ (t r ) 2 . From conservation of energy, at this time where Φ (t f ) is the amplitude the inflaton would have in the absence of rescattering if preheating finished at t f when q 1 = 1/4 due to backreaction. These estimates allow a simple, approximate, discussion of the effects of rescattering in a model with two sectors. During the final stage of preheating χ 2 1 grows ∼ e 4πµ 1 N where N is the number of inflaton oscillations after t b,1 . If preheating continued until backreaction was reached, then χ 2 1 Φ (t f ) 2 . Therefore, using χ 2 1/256 Φ (t r ) 2 and Eq.(3.15), rescattering will stop preheating to χ 1 after approximately second stage oscillations of the inflaton, where N 1 is the number of second stage oscillations that would occur in the absence of rescattering. Using Eq. (3.12), this means that is g 1 /g 2 3, preheating to χ 2 will still be taking place when rescattering blocks preheating to χ 1 . Note however, the actual dynamics will be complicated and potentially significantly different to this expression. Rescattering can therefore have interesting effects on the relative energy densities of the two sectors. If g 2 is not much smaller than g 1 , then χ 2 will continue to be preheated after resonance to χ 1 is stopped, because the mass induced by rescattering is ∼ g i whereas k * ,i ∼ √ g i . It is plausible that in a complete system the total energy transferred to χ 2 could be comparable to that transferred to χ 1 in such a case, despite the smaller coupling. However, our approximate understanding is far from this level of precision. Even if g 2 is substantially smaller than g 1 , if g 1 0.01 preheating will be ended by rescattering significantly earlier than backreaction would have stopped it [28]. At this time most of the energy will be in the form of inflaton-χ i interaction energy, and the final ratio of energy transferred to each sector depends on subsequent complicated dynamics. Further, rescattering will first block resonances with momentum k * , and efficient energy transfer could continue to higher momentum resonances. In Section 3.4, we use a lattice simulation to study processes at later times. If rescattering is effective, then at the end of preheating, any preheated states χ i are likely to have a large mass ∼ k * ,i (t r ). This will change due to later interactions of the inflaton, and would also redshift away, but could be important for subsequent dynamics. Using expressions for Φ (t r ) from [28], immediately after preheating ends m χ i ∼ √ g i × 10 15 GeV.
For reheating to complete, the energy remaining in the inflaton must be transferred to the sectors through perturbative decays, otherwise the universe will be matter dominated at late times. This requires trilinear couplings of the inflaton to χ i , in a full model. There are two especially motivated sizes of these couplings, given the quartic interactions of Eq. (3.2). If the inflaton potential is such that φ gets a VEV of the natural size φ ∼ m φ , the quartic interactions with χ will lead to trilinear interactions V ⊃ g 2 i m φ φχ 2 i . Alternatively, if the underlying theory is supersymmetric, and the potential comes from a superpotential, then there will automatically be terms ∼ g i m φ φχ 2 i [46]. In passing we note that if the theory comes from a superpotential W = m φ Φ 2 +g 1 ΦX 1 X 1 + g 2 ΦX 2 X 2 (where Φ, X i are chiral superfields with φ, χ i their θ = 0 components), then there will automatically be a large coupling between χ 1 and χ 2 that is ∼ g 1 g 2 χ 2 1 χ 2 2 . For all the values of couplings we consider, this will lead to fast thermalisation between the two sectors after preheating. Therefore in such models it is very hard to have a theory with separate sectors that both couple to the inflaton, but have different temperatures. Despite this, we continue to use the supersymmetry inspired relation between the quartic and trilinear couplings as a reasonable toy model for a theory with large trilinear couplings.
The presence of trilinear couplings can also affect preheating, or even be the dominant source of it. This has been studied analytically and numerically in [46], and can lead to interesting difference to the quartic case, since the states χ i can now get negative masssquared parameters (preheating is also altered if couplings between the inflaton and its decay products have negative sign [17]). In the case of a trilinear with size ∼ g 2 i , the quartic coupling is typically dominant throughout preheating. For a larger trilinear ∼ g i , the beginning of preheating is driven by the quartic interaction, but as the inflaton oscillation amplitude drops, the trilinear becomes increasingly important, eventually dominating the dynamics. Rather than attempt to analytically study a model containing multiple sectors coupled with large trilinear couplings to the inflaton, we instead simulate such a model numerically.
Finally, while we have focused on a simple model with well motivated parameter choices, it would also be interesting to study models with different inflaton masses. We expect the generic picture of backreaction and rescattering potentially leading to cross-sector effects to remains, but the details could be different. Other inflaton potentials might also lead to significant differences, though resonance is mostly sensitive to the dynamics of the inflaton around the minimum of its potential. More dramatic modifications are possible if the assumption that the two sectors are identical apart from the value of their coupling to the inflaton is relaxed, for example if there is a sector in which only fermions are coupled strongly to the inflaton. As mentioned, this can still undergo a form of preheating, and the dynamics of combining this with a preheated bosonic sector could be interesting.

Numerical analysis
Numerical simulations of preheating and the immediate aftermath can provide valuable information about a dynamical regime that cannot be easily studied analytically, for example on the energy transfer from the inflaton once backreaction and rescattering are important. Immediately after preheating the distribution of energy in both the inflaton and decay products sectors is dominantly in low momentum modes, compared to a thermal distribution with the same energy density. Scattering will move the system towards being thermalised, but the processes involved are complex [58][59][60][61][62][63][64][65]. Also a significant fraction of the total energy remains in the inflaton zero mode, and it is important for subsequent dynamics (discussed in Section 4) if this remains the case, or if scattering moves it into high momentum inflaton modes.
Previously, there have been many developments in understanding of the relevant physics utilizing lattice simulations, for example [46,[66][67][68][69][70][71][72]. 8 Further, recent progress has allowed for studies of dynamics at later times, and with finer grid spacing giving information on the behaviour of higher momentum states [54,73]. In order to determine the early stages of thermalisation, relatively large simulations are needed, and we do not carry out a full analysis here (in Section 4 we use information from [46,71], assuming the results will be similar in models with multiple sectors). Instead we implement models with multiple sectors in the code LATTICEEASY [74], to determine whether the analytic arguments based on the initial stages of preheating in Section 3.2 hold once the full dynamics of the system are included. The simulation is carried out on a grid of spatial size 128 3 (which is fairly small compared to the current best numerics). We run to a time t = 300/m φ by which point the fraction of energy in each sector is stable.
In Fig. 2 left panel, the model Eq. (3.2) is extended with a third scalar χ 3 coupled similarly to the inflaton. The total energy in each sector is plotted during preheating, which continues to times ∼ 100/m φ , including the kinetic energy, gradient energy, and also the energy due to interactions with the inflaton for each χ i . The interaction energy is on slightly different footing, since it will decrease as the inflaton expectation value drops. However, the amount of energy in this form is comparable to that in kinetic and gradient energies, so does not affect the results. At the start of preheating, transfer of energy happens to all three fields at a similar rate despite their different couplings, until a time ∼ 60/m φ . At this point χ 3 , which has the smallest coupling, stops preheating, while the other two sectors continue.  Left: An example of the evolution of the energy (compared to the initial energy of the system E init ) in each sector during preheating, for a model in which the inflaton can decay to three sectors. The inflaton potential is quadratic and couples to the other states by 1 2 g 2 i φ 2 χ 2 i . Sector 1 is fully preheated, and despite its smaller coupling sector 2 also gets a significant fraction of the energy. The coupling to sector 3 is just below threshold for efficient preheating and the fraction of energy transferred is much smaller. For clarity, the plotted energy is averaged over 3 inflaton oscillations. Right: The fraction of a system's total energy transferred to each sector in a two sector model during preheating, as the inflaton couplings to the two sectors are varied. The inflaton potential is the same form as in the left panel. In the purple shaded region, one or both sectors is efficiently preheated. Since the contours are close to straight lines, the presence of another sector does not significantly block preheating, even if it is more strongly coupled to the inflaton.
The coupling g 2 is such that χ 2 would get approximately 1/10 of the total energy during preheating if there were no other sectors present. This is still the case in the model including χ 1 shown in Fig. 2 left (and remains true even for theories with larger values of g 1 ). At the end of the simulation, the difference in energy between χ 1 and χ 2 is slightly smaller than that which would come from perturbative reheating. For larger g 2 , the energy transfer to χ 2 is increased, again approximately independent of the value of g 1 . Meanwhile the energy in χ 3 is highly suppressed, which is also the case in models containing only χ 3 given the value of g 3 . The low proportion of the energy in χ 3 , despite its coupling to the inflaton not being much smaller than that of χ 2 , shows the sharp threshold between efficient and inefficient preheating.
In the right panel, we plot contours for which the states χ 1 and χ 2 get a fixed fraction of the total energy, in a two sector model. The contours are almost horizontal and vertical over the parameter range accessible to our simulation, indicating that the presence of a second sector has a limited effect on the preheating of the first, compatible with the left panel. There is a minor increase in the value of coupling needed to get an order 1 fraction of the total energy density when a second sector is efficiently preheated. However this is significantly smaller than that predicted in the previous section assuming preheating is terminated by backreaction. In the limit that the second sector is decoupled from the inflaton, the values of q 1 (t 0 ) that lead to order 1, and suppressed, energy transfer are compatible with the analytic arguments, and also recent careful numerical analysis with large simulations [54,71]. Although the plot is shown up to couplings of g i = 10 −2 , close to this value numerical errors in the simulations become more problematic because typical resonance modes are higher momentum and further from the IR dynamics (more discussion can be found in [54]).
In the left panel, a significant fraction of the energy transfer from the inflaton to χ 2 happens at relatively late times, once the more strongly coupled χ 1 has almost finished preheating (at a time around 100/m φ ), and prior to this there is a short period where energy transfer halts. These features are common when a sector is relatively weakly coupled but still preheated, and can also appear in models containing a single decay product. A possible explanation for the smallness of the cross-sector effects is that a significant fraction of energy transfer to a fairly weakly coupled sector during preheating happens after rescattering has become important. Even though a more strongly coupled sector causes rescattering to happen earlier, this late time energy transfer might not be disrupted. As mentioned, rescattering will affect the more weakly coupled sector less than the more strongly coupled sector, and there is likely to be a significant proportion of energy remaining in the inflaton zero mode. Further investigation of the dynamics, for example by analysing the spectrum of the inflaton and its decay products, would be interesting.
We also consider models with relatively large trilinear couplings, and a potential where all fields are real scalars. The relation between the couplings follows the supersymmetry inspired pattern (but as discussed without a cross sector coupling). Repeating the numerical study we find that the cross-sector effects are very similar to those shown in Fig. 2 for the model with only quartic couplings. The presence of a second preheated sector slightly increases the coupling needed to get an order 1 fraction of the total energy, but not dramatically. We have also repeated the numerical simulation for a model with only trilinear inflaton χ i couplings, that is a potential Eq. (3.17) but without the g 2 i φ 2 χ 2 terms. In this case the results are again similar, even though the dynamics of preheating are different [46].

Post-preheating evolution
After the efficient energy transfer from the inflaton finishes, the remaining inflaton quanta will decay perturbatively, and effects at these times are crucial for the final relative temperatures of the hidden and visible sectors. In this section we discuss potentially relevant processes and estimate the parts of parameter space where they are significant. We stress that due to large uncertainty in the dynamics of this period even in the case of a single sector model, some of our results are necessarily rough approximations.
i . Models to the right or above the purple dashed lines have at least one sector efficiently preheated. Inflaton mediated thermalisation can happen at late times, when the temperature is m χ . For small g i reheating is purely perturbative, or there is a period of matter domination leading to a temperature ratio close to perturbative expectation. Early time thermalisation is possible if both sectors are preheated. If only one sector is preheated, and is not thermalised with the other, Bose enhancement of the late time inflaton decays can have a significant effect on the relative reheat temperatures. Right: Plot of the final temperature ratio between sectors in the part of parameter space where only one sector is preheated (solid lines). The model is the same as the left panel, and we assume that the preheated sector internally thermalises at t ∼ 10 4 /m φ , and that inflaton decays are perturbative after this and not efficient before. The Bose enhancement of inflaton decay can cause significant deviations from the perturbative prediction (dotted lines).
Aside from the dynamics we discuss, the relative temperatures of separate sectors could also be altered by changes in the number of relativistic degrees of freedom g at relatively late times, though the resulting effects are ∆T /T ∼ (∆g/g) 1/3 and typically not very large [75]. Alternatively, a non-relativistic scalar could dominate the energy of the universe at a time after reheating, and dominantly decay to a particular sector. This is a plausible scenario in models with flat directions, for example due to supersymmetry.
As an example we focus on the quadratic inflation model, Eq. (3.17), with SUSY inspired trilinear and quartic couplings to the other fields. Throughout we take m φ = 10 13 GeV, and χ 1 and χ 2 to have the same mass m χ , detailed study of other models would also be interesting, but is left to future work. A summary of the results for this model, with m χ = 10 6 GeV is shown in Fig. 3.

Late time thermalisation
For sufficiently large inflaton trilinear couplings, separate sectors will thermalise with each other at temperatures below the perturbative reheat temperature, mediated by inflaton scattering. If this is the case, at late times the two sectors will have comparable energy densities, regardless of earlier effects. The details of this computation have been studied carefully in [29], where explicit formula for the scattering rates are given. The relative energy transfer rate from sector 1 to sector 2 via inflaton scattering is defined as where ρ 1 is the energy density in sector 1, which is taken to be at higher temperature than sector 2. If Γ 1→2 (T 1 , T 2 ) exceeds the Hubble parameter at any time, the two sectors will thermalise.
In the model Eq. (3.17), the trilinear couplings are dimensionful, so thermalisation will happen most easily at the lowest temperatures satisfying T m χ . Provided m χ is not close to the predicted perturbative reheat temperature, both sectors will be close to a thermal distribution at this time, and at temperatures below the inflaton mass an approximate scaling is expected [29]. This dependence is only a rough estimate, and we obtain numerical results by computing the relevant scattering processes. In models with other types of couplings, for example an inflaton coupled to a fermion, thermalisation will be dominated by high temperatures. In Fig. 3 we show the parts of parameter space thermalised when the decay products both have a mass of 10 6 GeV (and the inflaton mass is m φ ∼ 10 13 GeV). 9 Because the dominant energy transfer happens at low temperatures, Bose enhancement of the cross section is not important. Consequently, the results are effectively the same if the constraint Γ 1→2 (T 1 , T 2 ) > H is imposed with T 1 = T 2 directly, rather than starting at early times with the post-reheating value of T 2 , and tracking its evolution with time. For the model we consider, a large part of the parameter space with at least one sector preheated is thermalised. The effect of changing the decay product masses in a full computation is well reproduced by the scaling Eq. (4.2) after setting T 1 = m χ (provided the reheat temperature is significantly below the inflaton mass), and the minimum value of the product of the coupling constants for thermalisation (with m φ = 10 13 GeV still fixed) is approximately If one of the sectors has typical mass scale much higher than the other, thermalisation is expected to dominantly occur when the hotter sector has a temperature close to the heavier mass scale (assuming scalars). However, a full computation is required for reliable results.
There may be interesting phenomenological possibilities if a sector which starts off relatively cold, but is warmed up by thermalisation at late times. In the model with m χ = 10 6 GeV, this only happens for the small part of parameter space in which only one sector is preheated, but thermalisation is efficient. However, it is possible over large parts of parameter space for smaller decay product masses.

Matter domination
If late time thermalisation does not occur, a possible period of matter domination will be the most important effect for the final relative temperatures of the two sectors. This happens if, after preheating, a significant population of non-relativistic inflaton states forms, and exists for an extended time before decaying perturbatively [18,71]. As a result, the energy transferred to the χ i sectors is redshifted away, and the late time temperatures are just set by the perturbative prediction. In particular, any mass the χ i states have from an expectation value at the end of preheating is expected to decrease as ∼ 1/a(t) where a (t) is the scale factor of the universe, and redshifting alone will not cause the χ i states to become more nonrelativistic if their bare mass is small. As discussed in Section 3, at the end of preheating the inflaton mass could also dominantly come from expectation values of thermal contributions if rescattering is important, but its bare mass is assumed larger than that of the χ i , so finite momentum quanta will become non-relativistic faster. Also, if interactions are not efficient at depleting the inflaton zero mode, this will immediately act as a matter component of the energy density. The trilinear inflaton couplings in the model we study, Eq. (3.17), are relatively large (for example compared to those in a model with couplings ∼ g 2 i m φ φχ 2 i ). This has two important effects. First it changes the perturbative inflaton decay rate, and the part of parameter space that is thermalised at late time, relative to the region that is preheated. The second, more subtle, effect is that numerical studies indicate that a large trilinear coupling causes the inflaton states after preheating to be far more relativistic than if the trilinear was absent [46].
In the case of a small trilinear, ∼ g 2 , lattice simulations show that from the time when preheating ends at t ∼ 100/m φ until times of 1000/m φ the fraction of χ states that are relativistic increases steadily to 1/3. Meanwhile, 1/5 of φ states are relativistic at t = 200/m φ , and this subsequently drops to 1/20 at t = 1000/m φ [46,71]. At these times the equation of state parameter w = p/ (where p is the pressure and the energy density), is w 0.15 and is decreasing towards the matter dominated value. More recent simulations run for longer, show that in such a model the fraction of energy in χ states drops approximately as ∼ 1/a (t) at late times. This shows that matter domination is reached relatively quickly, and barring extremely large values of g i at least some period of matter domination is inevitable. This is compatible with recent studies, which show that the matter dominated equation of state is quickly reached in a quadratic inflaton potential model without large interactions with decay products [76]. The dilution of the energy density in the χ i produced by preheating is given by the (t m /t) 2/3 , where t m is the time at which matter domination sets in, estimated as t m ∼ 250/m φ .
In contrast, models with large trilinear couplings ∼ g have different behaviour. Up to the times of the longest simulations carried out so far, t = 1000/m φ , the fraction of both inflaton and χ quanta that are relativistic are very similar, with both increasing to 1/2 steadily over the simulations [46]. 10 Meanwhile the equation of state is constant at w ∼ 0.25, and there is a relatively fast flow of the energy and number density spectrum towards UV modes, and kinetic equilibrium. However, the distribution of modes is still fairly IR localised compared to a thermalised system, which would have temperature 10 14 GeV. At these late times, the energy is still approximately equally distributed between the inflaton and an efficiently preheated sector, and therefore matter domination has not set in. We stress that by focusing on a model with large trilinears, we are choosing an option for which matter domination will have a relatively small effect, and in other cases it will be important over larger parts of parameter space.
Considering the model Eq. (3.17), an estimate for the parameter range for which matter domination is important can be made by assuming that both the inflaton and χ states continue moving towards a thermal distribution (numerical simulations would be useful to verify if this is accurate). Consequently, matter domination will occur if the timescale for perturbative inflaton decays exceeds the time taken for the energy in the inflaton sector to redshift to an effective temperature m φ , which happens at a time 10 5 /m φ (based on an equation of state close to the radiation dominated value). We make the simple assumptions that perturbative inflaton decay rate (including a possible Bose enhancement) is accurate, and that any significant expectation value induced and thermal masses are negligible, which is plausible given that this is happening at relatively late times. Energy could also be extracted from the inflaton by scattering processes, which would need to be included in a complete analysis, along with a study of the quasiparticle structure of the plasma [77].
Taking the potential of Eq. (3.17), and assuming the χ are thermalised, the inflaton decay rate is where the factor in brackets is the Bose enhancement from the χ number density [29]. This expression is obtained from the decoherence rate of a single inflaton mode [94], assuming that its spatial momentum is small [84] (which is a reasonable approximation since the majority of the quanta produced by rescattering have relatively small momentum). However, it does not include effects such as hydrodynamic slowdown, discussed in [56,57], which could lead to significant corrections, but are beyond the scope of our present analysis. As the temperature of the hottest χ i sector drops after preheating, there will be a time when Γ φ→χ i χ i (T i ) H (T i ) (the corresponding temperatures are significantly below the temperature immediately after preheating for typical values of g i ). For this to happen before matter domination, requires a coupling g i 0.005.
An alternative plausible assumption is that the distribution of decay products could remain close to that at the end of preheating (at a time t p ), apart from the effect of the expansion of the universe, corresponding to no significant flow of modes to the UV. If this is the case, the enhancement relative to the decay rate when there is no Bose enhancement, is a factor of f χ m φ 2 a(t) a(tp) , where f χ (k) is the occupation number distribution at a time t p . From lattice results, in this scenario the enhancement can be large. For example at t 1000/m φ , the occupation number distribution is approximately where k c is the comoving momentum of a mode defined relative to the beginning of preheating at a time t 0 m φ . In this case, decays are fastest compared to the Hubble parameter at early times due to the power dependence in Eq. (4.5). If the inflaton modes are thermalised, then provided decays are possible by a time 10 4 /m φ , matter domination will not occur if g i 0.001. In contrast, if there is no further flow of inflaton modes towards the UV, matter domination will happen earlier, within a time corresponding to a (t) /a (1000/m φ ) ∼ 2. Provided inflaton decays are not blocked, they could still be fast at these early time, however more complex dynamics, not well reproduced by perturbative calculations, might still lead to a period of matter domination. Another possibility is that the inflaton has a quartic potential in the part of field space where reheating happens, although a pure quartic potential is ruled out by observations (it may be possible to improve agreement with non-minimal couplings to gravity [78]). In this case, matter domination does not occur because the inflaton sector energy density evolves as radiation. This is supported by [54] who carry out a numerical simulation to relatively late times, finding an equation of state close to radiation and no relative decrease in the fraction of energy in inflaton decay products. Therefore, in models with a quartic potential, the final ratio of the temperatures of the sectors is determined by a combination of the energy distribution at early times, and the perturbative expectation from the order 1 fraction of energy remaining in the inflaton, if thermalisation does not take place. Recent studies [76] confirm this behaviour for any model in which the inflaton potential has form φ n with n 2.

Early time thermalisation
Another potentially important effect is inflaton mediated thermalisation at early times, enhanced by the high number densities immediately after preheating. Off-shell inflaton scattering in not included in the simulations of Section 3, but can be calculated directly using the distributions from numerical simulations, assuming perturbative calculations of the scattering rate are reasonably accurate (thermalisation by production of on-shell inflaton states will be seen in the numerics if it happens on short enough timescales). The computation is very similar to [29], with the difference that the energy distributions involved are very far from thermal, so that more integrals must be evaluated numerically. As before we compare the rate Γ 1→2 (T 1 , T 2 ) to the Hubble parameter (the rate at which the decay product's number distribution is altered by internal thermalisation could also give an important timescale, but this is unknown). If matter domination does not subsequently occur, early thermalisation will lead to similar final temperatures in the two sectors.
Since Bose enhancement is potentially large, we calculate the rate of energy transfer using occupation number distributions from the numerics [46], rather than assuming the two sectors already have equal energy density. A deficiency of this is that the number distributions are only giving in [46] for a model g i = 10 −3 , however many properties of preheating are not highly sensitive to the actual values of couplings provided they are above the threshold for efficient preheating. While we use the full distributions in our computations, as an indication the occupation numbers is approximately given by Eq.(4.5) at t 1000/m φ . We again neglect any expectation value sourced masses, and there remains the possibility that effects due to the high occupation number of χ states could lead to large deviations, so our results should be taken as indicative.
In parts of parameter space where one sector is preheated, the enhancement from the relatively high number densities is not enough to overcome the much larger value of the Hubble parameter at these times. As a result, in this scenario, early time thermalisation does not happen apart from for very large couplings, far inside the region where late time thermalisation is efficient. However, if both sectors are preheated, there in a potentially enormous Bose enhancement, even though the normal cancellation in the rate equation means it is only by one factor of the decay product occupation number not two [29]. Assuming the dynamics of the system are such that the thermalisation rate is well approximated by a perturbative calculation at t ∼ 100/m φ , this enhancement is as large as 10 6 for modes with physical momentum around k m φ /20 at this time (that is, modes with momentum around m φ at the beginning of preheating). This is large enough to cause thermalisation over all of the parameter space where both sectors are preheated efficiently.
Meanwhile, if thermalisation is not possible until a later time, for example due to large expectation value sourced masses, it will be suppressed. This is because, although the Hubble parameter decreases, the number density of modes with momentum not far below the mass of the inflaton decreases fast due to redshifting as well. For example, assuming Eq. (4.5) is valid over a range of times, the occupation number of modes with physical momentum around m φ drops as a (t) −3 . Supposing that the dynamics are such that thermalisation can first happen at a time t 1000/m φ , and both sectors get an order 1 fraction of the energy density, it occurs provided g 1 g 2 10 −3 . This is inside the late time thermalisation region, for m χ = 10 6 GeV.
As an aside, we note that in the case of a preheated bosonic sector potentially thermalising with a non-preheated fermionic sector, a computation shows that the thermalisation rate at the end of preheating is enhanced compared to a thermal distribution with the same energy density. This leads to an especially large change in the thermalised parameter space, as in such models thermalisation typically requires temperatures m φ [29]. Consequently, preheated thermalisation is stronger than thermalisation with same energy density, which itself is much stronger than thermalisation at the (relatively much lower) perturbative reheat temperature.

Bose enhanced perturbative reheating
Finally if none of these effects matter, the distribution of energy at the end of preheating can have an impact on the final ratio of reheat temperatures, beyond facilitating thermalisation. One possibility is that preheating could be the dominant source of energy for a relatively weakly coupled sector if perturbative inflaton decays are dominantly to another sector. This can happen because, provided its inflaton coupling is above a threshold, a sector typically gets an order 1 fraction of the total energy during preheating, even in the presence of a more strongly coupled sector. Alternatively, a sector could be preheated, but not have significant trilinear couplings to the inflaton, preventing perturbative inflaton decays (or the inflaton could dominantly decay to a non-preheated fermionic sector).
Another possibility is that this high number density after preheating could alter the relative rate of perturbative inflaton decays to the two sectors. 11 The energy in the inflaton could also be distributed differently from the case of pure perturbative reheating, since after preheating a substantial proportion is in higher momentum modes rather than the zero mode, which could change the dynamics although we do not investigate this. Even in perturbative reheating scenarios, the dynamics can be modified significantly by thermal and plasma effects [77,[79][80][81][82][83][84][85][86][87][88][89][90][91][92][93][94]. Thermal masses produced as a sector heats up could modify or block energy transport [81], and energy transfer by scattering can be important, as can the quasiparticle and collective excitation dynamics [77]. Quantifying these late time effects in a system with a complicated energy distributions left over from preheating is extremely challenging, even in the single sector case. A dedicated and computationally expensive numerical study is probably required, and this may not be within reach of current simulation power.
Rather than attempting a full analysis of preheated models with multiple sectors, we instead focus on the simple possibility that a Bose enhancement of the inflaton zero mode decays can change the relative reheat temperatures of two sectors. While not likely to be representative of the full dynamics of a model, our calculation shows that the ratio of final temperatures can be significantly different to the simple perturbative expectation g 1 /g 2 . In particular, we consider models in which only one sector is preheated, and therefore has a high number density. Unfortunately, even only considering the Bose enhancement, this depends on the unknown late time number distributions. However, it can be calculated if we assume that the states χ i reach a thermal distribution fairly quickly. Whether this is accurate will depend on the other unspecified interactions in a sector.
Assuming internal thermalisation, the temperature of the preheated sector will drop until the inflaton decay rate is equal to the Hubble parameter, at a temperature T d . The final temperature ratio is approximately given by the relative inflaton decay rates at this time 11 Thermal and plasma dynamics can alter the final temperature ratio from T1/T2 ∼ g 1/2 1 /g 1/2 2 even in the absence of preheating, however the effect is smaller.
In Fig. 3 (right), we plot the final ratio of the sectors temperatures over the part of parameter space in Fig. 3 (left) for which preheating happens to only one sector. Since sector 1 will not internally thermalise instantly after preheating, results are shown assuming this happens over a timescale 10 4 /m φ . We further assume that χ 1 has a large mass from an expectation value such that the inflaton does not decay to it before this time, and that by this time the inflaton energy is in modes that are fairly close to non-relativistic so the Bose enhancement is f χ1 (m φ /2). The simple perturbative prediction g 1 /g 2 for the relative temperatures of the two sectors is also shown for comparison.
For a particular point of parameter space, the ratio of temperatures can differ from g 1 /g 2 by a factor of ∼ 5 (that is, up to ∼ 1000 in the energy). The effect is most significant when the coupling g 1 is largest, and therefore the decay time is early when the first sector has a high temperature. We emphasise however that this is assuming the perturbative inflaton decay rate Eq. (4.4) is valid, and the high occupation density of inflaton modes may lead to complex modifications not capture by Eq. (4.6). Rather than claiming precise results, our point instead is simply that this is an effect that can cause large deviations from the expected temperature ratio.
If the comoving distribution of χ remains close to its form at the end of preheating, the relative change will depend on the occupation number at the time of decay. In this scenario, decays are more efficient at early times, due to the redshifting of the momentum distribution. The relative temperatures are fixed by the timescale over which decays from the inflaton are possible, for example due to a large expectation value source mass for the χ decreasing. This depends on the details of the complex dynamics, and could lead to larger temperature ratios than assuming thermalisation.

Non-Oscillatory models and large temperature asymmetries
In the models studied so far, large temperature ratios still require dramatic differences in sectors' couplings to the inflaton, although the ratio can be significantly different from the perturbative prediction. This is simply because backreaction and rescattering ends preheating when there is an order 1 fraction of the total energy left in inflaton modes. However, in Non-Oscillatory inflation models [31,[95][96][97] the inflaton potential does not have a minimum, and instead is constant at large field values. As a result, the inflaton does not oscillate and perturbative decays are absent. Early Non-Oscillatory models had inefficient reheating, leading to unrealistic cosmologies, however the Instant Preheating mechanism can lead to sufficient reheating [30,31]. In this scenario a state χ, with mass due to the inflaton expectation value, is resonantly produced. As the inflaton moves to the larger field values, the mass of χ increases, extracting more energy from the inflaton. If χ decays to other lighter states while heavy, sufficient energy can be transferred for successful reheating.
In this section we show that such models can lead to very large temperature differences from order 1 differences in the coupling constants. Analogously to Fig. 2, one sector could be above the threshold of resonance and reheated efficiently, whereas another similar sector just below threshold will have exponentially suppressed energy transfer. Of course, we do not claim this is the only, or the simplest, way to get large temperature asymmetries. For example perturbative inflaton decay to a sector could be forbidden if the inflaton is not coupled to any lighter states in that sector (although this constrains model building) or if it has no trilinear interactions [98]. In these scenarios the possibility of a sector being strongly preheated but not reheated might allow for interesting phenomenology [16].
We study a simple example model in which the inflaton couples to two scalars χ i in different sectors. These are themselves coupled to fermions ψ i in their own sectors, so that the energy transfer is φ → χ i → ψ i . The χ i must only couple strongly to the corresponding ψ i . The advantage, compared to having the inflaton itself coupled strongly to only one sector, is that χ i can easily be charged, for example under an unbroken (or weakly broken) gauge symmetry. This makes model building strong decays of χ i to only one sector straightforward. In contrast, attempting to make the inflaton part of the visible sector, in such a way that prevents it decaying to a generic hidden sector is more difficult [99].
For definiteness we assume the simple Non-Oscillatory model studied in [30,31] The inflaton is assumed to start at large negative field values, and moves towards φ = 0. At this point, the dynamics can become non-adiabatic allowing χ i to be efficiently produced. Analogously to normal preheating, the number density generated is whereφ 0 is the inflaton velocity near the origin. Following [30], in a model where inflation ends with φ −0.3M Pl , the inflaton speed at the origin is From Eq. (5.3), the production of χ i states is exponential suppressed if where the exponential suppression is because the theory is close to evolving adiabatically if this condition is satisfied. For a coupling of g i ∼ 1 the threshold in this model is m i ∼ 10 15 GeV.
Since particle production only happens once, when φ passes through the origin, this result is not affected by the presence of a second sector. Subsequently, φ evolves to larger values. Although n χi is constant, the mass of χ i grows ∼ g i φ(t), and the χ i energy density increases Meanwhile, energy remaining in the inflaton redshifts as a(t) −6 , because it is in the form of kinetic energy [95], so does not dominate the universe at late times. The produced χ i states affect the inflaton equation of motion, potentially making it move back towards φ = 0, which would cause ρ χi to decrease. We consider the simplest possibility, which is that the χ i decay to other lighter states before this happens. For definiteness, if χ i is coupled to significantly lighter fermions through an interaction κ i χ i ψ iψi , the decay rate is It can be shown that provided decays happen sufficiently fast that the inflaton does not change direction [30] (quintessence style models where the potential is not exactly zero for φ > 0 may be viable, and can relax this requirement). The total relative energy transferred to each sector (evaluated at late times) is therefore given by ρ 1 ρ 2 9) assuming matter domination between the decays of χ 1 and χ 2 , and making the approximation thatφ is constant until both the χ 1 and χ 2 have decayed, for simplicity. 12 It can be seen from Eq. (5.9) that if both sectors are in a part of parameter space where the number density Eq. (5.3) is not exponentially suppressed, they will get comparable energy density, with differences only from powers of coupling constants. In contrast, if one of the sectors is in the suppressed region it is colder with the energy ratio an exponential function of its coupling to the inflaton. The transition between these regimes only requires order 1 changes in the parameters of the two models, and also does not need tuning onto particular points in parameter space. 13 Given that these models typically involve relatively large couplings, both between the inflaton and the χ i , and between the χ i and ψ i , there is the potential for sectors to thermalise 12 The speed of φ will actually decrease, and φ = M Pl / 2 √ 3 log (t/t0) with t0 = 5/ √ 3πm φ due to the expansion of the universe [31]. For relatively large gi and κi, the effect of this will be small on the timescale of χi decays. For smaller coupling constants, it will be an order 1 difference to Eq. (5.9). However, the overall form of a relatively weakly varying function of the couplings multiplied by an exponential suppression below threshold remains. 13 The more strongly reheated sector could also be in the exponentially suppressed region, however in this case it must be fairly close to the boundary, so that sufficient energy is transferred to be reheated to a high enough temperature for acceptable phenomenology.
through inflaton scattering. However, thermalisation is only possible at very early times before the χ i decay, or when the typical energy of ψ i states is high enough to produce a population of χ i . Imposing that Eq. (5.8) is satisfied, a negligible fraction of the energy will be transferred to the other sector on these timescales. In particular, the mass of χ 1,2 is bigger than that of the inflaton, so any on-shell inflaton states produced cannot decay and transfer energy to these. The fraction of energy transferred by off-shell inflaton scattering on the relevant timescales, defined analogously to Eq. (4.1), is negligible as well. This is because the number density Eq. (5.3) is relatively low reducing scattering events compared to decays, and also the matrix element for s-channel scattering via an off-shell inflaton is suppressed since the center of mass energy is very large relative to the inflaton mass, while t-channel scattering is suppressed by the low number density of χ i in a sector for which reheating is exponentially small. Meanwhile, the inflaton continues moving to larger field values, so the masses of the χ i increase. Once these have grown by a factor ∼ 2, production of χ i from scattering of ψ i is no longer possible. This also happens on a timescale Γ −1 since the inflaton is moving at an approximately constant speed. Consequently, as long as the inflaton only has large couplings to the χ i , inflaton thermalisation mediated does not occur at later times either.

Discussion and phenomenological implications
For many beyond the SM scenarios the relative temperatures of two sectors at very late times is crucial. As an example we showed that for the DM relic abundance to be set purely by freeze-in, a temperature ratio of 10 is needed between the visible and dark sectors. Possible temperature asymmetries between sectors are also important for understanding cosmological bounds on new light states (for example, bounds on ∆N eff or energy injection to the visible sector). Low temperature hidden sectors might contain exotic states that could be detectable with future experiments, but would be ruled out by cosmological bounds in the absence of a temperature asymmetry [100,101], although thermalisation through the coupling to the Standard Model must also be considered [102]. Some formulations of the relaxion mechanism [103], which could solve the Electroweak hierarchy problem, also need large temperature differences between the visible and dark sectors [104].
In other models, the details of a thermal history of a sector, rather than its late time relative temperature, can have implications. Possibilities include restoration of a broken symmetry at high temperatures, which might only happen if a sector is preheated [105]. In this case, an understanding of the effects of multiple sectors on preheating is important. Similarly, whether preheating occurs can affect the production of very heavy, or very weakly coupled states [19], and potentially baryogenesis [106,107]. In another direction, interesting phenomenology might come from the possibility that a sector is initially cold after reheating, but is heated up by inflaton mediated thermalisation at late times.

Implications for dark matter
Preheating and related effects can have a significant impact on the DM relic abundance. This depends on the details of the particular DM candidate and production mechanism, and we simply consider a few simple scenarios, which by no means cover all the possibilities.
As a example suppose the DM is in a hidden sector uncoupled to the visible sector (except for inflaton mediated interactions), and that there is a light hidden sector state that the DM can annihilate to, similarly to the models of Section 2.2. The DM relic abundance is then given by Eq. (2.6), and is approximately proportional to the temperature ratio between the hidden and visible sectors at late times. This is determined by preheating, perturbative reheating, and possible inflaton mediated thermalisation. We also assume that the inflaton decays to scalars in the visible and hidden sectors through couplings of the form Eq. (3.17), and fix the coupling to the visible sector to g 1 = 0.03 so that the visible sector undergoes preheating.
The DM relic abundance in such a model is plotted in Fig. 4 as a function of the coupling of the inflaton to the hidden sector g 2 , for a DM candidate with mass m DM = 3 TeV and an s-wave annihilation cross section (σv) = 0.01/m 2 DM . The results are shown with the effects of preheating included (solid line), and for comparison the relic density that would be obtained if the effects of preheating was neglected is also plotted (dashed line). The energy transfered during preheating is calculated using the LATTICEEASY simulations of Section 3.4. In this model Bose enhancement of perturbative inflaton decays to the visible sector, combined with the large trilinear couplings, means that perturbative reheating happens relatively fast, at times t ∼ 1000/m φ . Consequently there is no period of matter domination. In computing the relic density we assume for definiteness that the energy in inflaton states redshifts as radiation between preheating and the time of perturbative decay, and the expansion of the universe is exactly radiation dominated at these times. However the results obtained are not very sensitive to these assumptions.
For very small values of g 2 less than approximately 10 −4 , the hidden sector is not preheated. In this case including preheating lowers the DM relic density relative to if its effects had been ignored. This is dominantly because the Bose enhancement increases the perturbative decay rate to the visible sector, decreasing the fraction of the energy transfered to the DM sector. The relative decrease in the hidden sector energy density is by a factor of approximately where Γ φ→visible (T ) is the inflaton decay rate at a temperature T , and T d is the visible sector temperature when perturbative decays dominantly occur, similarly to Eq. (4.6). The numerical result 0.01 is dependent on the particular model, and the magnitude of the Bose enhancement, and can vary from 1 ÷ 0.001 for 0.001 g 1 1 . Consequently the DM relic abundance is reduced by a factor of approximately ∼ 0.01 1/4 in the particular model we . The DM relic abundance Ω compared to the observed value Ω DM in the example model described in the text, as the inflaton coupling to the hidden sector g 2 is varied. The solid line shows the result including the effects of preheating, while the dashed lines show the relic density that would be found if preheating was ignored and only perturbative inflaton decays considered. The inflaton coupling to the hidden sector is taken to have the form Eq. (3.17), with a coupling constant g 2 . Meanwhile the inflaton coupling to the visible sector has the same form, with a fixed coupling constant g 1 = 0.03, so that this sector is preheated. Additionally, we assume the DM can annihilate to light hidden sector states, with an s-wave cross section (σv) = 0.01/m 2 DM , and that is has a mass of 3 TeV. For a given coupling to the inflaton, preheating can change the DM relic abundance by factors of order 10. consider (the actual decrease is slightly larger since preheating transfers an order 1 fraction of the inflaton energy to the visible sector, which is then unavailable for perturbative decays).
As g 2 is increased preheating to the hidden sector begins to occur and transfers a significant proportion of the inflaton energy, increasing the hidden sector temperature and the DM relic abundance. This corresponds to the sharp rise in Fig. 4, and leads to a larger relic abundance than if preheating was ignored. Additionally, once the hidden sector gets an order 1 fraction of the inflaton energy early time inflaton mediated thermalisation becomes efficient due to the high number densities in the two sectors, as studied in Section 4.3. When this happens the visible and hidden sector have the same temperature, apart from a small effect due to the differing perturbative inflaton decay rates (a fraction of which occur after the two sectors are not longer in thermal contact).
As a result, for g 2 0.003 the hidden and visible sectors have approximately the same temperature, and any further increase in g 2 has no effect on the DM relic abundance. Meanwhile, if preheating was neglected the hidden sector temperature would still be significantly below that of the visible sector, leading to a suppressed relic density. Finally for g 2 0.001 late time thermalisation is efficient, and the hidden and visible sectors have the same final temperatures even if preheating is neglected. This leads to a sharp increase in the DM relic density for the dashed curve for which preheating is ignored.
Preheating will also have a significant effect on the relic abundance in other dark matter scenarios. One simple case is if the dark matter sector has no number changing interactions, similarly to Section 2.1. If the dark matter sector is not preheated but the visible sector is, then the Bose enhancement of the perturbative decays to the visible sector results in production of fewer DM states. Analogously to Fig. 2 right and Eq. (6.1), Bose enhancement increases the inflaton decay rate to the visible sector by factors of up to 1÷1000 for couplings to the visible sector in the range g 1 0.001 (assuming trilinear couplings of the form Eq. (3.17)). In such models the DM relic abundance is simply proportional to the branching fraction of the inflaton to the hidden sector, and consequently is reduced by the same factor due to the Bose enhancement.
Another class of models in which preheating can have a dramatic effect on the relic abundance is if the DM is very heavy, known as "wimpzillas". One possibility is that superheavy DM is generated directly during preheating through its coupling to the inflaton [16][17][18]. In this case, in the absence of preheating there will be no DM relic abundance in parts of parameter space where the DM is too heavy to be produced by perturbative inflaton decays (production from the thermal bath is also extremely suppressed for such large DM masses).
Alternatively, DM could be very heavy and not be coupled to the inflaton, but instead produced from the thermal bath after preheating or reheating. In this case preheating still has a dramatic effect. For a DM mass above the maximum temperature the universe reaches after inflation T max , its relic density is parametrically suppressed by exp (−m DM /T max ) [109]. After preheating the effective temperature of the universe is very large, typically of order 10 15 GeV if the inflaton mass is ∼ 10 13 GeV. This can lead to significant production of even very heavy DM at this time, although the exact relic abundance depends on the details of interactions and thermalisation, which determine when production of DM from the thermal bath becomes efficient.
Meanwhile, in the absence of preheating the maximum temperature of the universe is parametrically T max ∼ H RH where H I is the Hubble parameter during inflation and T RH is the perturbative reheat temperature [108] (this temperature is reached before reheating completes). For a quadratic inflaton potential with energy density during inflation of m 2 φ M 2 Pl , and inflaton decay through a coupling to scalars gm φ φψ 2 , perturbative decays lead to a T max ∼ g 1/2 / (32π) 1/4 m Pl . To obtain the correct DM relic abundance in such models the exponential suppression must typically satisfy exp (−m DM /T max ) 1. As a result, including the effects of preheating leads to an exponentially large change in the relic density. Equivalently, the dark matter mass required to obtain the correct relic density is changed by a relative factor T PH /T pert , where T PH is the maximum effective temperature after preheating, and T pert the maximum temperature including only perturbative decays. This shift can easily be several orders of magnitude depending on the inflaton couplings in a particular model.

Summary
In this paper we have studied preheating and reheating when the inflaton is coupled to two otherwise decoupled sectors, and related effects such as inflaton mediated thermalisation between sectors. In particular, we have attempted to track the energy in separate sectors during each stage of the evolution of the universe, highlighting potentially important physical processes. If two sectors are both preheated this can result in final temperatures much closer than the perturbative prediction, either directly or due to efficient thermalisation, provided there is not a long period of matter domination before perturbative reheating. Meanwhile, if only one sector is preheated, thermal effects such as Bose enhancement of perturbative inflaton decays can alter the ratio of final temperatures.
Our analysis has two clear shortcomings. First, for many of the effects of interest, the underlying physics is complex. To make progress, we have had to make strong simplifying assumptions, and these are unlikely to fully reproduce the true dynamics of the system. However, despite being unable to obtain precise results, we have pointed out possibly interesting processes, the role of which could be clarified with further work. Second, we have restricted ourselves to very simple toy models, which are unlikely to be representative of the visible sector, or hidden sectors in interesting beyond SM physics scenarios. Although we have concentrated on scalars coupled to the inflaton, models coupled via fermions or gauge fields would also be worthwhile to explore [110]. Further, we have considered only sectors that are almost identical copies of each other, differing only in their coupling to the inflaton. Even with these simplifications, mapping out the dynamics in the g 1 vs g 2 plane is challenging. Again, further work to consider more general scenarios would be worthwhile.
Finally, we have noted that very large temperature asymmetries are possible in particular models of inflation and reheating, without requiring the inflaton coupling to different sectors to be vastly different. Large temperature differences could be interesting, for example in allowing heavy DM candidates that would have too large relic density in a warm sector, and there could be other model building possibilities.