# Axion Dark Matter from Heavy Quarks

Mohammad Aghaie  Dipartimento di Fisica E. Fermi, Università di Pisa and INFN-
Pisa, Largo B. Pontecorvo 3, I-56127 Pisa, Italy  Giovanni Armando
Dipartimento di Fisica E. Fermi, Università di Pisa and INFN-Pisa, Largo B.
Pontecorvo 3, I-56127 Pisa, Italy  Angela Conaci  Dipartimento di Fisica E.
Fermi, Università di Pisa and INFN-Pisa, Largo B. Pontecorvo 3, I-56127 Pisa,
Italy  Dipartimento di Fisica, Università della Calabria and INFN-Cosenza,
Arcavacata di Rende, I-87036 Cosenza, Italy  Alessandro Dondarini
Dipartimento di Fisica E. Fermi, Università di Pisa and INFN-Pisa, Largo B.
Pontecorvo 3, I-56127 Pisa, Italy  Institute for Theoretical Particle Physics,
KIT, 76128 Karlsruhe, Germany  Peter Maták  Department of Theoretical Physics,
Comenius University in Bratislava, Mlynská dolina, 84248 Bratislava, Slovak
Republic  Paolo Panci  Dipartimento di Fisica E. Fermi, Università di Pisa and
INFN-Pisa, Largo B. Pontecorvo 3, I-56127 Pisa, Italy  Zuzana Šinská
Department of Theoretical Physics, Comenius University in Bratislava, Mlynská
dolina, 84248 Bratislava, Slovak Republic  Robert Ziegler  Institute for
Theoretical Particle Physics, KIT, 76128 Karlsruhe, Germany  Institute for
Theoretical Physics, Heidelberg University, 69120 Heidelberg, Germany

###### Abstract

We propose simple scenarios where the observed dark matter abundance arises
from decays and scatterings of heavy quarks through freeze-in of an axion-like
particle with mass in the $10{\rm\,keV}-1{\rm\,MeV}$ range. These models can
be tested by future X-ray telescopes, and in some cases will be almost
entirely probed by searches for two-body decays $K\to\pi+{\rm invis.}$ at
NA62. As a byproduct, we discuss the cancellation of IR divergencies in
flavor-violating scattering processes relevant for thermal axion production,
and derive the general contribution to axion-photon couplings from all three
light quarks.

††preprint: TTP24-009, P3H-24-026

##  I Introduction

QCD Axions and axion-like particles (ALPs) with masses below the MeV scale are
excellent Dark Matter (DM) candidates, provided that the associated Peccei-
Quinn (PQ) breaking scale $f_{a}$ is sufficiently large in order to ensure
stability on cosmological scales. These particles are light enough to be
produced in stellar plasmas, and constraints from star cooling typically
require $f_{a}\gtrsim 10^{9}{\rm\,GeV}$, seemingly rendering axion production
at particle colliders hopeless. However, collider searches are actually
sensitive to such large scales, if the axion has flavor-violating (FV)
couplings to SM fermions Feng:1997tn ; Kamenik:2011vy ; Bjorkeroth:2018dzu ;
MartinCamalich:2020dfe ; Calibbi:2020jvd ; Ziegler:2023aoe . Precision flavor
experiments then allow to probe scales of the order of $f_{a}\sim 10^{12}$ GeV
by searching for $K\to\pi+{\rm invis.}$ at NA62 Goudzovski:2022vbt ,
$10^{10}{\rm\,GeV}$ with $\mu\to e+{\rm invis.}$ at MEG-II Calibbi:2020jvd ;
Jho:2022snj , Mu3e Knapen:2023zgi , Mu2e or COMET Hill:2023dym , and
$10^{9}{\rm\,GeV}$ with $B\to K+{\rm invis.}$ at Belle II
MartinCamalich:2020dfe .

As flavor-violating axion couplings are determined by the misalignment of PQ
charges and SM Yukawas, their prediction from UV scenarios requires a theory
of flavor. Particularly economic models of this kind can be constructed when
PQ acts as a flavor symmetry explaining Yukawa hierarchies Davidson:1981zd ;
Wilczek ; Berezhiani:1989fp ; Flaxion ; Calibbi:2016hwq , although the
resulting size of flavor-violating couplings largely depends on the particular
scenario. Here instead we link the size of flavor-violating axion couplings to
the observed DM relic abundance, requiring thermal production of DM axions in
the right amount via freeze-in of decays (and scatterings) of SM fermions. As
$\Omega_{a}\propto m_{a}\Gamma(f_{i}\to f_{j}a)$, this fixes the rate of these
decays for a given axion mass, which indeed is in the reach of current
experiments. This idea has originally been proposed in Ref. Panci:2022wlc  in
the context of lepton flavor-violating (LFV) decays, it is the purpose of this
article to extend the analysis to quarks.

The quark scenario differs from the lepton scenarios in several aspects.
Thermal axion production has to respect the Warm DM bound, $m_{a}\gtrsim
10{\rm\,keV}$, and axion decays into photons have to be sufficiently
suppressed in order to satisfy stringent constraints from X-ray and low-energy
$\gamma$-ray line searches. This requires the absence of EM and color
anomalies, so that the decay rate is additionally suppressed by power of
$m_{a}^{4}/m_{f}^{4}$. In the case of LFV decays $f=e,\mu$, so that some
hierarchy between diagonal and off-diagonal couplings is needed in order to
ensure sufficient stability. In contrast in the quark case the mass
suppression is at least $m_{a}^{4}/m_{\pi}^{4}$, thus improving axion
stability and reducing the need of coupling hierarchies. Another important
difference is the relative size of axion production rates from decays and
scattering processes. While in the LFV scenarios production from diagonal
scattering is suppressed with respect to decays by a factor of $\alpha_{\rm
em}$, in the quark case this becomes a factor $\alpha_{s}$, so that decays and
scattering are almost equally relevant for couplings of similar size. This
also implies that NLO corrections are sizable, and we will discuss the
corrections from flavor-violating scattering processes, which naively involves
IR divergences. However, we will demonstrate that such terms are cancelled in
the relevant temperature regime by thermal and virtual corrections, partially
reproducing results in Ref. Czarnecki:2011mr .

We now proceed by defining the basic framework in Section II and introduce two
classes of simple benchmark scenarios. We then discuss axion stability, axion
production in the early universe and constraints from structure formation and
astrophysics. We use these results to project present constraints and future
sensitivities on the 2-dimensional parameter of our benchmark models in
Section III, before concluding in Section IV.

##  II Framework

We consider an anomaly-free ALP $a$ that only couples to SM quarks

$\displaystyle\mathcal{L}=\frac{1}{2}(\partial_{\mu}a)^{2}-\frac{m_{a}^{2}}{2}a^{2}+\frac{\partial_{\mu}a}{2f_{a}}\overline{q}_{i}\gamma^{\mu}\left(C^{V}_{q_{i},q_{j}}+C^{A}_{q_{i},q_{j}}\gamma_{5}\right)q_{j},$
(1)

where $C^{V,A}_{q_{i},q_{j}}$ are traceless hermitian matrices in flavor
space. These couplings arise from the misalignment of PQ charges and quark
Yukawa matrices

$\displaystyle\begin{split}&C^{V,A}_{q}=U^{\dagger}_{q_{R}}X_{q_{R}}U_{q_{R}}\pm
U^{\dagger}_{q_{L}}X_{Q_{L}}U_{q_{L}}\,,\end{split}$ (2)

where $q=u,d$, and $X_{Q_{L}},X_{u_{R}},X_{d_{R}}$ are traceless diagonal
matrices containing the PQ charges of $Q_{L},u_{R},d_{R}$, respectively, while
$U_{q_{L},q_{R}}$ are unitary matrices that diagonalize the quark Yukawas
according to $Y_{q}^{\text{diag}}=U_{q_{L}}^{\dagger}Y_{q}U_{q_{R}}$, where
$V_{\rm CKM}=U_{u_{L}}^{\dagger}U_{d_{L}}$ is the CKM matrix.

Different scenarios can occur depending on the specific choices for flavor
rotations $U_{q}$ and PQ charges $X_{q}$. In the following we consider two
classes of benchmark scenarios. In the first class we take into account a
single flavor transition at a time, so only two charges of right-handed (RH)
quarks are different from zero, e.g. $X_{d_{R}}={\rm
diag}(0,1,-1),X_{u_{R}}=X_{Q_{L}}=0$. The corresponding unitary matrix is
restricted to a rotation in the same sector, i.e. is a rotation in the 2-3
plane by some angle $\alpha$ with $0\leq\alpha\leq\pi/2$. This gives

$\displaystyle C_{d}^{V}$ $\displaystyle=C_{d}^{A}=\begin{pmatrix}0&0&0\\\
0&\sin{\alpha}&\cos{\alpha}\\\ 0&\cos{\alpha}&-\sin{\alpha}\end{pmatrix}\,,$
$\displaystyle C^{V}_{u}$ $\displaystyle=C^{A}_{u}=0\,.$ (3)

We call this scenario the “$bs$ scenario”, analogously we define the $bd$,
$cu$, $sd$, $tu$ and $tc$ scenarios. These benchmarks scenarios have only
three free parameters: the ALP mass $m_{a}$, the decay constant $f_{a}$, and
the rotation angle $\alpha$ that controls the ratio of flavor-diagonal and
off-diagonal ALP couplings. We will fix one of these parameters ($f_{a}$) by
demanding that ALPs are thermally produced in the right abundance via thermal
freeze-in. As a consequence, we will obtain a two-dimensional parameter space
in the plane ($m_{a},\alpha$), which is subject to various constraints from
direct searches, astrophysics and cosmology. As we are going to discuss in
Section III, only few of the six possible scenarios are viable and give rise
to a distinct phenomenology.

