Revisiting signatures of thermal axions in nonstandard cosmologies

We revisit the formation of a thermal population of hadronic axions in nonstandard cosmologies, in light of the recent developments in obtaining continuous and smooth interaction rates for both the gluon and photon couplings. For certain cosmological histories, such as low-temperature reheating (LTR) and kination-like scenarios, the thermalization of the axion can be severely delayed to higher masses. In the case that thermal equilibrium is achieved, we improve the constraints on LTR for axion masses around the eV scale with respect to previous works and we constrain for the first time early matter-dominated cosmologies. We also point out the possibility of having the coexistence of cold and warm dark matter populations of axions in kination-like scenarios in the eV mass range.


Introduction
Axions are indisputably excellent cold dark matter (DM) candidates.They can be produced by the so-called misalignment mechanism, see, (e.g., Refs.[1][2][3]) or by the decay of topological defects such as cosmic strings and domain walls [4][5][6][7], in a wide range of masses in the sub-eV ballpark.In the standard cosmological scenario (SC), axions can make up the DM density in the range of ∼ µeV.However, this can be extended to smaller masses in scenarios such as kinetic misalignment [8][9][10], or in nonstandard cosmologies (NSC) featuring periods with an injection of entropy [11], typically due to the decay of a heavy long-lived particle [12][13][14][15][16][17][18][19][20][21][22][23][24][25], or by Hawking evaporation of primordial black holes [26][27][28][29].However, if the misalignment mechanism occurred during a non-standard period such as kination domination (KD) (or in a more general scenario with an equation of state parameter ω > 1/3), the observed DM relic abundance can be reproduced with masses near or around the eV scale [17,24].Furthermore, a thermal population of axions can also be produced from scatterings and/or decays of particles from the Standard Model (SM) to which they couple.In particular, a hot thermal population of axions forms because of the existence of the model-independent axiongluon coupling.The effect of such a population depends on the mass of the axion.Axions with masses smaller than ∼ 0.1 eV decouple from the SM thermal bath at temperatures above the QCD phase transition (QCDPT) and therefore experience an important entropy dilution due to the change in the relativistic degrees of freedom, or they may never reach thermal equilibrium at all, making them difficult to detect.Axions with higher masses that later decouple are more likely to have abundances that can be detected in our observations.In addition to the axion-gluon coupling, there is a coupling to photons; even though it is model-dependent, it is widely exploited to search for axions and axion-like particles.
The hot axion population and its impact on observable cosmology attracted attention long ago.One of the first approaches to study their phenomenological implications and potential sensitivity was made in Ref. [30], where the production of thermal axions above the electroweak scale via quark and gluon interactions was considered.In Refs.[31,32] thermal effects were reconsidered, with the finding that the coupling of axions to the top quark had the most significant contribution to the production rate [32].Studies below the electroweak scale (although far from the QCDPT) were performed in Refs.[33,34].Below the QCDPT, hadron processes are relevant, and they were studied in Refs.[35][36][37][38][39][40].That is, a hot axion population has been computed for a certain range of temperatures, and thus for a given mass range.A first approach to smoother calculations was made in Ref. [41], where the thermal production was calculated across the electroweak phase transition.A further improvement that allows consideration of a wider range of axion coupling to gluons across thresholds was made in Refs.[42,43].This complete interaction rate has been used in Refs.[44][45][46][47] to further constrain axion-like particles and axions.Ref. [44] also considered the coupling to photons and made a similar treatment to the one from Ref. [43], to smoothly extend the interaction rate to high temperatures.
Once the axion population in full thermal equilibrium is established, it is expected to decouple due to its feeble interactions while still relativistic.A relic that decouples in the early Universe while being relativistic can impact, on the one hand, Big Bang Nucleosynthesis (BBN) through the extra contribution to the effective number of neutrinos N ν ≡ 3 + ∆N ν , where the recent constraints set N ν = 2.880±0.144[48,49].On the other hand, if they are still relativistic before the photon decoupling epoch, they will leave their imprint on the cosmic microwave background (CMB) power spectrum.Here, the impact is also parameterized as the effective number of neutrinos N eff = 3.044 + ∆N eff , measured as N eff = 2.99 ± 0.17 by Planck 2018 [50].The CMB stage 4 experiment is expected to lower the sensitivity to ∆N eff = 0.06 at 95% CL [51].
Besides their impact on the CMB, hot/warm relics impact the Universe's current energy density, affecting the large-scale structure (LSS) spectra.They can potentially cause a reduction in amplitude for small scales in the matter power spectrum because of their freestreaming behavior.The suppression of small-scale matter fluctuations is controlled by the comoving free-streaming length of these relics when they transition from relativistic to nonrelativistic states.
The settlement and evolution of both thermal and non-thermal axions crucially depend upon the cosmological conditions at the moment of and after their generation.An intriguing scenario arises when the Universe undergoes a nonstandard cosmological era preceding BBN [11].A particularly interesting case corresponds to an early matter dominated (EMD) period, where the Hubble expansion rate was dominated at early times by a species with an energy density that scales as nonrelativistic matter.Such a species could correspond to a long-lived heavy particle (ubiquitous in UV-complete models) [52][53][54][55][56] or to very light primordial black holes that fully evaporate before the onset of BBN.However, examples with more general equation-of-state parameters can occur for oscillating scalar fields with certain potentials.A special case is that of kination in which a "fast-rolling" field whose kinetic term dictates the expansion rate of the post-inflation Universe, implying an equation of state ω = 1 [57][58][59].Moreover, the transition from an inflaton-dominated to a radiationdominated Universe, that is, the reheating era, is often assumed to be instantaneous and to occur at a very early time (i.e. a very high temperature).This picture cannot be taken for granted, and scenarios with low-temperature reheating can occur in which the reheating era is prolonged.We refer the reader to Ref. [11] for a review on non-standard cosmological scenarios.Previous studies have extensively explored the NSC scenario within the framework of axion physics.Production of cold DM axions has been studied in different contexts in Refs.[12-15, 17, 20, 21, 23-26].Concerning lightweight thermally produced axions, research has predominantly focused on LTR and kination cosmologies, as discussed in Refs.[16,22].
The aim of this work is, on the one hand, to complement and extend the studies done in Refs.[16,22], about the constraints on the production of an axion population in full thermal equilibrium, under the assumption that the Universe had a different expansion period than that due to radiation before BBN.On the other hand, we also point out interesting features that have only been briefly discussed in the literature before.
The axion framework with which we will work is the so-called "hadronic axion models" that have been extensively used in the literature due to the model-independent axion-gluon coupling, the prime example being the Kim-Shifman-Vainshtein-Zakharov (KSVZ) model [60,61].In this model, new heavy singlet quarks carry U (1) PQ Peccei-Quinn charges, leaving ordinary quarks and leptons without tree-level axion couplings.The most stringent constraint for the KSVZ scenario sets m a < 0.28 eV at 95% CL [45], using the complete gluon interaction rate from Ref. [42] and combining the CMB and LSS data.Therefore, to account for the axion-gluon coupling, we also make use of the complete interaction rate from Ref. [42].We also include in our analysis the coupling to photons, which, even though it is model-dependent, has been extensively exploited for QCD axion searches.This interaction rate has also been extended to high temperatures in Ref. [44], for hadronic axions.We will treat both interactions independently, assuming that one of them dominates over the other in producing the thermal population.In this work, we will impose restrictions on the thermal population in two nonstandard scenarios: low-temperature reheating (LTR) and early matter domination (EMD).The constraints will rely on the results on light massive relics presented in Ref. [62].Their analysis considers CMB, LSS, and weak-lensing effects.Our findings confirm that the population is dependent on the cosmological history; therefore, it could evade current and future observations.We comment on the thermalization of the population in different cosmologies and about the distinctive free-streaming signatures and their potential to further explore these nonstandard cosmological scenarios.Finally, we comment on the possibility of a coexistence of the axion cold DM population produced by misalignment and the population from the thermal bath produced via freeze-in, which is within the reach of experimental observation.In this work, we do not consider axion generation by other possible mechanisms, such as decays of topological defects.Such contributions can in principle be significant, but their actual value is still in dispute (see Ref. [63] and references therein).
The manuscript is organized as follows: In Section 2 we review the coupling of axions to gluons and photons, focusing on the resultant interaction rates with the thermal bath.In Section 3 we discuss the thermal axion population produced in a standard radiationdominated cosmology along with the constraints that we will use.In Section 4 we introduce the features of the NSC to be worked with and address the conditions for thermalization in different cosmological scenarios.Finally, we present our results and analysis for the thermally-produced axion population in two NSC scenarios, LTR and EMD.We point out differences with previous works and comment on their origin.Finally, in Section 5 we explore the possibility of the co-existence of a hot and cold axion population produced from the bath and the misalignment mechanism, respectively, in a kination-like history.We conclude in Section 6.

