Testing Axionic Dark Matter during Gravitational Reheating

Assuming axions are potential dark matter (DM) candidate that make up all of the DM abundance, we discuss production of axions via (i) standard misalignment mechanism during the period of reheating and (ii) graviton-mediated 2-to-2 scattering of the inflaton and bath particles, where the inflaton $\phi$ oscillates in a monomial potential $V(\phi)\propto\phi^k$ with a general equation of state. Considering reheating takes place purely gravitationally, mediated by massless gravitons, we explore the viable region of the parameter space that agrees with the observed DM relic abundance, satisfying bounds from big bang nucleosynthesis (BBN) and cosmic microwave background radiation (CMB). We also discuss complementarity between dedicated axion search experiments and futuristic gravitational wave search facilities in probing the viable parameter space.


Introduction
The QCD axions [1,2], pseudo-Nambu-Goldstone boson of the Peccei-Quinn (PQ) solution to the strong CP problem [3][4][5] and the axion like particles (ALP), that could also arise from the spontaneous breaking of a global U (1) symmetry, stand out as especially well motivated candidates for cold dark matter (DM) [6][7][8][9].QCD axions and ALPs arise in various extensions of the Standard Model (SM) through spontaneous symmetry breaking (SSB) or from string compactification [10,11], with a potential of getting discovered in the next decade (see, for example, Ref. [12]).
In the standard scenario, axions (by "axions" we collectively refer to QCD axions and ALPs) can be produced in the early universe via the "misalignment mechanism," in which the QCD axion or the ALP field modeled as the classical scalar field due to its bosonic nature and high occupation numbers, has a non-zero initial field value and non-zero potential energy, leading to oscillations of the field 1 .For an O(1) initial misalignment angle θ i , the allowed mass window for QCD axion that produces the observed DM abundance turns out to be 10 −6 ≲ m a ≲ 10 −5 eV, when the oscillation begins during radiation dominated Universe.For ALPs, on the other hand, the relic abundance depends on three parameters: the decay constant f a , ALP mass m a , and θ i , leading to strong bounds on the viable parameter space for θ i ∼ O (1).It has been shown that deviations from standard cosmological histories can significantly broaden the parameter space for both axions and ALPs [16][17][18][19].Very recently, Ref. [20] has pointed out that such a conclusion also holds true if misalignment happens during the epoch of reheating if the inflaton2 ϕ oscillates in a monomial potential of the form V (ϕ) ∝ ϕ k , that provides a general equation of state (EoS) 0 ≤ w = (k − 2)/(k + 2) ≲ 1 for the inflaton.
On the other hand, the irrefutable coupling that one can imagine between matter particles (irrespective of dark and visible sector), is of gravitational origin, mediated by massless graviton [24,25].Such an interaction inevitably exists between the inflationary sector and axions.A graviton portal between the inflaton and the SM sector can also produce the radiation bath, completing the process of reheating even in the absence of direct coupling between the inflaton and the SM fields.This is dubbed as the "gravitational reheating" scenario, that has been shown to be efficient for k ≳ 9(4) [26,27], considering (non-)minimal contribution to radiation.Now, k > 4 implies a stiff EoS for the inflaton, that results in a significantly blue-tilted primordial gravitational wave (GW) spectrum, originating from the tensor perturbations during inflation [28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45].Although a stiff period during reheating makes the GW signal potentially observable by GW detectors, but the very same enhancement also results in overproduction of GW energy density, that violates standard BBN and CMB predictions.
Motivated by these, in this paper we have discussed a scenario where axion misalignment takes place in an inflaton dominated background during reheating, supplemented by an attractor potential for the inflaton.The reheating takes place purely gravitationally via scattering of the inflaton condensate into a pair of Higgs bosons, mediated by massless gravitons, in contrast to [20], where reheating occurs via perturbative decay of the inflaton condensate into bosonic and fermionic final states.Interestingly, the oscillation temperature T osc (and hence the relic abundance) becomes sensitive to the choice of k (along with m a and θ i ) in case of minimal gravitational reheating, and to both {k, T rh } in case where reheating takes place via a non-minimal coupling between the Higgs and gravity (Ricci scalar).We find that the standard misalignment mechanism during reheating opens up more parameter space for axions, compared to misalignment during radiation-domination, making up all of the DM abundance.The overproduction of the primordial GW energy density around the time of BBN rules out the minimal reheating scenario, while the non-minimal reheating remains within the sensitivity of future GW and axion search experiments.We thus find a complementarity between axion search and GW experimental facilities in constraining the allowed parameter space.
The paper is organized as follows.In Sec. 2 we discuss the details of gravitational reheating and the generation of primordial gravitational wave.The mechanism of standard misalignment during reheating producing viable parameter space is discussed in Sec. 3. Pure gravitational production of axions, mediated by graviton, is elaborated in Sec. 4. In Sec. 5 we discuss the discovery potential of this framework.Finally, we summarize our findings in Sec. 6.

