Resonant Reheating
Basabendu Barman Nicolás Bernal Yong Xu



    MITP-24-045


Satellite group infall into the Milky Way: exploring the Crater-Leo case with new HST proper motions
    Mariana P. Júlio 1, 2 Marcel S. Pawlowski 1 Sangmo Tony Sohn 3 Salvatore Taibi 1 Roeland P. van der Marel 3,4 Stacy S. McGaugh 5
    Received 14 March 2024 / Accepted 23 April 2024
====================================================================================================================================






§ INTRODUCTION


Cosmic inflation is an elegant paradigm to solve several problems in cosmology <cit.>. In the conventional cosmological narrative, the initially (quasi) exponential expansion leads the Universe to be cold at the end of inflation. Understanding the transition from this state to a hot, thermal, radiation-dominated Universe, laying the foundation for Big Bang nucleosynthesis (BBN), is vital. This heating process not only elucidates the cosmic origins of the matter we observe, but also accounts for relics such as photons, neutrinos, and primordial gravitational waves (GW), along with exotic particles beyond the Standard Model (SM), for instance dark matter (DM).

The basic idea of reheating after inflation is as follows: it occurs due to the production of particles by the oscillating scalar field ϕ <cit.>. In the most basic inflationary models, this field, responsible for driving inflation, subsequently undergoes oscillations near the minimum of its effective potential V(ϕ) ∝ϕ^n, thus inducing particle production. Reheating can be performed through interactions between the inflaton and some light degrees of freedom. Such interactions, which result in energy transfer from the inflaton to the SM fields, can proceed through inflaton decays, for example, as discussed in Refs. <cit.> or inflaton scattering <cit.>.[The effect of time-dependent inflaton decay on reheating dynamics has been discussed in, for example, Refs. <cit.>.]

In the case of reheating through inflaton annihilations, contact interactions between a pair of inflatons and a pair of SM particles are typically assumed. Inflaton annihilations through a mediator have been studied in the framework of gravity, where the graviton plays the role of mediator. However, it has been shown that rather steep inflaton potentials with n > 9 are needed for successful reheating <cit.>. Interestingly, this bound can be relaxed to n > 4 if one introduces a non-minimal coupling between gravity and a pair of inflatons <cit.>, or even to n ≥ 2 if gravity couples non-minimally to a single inflaton <cit.>.

In this work, however, we investigate a previously unexplored novel possibility of reheating through scalar-mediated s-channel scattering.  In particular, we consider a model in which a pair of inflatons (ϕ) annihilate into a pair of daughter particles through the s-channel exchange of a massive spin-0 mediator (S). The daughter particles could be Higgs bosons in the SM or particles beyond the SM like right-handed neutrinos. If the inflaton oscillates in a potential steeper than that of the quadratic, the inflaton mass is not a constant but decreases with time. Consequently, the existence of higher modes can give enough energy to the inflaton so that mediators are produced resonantly. Once these mediator fields are produced, they quickly decay into SM particles that make up the thermal bath. Because of the resonant production of the mediator, the evolution of the SM temperature features non-trivial behavior (typically a bump) during reheating, which is very different from the standard cases characterized by a simple power law. We find that in the presence of a resonance, the maximum temperature during reheating is controlled by the location of the resonance relative to the end of the reheating. In particular, when the resonance occurs close to the end of reheating, the reheating temperature would be the maximum temperature and can be larger than that where the resonance is absent.

Interestingly, the mediator S could also play a fundamental role in the production of particles beyond the SM, such as the DM. In particular, if the DM has very suppressed couplings with the SM so that it never reaches thermal equilibrium, it can be produced in the early Universe through the FIMP mechanism <cit.>. The dynamics of reheating determines the abundance of DM if it is produced during reheating, and in that case it is called an ultraviolet (UV) FIMP <cit.>. We further explore the impact of such a resonant reheating on DM production from the visible sector or inflaton, mediated by S or graviton. We find that new and interesting freeze-in behaviors show up mainly due to the presence of the resonance feature during reheating.

Primordial gravitational waves (GW) serve as a robust prediction for inflationary models. In the course of inflation, quantum fluctuations inevitably lead to a spectrum of tensor-metric perturbations on scales beyond the Hubble horizon. In a conventional post-inflationary scenario, these tensor modes transition to subhorizon scales during the radiation-dominated epoch. However, it is widely recognized that the presence of nonstandard cosmological conditions, such as an early stiff era preceding radiation domination, disrupts this scale invariance <cit.>. Under such circumstances, the GW spectrum experiences a significant blue tilt within the frequency range corresponding to the modes that cross the horizon during the stiff period <cit.>. Such a blue-tilted spectrum turns out to be well within the reach of several GW detector facilities. Our scenario turns out to be testable on several next-generation GW detectors.

The paper is organized as follows. In Section <ref>, we offer the model setup and explore the features of resonance reheating with both numerical results and analytical approximations. In Section <ref>, we investigate the implications for the production of DM due to the presence of resonance. Section <ref> is devoted to discussions of the testability of the present framework with the primordial GW spectrum. We sum up our findings in Section <ref>.




§ POST-INFLATIONARY REHEATING DYNAMICS
 

After the end of cosmic inflation, inflatons oscillate around the minimum of their potential, transferring energy to SM particles that subsequently thermalize, forming the SM thermal bath. We consider the post-inflationary oscillation of the inflaton ϕ at the bottom of a monomial potential V(ϕ) of the form

    V(ϕ) = λ ϕ^n/M_P^n - 4 ,

where λ is a dimensionless coupling and M_P ≃ 2.4 × 10^18 GeV the reduced Planck mass. This type of potential could arise from α-attractor T- or E-models <cit.>, or the Starobinsky inflationary model <cit.>. Defining the energy density and pressure of ϕ as ≡1/2 ϕ̇^2+ V(ϕ) and p_ϕ≡1/2 ϕ̇^2 - V(ϕ), the background equation-of-state parameter is defined as w ≡ p_ϕ/ = (n - 2) / (n + 2) <cit.>. For example, for n=4, w = 1/3 corresponding to a radiation-like background.