Relevant axion couplings
Axions appear as pseudo-Nambu-Goldstone bosons when the so-called Peccei-Quinn (PQ) symmetry is spontaneously broken at an energy scale f a .Assuming f a is higher than the electroweak scale (the so-called invisible axion model), axion scattering in the early Universe may lead to a dark population, either in or out of thermal equilibrium with the SM bath.For temperatures smaller than f a , the effective Lagrangian density for the axion a contains couplings to gluons and photons and is given by A distinctive feature of the axion is that its mass has a definite relationship with the PQ scale through the topological susceptibility of QCD, which has been evaluated in the chiral limit [64,65], NNLO in chiral perturbation theory [66], and directly via QCD lattice simulations [67], giving where m a corresponds to the mass of the axion at zero temperature.

Coupling to gluons
Let us first examine the coupling to gluons, written as where α s is the strong coupling and G and G are the gluon field strength and its dual, respectively.For temperatures above the QCDPT T ≫ T QCD , the main channels that can thermalize an axion population are the 2-to-2 scattering processes g+g → g+a, q+ q → g+a, q +g → q +a and q +g → q +a, with g being a gluon, q a quark and q an antiquark.The total axion production rate was obtained in Refs.[30,31] by a hard thermal loop approximation.An improvement including higher-order effects was made consistently in Ref. [32], where couplings to quarks were also included.In general, the interaction rate density of axions with gluons γ gg (T ) is proportional to As the Universe evolves and the temperature decreases, quarks hadronize, and therefore other processes become responsible for thermalization.For temperatures T ≪ T QCD , the main channels are a + π → π + π and a + N → N + π, where π stands for π 0 and π ± .Since pions are more abundant than nucleons, the former process dominates over the latter.An analytical expression for such an interaction rate is [37] Γ ππ (T ) ≡ γ ππ (T ) where m π and f π are the mass and the coupling constant of the pion, respectively, A ≃ 0.215, and h(m π /T ) is a monotonically decreasing function for m π /T > 1, as h(0) = 1 [37].C aπ is the dimensionless axion-pion coupling constant given by where r ≡ m u /m d ≃ 0.56 is the ratio of up-and down-quark masses.The validity of this expression has been questioned in Ref. [68] for decoupling temperatures above ∼ 62 MeV, where the effective field theory is expected to break down.Instead, a unitarized next-to-leadingorder chiral perturbation theory can be used, extending the reliability of the interaction rate to decoupling temperatures up to ∼ 155 MeV [69].Alternatively, the axion-pion rate can be extracted directly from the experimental data on pion scattering, by rescaling the corresponding cross sections [70].
A step forward in accounting for this axion rate was made in Refs.[42,43], where a smooth transition of the interaction rate of axions to gluons across the QCDPT was achieved, by computing it above and below the confinement scale and interpolating between the two regimes.They also treat the threshold at the mass of the heavy PQ fermion Ψ, above which the binary collisions of Ψ become the dominant production channel for the axions.Their result and the corresponding comparison with the interaction rates mentioned above can be seen in Fig. 1a for the KSVZ model.The vertical dotted blue lines correspond to T = 62 MeV and T = m Ψ .In the following, we will use the continuous smooth rate of Ref. [42] for convenience.To account for the importance of using the continuous rate across the QCDPT we present in Appendix A a comparison of the bounds presented in this work with those that would be obtained if only the axion-pion interaction was used; cf.Eq. (5).Although the interpolated rate suffers from greater uncertainties, we have explicitly checked that using the rates from Refs.[69,70] does not significantly alter the limits nor the features presented in the next sections.

Coupling to photons
The coupling of axions to photons is undoubtedly the most exploited one to search for these particles, even though the strength is model dependent.It reads where F is the electromagnetic field tensor and F its dual.The coupling constant g aγγ is given by Here α em is the fine-structure constant, while E and N are the electromagnetic and color anomalies of the axial current associated with the axion, respectively.There are two benchmark axion models.In a model known as the Dine-Fishler-Srednicki-Zhitnitsky (DFSZ) model [71,72] at least two additional Higgs doublets are needed, and ordinary quarks and leptons carry PQ charges.The KSVZ model [60,61] instead features new heavy singlet quarks carrying U (1) PQ charges, leaving normal quarks and leptons without tree-level couplings to axions.For the KSVZ models E/N = 0, while for the DFSZ models, E/N = 8/3.There are, however, extensions to these models, in particular to the KSVZ case, that can spawn several other values of E/N (see, e.g., Ref. [73]).In this work, we assume that the axion comes from a KSVZ-like model and therefore does not couple to SM fermions directly.The axion-photon coupling contributes to the formation of a thermal population of axions via the Primakoff effect.The rate of axion production due to scattering in a multicomponent charged plasma is given by [74,75] where is the effective number density of charged particles, with Q i the charge of the i th particle species.g Q (T ) accounts for the effective number of charged relativistic degrees of freedom.In the hot early Universe, the photon has an effective mass (plasmon mass) given by m γ (T ) = g 1/2 Q T /(6 α em ).Fig. 1b shows the interaction rate density defined as γ Q ≡ n eq a Γ Q (T ).
Before closing this section, we note that the coupling of axions to photons allows for axion decay, with a rate given by [63] which implies that axions above m a ≳ 20 eV have decayed before the present time.
From the above analysis, it can be seen that at very high temperatures, above the new heavy fermion mass, the Primakoff interaction is efficient in thermalizing axion interactions, but below that scale, processes involving the gluon coupling take over up to temperatures similar to the pion mass, where the interaction ceases.Thus, for masses below the eV scale, in a radiation-dominated Universe, the gluon coupling is responsible for keeping thermal equilibrium to lower temperatures.
3 Thermal axions in standard cosmology

The population of thermal axions
If the interaction rate with the SM particles is strong enough, axions reach chemical equilibrium with the primordial plasma in the early Universe1 .However, when the interaction rate becomes smaller than the expansion rate of the Universe, axions decouple from the thermal bath and freeze out.An estimation of the decoupling temperature T d can be obtained from where H R is the Hubble expansion rate in a radiation-dominated Universe, where the SM radiation energy density is with M P ≃ 2.4 × 10 −18 GeV being the reduced Planck mass and g ⋆ (T ) the number of relativistic degrees of freedom contributing to ρ R [76].The decoupling temperature as a function of the axion mass is shown in Fig. 2 with solid gray lines in standard cosmology for axion couplings to gluons (left) and photons (right).The rapid change in slope when m a ≃ O (10 −2 ) eV corresponds to a decoupling temperature around the QCDPT.For the interaction with photons, the decoupling temperature can be found analytically to be GeV, (14) where g ⋆ (T ) corresponds to the number of relativistic degrees of freedom contributing to the SM energy density.At the moment of their decoupling,2 axions have the same temperature as the SM bath, and therefore their number and energy densities are given by since they decouple while being relativistic.As their total number after decoupling is conserved, their number density today is where the conservation of the SM entropy in the standard cosmological scenario was used, g ⋆S (T ) is the number of relativistic degrees of freedom contributing to the SM entropy [76], and T 0 is the photon temperature today.If axions become nonrelativistic before today, their energy density is simply ρ a (T 0 ) = m a n a (T 0 ), therefore, their contribution to the total energy density of the Universe is Equivalently, it can also be expressed in terms of today's axion temperature T a,0 , by noticing that after their decoupling, the axions leave thermal equilibrium, and their temperature T a redshifts until today in terms of the photon temperature as3 where R is the cosmic scale factor and R d ≡ R(T d ).
Alternatively, if axions remain relativistic until today, their energy density at present is By equating Eqs. ( 18) and ( 21), we find that the transition between axions being relativistic and nonrelativistic occurs at a mass of which shows that axions of masses roughly below 10 −3 eV are still relativistic today.Due to their early detachment from the thermal bath, particles with masses below the meV undergo significant entropy dilution from the SM, resulting in a negligible impact on the relic abundance and, consequently, on ∆N eff , as discussed in the subsequent subsection.This suppression is further pronounced during non-adiabatic expansion because of entropy dilution.Hence, the mass target for this study is set above the meV scale.

