Observable windows for the QCD axion through $N_\text{eff}$

We show that when the QCD axion is directly coupled to quarks with $c_q/f \, \partial_\mu a \, \bar{q} \gamma^\mu \gamma^5 q$, such as in DFSZ models, the dominant production mechanism in the early universe at temperatures $1 \, {\rm GeV}\lesssim T \lesssim 100 \,{\rm GeV}$ is obtained via $q \bar{q} \leftrightarrow g a$ and $q g \leftrightarrow q a$, where $g$ are gluons. Different heavy quarks $q_i$ can produce a thermal axion background that decouples at a temperature $T_i$: (1) top quark at $T_t \lesssim 100 \, {\rm GeV}$ for $f/c_t \lesssim 3\times 10^8 {\rm GeV}$; (2) bottom quark at $T_b \lesssim m_b$, for $f/c_b\lesssim 8\times 10^{7} {\rm GeV}$; (3) charm quark at $T_c \lesssim m_c $ for $f/c_c \lesssim 5\times 10^{7} {\rm GeV}$. Each of these cases corresponds to a contribution to the effective number of relativistic degrees of freedom, in the windows given by $0.027 \leq \Delta N_\text{eff}\leq 0.031$, $0.037 \leq \Delta N_\text{eff} \leq 0.039$ and $0.039 \leq \Delta N_\text{eff}$, respectively. These contributions are larger than the one obtained when thermalization happens only above the electroweak phase transition, $\Delta N_{\rm eff}\lesssim 0.027$, and are within reach of future CMB S4 experiments, thus opening an alternative window to detect the axion and to test the early universe at such temperatures.