The evolution of the energy densities of the inflaton  and the SM radiation  can be tracked using the Boltzmann equations[We also assume that the thermalization of the SM plasma occurs instantaneously. We note, however, that the SM particles do not necessarily thermalize suddenly, and thus they could be initially distributed with smaller occupation numbers and harder momenta <cit.>.]

    d/dt + 6  n/2 + n  H    = - 2  n/2 + n Γ^2→ 2 
    d/dt + 4  H    = + 2  n/2 + n Γ^2→ 2  ,

where 

    H = √(+/3 M_P^2) ,

is the Hubble expansion rate, and Γ^2→ 2 the 2-to-2 annihilation rate of a pair of inflaton condensate into a pair of SM states. The detailed expression of Γ^2→ 2 will be discussed in a moment.

During reheating, that is, when a_I ≤ a ≤, where a is the cosmic scale factor, the expansion H dominates over the interaction Γ^2→ 2, and then

    (a) ≃ (a_I) (a_I/a)^6  n/2 + n.

Here, a_I and  correspond to the scale factor at the end of inflation (that is, at the beginning of reheating) and at the end of reheating, respectively. Since the Hubble rate during reheating is dominated by the inflaton energy density, it follows that

    H(a) ≃ H() ×(/a)^3  n/n + 2    for  a_I ≤ a ≤ ,
    (/a)^2     for ≤ a ,

taking into account that after the end of reheating the Hubble expansion rate is dominated by free SM radiation.

The effective mass m_ϕ of the inflaton can be obtained from the second derivative of Eq. (<ref>), which reads <cit.>

    m_ϕ^2(a) ≡d^2V/dϕ^2 = n  (n - 1) λ  M_P^2 (ϕ/M_P)^n - 2≃ n  (n-1) λ^2/n  M_P^2 ((a)/M_P^4)^n-2/n.

In the last step, we have considered ≃ V(ϕ) and utilized the relation between the inflaton field ϕ and the potential given in Eq. (<ref>). The inflaton mass at a scale factor a can be related to its value m_I ≡ m_ϕ(a_I) at the beginning of the reheating by

    m_ϕ(a) ≃ m_I (a_I/a)^3 (n-2)/n+2.

It is interesting to note that for n ≠ 2, m_ϕ has a field dependence that will be inherited by the interaction rate Γ^2→ 2. Additionally, for n = 2, radiation cannot exceed the inflaton energy density, and therefore reheating is not viable if it proceeds through annihilations with n < 5/2 <cit.>.[Reheating through annihilations and n=2 could become viable in the case of an inflaton nonminimally coupled to gravity <cit.>.] In the following, we will therefore focus on the case n ≥ 4.




 §.§ Bosonic reheating


As advocated in the beginning, here we are interested in the 2-to-2 scattering of the inflaton ϕϕ→ S → h h, mediated by a real singlet scalar S in the s channel.[We note that if ϕϕ→ S → h h is present, one can not avoid the annihilation ϕϕ→ SS through t and u channels. Subsequently, the mediators S decay into a pair of SM particles S → hh. However, in the present setup, this channel is subdominant, as explored in Appendix <ref>.] This can arise from an effective Lagrangian density interaction of the form

    ℒ⊃1/2 μ_ϕ ϕ^2 S + 1/2 μ_h  S  h^2,

where ϕ is the inflaton field, S is the mediator and h a SM-like field.[We assume quartic interaction terms, e.g., ϕ^2 S^2, S^2 h^2, ϕ^2 h^2 or trilinear interaction between the inflaton and a pair of SM Higgs ϕ h^2, are absent. We remain agnostic about the UV completion of our effective toy model.] Note that both couplings have a mass dimension. Since the inflaton oscillates in a monomial potential, one must take care of the Fourier modes of the annihilating inflaton field. In that case, after averaging over several oscillations, the scattering rate for ϕϕ→ S → hh turns out to be

    Γ^ϕϕ→ hh≃ρ_ϕ/m_ϕ μ_ϕ^2 μ_h^2/32 π  m_ϕ^2 [ (4 m_ϕ^2 -m_s^2)^2 + Γ_s^2  m_s^2]√(1-m_h^2/m_ϕ^2) Θ(m_ϕ - m_h) ,

where S has a mass m_s and a total decay width 

    Γ_s = 1/8 π μ_h^2/m_s√(1-4 m_h^2/m_s^2) Θ(m_s - 2 m_h) + 1/8 π μ_ϕ^2/m_s√(1-4  m_ϕ^2/m_s^2) Θ(m_s - 2 m_ϕ) ,

with Θ being the Heaviside step function. Here, the effective coupling μ_ϕ is proportional to μ_ϕ and encodes the summation over several oscillations of Fourier modes <cit.>. It is interesting to note that the rate in Eq. (<ref>) features a pole when m_ϕ(a) = m_s/2. As the inflaton mass varies during reheating, this resonance could be reached. It occurs at the scale factor a_p given by

    a_p ≡ a_I (2  m_I/m_s)^2+n/3(n-2).

We emphasize that here we are interested in the case where n ≥ 4. Additionally, we will focus on the case where the resonance is reached during the reheating era.

This system of Boltzmann equations (<ref>) and (<ref>) can be solved analytically. Taking the initial conditions at a = a_I to be (a_I) = 0 and (a_I) = 3  H_I^2  M_P^2, with H_I being the inflationary scale,[The recent BICEP/Keck bound on tensor-to-scalar ratio implies that  H_I≤ 2.0× 10^-5 M_P <cit.>. This limit can also be used to derive an upper bound on the reheating temperature: ≲ 2.5 × 10^15 GeV under the assumption that the reheating is instantaneous <cit.>.] the SM energy density can be approximated by

    (a)    ≃H_I^3  M_P^4 μ_h^2 μ_ϕ^2/m_I^7
           ×9/512 πn/8 n - 17(a/a_I)^6(2n-7)/2+n[1 - (a_I/a)^2 (8 n-17)/n + 2]    for  a_I ≤ a ≪ a_p ,
    3 π/2n/n - 2(a_I/a)^4 (2 m_I/m_s)^n - 4/3 (n - 2)m_I^5/m_s^3 μ_h^2   for  a_p ≲ a ≲ a_t  ,
    9/32 πn/2 n - 5(a_I/a)^18/n + 2[1 - (a_I/a)^2(2 n - 5)/2 + n] (m_I/m_s)^4    for  a_t ≲ a ≤ .