The second class of scenarios is obtained by assuming that the unitary flavor
rotations are given by the CKM matrix, while PQ charges in the quark sector
are either vanishing or taken to be the most general assignment, $X_{q}={\rm
diag}(1,X,-1-X)$, where the PQ charge $X$ is a real number. Below we consider
two explicit benchmark scenarios: either only left-handed quarks are charged
under PQ, i.e., $X_{u_{R}}=X_{d_{R}}=0,X_{Q_{L}}={\rm diag}(1,X,-1-X)$, and
the CKM is coming entirely from the down-quark sector,
$U_{u_{L}}=1,U_{d_{L}}=V_{\rm CKM}$, or only right-handed down quarks are
charged, $X_{Q_{L}}=X_{u_{R}}=0,X_{d_{R}}={\rm diag}(1,X,-1-X)$ and the
relevant rotation is CKM-like, $U_{d_{R}}=V_{\rm CKM}$. We call these scenario
the “CKM${}_{Q_{L}}$ scenario” and the “CKM${}_{d_{R}}$ scenario”,
respectively. These scenarios are considered to be representative for the
phenomenology of more realistic models, where flavor-rotations are determined
by the same dynamics that explain fermion mass hierarchies, which may be the
PQ symmetry itself Flaxion ; Calibbi:2016hwq ; Linster:2018avp . As in the
first class, these two scenarios have just three parameters, where again
$f_{a}$ is determined by the observed relic abundance, leaving a two-
dimensional parameter space in the plane ($m_{a},X$), which is subject to
phenomenological constraints.

###  II.1 Dark Matter Stability

To be stable on cosmological scales, axion decays into SM particles must be
sufficiently suppressed. We will take $m_{a}\ll m_{\pi}$, so that only decays
into photons are possible. The dominant constraints on the decay rate come
from X-ray telescopes, and are of the order
$\tau_{\gamma\gamma}\gtrsim(10^{26}\div 10^{28})\,\text{sec}$, depending on
$m_{a}$, which is roughly 10 orders of magnitude larger than the age of the
Universe.

In our benchmark models, the decay $a\to\gamma\gamma$ takes place through
quarks loops. For heavy quarks one can use perturbative results, while for
lights quarks ($u,d,s$) one has to rely on chiral perturbation theory since
$m_{a}\ll\Lambda_{\rm QCD}$. In the following we use and extend the results of
Ref. Bauer:2017ris . The decay rate into photons is given by

$\displaystyle\Gamma_{\gamma\gamma}=\frac{\alpha_{\rm
em}^{2}m_{a}^{3}}{64\pi^{3}f_{a}^{2}}\left|C^{\rm heavy}_{\gamma\gamma}+C^{\rm
light}_{\gamma\gamma}\right|^{2}\,,$ (4)

where the effective photon couplings receive contributions from heavy and
light quarks. The heavy quark contribution is given by

$\displaystyle C^{\rm
heavy}_{\gamma\gamma}\approx\sum_{i=c,b,t}Q_{i}^{2}C_{i}\frac{m_{a}^{2}}{4m_{i}^{2}}\,,$
(5)

where $C_{i}\equiv C_{q_{i}q_{i}}^{A}$ and we have neglected terms of order
$m_{a}^{4}/m_{i}^{4}$. Using the results detailed in Appendix A, the light
quarks contribute dominantly through axion-$\pi$, axion-$\eta$ and
axion-$\eta^{\prime}$ mixing

$\displaystyle C^{\rm light}_{\gamma\gamma}$
$\displaystyle\approx\frac{C_{u}-C_{d}}{2}\frac{m_{a}^{2}}{m_{\pi}^{2}}+\frac{\sqrt{2}}{6}(C_{u}+C_{d}-C_{s})\frac{m_{a}^{2}}{m_{\eta}^{2}}$
$\displaystyle+\frac{\sqrt{2}}{3}(C_{u}+C_{d}+2C_{s})\frac{m_{a}^{2}}{m_{\eta^{\prime}}^{2}}\,,$
(6)

with $(m_{\pi},m_{\eta},m_{\eta^{\prime}})=(135,548,958){\rm\,MeV}$, and we
have neglected multiplicative corrections of order
$m_{a}^{2}/m_{\pi,\eta,\eta^{\prime}}^{2}$, besides small corrections from
isospin breaking.

Thus the effective coupling to photons is suppressed by at least
$m_{a}^{2}/m_{\pi}^{2}$, since there is no color nor electromagnetic anomaly
Nakayama:2014cza ; Takahashi:2020bpq ; Han:2020dwo ; Han:2022iig ;
Sakurai:2022roq . As a result the axion lifetime is given by (assuming that
axion-pion mixing is the dominant contribution)

$\displaystyle\tau_{a}\approx 3\times 10^{26}{\rm
sec}\left(\frac{0.1{\rm\,MeV}}{m_{a}}\right)^{7}\left(\frac{f_{a}/(C_{u}-C_{d})}{10^{9}{\rm\,GeV}}\right)^{2}\,,$
(7)

so that for parameters consistent with freeze-in production and WDM bounds
(see below) the axion lifetime easily exceeds the age of the universe, and can
be sufficiently large in order to satisfy the stringent limits from X-ray
telescopes.

For axion masses in the keV-MeV range we use the constraints summarized in
Appendix A of Ref. Panci:2022wlc , where the strongest bounds are set by
different X-rays and low energy gamma rays line searches: Chandra
Watson:2011dw ; Horiuchi:2013noa , Newton-XMM Foster:2022ajl , NuStar
Perez:2016tcq ; Roach:2019ctw ; Ng:2019gch ; Roach:2022lgo , and INTEGRAL
Laha:2020ivk . For heavier masses $1\,\text{MeV}\lesssim m_{a}\lesssim 1000$
MeV the most stringent limits on the $a\to\gamma\gamma$ decay rate come from
COMPTEL and EGRET, and we take the constraints presented in Ref. Essig:2013goa
. Further limits are provided by the optical depth since recombination, which
is measured by the Planck collaboration Planck:2018vyg . Fast DM decay into
photons would significantly modify the fraction of free electrons after
reionization and, consequently, would attenuate the small-scale acoustic peaks
of the CMB power spectrum. The model-independent bounds for the optical depth
can be found in Ref. Cirelli:2009bb ; Liu:2016cnk . These limits constrain
rates of the order $\tau_{\gamma\gamma}\approx 10^{24}$ sec, and are therefore
less constraining than the X-ray telescopes in the relevant parameter region
(cf. Fig. 2). Various future X-ray missions are designed to further extend the
limits, and we use the optimistic projections collected in Ref. Panci:2022wlc
for GECCO Coogan:2021rez , THESEUS Thorpe-Morgan:2020rwc  and Athena
Neronov:2015kca ; Dekker:2021bos ; Ando:2021fhj .

(a)

(b)

(c)

(d)

(e)

(f)

Figure 1:  Diagrams contributing to axion production from flavor-violating
quark and gluon scatterings.

###  II.2 Dark Matter Production

For sufficiently large decay constants, $f_{a}\gtrsim 10^{8}{\rm\,GeV}$, the
axion was never in thermal equilibrium with the SM bath. Thermal axions are
then produced via $2\to 2$ scattering and decay processes of quarks in the
thermal bath, which allows to explain the observed DM relic density through
the freeze-in mechanism Hall:2009bx . The total relic axion abundance is given
by $\Omega_{a}h^{2}=\Omega_{a}h^{2}|_{\rm dec}+\Omega_{a}h^{2}|_{\rm scatt}$,
where $\Omega_{a}h^{2}|_{\rm dec}$ and $\Omega_{a}h^{2}|_{\rm scatt}$ are the
contributions from (flavor-violating) quark decays $q_{i}\to q_{j}a$ and
flavor-diagonal quark scattering processes $q_{i}g(\gamma)\to q_{i}a$ and
$q_{i}\overline{q}_{i}\to g(\gamma)\,a$, respectively (we will comment on
flavor-violating scattering processes in Section II.1). The corresponding
cross-sections read

$\displaystyle\sigma_{q_{i}\gamma\to q_{i}a}$ $\displaystyle=$
$\displaystyle\frac{\alpha_{\rm
em}Q_{i}^{2}}{8f_{a}^{2}}|C^{A}_{q_{i}q_{i}}|^{2}\frac{x\left(-2\ln
x-3+4x-x^{2}\right)}{1-x}\,,$
$\displaystyle\sigma_{q_{i}\overline{q}_{i}\to\gamma a}$ $\displaystyle=$
$\displaystyle\frac{\alpha_{\rm
em}Q_{i}^{2}}{f_{a}^{2}}|C^{A}_{q_{i}q_{i}}|^{2}\frac{x\tanh^{-1}(\sqrt{1-4x})}{1-4x}\,,$
(8)

where $x=m_{q_{i}}^{2}/s$ and $Q_{i}$ is the electric charge of $q_{i}$. The
corresponding gluon scattering processes are obtained from these results by
replacing $\alpha_{\rm em}Q_{i}^{2}\to\alpha_{s}/6$ in $\sigma_{q_{i}\gamma\to
q_{i}a}$ and $\alpha_{\rm em}Q_{i}^{2}\to 4\alpha_{s}/9$ in
$\sigma_{q_{i}\overline{q}_{i}\to\gamma a}$, in agreement with e.g. Ref.
Arias-Aragon:2020shv . The decay rate is given by

$\displaystyle\Gamma_{q_{i}\to q_{j}a}$ $\displaystyle=$
$\displaystyle\frac{m_{q_{i}}^{3}}{64\pi
f_{a}^{2}}\left|C_{q_{i}q_{j}}\right|^{2}\left(1-\frac{m_{q_{j}}^{2}}{m_{q_{i}}^{2}}\right)^{3}\,,$
(9)

where
$C_{q_{i}q_{j}}\equiv\sqrt{|C^{V}_{q_{i}q_{j}}|^{2}+|C^{A}_{q_{i}q_{j}}|^{2}}$
and we have neglected the ALP mass. Following Refs. Hall:2009bx ;
Belfatto:2021ats , one can use these results to derive analytical estimates
for the corresponding contributions to the freeze-in abundance, assuming that
the effective number of relativistic degrees of freedom in the SM bath is
approximately constant and that the axion production takes place during
radiation domination. Under these assumptions one obtains, including charge
multiplicities (cf, Appendix C in Ref. Badziak:2024szg )

