Thermal QCD Axions across Thresholds

Thermal axion production in the early universe goes through several mass thresholds, and the resulting rate may change dramatically across them. Focusing on the KSVZ and DFSZ frameworks for the invisible QCD axion, we perform a systematic analysis of thermal production across thresholds and provide smooth results for the rate. The QCD phase transition is an obstacle for both classes of models. For the hadronic KSVZ axion, we also deal with production at temperatures around the mass of the heavy-colored fermion charged under the Peccei-Quinn symmetry. Within the DFSZ framework, standard model fermions are charged under this symmetry, and additional thresholds are the heavy Higgs bosons masses and the electroweak phase transition. We investigate the cosmological implications with a specific focus on axion dark radiation quantified by an effective number of neutrino species and explore the discovery reach of future CMB-S4 surveys.


Introduction
The unexpected invariance of strong interactions under transformations flipping the arrow of time is one of the most challenging puzzles in fundamental physics. This remarkable experimental fact is equivalent to state that Quantum ChromoDynamics (QCD) is invariant, within experimental uncertainties, under the combination of parity and charge conjugation (CP) in agreement with the CPT theorem. Such a CP violation by strong interactions is parameterized by an effective dimensionless parameterθ that is expected to be of order one, but the lack of observation of a fundamental neutron electric dipole moment put the spectacular constraintθ 10 −10 [1][2][3]. Anthropic explanations are not viable [4,5], and understanding this severe inequality is known as the strong CP problem. The Peccei-Quinn (PQ) mechanism [6,7] is one of the most appealing solutions. A new Abelian U (1) PQ symmetry, anomalous under strong interactions and spontaneously broken, plays the role of the main character. At energies much lower than the PQ breaking scale f a , the only residual degree of freedom is an approximate Nambu-Goldstone boson a known as the axion [8,9] that acquires the anomalous coupling to gluons Here, we denote the QCD fine structure constant by α s = g 2 s /(4π), the gluon field strength by G A µν with the index A = 1, . . . , 8 running over the adjoint indices of the color gauge group, and its dual by G Aµν ≡ µνρσ G A ρσ /2. Once strong interactions confine, QCD nonperturbative effects generate a potential that leads to an axion mass [10,11] m a 5.7 µeV 10 12 GeV f a . (1. 2) The PQ breaking scale f a suppresses axion couplings as well. Such a scale is constrained by terrestrial and astrophysical axion searches, and accounting only for the axion coupling to gluons in Eq. (1.1) leads to the rather conservative bound f a 10 8 GeV: the axion must be light and weakly-coupled. The field evolution in the early universe goes through two main phases: the axion is initially stuck by Hubble friction, and once its mass becomes comparable to the expansion rate it begins oscillating around the minimum of its potential. The oscillation amplitude gets damped by the Hubble friction, and the axion settles down at its minimum which is CP-conserving as ensured by the Vafa-Witten theorem [12]: QCD dynamics itself solves the strong CP problem. The energy density stored in the field oscillations can account for the observed dark matter abundance [13][14][15]. Furthermore, axion interactions with standard model (SM) fields are responsible for a plethora of phenomena in the early universe [16] and the target of an intensive experimental effort [17][18][19]. This work investigates a distinct cosmological signal of PQ theories: the production of relativistic axions from scatterings and/or decays of particles belonging to the primordial thermal bath [20]. Given their thermal origin, such hot axions are produced with energies of the size of the bath temperature, and their typical energy will stay of the size of the one for photons as long as they are relativistic. This statement holds regardless of whether they thermalize or not, and the axion mass in Eq. (1.2) ensures that axions produced thermally are still relativistic as late as at recombination.
Such hot axions manifest themselves experimentally as an additional contribution to radiation in the early universe. How do we have access to this quantity? Two key events in the expansion history allow us to bound the energy density stored in relativistic particles. Following a chronological order, the first is Big Bang Nucleosynthesis (BBN) when the thermal bath synthesized light nuclei. The successful agreement between predictions and observations gives us information about the expansion rate at BBN and it bounds additional radiation. This effect is parameterized by an effective number of neutrino species on top of the SM contribution, N ν = 3 + ∆N ν . After the recent measurement of the deuterium burning rate by LUNA [21], Ref. [22] found the constraint N ν = 2.880 ± 0.144.
Another important event is the formation of the Cosmic Microwave Background (CMB) since additional radiation alters the CMB anisotropy spectrum at small angular scales. This effect is also parameterized in terms of an effective number of additional neutrino species N eff = N SM eff + ∆N eff . The SM naive prediction N SM eff = 3 does not hold because neutrino decoupling is not instantaneous [23], and we have N SM eff 3.0440 [24][25][26]. The most stringent constraint, N eff = 2.99 ± 0.17, comes from the Planck collaboration [27].
Future CMB surveys forecast an extraordinary improvement in measuring this quantity, and conservative configurations of CMB-S4 can reach ∆N CMB−S4 eff (1σ) 0.02−0.03 [28,29]. What does this value imply for fundamental physics? Let us consider a scalar field, such as the axion, and let us assume that it reaches thermal equilibrium in the early universe. Even if we take the most pessimistic hypothesis that decoupling happened well above the weak scale, the resulting contribution would be ∆N eff 0.027. If decoupling happened at lower temperatures and/or if the field has a larger number of internal degrees of freedom the expected value is larger. Thus future CMB-S4 surveys are sensitive to any relic light particle that was once in equilibrium with the standard model thermal bath [30,31]. with g SM * s (T D ) the effective number of SM entropic degrees of freedom at T D . We employ two different sets of data for g SM * s (T D ), Refs. [32] (dashed lines) and [33] (dotted lines), and we notice how the treatment of thermal bath has a tiny effect on the final predictions (see App. D for more discussion). This result is valid for a standard thermal history with no significant releases of entropy that would dilute the expected amount population of relativistic axions. Exceptions are possible, as for example around the time of the QCD phase transition (QCDPT) [34] or the electroweak phase transition (EWPT) [35].
An important message from Fig. 1, which is quantified by the expression in Eq. (1.3), is that the later the relic decouples the larger is its contribution to ∆N eff . Planck data already exclude dark radiation that decoupled around the time of the QCDPT, and future experiments will probe hot relics that decoupled earlier. Production of dark radiation around or below the QCDPT is particularly timely and promising for future experiments. This is true for any dark radiation candidate, and in particular for the QCD axion which is the subject of our investigation.
These spectacular projections make rigorous calculations a top priority. Recently, several collaborations revised predictions for the hot axion abundance. Production above the weak scale via quark and gluon scattering was investigated by Ref. [36] with the inclusion of Debye screening effects. Later on, Refs. [37,38] provided rigorous treatments of thermal effects, and Ref. [38] extended the analysis to production via electroweak gauge fields and  Figure 1. Contribution to ∆N eff from light relics that were once in thermal equilibrium as a function of the decoupling temperature T D . Different colors represent different particle spins, and the width of each line corresponds to the two different treatments of the relativistic degrees of freedom provided in Refs. [32] (dashed lines) and [33] (dashed lines). The shaded region is excluded by Planck [27], the dot-dashed lines show the CMB-S4 discovery reach [28,29].
top quark. These studies considered production above the weak scale. Scatterings of heavy quarks below the electroweak scale were analyzed by Ref. [39,40], and this analysis is valid only well above the QCDPT. At lower temperatures, hadron scatterings were considered by Refs. [41][42][43][44][45][46]. The analysis in Ref. [47] bound the axion mass to m a < 7.46 eV and m a < 0.91 eV for thermalization with gluons and pions, respectively. Ref. [48] investigated production via lepton scatterings and decays, which is immune to QCD complications, and it pointed out a possible connection with the so-called Hubble tension [49,50]. The presence of mass thresholds across which production rates change dramatically is an issue rarely addressed in the literature. As a step forward toward the proper treatment of the EWPT, Ref. [51] provided a smooth connection between rates within an effective field theory framework containing only one Higgs doublet. Although it is not the most general case for PQ theories, it is always a good approximation in the so-called decoupling limit where the heavy Higgs bosons are much more massive than the weak scale. A threshold common to all axion models is the QCDPT where axion interactions with quarks and gluons become non-perturbative, strong interactions confine and one must resort to nonperturbative techniques. Recently, Ref. [52] provided the first smooth treatment of the QCDPT for the axion coupling to gluons in Eq. (1.1).
In this work, we provide predictions for the amount of axion dark radiation within UV complete models. Our methodology features two key steps. First, we evaluate the axion production rate at any temperature with smooth treatments of all mass thresholds. With the rate in hand, we solve the Boltzmann equation for the axion abundance and translate the resulting amount into a correspondent ∆N eff .
Notwithstanding the broad landscape of PQ theories [53,54], we can divide them into two main classes according to the origin of the color anomaly.
• KSVZ axion [55,56]. SM fields are PQ-neutral, and the color anomaly is due to a new heavy-colored fermion Ψ that gets mass from PQ breaking. Following the chronology of the expansion history, the first mass scale that we encounter is the fermion Ψ. Binary collisions involving Ψ are the main production channel above such a threshold. Gluon scatterings mediated by the operator in Eq. (1.1), which is generated once we integrate out Ψ, control the production rate at temperatures below the Ψ mass. These processes are the main production channel until we reach the scale where strong interactions confine. This is the second and the last threshold that we need to treat for this framework, and we connect the production rate above confinement with the one where the relevant degrees of freedom become hadrons.
• DFSZ axion [57,58]. SM quarks are responsible for the color anomaly and the Higgs sector is extended with another weak doublet. We work in the decoupling limit where the heavy Higgs bosons have a mass m A substantially larger than the weak scale. Such a high mass scale is the first threshold: we have a two Higgs doublet model (2HDM) above and the SM below, respectively. The second threshold is the EWPT, not present for the KSVZ case because the axion did not couple to the Higgs field; PQ charges of the Higgs and SM fermions make it relevant for this case. Finally, the QCDPT is also something to account for as we did for the previous case, although the details of the matching are different as a consequence of different axion couplings.
We compute the production rate at any temperature for the KSVZ and DFSZ axion in Secs. 2 and 3, respectively. We employ these results in Sec. 4 to quantify how many axions are produced thermally in the early universe and to predict the resulting contribution to ∆N eff . Sec. 5 contains our conclusions, and we defer all technical details to appendices.