Finally, after the end of the reheating, when a >, the SM radiation scales as free radiation and therefore (a) ∝ a^-4. Some comments on the analytical solution in Eq. (<ref>) are in order: i) before the resonance, that is, for a_I ≤ a ≪ a_p, the terms proportional to  and m_s in the denominator can be ignored, ii) near the resonance, the so-called narrow width approximation

    1/(s-m^2)^2 + m^2 Γ^2→π/m Γ δ(s - m^2)

is valid. It is worth mentioning here that although the resonance occurs at a particular epoch a=a_p, it persists over a finite time. And iii) after the resonance, one can safely take = m_I =0 in the denominator.

By equating the second and third solutions of Eq. (<ref>), one finds that the transition after the pole occurs at the scale factor a_t defined as

    a_t ≡ a_I [32 π^2/3 2n-5/n-2(m_I/μ_h)^2 (m_s/2 m_I)^2/3n-1/n-2]^2+n/2(2n-5).

Also, the end of reheating, defined as the onset of the radiation-dominated era, can be found by equating the third solution of Eq. (<ref>) with Eq. (<ref>), and the corresponding scale factor  is

    ≃ a_I (32π/32n - 5/nm_I^3  m_s^4/H_I  M_P^2 μ_h^2 μ_ϕ^2)^2+n/6(n-3).


Finally, the SM temperature T can be extracted from the SM radiation energy density  by the fact that

    (T) = π^2/30 (T)  T^4,

where (T) corresponds to the number of relativistic degrees of freedom contributing to . At a = a_p the SM bath reaches a maximum temperature  <cit.> given by

    ≃[45/π n/n-2(m_s/2 m_I)^n+4/n-2H_I^3  M_P^4 μ_ϕ^2/m_I^2  m_s^3]^1/4 .

Furthermore, at the end of reheating, when a =, the SM temperature  is

    ≃[90/π^2  H_I^2  M_P^2 (3/32 πn/2 n - 5H_I  M_P^2 μ_h^2 μ_ϕ^2/m_I^3  m_s^4)^n/n - 3]^1/4 .


In general, from Eq. (<ref>) the scaling of the SM temperature can be extracted, and it is given by

    T(a) ∝
            a^+3(2n-7)/2(n+2)   for  a_I ≤ a ≪ a_p ,
    
            a^-1   for  a_p ≤ a ≪ a_t ,
    
            a^-9/2(n+2)   for  a_t ≤ a ≪ .

We note that before the pole the SM temperature continuously increases with a slope T(a) ∝ a^+3(2n-7)/2(n+2), reaching the maximum T = at a=a_p. Interestingly, near the resonance, the temperature decreases as T(a) ∝ a^-1, independent of the value of n. In this case, the entropy injection and the Hubble expansion are counterbalanced, so that the SM radiation scales as free radiation. Well after the pole, the SM temperature continues to decrease, but more slowly, as T(a) ∝  a^-9/2(n+2) until the end of reheating at a =. After reheating, when the Universe becomes dominated by SM radiation, the SM entropy is conserved and therefore T(a) ∝ a^-1. Before moving on, we clarify that in the present setup there are six free parameters that dictate the reheating dynamics, and can be conveniently chosen to be

    n,  m_I,  H_I,  m_s, μ_h, and μ_ϕ ,

or

    n,  m_I,  H_I, a_p/a_I, a_t/a_I, and /a_I.

However, one can connect one set of parameters to the other using Eqs. (<ref>), (<ref>) and (<ref>). For simplicity, in all numerical evaluations of this section, we will fix m_I = 4 × 10^13 GeV and H_I = 2 × 10^13 GeV, as inspired by α-attractor models <cit.>.

The fully numerical solution of the system of Boltzmann equations (<ref>) and (<ref>), using the varying inflaton mass in Eq. (<ref>) together with the scattering rate in Eq. (<ref>), is presented in Fig. <ref> for n = 4, m_s = 10^8 GeV, μ_ϕ = 10^-2 GeV and μ_h = 10^6 GeV. The left panel shows the energy densities for inflaton (blue) and SM radiation (black), while the right panel shows the evolution of the SM bath temperature T as a function of the scale factor. The red dotted vertical lines correspond to a=a_p, a=a_t, and a=, while the red dotted horizontal lines (in the right panel) correspond to T=, T=T_t, and T=. Furthermore, the three dotted black lines in the left panel show the analytical solutions in Eq. (<ref>), which are in very good agreement with the numerical solution. 

 



 

In Fig. <ref> the dependence of the SM temperature on n is explored.  We show the evolution of T as a function of the scale factor for n=4 (solid line), n=6 (dashed line), or n=8 (dashed-dotted line). In the left frame, we fix the scale factors a_p/a_I = 10^5, a_t/a_I = 10^10, /a_I = 10^15 (we therefore change the mediator mass and the couplings) so that the different slopes of the SM temperature are clear. Their behavior fits well with Eq. (<ref>). In the right frame, we fix m_s = 10^8 GeV, μ_ϕ = 10^-2 GeV, and μ_h = 10^6 GeV.  In this panel, it is clear that, with an increase of n, the pole occurs earlier (as expected from Eq. (<ref>)), and consequently  is further enhanced due to the resonance.

In Fig. <ref> we show the evolution of the inflaton (blue) and SM radiation (black) energy densities for n=4 in the left panel, while the right panel shows the variation of the SM bath temperature T as a function of the scale factor.  For a mediator that is always heavier than the inflaton (solid black lines), cf. m_s ≫ m_I,  no resonance develops during reheating. Essentially, if the mediator is heavy, it can be integrated out, and the reheating dynamics corresponds to the scenario in which the inflaton annihilates through a contact interaction <cit.>. It corresponds to the third case explored in Eqs. (<ref>) and (<ref>). Similarly, if the mediator is always lighter than the inflaton (dashed black lines), the resonance cannot develop either. However, it is interesting to note that in this case, with a light mediator, the temperature increases monotonically until the end of the reheating and consequently the reheating temperature =. It corresponds to the first case in Eqs. (<ref>) and (<ref>). Note that the result presented in Fig. <ref> corresponds to a transition between the limit cases of the heavy mediator and the light mediator.

 





 §.§ Fermionic reheating