Constraints on thermal axions
In recent years, there has been extensive effort focused on constraining light, massive relics beyond the SM by integrating observations across various cosmic epochs.One approach has been to combine data from the early Universe, such as the Cosmic Microwave Background (CMB), with measurements from the late Universe, including galaxy clustering and gravitational lensing [77,78].This integration has attracted increased attention due to some preliminary observations suggesting the potential necessity of a certain amount of warm DM to address existing tensions, such as the H 0 and σ 8 tensions [79,80].
Restrictions on the population of any thermal relic are usually separated depending on when in the cosmological timeline they transitioned to become nonrelativistic.Let us estimate this transition by comparing the momentum of the axions to their mass where T a,nr is the temperature of the axion at the moment it becomes non-relativistic.As the temperature of the axion only redshifts, we can trade the temperature T a,nr for the corresponding scale factor at that time, T a,nr = T a,0 /R nr , and therefore For coupling to gluons, axions of mass m a ≳ 1.4 eV turn non-relativistic around the moment of matter-radiation equality, which corresponds to a decoupling temperature of around T d ≲ 50 MeV (see Fig. 2).Instead, for the Primakoff interaction, this happens for axions of m a ≳ 0.8 eV, that decouple at a temperature T d ≲ 400 GeV.Relativistic axions at matter-radiation equality contribute to the number of effective neutrinos present at the CMB decoupling.On the other hand, for those that transition already into a nonrelativistic state, their velocity dispersion can still be high enough such that there is an effective freestreaming scale below which perturbations in the relic component are suppressed.If they constitute a significant fraction of the energy budget, this suppression hinders the growth of matter perturbations on small scales, and they can still be constrained.

Hot axions
The relativistic and semi-relativistic populations of axions at CMB decoupling behave as dark radiation and contribute to the so-called effective number of neutrinos at that time, which can be accounted for using where the contribution from a dark relic is usually parametrized as with In the case of a relativistic relic that was once in full thermal equilibrium, its contribution can be written as Thus, for axions that decoupled well before the QCDPT, the number of relativistic degrees of freedom reaches a plateau at g ⋆S (T d ) = 106.8,and their contribution to the effective number of neutrinos is always ∆N eff ∼ 0.027.The recent value of N eff = 2.99 ± 0.17 has been found, with a 95% CL upper limit from Planck 2018 of ∆N eff ≲ 0.35 [50,81].In Ref. [45], Planck temperature and polarization data was used to set the bound ∆N eff < 0.3, which restricts the mass of the axion to m a < 1.04 eV.On the contrary, by introducing LSS data from Planck lensing and BAO (6dFGS, SDSS MGS and BOSS DR12 surveys) the bound drops to ∆N eff < 0.23 yielding m a < 0.28 eV, reinforcing the strength of the LSS data to constrain light relics.Axions lighter than ∼ 0.1 eV escape detection because they decouple above the QCDPT and, therefore, their abundance is strongly diminished.The coupling to two photons, on the other hand, is not effective enough to place a constraint on relativistic axions during the CMB decoupling epoch.
Upcoming CMB experiments, such as SPT-3G [82] and the Simons Observatory [83], will soon improve Planck's precision in N eff .In particular, CMB-S4 [84] and CMB-HD [85] will be sensitive to a precision of ∆N eff ∼ 0.06 and ∆N eff ∼ 0.027 at 95% CL, respectively.As calculated in Ref. [86], a combined analysis from BBN and CMB results in N eff = 2.880 ± 0.144.The next generation of satellite missions, such as COrE [87] and Euclid [88], shall impose limits at 2σ on ∆N eff ≲ 0.013.Furthermore, as mentioned in Ref. [89], a hypothetical cosmic-variance-limited CMB polarization experiment could presumably reduce the limit to as low as ∆N eff ≲ 3 × 10 −6 , although this seems experimentally challenging.

Warm axions
A second restriction on the mass of hot relics comes from their impact on the total energy density of the Universe today.The hot massive relic components in the Universe are included in the neutrino density of the Universe, which is given by Ω ν = m ν /(93.14 h 2 eV).Planck 2018+lensing+BAO data have set a stringent limit of m ν < 0.12 eV.Subtracting the contribution of active neutrinos considering the minimal mass allowed by neutrino flavor oscillation experiments m ν > 0.06 eV [50,63], we obtain that the contribution from a different kind of hot relic must satisfy From Eq. ( 18) this yields m a < 0.35 eV.
The hot DM component streams away from structures below its free streaming length, reducing the growth of matter fluctuations at smaller distances.The size of the suppression depends on the present hot DM abundance [90].The free-streaming length λ fs (sometimes called the free-streaming horizon) gives the distance traveled by the hot particle at a given time.Physically, this sets the scale below which collisionless particles can not stay confined in gravitational potential wells.It is defined by While axions are relativistic, they travel at the speed of light and their free-streaming length is simply equal to the Hubble radius.When they become non-relativistic, their thermal velocity drops to for z < z nr .In standard cosmology, sub-eV axions become non-relativistic in the matterdominated era.Their free-streaming length peaks around the moment they become nonrelativistic and then reaches a steady state.Therefore, it is customary to approximate λ fs ≃ λ nr fs ∼ (R nr H nr ) −1 which is given by Axions with larger masses (≳ eV), will become non-relativistic before the radiation-matter transition, assuming a standard cosmology.Integrating the above for t > t eq > t nr we obtain The maximum of this expression occurs around t ≃ t eq , and then it reaches a steady value.Using Eq. ( 24), the free-streaming at the moment of matter-radiation equality, in terms of the mass of the particle and today's temperature (related to today's thermal velocity) can be written as The free-streaming for particles in this mass range is important up to the point when the particles are indistinguishable from a cold DM relic and their free-streaming distance is negligible.

Cold axions
Thermal axions will behave as a cold relic at matter-radiation equality -thus contributing to the cold DM energy density -roughly when their velocity satisfies ⟨v a ⟩| zeq ≪ 1, i.e. for masses well above the eV scale.They could overclose the Universe if their relic density is higher than the cold DM one, Ω c h 2 ∼ 0.12.Assuming a standard cosmological history, this occurs for the axion-gluon coupling, for masses of m a ≳ 15 eV, with a free-streaming length of λ fs ≲ 44 Mpc.For the photon coupling, the bound saturates at m a ≳ 54 eV, with λ fs ≲ 12 Mpc.