$\displaystyle\Omega_{a}h^{2}|_{\rm dec}$ $\displaystyle\approx$
$\displaystyle
0.12\left(\frac{mx_{a}}{0.1{\rm\,MeV}}\right)\left(\frac{9.7\times
10^{9}\text{GeV}}{f_{a}/C_{q_{i}q_{j}}}\right)^{2}\left(\frac{m_{q_{i}}}{\text{GeV}}\right)\left(\frac{70}{g_{*}(m_{q_{i}})}\right)^{3/2}\qquad\mbox{for
decays}\,,$ (10) $\displaystyle\Omega_{a}h^{2}|_{\rm scatt}$
$\displaystyle\approx$ $\displaystyle
0.12\left(\frac{m_{a}}{0.1{\rm\,MeV}}\right)\left(\frac{1.4\times
10^{10}\text{GeV}}{f_{a}/C^{A}_{q_{i}q_{i}}}\right)^{2}\left(\frac{m_{q_{i}}}{\text{GeV}}\right)\left(\frac{70}{g_{*}(m_{q_{i}})}\right)^{3/2}\left(\frac{\alpha_{s}(m_{q_{i}})}{0.48}\right)\qquad\mbox{for
scattering}\,,$ (11)

where we have omitted the sub-dominant contribution from photon scattering.

It is clear from Eqs. (10) and (11) that the scattering contribution is only
slightly smaller than the contribution from quark decays, as a result of the
large size of the strong coupling close to the GeV scale. This also implies
that omitting higher-order QCD corrections is not a good approximation, so we
consider our leading-order results to be valid only up to ${\cal O}(1)$
corrections, which however only has a mild impact on the relevant model
parameter $f_{a}$. Keeping in mind this uncertainty, we can still obtain more
accurate expressions by solving the Boltzmann equation numerically, which
leads to the results presented in Section III. We stress that for this
procedure we only use the temperature dependence of energy and entropy degrees
of freedom $g_{(s)*}(T)$, and unlike Ref. DEramo:2021usm  we neither consider
thermal masses nor flavor off-diagonal scattering processes. Indeed both
effects represent only a subset of the full (and unknown) NLO corrections to
the leading order effects to which we restrict here, as explained in more
detail in Section II.1.

In addition to the purely IR contribution to the DM abundance discussed above,
freeze-in scenarios are potentially sensitive also to processes that are
dominated by high temperatures. In particular, at energies above the
electroweak scale, we have to take into account also operators like

$\displaystyle\mathcal{L}_{\rm
eff}=-C^{A}_{q_{i}q_{j}}\frac{ia}{f_{a}}\frac{m_{q_{i}}}{v}H\overline{Q}_{i}q_{Rj}\,,$
(12)

where $Q_{i}$ ($H$) denotes the quark (Higgs) doublet field. This operator can
be obtained from Eq. (1) upon integrating by parts and using the equations of
motion in the unbroken phase, and for simplicity we have set
$C^{V}_{q_{i}q_{j}}=C^{A}_{q_{i}q_{j}}$. This gives rise to scattering
processes like $\overline{q}_{i}q_{j}\to ha$, which lead to axion production
rates that are UV sensitive and thus depend on the reheating temperature
$T_{R}$. The corresponding UV contribution to the relic abundance can be
related to the decay contribution as Hall:2009bx

$\displaystyle\Omega_{a}h^{2}|_{\rm
UV}\approx\frac{m_{q_{i}}T_{R}}{3\pi^{3}v^{2}}\times\Omega_{a}h^{2}|_{q_{i}\to
q_{j}a}\,.$ (13)

While one could take into account such UV sensitive contributions on the price
of introducing $T_{R}$ as an additional parameter of the models, here we want
to stick to the minimal number of parameters and thus take $T_{R}$
sufficiently small such that the IR contributions always dominate the relic
abundance. As we will see below, this procedure also suppresses the
misalignment contribution, which also depends on additional parameters (the
original misalignment angle). Hence, we establish an upper bound on $T_{R}$ by
requiring that the axion abundance generated from UV sensitive processes is
smaller than the one from decays, giving
$T_{R}<3\pi^{3}v^{2}/m_{q_{i}}=3\times
10^{6}{\rm\,GeV}({\rm\,GeV}/m_{q_{i}})$. Equivalently we can consider an upper
bound on the Hubble parameter at reheating

$\displaystyle
H_{R}<11{\rm\,keV}\left(\frac{{\rm\,GeV}}{m_{q_{i}}}\right)^{2}\,,$ (14)

where have assumed $g_{*}(T_{R})\approx 106.75$.

We can now discuss possible sources of non-thermal production. The most
relevant is the misalignment mechanism for ALPs Arias:2012az ; Blinov:2019rhb
. Also this contribution depends on the reheating temperature, as the onset of
axion oscillations (defined by111For the numerical values below we have used
$m_{a}=1.6H(Tosc)$ as suggested in Ref. Blinov:2019rhb . $m_{a}\simeq H$)
occurs prior to reheating in the axion mass range under consideration. Todays
misalignment abundance is then suppressed due to the dilution that occurred
during an initial period of matter domination222If inflation ends in a period
of kination instead the misalignment contribution would be enhanced. that took
place between the onset of oscillations and $T_{R}$. The resulting ALP
abundance in terms of the misalignment angle $\theta_{0}$ is then given by
Blinov:2019rhb ; Visinelli:2009kt ; Arias:2021rer

$\displaystyle\Omega_{a}h^{2}|_{\rm mis}\approx 4\times
10^{-3}\left(\frac{H_{R}}{11{\rm\,keV}}\right)^{1/2}\left(\frac{f_{a}\theta_{0}}{10^{10}\,\text{GeV}}\right)^{2}.$
(15)

We notice that the misalignment contribution can be somewhat larger than in
models where the ALP is only coupled to leptons Panci:2022wlc , as heavy
quarks require sizable $f_{a}$ scales. Still, misalignment production is never
relevant in the interesting parameter region compatible with astrophysical
bounds, as we will discuss in Sec. III.

###  II.3 Flavor-violating Scattering Processes and Infrared Finiteness

In addition to $q_{i}\rightarrow q_{j}a$ decays, axion production via freeze-
in is also affected by flavor-changing quark and gluon scatterings Arias-
Aragon:2020shv ; DEramo:2021usm ; DEramo:2023asj . For the special case of
$b\bar{q}\rightarrow ag$, $bg\rightarrow aq$, and $qg\rightarrow ab$
reactions, with $q=s,d$, the diagrams are shown in Fig. 1. The squared
amplitude for any of these three processes depends both on the bottom and the
light quark masses. Neglecting the latter, an infrared divergence occurs,
while using non-zero light-quark mass may still lead to an unphysical
enhancement of the cross section due to the sizeable quark-mass hierarchy. Let
us consider the square of the diagram in Fig. 1a. If $q$ is massless and its
momentum is collinear to the momentum of the gluon, a singularity occurs. In
order to isolate the divergent terms, we assign the light quark a small mass
$m_{q}$ as a regulator, and take the limit $m_{q}\to 0$ when possible. The
square of the diagram summed over spins and integrated over the final state
momenta can be expressed as a unitary cut of the forward-scattering diagram

$\displaystyle\includegraphics[scale={1},valign={c}]{math1.pdf}\propto\hskip
2.84526pt2s+\frac{4m^{4}_{b}}{s}+\frac{2sm^{2}_{b}}{s-m^{2}_{b}}\ln\frac{(s-m^{2}_{b})^{2}}{sm^{2}_{q}}\,,$
(16)

where $s=(p_{b\vphantom{\bar{q}}}+p_{\bar{q}})^{2}$ is the squared centre-of-
mass energy. A logarithmic divergence occurs in the last term, and the same
happens in all contributions involving the diagrams in Figs. 1a and 1c. To
deal with these singularities, we follow the procedure introduced in Ref.
Racker:2018tzw  based on the Kinoshita-Lee-Nauenberg (KLM) theorem
Kinoshita:1962ur ; Lee:1964is ; Frye:2018xjj . The forward-scattering diagram
in Eq. (16) allows for two other unitary cuts corresponding to so-called
anomalous thresholds Mandelstam:1960zz ; Cutkosky:1961 ; Goddard:1969ci ;
Hannesdottir:2022bmo , which evaluate to

$\displaystyle\includegraphics[scale={1},valign={c}]{math2.pdf}+\includegraphics[scale={1},valign={c}]{math3.pdf}\propto\hskip
2.84526pt-4m^{2}_{b}-\frac{4m^{4}_{b}}{s}-\frac{2sm^{2}_{b}}{s-m^{2}_{b}}\ln\frac{(s-m^{2}_{b})^{2}m^{2}_{b}}{s^{2}m^{2}_{q}}\,,$
(17)

yielding a finite result for $m_{q}\rightarrow 0$ when added to Eq. (16). We
apply this procedure whenever an infrared divergence occurs in the total cross
sections of the processes listed in Fig. 1, which allows us to obtain well-
defined expressions. Still, another sort of singularity persists.

The square of the diagram in Fig. 1d leads to a finite total cross section,
but diverges for small gluon energy in thermal averaging. The problem has been
resolved in Ref. Czarnecki:2011mr  for electromagnetic corrections to charged
particle decays in a thermal medium. In our case, we must include the gluon-
induced bottom thermal mass and wave-function renormalization factor in the
$b\rightarrow qa$ decay. Furthermore, in analogy to Eq. (9) of Ref.
Czarnecki:2011mr , we add the temperature-dependent part of the gluon Bose-
enhancement factor in the $b\rightarrow qag$ decay. As a result, the
$s$-channel $bg\rightarrow aq$ cross section is found to be completely
canceled for values of $s<2m^{2}_{b}$. Putting everything together, the ratio
of temperature-dependent axion production rates from scattering $\gamma_{S}$
and decays $\gamma_{D}$ is finite in the limit $m_{q}\to 0$ and given by

$\displaystyle\frac{\gamma_{S}}{\gamma_{D}}=$
$\displaystyle\frac{\alpha_{s}}{\pi}\int\displaylimits^{\infty}_{1}d\tilde{s}\frac{\tilde{s}-1}{\sqrt{\tilde{s}}}\frac{K_{1}\big{(}x_{b}\sqrt{\tilde{s}}\big{)}}{K_{1}\big{(}x_{b}\big{)}}\bigg{\\{}-\tilde{s}-4-\frac{3}{\tilde{s}}+\frac{12}{\tilde{s}-1}$
(18)
$\displaystyle+\bigg{[}8-\frac{8}{\tilde{s}}+\frac{10}{\tilde{s}-1}+\frac{4}{(\tilde{s}-1)^{2}}\bigg{]}\ln\tilde{s}+\bigg{[}\tilde{s}+2-\frac{8}{\tilde{s}-1}+\frac{4}{(\tilde{s}-1)^{2}}\bigg{]}\theta(\tilde{s}-2)\bigg{\\}}\,,$