A successful reheating could also occur through inflaton annihilation into a pair of vector-like fermions ψ or e.g. right-handed neutrinos with mass m_ψ, mediated by S in the s channel: ϕϕ→ S →ψψ. This can arise from the effective Lagrangian density

    ℒ⊃1/2 μ_ϕ ϕ^2 S + y_f  S ψ̅ ψ .

The corresponding annihilation rate reads

    Γ^ϕϕ→ψψ≃ρ_ϕ/m_ϕ1/8πμ_ϕ^2  y_f^2/(4 m_ϕ^2 - m_s^2)^2 + m_s^2 Γ_s^2(1-m_ψ^2/m_ϕ^2)^3/2,

and the total decay width for S

    Γ_s = y_f^2  m_s/8 π(1-4 m_f^2/m_s^2)^3/2Θ(m_s - 2 m_f) + 1/8 π μ_ϕ^2/m_s√(1-4  m_ϕ^2/m_s^2) Θ(m_s - 2 m_ϕ) .

The system of Boltzmann equations (<ref>) and (<ref>) can be analytically solved following the same procedure as in the previous section. Consequently, the solution for the SM radiation energy is given by

    (a) ≃9/128 π H_I^3  M_P^4  y_f^2 μ_ϕ^2/m_I^5×n/5n-11(a/a_I)^6(n-5)/n+2[1 - (a_I/a)^2(5n-11)/2+n]   for  a_I ≤ a < a_p,
    n/n-1(2 m_I/m_s)^4 (a_I/a)^4 [1 - (a_I/a)^2(n-1)/n+2]   for  a_p < a ,

where a_p is again the scale factor where the pole is reached (i.e. m_ϕ(a_p) = m_s/2), given in Eq. (<ref>), and corresponds to a maximal temperature of the SM plasma

    ≃[135/64 π^3  n/n-1 H_I^3  M_P^4  y_f^2 μ_ϕ^2/m_I^5(2 m_I/m_s)^8(n-4)/3(n-2)]^1/4 .

The analytical solution for the SM energy density in the fermionic case in Eq. (<ref>) can be divided into two cases (instead of three as in the bosonic case; cf. Eq. (<ref>)), as the scaling of  during and after the pole is the same, (a) ∝ a^-4.

The end of reheating, defined as the moment at which the Universe transitions to a radiation-dominated epoch, happens at a scale factor

    ≃ a_I (8 π/3 n-1/n m_I  m_s^4/H_I  M_P^2  y_f^2 μ_ϕ^2)^n+2/2(n-4),

corresponding to the temperature

    ≃[90/π^2   H_I^2 M_P^2 (3/8 π n/n-1 H_I  M_P^2  y_f^2 μ_ϕ^2/m_I  m_s^4)^3n/n-4]^1/4 .

Finally, we observe from Eq. (<ref>) that the SM temperature scales as

    T(a) ∝
            a^3(n-5)/2(n+2)   for  a_I < a < a_p ,
    
            a^-1   for  a_p < a .


The fully numerical solution of the system of Boltzmann equations for the fermionic case is presented in Fig. <ref> for n = 6, m_s = 10^8 GeV, μ_ϕ = 10^-2 GeV and y_f = 10^-2. The left panel shows the energy densities for the inflaton (blue) and SM radiation (black), while the right panel shows the evolution of the SM bath temperature T as a function of the scale factor. The red dotted vertical lines correspond to a=a_p and a=, while the red dotted horizontal lines (in the right panel) correspond to T= and T=. Furthermore, the two dotted black lines in the left panel show the analytical solutions in Eq. (<ref>), which is in very good agreement with the numerical solution. 

 





§ DARK MATTER FREEZE-IN DURING RESONANT REHEATING
 

The DM number can then be tracked by solving the Boltzmann equation for its number density 

    d/dt + 3  H  = γ ,

where γ is the DM production reaction rate density. As the SM entropy is not conserved during reheating due to the annihilation of the inflaton in SM particles, it is convenient to introduce a comoving number density ≡  a^3, and therefore Eq. (<ref>) can be rewritten as

    d/da = a^2 γ/H ,

which has to be numerically solved together with Eqs. (<ref>) and (<ref>), taking the initial condition (a_I) = 0. To fit the whole observed DM relic density, it is required that

    Y_0  = Ω h^2  1/s_0 ρ_c/h^2≃ 4.3 × 10^-10 GeV,

where Y_0 ≡ Y(T_0) and Y(T) ≡(T)/s(T), with s being the SM entropy density defined by

    s(T) = 2π^2/45 (T)  T^3,

and (T) being the number of relativistic degrees of freedom contributing to the SM entropy. Furthermore, ρ_c ≃ 1.05 × 10^-5  h^2 GeV/cm^3 is the critical energy density, s_0≃ 2.69 × 10^3 cm^-3 the present entropy density <cit.>, and Ω h^2 ≃ 0.12 the observed abundance of DM relics <cit.>.


In the early Universe, DM can be produced from 2-to-2 scattering of SM thermal bath particles or inflatons, mediated by S.[We note that the inflaton annihilation into a pair of S particles than subsequently decay into DM states is also viable, however it is subdominant in our scenario as discussed in Appendix <ref>.] Additionally, it can also be mediated by the exchange of gravitons. In the following subsections, the DM production will be discussed in detail.




 §.§ From the thermal bath


In this section, we investigate the impact of resonant reheating, presented in the previous section, on dark matter (DM) yield. We consider the simplest scenario, where the DM is a real singlet scalar φ of mass , with an interaction Lagrangian that reads

    ℒ_ dm⊃1/2 ^2 φ^2 + 1/2  S φ^2 ,