The KSVZ Axion
The minimal ingredients for the KSVZ framework are an electroweak singlet complex scalar ϕ and a vector-like colored fermion Ψ. Electroweak charges for the fermion are allowed but not mandatory, and we work in the scenario where it is only charged under the fundamental of the SU (3) c gauge group. The Lagrangian for this case reads The theory features a global symmetry that acts on the fields as follows For any value of the transformation parameter α, the Lagrangian in Eq. (2.1) is invariant as long as the scalar potential V KSVZ (ϕ) does not change and the charges satisfy q ϕ = q R − q L . Furthermore, two crucial ingredients must be satisfied for this to be a viable PQ symmetry: broken in the vacuum state and anomalous under strong interactions. The potential where λ ϕ is the quartic self-coupling for the field and v ϕ its vacuum expectation value (vev), satisfies the first requirement. The condition on the anomaly is satisfied as long as the global symmetry is not vector-like, q L = q R (i.e., q ϕ = 0). Thus the phase of the complex scalar ϕ, which corresponds to the KSVZ axion, appears in the gluon anomaly operator in Eq. (1.1) and it eventually leads to a natural solution of the strong CP problem. The KSVZ axion originates solely from the phase of the complex scalar ϕ. As long as the quartic coupling is λ ϕ ∼ O(1), the radial mode of ϕ is rather heavy with a mass O(v ϕ ). We neglect fluctuations along the radial direction, and we identify the axion a as the phase of the complex field, ϕ → v ϕ / √ 2 e ia/vϕ . The Yukawa coupling is responsible for a fermion mass, m Ψ = y Ψ v ϕ / √ 2, which can be smaller than the symmetry breaking scale v ϕ if y Ψ is small, but not in conflict with collider searches for heavy colored states (m Ψ TeV).
The effective Lagrangian below the symmetry breaking scale reads We find this Lagrangian convenient to compute the axion production rate at temperatures above m Ψ . At lower temperatures, it is preferable to employ a different field basis. Let us describe in detail the difference between these two choices to exploit the interplay between the axion and PQ-charged fields. On one hand, the PQ symmetry can be linearly realized and the axion appears as the phase of the PQ breaking scalar as in Eq. (2.1). On the other hand, the PQ symmetry can be non-linearly realized and the axion shifts under a PQ transformation, a → a + const. The second option can be reached by performing the axion-dependent chiral rotation Ψ → exp i a 2vϕ γ 5 Ψ, and the resulting Lagrangian contains the changes at the classical level as well as the effects of the anomaly through Eq. (A.31) We introduce the axion decay constant f a and we set it to f a = v ϕ so we reproduce the normalization in Eq. (1.1). Although the field basis to describe axion interactions is not unique, the scattering cross sections calculated in App. C are independent on such a choice. The production of the KSVZ axion goes through three main cosmological phases that are separated by two mass thresholds: (i) the mass of Ψ; (ii) the confinement scale.

Matching at the heavy PQ fermion threshold
Above the heavy fermion mass m Ψ , axion production is driven by the scatterings where g is a gluon. Below m Ψ , the rate is controlled by quark (q) and gluon scatterings The long range nature of gluon interactions leads to IR divergences that need some care. The prescription to regularize such divergences of Ref. [59], which holds for soft external momenta (p g s T ), works only in the weak-coupling regime for axion production [37]. Once the QCD coupling g s gets stronger one needs to go beyond the hard thermal loop (HTL) approximation. Ref. [38] parameterized the rate in such a regime as follows 1 Here,c Ψ g (T ) denotes the effective gluon anomaly coefficient discussed in the next paragraph. The numerical factors are d g = 8 and ζ (3) 1.2 for the dimension of the SU (3) adjoint representation and the Riemann zeta function, respectively. The result for F 3 (T ) provided by Ref. [38] for temperatures well above the weak scale allows us to deal with the heavy PQ fermion threshold, but it is not enough to approach the QCDPT. We evaluate F 3 (T ) in App. B at any temperature in the QCD perturbative regime, and we keep into account the decoupling of heavy quarks. We evaluate the rate with the aid of the 'RunDec' [60] code that accounts for the running of the strong coupling constant α s up to four loops.
The UV origin of gluon scatterings is due to a heavy PQ-charged colored fermion. At temperatures much larger than its mass the effect is negligible, and it becomes relevant only once we integrate out the fermion for physical processes with typical energies smaller than the mass of the fermion itself. We can make this statement quantitative by evaluating the 1PI effective action. 2 For a generic colored fermion χ charged under PQ, we can perform a chiral rotation to induce the trilinear axion anomalous coupling to gluons Here, the Wilson coefficient c χ g is a constant number that depends on the quantum numbers of the fermion χ. However, if one cares about the axion production rate the relevant quantity to consider is the 1PI effective action that we parameterize as follows . KSVZ axion production rate across the heavy colored PQ fermion mass m Ψ = 10 5 GeV.
The black solid line indicates the total rate, and it is the sum of Ψ scatterings (solid blue) and thermal gluon scattering (solid red). We show for comparison the gluon scattering rate obtained with only the Wilsonian action contribution (dashed red).
Contrarily to the previous case, the effective couplingc χ g (q 2 ) depends on the momentum exchanged in the physical process under consideration. The relation between the Wilson coefficient and the 1PI effective coupling reads where we express the momentum dependence in terms of the dimensionless τ χ = q 2 /4m 2 χ . After this general discussion, we get back to the KSVZ axion (c Ψ g = 1) and set q 2 ≈ T 2 in Eq. (2.11). When Ψ is a degree of freedom of the thermal bath, with the universe much hotter than m Ψ , one-loop 1PI corrections effectively cancel out the Wilson coefficient leading The diminished gluon scatterings at T m Ψ ensures the dominance of Ψ collisions. At temperatures below m Ψ , the effective coefficient c Ψ g (T m Ψ ) becomes nearly unity and this is the radiative remnant of Ψ. We set m Ψ = 10 5 GeV (vertical green line), and we plot the combination γ a f 2 a as a function of the temperature in Fig. 2. Axion interactions below m Ψ mediated by the dimension 5 operator with gluons lead to the scaling γ a ∝ f −2 a . In the opposite regime, the axion interacts via a renormalizable Yukawa coupling y Ψ = √ 2m Ψ /f a , and the rate scales as γ a ∝ y 2 Ψ ∝ f −2 a also at high temperatures once we fix the fermion mass. The solid black line denotes the total rate that is the sum of Ψ (solid blue line) and gluon (solid red line) scatterings. In order to appreciate the difference between Wilsonian and 1PI descriptions, we also show the gluon scattering rate that we would get by accounting for the Wilsonian contribution only (dashed red line); this corresponds to settingc Ψ g (T ) = 1 in Eq. (2.8) at any temperature. Consistently with our picture, the solid and dashed red lines are in exact agreement for T < m Ψ , and they differ significantly for larger temperatures. In particular, γ w/o 1PI gg ∝ T 6 dominates even at large temperatures while the correct functional dependence for the production rate in the UV scales as γ Ψ ∝ T 4 .

Matching at the QCD threshold
The picture where the axion field interacts with quarks and gluons breaks down once we approach the scale Λ N ∼ 2 GeV and strong interactions become non-perturbative. Quarks are confined within hadrons at lower energies, and one must resort to non-perturbative techniques. Here, we employ the ones of chiral perturbation theory (ChPT) to compute the axion production rate from hadron collisions, and we provide a smooth result across the QCDPT. As we review in App. A, the correct prescription to determine axion couplings to hadrons is to match currents with the same symmetry properties between the UV and IR theories [62,63]. Such a procedure is straightforward within the KSVZ framework since the axion interacts with the strong sector only through the anomalous coupling to gluons. We find it convenient to rotate it away due to the large instanton effects at low energy, and this is done via the axion-dependent field redefinition of the light quarks The matching conditions provided in App. A allow us to find axion couplings to mesons (for coupling to baryons see Ref. [42]). Up to what UV cutoff Λ ChPT can we push ChPT? We treat the primordial plasma within the hadron resonance gas (HRG) approximation [64][65][66] which is inconsistent with lattice QCD results above T > 150 MeV [67]. Furthermore, leading order ChPT for axion production suffers perturbativity problems at T > 62 MeV [68]. Notwithstanding both values of Λ ChPT being smaller than Λ N , we find it plausible that the axion production rate between Λ ChPT and Λ N is connected smoothly since the QCDPT is a crossover where thermodynamic variables are continuos [69,70]. Such a connection for the anomalous coupling to gluons, and in particular for the KSVZ axion, was provided recently by Ref. [52].

Summary: production rate for the KSVZ axion
We summarize the KSVZ axion production rate in the whole temperature range in Fig. 4. Renormalizable interactions at large temperatures give γ Ψ ∝ T 4 , gluon scattering mediated by a dimension 5 operator below m ψ give γ gg ∝ T 6 . The rate drops below the QCDPT because of an exponential Maxwell-Boltzmann suppression for the pion number density.