Constraints
To study the phenomenology of the formation of thermal axions, we will assume each interaction acts independently of the other.To derive constraints, we make use of the results presented in Fig. 2 of Ref. [62] for Weyl fermions, which are based on cosmic microwave background (CMB), large-scale structure (LSS) and weak-lensing data.The analysis incorporates the full likelihoods of TT, TE, EE, low E, and lensing from Planck 2018 data [50].
Additionally, weak-lensing data from the Canada-France-Hawaii Telescope (CFHTLens) collaboration [91] are included, consisting of 2-point correlation functions of galaxy ellipticities.Furthermore, galaxy power-spectrum data from the Baryon Oscillation Spectroscopic Survey (BOSS) data release 12 [92] are incorporated.This dataset contains spectroscopic information from a large number of galaxies, divided into two redshift bins (z = 0.38 and z = 0.61).The analysis focuses on the linear regime, up to a maximum wave number of k max = 0.25 h Mpc −1 , which is within the validity range of the perturbative approach.The free streaming of light relics on small scales suppresses the growth of CDM+baryon fluctuations at the linear level, with a scale-dependent suppression characterized by k fs = 2π/λ fs , where the amplitude is set by the relic abundance, Ωh 2 .The analysis of Ref. [62] is distinguished from previous works on light massive relics in two key aspects: first, the inclusion of the full-shape galaxy data shows that they improve their constraints concerning the usual CMB+BOSS-BAO by approximately 30%.Second, the inclusion of weak-lensing data breaks the degeneracy between the CDM and the thermal relic abundance that would otherwise be present.That is, weak-lensing data prefer a fixed value of CDM (close to the ΛCDM one), pushing the mass of the light relic to smaller values.
In their analysis, it is assumed that the relic was once in full thermal equilibrium with the primordial bath.A standard cosmological history is implicitly assumed, and a prior is set with the lowest possible temperature for a scalar relic at 0.91 K.The constraints obtained are currently the tightest, with lower mass limits of m X ≥ 2.3, 11, 1.1, 1.6 eV for Weyl fermions, scalars, vectors, and Dirac fermion relics, respectively.
To map the constraints from Weyl fermions to axions, it is assumed that both behave indistinguishably (at least on linear scales) in terms of relic density, contribution to ∆N eff , and today's average velocity.This is ensured by the following mappings: where m W and T W are the mass and temperature of the Weyl fermion thermal relic.All in all, the results of Ref. [62] for light massive scalars once in full thermal equilibrium -referred to as Planck18+BOSS-FS+WLens -are depicted in Fig. 3. Data support everything within the blue region at a 95% CL.
Ref. [62] analyzed masses up to ∼ 10 eV, since higher masses lead to a free-streaming scale beyond the validity of the linear regime (k max = 0.25 h Mpc −1 ).To connect with the results of previous studies, Refs.[16,22], we extrapolate to larger masses, by assuming that the 95% contour converges to a line determined by the axion free-streaming set by λ min = 2π/k max ≃ 35 Mpc, this is shown by the dashed-dotted line in Fig. 3.
To gain deeper insights into the data presented in Ref. [62] and their implications for light relics, in Fig. 3 we have added two lines representing the limits imposed by ∆N eff < 0.3 and the hot/warm relic abundance Ω h h 2 ≲ 0.002, labeled '∆N eff = 0.3' and 'Ω = Ω h ', respectively.
The blue region corresponds to the 95 % CL allowed parameter space for light massive scalar relics obtained in Ref. [62], and the white region is excluded.The black dotted line represents the axion temperature today as a function of the mass for the QCD axion population produced in standard cosmology via their interaction with gluons.The dashed green line is the analog for the coupling to photons.Above m ∼ 10 eV we have assumed an asymptotic value for the free-streaming length of λ fs = 25 h −1 Mpc; see the text for details.
These limits are commonly referenced when analyzing Planck 2018 temperature and polarization data.These two constraints are similar to the allowed region in Fig. 3, however, they are weaker for masses around a few eV as well as below 1 eV.
Additionally, we have introduced the black-dotted and green-dashed lines illustrating the axion temperature today relative to its mass, considering interactions with gluons and photons, respectively, under the assumption of standard cosmological history.Masses ranging from 0.2 eV to 18 eV are ruled out for the gluon interaction, while for the Primakoff interaction, the excluded range narrows to 2.5 eV to 10.5 eV.
In particular, the integration of robust LSS data enables more precise constraints, particularly in the low-mass parameter space below the eV scale.Furthermore, it improves our understanding of the mass region between 2 eV and 10 eV, where the data exhibit a pronounced slope, and is much better than the naive estimate of the hot relic abundance.

Nonstandard cosmologies
The cosmological history of our Universe prior to BBN is so far totally unknown.In the so-called standard cosmological scenario, it is assumed that between the end of inflationary reheating and the onset of BBN, the energy density of the Universe was dominated by SM radiation.Additionally, the end of reheating is also taken at a very high scale, well above the typical scales of the processes studied.However, the energy density of the post-inflationary Universe could have been dominated by something other than SM radiation, resulting in a period of expansion that deviates from standard cosmology [11].The impact on axion DM production during this period has been extensively studied [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29].Here, we aim to take advantage of the progress made to obtain smooth decay rates for the coupling of axions to gluons and photons and to re-examine the status of EMD, LTR, kination, and kination-like scenarios.

Early matter domination
In this setup, the evolution of the background is governed by the system of Boltzmann equations where ρ R and ρ ϕ denote the SM radiation and the NSC-driving field energy densities, respectively, and Γ ϕ is the total decay width of ϕ into SM radiation.The Hubble expansion rate H is given by For an EMD, ϕ dominates the energy density of the Universe between T eq > T > T end , with T end > T BBN in order to not spoil the successful predictions of BBN.The temperature at the end of the NSC is defined as the stage at which the field has mostly decayed away, i.e., when 3 H(T end ) ≃ Γ ϕ , and therefore [15, 93] The nonstandard expansion has two different phases, an adiabatic one, where the ϕ field dominates the expansion without effectively decaying, and then a nonadiabatic one, where entropy is injected into the SM thermal bath due to the decay of the field.We denote the transition temperature between these periods as T c .The evolution of the Hubble parameter can be conveniently expressed as [15,24,94] for This allows us to find the relationship between the scale factor and temperature during the different phases, given by where R c ≡ R(T c ) and R end ≡ R(T end ).The total entropy injected into SM radiation can be estimated to be We note that if a process such as axion decoupling occurs during the nonadiabatic phase, only a portion of the total entropy injection is felt, with T c in this expression being replaced by the decoupling temperature.The evolution of the energy densities for the SM radiation and ϕ as a function of the scale factor is shown in Fig. 4a.It was obtained numerically by solving Eqs. ( 36) and ( 37), for T eq = 6.6 × 10 4 GeV and T end = 0.02 GeV (which implies T c = 0.2 GeV).The deviation of the scaling ρ R (R) ∝ R −4 for free radiation between R c and R end is due to the effective decay of ϕ that acts as a source term for SM radiation.

Low-temperature reheating
The low-temperature reheating (LTR) scenario can be recovered from the above in the limit of large T c , with ρ R (T c ) = 0, and identifying T end as the inflationary reheating temperature.Therefore, before the onset of the radiation-dominated era, only the nonadiabatic phase is present, corresponding to the inflationary reheating period.In this case, ϕ in Eqs.(36) and ( 37) is identified as the inflaton.

Kination-like scenarios
In this case, we assume that the state responsible for the NSC has an energy density that redshifts faster than radiation, having an equation-of-state parameter ω > 1/3.As its energy density eventually becomes subdominant with respect to the SM radiation, it does not have to decay.The Boltzmann equations describing the system can be written as As ϕ does not decay into SM particles, the SM entropy is conserved, and therefore the bath temperature scales as where R end ≡ R(T end ) corresponds to the scale factor at the moment of the equality ρ R (R end ) = ρ ϕ (R end ).In this scenario, the Hubble expansion rate can be approximated by A typical example of this scenario corresponds to kination [59,95], where ω = 1.However, larger equation-of-state parameters are also possible.In cosmologies with ω = 5/3, the Figure 5: Parameter space (below the lines) where axions reach chemical equilibrium with the SM plasma, for the SC, EMD (with T eq = 10 7 GeV), LTR, KD and ω = 5/3, and couplings to gluons (left) or photons (right).In the right panel, the SC, EMD, and LTR cases all coincide as described in the text, and therefore we only show LTR to represent these three, while KD and ω = 5/3 are shown separately.energy density of the state generating the NSC scales like R −8 , and therefore the Hubble expansion rate in the NSC phase scales as H(R) ∝ R −4 [96].This equation of state appears in the context of ekpyrotic [97,98] or cyclic scenarios [99][100][101][102]; see also Ref. [103].Fig. 4b shows the energy densities for ϕ and the SM radiation energy density as a function of the scale factor for ω = 1 (KD) and ω = 5/3, for T end = 0.02 GeV.
While the analytical estimates provided above offer valuable insights into the role of the NSC in the Universe's expansion, we emphasize that all subsequent calculations have been done by numerically solving the corresponding Boltzmann equations ( 36) and (37), and using the Hubble parameter of Eq. ( 38).