with $\tilde{s}=s/m^{2}_{b}$ and $x_{b}=m_{b}/T$.

Axion production via freeze-in is dominated by temperatures slightly below the
heavy quark mass. From Eq. (18), taking $x_{b}=3$ leads to
$\gamma_{S}/\gamma_{D}\approx 0.7$. This contrasts with the findings of Ref.
DEramo:2021usm , where thermal masses are used in the leading order cross-
sections to handle divergencies. This leads to an enhancement of the
scattering production by two orders of magnitude relative to decay at the
relevant temperature slightly below the heavy quark mass. In our procedure,
this enhancement is cancelled by the contributions of the anomalous
thresholds, as in Eq. (17). Those are indeed related to thermal-mass effects,
but instead of the scattering, they enter the leading-order $b\rightarrow qa$
decay kinematics through mass-derivative relations (see Eq. (35) in Ref.
Blazek:2021gmw  for an example). However, we note that our approximation of
thermal effects is not complete and only represents a minimal set of
contributions needed for infrared finiteness. A more complete treatment of
thermal corrections may be considered in future work.

###  II.4 Warm Dark Matter

Soon after its production, DM free-streams and suppresses the primordial
fluctuations related to the matter power spectrum. DM free-streaming leaves
its footprints on large-scale structures and can be constrained by looking at
the absorption features of the spectra of distant quasars through the
Lyman-$\alpha$ forest (Ly-$\alpha$) Boyarsky:2008xj . In particular, one can
set a “warmness bound” on the DM mass to avoid large free-streaming
Viel:2013fqw ; Baur:2015jsy ; Irsic:2017ixq . The Ly-$\alpha$ limits have been
recasted for different freeze-in processes by computing the exact DM velocity
distribution, which results in the “Warm Dark Matter” (WDM) constraint
DEramo:2020gpr ; Ballesteros:2020adh ; Decant:2021mhj

$\displaystyle m_{a}\gtrsim 0.01{\rm\,MeV}\left(\frac{m_{\rm
WDM}}{3.5\,\text{keV}}\right)^{4/3}\left(\frac{70}{g^{*}(m_{q})}\right)^{1/3}\,,$
(19)

where $m_{\rm WDM}\approx 3.5{\rm\,keV}$ or $5.3{\rm\,keV}$ for the
conservative and stringent bounds, respectively.

###  II.5 Other Astrophysical Bounds

Sufficiently light ALPs coupled to SM fermions can efficiently extract energy
from stellar objects and are subject to limits from star cooling
Raffelt:1996wa . Quark couplings induce axion couplings to nucleons, which
allow for efficient axion production in hot stellar plasmas, such as in the
proto-neutron star formed during core-collapse supernovae. Sufficiently light
axions ($m_{a}\lesssim 100{\rm\,MeV}$) would extract energy from the proto-
neutron star, which is constrained by the usual energy loss argument for
SN1987A Raffelt:1996wa . Lighter axions ($m_{a}\lesssim 0.4{\rm\,keV}$) are
also constrained by measurements of the White Dwarf (WD) luminosity function,
which primarily limits electron couplings at the order of $f_{a}/C_{e}\geq
2.5\times 10^{9}{\rm\,GeV}$ MillerBertolami:2014rka , but due to
Renormalization Group evolution this also puts constraints on axion couplings
to top quarks Feng:1997tn  at the level MartinCamalich:2020dfe

$\displaystyle f_{a}/C_{t}\geq 1.7\times 10^{9}{\rm\,GeV}\,,$ (20)

where we ignored the mild logarithmic dependence on the UV scale by setting
$f_{a}=10^{10}{\rm\,GeV}$ for simplicity.

In order to extract the resulting SN1987A limits on axion-quark couplings, we
match the Lagrangian in Eq. (1) to the axion-nucleon effective Lagrangian.
Following Ref. GrillidiCortona:2015jxo , we obtain the following effective
Lagrangian in the non-relativistic limit, which should be reliable as long as
the ALP mass and the relevant energies are smaller than the QCD mass gap
$\Delta\approx 100$ MeV

$\displaystyle\mathcal{L}_{aN}$
$\displaystyle=\overline{N}v^{\mu}\partial_{\mu}N+\frac{\partial_{\mu}a}{f_{a}}\frac{C_{u}-C_{d}}{2}\Delta_{u-d}\overline{N}S^{\mu}\sigma^{3}N$
(21)
$\displaystyle+\frac{\partial_{\mu}a}{f_{a}}\left[\frac{C_{u}+C_{d}}{2}\Delta_{u+d}+\sum_{q=s,c,b,t}C_{q}\Delta
q\right]\overline{N}S^{\mu}N\,,$

where $N=(p,n)$ is the nucleon isospin doublet, $v^{\mu}$ is the four-velocity
of the nucleon, $2S^{\mu}\equiv\gamma^{\mu}\gamma^{5}$ is the spin operator
and $\Delta_{u\pm d}\equiv\Delta u\pm\Delta d$. The coefficients $\Delta
q,q=u,d,s,c,b,t$ are extracted from lattice QCD studies and low-energy
experiments, and can be found in Ref. GrillidiCortona:2015jxo . Here we use
the recent analysis in Ref. Badziak:2023fsc , giving

$\displaystyle\begin{split}&C_{p}\approx
0.82\,C_{u}-0.45\,C_{d}-0.052\,C_{s}\,,\\\ &C_{n}\approx
0.82\,C_{d}-0.45\,C_{u}-0.052\,C_{s}\,,\\\ \end{split}$ (22)

where we have neglected the contributions from heavy quarks. Bounds on these
couplings can be obtained from the burst duration of the neutrino emission of
SN1987A, which yields Carenza:2019pxu

$\displaystyle 0.61g_{ap}^{2}+g_{an}^{2}+0.53g_{an}g_{ap}<8.26\times
10^{-19}\,,$ (23)

where $g_{ai}\equiv C_{i}m_{i}/f_{a}$. These constraints are roughly
comparable to limits that can be derived from observations of neutron star
cooling rates Buschmann:2021juv , and give for $C_{p}\approx C_{n}\approx
C_{N}$

$\displaystyle f_{a}/C_{N}\gtrsim 1.5\times 10^{9}{\rm\,GeV}\,.$ (24)

Figure 2: Parameter space in the $(m_{a},\alpha)$ plane for selected benchmark
models defined in Section II, which reproduce the observed DM relic abundance
everywhere. The $tu$ scenario (left panel) is representative for the $tc$
scenario, the $bs$ scenario (right panel) for the $bd$ and $cu$ scenarios. The
angle $\alpha$ controls the ratio of flavor-diagonal to flavor-off-diagonal
couplings, see Eq. (3). The shaded blue region shows CMB constraints on the
axion lifetime, while the shaded gray bound indicates to the region excluded
by X-ray telescopes. The dark gray region is excluded by SN (WD) cooling,
while the green shaded region is probed by collider searches. The stringent
WDM bound (the conservative is weaker by a factor 1.7) is denoted by a dashed
black line along with the corresponding target branching ratio for $B\to Ka$,
and green (orange) dashed lines indicate the prospective limits from future
laboratory searches (X-ray telescopes). The upper axis indicates the value of
$f_{a}$ needed to reproduce the relic abundance for $\alpha=0$, with a mild
dependence on $\alpha$.

##  III Results

It is clear from Eqs. (10) and (11) that for heavy quarks
$m_{q_{i}}\gtrsim{\rm\,GeV}$ and axion masses $m_{a}\sim 0.1{\rm\,MeV}$
satisfying the WDM bound in Eq. (19), the observed DM relic abundance can be
obtained for $f_{a}\sim 10^{10}{\rm\,GeV}$, while respecting the limits from
X-ray telescopes in Eq. (7) and the supernova bounds in Eq. (24). This also
implies that axion couplings to light quarks, i.e., the $sd$ scenario, are not
viable, since the strong constraints from $K^{+}\to\pi^{+}+{\rm inv.}$
searches NA62:2021zjw ; Goudzovski:2022vbt  essentially exclude the whole
parameter space, as the relic abundance require too low $f_{a}$ values. In
contrast freeze-in via heavy-quark couplings ($c,b,t$) gives larger values of
$f_{a}$, and there is not much difference between axion production via flavor-
diagonal or flavor-violating couplings, as $\alpha_{s}$ is sizable333This is
in contrast to the $\mu e$ scenario considered in Ref. Panci:2022wlc , where
axion production via flavor-diagonal scattering is suppressed by $\alpha_{\rm
em}$, and requires very low values of $f_{a}$ that are already excluded by
X-ray searches. These limits disappear in the limit where axion couplings are
mainly flavor-violating, which provides a scenario compatible with present
laboratory searches for $\mu\to e+{\rm invis.}$ and in the reach of near-
future experimental proposals.. It turns out that all five scenarios in the
first class discussed in Section II are indeed viable for all values of
$\alpha$, which controls the ratio of diagonal to off-diagonal couplings. This
is because for values of $m_{a}$ that respect the WDM bound, not only the
constraints on flavor-diagonal couplings from star cooling are satisfied, but
also the stringent laboratory limits on flavor-violating decays with missing
energy are respected, even taking into account near-future projections. This
regards $D\to\pi+{\rm inv.}$ searches at CLEO CLEO:2008ffk ;
MartinCamalich:2020dfe  ($cu$ scenario), $B\to\pi+{\rm inv.}$ searches at
BaBar BaBar:2004xlo ; MartinCamalich:2020dfe  ($bd$ scenario) and $B\to K+{\rm
inv.}$ searches at Belle II Belle-II:2023esi  and BaBar BaBar:2013npw  ($bs$
scenario), where we used the combined limit on the two-body decay recently
provided in Ref. Altmannshofer:2023hkn , ${\rm BR}(B\to Ka)<8.0\times 10^{-6}$
at 95% CL for $m_{a}\ll 100{\rm\,MeV}$. Note that the $tc$ and $tu$ scenarios
are constrained only mildly by WD cooling (Eq. (20)) and SN1987A (Eq. (24)),
and SM loop contributions to $K\to\pi a$ are absent in these cases as the
axion only couples to RH quarks444Otherwise one would obtain constraints on
$f_{a}$ of the order of few$\times 10^{8}{\rm\,GeV}$ MartinCamalich:2020dfe ..