Post-inflationary evolution of the Universe
The interaction between all matter fields and the gravitational field can be found by expanding the metric around Minkowski space-time η µν as g µν ≃ η µν + 2hµν M P , where h µν represents the canonically normalized quanta of the graviton.Consequently, one obtains possible gravitational interactions from the Lagrangian [24,25] where "SM" denotes the SM fields, ϕ is the inflaton while X represents other beyond the SM (BSM) fields, which in our case is an axion.Here M P ≃ 2.45 × 10 18 GeV is the reduced Planck mass.The form of stress-energy tensor T µν is dictated by the spin of the fields.For a generic spin-0 scalar S, such as the Higgs bosons, the inflaton or an axion3 , it can be expressed as where V (S) represents the potential of the respective scalar.
For the inflationary potential, we adopt the following α-attractor T-model [46] that provides a plateau region in the large field limit, leading to quasi-de Sitter expansion consistent with observation 3) The overall scale of the potential, parameterized by the coupling λ, can be expressed in terms of the amplitude of the scalar perturbation power spectrum , where N ⋆ is the number of e-folds measured from the end of inflation to the time when the pivot scale k ⋆ = 0.05 Mpc −1 exits the horizon.Here onward we will also fix α = 1/6.The ending of inflation is marked when ä = 0, at which the inflaton field has a magnitude Furthermore, one can compute the effective mass of the inflaton field by taking the second derivative of the scalar potential, i.e., m 2 ϕ = ∂ 2 V ∂ϕ 2 that at the end of inflation turns out to be m ϕ (a = a end ) ≃ 10 13 GeV for above choice of the CMB observable, with a very mild dependence on the exponent k.With this set-up we will now move on to the computation of gravitational reheating temperature.