The DFSZ Axion
A complex scalar singlet ϕ is also the starting point for the DFSZ framework. However, instead of introducing a new colored fermion we extend the SM Higgs sector with another weak doublet. We consider a two Higgs doublet model (2HDM) where the Higgs fields H u and H d carry opposite hypercharges ±1/2 and couple to up-and down-quarks, respectively. The essential 2HDM features needed for our discussion are summarized in App. A. The Lagrangian for scalars within the DFSZ framework takes the schematic form Unlike the previous case, SM fields are PQ-charged and in particular the combination H u H d carries a non-vanishing PQ charge. The specific scalar potential is model-dependent, but it must ensure the spontaneous breaking of two symmetries: PQ at the scale v ϕ , and electroweak at the scale v. The latter is due to the vevs of the Higgs field H u and H d that we parameterized as v u / √ 2 and v d / √ 2, respectively. Following a standard convention in the literature, we parameterize their ratio as tan β = v u /v d . Moreover, scalar potential interactions must couple the PQ breaking field ϕ with the two Higgs doublets to have a solution to the strong CP problem. Options for this latter constraint include the renormalizable coupling ϕ †2 H T u iσ 2 H d and the super-renormalizable coupling ϕ † H T u iσ 2 H d . We keep our discussion general and we parameterize this coupling as follows where we introduce the vev of the complex scalar field ϕ = v ϕ / √ 2. Upon appropriate field redefinitions, it is always possible to set the B parameter to be real and positive. The exponent r is connected to the so-called domain wall number via the relation N DW = 3r, and N DW corresponds to the number of degenerate vacua of the axion potential. Finally, SM fermions have the following Yukawa interactions with the Higgs fields whereH α = iσ 2 H * α and y (α) , with α = u, d, e, are Yukawa matrices for the type-II 2HDM. At energies below PQ breaking, the phase of the scalar field ϕ corresponds to the axion. As in the KSVZ scenario, there are the two typical ways to delineate the effective axion interactions to other fields. When the PQ symmetry is linearly-realized, ϕ → v ϕ / √ 2 e −ia/vϕ , the effective axion Lagrangian is given by On the contrary, we can realize the PQ symmetry non-linearly via axion-dependent field redefinitions ξ → exp iq ξ a vϕ ξ, where we rotate all fields ξ carrying a non-vanishing PQ charge q ξ . The PQ invariance of the scalar potential in Eq. (3.2) imposes the constraint q Hu + q H d = r = N DW /3. Likewise, the invariance of the Yukawa couplings in Eq. (3.3) imposes the relations among global charges: q Hu = −q Q L + q u R , q H d = −q Q L + q d R , and q H d = −q L L + q e R . After these rotations, we find the Lagrangian where f denotes SM fermions and we introduce the spin-one Higgs currents The anomaly coefficients after these rotations can be determined through the general result in Eq. (A.31), and they explicitly read As we discuss later, generic values of q Hu and q H d induce after electroweak symmetry breaking a mixing between the axion and the longitudinal Z weak gauge boson. Finally, we parameterize the axion anomalous coupling to photons below the weak scale in the standard form as follows E (a/v ϕ ) e 2 /32π 2 F µν F µν . We extract if from the couplings in Eq. (3.5) and we find E = (8/3)N DW .
We work in the so-called decoupling limit where the extra Higgs bosons are much heavier than the weak scale. Indeed, the 2HDM is phenomenologically constrained to be in such a region to respect LHC bounds [74][75][76][77]. There are three thresholds in this case. Two of them are analogous to the KSVZ scenario: the heavy Higgs bosons m A , and the QCD non-perturbative scale Λ N ∼ 2 GeV. An additional threshold is the EWPT.

Matching at the heavy Higgs bosons threshold
Above the EWPT, the axion field lives entirely inside the phase of ϕ. The linear realization of the PQ symmetry, with interactions as in Eq. (3.4), is the most convenient option to perform the rate calculation in this phase. The single axion coupling reads where the factor of (1/3) comes from the normalization of the Wilson coefficient of the gluon anomaly operator, v ϕ = N DW f a , to reproduce the convention in Eq. (1.1).
Above the mass scale m A , axion production is controlled by scatterings of Higgs bosons mediated by the interactions in Eq. (3.6). At temperatures below m A , the number density of heavy Higgs bosons gets Maxwell-Boltzmann suppressed and axion production is due to scatterings of SM particles (including the lighter Higgs doublet corresponding to the SM-like Higgs). The interactions mediating scatterings at low temperatures can be found by integrating out the heavy scalars with H SM the SM-like Higgs doublet. The Yukawa matrices y (α) for the 2HDM appearing in Eq. (3.3) and the correspondent Y (α) defined in Eq. (A.7) are related as follows The mixing angle α between the two doublets is a temperature dependent quantity and it is defined in Eq. (C.8). As the temperature drops below m A , thermal corrections to the Higgs mass matrix become sub-dominant with respect to the overall mass scale √ B. Hence the mixing angle α is approximated by β and the mass eigenstates coincide nearly with those in the vacuum defined in Eq. (A.23)-(A.25).
The interactions in Eq. (3.7) are equivalent, via appropriate field redefinitions, to the commonly used DFSZ axion interactions with SM fields parameterized as follows We point out how working with a linearly realized PQ symmetry and with axion interactions in Eq. (3.7) prevents any axion mixing with the Z boson. The lack of such a mixing, which for the nonlinear realization in Eq. (3.9) must be achieved by hand, is automatic with our procedure. The scatterings producing final state axions and their relative cross sections are summarized in App. C, and they lead to the rate shown in Fig. 5. Consistently with our choice to work in the decoupling limit, we set √ 2B = 10 5 GeV and the resulting heavy Higgs bosons mass is around the same scale. We visualize this mass threshold with a vertical green line. The total rate is given by the solid black line. At temperatures larger than m A , scatterings of heavy Higgs bosons dominate the total rate, and this is the contribution γ A that we denote with a solid magenta line. As expected, the magenta line drops exponentially at temperatures below m A . Scatterings of SM particles control axion production below m A . At temperatures much smaller than m A , the rate can be evaluated either with the interactions in Eq. (3.7) or the ones in Eq. (3.9). As explained, once the temperature is much smaller than m A the temperature dependent angle α reaches the constant value β and the two Lagrangians are equivalent. However, once we are not too far from m A , the correct prescription is to evaluate axion production via Eq. (3.7). For comparison, we report the rate computation obtained by using the non-linear realization in Eq. (3.9) at all temperatures (γ SM , dashed gray line). As expected, it agrees with the full result at temperatures below m A but it is substantially different at large temperatures.

Matching at the electroweak threshold
We work in the decoupling limit and therefore we can match across the electroweak threshold. Axion interactions in Eq. (3.9) are valid both above and below the EWPT and this is the field basis we employ to go across this threshold. We only consider the 2HDM parameter space with a smooth EWPT. If the mass of the heavy neutral pseudo-scalar m A defined in Eq. (A.24) is much heavier than the SM Higgs and Z boson, as it is the case for the decoupling limit, then the EWPT would be second order [78,79].
Similarly to the KSVZ scenario, the anomalous coupling to gluons mediates axion production, and the rate is given again by the expression in Eq. (2.8). The UV origin for this interaction in the KSVZ scenario was the Yukawa operator of the heavy colored fermion Ψ. On the contrary, for the DFSZ scenario this operator originates from the Yukawa operators of SM quarks. The couplingc q g (q 2 ) in the 1PI effective action receives threshold corrections from each quark as prescribed by Eq. (2.11). The low-energy remnant once we integrate each quark can be read off Eq. (3.9): c u g = cos 2 β/3 and c d g = sin 2 β/3 for the each family of up-type quarks and down-type quarks, respectively. The production rate through the gluon anomaly vanishes above the EWPT, and it subsequently becomes more and more significant due to the accumulated effective 1PI gluon coupling.
Furthremore, quark scatterings via the couplings given in Eq. (3.9) also contribute to axion production. Their cross sections are provided in App. C. The matching across the EWPT for this class of interactions was spelled out in Ref. [51]. Axion production via SM fermion scatterings requires a chirality flip. Above the EWPT, a chiral flipping can occur only via the Yukawa interactions in Eq. (3.7) so that only fermion scatterings with Figure 5. DFSZ axion production rate across the heavy Higgs bosons thresholds. We set tan β = 10 and √ 2B = 10 5 GeV m A . The solid black line denotes the total rate. The solid magenta line (γ A ) is the partial rate from scatterings of heavy Higgses. The dashed gray line (γ SM ) corresponds to the rate due to only SM particle scattering processes.
components of the Higgs doublet contribute to the axion production. On the contrary, after spontaneous electroweak symmetry breaking, the same Yukawa interactions provide quark masses that allow for chirality flips also for scatterings with gauge bosons, with gluons dominating the rate because of the hierarchy among the gauge coupling constants.

Matching at the QCD threshold
The procedure to investigate DFSZ axion production below the QCDPT is analogous to the one discussed in Sec. 2.2 for the KSVZ scenario. The leading order axion coupling to the strong sector in the KSVZ scenario originates only in the gluon anomalous term, whereas there are the additional axion interactions to the quark currents in the DFSZ scenario as given by Eq. (3.9). In other words, the effective axion interactions to the current of the light quarks (u, d, and s below Λ N ) can be written as Eq. (2.13) with the replaced coefficients (3.10) Figure 6. DFSZ axion production rate across the QCDPT. We set tan β = 10. QCD is perturbative for T > Λ N = 2 GeV and the rate is controlled by both gluon and quark scatterings. Right below confinement, but still above T ∼ 10 MeV, pion scatterings (γ ππ ) dominate. Below T 10 MeV, lepton scatterings (γ e,µ ) is the only production channel available. We interpolate for the two values Λ ChPT = 62 MeV (dashed red) and 150 MeV (dashed blue). The brown region identifies the EWPT.
Here, the first element for each coefficients comes from the gluon anomaly in common with the KSVZ scenario and the second one comes from the PQ charge of SM quarks.
Through the same matching procedure discussed in App. A, we find the effective axion couplings to hadrons. We report here the ones to pions, which dominate the rate, and they are still given by the operator in Eq. (2.14) but with the replaced coefficient (3.11) Fig. 6 shows the numerical result for the axion production rate across the QCDPT. The total rate is denoted by the solid black line. At temperatures right below the confinement scale, pion scatterings (γ ππ ) dominate axion production. As discussed already for the KSVZ scenario, this evaluation for the rate is trustworthy only up to the cutoff Λ ChPT , and we interpolate the axion production rate between Λ ChPT and Λ N = 2 GeV. The dashed red and dashed blue lines correspond to the interpolations for Λ ChPT = 62 MeV and 150 MeV, respectively. Unlike the KSVZ scenario, these two interpolations give slightly different results for the DFSZ case. We will discuss the impact of the interpolation on cosmological observables in the next section.
they do not carry color charge. The relevant scattering processes together with their cross sections are summarized in App. C. As shown in Fig. 6, when the universe cools down much below Λ ChPT (i.e., T m π ), the pion contribution to the axion production rate diminishes exponentially and lepton scatterings (γ e,µ ) become eventually dominant.
The bump arising near the EWPT (brown region) in the high temperature region of Fig. 6 is the combination of several effects. Below the EWPT, the production rate is the sum of two contributions: thermal gluon scatterings via the axion anomalous coupling with a rate γ gg ∝ T 6 , and bottom quark scatterings with gluons with a rate scaling as γ b ∝ T 4 (the different scaling is because the bottom mass provides the chirality flip). As we go above the weak scale, the thermal gluon scattering rate γ gg switches off, and top quark scatterings become available. However, we do not have fermions masses and therefore bottom and top quark scatterings (comparable since of tan β = 10) lead to the scaling as γ b,t ∝ T 6 .