The population of thermal axions
In the standard radiation-dominated period, both axion interactions under study (those due to gluon and photon couplings) can be strong enough to produce an axion population in thermal equilibrium.In Fig. 5 we explore this possibility for a set of non-standard cosmological histories including LTR, EMD, KD, and, for the sake of completeness, a cosmology with ω = 5/3.For parameter values in the regions below each curve, thermal equilibrium is achieved at some temperature T ≥ T end when T end < m Ψ , or at some temperature T ≥ m Ψ when T end > m Ψ .In the left panel, corresponding to the axion coupling to gluons, we can see from the solid black line that, for a standard radiation-dominated Universe, thermal equilibrium is guaranteed for f a < 3 × 10 10 GeV when the mass of the heavy colored PQ fermion is m Ψ = 10 5 GeV.
For the gluon interaction, assuming the validity of the interpolated rate presented in Sec-tion 2.1, it is important to note that the highest f a at which equilibrium can be achieved is determined by the mass of the PQ fermion, which we have fixed to m Ψ = 10 5 GeV.For temperatures higher than this mass, the rate scales approximately to T .Consequently, for T end > m Ψ , all curves in the left panel converge, since thermalization is achieved in the standard RD period that follows each particular history.For lower ending temperatures, T end ≲ m Ψ , the curves depart from standard radiation-dominated cosmology.In the LTR case, with H ∝ T 4 , the features of the interaction rate at the QCDPT and the heavy fermion mass threshold are clearly seen as the two "bumps"at the corresponding temperatures, indicating which part of the rate is responsible for ensuring equilibrium.The change in slope at T end ≃ 1 TeV corresponds to a case where equilibrium is instead being achieved by the gluon scattering rate at T ≃ T end < m Ψ , rather than the rate near the QCDPT or the heavy fermion mass during the nonadiabatic phase.The case of KD, which has a milder dependence of H ∝ T 3 also displays these features, although to a lesser degree, whereas the ω = 5/3 case, which also has H ∝ T 4 essentially mimics LTR.Finally, the case of EMD is less clear, due to the prior period of RD before the onset of EMD.This results in an additional dependence on the temperature at the beginning of the EMD, which was taken to be 10 7 GeV in the figure.For this starting temperature, the EMD curve follows along the LTR case at high T end , corresponding to H ∼ Γ occurring in the nonadiabatic phase with H ∝ T 4 , before becoming constant for T end < 3 × 10 4 GeV in which equilibrium is established in the adiabatic phase, with H ∝ T 3/2 .If the start of the EMD is taken to be at temperatures larger than 10 7 GeV, the curve would continue to follow the LTR case down to smaller f a , while if the start temperature is smaller, the EMD curve approaches the RD line, fully replicating it for T eq < m Ψ .
Based on the Γ ∝ T 3 /f 2 a scaling of the gluon interaction rate for T ≲ m Ψ , and our f a < 3×10 10 GeV value for m Ψ = 10 5 GeV, we can estimate the behavior of the curves for different values of the PQ fermion mass.The limiting value of f a below which equilibrium is achieved for a standard RD history is approximately given by f a ≈ 3 × 10 10 GeV m Ψ /(10 5 GeV), which increases with m Ψ .Therefore, we expect that the upper right portions of the curves will shift upward and to the right as m Ψ increases, maintaining their relative positions, while the lower left portions remain unchanged.In principle, if m Ψ is larger than all inflationary reheat temperatures being considered, one could integrate out the heavy fermions so that they decouple from the low-energy theory.This would result in a continuation of the T 3 scaling of the interaction rate from gluon scattering (as seen in Fig. 1a) to temperatures larger than T end for all T end considered.The feature at T end ∼ m Ψ , where all curves converge, would then be replaced by a scaling similar to that seen in the photon interaction discussed below, where larger f a can establish equilibrium as long as T end is sufficiently large.However, we find it informative to study the case where m Ψ is accessible precisely because it introduces an upper limit on the scale f a above which thermalization from the gluon interaction is no longer guaranteed for any cosmological history, independent of the value of T end as long as it is larger than m Ψ .
For the Primakoff interaction, shown in the right panel of Fig. 5, the interaction rate maintains the scaling Γ Q ∝ T 3 /f 2 a throughout the cosmological history.For the purposes of determining the largest f a that can accommodate equilibrium in a standard RD history, this introduces a dependence on the highest temperature reached in the RD period, namely the inflationary reheat temperature.Therefore, the standard RD case is essentially the same as LTR, albeit with a higher reheat temperature.The largest f a that allows for equilibration can then be estimated by evaluating H R = Γ Q at the inflationary reheat temperature using Eqs.( 9) and ( 12), giving f a ≃ 1.5 × 10 9 GeV T end /(10 7 GeV).In the figure, we show a single curve for both the RD and LTR cases, with T end corresponding to the inflationary reheat temperature.In the case of EMD, where the matter period is preceded by a radiation phase, the curve is the same as in the case of RD/LTR, with the understanding that T end is again the inflationary reheat temperature rather than the end of EMD.The cases of KD and ω = 5/3 also result in essentially the same curve, but in this case both the Hubble rate and the photon interaction rate scale as T 3 for KD, while H ∝ T 4 for ω = 5/3.Instead of the reheat temperature, T end now indicates the end of the kination-like period, which still corresponds to the maximum temperature reached in the radiation phase.We should note that this discussion applies only to finding the maximum value of f a below which equilibrium is guaranteed.Here, we have not treated the differences between the histories when it comes to the decoupling temperature.

Signatures of thermal axions formed during NSCs
In this section, we discuss the changes in observables related to axion thermal production under two different cosmological scenarios: LTR and EMD.We focus on these cosmologies, on the one hand, to make contact with previous works [16,22].On the other hand, cosmological periods with ω < 0 will differentiate from EMD mainly in the total amount of entropy injected into the primordial bath.That is, they further dilute the thermal population.Kination and similar cosmologies (ω > 1/3) do not deposit extra entropy, as they dilute faster than radiation.Therefore, the only change from SC is a higher decoupling temperature, reaching almost the same constraints as in the standard scenario.
We will constrain the LTR and EMD periods using the data in Fig. 3.When the axion population undergoes thermalization and decouples during an NSC, the decoupling temperature is higher than that of the RD scenario.This is due to the increased Hubble parameter, resulting in a higher dilution caused by the change in the relativistic degrees of freedom in the SM.Fig. 2 illustrates the comparison of decoupling temperatures between an RD Universe (solid gray) and different EMD histories for the gluon coupling (left) and photon coupling (right).
However, in the case of EMD and LTR there is also an entropy dilution from the new field ϕ when it begins to decay.To incorporate this effect into the observables discussed in Section 3, we focus on the axion temperature today, since we can express all observables in terms of this variable; see Eqs. ( 19) and (28).
We assume that the decoupling occurs at bath temperature T d,nsc > T end .After decoupling, the axion temperature will then continue to redshift as To relate this temperature to the current temperature, we consider the nonadiabatic period between T end < T < T c using the injection of entropy in Eq. ( 42).The axion temperature today is given by The second case in the above equation also applies to an LTR cosmology, where T c is a very high scale.For decoupling temperatures smaller than T end , the expression coincides with that of RD, Eq. (20).

Observables in NSC
We first write down the abundance of axions today, in the case where the decoupling happens during a NSC.We go back to Eq. ( 19) and replace today's temperature for the one in an NSC scenario, Eq. (48).Therefore, where in the last case, it is implied T d,nsc = T d .As a result of the smaller energy, the bounds relax with respect to the standard scenario.
The same strategy can be used for the contribution to the relativistic degrees of freedom during the CMB decoupling and the velocity of thermal axions after they become nonrelativistic.In the first case, we obtain while for the velocity, we use Eq. ( 31), It is interesting to note that the thermal velocity is the observable that is least affected by the dilution.Therefore, we expect that even though the constraints on axions of small mass (dark radiation) can be severely lifted in non-standard cosmological histories, they can still leave a distinctive imprint through their impact on LSS observations.In effect, we can easily find the redshift at the approximate moment when the decoupled axion becomes non-relativistic -the analog of Eq. ( 24) -because it can also be written in terms of today's temperature Note that in the above expression, the entropy injection appears, which is always ≳ 1.Therefore, the moment the axions become non-relativistic happens earlier if there is a nonadiabatic phase, leading to an earlier transition into a non-relativistic state.This trend has also been observed in Ref. [22].As an example, in Fig. 6 we show the free-streaming length at the moment of matter-radiation equality as a function of the axion mass for the gluon coupling.We can see that small axion masses, which decouple very early, have a much larger entropy dilution and thus do not follow the standard RD cosmology (solid gray line) because they became nonrelativistic prior to matter-radiation equality.For T end = 100 MeV (dotted blue) this behavior occurs for masses up to m a ≲ 0.2 eV and higher masses already decouple after the end of the NSC.For cosmologies with a higher entropy dilution, such as the dash-dotted red line, masses below 10 eV are already non-relativistic at the equality, with a much reduced free-streaming length.Therefore, the axions become CDM-like for observational purposes.
For several NSC histories, the thermal axion population can have more warm DM features than hot DM.Therefore, it is expected that they impact the matter power spectrum at high wave numbers k fs ≡ 2π/λ fs .As a consequence, most of their effects appear at the nonlinear scale, where the perturbative analysis breaks down.This distinctive feature can be the key to distinguishing hot relics formed during an RD history from those formed during NSC.
Let us now return to the general picture of the impact of an NSC on the population of relic thermal axions.The energy dilution of the population is the most significant feature of the NSC analyzed here.In the case of LTR, a larger loosening of the constraints will occur as  long as the entropy injection is also larger.This situation arises if T d,nsc is as high as possible compared to T end .Therefore, we expect the highest relaxation of the bounds for small values of the reheating temperature and small masses.In the case of EMD cosmologies, the decoupling can occur at three different stages of the NSC history: before the decay of the ϕ field, during the decay, or afterward.In the first case, the injection of entropy could be the highest possible (and thus loosen the constraints), especially for higher T c (or equivalently T eq ) and smaller T end .This will happen for long NSC periods and small masses (high decoupling temperatures).The second possibility has the same outcome as the LTR scenario; therefore, the less constrained scenario comes from small masses and small T end .The third possibility corresponds to decoupling during an RD scenario, so all constraints from Section 3 apply.With all these considerations, we will now explore the parameter space that can be constrained for EMD and LTR cosmologies using the same data of Ref. [62] we used previously for the standard cosmology.