Figure 3: Parameter space in the $(m_{a},X)$ plane for the CKM${}_{Q_{L}}$
(left panel) and CKM${}_{d_{R}}$ (right panel) benchmark models defined in
Section II, which reproduce the observed DM relic abundance everywhere. The
parameter $X$ controls the ratio of PQ charges for LH quarks (CKM${}_{Q_{L}}$)
or RH down-quarks (CKM${}_{d_{R}}$), which are rotated to the quark mass basis
by CKM rotations, see Eq. (3). The shaded blue region shows CMB constraints on
the axion lifetime, while the shaded gray bound indicates to the region
excluded by X-ray telescopes. The green shaded region is excluded by $K\to\pi
a$ searches at NA62. The stringent WDM bound (the conservative is weaker by a
factor 1.7) is denoted by a dashed black line along with the corresponding
target branching ratio for $K\to\pi a$, and green (orange) dashed lines
indicate the prospective limits from future laboratory searches (X-ray
telescopes). The upper axis indicates the value of $f_{a}$ needed to reproduce
the relic abundance for $X=1$, with a mild dependence on $X$.

Thus we obtain five simple benchmark models that generate the DM relic
abundance and are compatible with all present constraints. In Fig. 2 we
display the 2-dimensional parameter space for two of such models, the $tu$
scenario and the $bs$ scenario. Analogous figures for the $tc$ scenario and
the $bd$ and $cu$ scenarios are not shown since they are very similar to the
$tu$ and $bs$ scenario, respectively. All scenarios have in common that even
future laboratory searches for two-body flavor-violating decays with missing
energy will not probe the interesting region of axion masses satisfying the
WDM bound, which is right of the vertical dashed black line. For the $bs$
scenario we have indicated the maximal size of the $B\to Ka$ target branching
ratio compatible with the WDM bound, which is of the order of ${\rm BR}(B\to
Ka)_{\rm target}\approx 1\times 10^{-9}$, clearly beyond the reach of running
or near-future $B$-factories MartinCamalich:2020dfe  indicated by the green
dashed line. The same conclusions are valid for $D\to\pi$, with ${\rm
BR}(D\to\pi a)_{\rm target}\approx 2\times 10^{-10}$, and $B\to\pi$
transitions with ${\rm BR}(B\to\pi a)_{\rm target}\approx 6\times 10^{-10}$,
such that all scenarios will be tested only by future X-ray telescopes. Note
however that colliders are in principle better suited to probe the remaining
parameter space as compared to X-ray line searches, as the constrained axion
mass scales rather weakly with the axion decay rate into photons,
$m_{a}\propto\Gamma_{a\to\gamma\gamma}^{1/6}$ (Eq. (7)), but strongly with
flavor-violating decay rates $m_{a}\propto\Gamma_{q_{i}\to q_{j}a}^{1/2}$ (Eq.
(9)), after fixing $f_{a}^{2}\propto m_{a}$ with the relic abundance (Eq.
(10)).

We finally discuss benchmark models where present collider constraints exceed
the WDM bound, such that the remaining parameter space will be complementary
probed by precision flavor experiments and X-ray telescopes. This is the case
for the other two benchmark scenarios discussed in Section II, where the PQ
charge is taken as a free parameter and the rotation to the quark mass basis
is fixed by the CKM matrix. The resulting parameter space is shown in Fig. 3
for the CKM${}_{Q_{L}}$ and the CKM${}_{d_{R}}$ scenario, which reproduce the
observed DM relic abundance for the indicated values for $m_{a}$ and the PQ
charge $X$. The dominant contribution to axion production comes from
unsuppressed processes involving the heaviest quarks, which is $tt$ scattering
in the CKM${}_{Q_{L}}$ scenario and $bb$ scattering in the CKM${}_{d_{R}}$
scenario, since quark mixing in the CKM is small and cannot compensate the
mild $\alpha_{s}$ suppression in scatterings compared to decays. Still the CKM
involves a rather large rotation in the $sd$ sector of order $\lambda\approx
0.23$, which induces a sizable coupling of the axion to $sd$ quarks in both
scenarios for generic values of $X$, unless the first two generation have the
same PQ charge, i.e. $X=1$, leading to an approximate $SU(2)$ symmetry and
$C^{V}_{sd}$ involves additional CKM suppression555See Ref. Linster:2018avp
for a motivated scenario where this situation arises by identifying PQ as a
subgroup of a horizontal $U(2)$ symmetry explaining Yukawa hierarchies.. Close
to this value the stringent limits on $K\to\pi a$ from NA62 NA62:2021zjw  are
relaxed, which otherwise give constraints of order $f_{a}/C^{V}_{sd}>4\times
10^{11}{\rm\,GeV}$ Goudzovski:2022vbt . For $|X-1|\gtrsim 0.44$
(CKM${}_{d_{R}}$) or $|X-1|\gtrsim 0.24$ (CKM${}_{Q_{L}}$) the resulting
limits exceed the WDM bound, so that NA62 will probe the remaining parameter
space in the near future complementary to future X-ray telescopes.
Interestingly, the parameter space of the CKM${}_{d_{R}}$ scenario will be
almost entirely probed by $K\to\pi a$ and $a\to\gamma\gamma$ searches, leaving
only a narrow region between axion masses $10\lesssim m_{a}\lesssim
100{\rm\,keV}$ and PQ charge $|X-1|\lesssim 0.15$.

##  IV Conclusions

To summarize, we have explored the production of axion DM from decays and
scatterings of heavy quarks via thermal freeze-in. This gives rise to very
simple scenarios with few parameters able to explain the observed DM
abundance, which are subject to various constraints from precision flavor
experiments, star cooling, X-ray telescopes and structure formation. Similar
to the lepton case explored in Ref. Panci:2022wlc , we have focussed on two
classes of models with only two parameters after fixing the axion decay
constant to values that reproduce the observed DM relic abundance. Apart from
the axion mass, the free parameter is the ratio of flavor-diagonal couplings
to flavor-violating couplings in the first class (effective 2-flavor
scenarios), and the overall PQ charge in the second, with flavor violation
controlled by the CKM matrix.

Compared to the case of an ALP coupled to electrons Panci:2022wlc , in the
quark scenarios the axion decay rate into photons is additionally suppressed
by at least a factor $m_{e}^{4}/m_{\pi}^{4}$, which enhances axion stability
and eases constraints on flavor-diagonal couplings from X-ray line searches.
As axion production from quark scattering is only mildly suppressed with
respect to quark decays as a result of large values of $\alpha_{s}$ close to
the GeV scale, we find that there is not much difference between scenarios
with flavor-violating coupling and flavor-diagonal couplings of the same size,
in stark contrast to LFV models Panci:2022wlc . This also implies that next-
to-leading order corrections to axion production are sizable, which we have
calculated here for the first time for the case of flavor-violating $2\to 2$
scattering processes. We showed that in concordance with the KLN theorem IR
divergencies in these processes are cancelled by taking into account
contributions from anomalous thresholds (partially related to thermal
corrections), at least for energies below the heaviest quark involved. A more
complete analysis of NLO corrections is left for future work.

The main results of our analysis are summarized in Fig. 2 for the effective
2-flavor model and the CKM scenarios in Fig. 3. The allowed parameter space of
all 2-flavor models have similar shapes (thus we only show the $tu$ and $bs$
scenarios as representatives), and are viable except for the $sd$ model, which
is essentially ruled out by present $K\to\pi a$ constraints. These scenarios
will only be probed by future X-ray telescopes, as the sensitivities of future
flavor factories will still be weaker than the constraints on Warm Dark
Matter. On the other hand in the CKM scenarios the strongest limits in the low
axion mass regime arise from searches for $K\to\pi a$ at NA62. The expected
sensitivity together with X-ray line searches will allow to probe large
portions of the remaining parameter space, giving excellent prospects to
explore a very simple class of axion DM models at the high-intensity frontier.

###### Acknowledgements.

We would like to thank Francesco D’Eramo, Kirill Melnikov, Uli Nierste and
Diego Redigolo for useful discussions. Peter Maták and Zuzana Šinská were
supported by the Slovak Grant Agency VEGA, project No. 1/0719/23. Peter Maták
also received financial support from Slovak Education Ministry contract No.
0466/2022. The work of Robert Ziegler has received support from the European
Union’s Horizon 2020 research and innovation programme under the Marie
Skłodowska-Curie grant agreement No 860881-HIDDeN and is partially supported
by project B3a and C3b and of the DFG-funded Collaborative Research Center
TRR257 “Particle Physics Phenomenology after the Higgs Discovery”. The
research conducted by Mohammad Aghaie, Giovanni Armando, Alessandro Dondarini,
Angela Conaci and Paolo Panci receives partial funding from the European
Union–Next generation EU (through Progetti di Ricerca di Interesse Nazionale
(PRIN) Grant No. 202289JEW4).

##  Appendix A Axion-Photon Coupling

In this appendix we present the details of the light quark contribution to the
axion-photon coupling in Eq. (II.1) due to axion mixing with $\pi,\eta$ and
$\eta^{\prime}$ in leading order chiral perturbation theory ($\chi$PT). Our
analysis complements the results of Ref. Ertas:2020xcc , where the
contribution to axion-meson mixing from the axion-gluon coupling was
calculated. Here instead we provide the contribution from axion couplings to
all three light quarks.

We start by matching the Lagrangian in Eq.(1) to 3-flavor $\chi$PT. After
integrating out the heavy quarks, we define the effective axion couplings as
diagonal $3\times 3$ matrices $k_{R,L}=1/2\,{\rm diag}(C^{V}_{uu}\pm
C^{A}_{uu},C^{V}_{dd}\pm C^{A}_{dd},C^{V}_{ss}\pm C^{A}_{ss})$, and the
Lagrangian reads

$\displaystyle\mathcal{L}_{\text{light}}$
$\displaystyle=\frac{1}{2}(\partial_{\mu}a)^{2}-\frac{m_{a}^{2}}{2}a^{2}+\bar{\Psi}(i\not{D}-M_{q})\Psi$
(25)
$\displaystyle+\frac{\partial_{\mu}a}{f_{a}}\bar{\Psi}\gamma^{\mu}\left(k_{L}P_{L}+k_{R}P_{R}\right)\Psi\,,$
(26)