I. INTRODUCTION
The axion (a) is a field postulated for theoretical reasons [1,2], that possesses a very rich phenomenology [3].It has a coupling to gluons of the form α s /(8πf ) aG G, that can provide a compelling solution to the strong CP problem [4], i.e. the fact that the total coefficient of the CP violating operator θ/(8π)G G has a coefficient consistent with zero in vacuum with high precision, |θ| 10 −10 [5,6].Experimentally the possible range for f is bounded by f 10 6 − 10 9 GeV, via direct detection experiments [3] or stellar cooling bounds [7,8] 1 , although different bounds refer to coupling to different particles, namely photons, nucleons, quarks or leptons.From Supernova SN1987A [3,[9][10][11] one gets f O(10 8 ) GeV, but such bound should be considered less robust due to complicated SN physics.Moreover, the SN bound relies specifically on the axion-nucleons couplings and so is somehow model dependent.In fact the low energy effective Lagrangian for the QCD axion can arise from several UV completions, with different properties, leading thus to different effective couplings to the various Standard Model (SM) fields, including quarks and nucleons.The coherent oscillations of the axion can lead also to a good candidate for Dark Matter, for f 10 11 GeV [3,[12][13][14].
Axions also naturally constitute a population of detectable hot relics from the early universe [15][16][17][18][19][20][21] if the interactions that can produce them are fast enough at some stage of the radiation dominated universe.A rather generic expectation is that axions can be produced in thermal abundance either via interactions with gluons [18,19], or via direct interactions with the top quark [15,20], and they subsequently decouple and freeze out at a high temperature T D , before the Electroweak phase transition (EWPT).
Their relic abundance in this case would be proportional to the inverse of the effective number of degrees of freedom at the decoupling temperature g * ,D and this would affect the energy density of relativistic particles at recombination, traditionally parameterized by an effective number of neutrinos N eff , which can be tested by the CMB.Assuming only the existence of the SM plus the axion above the EWPT, the predicted change in N eff due to axions would be given by ∆N eff = 4/7 (43/(4g * ,D )) 4/3 ≈ 0.027 [22].Remarkably this is comparable to the forecasted sensitivity of future CMB-S4 experiments [23,24].If there are additional degrees of freedom beyond the SM at such high temperature T D , this number would be smaller, so that this value is actually an upper bound.A possible thermalization at much later stages, below the QCD phase transition (QCDPT), via scattering with pions and nucleons, has been studied in [16,17], concluding that this would require f O(10 6 ) GeV, which is already in conflict with the astrophysical bounds [21].
In this Letter we point out that if the axion is produced instead between the EWPT and the QCDPT, it can erase the previous prediction due to the high temperature thermalization before the EWPT, leading instead to larger values of ∆N eff , thus slightly easier to detect.Moreover, although there is room to have extra relativistic species at temperatures between the QCDPT and the EWPT, e.g.light relics or other weakly coupled particles, so far the Standard Model has passed precision tests with great success and so in this window of temperatures the values of ∆N eff can be seen more confidently as true predictions, as a function of f , instead of upper bounds.This happens due to direct couplings to various heavy quarks which are present in many axion models, as estimated already with dimensional analysis in [15].Similar couplings were also studied in [25] and [23], although this was applied mostly to leptons, whose cross-sections are suppressed due to the absence of coupling to gluons.We explicitly compute here the relevant cross-sections, with t-b and c quarks, at temperatures sufficiently above the QCDPT, so that we can use a perturbative approach, which should be roughly correct as long as α s 1.This calculation is timely because of the planned CMB-S4 experiments, which should be able to go close to the necessary sensitivity.

II. AXIAL COUPLINGS
The axion has a Lagrangian that includes couplings with gauge bosons and can also include direct derivative couplings with fermions.We are in particular interested in the direct interaction with the different quarks q i of the SM, The coefficients c i are non-zero already at tree level, for example, in DFSZ models [26,27] and are typically of order unity2 .Note also that in UV-completed axion models the c i are proportional to the Peccei-Quinn charges over the domain wall number N DW [22].At high temperature, above the EWPT, this type of couplings are the most efficient ones that can thermalize the axion [20], through the scatterings t + h → t + a, t + h → t + a, t + t → h + a, where h is the Higgs particle.This is even larger than the scattering with gluons because it is not suppressed by α s /(2π), contrary to the axiongluon coupling, jointly with the fact that the top Yukawa coupling is y t 1.By dimensional analysis the rate of the axion-top-Higgs process goes as Γ t−h ∝ y 2 t c 2 t T 3 /f 2 , and so can be larger than the Hubble rate H = g * /90πT 2 /M p , at sufficiently high temperature, leading to an equilibrium distribution and subsequent decoupling when Γ t−h H. Here, M p = 2.4 × 10 18 GeV.Such a rate is larger than the one through gluons, which is given parametrically by Γ g ∝ (α s /(2π)) 2 g 2 s T 3 /f 2 , where g s ≡ √ 4πα s .In this Letter, we will instead look at temperatures below the EWPT.At these energies the Higgs particles are Boltzmann suppressed, so a much more relevant scattering is obtained replacing the Higgs with a gluon, jointly with the fact that g s 1.So we will consider the processes A crucial feature of these cross sections is that, similarly to the case of Γ t−h , they are proportional to the Yukawa coupling squared, y 2 i , which means to the squared mass of the quark, since the Higgs takes a vev v below the EWPT.This can be understood performing a field redefinition, as in [20], which erases the coupling of eq. ( 1) and replaces it with a coupling of the type y i c i /f ahq i γ 5 q i .When the Higgs has a vev, one can set h = v and so the quark mass (y i v) 2 /2 = m 2 i will actually appear in the cross sections.We work in the present letter without performing such a redefinition, but we verify by explicit calculation that all cross sections indeed go as m 2 i .For T m i this represents a suppression.However, when T m i that is no longer the case, and such interactions dominate, e.g., over the ones with gluons.
For temperatures above the QCDPT, the cross sections σ i can be evaluated straightforwardly and they are rescaled for the various quarks q i simply by the respective m 2 i .In the end we are interested in computing standard thermally averaged rates (see e.g.[18,22]), where we integrate over the 3-momenta of the incoming states, f a,b are their thermal distributions (Bose-Einstein or Fermi-Dirac) and n eq is the axion equilibrium number density.Parametrically [15] one has If, instead, T m i a Boltzmann suppression arises, due to f a and f b , leading to In order to find the precise numerical prefactor we explicitly compute the cross-sections.We simply use tree-level diagrams, adding also thermal masses in the kinematics; for a more refined treatment one could evaluate in thermal field theory the axion self-energy, as in [19,20], although we expect this will lead only to minor corrections.Comparing the scattering rates Γ i with H allows us to know, at a given temperature, at which value of f the axion is in thermal equilibrium.More precisely, we check whether Γ i /(3H) > 1.At a given T the most important channel is the scattering with the heaviest quark which is not yet Boltzmann suppressed.Thus, the rates Γ i are maximal around T m i .
For the first two processes in eq. ( 2) the cross section can be quite simply expressed under the approximation of massless gluons (i.e.neglecting thermal masses) as where p and E are the momentum and energy of the quark, respectively, in the center of mass frame.Instead, the cross section for pair annihilation takes a simple form even including the gluon thermal mass m g,th , 4σ In our final results we will include gluon thermal masses [28], where N c and N f are respectively the number of colors and flavors.We will also use in the quark dispersion relations an overall effective mass which includes both the standard mass m i ≡ y i v/ √ 2 and the thermal mass [28] Note however that, while the thermal masses modify the kinematics, they cannot enter in the overall prefactor y 2 i v 2 that appears in the crosssections, since the thermal masses change the dispersion relations, but preserve chirality [28] 3 .This can also be seen by performing again the field redefinition of [20], which shows that the cross section must vanish for y i → 0, and therefore there cannot be a contribution to the prefactor proportional to m i,th , which is independent of y i .The inclusion of gluon thermal masses does not change dramatically the result, but it is important conceptually, because it kinematically forbids processes with emission of many gluons, which could be otherwise dominant since we have g s > 1.Finally, we evaluate the couplings inside the integral at the center-ofmass energy g s (2E) by using a simple one-loop RGE running.
We find that the rates are well approximated by the following fitting functions , where M Z is the Z-boson mass and

III. SCATTERING RATES AND N eff : RESULTS
In fig. 1 we plot Γ i /(3H) as a function of T , for the scatterings involving t, b and c quarks4 , plotting each curve at two different values of f /c i .Note that what determines the final axion abundance is the decoupling temperature.In all our plots we use the one-loop analytical expression for the temperature dependence of the Higgs vev v(T ), which is relevant close to 100 GeV, especially in the case of the t-quark.At large values of f , say f /c t 10 9 GeV, the axion can become thermal via the top-Higgs scatterings5 already at T O(100) GeV.However, the axion Γi/(3H) vs temperature associated with the axion production via scatterings with different quarks and gluons.Each process is evaluated at values of f /ci chosen to be the minimal ones such that Γi/(3H) = 1.The orange line uses the result derived in [20].
can also be produced below the EWPT through the scatterings with the t-quark, eq. ( 2), for f /c t 10 9 GeV or, for even smaller f , through scatterings with the b or c quarks.
The important consequence of axion production is the observable effect on the effective number of relativistic degrees of freedom N eff .In the SM this is predicted to be [30,31] N eff = 3.045, which tells us that the only relativistic species at recombination, apart from photons, are the 3 neutrinos.Any deviation from this would tell us something new about the universe, in particular, about light relics and/or additional degrees of freedom.
A simple way to estimate ∆N eff would be simply to read the decoupling temperature from fig. 1, and to evaluate g * D , i.e. from [32].In order to compute ∆N eff more precisely we solve the Boltzmann equation 6 for the axion number density n a , This equation can also be rewritten in terms of the abundance Y ≡ n a /s, where z ≡ m i /T , Y eq ≡ n eq /s, γ i ≡ Γ i n eq and s ≡ 2π 2 g * T 3 /45 is the entropy density.Assum- This equation is valid in the Boltzmann approximation. .∆N eff as a function of f /ci.The green, red and blue label correspond to the predictions for the charm, bottom and top particle, respectively.The orange bands represent the 1σ and 2σ forecasted contours of CMBS4, plus a more futuristic 1σ band, according to [24].See also [33] for similar forecasts combining CMB-S4 and large scale structure. .
ing a constant g * there is an analytical solution.
Parameterizing γ i ≡ γ mi z −4 e −z , H ≡ H mi z −2 and s ≡ s mi z −3 and imposing Y = 0 at z → 0 we get: where ) which corresponds to a change in N eff given by However, g * has a non-trivial dependence with T [32], so we need to solve the Boltzmann equation numerically 7 .We stop the simulation at T = 0.5 GeV, when α s = 1.In the case of t the results are insensitive to this cutoff.In the case of b the results are also insensitive to it, provided f /c b 10 7 GeV.However, if the axion production is still non-zero at T 0.5 GeV, as happens for the b, when f /c b 10 7 GeV, and for the c, when f /c c 2 × 10 8 GeV, that can affect ∆N eff , given that g * (T ) varies rapidly in 7 Other temperature dependences that we also included numerically are: the running of gs, the Higgs vev v(T ) and the logarithmic part of the rates Γ i .
that region.To get a good estimate one would need more refined techniques, which we leave for future work.We plot the results in fig.2, showing the predicted ∆N eff as a function of f /c i for i = t, b, c.As anticipated from fig. 1 thermalization via top scatterings happens when f /c t ≈ 4 × 10 8 GeV, leading to ∆N eff ≈ 0.03.For smaller f /c t the axion remains thermal for longer time, leading to a larger ∆N eff .For example, at 10 6 GeV ≤ f /c t ≤ 10 7 GeV, ∆N eff = 0.035 − 0.036.Note, however, that such values of the axial couplings to the top might already be excluded by stellar cooling constraints, f /c t 2.4 × 10 9 GeV [8].
Similarly, the scatterings with the bottom are able to thermalize the axion at f /c b ≈ (2 − 3) × 10 8 GeV, leading to ∆N eff = 0.038.In this case, at f /c b = 10 7 GeV, one has ∆N eff = 0.047 which is interestingly close to the forecasted 2σ contour of CMB-S4 [24].For f /c b below 10 7 GeV we encounter the problem that we have stressed before: the axion decouples at T 0.5 GeV, where the methods used here break down.To see this sensitivity we have plotted 4 red lines.The lowest one corresponds to a case where we stop the simulation at a final temperature T F = 0.5 GeV.Such line should be seen as a lower bound on ∆N eff at any given value of f /c b .The other 3 lines correspond to cases where we fix by hand α s = 1 for T < 0.5 GeV and then run the simulation down to lower temperatures T F = 0.4, 0.3, 0.2 GeV, respectively.The lowest red line gives a lower bound on ∆N eff , which at at 10 6 GeV ≤ f /c b ≤ 10 7 GeV, is ∆N eff > 0.049.
Regarding the scatterings with the charm they are able to thermalize around f /c c 10 8 GeV.This case is affected, to a greater extent, by the same uncertainty due to strong coupling and the final predictions turn out to be very sensitive to what happens below 0.5GeV.Also in this case we plot 4 lines to show the sensitivity to T F and again, from the lowest green line we can give a lower bound on ∆N eff which at f /c c < 10 8 GeV is ∆N eff > 0.045.Such values, obtained close to the strong coupling region, are of great interest, since they could be potentially measurable at 2σ by CMB-S4.However, a dedicated study of such cases would be needed, using more adequate methods to treat the axion decoupling around the QCDPT.
Finally, note that in all cases even if the axion does not thermalize completely it can still have a non-negligible abundance, and even give ∆N eff > 0.027.
In conclusion, we have shown that, in addi-tion to the possible thermalization at f 10 9 GeV, before the EWPT via gluons [18,19] or axion-top-Higgs scatterings [20], the axion can be produced at intermediate energy scales between the EWPT and the QCDPT via direct couplings to heavy quarks.Such couplings are present in various models, such as DFSZ, at tree level and thus are enhanced compared to interactions with gluons, which are generated at one loop.Clearly, because we are probing lower energies, this happens in a lower range of f , which by solving the associated Boltzmann equation turns out to be f /c i 10 9 GeV.This window can be viable if SN bounds are not so robust [11] or, as pointed out in [34], in DFSZ-like axion models which can avoid SN bounds, because of suppressed coupling to light quarks and thus to nucleons.Interestingly in the same window one could also explain the hints of extra stellar energy losses, as suggested in [11] 8 .
Axion production and decoupling happen at different thresholds T i , close to the heavy quark masses m i , leading to different predictions on ∆N eff as a function of f /c i , shown in fig.2, which contains our main results.We stress that such predictions are not upper bounds, but can be precisely given as far as we know the number of degrees of freedom g * ,D and the dynamics at such energies.For the top and bottom our main predictions (0.027 ≤ ∆N eff ≤ 0.036 for 10 6 GeV f /c t 4 × 10 8 GeV and 0.027 ≤ ∆N eff ≤ 0.047 for 10 7 GeV f /c b 3×10 8 GeV) are quite under control, while for the c-quark, since axion production is too close to the QCD phase transition, we only have a lower bound ∆N eff ≥ 0.045 (at f /c c 10 8 GeV), which nonetheless might represent the largest signal for a detection.
These observations open a new window to test the QCD axion.At the same time, if the axion exists in this range of f , this would allow us to probe the Early universe at energies so far unexplored in cosmology, between the QCD and the EW phase transitions.If the axion could 8 Note also that, although this region does not overlap with the standard window of axion dark matter, it may be possible to obtain QCD axion DM even at such small f , via domain wall decay [11,35,36] or via resonant production of axions at the PQ phase transition [37], although in such cases one should check whether the scatterings we consider can destroy the condensate.Furthermore our window of f can also be relevant in the context of Axion-like particles (ALPs), as in [38].
be detected directly within such window for f , e.g. with the CAST [39] or IAXO [40][41][42] experiments, our conclusions would make possible to confirm if our expectations about the total particle content and the thermal history at such temperatures are indeed correct.

t+g→t+a (f/ct=3x10 8
Figure1.Γi/(3H) vs temperature associated with the axion production via scatterings with different quarks and gluons.Each process is evaluated at values of f /ci chosen to be the minimal ones such that Γi/(3H) = 1.The orange line uses the result derived in[20].

Figure 2
Figure2.∆N eff as a function of f /ci.The green, red and blue label correspond to the predictions for the charm, bottom and top particle, respectively.The orange bands represent the 1σ and 2σ forecasted contours of CMBS4, plus a more futuristic 1σ band, according to[24].See also[33] for similar forecasts combining CMB-S4 and large scale structure..