assuming that the only mediator between the SM, the inflaton and DM is S, and again ignoring quartic couplings. Here, we consider any explicit coupling between the DM and the inflaton, and quartic interactions between the mediator S or the SM Higgs doublet and the DM are absent by construction. We again remain agnostic about the UV completion of such scenarios and focus mainly on the aftermath of resonant reheating on DM production. The only interaction that the DM can have with the visible sector is mediated by S, as can be seen in Eq. (<ref>). Therefore, DM can be produced in the early Universe from the bath via the scattering of SM Higgs, mediated by S. This process can occur during (inflaton domination) and after (radiation domination) reheating; however, the effect of resonant reheating can only be seen in the former scenario. Now, if the interaction strength between the DM and the visible sector is sufficiently weak to ensure that the DM never thermalizes with the radiation bath, then the production of DM takes place through the freeze-in paradigm. This is the approach that follows.

As mentioned above, the DM yield is governed by the 2-to-2 scattering of the bath particles (Higgs), mediated by the scalar S. The corresponding cross-section for this process reads

    σ(s)_hh→φφ≃1/8 π s( μ_h)^2/(s-m_s^2)^2+Γ_s^2  m_s^2√(1-4^2/s) ,

where Γ_s is given by (following Eqs. (<ref>) and (<ref>))

    Γ_s = 1/8 π μ_h^2/m_s√(1-4 m_h^2/m_s^2) + 1/8 π μ_ϕ^2/m_s√(1-4  m_ϕ^2/m_s^2) + 1/8 π ^2/m_s√(1-4 ^2/m_s^2) ,

where the last term corresponds to S decay into a pair of DM particles.
In the left panel of Fig. <ref> the production cross section is shown as a function of the center-of-mass energy s, for m_s = 10^8 GeV, = 10^6 GeV, μ_h = 10^6 GeV and μ_ϕ = 10^-2 GeV. The resonance at √(s)≃ m_s is clearly visible as a peak in the cross section.

 


One can further define the interaction rate density as

    γ = γ_hh→φφ = T/2π^4 ∫_4^2^∞ ds s^3/2 σ(s)_hh→φφ K_1(√(s)/T),

which takes the approximate form for ≪ m_s

    γ_hh→φφ≃T/2π^4 1/8 π( μ_h)^2/1×4 ^2  T/m_s^4  K_1^2(/T)     for  T ≪m_s/2,
    π/  K_1(m_s/T) √(1 - 4 ^2/m_s^2)    for  T ≫m_s/2 ,

or in the opposite case where m_s ≪

    γ_hh→φφ≃T/2π^4 1/8 π( μ_h)^2/1×π  T^2/8 ^3  e^-2/T    for  T ≪,
    T/6 ^2    for  T ≫ ,

with K_1(x) being the modified Bessel function of the first kind. Here, we have ignored the masses of the initial states, a reasonable approximation during reheating, which occurs at a very high temperature, much before the electroweak symmetry is broken. The right panel of Fig. <ref> shows the production rate density as a function of temperature T, for m_s = 10^8 GeV, = 10^6 GeV, μ_h = 10^6 GeV and μ_ϕ = 10^-2 GeV. The black dotted lines correspond to the analytical solution in Eq. (<ref>), in good agreement with the fully numerical solution. Interestingly, the sharp peak in the cross section appears as a smooth bump in the rate density.

In Fig. <ref> examples of the evolution of the comoving DM abundance as a function of the scale factor a are shown for n=4, a_p/a_I = 10^5, a_t/a_I = 10^10, /a_I = 10^15, and = 10^-4 GeV and = 2× 10^-3 GeV (top left), = 10^0 GeV and = 5× 10^-2 GeV (top right), or = 10^4 GeV and = 3× 10^-1 GeV (bottom). In this section, we take m_I = 10^12 GeV and H_I = 5 × 10^8 GeV, to avoid gravitational overproduction, as will be seen in Section <ref>. The vertical red dotted lines correspond to the scale factors for which T=, T=T_t, or T=, while the dashed red lines correspond to T =. All there benchmark points fit well the observed relic abundance, represented by a horizontal blue band. The detailed evolution of these benchmark points will be described in the following. 

 


For light DM particles, with masses below  (and m_s), the heavy mediator S can be safely integrated out so that the cross section in Eq. (<ref>) effectively scales as σ(s) ∝μ_h^2^2/(s  m_s^4), and therefore the interaction rate density γ(T) ∝μ_h^2 ^2  T^4/m_s^4; cf. the first line of Eq. (<ref>) in the limit ≪ T. In this case, DM is a standard IR FIMP, produced at low temperatures, mainly when T ∼, as shown in the upper left panel of Fig. <ref>. The approximate analytical solution of Eq. (<ref>) reads

    Y_0 ≃135/64π^7 √(10/) ^2 μ_h^2  M_P/  m_s^4

and, as expected, is independent on the reheating dynamics. Interestingly, to fit the total observed relic abundance of DM, cf. Eq. (<ref>), the required inflaton coupling to DM is independent of the mass of DM. This can be seen in Fig. <ref>, where the thick black line shows the parameter space in the plane [, ] that fits the entire observed abundance of DM.

 


The second case corresponds to intermediate masses of DM in the range ≲≲ T_t. This case resembles the previous one, with the difference that in the range a_t < a < the temperature scales as T(a) ∝ a^-9/2(n+2), cf. Eq. (<ref>). DM is constantly produced in the range T_t ≳ T ≳, and stops when the temperature is of the order of its mass, due to kinematical reasons. After production stops, during the range ≳ T ≳, the injection of entropy due to the decay of the inflaton into the SM states dilutes the abundance of DM. In this regime, the approximate analytical solution of Eq. (<ref>) reads

    Y_0 ≃135/8π^8 √(10/) Γ[6(n-2)/n+2]/Γ[7n-10/n+2]^2 μ_h^2  M_P/m_s^4 (/)^4(n-2)/3,

showing a strong dependence on the reheating dynamics, as expected in the case of a UV FIMP. To fit the observed DM abundance, the inflaton coupling to DM has to scale as

    ∝^4n-11/6,

as shown in Fig. <ref> in the range ≲≲ T_t.