where $\Psi\equiv(u,d,s)^{T}$ and $M_{q}={\rm diag}(m_{u},m_{d},m_{s})$. The
chiral Lagrangian is written in terms of the unitary $3\times 3$ matrix
$\Sigma$ containing the Goldstone boson octet $(\pi,K,\eta_{8})$ and the
singlet $\eta_{0}$ as $\Sigma=\exp{\left(i\sqrt{2}\Phi/f_{\pi}\right)}$, where

$\displaystyle\Phi=\begin{pmatrix}\pi^{0}+\frac{\eta_{8}}{\sqrt{3}}&\sqrt{2}\pi^{+}&\sqrt{2}K^{+}\\\
\sqrt{2}\pi^{-}&-\pi^{0}+\frac{\eta_{8}}{\sqrt{3}}&\sqrt{2}K^{0}\\\
\sqrt{2}K^{-}&\sqrt{2}\bar{K}^{0}&-\frac{2}{\sqrt{3}}\eta_{8}\end{pmatrix}+\sqrt{\frac{2}{3}}\eta_{0}\mathbb{1}\,.$
(27)

At leading order the $SU(3)_{\rm L}\times SU(3)_{\rm R}$ symmetry gives

$\displaystyle\mathcal{L}_{\rm\chi PT}$
$\displaystyle=\frac{1}{2}(\partial_{\mu}a)^{2}-\frac{m_{a}^{2}}{2}a^{2}+\frac{f_{\pi}^{2}}{8}\text{Tr}\left[D_{\mu}\Sigma
D^{\mu}\Sigma^{\dagger}\right]$
$\displaystyle+\frac{f_{\pi}^{2}}{4}B_{0}\text{Tr}\left[M_{q}\Sigma^{\dagger}+{\rm
h.c.}\right]-\frac{1}{2}M^{2}_{0}\eta_{0}^{2},$ (28)

where the explicit mass term $M_{0}$ takes into account the explicit breaking
of the anomalous $U(1)_{A}$ symmetry and the covariant derivative reads

$\displaystyle
D_{\mu}\Sigma=\partial_{\mu}\Sigma+ieA_{\mu}\left[Q,\Sigma\right]+i\frac{\partial_{\mu}a}{f_{a}}(k_{L}\Sigma-\Sigma
k_{R})\,,$ (29)

where $Q={\rm diag}(2/3,-1/3,-1/3)$ is the electric charge matrix of the light
quarks.

Notice that the axion enters the chiral Lagrangian only through derivative
terms, since there is no axion-gluon coupling. This gives the kinetic mixing
with the mesons in the diagonal entries of $\Phi$, apart from the usual meson
mass matrix. Defining $\phi=(a,\pi^{0},\eta_{8},\eta_{0})$, one obtains for
the relevant quadratic Lagrangian ${\cal L}\supset
1/2K_{ij}\partial_{\mu}\phi_{i}\partial^{\mu}\phi_{j}-1/2M^{2}_{ij}\phi_{i}\phi_{j}$,
with

$\displaystyle K_{ij}$
$\displaystyle=\delta_{ij}+K_{i1}\delta_{1j}+K_{j1}\delta_{1i}\,.$ (30)

Here

$\displaystyle K_{i1}$
$\displaystyle=-\frac{\epsilon}{2\sqrt{6}}\begin{pmatrix}{\cal
O}(\epsilon^{2})\\\ \sqrt{3}(C_{u}-C_{d})\\\ C_{u}+C_{d}-2C_{s}\\\
\sqrt{2}(C_{u}+C_{d}+C_{s})\end{pmatrix}\ ,\qquad
M^{2}=\begin{pmatrix}m_{a}^{2}&0&0&0\\\
0&2B_{0}\hat{m}&-\frac{B_{0}}{\sqrt{3}}\Delta&-\sqrt{\frac{2}{3}}B_{0}\Delta\\\
0&-\frac{B_{0}}{\sqrt{3}}\Delta&\frac{2}{3}B_{0}(\hat{m}+2m_{s})&\frac{4}{3\sqrt{2}}B_{0}(\hat{m}-m_{s})\\\
0&-\sqrt{\frac{2}{3}}B_{0}\Delta&\frac{4}{3\sqrt{2}}B_{0}(\hat{m}-m_{s})&\frac{2}{3}B_{0}(2\hat{m}+m_{s})+M^{2}_{0}\end{pmatrix}\
,$ (31)

where we have defined $\hat{m}\equiv(m_{u}+m_{d})/2$, $\Delta=m_{d}-m_{u}$ and
$\epsilon\equiv f_{\pi}/f_{a}$.

We continue by taking the isospin limit $\Delta=0$, so that only
$\eta_{8}-\eta_{0}$ mass mixing takes place. The mass basis is thus defined by
a single rotation

$\displaystyle\begin{pmatrix}\eta_{8}\\\ \eta_{0}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\cos{\theta}&\sin{\theta}\\\
-\sin{\theta}&\cos{\theta}\end{pmatrix}\begin{pmatrix}\eta\\\
\eta^{\prime}\end{pmatrix}\,,$ (32)

with the rotation angle $\theta$ given by

$\displaystyle\tan\theta=\frac{4}{3\sqrt{2}}\frac{B_{0}(\hat{m}-m_{s})}{m_{\eta^{\prime}}^{2}-2/3B_{0}(\hat{m}+2m_{s})}\,.$
(33)

In the limit where axion-meson mixing can be neglected, $f_{a}<<f_{\pi}$, the
entries of the meson mass matrix can be obtained as usual, giving in the
isospin limit $m_{\pi}^{2}=2B_{0}\hat{m}$ and
$m_{K}^{2}=B_{0}(\hat{m}+m_{s})$. The explicit $U(1)_{A}$ breaking term
$M_{0}^{2}$ is determined by the $\eta^{\prime}$-mass, so that not only the
mixing angle, but also the $\eta$-mass are predicted, at least in the leading-
order (LO) level we are considering here. We obtain for the $\eta$-mass (in
agreement with e.g. Ref. Leutwyler:1997yr )

$m_{\eta}^{2}=m_{\eta_{8}}^{2}-\frac{8}{9}\frac{(m_{\pi}^{2}-m_{K}^{2})^{2}}{m_{\eta^{\prime}}^{2}-m_{\eta_{8}}^{2}},$
(34)

where $m_{\eta_{8}}^{2}\equiv(4m^{2}_{K}-m_{\pi}^{2})/3\approx 566{\rm\,MeV}$
denotes the $\eta$-mass one obtains using the Gell-Mann-Okubo formula, that
is, decoupling the $\eta^{\prime}$. Instead Eq. (34) gives $m_{\eta}\approx
494$ MeV, which is in mild tension with the measured $m_{\eta}\approx 548$
MeV. It is well known that $\chi$PT at leading order is not adequate to
describe $\eta-\eta^{\prime}$ mixing Georgi:1993jn ; Gerard:2004gx , and
${\cal O}(p^{4})$ give important corrections to mixing angles and masses
Leutwyler:1997yr ; Beisert:2001qb ; Alves:2017avw . For our purposes however
the LO result for the mixing angle suffices, keeping in mind that
uncertainties from higher-order corrections are large. This gives in agreement
with Ref. Leutwyler:1997yr

$\displaystyle\tan\theta=\frac{4}{3\sqrt{2}}\frac{m_{\pi}^{2}-m_{K}^{2}}{m_{\eta^{\prime}}^{2}-m_{\eta_{8}}^{2}}\,,$
(35)

and numerically $\theta\approx-20^{\circ}$. Other methods give values ranging
from $-13^{\circ}$ to $-22^{\circ}$ Beisert:2001qb , so that in the following
we work with the choice $\sin\theta\approx-1/3$, as frequently done in the
literature Aloni:2018vki ; Cheng:2021kjg . The field redefinition in Eq. (32)
modifies kinetic mixing, resulting in a rotation acting on $K_{i1}$ in Eq.
(31), which becomes after setting
$s_{\theta}\equiv\sin\theta\approx-1/3,c_{\theta}\equiv\cos\theta\approx
2\sqrt{2}/3$

$\displaystyle K_{21}$
$\displaystyle\to-\frac{\epsilon}{2\sqrt{6}}\left[(C_{u}+C_{d})(c_{\theta}-\sqrt{2}s_{\theta})-2C_{s}(c_{\theta}+\frac{s_{\theta}}{\sqrt{2}})\right]$
$\displaystyle=-\frac{\epsilon}{2\sqrt{3}}(C_{u}+C_{d}-C_{s})\,,$
$\displaystyle K_{31}$
$\displaystyle\to-\frac{\epsilon}{2\sqrt{6}}\left[(C_{u}+C_{d})(s_{\theta}+\sqrt{2}c_{\theta})-2C_{s}(s_{\theta}-\frac{c_{\theta}}{\sqrt{2}})\right]$
$\displaystyle=-\frac{\epsilon}{2\sqrt{6}}(C_{u}+C_{d}+2C_{s})\,.$ (36)

Finally we canonically normalize kinetic terms, and re-diagonalize the mass
matrix. At linear order in $f_{a}/f_{\pi}$ this is straighforward, and gives
the following relation between the fields in the original basis
$\phi=(a,\pi^{0},\eta_{8},\eta_{0})$ and canonically normalized mass
eigenstates $\phi_{\rm phys}=(a_{\rm phys},\pi^{0}_{\rm phys},\eta_{\rm
phys},\eta^{\prime}_{\rm phys})$

$\displaystyle\pi^{0}$ $\displaystyle\approx\pi^{0}_{\rm
phys}+\epsilon\frac{C_{u}-C_{d}}{2\sqrt{2}}\frac{m_{a}^{2}}{m_{a}^{2}-m_{\pi}^{2}}a_{\rm
phys}\,,$ (37) $\displaystyle\eta_{8}$ $\displaystyle\approx\eta_{\rm
phys}+\epsilon\frac{C_{u}+C_{d}-C_{s}}{2\sqrt{3}}\frac{m_{a}^{2}}{m_{a}^{2}-m_{\eta}^{2}}a_{\rm
phys}\,,$ (38) $\displaystyle\eta_{0}$ $\displaystyle\approx\eta^{\prime}_{\rm
phys}+\epsilon\frac{C_{u}+C_{d}+2C_{s}}{2\sqrt{6}}\frac{m_{a}^{2}}{m_{a}^{2}-m_{\eta^{\prime}}^{2}}a_{\rm
phys}\,.$ (39)