Constraints from light massive relics
We commence obtaining constraints on the formation of a thermal axion population in the early Universe for LTR cosmologies for the axion-gluon coupling.To do so, we compute for each mass m a and reheating temperature T end , today's axion temperature from Eq. ( 48) and compare with the data shown in Fig. 3, which we remind the reader includes Planck, fullshape LSS data from BOSS DR12 and weak lensing from CFHTLens.The results appear Figure 7: Parameter space as a function of the axion mass/decay constant and the NSC end temperature, T end , for LTR cosmologies for the axion-gluon coupling.The blue region shows the parameter space where the axion is overproduced with respect to CDM.The red region labeled "Planck18+BOSS-FS+WLens", corresponds to the constraints obtained using the CMB, full-shape galaxy and lensing data from Ref. [62], and is excluded with a 95% CL.
The parameter space where the axions do not achieve full thermal equilibrium is marked in gray and hatched where they decay before today.See the text for more details.
in Fig. 7, where the red area (marked "Planck18+BOSS-FS+WLens") is excluded from the data.The gray area (marked "No thermalization") is the parameter space where the axions do not fully thermalize, in correspondence with Fig. 5a.The blue region corresponds to the parameter space where the abundance of the thermal population exceeds the abundance of CDM today.The streaked region shows the masses above which the axion has a lifetime smaller than the age of the Universe.As was previously anticipated, the bounds on small axion masses are lifted for cosmologies with a high entropy dilution, resulting in a LTR with small T end .The condition T d,nsc = T end is shown as a solid black line.Nonetheless, decoupling temperatures slightly smaller than T end still happen during the NSC period, as the decay of the new field is not instantaneous.For this, we have also included the black dot-dashed line labeled 'T d,nsc = T d ', where the decoupling occurs just at the moment when the radiation content is no longer affected by the NSC.That is, to the right of that line we have decoupling temperatures in SC.
The shape of the regions in the figure can be understood as follows: Higher reheating  temperatures correspond to the decoupling during SC; therefore, the bounds are independent of T end and are the same as in SC, i.e. m a ≲ 0.2 eV.As the reheating temperature decreases, it eventually reaches the line T end = T d,nsc , first at smaller masses that have higher decoupling temperatures.Decreasing T end further leads to no constraints for the smallest axion masses constrained in SC, due to the high dilution lowering the axion's temperature below the data's reach for those masses (for smaller T end , we also see that the small-mass region does not thermalize).For masses around and above an eV, the bounds of the data shown in Fig. 3 become stronger, due to the variety of LSS data used.This allows us to constrain reheating temperatures down to 25-30 MeV for that mass range.Let us point out that the constrained red-colored region lies very close to the line T d,nsc = T QCD , meaning reaching masses where the decoupling temperature is only just below the QCD temperature T QCD ≃ 150 MeV.This reinforces the relevance of using a complete and smooth interaction rate around that temperature.
For EMD cosmologies, we have chosen to show our constraints according to the total amount of dilution possible from the EMD, characterized by S(T c )/S(T end ) given in Eq. (42).
Therefore, to keep that ratio fixed, we vary T end and also T c (or equivalently T eq ) for each axion mass.In Fig. 8 we show three different scenarios for the gluon coupling: the red and blue regions bounded by solid lines correspond to EMD cosmologies with a dilution factor of S(T c )/S(T end ) = 2 × 10 −4 , while the regions bounded by dot-dashed and dotted lines correspond to S(T c )/S(T end ) = 0.2 and S(T c )/S(T end ) = 0.7, respectively.We see that in the first high-diluted scenario, we recover the same results as in LTR Fig. 7.This means that the decoupling for all masses that can be probed happens either during the nonadiabatic phase or in SC.Masses that decouple at higher temperatures are not constrained because of the strong dilution.The second and third scenarios, which correspond to shorter EMD cosmologies with less total entropy dilution, can be achieved with T c and T end being close to each other. 4The shape of the parameter space in these scenarios in Fig. 8 can be understood as follows.For high T end , where all three cases converge, the bounds from SC are recovered.As the temperature decreases, it eventually reaches T end ∼ T d,nsc , and the dilution begins to become important.As T end keeps decreasing and T d,nsc < T c , the dilution factor is given by S(T d,nsc )/S(T end ), which is smaller than the total dilution factor S(T c )/S(T end ), so the constraints are stronger.This can be seen for T end ≳ 35 − 40 MeV in the second case.For smaller T end the decoupling -especially for smaller masses -occurs during the adiabatic phase of the EMD.In that scenario, the dilution factor is the highest possible for all T end , and the bound flattens.
Finally, in Fig. 9 we show the corresponding plots for the axion-photon coupling.The left panel shows the case of LTR, where it was already anticipated that thermalization will be highly delayed for this interaction.Therefore, for high reheating temperatures, the bounds are the same as those of SC, as the decoupling happens during that period.As the reheating temperature decreases, it reaches the T d,nsc = T end line (which occurs at temperatures much higher than the gluon counterpart).Smaller reheating temperatures do not allow for axion thermalization, so they cannot be constrained by our analysis.In the right panel, we again show EMD cosmologies with a fixed total dilution factor.The photon coupling has much higher decoupling temperatures than the gluon coupling.For that reason, the smallest dilution taken is S(T c )/S(T end ) = 0.2, enclosed by solid red and blue lines.Smaller dilution factors lead to decoupling in SC and those bounds apply.In this case, the result is not identical to the LTR case due to the difference in thermalization.For this scenario, we can constrain a slim parameter space where the decoupling occurs during the EMD era, and it is better for masses around the eV scale, thanks to the data used.By reducing the dilution, we obtain the dot-dashed and dotted lines, with total entropy dilution factors of 0.4 and 0.7, respectively.The same feature as in the case of the gluon coupling is obtained, in the sense that for T end below T d,nsc , the dilution factor is given by S(T d,nsc )/S(T end ), which is smaller than the total factor, leading to a stronger constraint.As T end continues to decrease, the decoupling occurs in the adiabatic phase, and the bound becomes independent of T end .