The third case corresponds to heavy DM, with mass in the range T_t ≲≲ m_s. The interaction rate density can be approximated as γ(T) ∝^2  m_s^2  (T/m_s)^3/2  e^-m_s/T, cf. the second line of Eq. (<ref>) in the limit T ≪ m_s. Several comments are in order: i) the exponential factor stops DM production for temperatures below the mass of S, ii) since in this period the SM radiation behaves as free radiation (cf. Eq. (<ref>)), there is effectively no entropy dilution and, therefore, the DM yield is constant, and iii) there is, however, a dilution of the abundance of DM in the era T_t > T >. An example of the evolution of the DM yield for this case can be seen in the lower panel of Fig. <ref>. An approximate analytical solution of Eq. (<ref>) is

    Y_0 ≃135/4√(2) π^11/2 √(10/) ^2  m_s^2  M_P/^5(a_p/ /m_s)^6(n+1)/n+2.

We note that as the rate density is independent of the DM mass, it is required that 

    ∝^-1/2 ,

to fit the observed DM abundance, as shown in Fig. <ref>

Finally, we also note that when > m_s/2, the DM production becomes Boltzmann suppressed, requiring a sudden increase of , as shown in Fig. <ref>.
Overall, it is interesting to see that this behavior significantly differs from the freeze-in during standard reheating scenarios; see, e.g. Refs. <cit.>.




 §.§ From inflaton annihilations


In our setup, another possible DM production channel occurs through the inflaton annihilation via an s-channel exchange of the mediator S, i.e. ϕϕ→ S →φφ. The production rate density is given by[We note that a factor 1+ω≡2n/2+n shall be included on the right handed side of the Boltzmann equation for the DM number density.]

    γ = γ_ϕϕ→φφ = + 2  n/n+2 Γ_ϕϕ→φφ /m_ϕ ,

with

    Γ_ϕϕ→φφ≃ρ_ϕ/m_ϕ μ_ϕ^2 ^2/32 π  m_ϕ^2 [ (4 m_ϕ^2 -m_s^2)^2 + Γ_s^2  m_s^2]√(1-^2/m_ϕ^2) Θ(m_ϕ - ) .



 

Examples of the evolution of the comoving DM abundance N (left) and Y (right) as a function of the scale factor a are shown in Fig. <ref> for n=4, a_p/a_I = 10^5, a_t/a_I = 10^10, /a_I = 10^15, and = 10^-5 GeV and = 2 × 10^2 GeV (top), = 10^-1 GeV and = 25 GeV (second from top), = 10^4 GeV and = 2 GeV (third from top), or = 10^9 GeV and = 7× 10^5 GeV (bottom). The vertical red dotted lines correspond to the scale factors for which T=, T=T_t, or T=, while the dashed red lines correspond to m_ϕ =. All four benchmark points fit well the observed relic abundance, represented by a horizontal blue band (right). The detailed evolution of these four benchmark points will be described in the following. 

In this scenario, the abundance of DM can be analytically estimated in different regimes, depending on its mass with respect to the varying inflaton mass. Light DM, with mass < m_ϕ(), is produced mainly at T ≃ (top panel of Fig. <ref>). This corresponds to a UV freeze-in during radiation dominance, with an interaction driven by a heavy mediator. An approximate analytical solution of Eq. (<ref>) is

    Y_0 ≃135/64 π^3  n/n-3 ^2 μ_ϕ^2  H_I^3  M_P^4/m_I^4  m_s^4 ^3(a_I/)^3(8-n)/2+n.

To fit the total observed relic abundance of DM, cf. Eq. (<ref>), the required inflaton coupling to DM has to scale as

    ∝^-1/2,

as shown in Fig. <ref> in the plane [, ].

 


The second case corresponds to DM with mass m_ϕ() <  < m_ϕ(T_t). The bulk of DM density starts to be produced at T=T_t and stops when m_ϕ =, due to kinematics, at a scale factor a_dm given by

    a_dm≡ a_I (m_I/)^n+2/3(n-2),

as shown in the second from top panel in Fig. <ref>. Between m_ϕ = and T = large entropy injection dilutes the DM abundance; the final DM yield is estimated as

    Y_0 ≃135/64 π^3  n/n-3 μ_ϕ^2 ^2  M_P^4  H_I^3/m_I^4  m_s^4 ^3(a_I/)^3 (m_I/)^2(n-3)/n-2,

which implies that, to fit the observed DM abundance, the inflaton coupling to DM has to scale as

    ∝^n-4/2(n-2).

For the case n=4,  is independent of the DM mass, as shown in Fig. <ref>.

The third case corresponds to DM with mass m_ϕ(T_t) <  < m_ϕ(T_p) = 2  m_s produced mainly at = 2  m_s, when T =, which is a clear example of a UV freeze-in during reheating; cf. the third panel from the top in Fig. <ref>. Here, the interaction rate density in Eq. (<ref>) is dominated by the pole. The DM abundance is then diluted during the period T_t > T >, the final yield being 

    Y_0 ≃135/16 π^2  n/n-2 μ_ϕ^2 ^2  M_P^4  H_I^3/m_I^2  m_s^5 ^3 (a_I/)^3 (m_s/2 m_I)^2/n-2,

implying that

    ∝^-1/2

to fit the observed DM abundance, as seen in Fig. <ref>.

Finally, the fourth case corresponds to the heavy DM with mass 2  m_s <  < m_I, which is produced mainly around = m_ϕ, when the scale factor is given by Eq. (<ref>). The abundance of DM suffers two periods of dilution, from m_ϕ = to T=, and then during the interval T_t > T >; see the lower panel of Fig. <ref>. The final DM yield is therefore 

    Y_0 ≃135/1024 π^3  n/3 n - 7(a_I/)^3 H_I^3  M_P^4 ^2 μ_ϕ^2/m_I^8 ^3(m_I/)^2(3n-7)/n-2,

implying that

    ∝^5n-12/2(n-2)

to fit the observed DM abundance, as seen in Fig. <ref>. We note that DM particles heavier than the inflaton (> m_I) are kinematically forbidden in this channel.




 §.§ Gravitational production
 

Gravity is the only interaction that is guaranteed to mediate between the DM and the visible sector, and therefore it has to be taken into account. DM can indeed be produced by the scattering of inflatons or SM particles, mediated by the s-channel exchange of gravitons. Below we briefly discuss the DM gravitational production.




  §.§.§ From the thermal bath


DM is unavoidably produced from the UV freeze-in mechanism via 2-to-2 annihilation of the SM particles mediated by s-channel exchange of gravitons <cit.>. The interaction rate density for such a process reads

    γ = α T^8/M_P^4

