New Bounds on Heavy QCD Axions from Big Bang Nucleosynthesis

TL;DR

This study develops a comprehensive BBN-based framework to constrain heavy QCD axions that decay hadronically. By solving a coupled set of Boltzmann equations for the neutron fraction, axions, and hadrons, and by incorporating updated hadronic cross sections, energetic $K_L$ treatment, and secondary hadrons, the authors derive a robust bound $\tau_a \lesssim 0.02\,\mathrm{s}$ for $m_a \gtrsim 0.3\,\mathrm{GeV}$, stronger than projected $N_{\rm eff}$ limits from the CMB. The analysis spans two decay regimes, employing a data-driven method for $m_a<2$ GeV and pQCD with hadronization models for $m_a>2$ GeV, and demonstrates the bound’s insensitivity to branching fractions, cross sections, and initial axion abundance. The results have broad implications for hadronically decaying long-lived particles and provide a robust methodology for including kaon dynamics and secondary hadrons in early-Universe studies. Overall, the work showcases BBN as a powerful probe of beyond-Standard-Model physics with long-lived, hadron-injecting particles and yields a competitive and often superior constraint compared to future CMB measurements.

Abstract

We study Big Bang Nucleosynthesis (BBN) constraints on heavy QCD axions. BBN offers a powerful probe of new physics that modifies the neutron-to-proton ratio during the process, thanks to the precisely measured primordial Helium-4 abundance. A heavy QCD axion provides an attractive target for this probe, because not only is it a well-motivated hypothetical particle by the strong CP problem, but also it dominantly decays to hadrons if kinematically allowed. A range of its lifetime is thus excluded where the hadronic decays would significantly alter the neutron-to-proton ratio. We compute axion-induced modification of the neutron-to-proton ratio, and obtain robust upper bounds on the axion lifetimes, as low as 0.017 s for the axion mass higher than 300 MeV. Remarkably, this is stronger than projected future CMB bounds via $N_{\rm eff}$. Our bounds are largely insensitive to uncertainties in hadronic cross sections and the axion's branching fractions into various hadrons, as well as to the precise value of the initial axion abundance. We also incorporate, for the first time, several key improvements, such as scattering processes by energetic $K_L$ and secondary hadrons, that can also be important for studying general hadronic injections during BBN, not limited to those from axion decays.

Figures (27)

• Figure 1: Schematic depiction of the thermal history of the universe as well as steps taken in our analysis as summarized in Sec.$\,$\ref{['sec:overall']}. We adopt the 2$\sigma$ bound of $^4{\rm He}$ abundance $R_{{\rm Y}_{\rm p}} = 2.45\%$.
• Figure 2: ${\cal F}_\pi(m_a, T)$ used in Eq. \ref{['Eq:a_pi']} for different values of $m_a$. The dots indicate our numerical estimation, and the lines are interpolations of the dots.
• Figure 3: $Y_a^{\rm (min)}$ (solid lines) and $Y_a^{\rm (max)}$ (dashed lines) as functions of $f_a$ for various choices of $m_a$ represented by different colors. The sharp drops at $f_a \sim 10^8$ and $10^9\,\mathrm{GeV}$ for the $m_a=1$ and $10\,\mathrm{GeV}$ cases, respectively, are our truncations by hand because their respective temperature $T_{\rm decay}$ becomes above $T_{\pi}$ in the region indicated by "$\leftarrow$".
• Figure 4: The left (right) plots show $Y_a^{\rm (min)}$ ($Y_a^{\rm (max)}$) in the $(m_a,\,f_a)$ space (top) and in the $(m_a,\,\tau_a)$ space (bottom). The blue, red, and cyan dashed lines in the upper plots correspond to $\tau_a=1\,\mathrm{s}$, $0.02\,\mathrm{s}$, and $0.01\,\mathrm{s}$, respectively.
• Figure 5: The total decay width $\Gamma_a$ as a function of $m_a$. Here, $f_a$ is fixed to $1\,\mathrm{TeV}$, but one can re-scale the rate by $(\mathrm{TeV}/f_a)^2$ for a different value $f_a$. Estimations based on pQCD and the data-driven method are shown in blue and red, respectively. The bands represent the uncertainties of the estimations, where the blue band is determined by varying the renormalization scale from $\mu=m_a/2$ to $\mu=2m_a$, while we assign a factor-of-$2$ uncertainty to the data-driven estimation for the red band.