Summary: production rate for the DFSZ axion
We summarize the production rate for the DFSZ axion in the whole temperature range in Fig. 7. Besides providing results for tan β = 10 (solid black), as done already for the previous figures, we show the rate also for tan β = 3 (dashed black). To ease the comparison, we set B to reproduce the same heavy Higgs boson mass in the two cases, Similarly to the KSVZ case, axion production is controlled by renormalizable interactions above the heaviest threshold, the heavy Higgs boson mass in the DFSZ scenario, and the rate consequently scales as γ A ∝ T 2 . At temperatures below m A but still above the EWPT, axion production processes proceed via dimension 5 operators coupling the axion to SM fermions (with top and bottom dominating), and therefore the rate scales as γ t,b ∝ T 6 . Axion couplings to top quarks exhibit the tan 2 β suppression discussed previously, and as a consequence the production rate with tan β = 3 at temperatures between m A and the EWPT is relatively larger than the one for tan β = 10. Below the EWPT, quark scatterings with gluons are dominant but top quarks quickly disappear from the bath, and the production rate is almost independent on tan β. Contrarily to the previous case, below the weak scale the SM fermion scattering rate scales as γ t,b ∝ T 4 , and in this region it dominates over the thermal gluon scattering that becomes active below the top quark mass with scaling γ gg ∝ T 6 . Pion scatterings dominate below the QCDPT and before we hit the Maxwell-Boltzmann suppression. Another important difference with respect to the KSVZ scenario is that the production is active even below the QCDPT since the axion couples to leptons, and interactions with muons and electrons give a rate with the scaling γ e,µ ∝ T 4 .

QCD Axion Dark Radiation
Scatterings of thermal bath particles produce axions in the final state, and the typical energy involved in each process is the bath temperature: the produced axions carry a kinetic energy much larger than their mass and therefore they are ultra-relativistic. What happens next? Initially, there are not enough axions to give the inverse (axion destruction) process and to ensure ultimately thermal equilibrium; axions just free-streams with their momentum decreasingly as the inverse scale factor. If axion production is efficient, we produce enough of them to thermalize with the primordial bath until the universe gets too cold and diluted, and they decouple with a relativistic thermal abundance exactly as neutrinos do.
A natural and useful application of the rates computed in Secs. 2 and 3 is keeping track of the axion abundance. Our conceptual starting point is an early universe going through an inflationary expansion, and inflaton decays generate the thermal bath afterward. Our only assumption is that inflationary reheating ends at high scales, and in particular the primordial thermal bath of relativistic SM particles dominates the energy budget earlier than the EWPT. However, we do not commit to any explicit hypothesis about axion production during inflationary reheating, but we consider two opposite cases in our analysis: we end reheating with no axions whatsoever, or we begin the radiation dominated era with a full thermal axion abundance. These two extremes cover all the options in between.
The quantitative tool to track the axion abundance is the Boltzmann equation Here, n a is the axion number density and t is the cosmic time. The number density dilution due to the Hubble expansion is accounted for by the term on the left-hand side proportional to the Hubble parameter H; in the absence of interactions, this is the only effect changing the axion number density. If number changing processes happen at an appreciable rate, we have to include their effects through the collision term C a on the right-hand side. We focus here on collisions producing one axion in the final state. For this class of processes, which is by far dominant as a consequence of the tiny couplings, the general expression for the collision term takes the form C a = γ a 1 − n a n eq a .

(4.2)
Here, the function γ a is the total axion production rate which is the sum of several contributions, one for each process that we account for. If we consider thermal gluon scatterings the rate is given by the expression in Eq. (2.8). For a generic binary collision with B i a bath particle (SM or beyond the SM), the associated rate reads Initial state bath particles are in thermal equilibrium, the scattering cross section is multiplied by the Moeller velocity, and the brackets denote a thermal average over all possible initial states. We use the Maxwell-Boltzmann statistics for the equilibrium density of bath particles since quantum degeneracy effects lead to negligible corrections (4.5) The particle mass and internal degrees of freedom are denoted by m i and g i , respectively, and the second kind modified Bessel functions are denoted by K i (x). The explicit expression for the thermally averaged cross section reads [80] The integral accounts for all the possible squared center of mass energies s in the collision with cross section σ ij→ka (s). The lower integration extreme corresponds to the kinematical threshold s min = Max (m i + m j ) 2 , m 2 k , and the Källén function λ is defined as follows It is convenient to rewrite the Boltzmann equation in terms of dimensionless quantities. We trade n a with the comoving number density Y a = n a /s R , where s R is the entropy density of the thermal bath. Other than being dimensionless, the comoving number density is advantageous because it scales out the effect of the Hubble expansion and therefore it varies only if number changing processes are in action. Likewise, we replace the time evolution variable t with the dimensionless inverse temperature x = M/T . Here, the choice for the scale M is purely conventional. Upon using the general result provided in Eq. (D.8) of App. D, which is a consequence of entropy conservation, we trade t with x and rewrite the Boltzmann equation in terms of dimensionless quantities (4.8) Our ultimate goal is to quantify how axions contribute to ∆N eff . Regardless of the details of axion production, there will be a point where Hubble expansion takes over the production, and this can happen because of two reasons. Particles participating in axion production can be massive, and as the bath temperature decreases number densities get exponentially suppressed. Even if production is mediated by massless particles, the universe gets too cold and diluted eventually to give appreciable interactions within a Hubble time. Such a freeze-out of interactions happens when the bath temperature was T F.O. , long before the CMB formation, and the axion comoving density freezes to a constant value Finding such an asymptotic value is the goal of our Boltzmann equation analysis. Once we have it, we evaluate ∆N eff via the general relation in Eq. (D.23) that for the axion reads Here, g SM * s (T CMB ) is the SM contribution to the effective number of entropic degrees of freedom, and the second term in the denominator accounts for the axion contribution to the energy density. As explained in App. D, this correction can be at most 1/(1+g SM * s (T F.O. )) and therefore becomes more relevant for late axion production. The SM contribution cannot be less than approximately 4, hence the correction can be at most 25%. However, axions are produced well above the MeV scale for most of the parameter space we explore, and our complete results are well described by the approximated expression (4.11)

KSVZ axion
The production rate for the KSVZ axion, with a smooth treatment of the heavy colored PQ fermion and the QCDPT thresholds, is shown in Fig. 4 as a function of the temperature. We feed the Boltzmann equation in Eq. (4.8) with this rate, and we show results from numerical integrations in Fig. 8. We set the mass of the heavy PQ fermion to m Ψ = 10 5 GeV, as done already in Sec. 2, and we run our code starting from an initial temperature T i = 10 7 GeV. Thus we always go across the Ψ threshold. Furthermore, we set n a (T i ) = 0 as the initial condition to produce this figure, and this choice has no impact on our final results as we discuss below. We employ the dimensionless combination x = M/T as the evolution variable, and we set M = 1 GeV to have the QCDPT around the region x 1. Different colors denote different choices for the axion decay constant f a which is the only free parameter within the KSVZ framework. As a matter of fact, each solution is not a line but rather a band whose width is due to the two different datasets [32,33] that we employ for the temperature evolution of the effective relativistic bath degrees of freedom. Finally, the dashed black line denotes the axion equilibrium comoving density whose The three colors for the solid lines correspond to different axion decay constants f a . Each solution is presented as a band whose width corresponds to two different treatments of the thermal bath relativistic degrees of freedom. The dashed black line describes the equilibrium distribution. The vertical axis on the right identifies the contribution to ∆N eff .
analytical expression is rather simple, Y eq a (x) = 45 ζ(3)/(2π 4 g * s (x)), and whose temperature dependence is only due to the change in the effective entropic degrees of freedom of the thermal bath. The three solid lines reach the thermal equilibrium distribution rather quickly, as early as T 10 6 GeV, and therefore setting the initial condition n a (T i ) = 0 has no impact on the final results for the values of f a chosen in Fig. 8.
At small enough temperatures (i.e., large x), the number density reaches an asymptotic value in agreement with Eq. (4.9). We include a second vertical axis on the right of Fig. 8 to identify the corresponding value of ∆N eff as quantified by Eq. (4.10). This contribution increases as we go to lower values of f a , consistently with the picture that larger axion couplings keep physical processes efficient at lower temperatures (see Eq. (1.3) and Fig. 1).
Our choices for the axion decay constant correspond to generating a ∆N eff equal to the Planck bound at 2σ (red), 2σ and 1σ for CMB-S4 surveys (green and blue, respectively). Astrophysical constraints bound the axion decay constant from below [81][82][83]. We impose the bound from SN1987A provided by the recent Ref. [84] that for the KSVZ axion results in f a 1.4 × 10 8 GeV. Neutron star cooling provides bounds in the same ballpark [85,86]. Thus current Planck bounds on ∆N eff are sentitive to KSVZ axions with f a one order of magnitude below the stellar exclusion bound, and future CMB-S4 surveys will probe the range f a ∼ (10 8 , 10 9 ) GeV that is still not in conflict with any experimental constraint.
We investigate how ∆N eff depends on f a in Fig. 9. We solve the Boltzmann equation again with initial condition n a (T i ) = 0, but we consider a few different values for T i corre- sponding to the different solid colored band (whose width quantifies our uncertainties due to different treatments of the bath). The red line shows ∆N eff as a function of f a for any value of T i much larger than the PQ fermion mass that we keep m Ψ = 10 5 GeV as in the previous plots. Axion production at temperatures above m Ψ is controlled by a renormalizable coupling, and the rate normalized by the Ψ number density scales as γ a /n Ψ ∝ T as long as Ψ are relativistic. This has to be compared with the Hubble expansion rate that scales as H T 2 /M Pl . Thus axion production is most efficient at low temperatures, and most axions coming from Ψ scatterings are created at temperatures around m Ψ . 3 On the contrary, below the heavy PQ fermion we have a rate normalized by the bath number density scaling as γ a /n B ∝ T 3 , and this temperature behavior is stronger than the one for the Hubble rate: axion production is most efficient in the UV at the highest temperature available T i . This explains the different results for large values of f a : axions do not have enough interaction strength to thermalize in the early universe, and smaller initial temperatures lead to smaller ∆N eff because at low temperatures the production is less efficient. At low enough values of the axion decay constant thermalization is achieved, and all colored lines coincide: the resulting prediction for ∆N eff does not depend on T i .
We provide in Fig. 9 also the prediction for ∆N eff once we set the thermal equilibrium distribution as initial condition for the axion number density at temperatures above the weak scale (dashed gray line). Large axion decay constants, f a 3 × 10 9 GeV, lead to ∆N eff 0.027 which is the value associated with a spin-0 particle that was once in thermal equilibrium and decoupled above the weak scale (green line in Fig. 1 at large T D ). In this range of f a , we set the initial abundance to the equilibrium value by hand and interactions are completely harmless. Things are different as we approach lower f a since couplings get stronger and they can keep the axion in equilibrium below the weak scale. The resulting prediction for ∆N eff coincides with the solid red band for f a 3×10 9 GeV regardless of the initial value of T i as long as we keep T i 1 TeV (for lower values of T i the expected ∆N eff would be larger, see Eq. (1.3)). This result can be understood from the plot in Fig. 8: for the axion decay constant range we are interested in, axions always reach equilibrium long before the time when the bath temperature gets to the TeV scale. Thus in the physical region of our interest where the signal is detectable, 10 8 GeV f a 10 9 GeV, our predictions for ∆N eff do not depend on the initial condition for the axion number density.