with α≃ 1.9 × 10^-4 for real scalar DM <cit.>.
In the cosmological scenario presented here, DM is produced mainly at T = and then diluted during the period T_t > T >.
An approximate analytical solution of Eq. (<ref>) is

    Y_0 ≃45 α/4π^2  n+2/n+5 ^8/H_I  M_P^4 ^3(a_I/)^3 (a_p/a_I)^6(n+1)/n+2.





  §.§.§ From inflaton annihilations


During reheating, the whole observed DM abundance can be generated via 2-to-2 annihilations of inflatons, mediated by the s-channel exchange of gravitons. The interaction rate density for DM production reads <cit.>

    γ = 2 n/n+2 ^2/256π  M_P^4  f(/m_ϕ),

with
    f(x) = (x^2+2)^2 √(1-x^2)

for real scalar DM <cit.>. The bulk of the DM relic density is promptly produced at the beginning of the reheating era (a = a_I), and suffers two periods of dilution: a_I > a > a_p and a_t > a >. The final DM yield can be analytically estimated as

    Y_0 ≃135/256π^3  2n/n-1(H_I/)^3 (a_I/)^3 [1 - (a_I/)^6(n-1)/n+2].

Gravitation production from inflaton annihilation typically dominates over the annihilation from thermal-bath particles and can even overtake production through the exchange of S for very large values of H_I and .



 

Figure <ref> is an example of the combination of all previously described DM production mechanisms, for n=4, a_p/a_I = 10^5, a_t/a_I = 10^10, /a_I = 10^15, m_I = 10^12 GeV, and H_I = 5 × 10^8 GeV. It shows the production through annihilations of SM particles in the thermal bath (dashed blue; see Fig. <ref>), inflatons (dotted blue; see Fig. <ref>), gravitational production (sharp cut at ≃ 2× 10^8 GeV), and the combination of all of them (thick black). For this particular choice of parameters, the production from the SM thermal bath dominates for low masses. The production from inflaton scatterings is comparable near the pole and dominates in a small region around ∼ 10^8 GeV. For higher masses the gravitational production from inflaton annihilations dominates and overcloses the Universe for even higher DM masses. In summary, in this section, we find that, within the framework of resonance reheating, new features for DM UV freeze-in show up primarily due to the resonant temperature behavior near the pole.




§ PRIMORDIAL GRAVITATIONAL WAVES



The resonant reheating scenario described so far can have a possible signature at the next-generation gravitational wave (GW) detectors, courtesy to the blue-tilted primordial GWs with an inflationary origin. In this section we discuss such a prospect. 

The spectrum of GWs is described in terms of the
fraction of their energy density per logarithmic wavenumber (frequency) k interval

    (a,k) = 1/ρ_c d/dlog k ,

normalized to the critical density ρ_c=3 H^2  M_P^2. This expression can be recasted as <cit.>

    (a ,k) = 1/12(k/a  H(a))^2 𝒫_T,prim 𝒯(a,k) ,

where the primordial tensor power spectrum 𝒫_T,prim can be parametrized as <cit.>

    𝒫_T,prim = r 𝒫_ζ(k_⋆) (k/k_⋆)^n_T ,

with k_⋆ = 0.05 Mpc^-1 being the Planck pivot scale,  𝒫_ζ(k_⋆)≃ 2.1× 10^-9 <cit.> corresponding power spectrum of the scalar perturbation at the pivot scale, r is the tensor-to-scalar ratio and n_T represents the tensor spectral index. In single-field inflationary models, n_T ≃ -r/8. Considering the current bound on r<0.036 <cit.>, we set n_T =0 throughout our analysis, which corresponds to a scale-invariant primordial spectrum. In Eq. (<ref>) the transfer function 𝒯(a, k) connects primordial mode functions with mode functions at some later time as <cit.>

    𝒯(a, k) = 1/2(/a)^2,

where the prefactor 1/2 appears due to oscillation-averaging the tensor mode functions <cit.> and   is scale factor at the horizon crossing defined as H() = k. Different k modes might reenter the horizon at different epochs, i.e., during radiation domination or during reheating.

As is evident from Eq. (<ref>), the transfer function characterizes the expansion history between the moment of horizon crossing a = of a given mode k and some later moment a>. From Eq. (<ref>), one can see that the spectral GW energy density at present (that is, at a=a_0) reads

    (k) ≡(a_0, k) = 1/24(k/a_0  H_0)^2 𝒫_T,prim(/a_0)^2,

where H_0 is the Hubble expansion rated as measured today. Depending on the epoch of horizon re-entry for the perturbations, one can obtain the full GW spectral energy density as

    ()    ≃Ω_γ^(0) 𝒫_T,prim/24 (T_hc)/2((T_0)/(T_hc))^4/3
               ×()/(T_hc)((T_hc)/())^4/3(/)^2(n-4)/n+2   for  a_I < ≤ ,
    
            1    for ≤≤ a_ eq ,

where Ω_γ^(0)≡ρ_γ,0/ρ_c=2.47× 10^-5 h^-2 <cit.> is the fraction of the energy density of photons at the present epoch and  corresponds to the temperature at which the corresponding mode reenters the horizon. The first line holds for the modes that re-enter the horizon during reheating, while the second line applies for the modes re-entering after reheating, during the radiation-domination epoch. Note that, for ≤≤, we have assumed conservation of entropy from the moment of horizon crossing until today, implying T∝ a^-1 ^-1/3, with a = being the scale factor at the matter-radiation equivalence at T≡≃ 0.7 eV. The corresponding GW frequency f can be expressed as

    f() ≡k/2π  a_0 =    1/6√(()/10)((T_0)/())^1/3T_0 /M_P
           ×(/)^2(n-1)/n+2   for  a_I < ≤ ,
    /   for ≤≤ a_ eq ,

with T_0 ≃ 2.3 × 10^-13 GeV <cit.>. Combining Eqs. (<ref>) and (<ref>) we can rewrite the GR spectrum as a function of the frequency

    (f)    ≃Ω_γ^(0) 𝒫_T,prim/24 (T_hc)/2((T_0)/(T_hc))^4/3
       ×()/(T_hc)((T_hc)/())^4/3(f/)^n-4/n-1   for ≤ f < f_max ,
    
            1    for  f_eq≤ f ≤ ,