Gravitational reheating
In order to track the evolution of the inflaton (ρ ϕ ) and the radiation (ρ R ) energy densities during reheating, we solve the following set of coupled Boltzmann equations (BEQ) [48] ) together with Here, w ϕ ≡ [49] is the general equation of state parameter for the inflaton.Since during most part of the reheating, the total energy density is dominated by the inflaton, the expansion rate corresponding to the term 3 H (1 + w ϕ ) ρ ϕ dominates over the reaction rate (1 + w ϕ ) Γ ϕ ρ ϕ in Eq. (2.5).As a consequence, one can solve Eq. (2.5) analytically by neglecting the right-hand side and obtain with corresponding Hubble rate (2.9) Solution of Eq. (2.6) requires the information of the reaction rates of inflaton.Since we are interested in the purely gravitational reheating, i.e., no direct coupling between the inflaton and the SM states is present, the production rate of radiation in that case can be evaluated as [50][51][52] ( by considering the graviton propagator for momentum q as Here, we have considered the interaction of inflaton with only the SM Higgs field, neglecting its mass.consequently, N h = 4 is considered as the number of internal degrees of freedom for one complex Higgs doublet.While evaluating the interaction rate, the inflaton is treated as a time-dependent external and classical background field, which we parametrize as where ϕ 0 (t) is the time-dependent amplitude that includes the effects of redshift and T (t) describes the periodicity of the oscillation.We also expand the potential energy in terms of the Fourier modes as [49,50,[52][53][54][55]] where the frequency of oscillations of ϕ field reads [49]  where g * is the relativistic degrees of freedoms present in the thermal bath, and we assume instantaneous thermalization.From Eq. (2.5) and (2.6), we note that the thermalization process of the SM particles produced from the inflaton scattering helps the Universe to attain a maximum temperature T max right at the end of inflation.Subsequently the temperature falls to T rh , where equality between ρ ϕ and ρ R is achieved.As a result, reheating temperature can be evaluated as (2.17) One can further note that, for a end ≪ a ≪ a rh , the temperature evolves as [cf.Eq. (2.15)] Purely gravitational scattering process discussed above may not always lead to efficient reheating of our Universe.For example, with k = 2, the radiation energy density produced by inflaton scattering never comes to dominate the energy density and can not lead to reheating.For k > 4 reheating from gravitational scattering is possible.However it is very inefficient,i.e., with k = 6, from Eq. (2.17), we obtain T rh ≪ 1 eV < T BBN .One actually needs to go beyond k = 9, for which the gravitational reheating temperature can be found to be T rh ≃ 2 MeV.On top of that, as we shall explain, the minimal gravitational reheating scenario is completely ruled out from the overproduction of inflationary GW energy density around the time of BBN.This motivates us to go beyond the minimal scenario and introduce a non-minimal coupling.
The non-minimal coupling of Higgs bosons to gravity (via an interaction of the form ξ h H † H R, R being the Ricci scalar [51]) provides an additional channel to reheat with the rate [51,52] ( ) where numerical estimates of the co-efficient α for different values of k are reported in Table .1 of Ref. [27].The non-minimal reheating temperature in this case can be obtained as while the maximum temperature in this case is determined by  which is O(10 12 ) GeV, with mild dependence on k [27,50].With non-minimal contribution taken into account, we note, reheating can be completed before the onset of BBN for k > 4, by tuning the non-minimal coupling properly.We show reheating temperature as a function of k for different choices of the nonminimal coupling in the left panel of Fig. 1, where the black solid curve corresponds to the minimal gravitational reheating scenario (ξ h = 0).The shaded region, corresponding to minimal gravitational reheating, is ruled out from GW overproduction as we shall elaborate in subsection.2.2.The right panel shows corresponding maximum temperature T max for each ξ h .For higher ξ h , as expected, the ratio T max /T rh becomes smaller since T rh can be larger in those cases.Thus, for a given k, thanks to the non-minimal coupling, it is possible to have T rh much larger than that compared to minimal gravitational reheating scenario.Hereafter, for the non-minimal case, we will consider {k, T rh } as free-parameters.

Inflationary gravitational wave
In this section we briefly describe the background theory for computing the spectral energy density of primordial gravitational wave (GW), emerging from the tensor fluctuations during inflation (for a review, see, for example, Ref. [36]).GWs are transverse (∂ i h ij = 0) and traceless (h ii = 0) metric perturbations Their energy density spectrum per momentum mode (at sub-horizon scales) is defined as [36,56] where ∆ 2 h (t, k) is the tensor power spectrum at arbitrary times, defined as with ⟨...⟩ denoting an average over a statistical ensemble.One can factorize the tensor power spectrum as [57,58] with the transfer function [36,56,58] (1/2 appears due to oscillation-averaging the tensor mode functions) and ∆ 2 h,inf (k) representing the primordial tensor spectrum from inflation [36] where n t is the spectral tilt, k p denotes a pivot scale of the order Hubble rate at the time of CMB decoupling and H inf represents the Hubble rate when the mode k p exited the Hubble radius during inflation.For simplicity, we shall assume an exact scale-invariant inflationary spectrum or, in other words, n t = 0.The GW background modes produced from such tensor fluctuations can cross outside the horizon during inflation when k < a H holds and can be considered as classical modes.Subsequently, at a later time after inflation, these modes reenter the horizon (k > a H) and form a stochastic background.Now, any extra radiation component, in addition to those of the SM, can be quantified in terms of the ∆N eff .Since GW energy density scales the same way as that of free radiation, it is possible to put an upper bound on ρ GW as an extra radiation component at the time of BBN and/or CMB decoupling.This can be done by computing the total radiation energy density in the late Universe as where ρ γ , ρ ν , and ρ GW correspond to the photon, SM neutrino, and GW energy densities, respectively, with T ν /T γ = (4/11) 1/3 .Within the SM, taking the non-instantaneous neutrino decoupling into account, one finds N SM eff = 3.044 [59][60][61][62][63][64][65][66][67], while the presence of GW results in a modification (2.28) The above relation can be utilized to put a constraint on the GW energy density red-shifted to today via [32,36,68] 4 11 4 3 which leads to Ω (0) GW h 2 ≃ 5.62 × 10 −6 ∆N eff where, Ω γ h 2 ≃ 2.47 × 10 −5 is the relic photon abundance at the present epoch while f = k/(2πa 0 ) is the present day frequency of the physical wave number k.Here f max corresponds to maximum frequency that re-enter the horizon right after the end of inflation when k max = a end H end , and is given by while f BBN ≃ 2 × 10 −11 Hz, corresponds to the mode k BBN = a BBN H BBN at the time of BBN.Using the present CMB measurement from Planck legacy data [47], we find Ω (0) GW h 2 ≲ 2×10 −6 , with ∆N eff ≃ 0.34 4 .Now, the ratio of the gravitational wave (GW) energy density to that of the radiation bath is given by ρ [27,56], where k hc = a hc H hc is the momentum mode, calculated at the moment it re-enters the horizon ("hc" denotes horizon-crossing).If horizon crossing occurs during radiation domination k 2 hc ∝ ρ R , then the GW spectrum becomes scale invariant.On the other hand, if horizon crossing occurs during the inflaton-dominated era, the GW energy density is enhanced by a factor of ρ ϕ /ρ R evaluated at T hc .As a result, the largest enhancement occurs for the mode that reenters the horizon right after inflation at T max .For minimal gravitational reheating (ξ h = 0), such enhancement turns out to be ρ end /ρ R ≃ (4 − 6) × 10 13 for k ∈ [6,20], resulting in Ω GW h 2 ≃ (8 − 10) × 10 −6 .This is in clear conflict with the present bound from Planck on ∆N eff , as discussed above.The constraint can be relaxed by increasing the value of T max , since in that case the the energy density of GW relative to that of radiation at T max becomes smaller, compared to the minimal gravitational reheating scenario.The region labelled as 'inconsistent reheating' in Fig. 1 thus corresponds to ξ h < 0, which is forbidden from ∆N eff bound due to Planck on excessive GW production (shown in blue).
Before closing this section, it is important to highlight for k > 4, the inflaton energy density redshifts faster than radiation, and the equation of state rapidly evolves from some stiffer fluid w ϕ > 1/3 to w ϕ = 1/3, and the universe is dominated by massless inflaton particles (not a condensate anymore), with energy density redshifting as a −4 [74,75].These free particles are produced as a consequence of the self-interaction of the inflaton, a process known as fragmentation.A detailed study of such processes, taking into account the effect of parametric and tachyonic resonance, requires dedicated lattice simulations and has been studied, for example, in Refs.[74][75][76][77][78][79][80][81].This, however, is beyond the scope of the present draft.
The Lagrangian density for axion field reads [82] which leads to equation of motion of the zero modes as where θ(t) ≡ a(t)/f a and f a denotes the decay constant.Here, the temperature dependent mass of the axion is denoted by ma (T ), which for the QCD axion, is shown to be dependent on the topological susceptibility of QCD χ(T ) via [83] ma (T ) = χ(T )/f a . (3. 3) The lattice QCD simulations have provided an estimates of the zero-temperature value of χ(T ), which is given by χ 0 ≡ χ(0) ≃ 0.0245 fm −4 , in the symmetric isospin case.Subsequently, the form of the ma (T ) is found to be with m a representing the axion mass at zero temperature and is given by [84] m a ≃ 5.7 × 10 −6 10 12 GeV f a eV .
Following Eq. (3.2), axion begins to oscillate at the temperature T = T osc defined by 3 H(T osc ) ≡ ma (T osc ) [85].Assuming a radiation dominated Universe, the corresponding oscillation temperature can be evaluated as T osc ≥ T qcd , (3.6) where the Hubble expansion rate takes the form H(T ) = ρ R (T )/(3 M 2 P ).Below T osc , the axion can be considered as a non-relativistic particle.Considering the conservation of the axion number density and assuming conservation of SM entropy, the energy density for such non-relativistic axions ρ a at present is given by with T 0 representing the temperature today and the SM entropy density is defined as where g ⋆s (T ) denotes the corresponding number of relativistic degrees [86].The WKB approximation leads to ρ a (T osc ) , where θ i is the initial misalignment angle [84,87].Eventually, using Eq.(3.7), the axion abundance can be found to be for m a ≳ m qcd a , (3.9) GeV/cm 3 is the critical energy density and s(T 0 ) ≃ 2.69 × 10 3 cm −3 is the entropy density at present [88].
Next, we take up the scenario when the oscillation of the axion starts during the reheating period i.e., T osc > T rh .In this case, the axion energy density at present epoch reads where the last factor takes into account the dilution of the axion energy density due to the entropy injection as a result of energy transfer from the inflaton to the bath.This dilution factor can be determined as where we have used Eq.(2.15).Note that since the non trivial expansion of the Universe is the sole reason for the dilution of radiation, the factor accounting for entropy injection is simply ratio of DOFs, unlike the case in [20].Now, one can further consider two sub-cases depending on the hierarchy between T qcd and T osc .Below we discuss them one by one.

Scenario-I:
We first consider the case T qcd < T osc , for which the axion mass shows a temperature dependence ma = m a (T /T qcd ) −4 .One can then obtain the expression for oscillation temperature as (3.12) Note that, this scenario can be further classified as: T rh < T qcd < T osc and T qcd < T rh < T osc .
The first inequality provides while from the second inequality we get Scenario-II: In the second case, T qcd > T osc , the axion mass remains constant.In this case we can again derive an analytical expression for the oscillation temperature as for ma = m a .Using the inequality T rh < T osc one finds while T osc < T qcd further provides In Fig. 2 we show QCD axion oscillation temperature as a function of the scale f a .In the left panel we consider standard radiation dominated Universe, while in the right panel we consider an inflaton-dominated background.In either cases we see a bend around T osc ≃ T qcd due to the change in the temperature dependence of the axion mass.In the right panel we note, depending on the choice of k (that decides steepness of the inflaton potential) the oscillation temperature changes, following Eq.(3.12) and (3.15).The analytically derived expressions for T osc are denoted with the red dotted lines for one of the values of k, showing the agreement between analytical and numerical results.The gray shaded region in both plots is forbidden as it surpasses the Planck scale.While k > 9 is necessary to have successful BBN in case of minimal gravitational reheating, as k grows, T rh starts increasing till it reaches k ≃ 12, above which T rh > T osc , rendering the oscillation to take place during radiation domination.So, in case of minimal gravitational reheating our parameter space is confined within 9 < k ≲ 12 (this is more prominent from the right panel Fig. 3 which will be discussed in a moment),   for which misalignment happens during reheating.In case of gravitational reheating via non-minimal coupling, one can, however, go to higher k-values.
Utilizing Eq. (3.10), we can obtain the relic abundance for T rh < T osc as for m a ≳ m qcd a , where we have fixed k = 10, which also fixes T rh = 0.05 GeV in case of minimal reheating scenario.In the left panel of Fig. 3 we show parameter space in [k, f a ] plane that reproduces all of the Planck observed relic abundance (shown via black curve), considering misalignment happens during gravitational reheating.As we see, a larger initial misalignment angle shifts the parameter space to lower f a (or heavier axion mass) since the relic abundance varies inversely with m a for a fixed θ i , as seen from Eq. (3.18).We show contours corresponding to T osc = T rh and T osc = T qcd via blue and orange broken curves respectively.As a result, parameter space of our interest (T osc > T rh ) lies below the blue dashed contour.The shaded regions are disallowed from T rh < T BBN (or, equivalently, k ≲ 9) and f a > M P .For each k, we denote corresponding T rh (fixed for ξ h = 0) along the right vertical axis.The gray vertical dashed lines show required f a that can produce the observed relic abundance if the oscillation happens during radiation dominated background.For the same axion mass, gravitational reheating thus opens up larger parameter space, satisfying DM abundance.More precisely, for θ i = π/ √ 3, misalignment during RD produces right abundance for f a ≃ 10 12 GeV (equivalently, m a ≃ 4.2 × 10 −6 eV), whereas for minimal gravitational reheating this window becomes 10 10 ≲ f a ≲ 10 12 GeV (or, 10 −4 ≲ m a ≲ 10 −6 eV).In the right panel of Fig. 3, the red and black curves correspond to contours producing observed DM abundance for different choices of θ i , where we fix f a = {10 11 , 10 12 } GeV.Here we consider non-minimal contribution to the radiation, and show the resulting parameter space in [T rh − k] plane.The relic density satisfying contours correspond to θ i = {0.1,1.0}, shown via solid and dashed pattern respectively.As mentioned in Sec. 2, for each k, in case of non-minimal gravitational reheating, one now have the freedom choose appropriate ξ h that can provide the desired T rh (larger than that corresponding to minimal reheating scenario).The region below the blue dashed line, parallel to the horizontal axis, indicates T osc > T rh .Here we see, a larger θ i requires larger T rh , as one can infer from Eq. (3.18).Also, for a fixed θ i , lower f a (higher m a ) requires lower T rh in order to produce right abundance, again following Eq.(3.18).As one can notice, the ∆N eff bound plays a very crucial role in constraining the viable parameter space, ruling out θ i ≲ 0.1 for both f a 's.
Finally, we turn our attention towards ALPs, for which we consider the mass m a to be time independent.For oscillations during reheating T rh < T osc we find the oscillation temperature to be same as Eq.(3.15).This puts a bound on the reheating temperature, which in turn can be translated into a lower bound on the ALP mass via The ALP energy density at present day reads with ρ a (T osc ) ≃ (1/2) m 2 a f 2 a θ 2 i , this turns out to be In case of ALPs, f a and m a can vary independently, that makes following set of parameters free: {m a , f a , θ i , k} for minimal and {m a , f a , θ i , k, T rh } for non-minimal reheating scenarios.From Eq. (3.21) we find, in order to satisfy the observed DM abundance, f 2 a ∝ m 2−k k a , implying, a larger f a requires lighter ALP for a fixed k > 2 (and θ i ).This is what we see in the top panel of Fig. 4, where the black thick contours satisfy the observed relic abundance for different choices of k with fixed θ i 's.Corresponding to each k-values, we also show contours satisfying T osc = T rh in blue.Since in the minimal reheating scenario fixing k fixes the reheating temperature, hence we see simple straight line contours parallel to the f a axes.The change in slope for each contour occurs at T osc = T rh , indicated by the blue vertical straight lines.To the left of each vertical line oscillation occurs during RD, while to the right during reheating.The parameter space of our interest therefore lies to the right of each blue vertical line.In the bottom panel of Fig. 4 we show contours of correct relic density for different choices of f a , while fixing m a for a given θ i = 1, now considering non-minimal gravitational reheating.Again, the allowed parameter space lies below the blue dashed curve, for which T osc = T rh .As we see, for m a = 1 eV (bottom left panel), f a ≲ 10 9 GeV is completely ruled out from BBN bound on T rh , together with ∆N eff constraint on GW energy density, irrespective of the choice of k.This can be relaxed by considering a heavier m a = 1 GeV, as seen from the bottom right panel.The reason being, the relic density iso-contours satisfy f 2 a m a ∝ T rh for a fixed θ i , following Eq.(3.22).Therefore, for the same choice of f a and θ i , heavier ALPs require higher T rh to produce the right abundance.

Purely gravitational axions
To this end we have discussed axion production via standard misalignment mechanism in an inflaton-dominated background.However, due to the inevitable gravitational interactions, through Eq. (2.1), axions are also produced via gravity-mediated 2-to-2 process, which we refer to as purely gravitational production of axions.Such purely gravitational production can happen from (i) thermal bath and (ii) scattering of the inflaton condensate.Below we discuss the details of each such processes and obtain the resulting parameter space.

Thermal scattering
In case of production from thermal bath due to scattering of the SM particles mediated by gravity, the rate of production reads [51,89] R T a = 3997 π 3 41472000 Since we are interested in axion production during reheating, it is instructive to consider the comoving axion number N a = n a × a 3 as the entropy is not conserved during reheating.To track the comoving axion number, we solve the corresponding Boltzmann equation One can then find the axion number density n a at the end of reheating as [27] n T a (a rh ) = where we have made use of Eq. (2.15).Since the gravitational reheating temperature is generally quite low [cf.Sec.2.1], considering a ≫ 1 we obtain where a = a rh /a end ≡ (H end /H rh ) k+2/(3 k) [cf.Eq. (2.9)] and α ⋆ = π 2 /30 g ⋆ (T rh ).

Inflaton scattering
For gravitational production from inflaton condensate scattering, viz., ϕϕ → aa, the particle production rate per unit volume and unit time for arbitrary k reads [27,50,51] where The factor of 2 accounts for pair production.Here E n = n ω is the energy of the n-th inflaton oscillation mode.Note that, such a production is possible only during reheating and not after, unlike the production from radiation bath that can take place both during and post-reheating.The evolution of comoving axion number produced from inflaton scattering mediated by graviton reads Again, the number density of axions at the end of reheating can be computed as Taking both contributions into account, we determine the gravitationally produced axion relic abundance as where n a (T rh ) = j={ϕ, T } n j a (T rh ).It is interesting to note, the ratio of axion number density produced from thermal and from inflaton scattering reads since ρ end ≫ ρ rh .This implies that purely gravitational production of axions is always dominated by the scattering of the inflaton zero modes.Now, axions produced from inflaton scattering during reheating are relativistic as the average momenta they carry is of the order of the inflaton mass, which is ∼ 10 13 GeV at the beginning of reheating.Such a (semi-)relativistic population of axions at CMB decoupling behave as dark radiation and contribute to ∆N eff .Following the analysis in Refs.[90,91],the average thermal velocity of axions today ⟨v a,0 ⟩ is related to ∆N eff and relic abundance ⟨v a,0 ⟩ ≃ 5.62 × 10 −6 × ∆N eff /Ω a h 2 , where we assume the axions are relativistic/semi-relativistic at photon decoupling, a legitimate assumption for masses around the eV scale.Such a population of (semi-) relativistic axion is distinguishable from CDM if ⟨v a,0 ⟩ ≳ 1 km/s.Since the axion number density: n a ∝ T 2 rh /M 3 P , in the limit of large k [cf.Eq.(4.8)] and m a ≪ M P , for purely gravitational QCD axions Ω a h 2 ≪ 10 −10 .This is however not true in case of gravitationally produced ALPs, since their masses can be large enough, independent of the choice of f a .In Fig. 5 we show parameter space where the co-existence of misalignment and gravitationally produced population of ALPs can happen for a benchmark value of mass m a = 0.01 MeV for scale f a = 10 10 GeV.For these choice of parameters, ALPs produced via misalignment during reheating can produce the observed relic abundance for θ i = [0.5, π/ √ 3].We therefore show underabundant contours for ALPs produced purely from gravity.We thus see, two populations of QCD axions: one cold, via the standard misalignment mechanism during reheating, and the other originating purely gravitationally from scattering of inflaton can actually co-exist.

Experimental limit on the parameter space
In this section, we wish to explore the possibilities of probing the parameter space satisfying the observed relic abundance, for axions produced via standard misalignment during gravitational reheating, at the present and proposed axion search facilities.In order to do that, we examine the interaction of axions with two photons, a highly utilized channel for detecting signatures in both observational studies and experimental investigations.The Lagrangian for such an interaction has the following form [84,92] where the coupling strength g aγ is model dependent and is related to the PQ scale f a as [93] |g where z = m u /m d is the ratio of quark masses, and E and N are the electromagnetic and colour anomalies associated with the axion anomaly.For KSVZ models E/N = 0 [94,95], whereas for DFSZ models E/N = 8/3 [96,97].For ALPs, on the other hand, one can relate the ALP-photon coupling with the PQ scale f a via [98] where |c γ | can vary between 1 and 10.We map our viable parameter space on a bi-dimensional plane of [g aγ , m a ], on which we further project limits from different experimental facilities.In the left panel of Fig. 6 we show reach of existing and future experiments in probing the relic density allowed parameter space for QCD axions produced via misalignment during non-minimal gravitational reheating.In obtaining the parameter space, we have scanned over a range of k ∈ [6,20] and T rh ∈ [10 −4 , 10 14 ] GeV along with θ i ∈ [0.5, π/ √ 3].The vertical green shaded band corresponds to right relic density for QCD axions produced via standard misalignment mechanism during radiation domination.As mentioned before, oscillation during gravitational reheating broadens the window of QCD axion mass (equivalently, the PQ breaking scale) satisfying the observed relic abundance.This is what is reflected here as well.The thick red slanted band is the allowed parameter space that satisfies T osc > T rh > T BBN .However, when the bound Ω (0) GW h 2 ≳ 2 × 10 −6 is imposed on top of that, the parameter space is confined within the blue band.As we see, a part of the viable parameter space is already constrained from existing limits from Haloscope experiments like ADMX [99], CAPP [100], ORGAN [101], as shown by the gray shaded region, while the darker red shaded region is the future projection from Haloscope experiments 5 .Proposed experiments like DMRadio [102] (shown in orange) or future projection from broadband axion-search experiment ABRACADABRA [103] (shown in blue) are capable of constraining the parameter space further.
Finally, in the right panel we show relic density allowed parameter space for ALP for two representative values of k = {8, 12}, shown via cyan and orange bands, respectively.In this case we consider the entire abundance of ALP is produced through standard misalignment during gravitational reheating.For k = 8 scenario we choose T rh = T BBN , whereas for k = 12, T rh ≤ 1 GeV is disallowed from ∆N eff bound on GW energy density [cf.Fig. 1].We therefore choose T rh = 10 GeV corresponding to k = 12 scenario.In all cases we have considered |c γ | = 1, and scanned over the initial misalignment angle θ i ∈ [0.5, π/ √ 3].The first noticeable feature here is that, similar to the case of QCD axions, more parameter space becomes available for ALPs when the misalignment occurs during gravitational reheating.However, since m a and f a now can vary independently, hence the resulting parameter space is broader, compared to the QCD axion case.The important point here is to note that for smaller k, it is possible to produce the observed abundance with lighter ALPs, a feature we already noticed in the top panel of Fig. 4. As a result, here we see, with k = 12, it is possible to obtain the right relic density for m a ≳ 10 −15 eV, while for masses below this, oscillation happens during RD for k = 12 and overlaps with the blue dashed curves labelled as "Standard RD".Another point to note here is that, a lower T rh requires larger g aγ to produce the right abundance for a given mass.This is expected because from Eq. (3.22) we can already see that the relic density goes as Ω a h 2 ∝ 1/ g aγ T 2 rh , hence the k = 8 band (with T rh = T BBN ) lies above k = 12 band (with T rh = 10 GeV).Regarding the experimental reach, a part of the parameter space for k = 8 is already ruled out from present Haloscope experiment, while rest of the parameter space remains well within the reach of future sensitivities from DMRadio, ABRACADABRA and Haloscope experiments.Larger k, on the other hand, is still outside the reach of futuristic experimental facilities as in that case even smaller g aγ is required.The precise message here is crucial: any potential discovery of axions in this mass range through futuristic experiments shall not only imply a signature of new physics beyond the SM, but also hint towards gravitational reheating and therefore inflationary paradigm, that may have complementary validation from the detection of primordial gravitational waves proposed GW detection facilities.

Conclusions
In the simplest scenario, a lone scalar field drives inflation, and the interaction between that scalar field (namely, the inflaton) with the Standard Model (SM) fields is crucial for successful reheating at the end of inflation.However, even without an explicit inflaton-visible sector coupling, gravity-mediated processes can efficiently heat up the Universe after inflation.This gravitational reheating emerges as a minimal and inevitable mechanism for obtaining our current Universe.In addition to the SM particles forming the radiation bath, fields beyond the SM can also be produced through a similar graviton-exchange process, highlighting the democratic nature of gravitational interaction.
With this underlying motivation, in this work we have discussed a scenario where axions, that arise as an elegant solution to the strong-CP problem and can serve as a viable cold dark matter (DM) candidate, are produced via standard misalignment mechanism during the epoch of reheating.We consider the production of radiation bath purely gravitationally, i.e., from the scattering of inflaton condensate to Higgs final state, mediated by massless graviton.We find, misalignment during gravitational reheating offers a larger window for axion mass, for a natural choice of initial misalignemnt angle θ i ∼ O(1), compared to the scenario where misalignment takes place in a radiation-dominated Universe.As the inflaton ϕ oscillates in a general monomial potential ϕ k , its equation of state mimics that of a stiffer-than-radiation fluid for k > 4, a value that is anyway required to reheat the Universe prior to BBN via purely gravitational coupling.Such a stiff background equation of state results in a bluetilted primordial gravitational waves (GW), having inflationary origin, that rules out the minimal gravitational reheating scenario due to BBN bound on GW energy density, encoded in ∆N eff .For non-minimal reheating case, however, this puts a bound on the parameter space satisfying relic abundance, constraining typically smaller θ i 's (see Fig. 3 and 4).Due to irreducible gravitational interaction, apart from standard misalignment, axions are also produced via gravity-mediated scattering of the bath particles and inflaton.Such axions are semi-relativistic in nature, and can be distinguishable from cold axions produced via misalignment.However, these two populations of axions, namely, (semi-) relativistic and cold, can effectively co-exist (see Fig. 5).Once again, overproduction of primordial GW energy density plays an important role in constraining the resulting parameter space.
We finally discuss the discovery potential of our set-up at present and future axion search frontiers, and find out that Haloscope experiments are quite capable of ruling out and/or constraining the parameter space for both QCD axion and ALPs (see Fig. 6).A large part of the parameter space that is within the reach of (future) Haloscope experiments, can also be probed at several GW detectors, or ruled out from GW overproduction (∆N eff bound).In conclusion, our framework not only provides a complementary avenue for axion searches, but any potential discovery of axions at any of these experiments can also validate a non-standard cosmological epoch during reheating, prior to the onset of BBN.

Figure 1 .
Figure1.Left: Gravitational reheating temperature as function of k, where each curve corresponds to a fixed non-minimal coupling ξ h , as indicated in the plot.The shaded regions are forbidden from ∆N eff bound due to Planck from excessive production of primordial GW as discussed in subsection.2.2, ruling out the minimal gravitational reheating scenario ξ h = 0 (solid black curve).Right: Ratio of T rh to T max , as a function of k for the same choice of non-minimal couplings as in the left panel.

Figure 2 .
Figure 2. Left: QCD axion oscillation temperature as a function of axion mass in radiation dominated background.Right: Same as left, but in an inflaton dominated background for different choices of k (≡ T rh ), considering minimal gravitational reheating.The red dotted lines correspond to analytical solutions [cf.Eq. (3.12) and (3.15)].The brown horizontal lines correspond to T osc = T rh for k = {9, 10} from bottom to top.For each case the arrowheads show the region of parameter space where oscillation happens during reheating.

3 T
o sc .> T rh T o sc .< T rh T o sc .= T rh T osc .

Figure 3 .
Figure 3. Left: The black curves satisfy the observed relic abundance for θ i = 0.1 (thick) and θ i = π/ √ 3 (dashed), considering minimal gravitational reheating.The gray shaded regions are disallowed from BBN bound on T rh (k ≥ 9) and super-Planckian scale f a > M P .Along each vertical gray dashed line total relic abundance is produced if oscillation takes place during radiation domination.Right: Contours producing observed relic abundance, in case when the radiation has non-minimal contribution to radiation for f a = {10 11 , 10 12 } GeV depicted by black and red contours.The solid and dashed patterns represent θ i = 0.1 and θ i = 1 respectively.The red shaded region is ruled out from BBN bound on ∆N eff for overproduction of GW.The gray shaded region is ruled out from BBN bound on T rh ≳ 4 MeV.The gray dashed curves correspond to exclusion limits from future GW detectors.In both cases the region of parameter space of our interest (T osc > T rh ) is shown by arrowheads.

Figure 4 .
Figure 4. Top: The black contours satisfy the observed relic abundance in case of ALPs for different choices of k shown by different patterns for the minimal reheating scenario.Here we choose θ i = {0.1,1} in the left and in the right panel respectively.The blue vertical straight lines correspond to T osc = T rh for each k.The gray shaded region is discarded by super-Planckian bound on f a .Bottom: Contours producing observed relic abundance for θ i = 1, m a = 1 eV (left panel) and m a = 1 GeV (right panel), for different choices of f a shown with different patterns.Other constraints are same as those appearing in Fig. 3.In all cases the region of parameter space of our interest (T osc > T rh ) is shown by arrowheads.

Figure 5 .
Figure 5. Parameter space in [T rh − k] plane, where ALP produced via standard misalignment can co-exist with (semi-) relativistic ALPs produced from inflaton scattering during reheating.The gray dashed contours (labeled as Ω grav a ) correspond to underabundance for ALPs produced purely via gravitational scattering.The red dashed contours provide right relic abundance for ALPs produced entirely from misalignment during reheating for θ i = 0.5 (lower dashed contour) and θ i = π/ √ 3 (upper dotted contour).We take m a = 10 −2 MeV and f a = 10 10 GeV for all curves.The region of parameter space of our interest (T osc > T rh ) is shown by arrowhead.

8 k = 1 2 Figure 6 .
Figure 6.Left: Parameter space producing right relic abundance, considering gravitational reheating for QCD axion for E/N = {0, 8/3}, corresponding to KSVZ and DFSZ model, respectively.The red thick is disallowed from ∆N eff bound from Planck on Ω GW , while the blue thick band is the viable part of the parameter space.The green shaded band shows axion DM parameter space for standard misalignment during radiation domination.We project limits from a few proposed and existing axion search experiments and astrophysical bounds.Right: Same as left panel, but for ALPs, where the blue dashed lines correspond to right relic abundance produced from standard misalignment during radiation domination for |c γ | = 1 and θ i ∈ [0.5, π/ √ 3].The cyan and orange shaded bands show parameter space where observed DM abundance can be obtained for misalignment during reheating for k = {8, 12}, with T rh = {T BBN , 10 GeV} respectively.In all cases we scan over θ i ∈ 0.5, π/ √ 3 .