DFSZ axion
We now turn to the DFSZ framework. The production rate, this time with smooth treatments of three different mass thresholds, is shown in Fig. 7 as a function of the temperature. Exactly as we just did for the KSVZ axion, we feed the Boltzmann equation in Eq. (4.8) with this rate, find the asymptotic value of the axion comoving number density and quantity the correspondent contribution to ∆N eff . We fix the model parameters to the same values as in Sec. 3, tan β = 10 and m A √ 2B = 10 5 GeV, and we run our Boltzmann code again starting from an initial temperature T i = 10 7 GeV; this ensures that we pass again all the mass thresholds in the scenario under investigation. The numerical output of the differential equation integrations is shown in Fig. 10. The axion number density reaches Figure 11. Contribution to ∆N eff as a function of the axion decay constant f a for the DFSZ axion.
its equilibrium value rather quickly, and therefore setting its initial value to zero does not impact our final results. At low temperatures, consistently with our previous discussion, the comoving number density settles to a constant value. We choose again three numerical values for the axion decay constant leading to ∆N eff equal to the Planck bound and the projected sensitivities of future CMB experiments. They are in the same ballpark as the ones for the KSVZ scenario.
Astrophysical bounds are more severe for this case. Data from SN 1987A [83] constrain again the axion decay constant, but this time the numerical value associated to the bound is slightly different because the DFSZ axion couples also to quarks. For the tan β chosen in this analysis, we find f a 1.9 × 10 8 GeV. However, this is not the leading bound since the DFSZ axion couples to electrons as well. Studies of red giants [88] and white dwarfs [89] provide competitive bounds, with the one coming from the latter slightly stronger. For the DFSZ parameter chosen in our analysis, this corresponds to the bound on the axion decay constant f a 5.2 × 10 8 GeV. Thus future CMB-S4 surveys will probe a rather small region of the DFSZ parameter space that is still not excluded.
We quantify the last statement in Fig. 11 where we explore how ∆N eff depends on the axion decay constant. Solid colored bands provide the prediction obtained with vanishing initial axion abundance at various initial temperatures T i . For comparison, we report also the prediction for the case when we begin the Boltzmann equation evolution with axions already in thermal equilibrium. As it was the case for the KSVZ axion, the predictions Finally, we explore how our predictions depend on tan β in Fig. 12 where we show ∆N eff as a function of f a for tan β = 3. We set the heavy Higgs boson mass m A to the same value as for Fig. 11, and we update stellar bounds consistently to account for the different axion couplings. As explained in Sec. 3, the production rates for tan β = 3 and 10 differ only above the EWPT with the former enhanced by approximately one order of magnitude. If axions are in thermal equilibrium until the time when the bath temperature is of the order of the weak scale there will be no difference between the two cases. This is manifest from a comparison between Figs. 11 and 12: the resulting hot DFSZ axion abundances with f a O(10 9 ) GeV are very similar for tan β = 10 and 3. However, we notice an effect at larger values of f a where it takes more effort for the axion to thermalize, and the predicted amount of dark radiation is enhanced for tan β = 3. Thus ∆N eff does depend on tan β at large values of f a , and smaller tan β makes the signal detectable for a wider range of axion decay constant values.

Interplay with inflationary reheating
Our analysis so far relied upon the assumption that the energy density of the universe was dominated by a gas of relativistic particles at the time of axion production. This is the extrapolation of how we "look at" our universe at the time of BBN, and it is worth keeping in mind that it is an extrapolation not supported by any observation. Within the inflationary paradigm, this extrapolation must come to an end because back enough in time the energy budget was controlled by the vacuum energy of the inflaton field. Inflaton decays populate the thermal bath eventually, and the highest temperature T R ever achieved during the radiation dominated epoch is known as the reheating temperature. Thus throughout our numerical Boltzmann analysis we have always implicitly assumed the hierarchy T R > T i .
Our predictions are insensitive to the dynamics of inflation as long as the reheating temperature is high enough. For both of the frameworks under investigation, axion production is mediated by renormalizable couplings above the highest mass threshold; we have a Yukawa interactions with the Ψ fermion and a scalar potential cubic term with the Higgs fields H u and H d in the KSVZ and DFSZ framework, respectively. The renormalizability of the axion couplings ensures that production is most efficient in the IR, and therefore around the mass of the heavy particles. Thus all we need is a reheat temperature larger than the heavy thresholds: T R > m Ψ and T R > m A for the KSVZ and the DFSZ framework, respectively. As a matter of fact, inflaton decays can provide an additional source for axion dark radiation. However, we consider also this option since in our analysis we accounted for the two extreme situations, n a (T i ) = 0 and n a (T i ) = n eq a (T i ), and this ensures that we cover all possible options. As we have already explained, our final predictions in the parameter space region where the signal is detectable do not depend on the initial conditions.
What happens for lower values of T R ? Inflationary dynamics can play a relevant role only for T R < m Ψ or T R < m A . In this regime, we solve the coupled Boltzmann equations describing inflationary reheating 4 This system describes the evolution of the inflaton (φ) and the radiation bath (R) energy densities. Inflaton decays, with a rate Γ φ , deplete the former and enhance the latter. The Hubble expansion rate, which allows us to understand when axions are produced most efficiently once we compare it with the production rate, is given by the Friedmann equation and it has the following scaling (4.14) Below the reheat temperature T R we have the typical scaling for a radiation dominated epoch. The reheat temperature is connected to the inflaton decay width through the relation Although this is defined as the highest temperature ever achieved by the thermal bath during the radiation dominated epoch, this is not the highest temperature achieved by the thermal bath in general. The bath itself exists even for temperatures larger than T R as a sub-dominant component since the decaying inflaton is still dominating the energy budget. The highest temperature ever achieved is usually denoted by T MAX and it scales as T MAX M Pl Γ φ E 2 I 1/4 , with E 4 I the constant energy density driving the inflationary expansion. Thus the bath temperature spans a potentially large range between T MAX and T R before becoming the dominant energy component. During this phase, which corresponds to an early matter domination with the energy budget controlled by inflaton oscillations, the Hubble rate is proportional to T 4 .
If we consider renormalizable interactions, production of particles via scatterings is efficient at low temperatures for a radiation dominated universe. This is the case also during inflationary reheating given the higher power of the temperature appearing in the scaling for the Hubble parameter. One can be more quantitative and state that in the range 5 ≤ d ≤ 8, with d the mass dimension of the operator mediating scatterings, the production is dominated at small temperatures during reheating and therefore maximized at T R [91][92][93][94]. In particular, for d < 8 the final abundance does not depend on T MAX . Only for d > 8 particle production is UV dominated also during inflationary reheating, and therefore the resulting abundance is sensitive to T MAX and to the details of reheating.
For the axion frameworks studied in this work, we never go above mass dimension 5 and therefore we are never sensitive to the peculiarities of inflationary reheating. Whether we are below the mass of Ψ for the KSVZ or the mass of A for the DFSZ, axion interactions are mediated by dimension 5 effective operators and the production is maximed at T R . This is the scenario investigated by Ref. [38] with the gluon and top quark couplings dominating the production rate for the KSVZ and DFSZ axion, respectively.

Conclusions
The PQ mechanism, where the θ parameter of QCD is promoted to a dynamical field, is undeniably one of the most elegant solutions to the strong CP problem. A plenitude of UV complete candidate models provides viable realizations of PQ symmetry breaking, but they all share the same low-energy residual: an approximate Nambu-Goldstone boson. Such a field, known as the axion, features the model-independent coupling to gluons given in Eq. (1.1) as well as model-dependent interactions with other SM particles. Given the Nambu-Goldstone nature of the axion, its couplings to visible matter are suppressed by the large PQ breaking scale and this makes axion detection very challenging. Notwithstanding these difficulties, the field of axion experimental searches has been literally blossoming in the recent decade, and the present time is rather unique for the quest for axions.
In spite of the rich set of options for axion couplings, all the terrestrial searches are sensitive to a handful of them: the ones to light quarks and gluons that in turn describe coupling to nuclei, the one to electrons, and the one to photons. An effective low-energy theory with only these interactions, with any UV completion that can be matched onto it, is enough to capture the phenomenology of axion searches.
We focused on an experimental signature to which all axion couplings can potentially contribute. The physics is the one of thermal axion production in the early universe, and the experimental manifestation is the presence of additional radiation that we infer from the CMB anisotropy spectrum. The net signal is accumulated through the expansion history with axions produced from bath particle collisions possibly at all temperature scales. Trustworthy predictions are possible only upon knowing axion couplings to all SM particles and to the model-dependent beyond the SM degrees of freedom specific to each theory.
The effect we consider is quantified by an additional contribution to the effective number of neutrinos species ∆N eff . Bounds on the amount of axion dark radiation from the Planck data are already quite remarkable, and prospects provided by CMB-S4 surveys make this signal rather intriguing for the future. This population of relativistic axions can also leave an imprint on cosmological structure through baryon acoustic oscillations (BAO), and this effect provides an additional constraint on ∆N eff [95][96][97]. Future large scale structure surveys will provide an improved BAO measurement, and this determination will be complementary to the CMB anisotropy spectrum. Furthermore, Ref. [98] suggested recently how such a cosmic axion background could be detected even with experiments in our terrestrial laboratories, although the signal is more sensitive to non-thermal energy spectra. These impressive projections combined with the top-down motivation make reliable theoretical predictions for the amount of axion dark radiation of the utmost importance.
This work addressed the presence of mass thresholds through the expansion history. As we showed with two explicit examples, the KSVZ and the DFSZ frameworks, passing through them alters the production rate significantly. For the KSVZ scenario, the rate changes significantly across the heavy fermion mass because renormalizable axion interactions become non-renormalizable, and this changes the rate temperature dependence. An analogous threshold is due to the heavy Higgs bosons in the DFSZ scenario, and in such a case the EWPT is also an important threshold where axion production mediated by fermion scattering changes drastically its temperature behavior as a consequence of chirality flips induced by fermion masses [51]. Finally, the QCDPT is common to both frameworks, and the matching procedure is far from being straightforward [52].
The central results of our analysis are predictions for ∆N eff as a function of the axion decay constant f a . They are explicitly presented in Fig. 9 for the KSVZ framework, and in Figs. 11 and 12 for the DFSZ framework. Scatterings of thermal bath particles lead to a detectable signal in the future for both cases. For the KSVZ axion, the less severe stellar bounds allow for a stronger signal, as large as the 2σ sensitivity of future CMB-S4 surveys and for values of the axion decay constant f a 3 × 10 9 GeV. On the contrary, white dwarf bounds for the DFSZ axion allow for a signal only detectable at 1σ in the future. The range of testable axion decay constants for the KSVZ axion depends on the specific value of tan β. For the two representative cases we analyzed, tan β = 10 and 3, the signal is testable for axion decay constants satisfying the upper bounds f a 2 × 10 9 GeV and f a 6 × 10 9 GeV, respectively.
We assumed a radiation dominated universe through our analysis, but the production rates in Figs. 4 and 7 are independent of the cosmological history. They can be employed to investigate axion production for alternative scenarios such as late inflationary reheating.
Our methodology can also be extended to other microscopic realizations besides the two frameworks studied here. Flavor-violating axion couplings, with the production rate controlled by decays of bath particles instead of scatterings, are of particular interest. Plausible origins for the flavor violation can be loop corrections to axion couplings [99][100][101][102][103], or they can even be present at tree-level as a consequence of the PQ charge assignments [104,105]. Predicting ∆N eff for specific axion UV complete models, along the lines of the analysis presented here, would be a piece of useful information to discriminate among them.