where f_eq≡ f() and

    f() ≡1/6√(()/10)((T_0)/())^1/3T_0 /M_P .

Finally, the GW spectrum presents a cutoff in frequency corresponding to =a_I,

    f_max≡H_I/2 π a_I/a_0
        = H_I/2 π((T_0)/())^1/3a_I/ T_0/ .

From Eq. (<ref>) one can see that modes that reenter the horizon after reheating (that is f  <) keep the same dependence with the frequency as the original primordial spectrum 𝒫_T,prim. However, if they reenter during reheating (< f <), (f) gets a boost ∝ f^n-4/n-1, which implies that the spectrum becomes blue tilted with respect to the primordial one. Interestingly, this only happens for n > 4, as for n=4 the Universe scales as free radiation during (and after) reheating.


 

In Fig. <ref>, we show the GW spectral energy density at present epoch, as a function of the frequency f for /a_I = 10^3 (left) and /a_I=10^10 (right), and three choices of the inflaton potential: n=4 (solid black), n=6 (dashed black) and n=8 (dashed-dotted black). We also used H_I = 2 × 10^13 GeV and the maximum value for the tensor-to-scalar ratio r = 0.035 <cit.>. We overlaid in yellow projected sensitivity curves from the Big Bang Observer (BBO) <cit.>, ultimate DECIGO <cit.>, LISA <cit.>, μAres <cit.>, the cosmic explorer (CE) <cit.>, the Einstein Telescope (ET) <cit.>, and from resonant cavities <cit.>, see Ref. <cit.> for a recent review. In all cases, the GW energy density is safe from the constraint on  <cit.>, which quantifies the extra radiation-like degrees of freedom around BBN. Note that for a larger /a_I (right panel) the GW spectrum receives a larger boost because in that case the reheating lasts longer, resulting in a longer stiff epoch that contributes to the blue tilt. As a consequence, the spectrum falls within the sensitivity of not only the high-frequency detectors like resonant cavity but also low-frequency detectors like BBO and DECIGO. We emphasize that although such a spectrum is not unique to resonant reheating scenario (and can be generic to any stiff background in the early Universe), any potential signal of primordial GW either in the low or in the high frequencies can not rule out the possibility of resonant production of radiation bath.




§ CONCLUSIONS
 

In this study, we investigate a reheating scenario where the transfer of energy from the inflaton ϕ to the Standard Model (SM) particles proceeds through an s-channel mediated scattering. We consider the inflaton to oscillate around a generic monomial potential ∝ϕ^n, with the inflaton mass m_ϕ being field dependent and therefore varying over time if n>2. If the mediator mass m_s is of the order of the inflaton mass m_ϕ, a resonance can develop during reheating once the inflaton mass is m_ϕ(ϕ) ≃ m_s/2.  We validate these characteristics by solving a system of coupled Boltzmann equations for the inflaton and SM energy densities (Eqs. (<ref>) and (<ref>)) and observe a characteristic bump occurring both in the radiation and in the temperature profiles, as depicted in Fig. <ref>. Furthermore, we derive several analytical approximations, detailed in Eq. (<ref>), which nicely fits the numerical results. Our investigation reveals that resonance reheating exhibits more nuanced temperature scaling compared to conventional scenarios presented in the existing literature, which is typically characterized by a featureless power law, as given in Eq. (<ref>). In scenarios where the mediator is always heavier than the inflaton, precluding the development of resonance, our model successfully reproduces the results of the existing literature. Alternatively, when the mediator is always lighter than the inflaton, the SM temperature rises throughout the reheating phase, with the reheating temperature serving as the maximum temperature achieved; see Fig. <ref>.

Some implications for particle phenomenology and cosmology during resonant reheating were studied. First, the mediator could play a pivotal role in producing not only the SM states but also particles beyond the SM, such as dark matter (DM). If the interaction rates between the DM and the visible sector are very suppressed, DM could be produced during reheating through the annihilation of a pair of SM or inflaton particles. Similarly to radiation, DM production also has a very rich freeze-in phenomenology, as illustrated in Figs. <ref>, <ref>, and <ref>.

Second, we also proposed the prospect of finding a signature of this model via the primordial gravitational wave (GW) spectrum. For n>4, the inflationary GW spectrum experiences a significant blue tilt within the frequency range corresponding to the modes that cross the horizon during the stiff period. The resulting spectrum could be well within the reach of several next-generation low- and high-frequency GW detectors, as shown in Fig. <ref>.



NB received funding from the Spanish FEDER / MCIU-AEI under the grant FPA2017-84543-P. YX acknowledges the support from the Cluster of Excellence “Precision Physics, Fundamental Interactions, and Structure of Matter” (PRISMA^+ EXC 2118/1) funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) within the German Excellence Strategy (Project No. 390831469).







§ ANNIHILATION TO MEDIATORS
 

In our setup, the annihilation ϕϕ→ s s through a t- and u-channel exchange of a ϕ is also kinematically open. We assume that s particles have a sufficiently large decay width and rapidly decay into SM radiation once they are produced.
The total production rate considering the t and u channels is

    Γ^ϕϕ→ s s≃ρ_ϕ/m_ϕ 1/16π  m_ϕ^2(μ_ϕ^2/m_s^2 - 2 m_ϕ^2)^2 √(1-m_s^2/m_ϕ^2) .

The solution for the radiation energy density is

    ρ_R(a) ≃9/64π n/8n-17( a/a_I)^6 (2 n-7)/2+n[1 -( a_I/a)^2 (8 n-17)/2+n] H_I^3  M_P^4 μ_ϕ^4/M_I^7

in the limit m_s ≪ m_ϕ. Note that if μ_h = μ_ϕ, Eq. (<ref>) and the first line of Eq. (<ref>) are the same (up to a factor 8), making this channel potentially important. However, if one demands an extended reheating period, a long-lived inflaton is desirable, and hence a suppressed μ_ϕ. Under these circumstances, the inflaton annihilation into mediators is suppressed. We have checked that in all the benchmarks used in this work, this channel is completely subdominant.


JHEP