These results allow to compute the light quark contribution to the axion
couplings to photons, which the axion inherits from the meson couplings,
suppressed by mixing. The pseudoscalar couplings to photons induced by the
electromagnetic anomaly read

$\displaystyle\mathcal{L}_{\rm EMA}$
$\displaystyle=\frac{i}{2}\frac{N_{c}\alpha_{\rm
em}}{4\pi}F_{\mu\nu}\widetilde{F}^{\mu\nu}{\rm
Tr}\,[Q^{2}(\log{\Sigma}-\log{\Sigma^{\dagger}})]$
$\displaystyle=-\frac{\alpha_{\rm em}}{4\pi
f_{\pi}}F\widetilde{F}\left(\sqrt{2}\pi^{0}+\sqrt{\frac{2}{3}}\eta_{8}+\frac{4}{\sqrt{3}}\eta_{0}\right)\,.$
(40)

Plugging in the mixing relations in Eq. (37), and matching to the axion-photon
coupling defined as

$\displaystyle\mathcal{L}_{\rm a\gamma\gamma}$
$\displaystyle=C_{\gamma\gamma}\frac{\alpha_{\rm
em}}{4\pi}\frac{a}{f_{a}}F_{\mu\nu}\widetilde{F}^{\mu\nu}\,,$ (41)

we finally obtain the light quark contribution $C^{\rm light}_{\gamma\gamma}$
in Eq. (II.1), which is valid in the limit $m_{a}\ll
m_{\pi,\eta,\eta^{\prime}}$ and taking $s_{\theta}\approx-1/3$.