Coexistence of cold and hot populations
An intriguing question is whether axions can leave imprints in cosmology from both a cold dark matter (CDM) population and also a "thermally" induced one, in the sense that it emerges from the primordial bath.
In standard cosmology, the CDM population arises around the QCDPT, with a mass range of 10 −6 − 10 −4 eV [104] from the misalignment mechanism.However, a complete thermal population of axions can only emerge due to the gluon interaction for masses above m a ∼ 10 −4 eV, as we have seen in Fig. 5, for m Ψ = 105 GeV.Thus, to have a co-existence of a CDM population from misalignment and a thermal population from the gluon interaction at masses around µeV, the latter must be produced via freeze-in.But for such small masses the energy density is negligible (see e.g.Fig. 1 of Ref. [105]), which makes them unobservable.On the other hand, from the photon coupling an axion mass m a ∼ 10 µeV decouples from the bath at around T d ∼ 10 8 GeV, much earlier than the oscillation temperature of the axion field. 5 This is not the case in NSC scenarios with ω ≥ 1.On the one hand, from Fig. 5 we can see that the thermalization of the axions can be largely delayed to higher masses by decreasing T end .On the other hand, the misalignment production of CDM axions is also altered in cosmologies with ω > 1/3, shifting the correct relic abundance to higher masses.Those cosmologies do not feature entropy injection into the thermal bath; the only net effect is to delay the oscillation of the axion field, leading to an overproduction of the relic abundance for the classical axion CDM mass window [17,24].
In this section, we are interested in finding the parameter space where two populations of QCD axions can actually co-exist in sizable amounts: one cold, via the misalignment mechanism, and the other originating from interactions in the bath, from the coupling with gluons.We will use the results of Ref. [24] where it was found that for cosmologies with ω = 5/3, the misalignment mechanism produces the right amount of CDM observed today for masses smaller than or around the eV range for T end ≳ 4 MeV.In principle, the coexistence could also happen for the ω = 1 KD cosmology, but since the correct CDM abundance takes place for smaller masses, the abundance (and therefore impact) of hot axions will be fairly small to escape detection.We will come back to this later in the section.We will only focus on the axion-gluon coupling because, for the photon coupling, thermalization is largely delayed to masses well above the eV for the range of T end we are interested in, according to Fig. 5b.
In order to find the yield of axions produced through freeze-in, we numerically solve the Boltzmann equation with the change of variables Y ≡ n a /s.Then, the relic density is simply found as ρ a = Y ∞ s 0 m a , where Y ∞ is the yield long after the freeze-in.
For the axion-gluon coupling, co-existence can occur if axions do not fully thermalize due to their active interactions, i.e.Γ < H, otherwise this would lead to the disappearance of the condensate.The population built through freeze-in respects the constraints on hot/cold relics that we have analyzed in the previous section.However, it is not possible to use the results of Fig. 3, as they are only valid for thermal relics, which is not the case considered here.Nevertheless, in order to give an educated guess on whether the co-existence parameter space could be allowed from cosmology, we use constraints on non-thermal light sterile neutrinos presented in Ref. [107].There, a relic population of sterile neutrinos was assumed to be produced nonthermally through the Dodelson-Widrow mechanism [108].There it is used that despite our ignorance about the phase space distribution, the three observables, ∆N eff , Ω a h 2 , and ⟨v a,0 ⟩, satisfy a constraint equation given by Above it has been assumed that the particles are relativistic/semi-relativistic at photon decoupling, which is a fair assumption for masses around an eV.For a nonrelativistic relic today, which is still distinguishable from the CDM, the velocity should satisfy ⟨v a,0 ⟩ ≳ 1 km/s.From the Boltzmann equation, we get the relic abundance of today together with the contribution to ∆N eff .For their analysis, Ref. [107] used WMAP5 data plus small-scale CMB, SDSS LRG data, SNIa data from SNLS, and Lyman-α (conservative) from VHS.
In Fig. 10 we show the parameter space where the co-existence of a misalignment and a freeze-in population can happen for a non-standard cosmology with ω = 5/3 in terms of the axion mass and the temperature where the NSC ends, T end .The blue band corresponds to the parameter space where axions produced through the misalignment mechanism can account for the entire observed DM today, assuming a natural range of initial angles, θ ∈ [0.5, 1.8], see Ref. [24] for details.Masses above the solid black line achieve full thermal equilibrium and those below them freeze in.Contours of equal thermal abundances Ω a h 2 are indicated in dot-dash lines with the corresponding value.
As we can see, most of the parameter space where co-existence happens is unconstrained.However, axions with masses below the eV scale have seemingly small relic abundances, which could lead them to be undetected by current and future surveys.We have checked that all the parameter space scanned in Fig. 10 corresponds to velocities today higher than 1 km/s.Therefore, even though these axions contribute to the total CDM abundance, they are still distinguishable from CDM.The mass range where both populations coexist, with a significant abundance of axions from freeze-in (Ω a h 2 ≳ 10 −4 ), is within the sensitivity range of experiments such as IAXO [109], BREAD [110], and LAMPOST [111].
Finally, a comment on the possibility of generating a similar scenario during kination is worth mentioning.From Fig. 5a, it can be observed that for the smallest ending temperature allowed by BBN, T end ∼ 10 MeV, masses above m a ∼ 0.5 eV thermalize for the axion-gluon coupling in KD.On the contrary, a population of CDM axions will form for that scenario, for T end ∼ 10 MeV, for masses between 10 −3 − 10 −2 eV (see Fig. 6 from Ref. [24]).Hence, it seems possible that the CDM population can co-exist with a nonfully thermal one in KD.However, the relic abundance of the latter is expected to be similar to or slightly higher than the one for a ω = 5/3 cosmology.This is because, despite ω = 5/3 having a higher slope than KD, they converge in the vicinity of T end .This leads to the yield of KD being slightly higher than ω = 5/3.From Fig. 10 we see that the relic abundance is much smaller than Ω a h 2 ∼ 10 −8 for axion masses m a ≲ 10 −2 eV, for T end ∼ few MeV.Thus, even though co-existence is possible during KD, the thermal bath population is not abundant enough to be detectable in the near future.

Summary and conclusions
In this study, we have explored the generation of axions from the thermal bath by their interactions with gluons and photons during cosmological eras characterized by expansion rates different than the standard radiation-dominated scenario, the so-called non-standard cosmologies.Our focus was on KSVZ-like models, where axions lack direct couplings to fermions.
We started in Section 3 with a comprehensive analysis of the formation of a fully thermalized population in standard cosmology and the potential signatures that these particles could leave in cosmological data.Massless or nearly massless axions are well characterized by their contribution to ∆N eff and m ν .On the other hand, light but massive axion relics, while non-relativistic in the present epoch, still maintain non-zero temperatures.This characteristic feature makes it possible to distinguish them from the majority of cold relics, offering a novel approach to identifying tiny relics in cosmological data.Unlike CDM, these thermal relics possess thermal velocities that hinder their clustering beyond a characteristic free-streaming scale.Consequently, they exert a discernible influence on the growth of matter fluctuations, making it possible to detect them through the study of their impact on the large-scale structure of the Universe.Finally, we introduced the data to be used to constrain the fully thermal axion population.We relied on the latest constraints derived from the temperature-mass parameter space for massive light relics once in thermal equilibrium, as reported in Ref. [62].These constraints were obtained through a comprehensive analysis that combined BOSS DR12 full-shape galaxy data, Planck 2018 temperature polarization and lensing anisotropies, and CFHTLens galaxy-galaxy ellipticity correlations.We applied those restrictions to SC in Fig. 3, to find excellent agreement with previous results reported in the literature.
In Section 4 we started by making a detailed analysis of the thermalization of axions in standard and nonstandard scenarios for both couplings to gluons and photons.For the gluon coupling we found that the mass of the heavy PQ fermion sets an important upper limit on the maximum value of the scale f a that can lead to thermalization in a given history.This is in contrast to the photon coupling which becomes increasingly efficient at thermalization as the maximum temperature of the RD phase increases.
Subsequently, we moved towards analyzing the scenario where axions thermalize in LTR and EMD cosmologies for the gluon and photon couplings.We assessed the influence of the NSC scenarios on three key parameters: the relic abundance, Ω a h 2 , the number of extra relativistic degrees of freedom, ∆N eff , and the free-streaming length, characterized by the thermal velocity, ⟨v a ⟩.Although the impact of light massive relics is inherently complex, in many models with non-thermal distortions, the observable effects can be effectively parameterized using these three quantities with considerable accuracy, as discussed in Ref. [112].
Axion thermalization in LTR cosmologies for the gluon coupling has been addressed in the literature before, in Refs.[16,22] and our results are in good agreement.The parameter space where the thermal population can exist without restrictions opens up, especially between 0.2 eV ≲ m a ≲ 10 eV for cosmologies with reheating temperatures smaller than 100 MeV.Our study expands the previous findings in the following way: we use interaction rates for both gluon and photon couplings that are continuous across the QCDPT.The former allows us to scan smaller axion masses that have higher decoupling temperatures.However, the data used to constrain the NSC include the full-shape from BOSS and weak-lensing data as an additional component with respect to Refs.[16,22].Due to the above, in contrast to previous works, we are able to constrain LTR scenarios where the axion freeze-out occurs during the nonstandard expansion, for masses between 1 eV ≲ m a ≲ 15 eV and for reheating temperatures around 25-30 MeV.We present a deeper insight on the complementarity of our limit with previous ones in Appendix A.
A thermal population in EMD cosmologies has not been studied before, to the best of our knowledge.We have established our constraints based on the total amount of entropy injected into the thermal bath.First, for the axion-gluon coupling, we have found that in the case of high dilution, the EMD scenario has the same constraints as LTR.This is because due to the decreased abundance and velocity, the data can only constrain masses that decouple during the nonadiabatic phase.As the total entropy injection decreases, it is possible to probe masses that decouple during the adiabatic phase of the EMD, resulting in smaller T end temperatures.Our results allow us to easily extrapolate EMD with other entropy injection factors.
Next, in Section 4, we repeated the analysis for the axion-photon coupling.Our results show that it is only possible to constrain LTR cosmologies with T end > T d,nsc because otherwise, axions do not thermalize, as seen in Fig. 5b.On the other hand, for EMD cos-mologies, we have found that, due to the high decoupling temperature (see Fig. 2), highly diluted cosmologies hide the population completely, even at high T end .Only for factors S(T c )/S(T end ) ≳ 0.2, constraints appear.We emphasize that a key feature to search for signatures of NSC is their effect on the free-streaming length.In contrast to the severe effects on the contribution to the effective number of neutrinos and the relic abundance, the velocity of the thermal relic is less affected.This means that they can still play an important role in their impact on the formation of LSS.Such an impact is shifted to smaller masses than in SC because of the earlier transition to non-relativistic states.Upcoming large-scale structure surveys such as the Dark Energy Spectroscopic Instrument (DESI) [113], the Vera Rubin Observatory [114], and Euclid [88], together with the next generation of CMB experiments will play a major role in either discovering or placing severe constraints.
Finally, in Section 5 we explored the possibility of having axion signatures from both cold dark matter and hot/warm populations.Axions are naturally produced as a cold condensate by the misalignment mechanism that, in a SC, is well motivated in the mass range around the µeV.However, for that mass range, axions do not thermalize for the gluon interactions and decouple much earlier for the photon coupling.However, it has been pointed out in the literature that for cosmologies driven by ω > 1/3 the formation of a CDM condensate, with the right characteristics demanded by observations, occurs at higher masses.In particular, for ω = 5/3 it can reach the eV mass range.We have found that around that mass, a non-negligible population of hot/warm axions from the thermal bath that is not in conflict with cosmological constraints can also exist.A dedicated analysis could constrain a larger portion of the parameter space than the one we showed, as we used data for nonthermal sterile neutrinos that are outdated.The parameter space in which the two populations can coexist is in the range of laboratory experiments, such as IAXO, BREAD, and LAMPOST.Figure 11: Parameter space as a function of the axion mass/decay constant and the NSC end temperature, T end , for LTR cosmologies for the axion-gluon coupling.The continuous solid black line is the bound using the smooth interaction rate found in Ref. [42] and is labeled 'Γ across '.Instead, the bound using the axion-pion interaction rate from Eq. ( 5) is shown as the dashed black line, labeled 'Γ π π '.Both exclusion regions use the LSS data from Ref. [62].