A Conventions and Useful Results I: Particle Physics
We collect in this Appendix our notations and conventions for SM fields and Lagrangian. Starting from the electroweak symmetric phase, we describe the SM gauge group and the matter content. We present spectrum and interactions, above and below the weak scale, for the case of a minimal Higgs sector with just one scalar weak doublet. The DFSZ framework features two Higgs doublets, and we discuss spectrum and interactions for this case as well. We quantify the effect of anomalous chiral rotations that are necessary to perform changes of field basis. Finally, we provide basic notions of ChPT and we introduce the formalism that, among several applications, allows us to determine axion couplings to hadrons.

Standard Model with minimal Higgs sector
At energies above the Fermi scale, the theory enjoys a full SU (3) c × SU (2) L × U (1) Y gauge symmetry, and the Lagrangian takes the form L SM = L gauge + L fermion + L Higgs + L Yukawa . (A.1) The first term contains gauge boson kinetic terms and they are constructed by employing the field strengths defined as follows with f ABC and IJK the structure constants of the non-Abelian groups SU (3) c and SU (2) L , respectively. The indices A = 1, . . . , 8 and I = 1, 2, 3 run over the adjoint representations. Fermion and scalar kinetic terms are expressed via the gauge covariant derivative The non-Abelian generators are Gell-Mann (λ A ) and Pauli (σ I ) matrices, and they act on color and weak-isospin indices (if any), respectively. The Abelian part has instead a term proportional to the hypercharge Y of the field the covariant derivative acts on. The SM matter fields with their quantum numbers are listed in Tab. 1. Upper case fermions denote weak doublets whereas lower case fermions are singlet under the weak-isospin group, and the index i runs over the three fermion generations. Fermion kinetic terms read We use a compact notation where the sum over the three different generations is understood. The Lagrangian for the Higgs field has a canonically normalized kinetic term and the most general renormalizable scalar potential leading to electroweak symmetry breaking Finally, the Yukawa part of the Lagrangian (also in a compact matrix form) reads where H = iσ 2 H * . In the most general fermion basis, Y (u,d,e) are generic 3 × 3 matrices in flavor space, and they can be diagonalized upon performing bi-unitary rotations Here, ψ = u, d, e and the matrices Y (ψ) are diagonal. Hence we redefine fermion fields by performing the following unitary rotations in flavor space After these operations, the Yukawa Lagrangian in Eq. (A.7) takes the form

Electroweak symmetry breaking
The scalar potential in Eq. (A.6) has an electroweak symmetry breaking minimum where the Higgs field gets a vacuum expectation value (vev) Table 1. SM matter fields and gauge quantum numbers (i runs over the three generations).
As a consequence, electroweak gauge bosons acquire mass terms and mix among each other. The complete gauge boson mass spectrum reads where the photon remains massless and we define the weak mixing angle (s w , c w ) = 1 .
Here, we introduce the Pauli matrices σ ± ≡ σ 1 ± iσ 2 that are used to define the ladder operators for the weak isospin group. Furthermore, we define (minus) the electron charge e ≡ gg /g w > 0 and the electric charge generator Q ≡ σ 3 /2 + Y . The same covariant derivative describes gauge interactions for the radial model of the Higgs field, the Higgs boson h, and this is the only scalar appearing in the Yukawa interactions with fermions.

Two Higgs Doublet Model
The Lagrangian for a theory with two weak doublets H u and H d is richer than the one given in Eq. (A.6), and it takes the schematic form We consider the renormalizable and gauge invariant scalar potential (A. 16) The only operator that is not invariant under PQ is the quadratic one proportional to B.
Within the DFSZ framework, it arises once the PQ breaking scalar gets a vev, and the phase of B contains the axion field above electroweak symmetry breaking. Upon redefining the Higgs doublets, it is possible to take the PQ-breaking coefficient B to be real and positive.
We search for an EWSB minimum where only the neutral components acquire vevs The vevs satisfy the constraint v 2 = v 2 u + v 2 d = (246 GeV) 2 , and we define the angle β as tan β ≡ v u /v d . The minimum conditions are We expand around this EWSB vacuum and we determine the mass spectrum. The Higgs doublets can be decomposed in terms of neutral and charged scalar components For neutral fields, we distinguish between CP-even and CP-odd scalars. The squared mass matrices for the charged, pseudo-scalar and scalar fields read respectively (A.20) The first two matrices have vanishing determinant, and this ensures massless Goldstones to provide longitudinal components for the W and Z gauge bosons. The masses of the heavy charged (H ± ) and pseudo-scalar (A) Higgs bosons can be found from the trace . (A.21) The CP-even Higgs bosons, the SM-like h and the heavier H, are both massive. We provide here the mass eigenvalues in the decoupling limit which is valid when B v 2 The charged mass eigenstates result in where G + is the charged Goldstone eaten up by the W boson. Likewise, if we dub G 0 the Goldstone eaten by the Z boson, we have the pseudo-scalar mass eigenstates Finally, in the decoupling limit (B v 2 ), we have the CP-even mass eigenstates

Anomalous chiral rotations
Chiral, axion dependent, rotations on fermion fields can be useful to find a new field basis better suited for the specific framework under investigation. Here, we state our conventions for these field transformations and we quantify gauge anomaly effects. For a generic Dirac fermion χ, before accounting for any interaction, the theory has a U (1) V symmetry where both left-and right-handed Weyl components are rotated with the same phase, and we have a conserved Noether's vector current J (V ) µ =χγ µ χ. If the fermion is massless (m χ = 0) then the theory has also the U (1) A symmetry, where left-and right-handed Weyl components are rotated with opposite phases, and the resulting Noether's current results in J On the contrary, if the fermion field is massive, we have that the axial current has a nonvanishing divergence The result above is valid at the classical level. Once we include quantum corrections, the axial current can be non-conserved if the fermion carries gauge charges, even if the fermion itself is massless. It is possible to derive its divergence by different methods, like evaluating the Green function of the axial current with gauge bosons in perturbation theory (triangle diagrams) [106,107], or via the Jacobian of the path integral measure [108].
We start from the well-known QED result The generalization to a non-Abelian gauge theory, such as QCD, does not require any new calculation, all we need to do is adding a group theory factor to the QED expression in Eq. (A.27) (the triangle diagrams have the same Lorentz structure) where t C are generators of the color group normalized as Tr t C t D = δ CD /2. We have all the tools to quantify the effects of performing chiral rotations. For a generic massless Dirac fermion χ, charged under a representation R of the SU (3) c gauge group, we perform the local axial rotation As a consequence, the Lagrangian changes due to both classical and quantum effects. The former is straightforward whereas the latter can be found by computing the change in the path integral measure [108]. We do not need to reproduce the derivation since we know that this contribution must reproduce the equation of motion in Eq. (A.28). Thus we have Upon comparing we find how chiral rotations alter the Lagrangian at the quantum level This is valid for an axial rotation of a Dirac fermion as in Eq. (A.29). If we only rotate one Weyl component the result is half the one above and with the appropriate sign.