## References

  * (1) J. L. Feng, T. Moroi, H. Murayama and E. Schnapka, _Third generation familons, b factories, and neutrino cosmology_ , _Phys. Rev. D_ 57 (1998) 5875–5892, [hep-ph/9709411].
  * (2) J. F. Kamenik and C. Smith, _FCNC portals to the dark sector_ , _JHEP_ 03 (2012) 090, [1111.6402].
  * (3) F. Björkeroth, E. J. Chun and S. F. King, _Flavourful Axion Phenomenology_ , _JHEP_ 08 (2018) 117, [1806.00660].
  * (4) J. Martin Camalich, M. Pospelov, P. N. H. Vuong, R. Ziegler and J. Zupan, _Quark Flavor Phenomenology of the QCD Axion_ , _Phys. Rev. D_ 102 (2020) 015023, [2002.04623].
  * (5) L. Calibbi, D. Redigolo, R. Ziegler and J. Zupan, _Looking forward to lepton-flavor-violating ALPs_ , _JHEP_ 09 (2021) 173, [2006.04795].
  * (6) R. Ziegler, _Flavor Probes of Axion Dark Matter_ , _PoS_ DISCRETE2022 (2024) 086, [2303.13353].
  * (7) E. Goudzovski et al., _New Physics Searches at Kaon and Hyperon Factories_ , 2201.07805.
  * (8) Y. Jho, S. Knapen and D. Redigolo, _Lepton-flavor violating axions at MEG II_ , 2203.11222.
  * (9) S. Knapen, K. Langhoff, T. Opferkuch and D. Redigolo, _A Robust Search for Lepton Flavour Violating Axions at Mu3e_ , 2311.17915.
  * (10) R. J. Hill, R. Plestid and J. Zupan, _Searching for new physics at $\mu\rightarrow e$ facilities with $\mu^{+}$ and $\pi^{+}$ decays at rest_, 2310.00043.
  * (11) A. Davidson and K. C. Wali, _MINIMAL FLAVOR UNIFICATION VIA MULTIGENERATIONAL PECCEI-QUINN SYMMETRY_ , _Phys. Rev. Lett._ 48 (1982) 11.
  * (12) F. Wilczek, _Axions and Family Symmetry Breaking_ , _Phys. Rev. Lett._ 49 (1982) 1549–1552.
  * (13) Z. G. Berezhiani and M. Y. Khlopov, _Cosmology of Spontaneously Broken Gauge Family Symmetry_ , _Z. Phys. C_ 49 (1991) 73–78.
  * (14) Y. Ema, K. Hamaguchi, T. Moroi and K. Nakayama, _Flaxion: a minimal extension to solve puzzles in the standard model_ , _JHEP_ 01 (2017) 096, [1612.05492].
  * (15) L. Calibbi, F. Goertz, D. Redigolo, R. Ziegler and J. Zupan, _Minimal axion model from flavor_ , _Phys. Rev. D_ 95 (2017) 095009, [1612.08040].
  * (16) P. Panci, D. Redigolo, T. Schwetz and R. Ziegler, _Axion dark matter from lepton flavor-violating decays_ , _Phys. Lett. B_ 841 (2023) 137919, [2209.03371].
  * (17) A. Czarnecki, M. Kamionkowski, S. K. Lee and K. Melnikov, _Charged-particle decay at finite temperature_ , _Phys. Rev. D_ 85 (Jan, 2012) 025018.
  * (18) M. Linster and R. Ziegler, _A Realistic $U(2)$ Model of Flavor_, _JHEP_ 08 (2018) 058, [1805.07341].
  * (19) M. Bauer, M. Neubert and A. Thamm, _Collider Probes of Axion-Like Particles_ , _JHEP_ 12 (2017) 044, [1708.00443].
  * (20) K. Nakayama, F. Takahashi and T. T. Yanagida, _Anomaly-free flavor models for Nambu–Goldstone bosons and the 3.5keV X-ray line signal_ , _Phys. Lett. B_ 734 (2014) 178–182, [1403.7390].
  * (21) F. Takahashi, M. Yamada and W. Yin, _XENON1T Excess from Anomaly-Free Axionlike Dark Matter and Its Implications for Stellar Cooling Anomaly_ , _Phys. Rev. Lett._ 125 (2020) 161801, [2006.10035].
  * (22) C. Han, M. L. López-Ibáñez, A. Melis, O. Vives and J. M. Yang, _Anomaly-free leptophilic axionlike particle and its flavor violating tests_ , _Phys. Rev. D_ 103 (2021) 035028, [2007.08834].
  * (23) C. Han, M. L. López-Ibáñez, A. Melis, O. Vives and J. M. Yang, _Anomaly-free ALP from non-Abelian flavor symmetry_ , 2203.16376.
  * (24) K. Sakurai and F. Takahashi, _Anomaly-free axion dark matter in three Higgs doublet model and its phenomenological implications_ , _JHEP_ 07 (2022) 124, [2203.17212].
  * (25) C. R. Watson, Z.-Y. Li and N. K. Polley, _Constraining Sterile Neutrino Warm Dark Matter with Chandra Observations of the Andromeda Galaxy_ , _JCAP_ 03 (2012) 018, [1111.4217].
  * (26) S. Horiuchi, P. J. Humphrey, J. Onorbe, K. N. Abazajian, M. Kaplinghat and S. Garrison-Kimmel, _Sterile neutrino dark matter bounds from galaxies of the Local Group_ , _Phys. Rev. D_ 89 (2014) 025017, [1311.0282].
  * (27) J. W. Foster, S. Kumar, B. R. Safdi and Y. Soreq, _Dark Grand Unification in the axiverse: decaying axion dark matter and spontaneous baryogenesis_ , _JHEP_ 12 (2022) 119, [2208.10504].
  * (28) K. Perez, K. C. Y. Ng, J. F. Beacom, C. Hersh, S. Horiuchi and R. Krivonos, _Almost closing the $\nu$MSM sterile neutrino dark matter window with NuSTAR_, _Phys. Rev. D_ 95 (2017) 123002, [1609.00667].
  * (29) B. M. Roach, K. C. Y. Ng, K. Perez, J. F. Beacom, S. Horiuchi, R. Krivonos et al., _NuSTAR Tests of Sterile-Neutrino Dark Matter: New Galactic Bulge Observations and Combined Impact_ , _Phys. Rev. D_ 101 (2020) 103011, [1908.09037].
  * (30) K. C. Y. Ng, B. M. Roach, K. Perez, J. F. Beacom, S. Horiuchi, R. Krivonos et al., _New Constraints on Sterile Neutrino Dark Matter from $NuSTAR$ M31 Observations_, _Phys. Rev. D_ 99 (2019) 083005, [1901.01262].
  * (31) B. M. Roach, S. Rossland, K. C. Y. Ng, K. Perez, J. F. Beacom, B. W. Grefenstette et al., _Long-exposure NuSTAR constraints on decaying dark matter in the Galactic halo_ , _Phys. Rev. D_ 107 (2023) 023009, [2207.04572].
  * (32) R. Laha, J. B. Muñoz and T. R. Slatyer, _INTEGRAL constraints on primordial black holes and particle dark matter_ , _Phys. Rev. D_ 101 (2020) 123514, [2004.00627].
  * (33) R. Essig, E. Kuflik, S. D. McDermott, T. Volansky and K. M. Zurek, _Constraining Light Dark Matter with Diffuse X-Ray and Gamma-Ray Observations_ , _JHEP_ 11 (2013) 193, [1309.4091].
  * (34) Planck collaboration, N. Aghanim et al., _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6, [1807.06209].
  * (35) M. Cirelli, F. Iocco and P. Panci, _Constraints on Dark Matter annihilations from reionization and heating of the intergalactic gas_ , _JCAP_ 10 (2009) 009, [0907.0719].
  * (36) H. Liu, T. R. Slatyer and J. Zavala, _Contributions to cosmic reionization from dark matter annihilation and decay_ , _Phys. Rev. D_ 94 (2016) 063507, [1604.02457].
  * (37) A. Coogan, A. Moiseev, L. Morrison and S. Profumo, _Hunting for Dark Matter and New Physics with (a) GECCO_ , 2101.10370.
  * (38) C. Thorpe-Morgan, D. Malyshev, A. Santangelo, J. Jochum, B. Jäger, M. Sasaki et al., _THESEUS insights into axionlike particles, dark photon, and sterile neutrino dark matter_ , _Phys. Rev. D_ 102 (2020) 123003, [2008.08306].
  * (39) A. Neronov and D. Malyshev, _Toward a full test of the $\nu$MSM sterile neutrino dark matter model with Athena_, _Phys. Rev. D_ 93 (2016) 063518, [1509.02758].
  * (40) A. Dekker, E. Peerbooms, F. Zimmer, K. C. Y. Ng and S. Ando, _Searches for sterile neutrinos and axionlike particles from the Galactic halo with eROSITA_ , _Phys. Rev. D_ 104 (2021) 023021, [2103.13241].
  * (41) S. Ando et al., _Decaying dark matter in dwarf spheroidal galaxies: Prospects for x-ray and gamma-ray telescopes_ , _Phys. Rev. D_ 104 (2021) 023022, [2103.13242].
  * (42) L. J. Hall, K. Jedamzik, J. March-Russell and S. M. West, _Freeze-In Production of FIMP Dark Matter_ , _JHEP_ 03 (2010) 080, [0911.1120].
  * (43) F. Arias-Aragón, F. D’Eramo, R. Z. Ferreira, L. Merlo and A. Notari, _Production of Thermal Axions across the ElectroWeak Phase Transition_ , _JCAP_ 03 (2021) 090, [2012.04736].
  * (44) B. Belfatto, D. Buttazzo, C. Gross, P. Panci, A. Strumia, N. Vignaroli et al., _Dark Matter abundance via thermal decays and leptoquark mediators_ , _JHEP_ 06 (2022) 084, [2111.14808].
  * (45) M. Badziak, K. Harigaya, M. Łukawski and R. Ziegler, _Thermal production of astrophobic axions_ , 2403.05621.
  * (46) F. D’Eramo and S. Yun, _Flavor violating axions in the early Universe_ , _Phys. Rev. D_ 105 (2022) 075002, [2111.12108].
  * (47) P. Arias, D. Cadamuro, M. Goodsell, J. Jaeckel, J. Redondo and A. Ringwald, _WISPy Cold Dark Matter_ , _JCAP_ 06 (2012) 013, [1201.5902].
  * (48) N. Blinov, M. J. Dolan, P. Draper and J. Kozaczuk, _Dark matter targets for axionlike particle searches_ , _Phys. Rev. D_ 100 (2019) 015049, [1905.06952].
  * (49) L. Visinelli and P. Gondolo, _Axion cold dark matter in non-standard cosmologies_ , _Phys. Rev. D_ 81 (2010) 063508, [0912.0015].
  * (50) P. Arias, N. Bernal, D. Karamitros, C. Maldonado, L. Roszkowski and M. Venegas, _New opportunities for axion dark matter searches in nonstandard cosmological models_ , _JCAP_ 11 (2021) 003, [2107.13588].
  * (51) F. D’Eramo, _Thermal Axions: Production Mechanisms and Cosmological Signals_ , _LHEP_ 2023 (2023) 353.
  * (52) J. Racker, _Unitarity and CP violation in leptogenesis at NLO: general considerations and top Yukawa contributions_ , _Journal of High Energy Physics_ 2019 (Feb, 2019) 42, [1811.00280].
  * (53) T. Kinoshita, _Mass singularities of feynman amplitudes_ , _J. Math. Phys._ 3 (1962) 650–677.
  * (54) T. D. Lee and M. Nauenberg, _Degenerate systems and mass singularities_ , _Phys. Rev._ 133 (Mar, 1964) B1549–B1562.
  * (55) C. Frye, H. Hannesdottir, N. Paul, M. D. Schwartz and K. Yan, _Infrared finiteness and forward scattering_ , _Phys. Rev. D_ 99 (Mar, 2019) 056015, [1810.10022].
  * (56) S. Mandelstam, _Unitarity condition below physical thresholds in the normal and anomalous cases_ , _Phys. Rev. Lett._ 4 (Jan, 1960) 84–87.
  * (57) R. E. Cutkosky, _Anomalous thresholds_ , _Rev. Mod. Phys._ 33 (Jul, 1961) 448–455.
  * (58) P. Goddard, _Anomalous threshold singularities in S-matrix theory_ , _Il Nuovo Cimento A (1965-1970)_ 59 (Feb, 1969) 335–355.
  * (59) H. S. Hannesdottir and S. Mizera, _What is the i $\varepsilon$ for the S-matrix?_ SpringerBriefs in Physics. Springer, 1, 2023, 10.1007/978-3-031-18258-7.
  * (60) T. Blažek and P. Maták, _Mass-derivative relations for leptogenesis_ , _The European Physical Journal C_ 82 (Mar, 2022) 214.
  * (61) A. Boyarsky, J. Lesgourgues, O. Ruchayskiy and M. Viel, _Lyman-alpha constraints on warm and on warm-plus-cold dark matter models_ , _JCAP_ 05 (2009) 012, [0812.0010].
  * (62) M. Viel, G. D. Becker, J. S. Bolton and M. G. Haehnelt, _Warm dark matter as a solution to the small scale crisis: New constraints from high redshift Lyman- $\alpha$ forest data_, _Phys. Rev. D_ 88 (2013) 043502, [1306.2314].
  * (63) J. Baur, N. Palanque-Delabrouille, C. Yèche, C. Magneville and M. Viel, _Lyman-alpha Forests cool Warm Dark Matter_ , _JCAP_ 08 (2016) 012, [1512.01981].
  * (64) V. Iršič et al., _New Constraints on the free-streaming of warm dark matter from intermediate and small scale Lyman- $\alpha$ forest data_, _Phys. Rev. D_ 96 (2017) 023522, [1702.01764].
  * (65) F. D’Eramo and A. Lenoci, _Lower mass bounds on FIMP dark matter produced via freeze-in_ , _JCAP_ 10 (2021) 045, [2012.01446].
  * (66) G. Ballesteros, M. A. G. Garcia and M. Pierre, _How warm are non-thermal relics? Lyman- $\alpha$ bounds on out-of-equilibrium dark matter_, _JCAP_ 03 (2021) 101, [2011.13458].
  * (67) Q. Decant, J. Heisig, D. C. Hooper and L. Lopez-Honorez, _Lyman- $\alpha$ constraints on freeze-in and superWIMPs_, _JCAP_ 03 (2022) 041, [2111.09321].
  * (68) G. G. Raffelt, _Stars as laboratories for fundamental physics: The astrophysics of neutrinos, axions, and other weakly interacting particles_.  5, 1996.
  * (69) M. M. Miller Bertolami, B. E. Melendez, L. G. Althaus and J. Isern, _Revisiting the axion bounds from the Galactic white dwarf luminosity function_ , _JCAP_ 10 (2014) 069, [1406.7712].
  * (70) G. Grilli di Cortona, E. Hardy, J. Pardo Vega and G. Villadoro, _The QCD axion, precisely_ , _JHEP_ 01 (2016) 034, [1511.02867].
  * (71) M. Badziak and K. Harigaya, _Naturally astrophobic QCD axion_ , _JHEP_ 06 (2023) 014, [2301.09647].
  * (72) P. Carenza, T. Fischer, M. Giannotti, G. Guo, G. Martínez-Pinedo and A. Mirizzi, _Improved axion emissivity from a supernova via nucleon-nucleon bremsstrahlung_ , _JCAP_ 10 (2019) 016, [1906.11844].
  * (73) M. Buschmann, C. Dessert, J. W. Foster, A. J. Long and B. R. Safdi, _Upper Limit on the QCD Axion Mass from Isolated Neutron Star Cooling_ , _Phys. Rev. Lett._ 128 (2022) 091102, [2111.09892].
  * (74) NA62 collaboration, E. Cortina Gil et al., _Measurement of the very rare K ${}^{+}\to{\pi}^{+}\nu\overline{\nu}$ decay_, _JHEP_ 06 (2021) 093, [2103.15389].
  * (75) CLEO collaboration, B. I. Eisenstein et al., _Precision Measurement of B(D+ — $>$ mu+ nu) and the Pseudoscalar Decay Constant f(D+)_, _Phys. Rev. D_ 78 (2008) 052003, [0806.2112].
  * (76) BaBar collaboration, B. Aubert et al., _A search for the decay $B^{+}\to K^{+}\nu\bar{\nu}$_, _Phys. Rev. Lett._ 94 (2005) 101801, [hep-ex/0411061].
  * (77) Belle-II collaboration, I. Adachi et al., _Evidence for $B^{+}\to K^{+}\nu\bar{\nu}$ Decays_, 2311.14647.
  * (78) BaBar collaboration, J. P. Lees et al., _Search for $B\to K^{(*)}\nu\overline{\nu}$ and invisible quarkonium decays_, _Phys. Rev. D_ 87 (2013) 112005, [1303.7465].
  * (79) W. Altmannshofer, A. Crivellin, H. Haigh, G. Inguglia and J. Martin Camalich, _Light New Physics in $B\to K^{(*)}\nu\bar{\nu}$?_, 2311.14629.
  * (80) F. Ertas and F. Kahlhoefer, _On the interplay between astrophysical and laboratory probes of MeV-scale axion-like particles_ , _JHEP_ 07 (2020) 050, [2004.01193].
  * (81) H. Leutwyler, _On the 1/N expansion in chiral perturbation theory_ , _Nucl. Phys. B Proc. Suppl._ 64 (1998) 223–231, [hep-ph/9709408].
  * (82) H. Georgi, _A bound on m(eta) / m(eta-prime) for large n(c)_ , _Phys. Rev. D_ 49 (1994) 1666–1667, [hep-ph/9310337].
  * (83) J. M. Gerard and E. Kou, _eta-eta-prime masses and mixing: A Large N(c) reappraisal_ , _Phys. Lett. B_ 616 (2005) 85–92, [hep-ph/0411292].
  * (84) N. Beisert and B. Borasoy, _eta eta-prime mixing in U(3) chiral perturbation theory_ , _Eur. Phys. J. A_ 11 (2001) 329–339, [hep-ph/0107175].
  * (85) D. S. M. Alves and N. Weiner, _A viable QCD axion in the MeV mass range_ , _JHEP_ 07 (2018) 092, [1710.03764].
  * (86) D. Aloni, Y. Soreq and M. Williams, _Coupling QCD-Scale Axionlike Particles to Gluons_ , _Phys. Rev. Lett._ 123 (2019) 031803, [1811.03474].
  * (87) H.-C. Cheng, L. Li and E. Salvioni, _A theory of dark pions_ , _JHEP_ 01 (2022) 122, [2110.10691].