A Appendix
Changes in the exclusion limit of the axion-gluon coupling by considering a continuous smooth interaction rate In Subsection 4.4, we discuss the importance of employing a smooth rate around the QCDPT, as the data are sensitive enough to constrain decoupling temperatures nearing approximately 150 MeV.Here, we aim to refine this assertion by quantitatively comparing the constraint depicted in Fig. 7, utilizing the smooth rate derived in Ref. [42,43], with the one obtained using the rate presented in Eq. ( 5), the validity of which has been challenged due to unitarity concerns above T ∼ 62 MeV [68].This comparison is shown in Fig. 11.We observe that the differences between the two are approximately 10% − 15% for temperatures below 30 MeV, escalating as the temperature approaches T = 62 MeV, reaching around a 50% difference.This discrepancy is evident in Fig. 11, where the most significant deviation occurs for masses around the eV range, aligning with the decoupling temperatures approximately at that value for T end ≲ 100 MeV.From this figure, we can also stress the importance of the smooth rate in constraining small axion masses and higher T end temperatures.
Changes in the exclusion limit of the axion-gluon coupling by considering more LSS data We present a comparison of our findings with those in the literature, specifically Refs [16,22].They investigated the parameter space for thermal axions within the framework of KSVZ axions, focusing on the gluon coupling, within the context of LTR cosmology.Both works utilized the axion-pion interaction rate outlined in Ref. [37], as presented in Eq. ( 5), noting its applicability up to temperatures of 62 MeV.
To set their constrained parameter space, Ref. [16] used the saturation of the relic abundance from cold dark matter, the ISW effect on the CMB and data from WMAP first year and SDSS measurements of the galaxy power spectrum found in Ref. [37].On the other hand, Ref. [22] updated the parameter space for thermal axions in LTR cosmologies by using CMB observations from Planck 2018 legacy data release [50], together with BAO measurements from BOSS Data Release 12 [92], 6dFGS [115] and SDSS-MGS [116].
In the left panel of Fig. 12, we present a comparative analysis.The excluded parameter space from Ref. [16] is depicted in the hatched area, while the limit from [22] is illustrated by the dashed black line, considering only the segment where the rate's validity is acknowledged.Our excluded region, previously presented in Fig. 7, is represented by the solid black line.Thus, it becomes apparent that, on the one hand, the smooth rate allows the low-axion-mass region to be confined, where the decoupling temperature is higher than 62 MeV, and, on the other hand, the incorporation of more robust LSS data allows for a better constraint of the parameter space overall.In particular, the bound on ∆N eff is refined for sub-eV masses, and in the eV region the LSS data impose significant restrictions, as expected from Fig. 3. Figure 12: Parameter space as a function of the axion mass/decay constant and the reheating NSC end temperature, T end , for LTR cosmologies for the axion-gluon coupling.The hatched area corresponds to the limits imposed by Ref. [16].The dashed black line corresponds to the exclusion set by Ref. [22], where only masses whose decoupling temperature is less than 62 MeV are shown.Finally, the continuous black line is the bound found in this work, using the smooth rate found in Ref. [42] and the data from Ref. [62].

Figure 1 :
Figure 1: (a) KSVZ total interaction rate of axions across both the heavy colored Peccei-Quinn fermion mass m Ψ = 10 5 GeV and the QCDPT, represented by the dashed black line (adapted from Ref. [43]).Solid red and magenta lines show the pion and gluon scattering rates quoted in the text.(b) Interaction rate with photons (Primakoff effect).

Figure 2 :
Figure 2: Decoupling temperature as a function of axion mass considering the coupling of axions to (a) gluons and (b) photons.Gray lines correspond to the standard cosmology, whereas solid lines are for LTR and segmented colored lines correspond to different EMD scenarios.

Figure 4 :
Figure 4: Evolution of the energy densities for radiation (dashed black) and the ϕ field (red) as a function of the scale factor R. (a) In EMD with T eq = 6.6 × 10 4 GeV and T end = 0.02 GeV.(b) In KD (dot-dashed) and ω = 5/3 (solid) with T end = 0.02 GeV, for both cosmologies, the energy density for radiation evolves the same as in SC.

Figure 6 :
Figure6: Comoving free-streaming length evaluated at matter-radiation equality as a function of the axion mass.The continuous gray line corresponds to the standard cosmological scenario.The blue and red lines correspond to EMD cosmologies with T eq ≃ 7 × 10 3 GeV and T end = 0.1 GeV, T end = 0.004 GeV respectively.

Figure 8 :
Figure 8: Parameter space as a function of the axion mass/decay constant and the NSC end temperature, T end , for EMD cosmologies for the axion-gluon coupling.Constrained colored areas have the same meaning as in Fig. 7.We show three different scenarios: the regions enclosed by the solid red and blue lines have a fixed ratio of S(T c )/S(T end ) = 2 × 10 −4 , those enclosed by the dot-dashed lines have S(T c )/S(T end ) = 0.2, while those enclosed by the dotted lines have S(T c )/S(T end ) = 0.7.White regions are allowed.

Figure 9 :
Figure 9: Parameter space for the Primakoff interaction in LTR (left) and EMD (right), the colored areas have the same meaning as in Fig. 7. (a) The axion never decouples during LTR because of the Hubble parameter's temperature dependence.Once standard cosmology is regained (T < T end ), decoupling proceeds as in the standard case.(b) We compare three different EMD cosmologies: the regions enclosed by the solid red and blue lines have a fixed ratio of S(T c )/S(T end ) = 0.2, those enclosed by the dot-dashed lines have S(T c )/S(T end ) = 0.4, while those enclosed by the dotted lines have S(T c )/S(T end ) = 0.7.White and gray regions are allowed.

Figure 10 :
Figure10: Parameter space of T end vs m a where a misalignment and a thermal (freeze-in) population of QCD axions can co-exist for a cosmology with ω = 5/3.The blue region shows the parameter space where the CDM relic abundance can be produced with natural initial misalignment angles (see the text for details).The red region corresponds to the bounds on the sterile neutrinos of Ref.[107] applied to the QCD axion.To the right of the solid black line axions produced via the gluon coupling thermalize.The dot-dashed lines show isocontours of relic abundance Ω a h 2 , marked with the corresponding value.White regions are unconstrained.