Rudimental ChPT
We review basic notions of Chiral Perturbation Theory (ChPT) [109][110][111] needed to study axion couplings. Our starting point is the QCD Lagrangian with N f = 3 quark flavors For the ease of notation, we introduce the quark vector in flavor space q = (u d s) T where each entry is a Dirac field with left and right-handed components. Chiral projectors, defined in the usual way P L,R = (1 ∓ γ 5 )/2, extract the different quark chiralities: q L,R = P L,R q.
The quark mass matrix in the mass eigenbasis reads M q = diag (m u , m d , m s ).
If we neglect the quark mass matrix, left-and right-handed quarks are decoupled and the QCD Lagrangian in Eq. (A.32) is invariant under independent rotations of the two fermion chiralities: the theory has a global U (3) L × U (3) R symmetry. The vectorial part of the symmetry group where both chiralities are rotated by the same angle, namely the baryon number U (1) V and the isospin SU (3) V , are good approximate symmetries of Nature; the former is broken by gauge anomalies whereas the latter is only broken by quark mass differences and electroweak interactions. The axial part U (1) A × SU (3) A is spontaneously broken by the quark condensate, and we do not expect mixed parity multiplets in the hadronic spectrum but rather Goldstone bosons associated to the broken axial generators. However, there is no Goldstone boson associated to the broken U (1) A symmetry since m η m π . This was dubbed as the U (1) A problem of QCD [112] and it was solved only thanks to a complete understanding of the rich structure of the QCD vacuum [113,114]. We do not have a ninth Goldstone boson η in the spectrum because the U (1) A is not even an approximate symmetry of the QCD Lagrangian in the massless limit since it is anomalous.
Concerning the spontaneous breaking of the SU (3) A part, we employ a non-linear sigma model to describe the associated Goldstone bosons. Keeping only terms with two derivatives, which correspond to quadratic terms in the exchanged momentum O(p 2 ), we have the low-energy chiral Lagrangian where λ a are the SU (3) Gell-Mann matrices and f π 93 MeV. The Goldstone bosons π a , with a = 1, . . . , 8, enter through the unitary matrix U that under a generic chiral rotation transforms as U → L U R † . We pick the basis This choice is convenient since these fields are physical eigenstates once we introduce chiral symmetry breaking quark masses. The Lagrangian in Eq. (A.33) holds for exact chiral symmetry in the high-energy theory. Quark masses, which break chiral symmetry, are easily incorporated by applying the formalism of Refs. [110,111] for matrix elements of currents in the chiral effective theory. Furthermore, this method also allows us to derive axion couplings to the Goldstone octet in Eq.(A.34). We review this method starting from the QCD Lagrangian written as follows Here, we include four different external sources S: scalar s, pseudo-scalar p, vector left l µ , and vector right r µ . The QCD Lagrangian with quark mass terms in Eq. (A.32) is recovered for l µ = r µ = p = 0 and s = M q . Setting the spin-one sources to a non-vanishing value allows us to deal with coupling to vector bosons, such as the photon, as well as the spin-one axion currents. Finally, the pseudo-scalar current p also plays an important role to determine axion couplings. The Lagrangian in Eq. (A.35) has a local SU (2) L × SU (2) R symmetry where left-and right-handed quarks transform separately 36) and the sources also transform as follows (A.37) In the above equations, we restore the explicit dependence on the space-time location x to emphasize that the transformation is local. The spin-one sources transform as gauge fields so we define the covariant derivative and we can match the Lagrangian in Eq. (A.35) onto a low-energy chiral Lagrangian Here, µ is a dimensionful parameter that we determine by the requirement of reproducing the meson masses If we look at the pion mass, and we neglect the mixing with the η, we find m 2 π = µ(m u +m d ).

B Conventions and Useful Results II: Thermal Corrections
We discuss thermal corrections in this Appendix. We analyze axion production via thermal gluon scattering, and we compute thermal masses for the DFSZ Higgs sector.

Thermal corrections to axion production via gluon loops
The axion production rate can be expressed in terms of the axion self-energy [115,116] Here, f BE (E a ) is the Bose-Einstein distribution and Π < a denotes the non time-ordered axion two-point function. At finite temperature, the expansion parameter is not anymore α s /(4π) but rather the gauge coupling constant g s because of collinear enhancements [59], and this require in principle the resummation of infinite processes involving many particles. However, as explained by Ref. [38], such an enhancement is absent for the axion anomalous interaction given in Eq. (1.1), and we can safely consider only binary scatterings. Nevertheless, even if one restricts to binary collisions there are still IR divergences to take care of. The production rate with IR divergence properly accounted for can be found in Ref. [38] but only at high temperatures, T 10 4 GeV. We extend this study to lower temperatures.
Only one diagram contributes to the production rate, the one-loop axion two-point functions with virtual gluons. Within the context of the optical theorem, this contribution could be interpreted as the thermal gluon decay [38]; using the cutting rules, we could identify diagrammatically as g th → g th + a where g th indicates the thermal excitation of the gluon field in a medium. This is the reason why Ref. [38] dubbed it the 'decay" diagram. We evaluate the it with the resummed thermal gluon propagators in the loop, and we work at the leading order in g s because we can restrict to binary collisions. At this point, all we need to derive a proper thermal gluon propagator.
For a generic gauge theory, thermal effects induce the following correction to the gluon two-point function in momentum space The thermal self energy Π µν is a function of the external momentum K µ = (ω, kk). Here and thereafter, capital characters denote four-vectors whereas lower-case characters their components. The four-vector K has a spatial component of size k and directionk.
There are two 5 sources for the gluon self-energy at one-loop: gluon self-interactions and gauge interactions of colored fermions (i.e., quarks) The pure gauge contribution reads where C 2 (G) denotes the quadratic Casimir operator for the adjoint representation and f B.E. (p) is the Bose-Einstein distribution. The contribution from a quark q with mass m q results in where E p = p 2 + m 2 q , and f FD (E p ) andf FD (E p ) are the Fermi-Dirac distributions for the quark q and the anti-quarkq, respectively.
We decompose Π µν into its longitudinal (L) and transverse (T) components [118,119] The longitudinal and transverse gauge contributions read where we set E p = p and with ω ± ≡ (ω ± k)/2. We assume a negligible particle/antiparticle asymmetry for quarks, and consistently we set f FD (E p ) =f FD (E p ). We find where For a massless quark (m 2 q = 0), we recover the analytic expressions in the literature [38]. With the gluon self-energy at our disposal, we can evaluate the spectral densities The imaginary part must be extracted according to the rescription ω → ω +i0 + . In order to deal with technical difficulties in the numerical analysis, we take into account the spectral densities rewritten as where i = (T, L) and Z i indicates the residues at the poles of ω = ±ω i (k) which are located in the time-like region (|ω| > k), and we consider the contribution from the continuum parts ρ cont i only in the space-like region (|ω| < k) [117]. We express the axion two-point function in terms of spectral densities where d g = 8 is the dimension of the SU (3) strong gauge group and the four momentum of the gluons are K µ = (k 0 , kk) and Q µ = P µ a − K µ = (q 0 , qq). We integrate numerically Eq. (B.1) where we use Eq. (B.19) for the axion self-energy. We employ the 'RunDec' [60] code to account for the running of the strong coupling constant up to four loops. The numerical result of the control function F 3 defined in Eq. (2.8) is shown in Fig. 13; the left and right plots illustrate the value of F 3 in the function of temperature T and the strong coupling g s , respectively, including the decoupling of quarks at different temperatures. Figure 13. The control function F 3 defined in Eq. (2.8) as a function of the temperature (left) and the strong coupling g s (right). The dashed gray line on the right panel corresponds to the Hard Thermal Loop (HTL) approximation [37,59], and it is valid only for small couplings g s < 0.5.

Thermal masses for the electroweak sector
The tree-level scalar potential for the 2HDM is given in Eq. (A.16). The thermal evolution of electroweak sector in the 2HDM is the subject of Refs. [120][121][122][123], and one-loop thermal corrections to the potential read Here, we include contributions from bosons and fermions, and we denote their (Higgs fields dependent) masses m b (H i ) and m f (H i ), respectively. The relative minus sign is due to fermion fields running in the loop. The dimensionless quantities n b and n f count the number of internal degrees of freedom. The function J F and J B are defined as follows [118,119] At high temperatures, more specifically for m b,f /T < 1.8 [124], we can approximate where a f = π 2 exp(3/2 − 2γ E ) and a b = 16a f with the Euler-Mascheroni number γ E . In the electroweak symmetric phase, the global minimum is located with the Higgs fields at the origin, v u,d (T > T EWPT ) = 0. The thermal corrections to the quadratic fluctuations of the Higgs fields around the origin read The factors of cos β or sin β come from the Yukawa interactions to reproduce the SM fermion spectrum. One can obtain the thermally corrected mass matrices for the Higgs sector by denotes the zero-temperature and tree-level scalar potential in Eq. (A.16). These corrected masses will be used in the calculation of cross sections for axion production above the weak scale (see App. C for more details).

C Conventions and Useful Results III: Cross Sections
We present results for the cross section for each binary collisions producing hot axions in the final state. With the only exception of thermal gluon scatterings, which we evaluated in the previous Appendix, these are the processes that we have to account for. We express each cross section as a function of the (squared of the) center of mass energy, and we evaluate the thermal average as prescribed by Eq. (4.6). The interaction rate in Eq. (4.4) is what we need to incorporate into our Boltzmann equation analysis. We employ the FeynCalc package to check all analytical expressions for cross sections [125,126].

KSVZ axion above the heavy fermion threshold
At temperatures larges than m Ψ the dominant processes for axion production are scatterings of the heavy colored PQ fermion through the interaction in Eq. (2.5). We set the axion decay constant f a to normalize the gluon anomalous coupling, and we find the cross sections

DFSZ axion above the heavy Higgs bosons threshold
In this phase, we find it convenient to work with the linear realization of PQ symmetry with the only axion interaction in Eq. (3.6). Quark-antiquark annihilations have cross sections When Higgs fields appear in the initial state we have the cross sections for scatterings mediated by up-type and down-type quark Yukawa interactions, respectively. The expressions in Eqs. (C.4) and Eq. (C.7) describe lepton scatterings withŶ (d) →Ŷ (l) . The doublets H u and H d are not mass eigenstate but we provide a simple two-step procedure to convert the cross sections above into the ones for mass eigenstates. 6 1. We introduce the temperature dependent mixing angle α where H 1 and H 2 denote the lighter and the heavier physical states, respectively. As expected, we recover α ≈ β at low temperatures, T (2B/ sin 2β) 1/2 .
2. We can replace the interaction states H u,d in Eqs. (C.3)-(C.7) with the physical eigenstates H 1,2 through the rotation above. As an example Furthermore, there are the additional contributions to the axion production from gauge boson scatterings. In this case, since the Higgs doublets are identical in terms of the gauge charge assignment there are no mixing angles appearing in the cross section for physical states. After straightforward calculations, we find the cross section for the case of the initial gauge boson state 3 (C.11) with s i ≡ m 2 H i /s and g V the corresponding gauge coupling. If the gauge boson in the final state we find (C.12) DFSZ axion below the heavy Higgs bosons and above the EWPT As the universe cools down further below the heavy Higgs boson masses, such heavy degrees of freedom are integrated out and the bath contains effectively only SM particles. We employ here the non-linear realization of the PQ symmetry with axion couplings given in Eq. (3.9). Matrix elements of SM fermion scatterings depend only on the combinations [51] for the up-type quarks u = (u, c, t), the down-type quarks d = (d, s, b), and the chargedlepton e = (e, µ, τ ), respectively. The chirality flip mentioned in the main text is such that only processes with the components of the complex Higgs doublet contribute to the rate. Thus the rate will be dominated by third generation SM fermions since their interaction strength with the Higgs field is proportional to the Yukawa couplings. We parameterize the Higgs field H T = (χ + , χ 0 ), where each doublet component is a complex scalar field, and we also introduce χ − ≡ χ † + and χ c 0 ≡ χ † 0 . The scattering cross sections take a particular simple form once we ignore CKM factors, which lead only to few percent corrections since the rate is controlled by third generation fermions. If both initial state particles are SM fermions Here, f is a generic SM fermions and y f the associated Yukawa coupling in the basis where such a coupling is diagonal. The fermion f is the weak-isospin partner of f . If a scalar appears in the initial state we have (C.17) DFSZ axion below the EWPT and above the QCDPT Below the EWPT, SM fermions and gauge bosons acquire a finite mass. We perform calculations in this phase with the PQ symmetry non-linearly realized, and cross sections still depend only on the same combinations in Eq. (C.13). We report here explicit expressions for quark scattering cross sections, the lepton case is a straightforward generalization. Here, we generalize the results provided by Ref. [51] by accounting also for flavor-violating processes whose contributions lead to corrections proportional to CKM factors. We begin with quark/antiquark annihilations. For final state gluons we have If we replace the gluon with the SM Higgs boson, we find Quarks can also annihilate to weak gauge bosons. For final state Z bosons we have for up and down quarks, respectively. Quark/antiquark annihilations to the charged weak gauge boson W ± can be flavor-changing processes, and their cross sections read We switch to quark or antiquark scattering with SM bosons. For a gluon we have For a SM Higgs boson in the initial state we find (C. 24) In the case of incident Z bosons, cross sections read for up and down quarks, respectively. Likewise, flavor-chaging processes with charged weak gauge bosons W ± give the cross section (C.27) One can easily derive the cross section of the scatterings of d j + W + → u i + a (equivalently,d j + W − →ū i + a) from Eq. (C.27) with the exchange of m u i ↔ m d j andĉ u ↔ĉ d .

KSVZ and DFSZ axions below the QCDPT
The ChPT formalism describes low-energy axion interactions with the strong sector. As discussed in the main text, we trust calculations in this regime only up to Λ ChPT ∼ 100 MeV. Thus axion production is dominated by pion scatterings mediated by the Lagrangian ∂ µ a f a c aπππ f π π 0 π + ∂ µ π − + π 0 π − ∂ µ π + − 2π + π − ∂ µ π 0 . (C.28) where g * (T ) denotes the effective number of relativistic degrees of freedom contributing to the energy density. Another crucial property of the thermal bath is its entropy density Likewise, g * s (T ) are the effective number of entropic relativistic degrees of freedom.

Temperature as the evolution variable
The cosmic time t appearing in the FLRW metric in Eq. (D.1) is not the most convenient variable to describe the evolution of a radiation dominated universe. The presence of a thermal bath makes the temperature T of the bath itself the most natural variable to keep track of the expansion. For a radiation dominated universe the expansion is adiabatic and the entropy in a comoving volume s R a 3 does not change with time We plug the definition given in Eq. (D.4) into Eq. (D.5) and we find Given a generic function of time ξ, such as the axion number density n a appearing in the Boltzmann equation, we can trade easily time with temperature derivatives It is often convenient to employ the dimensionless evolution variable x ≡ M/T , with the overall mass scale M purely conventional. Thus we find another useful relation Temperature dependence of g SM * (T ) and g SM * s (T ) At large temperatures all the degrees of freedom are relativistic so g * (T ) and g * s (T ) are constant. However, we consider axion production at temperatures below the weak scale where these quantities change as SM particles become non-relativistic. It is worth noting thanks to Eq. (D.8) that not only the absolute values matter but also their temperature derivatives. This effect is particularly significant around the QCDPT. We employ in our analysis the two different choices for the SM effective relativistic degrees of freedom.
• Ref. [32]. At temperatures above the EW scale all particles are considered free and massless and respecting the Stephan-Boltzmann law for bosons and fermions due to the crossover nature of EW transition in the SM. For massive particles around and below the EW scale when the temperature reaches each particle mass one should follow Fermi and Bose statistics. For the strongly interacting fluid especially above g✶ -Ref. [32] g✶ -Ref. [33] 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 0 20 40 60 80
100 MeV including the crossover QCD transition at 150 MeV the result of lattice simulation is used for up, down and strange quarks (2 + 1 flavors) added to the result for the charm quark at 1 GeV. Then they matched to the free gas limit at very high temperatures. For temperatures below 100 MeV the hadron resonance gas result for equation of state is used that matches to the lattice simulation of equation of state below the QCD transition epoch. At temperatures around 1 MeV the result of evolution of neutrino temperature with respect to photon temperature that includes the effect of decoupling of different types of neutrinos is implemented. In this model the number effective neutrinos based on previous calculation is assumed as N eff 3.046. The rest of SM particles considered free. Considering all these effects improves the calculation for the extra number of relativistic particles for any given models. There are uncertainties on hadron resonance gas model, lattice simulation, and thermal effect of QCD at high temperatures, and electroweak transition.
• Ref. [33]. This study uses a different treatment for the electroweak and QCD transitions, hadron resonance gas model, and neutrino decoupling. Around the EW transition the thermal corrections on the Higgs field evolution including the perturbative and nonperturbative effects for the interaction in the EW sector are used. Since around the EW transition the change of d.o.f. is not abrupt like the QCD case, due to lesser interacting particles in the thermal bath, these corrections will have tiny effects on the final result. Below 120 MeV the hadron resonance gas model and above that the QCD equation of state from a different lattice simulation for 2+1+1 flavors are used. Then it is freely matched to the perturbative QCD equation of state above 1 GeV. Around 1 MeV the neutrino decoupling is considered assuming N eff 3.045. Also, the negligible effects of plasma on electrons and photons are illustrated.
We compare the two different treatments in Fig. 14 where we show the temperature evolution for g SM * (T ) (left panel) and g SM * s (T ) (right panel). As it is manifest from these results, theoretical uncertainties will cause at most 10% difference in our prediction for the energy density stored in axion dark radiation. Furthermore, there are additional contributions to the effective relativistic degrees of freedom at high temperatures for the frameworks studied in this paper: the heavy PQ fermion and the extra Higgs bosons for the KSVZ and the DFSZ models, respectively. We include their effects by treating them as free particles with contributions given by the integrals in Eqs. (2.9) and (2.10) of Ref. [32].

How to compute ∆N eff
We provide the definition for the effective numbers of neutrino species ∆N eff valid for a generic dark radiation candidate Φ. When the universe was approximately 380,000 years old, the plasma opacity to electromagnetic radiation suddenly dropped and photons freestreamed until they reached our detectors today. At this stage the bath temperature was approximately T CMB 0.3 eV, and the only relativistic SM degrees of freedom were photons and neutrinos. The total energy density stored in radiation reads ρ R (T CMB ) = ρ γ + ρ ν + ρ Φ = 1 + We evaluate ∆N eff for a dark radiation candidate Φ that reaches thermal equilibrium with the bath at early times and it decouples subsequently. Thermal equilibrium erases the memory of whatever happened at earlier times and we can neglect physics before decoupling. As long as Φ is coupled, the number and energy densities result in n Φ (T ) = g nΦ ζ(3) π 2 T 3 , g nΦ = g Φ 1 boson 3/4 fermion , (D.11) ρ Φ (T ) = g * Φ π 2 30 T 4 , g * Φ = g Φ 1 boson 7/8 fermion . (D.12) Here, g Φ is a constant number accounting for the internal degrees of freedom (e.g., spin) of the particle Φ, and the distinction between bosons and fermions is due to the different phase-space equilibrium distributions. The Riemann ζ function appearing in the number density is approximately ζ(3) 1.2. We get rid of the temperature in the equations above to find the relation between energy and number densities π 7/2 ζ(3) with Y Φ = n Φ /s R the Φ comoving number density. After decoupling, which happens at a temperature T D , Φ's just free-stream: the phase-space distribution keeps a thermal shape with temperature red-shifting with the scale factor as T Φ ∝ a −1 , and the number density gets diluted as n Φ ∝ a −3 . Thus the comoving number density stays constant because of entropy conservation throughout the expansion Unlike Eq. (D.17), this result contains only the SM contribution to the entropic degrees of freedom that we can read off the plots in Fig. 14. Furthermore, this relation is consistent with the temperature ratio at the CMB formation as dictated by entropy conservation The numerical result in Eq. (1.3) of the introduction is a consequence of Eq. (D.20), and the value g SM * s (T CMB ) = 2 + N SM eff × (7/11) 3.94 that we use accounts for non-instantaneous neutrino decoupling. The output of this analysis describes the curves in Fig. 1.
The case discussed above is not the most general one. Thermalization may not be achieved, and even if Φ's reach thermal equilibrium the temperature T D is not the most practical variable to employ. As we do in our analysis for the QCD axion, the standard procedure is to solve the Boltzmann equation and find the asymptotic density. We conclude this Appendix with the explanation of how to use such an asymptotic value to find ∆N eff . The starting point is still Eq. (D.15) since it does not rely upon any assumption about thermalization. The only unknown quantity in that expression is the number of effective entropic degrees of freedom at recombination g * s (T CMB ): the SM part is known, we need to quantify the contribution from Φ in terms of Y Φ (T CMB ) The full number of entropic relativistic degrees of freedom appearing after the last equality is given by the two contributions in Eq. (D.18). Thus the above equation allows us to solve for g Φ * s (T CMB ) and eventually for g * s (T CMB ), and we find our final result The second term in the denominator accounts for the entropy associated to the dark radiation particle Φ. We estimate its relevance by looking back at the case when Φ's decouple at the temperature T D , and we plug the explicit equilibrium comoving density as given in Eq. (D. 16). We find that the correction results in g * Φ /(g SM * s (T D ) + g * Φ ), and therefore it is relevant only if the dark radiation stays in thermal equilibrium until a time when its effective number of entropic degrees of freedom is comparable with the one of the SM bath.
Theoretical uncertainty on ∆N eff due to the interpolation We conclude this Appendix with a discussion of theoretical uncertainties associated with our interpolation across the QCDPT. In our work, we performed a smooth interpolation for the production rate between Λ ChPT and Λ N with the cubic 'spline' method, motivated by the fact that the QCDPT is a crossover leading to mild shifts of thermal properties. One may wonder how our predictions for ∆N eff are sensitive to the details of such an interpolating method. We take the KSVZ axion case for concreteness, and we modify the production rate as shown in the left panel of Fig. 15. The solid black line corresponds to the rate used in our analysis. We consider two extreme cases where around the temperature scale 500 MeV the actual rate is a factor of two larger (red line) or smaller (green line), and we make sure to match these lines with our results below Λ ChPT and above the mass of QCD resonances. The resulting predictions for ∆N eff are shown in the right panel of Fig. 15. For values of the axion decay constant not excluded experimentally, our predictions are quite insensitive to the detail of the interpolation and therefore utterly solid.