aainstitutetext: The Oskar Klein Centre, Department of Physics, Stockholm
University, 10691 Stockholm, Swedenbbinstitutetext: Kavli Institute for the
Physics and Mathematics of the Universe (WPI), Chiba 277-8583,
Japanccinstitutetext: International Center for Quantum-field Measurement
Systems for Studies of the Universe and Particles (QUP,WPI), High Energy
Accelerator Research Organization (KEK), Oho 1-1, Tsukuba, Ibaraki 305-0801,
Japanddinstitutetext: Theory Center, Institute of Particle and Nuclear Studies
(IPNS), High Energy Accelerator Research Organization (KEK), Tsukuba 305-0801,
Japaneeinstitutetext: Graduate University for Advanced Studies (SOKENDAI),
1-1 Oho, Tsukuba, Ibaraki 305-0801, Japan

# Diffuse Axion Background

Joshua Eby c,d,e,b and Volodymyr Takhistov<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Relativistic axions can be readily produced in a broad variety of transient
sources, such as axion star bosenova explosions or supernovae. We develop a
general framework describing the resulting persistent diffuse axion background
(D$a$B) due to accumulated axions from historic transient events. We derive
strong constraints on the D$a$B flux from light axions $m\lesssim
10^{-3}\,{\rm eV}$ emitted from sources with energies $\omega\gtrsim{\rm MeV}$
considering the non-observation of excess photons associated with axion-photon
coupling from experiments, including COMPTEL, NuSTAR, XMM-Newton, INTEGRAL,
EGRET and Fermi. Future searches in experiments such as SKA, JWST, XRISM, Vera
C. Rubin Observatory, AMEGO/e-ASTROGAM will allow probing D$a$B and associated
axion-photon couplings with unprecedented sensitivity covering a wide range of
possible source energies as low as $0.1\,\mu$eV and multiple decades in axion
masses. We highlight the differences between astrophysical and dark sector
sources of D$a$B. Further, we discuss complementarity with direct detection as
well as prospects for other D$a$B searches. Our analysis demonstrates that
D$a$B can act as a promising probe of populations of axion emission sources as
well as emission mechanisms.

††preprint: KEK-QUP-2024-0001, KEK-TH-2596, KEK-Cosmo-0337, IPMU24-0003

##  1 Introduction

The discovery of Higgs boson ATLAS:2012yve ; CMS:2012qbp  has highlighted the
central role of scalar spin-0 particles in Nature. Existence of pseudo-scalar
spin-0 quantum chromodynamics (QCD) axions would readily address the long-
standing strong CP problem of the Standard Model (SM) Peccei:1977hh ;
Weinberg:1977ma ; Wilczek:1977pj ; Kim:1979if ; Shifman:1979if ;
Zhitnitsky:1980tq ; Dine:1981rt  and provide a compelling candidate for the
dark matter (DM) Preskill:1982cy ; Abbott:1982af ; Dine:1982ah . Axions and
more generally axion-like particles111Here, we denote axion and axion-like
particles as _axions_., pseudo-scalars that generically arise from broken
shift symmetries, are expected to be ubiquitous from fundamental theory
Svrcek:2006yi ; Arvanitaki:2009fg . Significant efforts have been devoted to
explore the axion landscape (see e.g. Irastorza:2018dyq ; Agrawal:2021dbo ;
Adams:2022pbo  for an overview), primarily focusing on non-relativistic cold
axion DM.

Relativistic axions can appear from a variety of sources. Cosmological axion
production at redshifts prior to galaxy and star formation ($z\gtrsim 30$)
will result in new contributions to the cosmic axion background Conlon:2013isa
; Marsh:2013opc ; Dror:2021nyr . Another prominent possibility are thermal
axions, which are produced in the early Universe if axions were in thermal
equilibrium through interactions with SM constituents Turner:1986tb ;
Chang:1993gm ; Masso:2002np ; Hannestad:2005df ; Graf:2010tv ; Salvio:2013iaa
; Arias-Aragon:2020shv ; Dror:2021nyr . Efficient non-thermal cosmological
axion production is also possible, such as due to particle decays Chun:1995hc
; Asaka:1998xa ; Ichikawa:2007jv ; Kawasaki:2007mk  or non-linear dynamics of
field oscillations Kasuya:1996ns ; Ema:2017krp ; Co:2017mop . On the other
hand, continuous emission at low-redshifts of relativistic axions from
astrophysical sources like stars and the Sun has resulted in variety of
sensitive tests of possible axion parameter space, such as by considerations
of star cooling Raffelt:1996wa  or direct detection searches for solar axion
emission Sikivie:1983ip ; Moriyama:1995bz .

Recently, a fruitful new direction focusing on manifestations of relativistic
burst axion emissions from transient sources has emerged and is being
explored, including axion (boson) star explosions Eby:2016cnq ; Levkov:2016rkk
, black hole superradiance Arvanitaki:2010sy ; Arvanitaki:2014wva ;
Arvanitaki:2016qwi ; Baryakhtar:2020gao ; Unal:2020jiy ; Branco:2023frw , and
solar halos Banerjee:2019epw ; Banerjee:2019xuy ; Budker:2023sex  that can
repeatedly explode on astrophysical timescales. Direct detection of
relativistic bosons emitted from transient sources was found to result in a
rich variety of signatures in the context of axion searches Eby:2021ece ;
Chigusa:2023bga  as well as precision quantum sensors Arakawa:2023gyq , often
allowing to probe novel parameter space that is beyond the reach of
conventional cold DM searches. Furthermore, such axion manifestations hold
promising avenues for exploration of new physics in the context of
multimessenger astronomy Dailey:2020sxa ; Eby:2021ece ; Arakawa:2023gyq .

Accumulation of relativistic axions from historic transient sources will
result in persistent axion flux filling the Universe Raffelt:2011ft ;
Calore:2020tjw ; Calore:2021hhn ; Takhistov:2021qzv . We comprehensively
explore this diffuse axion background (D$a$B). An analogy can be drawn between
D$a$B and neutrinos emitted from historic astrophysical sources, such as
supernovae, which will contribute to a diffuse supernova neutrino background
(see e.g. Beacom:2010kk  for review), detection of which is a target of
ongoing efforts by Super-Kamiokande Horiuchi:2008jz ; Super-Kamiokande:2021jaq
and other neutrino experiments, or other possible sources such as supermassive
star explosions Shi:1998nd ; Linke:2001mq ; Munoz:2021sad .

In this work, we develop a general framework and systematically analyze D$a$B
originating from distinct sources. Depending on the origin, D$a$B
contributions can dramatically differ. We discuss variety of general
approaches for D$a$B detection, including with existing as well as future
experiments.

The format of this paper is as follows. In Sec. 2 we introduce the general
framework for D$a$B and discuss different possible contributing sources. In
Sec. 3 we focus on photon production from D$a$B. Then, in Sec. 4 we analyze
observable signatures with photons originating from D$a$B for different
existing as well as future experiments. Later, in Sec. 5 we focus on analyzing
D$a$B associated with axion star bosenovae emission as a characteristic source
example. We discuss other D$a$B search prospects in Sec. 6, including
different axion couplings as well as direct detection. We conclude in Sec. 7.

##  2 Diffuse axion background flux

The present-day accumulated D$a$B number flux spectrum that originates from
historic transient sources can be stated in generality in analogy with the
well-studied diffuse background flux of SM weakly interacting neutrinos
originating from sources like historic core-collapse supernovae Beacom:2010kk
; Ando:2023fcc  (or e.g. supermassive star explosions Shi:1998nd ; Shi:1998jx
; Munoz:2021sad ) as222The conventional units for emission spectrum $\phi$ are
cm-2s-1.

$\frac{d\phi}{d\omega}=\int_{0}^{\infty}dz(1+z)\frac{dN_{a}(\omega(1+z))}{d\omega}R_{\rm
burst}(z)\left|\frac{dt}{dz}\right|~{},$ (1)

where $\omega$ is the axion energy, $R_{\rm burst}(z)$ is the rate density of
transient “burst” events emitting axions due to a particular source (in units
of events per time and per volume) as a function of cosmological redshift $z$,
and $dN_{a}/d\omega$ is the differential spectrum of the number of axions
$N_{a}$ per unit energy $\omega$ characterizing emitted axions from a single
transient burst with the $(1+z)$ factor accounting for the energy redshifting
of axions emitted at $z>0$. The last term encodes the cosmology of flat
expanding Universe dominated by matter and cosmological constant
$\left|dt/dz\right|^{-1}=H_{0}(1+z)\sqrt{\Omega_{\Lambda}+\Omega_{M}(1+z)^{3}}$,
with $H_{0}=67.4$ km s-1 Mpc-1, and densities $\Omega_{\Lambda}=0.685$ for
cosmological constant and $\Omega_{M}=0.315$ for matter compatible with latest
measurements from Planck Planck:2018vyg . The flux for various possible
sources, including both astrophysical as well as dark-sector sources, is
illustrated in Fig. 1.

Figure 1: Sources of relativistic axion flux emission, which can contribute to
the D$a$B as described in Sec. 2.2. Axion star bosenova emission estimated
from numerical simulations of Ref. Levkov:2016rkk  (orange shaded region
corresponds to axion decay constant between $f_{a}=10^{14}$ GeV (top) and
$f_{a}=10^{10}$ GeV (bottom)), supernova (SN) emission estimated from Ref.
Calore:2020tjw  (blue dashed corresponding to brehsstrahlung with axion proton
coupling $g_{ap}=10^{-9}$, $g_{a\gamma}=0$ allowed by constraints from
excessive SN cooling Carenza:2019pxu , axion-photon coupling, blue solid for
Primakoff processes with $g_{a\gamma}=10^{-12}$ GeV-1, and $\pi N\to Na$
production using $g_{ap}=10^{-9}$ Carenza:2020cis ; Lella:2022uwi ), and
neutron star (NS) merger flux estimated from Ref. Harris:2020qim ;
Fiorillo:2022piv  (with the energy fixed by the temperature as
$\omega/3=T=(20-100)\,{\rm MeV}$, density between $n=(1-7)n_{0}$, and neutron
coupling fixed to $g_{an}=1.5\cdot 10^{-9}$). Expected energy range for NS
polar cap emission, astrophysical thermal emission, and dark star thermal
emission are denoted by arrows (see discussion of Sec. 2.2.3). Reference flux
from the local DM density inside an experimental volume of ${\rm m}^{3}$ is
given by the gray dotted line.

Looking ahead, relevant for D$a$B detection are axion interactions with SM,
which generally depend on specific theory. More so, for astrophysical
transient sources, axion-SM interactions are also generally responsible for
their production as energy is converted from SM constituents to axions. A
canonical coupling of pseudoscalar axion $a$ is the two-photon coupling
described by the Lagrangian contribution

${\mathcal{L}}\supset\frac{1}{4}g_{a\gamma}aF_{\mu\nu}\tilde{F}^{\mu\nu}=g_{a\gamma}a\textbf{E$\cdot$B}~{},$
(2)

where $F_{\mu\nu}$ is the electromagnetic field-strength tensor,
$\tilde{F}^{\mu\nu}$ is its dual, and $g_{a\gamma}$ the axion-photon
coupling333This coupling has typical dimension of GeV-1.. Interacting with
external magnetic field B, the coupling of Eq. (2) results in axion-photon
conversion that is a prominent source of phenomenological signatures.
Additional possible axion interactions include

${\mathcal{L}}\supset\frac{g_{a\chi}}{2m_{\chi}}\partial_{\mu}a\overline{\chi}\gamma^{\mu}\gamma_{5}\chi+\frac{g_{ag}}{4}aG_{\mu\nu}\tilde{G}^{\mu\nu}\,,$
(3)

where $G_{\mu\nu}$ is the gluon field strength tensor, $\tilde{G}_{\mu\nu}$ is
its dual, $\chi$ denotes fermions such as nucleons $N$ (protons, neutrons) or
electrons $e$, $m_{\chi}$ is their mass, and $g_{a\chi}$ and $g_{ag}$ are
coupling constants. For QCD axion that resolves the CP problem in the SM QCD
sector, the axion mass $m$ and couplings are related to the Peccei-Quinn
symmetry breaking scale, also called the axion decay constant, $f_{a}$. In
particular, $(m/1~{}{\rm eV})=(7\times 10^{6}\,{\rm GeV}/f_{a})=0.5\xi
g_{a\gamma}/(10^{-10}~{}{\rm GeV}^{-1})$ with $\xi$ being a model-dependent
constant, which in KSVZ QCD axion scenario Kim:1979if ; Shifman:1979if  is
$|\xi|\simeq 0.5$ and in DFSZ QCD axion scenario Dine:1981rt ;
Zhitnitsky:1980tq  is $|\xi|\simeq 1.4$.

###  2.1 Cosmological transient event rate

To capture a broad class of possible D$a$B sources, it is convenient to
parameterize the cosmological transient event rate as

$R_{\rm burst}(z)=\frac{A\rho_{\rm loss}H_{0}}{E_{\rm tot}(z)}f(z)~{},$ (4)

where $f(z)$ is a dimensionless function of redshift $z$ characterizing
transient source distribution. Here, the normalization constant $A$ is fixed
by the requirement that the energy density $\rho_{\rm loss}$ in axions emitted
from a particular transient source type, such as exploding axion stars,
constitutes a specified subdominant fraction of the total DM density444In case
of exploding supermassive stars, analogous normalization can be imposed by
considering relation with baryon abundance Munoz:2021sad .

$\rho_{\rm loss}=\int dtE_{\rm tot}(z)\frac{R_{\rm
burst}(z)}{(1+z)^{3}}=\frac{1}{H_{0}}\int dz\frac{E_{\rm tot}(z)R_{\rm
burst}(z)}{(1+z)^{4}\sqrt{\Omega_{\Lambda}+\Omega_{M}(1+z)^{3}}}~{},$ (5)

where $E_{\rm tot}(z)$ is the total energy emitted in each transient event
burst at redshift $z$. Since for the sources we consider the total emitted
energy for each event is predominantly independent of redshift, we take
$E_{\rm tot}$ to be constant in what follows. We can set the normalization
constant $A$ by inverting Eq. (5) as

$A=\left[\int
dz^{\prime}\frac{f(z^{\prime})}{(1+z)^{4}\sqrt{\Omega_{\Lambda}+\Omega_{M}(1+z)^{3}}}\right]^{-1}\,,$
(6)

which can be directly computed for a given source model $f(z)$. In general,
$f(z)$ could have non-vanishing contributions over the redshift range
$z\in\\{z_{\rm min},z_{\rm max}\\}$, where $z_{\rm max}$ is the redshift when
explosive transient axion emission events begin to occur and $z_{\rm min}$
when the transient events subside. We consider $z_{\rm min}\rightarrow 0$,
although theories of transient axion sources with non-negligible $z_{\rm min}$
can be readily incorporated in our analysis.

The energy density $\rho_{\rm loss}$ can be intuitively expressed as the
fraction ${\mathcal{F}}$ of the average DM energy density in the Universe,
$\rho_{\rm DM}\simeq 1.2\times 10^{-6}\,{\rm GeV/cm}^{3}$, as $\rho_{\rm
loss}=\mathcal{F}\rho_{\rm DM}$. We consider ${\mathcal{F}}=0.1$ as a
reference benchmark, which can evade various constraints. As we will show
below, some of our signatures are detectable even for significantly smaller DM
fractions.

The cosmological event rate $R_{\rm burst}(z)$ characterizes the epochs during
which the transient sources are most active. Employing the normalization of
Eq. (4), we analyze two characteristic phenomenological models of source
events, described by

$f(z)=\begin{cases}(1+z)^{p}\Theta\left(z-z_{\rm max}\right)&{\rm
power\,\,law}\\\ \exp\left[-(z-\bar{z})^{2}/\delta z^{2}\right]&{\rm
Gauss.}\end{cases}\,\,,$ (7)

where $p$ is an integer, and $\bar{z}$ with $\delta z$ denote the Gaussian
distribution mean and width, respectively. The power law is motivated by
scenarios where the transient source rate is smoothly varying with possible
additional redshift dependence as specified by parameter $p$, reaching a
maximum in the present-day for $p<0$ or at earlier times for $p>0$. The
Gaussian model instead characterizes general scenarios where transient events
predominantly occur during a relatively short cosmological epoch, which could
be associated with some dynamical processes, with a duration approximated by
$z\in\\{\bar{z}-\delta z,\bar{z}+\delta z\\}$. We illustrate several
distributions considered in this work in Fig. 2. Our analysis can readily be
carried over to other choices of source event rates, specific to particular
theories.

Figure 2: Dimensionless cosmological transient axion burst event rate
distribution $f(z)$ sourcing relativistic axions that are considered in this
work. The power-law and Gauss distributions are defined in Eq. (7), and we set
the input parameters as labeled. For illustration we normalize $f(z)$ to
unity.

With $f(z)$ given by Eq. (7), we can directly compute normalization $A$ of Eq.
(6). For the power-law model we find

$A_{1}=\left[\int_{0}^{z_{\rm
max}}dz^{\prime}\frac{(1+z^{\prime})^{p-4}\Theta(z-z_{\rm
max})}{\sqrt{\Omega_{\Lambda}+\Omega_{M}(1+z^{\prime})^{3}}}\right]^{-1}\simeq\frac{1}{\sqrt{\Omega_{M}}}\left(\frac{9}{2}-p\right)+{\mathcal{O}}(z_{\rm
max}^{p-9/2})~{},$ (8)

which is valid as long as $p<9/2$, and where we approximated
$\sqrt{\Omega_{\Lambda}+\Omega_{M}(1+z)^{3}}\simeq\sqrt{\Omega_{M}(1+z)^{3}}$
considering Universe is matter-dominated for the majority of redshift values
of interest555Corrections induced by this approximation are of order
${\mathcal{O}}(1)$ at most.. Further corrections are suppressed by powers of
$1/z_{\rm max}$, which we consider to be subdominant compared to unity. For
the Gaussian distribution

$A_{2}=\left[\int_{0}^{z_{\rm
max}}dz^{\prime}\frac{\exp\left(-\left(\frac{z^{\prime}-\bar{z}}{\delta
z}\right)^{2}\right)}{(1+z^{\prime})^{4}\sqrt{\Omega_{\Lambda}+\Omega_{M}(1+z^{\prime})^{3}}}\right]^{-1}\,.$
(9)

Evaluating Eq. (9) numerically, one finds $A_{2}\simeq 1-10^{6}$ upon varying
$\delta z\simeq 10^{-2}-1$ and $\bar{z}\simeq 0-20$. In this estimate we used
$z_{\rm max}=20$, but the results do not depend significantly on this choice
as long as $z_{\rm max}>\bar{z}$ is satisfied. In what follows, we focus on
the power-law model for simplicity and outline results for the Gaussian model
in App. A.

###  2.2 Emission sources

The flux of axions in the relativistic background, given in Eq. (1), depends
sensitively on the spectrum of emission from transient sources
$dN_{a}/d\omega$. Our general analysis framework can be applied to a wide
range of possible transient sources. A generic description of a peaked axion
burst can be captured by a Gaussian function

$\frac{dN_{a}(\omega)}{d\omega}\Bigg{|}_{\rm Gauss}=C\frac{E_{\rm
tot}}{m^{2}}\frac{\exp\left(-\left(\frac{\omega-\bar{\omega}}{\delta\omega}\right)^{2}\right)}{\sqrt{\pi}\,(\delta\omega/m)}\,,$
(10)

with a central peak energy $\bar{\omega}$ and a width $\delta\omega$. The
normalization of the emission spectrum $m\int d\omega(dN_{a}/d\omega)=E_{\rm
tot}$ implies that the normalization constant is given by

$\displaystyle
C=m\left[\int\frac{\exp\left(-\left(\frac{\omega-\bar{\omega}}{\delta\omega}\right)^{2}\right)}{\sqrt{\pi}\,(\delta\omega/m)}d\omega\right]^{-1}$
$\displaystyle=2\left[1+{\rm
erf}\left(\frac{\bar{\omega}}{\delta\omega}\right)\right]^{-1}\,.$ (11)

Note that $C$ is independent of $m$ and $E_{\rm tot}$.

We can employ the Gaussian function description to determine the axion
background density for a generalized emission source, considering the emission
follows an approximately peaked structure. In the limit $\delta\omega\to 0$,
$(dN_{a}/d\omega)_{\rm Gauss}\to(CE_{\rm tot}/m)\delta(\omega-\bar{\omega})$,
capturing the effect of a narrow burst spectrum in the $\delta\omega\ll m$
limit. Combining Eqs. (1), (4), (10), and (11), we obtain

$\frac{d\phi}{d\omega}\Big{|}_{\rm Gauss}=\frac{{\mathcal{F}}\rho_{\rm
DM}}{m\delta\omega}\frac{2}{\sqrt{\pi}}\frac{A}{\left[1+{\rm
erf}(\bar{\omega}/\delta\omega)\right]}\int\frac{dzf(z)}{\sqrt{\Omega_{\Lambda}+\Omega_{M}(1+z)^{3}}}\exp\left[-\left(\frac{(\omega(1+z)-\bar{\omega})}{\delta\omega}\right)^{2}\right]\,.$
(12)

We stress that the flux model of Eq. (12) captures key behavior while being
particularly straightforward to analyze. For a given particle mass $m$ and
cosmological rate function $f(z)$, one needs only to specify the central burst
energy $\bar{\omega}$, the emission width $\delta\omega$, and the total DM
fraction converted to relativistic particles ${\mathcal{F}}$.

We systematically discuss below several prominent examples of transient axion
burst sources resulting in D$a$B highlighting the generality and reach of our
analysis. In Fig. 1 we summarize axion emission flux from several
characteristic sources, as compared to reference local DM flux within a volume
of m3 (see Eq. (22)). Analogously to astrophysical multimessenger targets,
astrophysical sources such as supernovae can serve as prime production sites
of relativistic axions666We note that some astrophysical axion emission
mechanisms can persist over significant timescales. Here, we comment on them
as characteristic and reference sources in the context of potentially
contributing to D$a$B.. This exemplifies conversion of energy associated with
SM constituents (e.g. baryons) into dark sectors with large kinetic energy. In
this work, in addition to astrophysical sources, we emphasize the significant
potential of dark-sector dynamics to produce relativistic axions resulting in
sizable contributions to D$a$B, which we can describe using a general and
flexible framework. We further illustrate how this can give rise to photon
signatures across a broad variety of source classes, range of energy scales,
and distinct photon-production scenarios.

####  2.2.1 Supernovae

Astrophysical sources could serve as powerful production sites not only of
neutrinos, but also of other weakly interacting particles such as axions.
Core-collapse supernovae (SN) allow to sensitively probe axions through
variety of methods (e.g. Raffelt:1990yz ; Raffelt:1996wa ; Raffelt:2006cw ;
Fischer:2016cyd ). Aside first direct SN neutrino observations Kamiokande-
II:1987idp ; Bionta:1987qt ; Alekseev:1988gp , SN1987A that exploded in the
Large Magellanic Cloud at a distance of 50 kpc proved to be an essential
laboratory for axion physics (e.g. Turner:1987by ; Burrows:1988ah ;
Burrows:1990pk ; Raffelt:1987yt ; Raffelt:1990yz ; Keil:1996ju ;
Hanhart:2000ae ; Chang:2018rso ). In SN, efficient axion production can occur
through axion nucleon-nucleon bremsstrahlung process Turner:1987by
$NN\rightarrow NNa$, where $N$ denotes nucleons (protons, neutrons). Refs.
Carenza:2020cis ; Lella:2022uwi  have illustrated the importance of a pion
production channel $\pi N\rightarrow Na$ as well, which can dominate axion
emission above energies of order $\sim 150\,{\rm MeV}$. Requiring that SN1987A
neutrino emission is not significantly reduced by axion emission, we consider
axion-nucleon coupling777Exact value depends on specific considered axion
model. $g_{aN}\lesssim 10^{-9}$ Carenza:2019pxu . When axions couple only to
photons with coupling $g_{a\gamma}$, efficient production can occur through
Primakoff process when thermal photons are converted into axions within
electromagnetic field of stellar medium Raffelt:1985nk . This results in
strong limits on axions from SN energy loss arguments more stringent than
other limits above tens of keV Jaeckel:2010ni ; Lella:2023bfb . For light SN
axions constraints arise from conversion to photons in Milky Way Galactic
magnetic field Grifols:1996id ; Brockway:1996yr ; Payez:2014xsa . For heavy SN
axions constraints arise from decays along the emitted axion line of sigh and
non-observation of SN1987A coincidence signal Jaeckel:2017tud .

Contributions to D$a$B from SN Raffelt:2011ft ; Calore:2020tjw ;
Calore:2021hhn  can be effectively parameterized using a fit to the axion
spectrum arising due to a given process channel $x$ of the form Calore:2020tjw

$\frac{dN_{a}(\omega)}{d\omega}\Bigg{|}_{\rm
supernova}=c\left(\frac{g_{ax}}{g_{ax}^{\rm
ref}}\right)^{2}\left(\frac{\omega}{\omega_{0}}\right)^{\beta}\exp\left(-\frac{(\beta+1)\omega}{\omega_{0}}\right)\,,$
(13)

where $c$, $g_{ax}^{\rm ref}$, $\omega_{0}$, and $\beta$ are fit parameters
which depend on the production process (see e.g. Tab. 1 of Ref. Calore:2020tjw
).

Importantly, for astrophysical transients such as SN contributing to D$a$B the
event rate $R_{\rm burst}$ is determined by conventional astrophysical
processes (e.g. star formation, stellar nuclear fusion) such as in case of
diffuse neutrino background. However, this is distinct for D$a$B contributions
arising from dark sector transient sources such as axion stars explosions.

####  2.2.2 Axion star bosenovae

Axion stars are gravitationally-bound states formed from ultralight bosons
Kaup:1968zz ; Ruffini:1969qy ; BREIT1984329 ; Colpi:1986ye . They can appear
in the early Universe within the in cores of diffuse axion miniclusters
Kolb:1993zz ; Eggemeier:2019jsu . When their mass reaches a critical value
$M_{c}$ Chavanis:2011zi ; Chavanis:2011zm ; Barranco:2010ib ; Eby:2014fya ;
Schiappacasse:2017ham , defined by the strength of the axion self-interaction
coupling, these states undergo gravitational collapse and finally emit a large
fraction of their mass-energy in relativistic bosons Eby:2016cnq ;
Levkov:2016rkk . This explosive production of relativistic particles
signifies888The term was originally coined to describe an explosive process in
atomic physics Donley_2001 . _bosenova_.

Detailed numerical simulations of boson star collapse and explosive axion
emission were carried out in Ref. Levkov:2016rkk . The spectrum of emitted
axions $dE/dk$ has been analyzed up to $k/m\simeq 9$, where $E(\omega)=\omega
N_{a}(\omega)$ is the total energy emitted in axions of energy $\omega$, and
$k=\sqrt{\omega^{2}-m^{2}}$ is the corresponding momentum. The total energy
emitted is a linear function of the axion decay constant $f_{a}$, well-fitted
by $E_{\rm tot}\simeq M_{c}[0.3+10^{3}f_{a}/M_{\rm Pl}]$, which saturates to
unity near $f_{a}\simeq 10^{16}\,{\rm GeV}$, and where $M_{c}=10M_{\rm
Pl}f_{a}/m$. To analyze resulting D$a$B flux using Eq. (1), we translate

$\frac{dN_{a}(\omega)}{d\omega}\Bigg{|}_{\rm
bosenova}=\left[\frac{dE}{dk}-\frac{d\omega}{dk}\right]\frac{1}{\omega\displaystyle{\frac{d\omega}{dk}}}\simeq\frac{dE}{dk}\frac{1}{\sqrt{\omega^{2}-m^{2}}}~{},$
(14)

where in the last step we used $d\omega/dk=k/\omega$, and neglected the second
term in the bracket, which is suppressed by a factor $1/N_{a}\ll 1$ relative
to the leading term. In this case, Eqs. (1) and (4) give

$\displaystyle\frac{d\phi}{d\omega}\Bigg{|}_{\rm bosenova}$
$\displaystyle\simeq\frac{{\mathcal{F}}\rho_{\rm
DM}A}{\sqrt{\omega^{2}-m^{2}}}\int_{0}^{\infty}dz\frac{f(z)}{\sqrt{\Omega_{\Lambda}+\Omega_{M}(1+z)^{3}}}\frac{1}{E_{\rm
tot}}\frac{dE(\omega(1+z))}{dk}~{}.$ (15)

Figure 3: Spectrum of emitted axions as a function of energy, from numerical
simulations of axion bosenova explosion Levkov:2016rkk  (blue solid) with a
Gaussian fit (black dashed) using Eq. (10) to the leading relativistic peak at
$\bar{\omega}\simeq 2.2m$ and $\delta\omega\simeq 0.4m$.

Key properties of these axion star explosions can be obtained from the
emission spectrum Levkov:2016rkk , which is peaked at semi-relativistic
energies (i.e. momenta of order $\sim$few$\times m$). The peaks can be
understood as a consequence of number-changing processes (see e.g. Ref.
Eby:2015hyx ; Eby:2017azn ). In Fig. 3 we illustrate the full simulated
spectrum of Ref. Levkov:2016rkk , along with a Gaussian fit of the first
relativistic peak, which is well-described by Eq. (10) with
$\bar{\omega}\simeq 2.2m$ and $\delta\omega\simeq 0.4$ as shown. Explosive
bosenova emissions across astrophysical timescales constitute a prominent
contribution to D$a$B arising from dark-sector sources.

Recently, our analyses of direct detection signatures from bosenova explosions
occurring inside the Milky Way demonstrated that they constitute particularly
promising targets for current and near-future experimental searches
Eby:2021ece ; Arakawa:2023gyq . Ref. Du:2023jxh ; Escudero:2023vgv  suggested
that repeated mergers of DM subhalos can give rise to a large abundance of
near-critical boson stars, obtaining a DM fraction as large as $\sim 10^{-6}$.
These stars would collapse at redshift $z\lesssim 30$, with the resulting
relativistic axions contributing to D$a$B. Bosenovae occurring at earlier
epochs of redshifts $z\gtrsim 30$ can also contribute, albeit with
significantly redshifted axion energies, and can be constrained by
cosmological observations Fox:2023aat .

####  2.2.3 Neutron stars and other sources

Neutron stars (NSs) constitute another possible prominent astrophysical source
of relativistic axion emission. Axions can be thermally produced inside NSs
through various channels, such as nucleon-nucleon bremsstrahlung
$NN\rightarrow NNa$ and Cooper pair-breaking processes, cooling NSs999Cooling
of white dwarfs and red giants through axion emission can also result in
sensitive constraints Altherr:1993zd ; Raffelt:1994ry ; Corsico:2001be ;
MillerBertolami:2014rka . (e.g. Iwamoto:1984ir ; Keller:2012yr ;
Sedrakian:2015krq ). Such emission can be enhanced during NS mergers, a
process which has been studied using numerical Dietrich:2019shr  and analytic
Harris:2020qim ; Fiorillo:2022piv  approaches. The typical energy scale of
emission is near the temperature of a NS $T\simeq{\mathcal{O}}(1-100)\,{\rm
MeV}$, provided the axions are massless by comparison, i.e. $m\ll T$.
Furthermore, for typical nuclear couplings allowed by current constraints (see
e.g. Adams:2022pbo  for review), the axions easily escape the star and free-
stream away.

The axion emissivity in NS mergers as a function of temperature, $Q_{0}(T)$,
has been estimated in relativistic mean field theory Harris:2020qim  for an
axion-neutron coupling of $g_{0}=1.5\times 10^{-9}$. For other couplings, the
emissivity scales as

$Q(T)=\left(\dfrac{g_{an}}{g_{0}}\right)^{2}Q_{0}(T)\,.$ (16)

Following Ref. Harris:2020qim , we use the temperature $T$ as a proxy for the
typical axion emission energy as $\omega\simeq 3T$. Taking an approximate
merger duration of $t_{\rm merger}\sim 1\,{\rm s}$ and emitting volume of
$V\sim(15\,{\rm km})^{3}$ Fiorillo:2022piv , we estimate

$\frac{dN_{a}}{d\omega}\Big{|}_{\rm NS\,merger}\simeq\frac{V\,t_{\rm
merger}}{\omega^{2}}Q(T)\,.$ (17)

Here, the density of the Fermi fluid $n$ is varied between the nuclear
saturation density, $n_{0}\simeq 10^{-44}\,{\rm m}^{3}$, and $7n_{0}$. This
range determines the extent of the gray shaded region in Fig. 1. Some of the
emitted axions convert to gamma rays in the magnetosphere, leading to
detectable, transient photon signals Fiorillo:2022piv . Here, we consider the
possibility of such axions accumulating over time and contributing to D$a$B.

Recently, it was pointed out that axions can be produced directly via a
coupling to photons in the polar-cap region of NS Prabhu:2021zve ;
Noordhuis:2022ljw . While some become trapped around the NS, for axions with
masses $m\lesssim 10^{-7}\,{\rm eV}$ the vast majority escape the NS (as many
as $\sim 10^{50}$ axions per second), as the energy scale of emission is fixed
by dynamics in the polar cap in the range $10^{-7}\,{\rm
eV}\lesssim\omega\lesssim 10^{-4}\,{\rm eV}$ Noordhuis:2023wid . As we discuss
further below, such frequency ranges could be probed with future Square
Kilometer Array (SKA) observations.

A broad variety of other sources of relativistic axion emission that could
contribute to D$a$B are possible, such as:

  * •

Black hole superradiance: “clouds” of axions can be generated by sapping the
angular momentum and energy from rapidly-rotating black holes Press:1972zz ;
Arvanitaki:2010sy ; Arvanitaki:2014wva ; Arvanitaki:2016qwi ;
Baryakhtar:2020gao ; Unal:2020jiy  (see e.g. Brito:2015oca  for review).
Superradiance effects around early-Universe (primordial) black holes have also
been considered Ferraz:2020zgi ; Bernal:2022oha ; Branco:2023frw . Massive
axion fields that have quasi-bound states around black holes with
$\omega<l\Omega_{H}$, where $\Omega_{H}$ is the angular velocity of the black
hole horizon101010Here,
$\Omega_{H}=\frac{1}{r_{S}}\Big{(}\frac{a_{*}}{1+\sqrt{1-a_{*}^{2}}}\Big{)}$,
where $r_{S}=2GM_{\rm BH}$ is the black hole Schwarzschild radius and
$a_{*}=J/GM_{\rm BH}^{2}$ is the dimensionless spin of black hole with angular
momentum $J$.  and $l$ is the angular momentum around the black hole spin
axis, can experience exponentially unstable growth. The timescale for spin-
down of the BH is give by Baryakhtar:2020gao

$\tau_{\rm spin-down}\simeq 10^{7}\,{\rm
years}\left(\frac{0.02}{r_{S}m}\right)^{5}\left(\frac{10^{-12}\,{\rm
eV}}{m}\right)\left(\frac{0.9}{a_{*}}\right)^{3/2}\left(\frac{10^{15}\,{\rm
GeV}}{f_{a}}\right)^{2}\,.$ (18)

The growth continues until the black hole loses most of its angular momentum
or the superradiance is quenched by the axion self-interactions, either by
exciting axions out of superradiant modes leading to stable steady-state
configuration or by inducing gravitational collapse. In the latter case, an
explosive bosenova-like event occurs near the black hole with rapid emission
of relativistic axions Arvanitaki:2010sy . Axion emission can be modified in
case of binary black hole mergers Baumann:2018vus . We note that the event
rates of bosenovae-type events associated with black hole supperradiance
sensitively depends on a variety of factors such as black hole population
statistics, distributions of black hole masses and spins, as well as particle
physics effects like interaction couplings of the axions.

  * •

Stellar axion halos: axions can be captured directly in the gravitational
field of astrophysical bodies, such as stars. These configurations are akin to
gravitational atoms, having a spherically-symmetric density profile of the
form Banerjee:2019epw ; Banerjee:2019xuy

$\rho(r,t)\simeq\rho_{0}(t)\exp\left(-\frac{r^{2}}{R_{\star}^{2}}\right)\,,$
(19)

with radius $R_{\star}$ and central density $\rho_{0}(t)$ which can grow
exponentially in magnitude with time Budker:2023sex . Gravitational collapse
of these structures around stars across galaxies can potentially lead to
occurrence of yet another class of bosenova-like explosive events. The rate of
growth happens on a timescales similar to the known relaxation timescales for
ultralight particles111111See also Ref. Levkov:2018kau ; Chen:2020cef ;
Kirkpatrick:2020fwd ; Chen:2021oot ; Kirkpatrick:2021wwz ; Dmitriev:2023ipv ;
Jain:2023tsr  for discussion of relaxation of ultralight fields through
gravity and self-interactions. Budker:2023sex , which for $2\to 2$ self-
interactions takes the form

$\tau_{\rm rel}\simeq\frac{64f_{a}^{4}m^{3}v_{\rm lDM}^{2}}{\rho_{\rm
lDM}^{2}}\simeq 9\,{\rm Gyr}\left(\frac{f_{a}}{10^{8}\,{\rm
GeV}}\right)^{4}\left(\frac{m}{10^{-14}\,{\rm
eV}}\right)^{3}\left(\frac{0.4\,{\rm GeV/cm}^{3}}{\rho_{\rm
lDM}}\right)\left(\frac{v_{\rm lDM}}{10^{-3}}\right)\,,$ (20)

where $\rho_{\rm lDM}\simeq 0.4\,{\rm GeV/cm}^{3}$ and $v_{\rm lDM}\simeq
10^{-3}$ are the energy density and velocity dispersion, respectively, of
local DM in the solar neighborhood vicinity (see e.g. Gelmini:2015zpa  for
review). Thus, on astrophysical timescales for sufficiently-strong self
interactions, solar halos can collapse and explode multiple times for a given
star, potentially converting a significant fraction of axionic DM to
relativistic D$a$B constituents in the process.

  * •

Dark non-bosonic stars: Astrophysical objects formed from SM particles (e.g.
stars and NS) could have analogues in the dark sector, sometimes called _dark
stars_ Maselli:2017vfi ; Kamenetskaia:2022lbf , composed of fermions rather
than bosons. If their constituents are coupled to axions, they will be emitted
by dynamical processes over a broad range of energies dictated by the dark
star core temperature $T_{\rm dark}$, which can be drastically distinct from
$\sim{\rm MeV}$ energies associated with ordinary astrophysical stars. An
intriguing possibility is _mirror star_ , formed from dark baryons in the
context of complex dark sector models Curtin:2019lhm ; Hippert:2021fch ;
Hippert:2022huw . The constituents of dark stars and the interacting SM
particles could contribute to axion emission, and therefore to the D$a$B.
Population statistics and dynamics of such objects could be dramatically
distinct from that of conventional astrophysical stars.

Many relativistic axion emission mechanisms remain poorly understood without a
well studied axion emission flux. Our analysis highlights the necessity and
serves as an impetus for further detailed numerical studies and simulations
required to comprehensively describe and quantify the axion emission
mechanisms from distinct sources. A population of astrophysical sources that
produces energetic axions can contribute to D$a$B.

###  2.3 Comparison to local DM flux

A simple estimate of D$a$B flux can be obtained as follows. Assuming a total
density $\rho_{\rm loss}$ is converted into relativistic axions of average
energy $\bar{\omega}$, ignoring redshift effects, the differential flux can be
approximated as

$\frac{d\phi}{d\omega}\sim\frac{\rho_{\rm loss}}{\bar{\omega}^{2}}\,v_{\rm
burst}~{}.$ (21)

This scaling can be roughly obtained from Eq. (15) in the approximation
$dE/dk\to E_{\rm tot}/\bar{\omega}$,
$\sqrt{\omega^{2}-m^{2}}\simeq\bar{\omega}$ and simultaneously setting $v_{\rm
burst}\to 1$, and ignoring the remaining prefactor. Similarly we can obtain
this scaling from Eq. (12) in the limit121212Here we use the identity
$\sqrt{2\pi}\delta(x)={\rm lim}_{\epsilon\to
0+}\exp(-x^{2}/2\epsilon^{2})/\epsilon$. $\delta\omega\to 0$,
$\omega\to\bar{\omega}$.

Using Eq. (21), we can estimate a typical ratio of axion background flux to
that expected of cold axion DM. Given a local DM density of $\rho_{\rm
lDM}\simeq 0.4$ GeV/cm3 and velocity dispersion $v_{\rm lDM}\simeq 10^{-3}$
(see e.g. Gelmini:2015zpa  for review), the local DM flux is

$\left(\frac{d\phi}{d\omega}\right)_{\rm lDM}\sim\frac{\rho_{\rm
lDM}}{m^{2}}\,v_{\rm lDM}\simeq 10^{26}\,{\rm eV}^{-1}{\rm cm}^{-2}\,{\rm
sec}^{-1}\left(\frac{v_{\rm lDM}}{10^{-3}}\right)\left(\frac{10^{-10}\,{\rm
eV}}{m}\right)^{2}~{}.$ (22)

Comparing Eqs. (21-22), the expected ratio of D$a$B flux to DM flux to be

$\frac{\displaystyle{d\phi/d\omega}}{\displaystyle{\left(d\phi/d\omega\right)_{\rm
lDM}}}\sim\left(\frac{v_{\rm burst}}{v_{\rm
lDM}}\right)\left(\frac{m}{\bar{\omega}}\right)^{2}\left(\frac{\rho_{\rm
loss}}{\rho_{\rm loc}}\right)\simeq 3\cdot
10^{-3}{\mathcal{F}}\left(\frac{m}{\bar{\omega}}\right)^{2}\,.$ (23)

Compared to the local cold DM flux, the D$a$B flux is boosted by its
relativistic velocity, $v_{\rm burst}/v_{\rm lDM}\sim 10^{3}$. However, it is
suppressed by the small average DM density of the Universe relative to the
expected local DM density, and by ${\mathcal{F}}<1$, as $\rho_{\rm
loss}/\rho_{\rm loc}\sim 10^{-6}{\mathcal{F}}$. Finally, for a fixed total
energy output the D$a$B is suppressed by the larger energy per particle,
$\omega$, as $(m/\bar{\omega})^{2}\ll 1$. Overall, the D$a$B flux is expected
to be at least a few orders of magnitude below that of local cold axion DM in
the solar neighborhood. We stress that these naive estimates do not account
for redshift dependencies. Since the energy densities are larger at earlier
times, converting a fixed axion DM fraction ${\mathcal{F}}$ one could convert
more axion DM by populations of transient sources at higher redshifts. Hence,
Eq. (23) should be modified by $(1+\bar{z})^{3}$ scaling when the source
distribution is sharply peaked at some redshift $\bar{z}$ in order to properly
capture such behavior.

Figure 4:  The differential flux, $d\phi/d\omega$, of D$a$B relative to the
flux of cold axion DM in the local neighborhood, $(d\phi/d\omega)_{\rm lDM}$,
as a function of $\omega/m$. For the source of the axion background, the thick
(dashed) lines correspond to a Gaussian energy distribution at the source with
peak energy $\bar{\omega}=3m$ and width $\delta\omega=m$ ($m/10$). The
cosmological rate function $R_{\rm burst}$ is given by: (Left panel) a power
law $\propto(1+z)^{p}$ with $p=-1,1$ and $z_{\rm max}=5,20$ as illustrated;
(Right panel) a Gaussian $\propto\exp\left[-(z-\bar{z})^{2}/\delta
z^{2}\right]$ with $\bar{z}=1,10$ and $\delta z=0.1,3$ as illustrated.

We illustrate the relativistic axion flux considering reference parameters
(see Sec. 2.2), normalized to the local cold axion DM expectation as in Eq.
(22), in Fig. 4. Here, for our estimates we assume a DM fraction
${\mathcal{F}}=0.1$, and a Gaussian emission spectrum with a central energy
$\bar{\omega}=3m$ and with a width $\delta\omega=m$ (dashed) or
$\delta\omega=m/10$ (thick). The left panel displays the results for a power-
law distribution of sources with an exponent $p=1$ (red) or $p=-1$ (purple),
and $z_{\rm max}=5,20$. For comparison, we also show the result for a Gaussian
distribution with $\delta z=0.1$ (red) or $\delta z=3$ (purple), and
$\bar{z}=1,10$ as labeled. We note that the D$a$B flux is within a few orders
of magnitude of the flux expected from the local axion DM density, matching at
the order-of-magnitude level to the rough estimates of Eq. (23), with a spread
in energy that depends on the $\delta\omega$ as well as the cosmological
history encoded in $R_{\rm burst}(z)$.

Fig. 4 illustrates how the parameters of the cosmological rate function and
burst spectrum determine the present-day flux in the D$a$B. For example, in
the left panel (power-law rate function), extending the duration of the burst
to higher redshifts, from $z_{\rm max}=5$ (faint lines) to $z_{\rm max}=20$
(opaque lines) leads to a more significant low-energy component. The power-law
exponent $p$, in contrast, has a more subtle effect on the present-day
spectrum. For a Gaussian rate function, the peak energy after redshifting is
found near $\bar{\omega}/(1+\bar{z})$, and when the Gaussian width in redshift
is smaller more of the energy remains near the peak, as in the red lines
($\delta z=0.1$) compared to the purple ones ($\delta z=3$). As we have
demonstrated, at high redshift the DM density is higher than today, allowing a
D$a$B flux larger than the DM background for example in case of Gaussian
distribution at high redshift (shown by faint lines in the right panel of Fig.
4). In both panels, the dashed lines, which correspond to a larger width in
energy $\delta\omega$, spread significantly in energy away from the initial
peak at $\bar{\omega}=3m$ than the solid lines, which exhibit a sharper
cutoff.

Our analysis indicates that D$a$B flux could be sub-leading, and in some
scenarios even dominate, compared to that of local cold axion DM. Successful
detection of D$a$B flux would prove particularly valuable, providing unique
information regarding the evolution of the axion field as well as its SM
interactions in cosmological and astrophysical contexts dramatically distinct
from cold DM. We further comment on direct-detection possibilities for D$a$B
in Sec. 6.2.

##  3 Axion conversion to photons

Complementary to direct detection, the D$a$B can can be probed through
indirect signatures due to axion interactions with SM constituents. Prominent
signals appear from axion conversion through axion-photon interaction of Eq.
(2), either in the presence of magnetic fields or through axion decays. We
calculate the rates of axion-photon conversion for each scenario. The
resulting photon signals can be searched for in a wide range of telescope and
satellite experiments outlined in Sec. 4. Additional detectable D$a$B
signatures could arise from other axion couplings as well, as we briefly
discuss in Sec. 6.

###  3.1 Magnetic fields

Propagating axions, including those contributing to D$a$B from various sources
as discussed in the previous section, can be efficiently converted (oscillate)
to photons $a\to\gamma$ by astrophysical magnetic fields Mirizzi:2006zy . We
shall subsequently focus on Milky Way Galactic magnetic fields. Consider
axions traveling a distance $R$ in a homogeneous magnetic field with a
transverse component in the plane normal to direction of propagation $B_{T}$,
which is of order $\mu{\rm G}$ in the Milky Way Beck:2008ty . We assume Milky
Way’s magnetic field to be $B_{T}=1$ $\mu{\rm G}$ throughout. Then, the
probability of conversion is Raffelt:1987im

$P_{\gamma\to a}=\left(\Delta_{a\gamma}R\right)^{2}\frac{\sin^{2}(\Delta_{\rm
osc}R/2)}{(\Delta_{\rm osc}R/2)^{2}}~{},$ (24)

where the oscillation wavenumber is $\Delta_{\rm
osc}\equiv\sqrt{(\Delta_{a}-\Delta_{\rm pl})^{2}+4\Delta_{a\gamma}^{2}}$. We
have considered the ultra-relativistic approximation, $\omega\gg m$, which we
expect to be a generally reasonable approximation131313In other energy
regimes, we expect additional corrections.. There are three distinct
contributions to the conversion probability:

$\displaystyle\Delta_{a\gamma}$
$\displaystyle\equiv\frac{g_{a\gamma}B_{T}}{2}\simeq 1.5\cdot
10^{-4}\left(\frac{g_{a\gamma}}{10^{-12}\,{\rm
GeV}^{-1}}\right)\left(\frac{B_{T}}{\mu{\rm G}}\right){\rm kpc}^{-1}\,,$ (25)
$\displaystyle\Delta_{a}$
$\displaystyle\equiv-\frac{m^{2}}{2\omega}\simeq-7.8\cdot
10^{13}\left(\frac{m}{10^{-11}\,{\rm eV}}\right)^{2}\left(\frac{10^{-10}\,{\rm
eV}}{\omega}\right){\rm kpc}^{-1}\,,$ (26) $\displaystyle\Delta_{\rm pl}$
$\displaystyle\equiv-\frac{\omega_{\rm pl}^{2}}{2\omega}\simeq-1.1\cdot
10^{13}\left(\frac{n_{e}}{10^{-2}\,{\rm
cm}^{-3}}\right)\left(\frac{10^{-10}\,{\rm eV}}{\omega}\right){\rm
kpc}^{-1}\,,$ (27)

where $\omega_{\rm pl}=\sqrt{4\pi\alpha n_{e}/m_{e}}\simeq 4\cdot 10^{-11}{\rm
eV}\sqrt{n_{e}/{\rm cm}^{-3}}$ is the plasma frequency corresponding to
effective mass that photons acquire induced by the SM background, with $m_{e}$
being the electron mass and $n_{e}$ density of free electrons in medium.

For lower energy emitted axions, when the energy scales could be semi-
relativistic with $\omega$ are not very far from $m$, we find
$\Delta_{a\gamma}$ is negligible at all relevant $m$. The other two
contributing terms are equal, $\Delta_{a}\simeq\Delta_{\rm pl}$, when
$m=\omega_{\rm pl}$. As indicated by Eqs. (26-27), both
$|\Delta_{a}|,|\Delta_{\rm pl}|\gg{\rm kpc}^{-1}$ in this case. Hence,
averaging the resulting fast oscillations in Eq. (24) as $\sin^{2}\to 1/2$, we
can write

$\displaystyle P_{\gamma\to a}$
$\displaystyle\simeq(\Delta_{a\gamma}R)^{2}\dfrac{\sin^{2}(|\Delta_{a}-\Delta_{\rm
pl}|R/2)}{(\Delta_{a}-\Delta_{\rm pl})^{2}R^{2}/4}$
$\displaystyle\simeq\begin{cases}\displaystyle
2\left(\dfrac{\Delta_{a\gamma}}{\Delta_{\rm pl}}\right)^{2}\\\ \displaystyle
2\left(\dfrac{\Delta_{a\gamma}}{\Delta_{a}}\right)^{2}\end{cases}$ (28)
$\displaystyle=\begin{cases}2\cdot
10^{-32}\left(\dfrac{g_{a\gamma}}{10^{-12}\,{\rm
GeV}^{-1}}\right)^{2}\left(\dfrac{\omega}{10^{-10}\,{\rm
eV}}\right)^{2}\left(\dfrac{10^{-2}\,{\rm
cm}^{-3}}{n_{e}}\right)^{2}\left(\dfrac{B_{T}}{\mu{\rm
G}}\right)^{2}\,,&m\lesssim\omega_{\rm pl}\,,\\\ 4\cdot
10^{-34}\left(\dfrac{g_{a\gamma}}{10^{-12}\,{\rm
GeV}^{-1}}\right)^{2}\left(\dfrac{\omega}{10^{-10}\,{\rm
eV}}\right)^{2}\left(\dfrac{10^{-11}\,{\rm
eV}}{m}\right)^{4}\left(\dfrac{B_{T}}{\mu{\rm
G}}\right)^{2}\,,&m\gtrsim\omega_{\rm pl}~{}.\end{cases}$

Note that the result is essentially independent of $R$, since the argument of
the sine function in Eq. (24) is significantly larger than unity for all
astrophysical, and certainly cosmological, distances. In Fig. 5, we illustrate
the two limits of Eq. (3.1) (red and blue curves, respectively), and compare
to the full expression of Eq. (24) (purple). As indicated, the conversion
probability is proportional to $P_{\gamma\rightarrow a}\propto
g_{a\gamma}^{2}\omega^{2}$.

In the other limiting regime, when $m\ll 10^{-11}$ eV and $\omega\gg$MeV (as
in e.g. environment of supernovae), the first term dominates and $P_{\gamma\to
a}$ becomes independent of $m$ and $\omega$ Mirizzi:2006zy ; Calore:2020tjw .
In that case

$P_{\gamma\to a}\simeq(\Delta_{a\gamma}R)^{2}\simeq 2\cdot
10^{-6}\left(\frac{g_{a\gamma}}{10^{-12}\,{\rm
GeV}^{-1}}\right)^{2}\left(\frac{B_{T}}{\mu{\rm
G}}\right)^{2}\left(\frac{R}{\rm kpc}\right)^{2}\,,$ (29)

which approaches unity $P_{\gamma\to a}\to 1$ for sizable distances $R$. In
panel (c) of Fig. 5 we illustrate the conversion probability for $R=$ kpc and
for different axion energies $\omega=$ keV, MeV, and GeV (purple, brown, and
green lines respectively).

(a)

(b)

(c)

Figure 5:  Conversion probability for axion to photon $a\to\gamma$ as a
function of axion mass $m$, in the (transverse) Galactic magnetic field
$B_{T}=\mu{\rm G}$ considering an electron density of $n_{e}=10^{-2}\,{\rm
cm}^{-3}$. The probability is computed for fixed $\omega$ equal to $10m$
(panel a), $\rm meV$ (panel b), or $\gg{\rm eV}$ (panel c, energies keV, MeV,
GeV as labeled). The purple, brown, and green lines show the full solutions of
Eq. (24), whereas the dashed red and blue fit lines illustrate the limiting
cases at small and large $m$ (see text for details). The axion-photon coupling
is fixed in the thick (faint) lines to a value of $g_{a\gamma}=10^{-12}\,{\rm
GeV}^{-1}$ ($10^{-15}\,{\rm GeV}^{-1}$).

Using the above, in Sec. 4 below we estimate the sensitivity of experiments to
detect photon signals generated by axion-photon conversion. For this, the flux
of photons can be readily computed from D$a$B flux as

$\displaystyle\frac{d\phi_{\gamma}}{d\omega}\Big{|}_{B{\rm-
field}}=P_{\gamma\to a}\frac{d\phi}{d\omega}$ (30)

where the axion flux $d\phi/d\omega$ is given by Eq. (12). In what follows, we
consider characteristic $R=1\,{\rm kpc}$, as the typical coherence length for
Galactic magnetic fields is ${\mathcal{O}}({\rm few})\,{\rm kpc}$ Beck:2008ty
.

Figure 6: Axion-photon decay probability factor, $P_{\rm decay}$, as a
function of axion mass $m$ at fixed energy $\omega$ (left panel), and as a
function of axion energy $\omega$ at fixed mass $m$ (right panel).
Representative values of axion coupling $g_{a\gamma}=10^{-12}\,{\rm GeV}^{-1}$
and cosmological propagation distance $R={\rm Gpc}$ are considered.

###  3.2 Axion decay

Axions can decay to photons, $a\to\gamma\gamma$, through the axion-photon
coupling in Eq. (2). The decay rate of an axion in vacuum is given by

$\Gamma_{\gamma\gamma}=\frac{g_{a\gamma}^{2}m^{3}}{64\pi}~{}.$ (31)

Hence, axions of energy $\omega$ traveling a distance $R$ will decay with
probability $P_{\rm decay}(R,\omega)=1-\exp(-R/\ell(\omega))$, where the mean
free path is given by

$\ell(\omega)\simeq\frac{\gamma v_{\rm
burst}}{\Gamma_{\gamma\gamma}}\simeq\left(\frac{\omega}{m}\right)\frac{64\pi}{g_{a\gamma}^{2}m^{3}}\simeq{\rm
Mpc}\left(\frac{\omega}{{\rm MeV}}\right)\left(\frac{100\,{\rm
keV}}{m}\right)^{4}\left(\frac{10^{-12}\,{\rm
GeV}^{-1}}{g_{a\gamma}}\right)^{2}\,,$ (32)

with $\gamma$ is the relativistic Lorentz factor. As Eq. (32) illustrates,
axions in the D$a$B can travel cosmological distances before decaying, and
therefore the photon flux from decay must be computed self-consistently over
cosmological redshifts $z$ as

$\frac{d\phi_{\gamma}}{d\omega}\Big{|}_{\rm
decay}=\int_{0}^{\infty}dz(1+z)\frac{dN_{a}(\omega(1+z))}{d\omega}R_{\rm
burst}(z)\left|\frac{dt}{dz}\right|2P_{\rm decay}(R(z),\omega(1+z))~{},$ (33)

which accounts for the dependence of the source distance $R$ and emission
energy $\omega$ on redshift $z$. For a sharply-peaked distribution around some
redshift $\bar{z}$, we can treat $P_{\rm decay}$ as a constant in the
integration, and the decay flux takes the form

$\displaystyle\frac{d\phi_{\gamma}}{d\omega}\Big{|}_{\rm decay}\simeq 2P_{\rm
decay}\left(R(\bar{z}),\omega(1+\bar{z})\right)\frac{d\phi}{d\omega}~{},$ (34)

which is analogous to Eq. (30) for magnetic field conversion. Note that there
is an additional factor of 2 compared to magnetic field conversion case, since
each axion decay creates two photons.

For example for Galactic contributions to D$a$B, such as from astrophysical
supernovae, the distances to the sources $R$ could be related to their
astrophysical population distributions, e.g. localized around Galactic Center.
More generally, particularly for dark sector sources, the source distances
$R(z)$ can be cosmological in scales, e.g. for $z\gtrsim 1$ we have
$R(z)\gtrsim{\rm Gpc}$, but depend in detail on the population distribution of
axion emission sources. We define an effective redshift for the purpose of
estimating the typical source distance in generality, $z_{\rm eff}$, as

$z_{\rm eff}\equiv\dfrac{\int dt\dfrac{z\,R_{\rm burst}(z)}{(1+z)^{3}}}{\int
dt\dfrac{R_{\rm burst}(z)}{(1+z)^{3}}}\,.$ (35)

For a peaked distribution $R_{\rm burst}(z)$, such as the Gaussian introduced
in Eq. (7), Eq. (35) gives $z_{\rm eff}\simeq\bar{z}$ to a good approximation
if $\bar{z}\gtrsim 5$. For the power-law distribution in Eq. (7), which can be
broader in redshift, $z_{\rm eff}$ is somewhat smaller, and we can approximate
$z_{\rm eff}\simeq 2/(7-2p)$. To estimate the decay probability in what
follows, we use $R(z)\simeq 1/H(z_{\rm eff})$ where $H(z)$ is the Hubble
parameter.

As noted by Ref. Calore:2020tjw , for large axion masses and energies
$\omega,m\gtrsim{\rm MeV}$ (appropriate for e.g. heavy axions produced in
supernovae), nearly all axions are expected to decay on cosmological
distances, as the exponential in Eq. (34) approaches zero, or equivalently
when $R/\ell\gg 1$. On the other hand, when $R/\ell\ll 1$, we can approximate
the flux as $d\phi_{\gamma}/d\omega\simeq(2R/\ell)d\phi/d\omega$, which is
proportional to $\omega^{-1}m^{4}$ multiplied by the axion flux. We illustrate
the behavior of $P_{\rm decay}$ in different limits in Fig. 6.

##  4 Observable signatures with photons

The generation of photon flux from D$a$B, either through magnetic $B$-field
conversion or axion decays, allows to search for axions over a broad parameter
space with a range of detectors and telescopes. In principle, the
contributions from these processes are additive, stemming from the same axion-
SM interaction $g_{a\gamma}$. However, as we demonstrated, the two processes
are expected to contribute a significant flux of photons in distinct kinematic
regimes. Approximately, axion magnetic field conversion tends to dominate at
high energies when $\omega\gg m$ whereas axion decays are more prominent when
their energy $\omega$ is within a few orders of magnitude of the mass $m$.

Figure 7:  Photon flux from Galactic magnetic $\mu$G $B$-field conversion of
the D$a$B (black lines) assuming ${\mathcal{F}}=0.1$ (dashed) or $10^{-10}$
(thick), for axions with mass $m=10^{-10}\,{\rm eV}$ and coupling
$g_{a\gamma}=10^{-12}$ GeV-1 emitted from transient sources with cosmological
burst rate $R_{\rm burst}\propto(1+z)^{-1}$, ratio of central energy spread to
average magnitude $\delta w/\bar{\omega}=1/10$. Each pair of black lines
corresponds to burst sources with central energy $\bar{\omega}$ as labeled
along the x-axis. Contributing background fluxes for a wide range of photon
frequencies (from radio, microwave, infrared, optical, and UV) are reproduced
from Ref. Cooray:2016jrk  (see also references therein). At higher energies
$\omega\gtrsim 100\,{\rm eV}$, we display the photon flux measured by specific
experiments, which lead to constraints on D$a$B, including NuSTAR
NuSTAR:2013yza ; Roach:2022lgo , XMM-Newton Foster:2021ngm , INTEGRAL
Bouchet:2008rp , COMPTEL 2008ICRC….3.1085S ; Strong:2019yfj , EGRET
Strong:2004de , and Fermi-LAT Fermi-LAT:2014ryh ; Principe:2022ttq .

In Fig. 7 and Fig. 8, we illustrate the photon flux obtained from Galactic
magnetic $B$-field conversion as well as axion decays, respectively. For
comparison, we include the measured background flux from telescope experiments
across a wide range of frequencies, from gamma-ray to radio, spanning nearly
20 orders of magnitude in photon energy. We observe that over a very broad
parameter space, the photon flux obtained from D$a$B can greatly exceed other
contributing backgrounds, allowing to impose constraints and promising future
detection prospects.

Highlighting our results, we estimate in detail the sensitivity of several
distinct experiments below. For concreteness, we consider a power-law
distribution of cosmological transient sources $R_{\rm burst}\propto(1+z)^{p}$
taking $p=-1$ and assuming a maximum redshift of $z_{\rm max}=20$.
Computations of results for the case of power-law distribution of cosmological
transient sources but with $p=1$, as well as a Gaussian distribution, are
depicted in App. A. We find the results for these cases to be overall similar
to those presented below for power-law distribution with $p=-1$, suggesting
that dependence on the detailed cosmological assumptions about transient
source population is not very strong. Further, we find that the spread of the
axion energy distribution, $\delta\omega/\bar{\omega}$, determines the
position of the low-energy cutoff of each black curve.

###  4.1 Limits from X-ray and gamma-ray searches

Figure 8: Photon flux from decay of axions from D$a$B (black lines) assuming
${\mathcal{F}}=0.1$ (dashed) or $10^{-10}$ (thick), considering axion coupling
$g_{a\gamma}=10^{-12}$ GeV-1 and mass $m=10^{-7}\,{\rm eV}$ (left),
$10^{-2}\,{\rm eV}$ (center), and $10\,{\rm keV}$ (right). D$a$B is considered
to arise from transient sources with cosmological burst rate of $R_{\rm
burst}\propto(1+z)^{-1}$, assuming typical source distance of $R(z_{\rm
eff})\simeq 1/H(z_{\rm eff})$ with $z_{\rm eff}\simeq 2/9$. Each pair of black
lines corresponds to transient sources with central energy $\bar{\omega}$ as
labeled along the x-axis. (right) Photon fluxes measured by specific
experiments and leading to constraints on D$a$B are as described in the
caption of Fig. 7. At lower energies (below MeV, center, and below eV, left),
we illustrate sensitivity of future searches by Vera C. Rubin Observatory
(formerly LSST) Tyson:2002mzo ; Mao:2022fyx , JWST JWST:2023abc ;
Janish:2023kvi ; Roy:2023omw , XRISM XRISMScienceTeam:2020rvx  and SKA
Weltman:2018zrl ; Braun:2019gdo ; Ghosh:2020hgd  as described in Sec. 4.2 (see
also Fig. 10).

Variety of experiments can sensitively probe flux of energetic photons with
$\omega\gg{\rm eV}$ from D$a$B. We illustrate this in Fig. 7 and Fig. 8, for
selected input parameters, using the measured photon background flux of
several existing experiments including NuSTAR NuSTAR:2013yza ; Roach:2022lgo ,
XMM-Newton Foster:2021ngm , INTEGRAL Bouchet:2008rp , COMPTEL
2008ICRC….3.1085S ; Strong:2019yfj , EGRET Strong:2004de , and Fermi-LAT
Fermi-LAT:2014ryh ; Principe:2022ttq . As a characteristic example of Fermi-
LAT, we consider 8-yr dataset observations of isotropic diffuse gamma-ray
background141414ULTRACLEANVETO event class, with reduced cosmic-ray
contamination. See iso P8R3 ULTRACLEAN V3 v1.txt at
https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html., which
above can be modeled using a simple fit function Calore:2020tjw

$\frac{d\phi_{\gamma}(\omega)}{d\omega}\Bigg{|}_{\rm Fermi}\simeq 5.5\cdot
10^{-10}\left(\frac{\rm GeV}{\omega}\right)^{2.2}\,{\rm MeV}^{-1}{\rm
cm}^{-2}{\rm sec}^{-1}$ (36)

in the $\omega$ energy range $30\,{\rm MeV}-300\,{\rm GeV}$. Analogously,
COMPTEL observations 2008ICRC….3.1085S  of diffuse gamma-ray emission can be
fit as Calore:2020tjw

$\frac{d\phi_{\gamma}(\omega)}{d\omega}\Bigg{|}_{\rm COMPTEL}\simeq 5\cdot
10^{-3}\left(\frac{{\rm MeV}}{\omega}\right)^{2.4}{\rm MeV}^{-1}{\rm
cm}^{-2}{\rm sec}^{-1}$ (37)

in the $\omega$ energy range $0.8-30\,{\rm MeV}$. We summarize observations
employed in our analysis in Tab. 1. Overall, we analyze combined observations
in the $\omega$ energy range of approximately ${\rm MeV}-{\rm 300~{}GeV}$
showcasing their implications for D$a$B, and leave a more comprehensive study
of signals outside of this range to future work. Other measured fluxes,
including those of EGRET, INTEGRAL, XMM, NuSTAR, and low-energy telescopes,
are summarized in Tab. 1.

Experiment(s) | Energy Range | Flux $(d\phi_{\gamma}/d\omega)$ Fit [eV-1cm-2sec-1]
---|---|---
Radio* Cooray:2016jrk  | $4-340\,{\rm neV}$ | $6.2\cdot 10^{14}\left(\dfrac{\mu{\rm eV}}{\omega}\right)^{1.5}$
Microwave* Cooray:2016jrk  | $67\,{\rm neV}-730\,\mu{\rm eV}$ | $3.2\cdot 10^{17}\left(\dfrac{\omega}{\rm meV}\right)^{0.5}\exp\left(-\dfrac{15\omega}{\rm meV}\right)$
Infrared* Cooray:2016jrk  | $0.1-5\,{\rm meV}$ | $6.3\cdot 10^{13}\left(\dfrac{\omega}{\rm meV}\right)^{0.5}\exp\left(-\dfrac{1.6\omega}{\rm meV}\right)$
Optical* Cooray:2016jrk  | $10\,{\rm meV}-1.7\,{\rm eV}$ | $4.4\cdot 10^{6}\left(\dfrac{\rm eV}{\omega}\right)^{1.8}$
UV* Cooray:2016jrk  | $46-58\,{\rm eV}$ | $1.7\cdot 10^{4}\left(\dfrac{\rm eV}{\omega}\right)^{2}$
NuSTAR NuSTAR:2013yza ; Roach:2022lgo  | $0.3-200\,{\rm keV}$ | $10^{-4}\left(\dfrac{\rm keV}{\omega}\right)$
XMM-Newton Foster:2021ngm  | $1-8\,{\rm keV}$ | $8.9\cdot 10^{-6}\left(\dfrac{\rm keV}{\omega}\right)^{0.65}$
INTEGRAL Bouchet:2008rp  | $24\,{\rm keV}-2\,{\rm MeV}$ | $10^{-8}\left(\dfrac{\rm MeV}{\omega}\right)^{1.6}$
COMPTEL 2008ICRC….3.1085S ; Strong:2019yfj  | $0.8-30\,{\rm MeV}$ | $5\cdot 10^{-9}\left(\dfrac{{\rm MeV}}{\omega}\right)^{2.4}$
EGRET Strong:2004de  | $6-240\,{\rm MeV}$ | $4\cdot 10^{-16}\left(\dfrac{\rm GeV}{\omega}\right)^{2.3}$
Fermi-LAT Fermi-LAT:2014ryh ; Principe:2022ttq  | $30\,{\rm MeV}-300\,{\rm GeV}$ | $5.5\cdot 10^{-16}\left(\dfrac{\rm GeV}{\omega}\right)^{2.2}$

* Denotes combined observations from multiple experiments Cooray:2016jrk .

Table 1: Summary of experiments whose observations are used to analyze photon
signals from D$a$B derived in this work, along with their sensitivity energy
range. The fits to the total measured photon flux in the Milky Way by
experiments are also displayed. See Fig. 7 for further details.

In order to obtain constraints on axion interactions, we vary the coupling
$g_{a\gamma}$ to derive the sensitivity regions for Fermi-LAT and COMPTEL,
highlighting our analysis, where the peak photon fluxes from the D$a$B exceeds
the gamma-ray fluxes displayed in Tab. 1. We stress that our analysis can be
readily analogously carried out for other experimental observations (see e.g.
Tab. 1). The resulting limits and constraints are illustrated in Fig. 9. The
lines correspond to searches over the range of energies given in the legend,
for COMPTEL (blue and red lines) and Fermi-LAT (green, yellow, and brown
lines). As discussed previously, Galactic $B$-field conversion (see Eq. (30))
is most relevant for small axion masses $m\ll{\rm eV}$, whereas decays (see
Eq. (33)) dominate at larger axion masses. The change in the slope of the
lines correspond to shifts in the behavior of the photon conversion and decay
probabilities, seen in Fig. 5 and Fig. 6, respectively. For comparison, we
also show the band corresponding to the QCD axion in pink Kim:1979if ;
Shifman:1979if ; Dine:1981rt ; Zhitnitsky:1980tq .

The total flux is proportional to $\propto g_{a\gamma}^{2}{\mathcal{F}}$,
allowing for simple rescaling of our benchmarks for different DM fraction
${\mathcal{F}}$. We illustrate in Fig. 9 scenarios of ${\mathcal{F}}=0.1$ and
${\mathcal{F}}=10^{-20}$. We observe that even for a very small converted DM
fraction ${\mathcal{F}}\lesssim 10^{-20}$, we find the estimated sensitivity
for axions searches converting to photons due to Galactic magnetic fields far
exceeds existing constraints from astrophysics or conventional direct DM
detection searches using haloscopes (see Fig. 9 caption for details). At large
masses, we find signals from axion decays can be competitive with strong
existing cosmological limits Cadamuro:2011fd .

We note that our estimates are conservative, since the photon flux associated
with D$a$B could be detectable even if it is subdominant to background. Our
analysis thus motivates dedicated searches, which would include detailed
information about the distribution of D$a$B photons as a function of energy
and could potentially significantly further improve these results.

Figure 9: Sensitivity to D$a$B photons of COMPTEL (blue, red) and Fermi-LAT
(green, yellow, and brown), as a function of axion mass $m$ for average
transient burst energy $\bar{\omega}$. The emitted photons are produced from
either Galactic magnetic $\mu$G $B$-field conversion (left side) or axion
decays (right). D$a$B is considered to arise from transient sources with
cosmological burst rate of $R_{\rm burst}\propto(1+z)^{-1}$. For $B$-field
conversion we take $R={\rm kpc}$, and for axion decay we assume a typical
source distance of $R(z_{\rm eff})\simeq 1/H(z_{\rm eff})$ with $z_{\rm
eff}\simeq 2/9$. The upper (lower) lines correspond to ${\mathcal{F}}=0.1$
(${\mathcal{F}}=10^{-20}$). The thick (dashed) lines correspond to transient
burst with energy spread $\delta\omega=\bar{\omega}/10$
($\delta\omega=\bar{\omega}$). For enumeration of astrophysical limits (gray),
see Fig. 12. At large $m$ we also include limits from NS (NS) in the
Breakthrough Listen Galactic Center Survey Foster:2022fxn  and direct DM
detection searches using haloscopes (red) including ADMX ADMX:2018gho ;
ADMX:2018ogs ; ADMX:2019uok ; Crisosto:2019fcj ; ADMX:2021mio ; ADMX:2021nhd ,
CAPP Lee:2020cfj ; Jeong:2020cwz ; CAPP:2020utb ; Lee:2022mnc ; Kim:2022hmg ;
Yi:2022fmn ; Yang:2023yry , CAST-CAPP Adair:2022rtw , CAST-RADES CAST:2020rlf
, GrAHal Grenet:2021vbb , RBF DePanfilis , HAYSTAC HAYSTAC:2018rwy ;
HAYSTAC:2020kwv ; HAYSTAC:2023cam , QUAX Alesini:2019ajt ; Alesini:2020vny ;
Alesini:2022lnp , ORGAN McAllister:2017lkb ; Quiskamp:2022pks , and TASEH
TASEH:2022vvu , as well as leading constraints from cosmology (blue) from Ref.
Cadamuro:2011fd .

###  4.2 Sensitivity of future searches

Rich opportunities exist to probe photons from D$a$B over a broad range of
energies. Complementing observations at high energies $\omega\gg$ eV discussed
in previous section, at low energies $\omega\lesssim{\rm eV}$ there also exist
significant backgrounds from radio, microwave, infrared, and optical frequency
photons, as summarized in Fig. 7. However, many recently-launched as well as
proposed future experiments are expected to probe significantly smaller fluxes
over a wide range of energies. These include the Square Kilometer Array (SKA)
with sensitivity in the radio range $\omega<10^{-4}\,{\rm eV}$ Weltman:2018zrl
; Braun:2019gdo ; Ghosh:2020hgd , Vera C. Rubin Observatory151515Formerly,
Large Synoptic Survey Telescope, LSST. operating at optical energy range
$1-10\,{\rm eV}$ Tyson:2002mzo ; Mao:2022fyx , James Webb Space Telescope
(JWST) sensitivity in the infrared of $0.04-2\,{\rm eV}$ JWST:2023abc ;
Janish:2023kvi ; Roy:2023omw , recently launched XRISM with sensitivity for
X-rays of $0.3-12\,{\rm keV}$ XRISMScienceTeam:2020rvx , as well as the
e-ASTROGAM and AMEGO proposals with sensitivity for gamma rays in the range of
$0.1\,{\rm MeV}-300\,{\rm MeV}$, estimated as having similar sensitivity
e-ASTROGAM:2017pxr ; DeAngelis:2021abc ; AMEGO:2021abc ; Engel:2022bgx .

We estimate, see App. B for details, the expected photon background flux
observations from experiments and illustrate comparison with the D$a$B flux in
Fig. 10. We note that precise sensitivity estimates could depend on a
multitude of relevant factors, such as the source, telescope direction, noise
estimates as well as duration.

Figure 10: Same as Fig. 7, but illustrating the expected sensitivity of
recently-launched and future experiments (see App. B and references for flux
estimation) including SKA Weltman:2018zrl ; Braun:2019gdo ; Ghosh:2020hgd ,
Vera C. Rubin Observatory (formerly LSST) Tyson:2002mzo ; Mao:2022fyx , JWST
JWST:2023abc ; Janish:2023kvi ; Roy:2023omw , XRISM XRISMScienceTeam:2020rvx ,
and gamma-ray proposals e-ASTROGAM and AMEGO e-ASTROGAM:2017pxr ;
DeAngelis:2021abc ; AMEGO:2021abc ; Engel:2022bgx . Note that the blue arrows
labeled “Rubin Observatory” denote only a range of sensitivity in $\omega$,
not a photon flux sensitivity estimate.

To showcase different detection approaches, we highlight the sensitivity
estimate of SKA to the flux of photons from D$a$B. SKA is expected to have
unprecedented sensitivity reach at low radio energies $\omega<10^{-4}$ eV and
the methodology is qualitatively distinct from that of searches with high-
energy gamma-ray telescopes, such as Fermi-LAT as discussed in previous Sec.
4.1. Expected performance of SKA can be found in Ref. Braun:2019gdo . We
consider two planned realizations of the experiment: _SKA-low_ , with a
collecting surface area of $A\simeq 419,000\,{\rm m}^{2}$ and operating in the
$50-350~{}{\rm MHz}$ frequency range, and _SKA-mid_ , with a collecting
surface area of $A\simeq 10^{6}\,{\rm m}^{2}$ (i.e. 5659 antennas with 15 m
diameter) and detection range of $350~{}{\rm MHz}-15.4~{}{\rm GHz}$. For both
configurations we assume a detection efficiency of $\eta\simeq 0.8$.

The experimental sensitivity can be conventionally estimated via by comparing
the total power in the signal, $P_{\rm sig}$, to the power in the noise,
$P_{\rm noise}$. The signal power can be calculated from the differential flux
as

$P_{\rm sig}=\eta AS~{}.$ (38)

The integrated flux density, $S$, is defined by161616Conventionally, this is
measured in Janskys (Jy) per second. Recall that $1$ Jy $=10^{-23}$ erg
cm${}^{-2}\simeq 6.2\cdot 10^{-12}$ eV cm-2.

$\displaystyle
S\equiv\int_{\omega_{1}}^{\omega_{2}}d\omega\omega\,\frac{d\phi_{\gamma}}{d\omega}~{}.$
(39)

We estimate the flux density $S$ in Fig. 11 for several choices of average
transient burst emission energy $\bar{\omega}$.

Figure 11: (Left) Total integrated photon flux density $S$, as described by
Eq. (39). (Right) Detection signal-to-noise ratio SNR = $P_{\rm
sig}/\overline{P}_{\rm noise}$ (right panel) in the sensitivity range of SKA-
low (blue, red) as well as SKA-mid (green, yellow, brown), as a function of
axion mass $m$ for an average transient emission burst energy $\bar{\omega}$.
The thick (dashed) lines correspond to burst emission energy spread of
$\delta\omega=\bar{\omega}/10$ ($\delta\omega=\bar{\omega}$). D$a$B is
considered to arise from transient sources with cosmological burst rate of
$R_{\rm burst}\propto(1+z)^{-1}$ and with $\mathcal{F}=0.1$.

The system’s noise can be described by Dicke’s radiometric receiver formula
Dicke:1946glx

$P_{\rm noise}=T_{\rm sys}\sqrt{\frac{\Delta\nu}{t_{\rm obs}}}~{},$ (40)

where $\Delta\nu$ is the frequency bandwidth in the experiment, $t_{\rm obs}$
is the observation time, $T_{\rm sys}$ is the system noise temperature. We
consider the following contributions to the system noise temperature

$T_{\rm sys}=T_{\rm atm}+T_{\rm CMB}+T_{\rm bknd}+T_{r}~{},$ (41)

where $T_{\rm atm}\simeq 3\,{\rm K}$ Ajello:1995aaa  is the atmospheric sky
emission temperature, $T_{\rm CMB}=2.725\,{\rm K}$ is the cosmic microwave
background (CMB) radiation temperature, $T_{\rm bknd}$ is the sky brightness
temperature contributed by emission of all “background” radio sources and
$T_{r}$ is the noise contributed by radiometric receiver itself. For SKA-low
(mid), we take an approximate averaged contribution over frequency bands of
$T_{r}\simeq 20$ ($40$) K Braun:2019gdo . The background source emission
temperature can be estimated as Braun:2019gdo ; Ghosh:2020hgd

$T_{\rm bknd}\simeq 60\,{\rm K}\left(\frac{300\,{\rm
MHz}}{\nu}\right)^{2.55}~{}.$ (42)

Converting the noise power $P_{\rm noise}$ in Eq. (40) to an equivalent flux,
we arrive at the red lines in Fig. 7 and Fig. 8. As a simple noise model for
our study, we use the noise averaged over the frequency range of each
experiment

$\bar{P}_{\rm noise}\equiv\frac{1}{\Delta\nu}\int P_{\rm noise}(\nu)d\nu.$
(43)

For SKA-low (SKA-mid) we have $\bar{P}_{\rm noise}\simeq 0.3$ ($0.17$) eV/sec.

Figure 12: Sensitivity estimation for SKA-low (blue, red) and SKA-mid (green,
yellow, and brown) to detect D$a$B, as a function of axion mass $m$ for
average transient burst energy $\bar{\omega}$. The dominant axion-photon
conversion process in this region is through the Galactic magnetic B-fields
(see Eq. (30)). The thick (dashed) lines correspond to transient burst energy
spread of $\delta\omega=\bar{\omega}/10$ ($\delta\omega=\bar{\omega}$). D$a$B
is considered to arise from transient sources with cosmological burst rate of
$R_{\rm burst}\propto(1+z)^{-1}$ and with $\mathcal{F}=0.1$. Existing
constraints are overlaid in gray including astrophysical limits (only leading
limits shown, see e.g. Ref. AxionLimits ; Adams:2022pbo  for a more
comprehensive list) from Chandra (H1821-643 Reynes:2021bpe , M87 Marsh:2017yvc
, and NGC1275 Reynolds:2019uqt  sources), super star clusters Dessert:2020lil
, magnetic white dwarf (MWD) polarization Dessert:2022yqq , CAST CAST:2007jps
; CAST:2017uph , pulsar polar-cap cascades Noordhuis:2022ljw , and SN1987a
Payez:2014xsa ; Hoof:2022xbe , direct detection limits (red) from SHAFT
Gramolin:2020ict  and ABRA 10cm Salemi:2021gck  as well as cosmology-dependent
limit from axion star explosions inducing early-Universe heating of the
intergalactic medium shown with gray long-dashed line Du:2023jxh ;
Escudero:2023vgv .

Defining the signal-to-noise-ratio parameter ${\rm SNR}\equiv P_{\rm
sig}/\bar{P}_{\rm noise}$, we also illustrate the sensitivity region for SKA
for the same benchmarks in $\bar{\omega}$ in the right panel of Fig. 11. As
before, the slope of the lines changes when the behavior of the photon
conversion probability changes, near $10^{-12}-10^{-11}\,{\rm eV}$, as seen in
Fig. 5(b). In order to determine sensitivity to axions contributing to D$a$B,
as before we vary coupling $g_{a\gamma}$ values to find sensitivity regions
where SNR$>2$. The results are shown in Fig. 12 We observe that SKA can
efficiently detect D$a$B signatures for sources with central burst energies
$\bar{\omega}\simeq 10^{-8}-10^{-4}\,{\rm eV}$ over a wide range of axion
masses $m\lesssim 10^{-9}\,{\rm eV}$. We note that the precise sensitivity
depends on the assumed cosmological population distribution of transient burst
sources, as further discussed in App. A.

##  5 Case study: axion star bosenovae

As demonstrated, photons arising from D$a$B constitute promising targets for
future searches. Our analysis can hence be utilized to set constraints on and
probe the existence of transient sources of relativistic axion emission as
well as their population distributions originating from distinct theories of
new physics.

As an example, we outline the analysis steps for the scenario of relativistic
axions emitted from the bosenovae that are the endpoint evolution of
collapsing boson stars Eby:2016cnq ; Levkov:2016rkk , discussed in Sec. 2.2.2.
Importantly, the axion emission spectrum of the source in this case,
$dN_{a}(\omega)/d\omega$, has been precisely determined by dedicated numerical
simulations Levkov:2016rkk . We can model the leading relativistic peak as
shown in Fig. 3, which has central axion emission energy of
$\bar{\omega}\simeq 2.2m$ and energy spread of $\delta\omega\simeq 0.4m$,
using Gaussian distribution described in Eq. (10). On the other hand, the
cosmological history and transient rate of bosenovae remains poorly understood
and can depend on multiple complex processes. For illustration, we adopt a
power-law cosmological bosenovae source distribution of Eq. (7) with parameter
$p=-1$, and analyze the resulting constraint on ${\mathcal{F}}$ assuming non-
observation of signals in experiments.

With above approximations, one can compute the total flux in D$a$B using Eq.
(12). Assuming boson stars are formed from ultralight particles, $m\ll{\rm
eV}$, and the energy scale of emission is of order $\bar{\omega}\sim 2.2m$, we
will be interested in low-energy signals detectable by SKA. In this parameter
range, axion decays rather than Galactic magnetic field axion conversion
dominates the production of photons. Therefore, we employ Eq. (33) to
determine the total flux of photons, and Eq. (38) as well as Eq. (39) to
determine the total signal power.

Figure 13: Sensitivity estimate for photons from the D$a$B associated with
axion bosenovae, with emission spectrum modeled as Gaussian with average
energy $\bar{\omega}=2.2m$ and spread $\delta\omega=0.4m$, considering SKA-low
(blue lines) or SKA-mid (green). In this range, the flux is dominated by
photons produced from axion decay (see Eq. (33)). The solid lines correspond
to ${\mathcal{F}}=0.1$, while the dashed lines have ${\mathcal{F}}=10^{3}$.
D$a$B is considered to arise from bosenovae with cosmological burst rate of
$R_{\rm burst}\propto(1+z)^{-1}$. Constraints from astrophysics and haloscopes
are further described in the caption of Fig. 9.

We vary $g_{a\gamma}$ to determine the sensitivity region where ${\rm
SNR}\equiv P_{\rm sig}/\bar{P}_{\rm noise}>2$. The results are shown in Fig.
13 for SKA-low (blue) as well as SKA-mid (green). For the benchmark we used
above, ${\mathcal{F}}=0.1$, the best projected sensitivity in the SKA range is
$g_{a\gamma}\gtrsim 10^{-10}\,{\rm GeV}^{-1}$ at $m\simeq 10^{-5}\,{\rm eV}$,
which is weaker than existing astrophysical constraints (shaded gray region).
To reach the sensitivity in axion coupling of the order of $g_{a\gamma}\simeq
10^{-12}\,{\rm GeV}^{-1}$, D$a$B DM fraction would need to exceed
${\mathcal{F}}\simeq 10^{3}$, as can be readily seen from SNR scaling of
$\propto{\mathcal{F}}g_{a\gamma}^{2}$. For illustration, in Fig. 13 we also
display the sensitivity reach considering the unphysical benchmark of
${\mathcal{F}}=10^{3}$. The results are similar for other cosmological
bosenovae rate distribution functions, including Gaussian analyzed in App. A.

This example illustrates how D$a$B can be exploited to constrain the
distribution of relativistic transient axion bursts from a source population
with a known emission spectrum, in this case bosenovae. Such search also
represents a complementary approach to existing ones that can probe couplings
near the astrophysical limits. Based on the results above, we expect improved
sensitivity for more energetic axion emission bursts, with $\bar{\omega}\gg m$
(see Sec. 4.2), or for astrophysical sources giving rise to X-ray or gamma-ray
signals (see Fig. 7 and Fig. 8). Extrapolating this study to larger $m$, our
results also motivate high-sensitivity next-generation experiments at
frequencies above 15 GHz where SKA-mid loses sensitivity, which may lead to
novel constrains on the D$a$B. We leave detailed investigation of such sources
to future work.

##  6 Complementary search prospects

###  6.1 Axion-electron and other couplings

While we have focused on axion couplings to photons, other axion couplings
(see Eq. (3)) could give rise to different novel signals associated with D$a$B
and call for dedicated analyses.

As an example, consider the coupling to electrons of the form
${\mathcal{L}}\supset(g_{ae}/2m_{e})(\partial_{\mu}a)\bar{e}\gamma^{\mu}\gamma_{5}e$,
where $m_{e}$ and $e$ are the electron mass and field, and $g_{ae}$ is the
dimensionless axion-electron coupling. If the axions have mass greater than
twice the electron mass, $m\gtrsim 2m_{e}\simeq{\rm MeV}$, this creates a new
decay channel for the D$a$B and could be searched for in a wide range of
experiments. Note that if the electron coupling is present at tree level,
there is no analogue of magnetic-field conversion that we described for photon
couplings. However, photon couplings could be induced beyond tree level, e.g.
through a loop of electrons with tree-level coupling given by Eq. (3),
implying a potentially interesting interplay between the two Calore:2020tjw .

The typical decay length for the process $a\to\bar{e}e$ is given by
Altmann:1995abc ; Lucente:2021hbp

$\ell_{e}=\dfrac{\gamma v_{\rm burst}}{\Gamma_{e^{+}e^{-}}}=8\cdot
10^{-5}\,{\rm
pc}\left(\dfrac{10^{-13}}{g_{ae}}\right)^{2}\left(\dfrac{20~{}{\rm
MeV}}{m}\right)^{2}\left(\dfrac{\omega}{100~{}{\rm
MeV}}\right)\left[\dfrac{1-\left(\dfrac{m}{\omega}\right)^{2}}{1-\left(\dfrac{2m_{e}}{m}\right)^{2}}\right]^{1/2}\,,$
(44)

where we used171717This also matches Ref. Bauer:2017ris , taking
$g_{ae}=c_{\ell\ell}^{\rm eff}m_{\ell}/\Lambda$.
$\Gamma_{e^{+}e^{-}}\simeq(g_{ae}^{2}/8\pi)\sqrt{m^{2}-4m_{e}^{2}}$. Compared
with the case of axion decays to photons, c.f. Eq. (32), the decay length
associated with axion electron decays tends to be significantly shorter, but
is restricted to sizable $m\gtrsim 2m_{e}$. The flux of electrons and
positrons from axion decays is similar to the photon case in Eq. (33), and is
given by

$\frac{d\phi_{e}}{d\omega}\Big{|}_{\rm
decay}=\int_{0}^{\infty}dz(1+z)\frac{dN_{a}(\omega(1+z))}{d\omega}R_{\rm
burst}(z)\left|\frac{dt}{dz}\right|2P^{(e)}_{\rm decay}(R(z),\omega(1+z))~{},$
(45)

where $P_{\rm
decay}^{(e)}(R,\omega)=\left[1-\exp\left(-R/\ell_{e}(\omega)\right)\right]$ is
they decay probability. The flux reduces, as in Eq. (34), to

$\frac{d\phi_{e}}{d\omega}\Big{|}_{\rm decay}\simeq 2P_{\rm
decay}^{(e)}(R(\bar{z}),\omega(1+\bar{z}))\frac{d\phi}{d\omega}~{},$ (46)

when $P_{\rm decay}^{(e)}$ is independent of redshift, e.g. for a rapid epoch
of bursts at a fixed redshift $z=\bar{z}$. Here we sum the flux of electrons
and positrons together.

Figure 14: Axion-electron decay probability, $P_{\rm decay}^{(e)}$, as a
function of axion mass $m$ and a fixed energy $\omega$ (left panel), and as a
function of energy $\omega$ at fixed mass $m$ (right panel). Reference values
of coupling $g_{ae}=10^{-18}$ and distance $R={\rm Gpc}$ are considered.

For cosmological distances $R\sim{\rm Gpc}$, the decay probability generally
tends to unity $P_{\rm decay}^{(e)}\to 1$ when the axion-electron coupling
values are near the current experimental limits for light axions,
$g_{ae}\simeq 10^{-13}$ (see e.g. Capozzi:2020cbu ). Hence, one expects nearly
all sufficiently heavy axions contributing to D$a$B to decay to electrons. For
illustration, in Fig. 14 we demonstrate that even for a significantly smaller
coupling $g_{ae}=10^{-18}$ a significant fraction of the D$a$B is expected to
decay, providing opportunities for related signatures.

Electrons and positrons produced by the D$a$B can serve as novel efficient
sources of photons. For astrophysical sources inside the Galaxy, primary
mechanisms include interactions with interstellar medium gas, Galactic
magnetic fields, or interstellar radiation fields including the cosmic
microwave background (CMB). These can lead to production of X-ray and
$\gamma$-ray photons, serving as targets for experiments like XMM-Newton,
INTEGRAL and COMPTEL, as has been shown for positrons originating in decays of
astrophysical supernova axions DelaTorreLuque:2023nhh ; DelaTorreLuque:2023huu
. Our analysis also highlights the significance of distinct sources producing
axions contributing to D$a$B and their energetics. As we have demonstrated,
for example, axion decay probability can be dramatically different between the
ultrarelativistic source emission with $\omega\gg m$ and semi-relativistic
source emission with $\omega\sim m$.

However, the D$a$B flux need not originate from within the Galaxy,
particularly in case of dark sector sources, allowing for other search
approaches as we outline. In analogy to the scenario of decaying DM (e.g.
Cirelli:2012ut ; Ando:2015qda ; Blanco:2018esa ), D$a$B contributions to
extragalactic $\gamma$-ray background can arise from Inverse Compton
scattering off the CMB. Unlike decaying DM case, the origin of decaying
particles is not associated with a constant DM density but depends on the
cosmological redshift distribution of axion emission sources as specified by
$R_{\rm burst}$ model. We can thus estimate the flux of photons of energy
$\omega$ in generality as

$\displaystyle\frac{d\phi_{\gamma}}{d\omega}\Big{|}_{\rm decay}=$
$\displaystyle~{}\int_{0}^{\infty}dz(1+z)R_{\rm
burst}(z)\left|\frac{dt}{dz}\right|$ (47) $\displaystyle\times\int
d\omega_{e}\int d\omega_{\rm
CMB}\int_{\omega_{e}}^{\infty}d\omega^{\prime}P_{\rm
decay}^{(e)}(R(z),\omega^{\prime}(1+z))\frac{dN_{a}(\omega^{\prime}(1+z))}{d\omega^{\prime}}Q(\omega,\omega^{\prime},\omega_{e},\omega_{\rm
CMB})~{},$

where $Q(\omega,\omega^{\prime},\omega_{e},\omega_{\rm CMB})$ is a function
characterizing the photon production from electron and positron conversion
originating from D$a$B associated with Inverse Compton scattering against CMB
photons. Here, the integrations are taken over the distribution of CMB
energies $\omega_{\rm CMB}$, electron energies $\omega_{e}$, and axion
energies $\omega^{\prime}$. The flux of photons will also experience
attenuation. We leave detailed computations of realizations of such scenarios
and discussion of resulting constraints for future work.

Other potentially promising scenarios for observations can also arise from
axions coupled to gluons or quarks, as is the case for QCD axions. In such
case, axions can decay hadronically. For $m\lesssim 1$ GeV, the dominant
decays are $a\rightarrow 3\pi^{0}$ and $a\rightarrow\pi^{+}\pi^{-}\pi^{0}$
Bauer:2017ris , with neutral pions subsequently decaying primarily into
photons $\pi^{0}\rightarrow\gamma\gamma$ and charged pions into muons and
neutrinos. We leave detailed discussion of these possibilities for future
studies.

###  6.2 Direct detection

While we have focused on indirect D$a$B searches, direct searches are possible
and complementary. Here, we briefly describe prospects and advantages of D$a$B
direct detection.

As computed in Sec. 2.3, the total flux from D$a$B is generally sub-leading
compared to the flux of local DM in the solar neighborhood vicinity. However,
in some scenarios, we found that D$a$B flux can even be larger than the
expected local DM flux. A direct search for D$a$B in terrestrial axion
experiments could be viable for a number of reasons. First, the frequency
range of interest is shifted from around axion mass $m$, associated with non-
relativistic cold axion DM described by Maxwell-Boltzmann distribution with
velocity dispersion of $v_{0}\simeq 220$ km/s (e.g. Gelmini:2015zpa  for
review), to $\omega\gg m$ for relativistic axions, allowing for additional
sensitivity to very light fields. Further, the broad distribution of axion
energies contributing to D$a$B is also distinct from the narrow distribution
expected of cold DM, see Fig. 4. It is also worth noting that the very local
DM density, on scales of $\sim$AU-pc within our solar neighborhood, is not
well-known Adler:2008rq ; Pitjev:2013sfa , and can deviate from the
conventionally-assumed value of $\rho_{\rm loc}$.

D$a$B direct detection searches are complementary to direct detection searches
of relativistic axions from transient sources, as analyzed in Ref.
Dailey:2020sxa ; Eby:2021ece ; Arakawa:2023gyq . Importantly, the presence of
D$a$B constitutes continuous signals with expected constant density in
experiments. On the other hand, transient signals result in temporary flux
density increases of axions traversing the experiments, which depending on
effects such as wave-spreading can be nearly instantaneous or lasting even
months or years Eby:2021ece ; Arakawa:2023gyq . While the D$a$B flux depends
on the population distribution of axion emission sources, persistent
associated signals do not require existence of a particular transient source
emitting axions recently or near the Earth. We also note searches for D$a$B
are complementary to _cosmological_ sources of relativistic axions, the so-
called cosmic axion background Conlon:2013isa ; Marsh:2013opc ; Dror:2021nyr ,
which are limited by the number of relativistic degrees of freedom in the
early Universe, $N_{\rm eff}$, and also are diminished by significant
redshifting of axion energies when $z\gtrsim 1000$.

For axions of mass $m\lesssim{\rm eV}$, the integration of coherent
oscillations over long timescales poses a significant source of experimental
sensitivity. For such oscillations this is characterized by a coherence
timescale $\tau\simeq 2\pi(mv^{2})^{-1}$ for axion velocity $v$, which
corresponds to an oscillator quality factor $Q\equiv v^{-2}$. For cold DM in
the solar neighborhood vicinity, $Q_{\rm lDM}\simeq v_{\rm lDM}^{-2}\simeq
10^{6}$, and for transient bosenovae the effective coherence can be nearly as
long due to wave-spreading effects Eby:2021ece ; Arakawa:2023gyq . The D$a$B,
however, is expected to be largely composed of incoherent waves from sources
across the Universe, implying $Q\simeq{\mathcal{O}}(1)$. This is also the case
for the cosmic axion background Conlon:2013isa ; Marsh:2013opc ; Dror:2021nyr
. These effects may pose challenges for direct detection, and call for further
investigations.

We further stress that the distribution of axion energies in the D$a$B can
encode unique information about the cosmological distribution of sources, see
Fig. 4. In the event of a successful cold axion DM detection, a dedicated
D$a$B search carries the potential to unlock substantial novel information
related to the cosmological production, SM couplings, and self-interactions of
axions. We leave detailed study of these topics for future work.

##  7 Conclusions

While conventional axion searches have primarily focused on cold axions
constituting DM, relativistic axions offer novel rich opportunities to probe
fundamental physics and theories. We have developed a general framework to
describe D$a$B formed from accumulation of emitted relativistic axions from
populations of distinct historic transient sources. Utilizing a small number
of key input parameters that describe the cosmological population and emission
properties of the sources, depending on the theory, our framework allows to
easily extract the energy density and flux of the resulting axions in the
vicinity of terrestrial and space-based experiments. We show that our analysis
allows to comprehensively and systematically explore a wide range of D$a$B
sources, include those related to astrophysics such as supernovae or NS
mergers, as well as originating from dark sector dynamics such as axion star
bosenova explosions.

For axions coupled to photons, we established sensitive constraints on the
D$a$B flux associated with light axions of mass $m\lesssim 10^{-3}\,{\rm eV}$
originating from energetic sources with $\omega\gtrsim{\rm MeV}$ considering
axion-photon conversion and subsequent photon detection within a range of
experiments such as COMPTEL, NuSTAR, XMM-Newton, INTEGRAL, EGRET and Fermi. We
demonstrate that future experiments, including SKA, JWST, XRISM, Vera C. Rubin
Observatory and AMEGO/e-ASTROGAM have the capability to probe axion-photon
couplings with remarkable sensitivity over multiple orders of magnitude in
parameter space in axion masses and couplings for sources emitting axions with
energies as low as $0.1~{}\mu$eV.

Focusing on relativistic axion bursts originating during or after structure
formation $z\lesssim 20$, we find that the details of the cosmological burst
history in this parameter range play a less significant role than the total
converted fraction ${\mathcal{F}}$ and the axion emission spectrum
$dN_{a}/d\omega$, which largely drive the size and detection sensitivity of
signals. This allows to derive strong constraints from non-observation of
photons from the D$a$B that largely do not significantly rely on cosmological
assumptions. In the event of a detection, the spectrum of the D$a$B could
provide detailed information on the cosmological evolution of axions and
properties of the emission sources. When the axion emission bursts primarily
occur within a short period at redshift $z>1$, the flux can exceed the
expected background flux from local DM, providing novel prospects for both
direct and indirect detection that are worthy of further exploration.

Other D$a$B searches offer additional prospects and opportunities to probe
theories beyond SM. In particular, different signatures, such as novel
extragalatic $\gamma$-ray background contributions, can appear from non-photon
axion couplings, motivating further analyses. Direct detection of D$a$B will
manifest distinctly from cold axion DM and is complementary to direct searches
of relativistic axions from individual transient events.

We demonstrated that a broad variety of distinct sources can contribute to
D$a$B, offering rich prospects for studying new physics. Many relativistic
axion emission mechanisms and resulting axion emission flux sources, as well
as their cosmological population distributions and statistics remain poorly
understood. Our work highlights the necessity and serves as an impetus calling
for detailed numerical studies and simulations required to systematically
describe and quantify these scenarios.

## Acknowledgements

We would like to thank Jason Arakawa, John Beacom, Oindrila Ghosh, Shunsaku
Horiuchi, Marianna Safronova, and Muhammad Zaheer for discussions, as well as
Pierluca Carenza and Hyungjin Kim for comments on an early draft. The work of
JE was supported by the World Premier International Research Center Initiative
(WPI), MEXT, Japan and by the JSPS KAKENHI Grant Numbers 21H05451 and
21K20366, as well as by the Swedish Research Council (VR) under grants
2018-03641 and 2019-02337. VT acknowledges support by the World Premier
International Research Center Initiative (WPI), MEXT, Japan and JSPS KAKENHI
grant No. 23K13109. This article is based upon work from COST Action COSMIC
WISPers CA21106, supported by COST (European Cooperation in Science and
Technology). We acknowledge the use of the AxionLimits on GitHub AxionLimits
for compilation of existing constraints.

##  Appendix A Results for different cosmological rate functions

In the main text we analyzed D$a$B detection sensitivity for a power-law rate
distribution (see Eq. (7)) of cosmological transient sources, considering
$z_{\rm max}=20$ and power-law exponent $p=-1$. Here, we show results for
different cosmological assumptions, finding similar features. In Fig. 15 we
display the photon flux associated with magnetic field conversion as well as
decays of axions from D$a$B as in Fig. 7 (upper panel) and Fig. 8 (lower), but
considering sources distributed according to power-law rate function with an
exponent of $p=1$ (blue lines) and $p=-3$ (red). We observe that the change in
exponent modifies distribution of energies, with larger $p$ implying a greater
fraction of the flux lying at lower energies. Fig. 16 illustrates the photon
flux from D$a$B for sources distributed according to Gaussian rate function
with $\bar{z}=5$ (red lines) and $\bar{z}=15$ (black), which is more peaked
than in the power-law case.

Figure 15:  Analogous to Fig. 7 (upper) and Fig. 8 (lower), but considering
D$a$B arises from transient sources with cosmological burst rate of power-law
with exponent $p=1$ (red) and $p=-3$ (blue) in addition to $p=-1$ (black),
$R_{\rm burst}\propto(1+z)^{p}$. For $B$-field conversion we take $R={\rm
kpc}$, and for axion decay we assume a typical source distance of $R(z_{\rm
eff})\simeq 1/H(z_{\rm eff})$ with $z_{\rm eff}$ given by Eq. (35) for
different values of $p$.

Figure 16: Analogous to Fig. 7 (upper) and Fig. 8 (lower), but considering
D$a$B arises from transient sources with cosmological burst rate of gaussian
with width $\delta z=1$ and mean $\bar{z}=5$ (red) and $\bar{z}=15$ (black),
$R_{\rm burst}\propto\exp(-(z-\bar{z})^{2}/\delta z^{2})$. For $B$-field
conversion we take $R={\rm kpc}$, and for axion decay we assume a typical
source distance of $R(z_{\rm eff})\simeq 1/H(z_{\rm eff})$ with $z_{\rm eff}$
given by Eq. (35).

In Fig. 17, we display sensitivity to D$a$B as in Fig. 9 and Fig. 12, but
considering a power-law distribution of transient sources with a different
exponent parameter of $p=1$. While the sensitivity estimates are seen to be
modified, our conclusions regarding the reach of present and future searches
are not significantly affected. We also repeat our analysis using a Gaussian
distribution of axion emission sources as described by Eq. (7), with peak at
$z=\bar{z}$ and spread $\delta z$. The sensitivity estimates for this case are
shown in Fig. 18 and Fig. 19 for $\bar{z}=5$ and $\bar{z}=15$, respectively.
In both cases we consider $\delta z=1$, noting that that our results do not
depend sensitively on this choice. Importantly, the peak of the energy
spectrum is always located at $\bar{\omega}/(1+\bar{z})$, hence the same
experiments can be utilized to search for high-redshift (i.e. $\bar{z}\gg 1$)
sources. The detection of a peaked signal originating at $\bar{z}$ would be
optimized if the search occurs near energies of order
$\bar{\omega}/(1+\bar{z})$.

Figure 17:  Analogous to Fig. 9 (upper) and Fig. 12 (lower), but considering
D$a$B arises from transient sources with cosmological burst rate of power-law
with exponent $p=1$, $R_{\rm burst}\propto(1+z)$. For $B$-field conversion we
take $R={\rm kpc}$, and for axion decay we assume a typical source distance of
$R(z_{\rm eff})\simeq 1/H(z_{\rm eff})$ with $z_{\rm eff}\simeq 2/5$.

Figure 18: Analogous to Fig. 9 (upper) and Fig. 12 (lower), but using a
Gaussian cosmological transient rate function
$f(z)=\exp[-(z-\bar{z})^{2}/\delta z^{2}]$ instead of a power law. For
$B$-field conversion we take $R={\rm kpc}$, and for axion decay we assume a
typical source distance of $R(z_{\rm eff})\simeq 1/H(z_{\rm eff})$ with
$z_{\rm eff}$ given by Eq. (35). Here, we set $\bar{z}=5$ and $\delta z=1$.

Figure 19:  Analogous to Fig. 9 (upper) and Fig. 12 (lower), but using a
Gaussian cosmological transient rate function
$f(z)=\exp[-(z-\bar{z})^{2}/\delta z^{2}]$ rather than a power law. For
$B$-field conversion we take $R={\rm kpc}$, and for axion decay we assume a
typical source distance of $R(z_{\rm eff})\simeq 1/H(z_{\rm eff})$ with
$z_{\rm eff}$ given by Eq. (35). Here, we set $\bar{z}=15$ and $\delta z=1$.

##  Appendix B Photon flux estimates for future experiments

Here, we briefly summarize considerations for estimating differential flux
$d\phi_{\gamma}/d\omega$ sensitivity of future experiments as shown in e.g.
Fig. 10. These results highlight their significant potential to probe a wide
range of D$a$B contributions.

Square Kilometer Array (SKA) Weltman:2018zrl ; Braun:2019gdo ; Ghosh:2020hgd :
the signal and noise flux is estimated as described in Sec. 4.2.

Vera C. Rubin Observatory (Rubin, formerly LSST) Tyson:2002mzo ; Mao:2022fyx :
the expected range of sensitivity to photon energy between $(1-10)\,{\rm eV}$
is shown, without detailed sensitivity estimate.

James Webb Space Telescope (JWST) JWST:2023abc ; Janish:2023kvi ; Roy:2023omw
: the expected JWST background spectrum for Galactic Center observations with
polar angle $\theta_{\rm GC}=10^{\circ}$ and $45^{\circ}$ are shown in
supplemental material of Ref. Roy:2023omw . Fig. S3, depicts that while the
two differ in detail, both lie in the range ${\mathcal{B}}=(0.1-10^{3}){\rm
MJy}=6.2\cdot(10^{-19}-10^{-15})\,{\rm eV/cm}^{2}$ for wavelengths
$\lambda=(0.6-29)\,\mu{\rm m}$, which corresponds to energy range
$\omega=(0.04-2)\,{\rm eV}$. Considering this estimate along with a total
observing time of $t_{\rm JWST}=10^{3}\,{\rm sec}$ (see Tab. 1), we calculate
conversion ${\mathcal{B}}\simeq t_{\rm
JWST}\omega^{2}\,d\phi_{\gamma}/d\omega$ to generate results shown in Fig. 7.

X-Ray Imaging and Spectroscopy Mission (XRISM) XRISMScienceTeam:2020rvx :
simulated spectrum for XRISM observations of the core of M87 (Virgo A) galaxy
can be found in Fig. 7 of Ref. XRISMScienceTeam:2020rvx , which gives
${\mathcal{C}}=(0.1-10)$ counts/keV/s in the energy range $0.3-10\,{\rm keV}$.
The effective collecting area in XRISM is reported to be $>160\,{\rm cm}^{2}$
at $1\,{\rm keV}$ and $>210\,{\rm cm}^{2}$ at $6\,{\rm keV}$
XRISMScienceTeam:2022exn . For our estimates, we consider an area of
$A=200\,{\rm cm}^{2}$. To translate the relevant quantities to units of Fig.
10, we use ${\mathcal{C}}\simeq A\,d\phi_{\gamma}/d\omega$.

All-sky Medium Energy Gamma-ray Observatory (AMEGO) and A space mission for
MeV-GeV gamma-ray astrophysics (e-ASTROGAM) e-ASTROGAM:2017pxr ;
DeAngelis:2021abc ; AMEGO:2021abc ; Engel:2022bgx : These distinct proposals
for MeV-energy gamma-ray search experiments have been estimated to have
similar sensitivity, prompting us to consider them together. Estimate of the
sensitivity of ASTROMEV of ${\mathcal{S}}\simeq(2-4)\cdot 10^{-12}\,{\rm
erg\,cm}^{-2}{\rm s}^{-1}$ in the photon energy range $\omega=(0.1-300)\,{\rm
MeV}$ can be found in Fig. 1 of Ref. DeAngelis:2021abc . We use the relation
${\mathcal{S}}\simeq\omega^{2}d\phi_{\gamma}/d\omega$ to convert this into
flux employed in Fig. 10. A more detailed sensitivity estimate for e-ASTROGAM
is given in Ref. AMEGO:2021abc  (see e.g. their Fig. 8), which agrees with the
above within approximately a factor of two. See also discussion of Ref.
Engel:2022bgx .

## References

  * (1) ATLAS collaboration, G. Aad et al., _Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC_ , _Phys. Lett. B_ 716 (2012) 1 [1207.7214].
  * (2) CMS collaboration, S. Chatrchyan et al., _Observation of a New Boson at a Mass of 125 GeV with the CMS Experiment at the LHC_ , _Phys. Lett. B_ 716 (2012) 30 [1207.7235].
  * (3) R. Peccei and H. R. Quinn, _CP Conservation in the Presence of Instantons_ , _Phys. Rev. Lett._ 38 (1977) 1440.
  * (4) S. Weinberg, _A New Light Boson?_ , _Phys. Rev. Lett._ 40 (1978) 223.
  * (5) F. Wilczek, _Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons_, _Phys. Rev. Lett._ 40 (1978) 279.
  * (6) J. E. Kim, _Weak Interaction Singlet and Strong CP Invariance_ , _Phys. Rev. Lett._ 43 (1979) 103.
  * (7) M. A. Shifman, A. Vainshtein and V. I. Zakharov, _Can Confinement Ensure Natural CP Invariance of Strong Interactions?_ , _Nucl. Phys. B_ 166 (1980) 493.
  * (8) A. Zhitnitsky, _On Possible Suppression of the Axion Hadron Interactions. (In Russian)_ , _Sov. J. Nucl. Phys._ 31 (1980) 260.
  * (9) M. Dine, W. Fischler and M. Srednicki, _A Simple Solution to the Strong CP Problem with a Harmless Axion_ , _Phys. Lett. B_ 104 (1981) 199.
  * (10) J. Preskill, M. B. Wise and F. Wilczek, _Cosmology of the Invisible Axion_ , _Phys. Lett. B_ 120 (1983) 127.
  * (11) L. Abbott and P. Sikivie, _A Cosmological Bound on the Invisible Axion_ , _Phys. Lett. B_ 120 (1983) 133.
  * (12) M. Dine and W. Fischler, _The Not So Harmless Axion_ , _Phys. Lett. B_ 120 (1983) 137.
  * (13) P. Svrcek and E. Witten, _Axions In String Theory_ , _JHEP_ 06 (2006) 051 [hep-th/0605206].
  * (14) A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper and J. March-Russell, _String Axiverse_ , _Phys. Rev. D_ 81 (2010) 123530 [0905.4720].
  * (15) I. G. Irastorza and J. Redondo, _New experimental approaches in the search for axion-like particles_ , _Prog. Part. Nucl. Phys._ 102 (2018) 89 [1801.08127].
  * (16) P. Agrawal et al., _Feebly-interacting particles: FIPs 2020 workshop report_ , _Eur. Phys. J. C_ 81 (2021) 1015 [2102.12143].
  * (17) C. B. Adams et al., _Axion Dark Matter_ , in _Snowmass 2021_ , 3, 2022, 2203.14923.
  * (18) J. P. Conlon and M. C. D. Marsh, _The Cosmophenomenology of Axionic Dark Radiation_ , _JHEP_ 10 (2013) 214 [1304.1804].
  * (19) M. C. D. Marsh, _Hints of a Cosmic Axion Background_ , in _9th Patras Workshop on Axions, WIMPs and WISPs_ , pp. 159–163, 2013, DOI.
  * (20) J. A. Dror, H. Murayama and N. L. Rodd, _Cosmic axion background_ , _Phys. Rev. D_ 103 (2021) 115004 [2101.09287].
  * (21) M. S. Turner, _Thermal Production of Not SO Invisible Axions in the Early Universe_ , _Phys. Rev. Lett._ 59 (1987) 2489.
  * (22) S. Chang and K. Choi, _Hadronic axion window and the big bang nucleosynthesis_ , _Phys. Lett. B_ 316 (1993) 51 [hep-ph/9306216].
  * (23) E. Masso, F. Rota and G. Zsembinszki, _On axion thermalization in the early universe_ , _Phys. Rev. D_ 66 (2002) 023004 [hep-ph/0203221].
  * (24) S. Hannestad, A. Mirizzi and G. Raffelt, _New cosmological mass limit on thermal relic axions_ , _JCAP_ 07 (2005) 002 [hep-ph/0504059].
  * (25) P. Graf and F. D. Steffen, _Thermal axion production in the primordial quark-gluon plasma_ , _Phys. Rev. D_ 83 (2011) 075011 [1008.4528].
  * (26) A. Salvio, A. Strumia and W. Xue, _Thermal axion production_ , _JCAP_ 01 (2014) 011 [1310.6982].
  * (27) 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].
  * (28) E. J. Chun and A. Lukas, _Axino mass in supergravity models_ , _Phys. Lett. B_ 357 (1995) 43 [hep-ph/9503233].
  * (29) T. Asaka and M. Yamaguchi, _Hadronic axion model in gauge mediated supersymmetry breaking and cosmology of saxion_ , _Phys. Rev. D_ 59 (1999) 125003 [hep-ph/9811451].
  * (30) K. Ichikawa, M. Kawasaki, K. Nakayama, M. Senami and F. Takahashi, _Increasing effective number of neutrinos by decaying particles_ , _JCAP_ 05 (2007) 008 [hep-ph/0703034].
  * (31) M. Kawasaki, K. Nakayama and M. Senami, _Cosmological implications of supersymmetric axion models_ , _JCAP_ 03 (2008) 009 [0711.3083].
  * (32) S. Kasuya, M. Kawasaki and T. Yanagida, _Cosmological axion problem in chaotic inflationary universe_ , _Phys. Lett. B_ 409 (1997) 94 [hep-ph/9608405].
  * (33) Y. Ema and K. Nakayama, _Explosive Axion Production from Saxion_ , _Phys. Lett. B_ 776 (2018) 174 [1710.02461].
  * (34) R. T. Co, L. J. Hall and K. Harigaya, _QCD Axion Dark Matter with a Small Decay Constant_ , _Phys. Rev. Lett._ 120 (2018) 211602 [1711.10486].
  * (35) G. G. Raffelt, _Stars as laboratories for fundamental physics: The astrophysics of neutrinos, axions, and other weakly interacting particles_. 5, 1996.
  * (36) P. Sikivie, _Experimental Tests of the Invisible Axion_ , _Phys. Rev. Lett._ 51 (1983) 1415.
  * (37) S. Moriyama, _A Proposal to search for a monochromatic component of solar axions using Fe-57_ , _Phys. Rev. Lett._ 75 (1995) 3222 [hep-ph/9504318].
  * (38) J. Eby, M. Leembruggen, P. Suranyi and L. C. R. Wijewardhana, _Collapse of Axion Stars_ , _JHEP_ 12 (2016) 066 [1608.06911].
  * (39) D. G. Levkov, A. G. Panin and I. I. Tkachev, _Relativistic axions from collapsing Bose stars_ , _Phys. Rev. Lett._ 118 (2017) 011301 [1609.03611].
  * (40) A. Arvanitaki and S. Dubovsky, _Exploring the String Axiverse with Precision Black Hole Physics_ , _Phys. Rev. D_ 83 (2011) 044026 [1004.3558].
  * (41) A. Arvanitaki, M. Baryakhtar and X. Huang, _Discovering the QCD Axion with Black Holes and Gravitational Waves_ , _Phys. Rev. D_ 91 (2015) 084011 [1411.2263].
  * (42) A. Arvanitaki, M. Baryakhtar, S. Dimopoulos, S. Dubovsky and R. Lasenby, _Black Hole Mergers and the QCD Axion at Advanced LIGO_ , _Phys. Rev. D_ 95 (2017) 043001 [1604.03958].
  * (43) M. Baryakhtar, M. Galanis, R. Lasenby and O. Simon, _Black hole superradiance of self-interacting scalar fields_ , _Phys. Rev. D_ 103 (2021) 095019 [2011.11646].
  * (44) C. Ünal, F. Pacucci and A. Loeb, _Properties of ultralight bosons from heavy quasar spins via superradiance_ , _JCAP_ 05 (2021) 007 [2012.12790].
  * (45) N. P. Branco, R. Z. Ferreira and J. a. G. Rosa, _Superradiant axion clouds around asteroid-mass primordial black holes_ , _JCAP_ 04 (2023) 003 [2301.01780].
  * (46) A. Banerjee, D. Budker, J. Eby, H. Kim and G. Perez, _Relaxion Stars and their detection via Atomic Physics_ , _Commun. Phys._ 3 (2020) 1 [1902.08212].
  * (47) A. Banerjee, D. Budker, J. Eby, V. V. Flambaum, H. Kim, O. Matsedonskyi et al., _Searching for Earth/Solar Axion Halos_ , _JHEP_ 09 (2020) 004 [1912.04295].
  * (48) D. Budker, J. Eby, M. Gorghetto, M. Jiang and G. Perez, _A Generic Formation Mechanism of Ultralight Dark Matter Solar Halos_ , 2306.12477.
  * (49) J. Eby, S. Shirai, Y. V. Stadnik and V. Takhistov, _Probing relativistic axions from transient astrophysical sources_ , _Phys. Lett. B_ 825 (2022) 136858 [2106.14893].
  * (50) S. Chigusa, A. Ito, K. Nakayama and V. Takhistov, _Effects of Finite Material Size On Axion-magnon Conversion_ , 2310.17704.
  * (51) J. Arakawa, J. Eby, M. S. Safronova, V. Takhistov and M. H. Zaheer, _Detection of Bosenovae with Quantum Sensors on Earth and in Space_ , 2306.16468.
  * (52) C. Dailey, C. Bradley, D. F. Jackson Kimball, I. A. Sulai, S. Pustelny, A. Wickenbrock et al., _Quantum sensor networks as exotic field telescopes for multi-messenger astronomy_ , _Nature Astron._ 5 (2021) 150 [2002.04352].
  * (53) G. G. Raffelt, J. Redondo and N. Viaux Maira, _The meV mass frontier of axion physics_ , _Phys. Rev. D_ 84 (2011) 103008 [1110.6397].
  * (54) F. Calore, P. Carenza, M. Giannotti, J. Jaeckel and A. Mirizzi, _Bounds on axionlike particles from the diffuse supernova flux_ , _Phys. Rev. D_ 102 (2020) 123005 [2008.11741].
  * (55) F. Calore, P. Carenza, C. Eckner, T. Fischer, M. Giannotti, J. Jaeckel et al., _3D template-based Fermi-LAT constraints on the diffuse supernova axion-like particle background_ , _Phys. Rev. D_ 105 (2022) 063028 [2110.03679].
  * (56) V. Takhistov, _Neutrino and Axion Astronomy with Dark Matter Experiments_ , _J. Phys. Conf. Ser._ 2156 (2021) 012049 [2111.04088].
  * (57) J. F. Beacom, _The Diffuse Supernova Neutrino Background_ , _Ann. Rev. Nucl. Part. Sci._ 60 (2010) 439 [1004.3311].
  * (58) S. Horiuchi, J. F. Beacom and E. Dwek, _The Diffuse Supernova Neutrino Background is detectable in Super-Kamiokande_ , _Phys. Rev. D_ 79 (2009) 083013 [0812.3157].
  * (59) Super-Kamiokande collaboration, K. Abe et al., _Diffuse supernova neutrino background search at Super-Kamiokande_ , _Phys. Rev. D_ 104 (2021) 122002 [2109.11174].
  * (60) X.-D. Shi and G. M. Fuller, _Neutrinos and supermassive stars: Prospects for neutrino emission and detection_ , _Astrophys. J._ 503 (1998) 307 [astro-ph/9801106].
  * (61) F. Linke, J. A. Font, H.-T. Janka, E. Muller and P. Papadopoulos, _Spherical collapse of supermassive stars: Neutrino emission and gamma-ray bursts_ , _Astron. Astrophys._ 376 (2001) 568 [astro-ph/0103144].
  * (62) V. Munoz, V. Takhistov, S. J. Witte and G. M. Fuller, _Exploring the origin of supermassive black holes with coherent neutrino scattering_ , _JCAP_ 11 (2021) 020 [2102.00885].
  * (63) S. Ando, N. Ekanger, S. Horiuchi and Y. Koshio, _Diffuse neutrino background from past core-collapse supernovae_ , 6, 2023, 2306.16076.
  * (64) X.-D. Shi, G. M. Fuller and F. Halzen, _Observing the birth of supermassive black holes with the ’ICECUBE’ neutrino detector_ , _Phys. Rev. Lett._ 81 (1998) 5722 [astro-ph/9805242].
  * (65) Planck collaboration, N. Aghanim et al., _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6 [1807.06209].
  * (66) 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].
  * (67) P. Carenza, B. Fore, M. Giannotti, A. Mirizzi and S. Reddy, _Enhanced Supernova Axion Emission and its Implications_ , _Phys. Rev. Lett._ 126 (2021) 071102 [2010.02943].
  * (68) A. Lella, P. Carenza, G. Lucente, M. Giannotti and A. Mirizzi, _Protoneutron stars as cosmic factories for massive axionlike particles_ , _Phys. Rev. D_ 107 (2023) 103017 [2211.13760].
  * (69) S. P. Harris, J.-F. Fortin, K. Sinha and M. G. Alford, _Axions in neutron star mergers_ , _JCAP_ 07 (2020) 023 [2003.09768].
  * (70) D. F. G. Fiorillo and F. Iocco, _Axions from neutron star mergers_ , _Phys. Rev. D_ 105 (2022) 123007.
  * (71) G. G. Raffelt, _Astrophysical methods to constrain axions and other novel particle phenomena_ , _Phys. Rept._ 198 (1990) 1.
  * (72) G. G. Raffelt, _Astrophysical axion bounds_ , _Lect. Notes Phys._ 741 (2008) 51 [hep-ph/0611350].
  * (73) T. Fischer, S. Chakraborty, M. Giannotti, A. Mirizzi, A. Payez and A. Ringwald, _Probing axions with the neutrino signal from the next galactic supernova_ , _Phys. Rev. D_ 94 (2016) 085012 [1605.08780].
  * (74) Kamiokande-II collaboration, K. Hirata et al., _Observation of a Neutrino Burst from the Supernova SN 1987a_ , _Phys. Rev. Lett._ 58 (1987) 1490.
  * (75) R. M. Bionta et al., _Observation of a Neutrino Burst in Coincidence with Supernova SN 1987a in the Large Magellanic Cloud_ , _Phys. Rev. Lett._ 58 (1987) 1494.
  * (76) E. N. Alekseev, L. N. Alekseeva, I. V. Krivosheina and V. I. Volchenko, _Detection of the Neutrino Signal From SN1987A in the LMC Using the Inr Baksan Underground Scintillation Telescope_ , _Phys. Lett. B_ 205 (1988) 209.
  * (77) M. S. Turner, _Axions from SN 1987a_ , _Phys. Rev. Lett._ 60 (1988) 1797.
  * (78) A. Burrows, M. S. Turner and R. P. Brinkmann, _Axions and SN 1987a_ , _Phys. Rev. D_ 39 (1989) 1020.
  * (79) A. Burrows, M. T. Ressell and M. S. Turner, _Axions and SN1987A: Axion trapping_ , _Phys. Rev. D_ 42 (1990) 3297.
  * (80) G. Raffelt and D. Seckel, _Bounds on Exotic Particle Interactions from SN 1987a_ , _Phys. Rev. Lett._ 60 (1988) 1793.
  * (81) W. Keil, H.-T. Janka, D. N. Schramm, G. Sigl, M. S. Turner and J. R. Ellis, _A Fresh look at axions and SN-1987A_ , _Phys. Rev. D_ 56 (1997) 2419 [astro-ph/9612222].
  * (82) C. Hanhart, D. R. Phillips and S. Reddy, _Neutrino and axion emissivities of neutron stars from nucleon-nucleon scattering data_ , _Phys. Lett. B_ 499 (2001) 9 [astro-ph/0003445].
  * (83) J. H. Chang, R. Essig and S. D. McDermott, _Supernova 1987A Constraints on Sub-GeV Dark Sectors, Millicharged Particles, the QCD Axion, and an Axion-like Particle_ , _JHEP_ 09 (2018) 051 [1803.00993].
  * (84) G. G. Raffelt, _ASTROPHYSICAL AXION BOUNDS DIMINISHED BY SCREENING EFFECTS_ , _Phys. Rev. D_ 33 (1986) 897.
  * (85) J. Jaeckel and A. Ringwald, _The Low-Energy Frontier of Particle Physics_ , _Ann. Rev. Nucl. Part. Sci._ 60 (2010) 405 [1002.0329].
  * (86) A. Lella, P. Carenza, G. Co’, G. Lucente, M. Giannotti, A. Mirizzi et al., _Getting the most on supernova axions_ , _Phys. Rev. D_ 109 (2024) 023001 [2306.01048].
  * (87) J. A. Grifols, E. Masso and R. Toldra, _Gamma-rays from SN1987A due to pseudoscalar conversion_ , _Phys. Rev. Lett._ 77 (1996) 2372 [astro-ph/9606028].
  * (88) J. W. Brockway, E. D. Carlson and G. G. Raffelt, _SN1987A gamma-ray limits on the conversion of pseudoscalars_ , _Phys. Lett. B_ 383 (1996) 439 [astro-ph/9605197].
  * (89) A. Payez, C. Evoli, T. Fischer, M. Giannotti, A. Mirizzi and A. Ringwald, _Revisiting the SN1987A gamma-ray limit on ultralight axion-like particles_ , _JCAP_ 02 (2015) 006 [1410.3747].
  * (90) J. Jaeckel, P. C. Malta and J. Redondo, _Decay photons from the axionlike particles burst of type II supernovae_ , _Phys. Rev. D_ 98 (2018) 055032 [1702.02964].
  * (91) D. J. Kaup, _Klein-Gordon Geon_ , _Phys. Rev._ 172 (1968) 1331.
  * (92) R. Ruffini and S. Bonazzola, _Systems of selfgravitating particles in general relativity and the concept of an equation of state_ , _Phys. Rev._ 187 (1969) 1767.
  * (93) J. Breit, S. Gupta and A. Zaks, _Cold bose stars_ , _Physics Letters B_ 140 (1984) 329 .
  * (94) M. Colpi, S. Shapiro and I. Wasserman, _Boson Stars: Gravitational Equilibria of Selfinteracting Scalar Fields_ , _Phys. Rev. Lett._ 57 (1986) 2485.
  * (95) E. W. Kolb and I. I. Tkachev, _Axion miniclusters and Bose stars_ , _Phys. Rev. Lett._ 71 (1993) 3051 [hep-ph/9303313].
  * (96) B. Eggemeier and J. C. Niemeyer, _Formation and mass growth of axion stars in axion miniclusters_ , _Phys. Rev. D_ 100 (2019) 063528 [1906.01348].
  * (97) P.-H. Chavanis, _Mass-radius relation of Newtonian self-gravitating Bose-Einstein condensates with short-range interactions: I. Analytical results_ , _Phys. Rev. D_ 84 (2011) 043531 [1103.2050].
  * (98) P. Chavanis and L. Delfini, _Mass-radius relation of Newtonian self-gravitating Bose-Einstein condensates with short-range interactions: II. Numerical results_ , _Phys. Rev. D_ 84 (2011) 043532 [1103.2054].
  * (99) J. Barranco and A. Bernal, _Self-gravitating system made of axions_ , _Phys. Rev. D_ 83 (2011) 043525 [1001.1769].
  * (100) J. Eby, P. Suranyi, C. Vaz and L. Wijewardhana, _Axion Stars in the Infrared Limit_ , _JHEP_ 03 (2015) 080 [1412.3430].
  * (101) E. D. Schiappacasse and M. P. Hertzberg, _Analysis of Dark Matter Axion Clumps with Spherical Symmetry_ , _JCAP_ 01 (2018) 037 [1710.04729].
  * (102) E. A. Donley, N. R. Claussen, S. L. Cornish, J. L. Roberts, E. A. Cornell and C. E. Wieman, _Dynamics of collapsing and exploding bose–einstein condensates_ , _Nature_ 412 (2001) 295.
  * (103) J. Eby, P. Suranyi and L. Wijewardhana, _The Lifetime of Axion Stars_ , _Mod. Phys. Lett. A_ 31 (2016) 1650090 [1512.01709].
  * (104) J. Eby, M. Ma, P. Suranyi and L. Wijewardhana, _Decay of Ultralight Axion Condensates_ , _JHEP_ 01 (2018) 066 [1705.05385].
  * (105) X. Du, D. J. E. Marsh, M. Escudero, A. Benson, D. Blas, C. K. Pooni et al., _Soliton Merger Rates and Enhanced Axion Dark Matter Decay_ , 2301.09769.
  * (106) M. Escudero, C. K. Pooni, M. Fairbairn, D. Blas, X. Du and D. J. E. Marsh, _Axion Star Explosions: A New Source for Axion Indirect Detection_ , 2302.10206.
  * (107) P. J. Fox, N. Weiner and H. Xiao, _Recurrent Axinovae and their Cosmological Constraints_ , 2302.00685.
  * (108) T. Altherr, E. Petitgirard and T. del Rio Gaztelurrutia, _Axion emission from red giants and white dwarfs_ , _Astropart. Phys._ 2 (1994) 175 [hep-ph/9310304].
  * (109) G. Raffelt and A. Weiss, _Red giant bound on the axion - electron coupling revisited_ , _Phys. Rev. D_ 51 (1995) 1495 [hep-ph/9410205].
  * (110) A. H. Corsico, O. G. Benvenuto, L. G. Althaus, J. Isern and E. Garcia-Berro, _The Potential of the variable DA white dwarf G117 - B15A as a tool for fundamental physics_ , _New Astron._ 6 (2001) 197 [astro-ph/0104103].
  * (111) 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].
  * (112) N. Iwamoto, _Axion Emission from Neutron Stars_ , _Phys. Rev. Lett._ 53 (1984) 1198.
  * (113) J. Keller and A. Sedrakian, _Axions from cooling compact stars_ , _Nucl. Phys. A_ 897 (2013) 62 [1205.6940].
  * (114) A. Sedrakian, _Axion cooling of neutron stars_ , _Phys. Rev. D_ 93 (2016) 065044 [1512.07828].
  * (115) T. Dietrich and K. Clough, _Cooling binary neutron star remnants via nucleon-nucleon-axion bremsstrahlung_ , _Phys. Rev. D_ 100 (2019) 083005 [1909.01278].
  * (116) A. Prabhu, _Axion production in pulsar magnetosphere gaps_ , _Phys. Rev. D_ 104 (2021) 055038 [2104.14569].
  * (117) D. Noordhuis, A. Prabhu, S. J. Witte, A. Y. Chen, F. Cruz and C. Weniger, _Novel Constraints on Axions Produced in Pulsar Polar-Cap Cascades_ , _Phys. Rev. Lett._ 131 (2023) 111004 [2209.09917].
  * (118) D. Noordhuis, A. Prabhu, C. Weniger and S. J. Witte, _Axion Clouds around Neutron Stars_ , 2307.11811.
  * (119) W. H. Press and S. A. Teukolsky, _Floating Orbits, Superradiant Scattering and the Black-hole Bomb_ , _Nature_ 238 (1972) 211.
  * (120) R. Brito, V. Cardoso and P. Pani, _Superradiance: New Frontiers in Black Hole Physics_ , _Lect. Notes Phys._ 906 (2015) pp.1 [1501.06570].
  * (121) P. B. Ferraz, T. W. Kephart and J. a. G. Rosa, _Superradiant pion clouds around primordial black holes_ , _JCAP_ 07 (2022) 026 [2004.11303].
  * (122) N. Bernal, Y. F. Perez-Gonzalez and Y. Xu, _Superradiant production of heavy dark matter from primordial black holes_ , _Phys. Rev. D_ 106 (2022) 015020 [2205.11522].
  * (123) D. Baumann, H. S. Chia and R. A. Porto, _Probing Ultralight Bosons with Binary Black Holes_ , _Phys. Rev. D_ 99 (2019) 044001 [1804.03208].
  * (124) D. Levkov, A. Panin and I. Tkachev, _Gravitational Bose-Einstein condensation in the kinetic regime_ , _Phys. Rev. Lett._ 121 (2018) 151301 [1804.05857].
  * (125) J. Chen, X. Du, E. W. Lentz, D. J. E. Marsh and J. C. Niemeyer, _New insights into the formation and growth of boson stars in dark matter halos_ , _Phys. Rev. D_ 104 (2021) 083022 [2011.01333].
  * (126) K. Kirkpatrick, A. E. Mirasola and C. Prescod-Weinstein, _Relaxation times for Bose-Einstein condensation in axion miniclusters_ , _Phys. Rev. D_ 102 (2020) 103012 [2007.07438].
  * (127) J. Chen, X. Du, E. W. Lentz and D. J. E. Marsh, _Relaxation times for Bose-Einstein condensation by self-interaction and gravity_ , _Phys. Rev. D_ 106 (2022) 023009 [2109.11474].
  * (128) K. Kirkpatrick, A. E. Mirasola and C. Prescod-Weinstein, _Analysis of Bose-Einstein condensation times for self-interacting scalar dark matter_ , _Phys. Rev. D_ 106 (2022) 043512 [2110.08921].
  * (129) A. S. Dmitriev, D. G. Levkov, A. G. Panin and I. I. Tkachev, _Self-similar growth of Bose stars_ , 2305.01005.
  * (130) M. Jain, W. Wanichwecharungruang and J. Thomas, _Kinetic relaxation and nucleation of Bose stars in self-interacting wave dark matter_ , 2310.00058.
  * (131) G. B. Gelmini, _The Hunt for Dark Matter_ , in _Theoretical Advanced Study Institute in Elementary Particle Physics: Journeys Through the Precision Frontier: Amplitudes for Colliders_ , pp. 559–616, 2015, 1502.01320, DOI.
  * (132) A. Maselli, P. Pnigouras, N. G. Nielsen, C. Kouvaris and K. D. Kokkotas, _Dark stars: gravitational and electromagnetic observables_ , _Phys. Rev. D_ 96 (2017) 023005 [1704.07286].
  * (133) B. B. Kamenetskaia, A. Brenner, A. Ibarra and C. Kouvaris, _Proton capture in compact dark stars and observable implications_ , _JCAP_ 10 (2023) 027 [2211.05845].
  * (134) D. Curtin and J. Setford, _How To Discover Mirror Stars_ , _Phys. Lett. B_ 804 (2020) 135391 [1909.04071].
  * (135) M. Hippert, J. Setford, H. Tan, D. Curtin, J. Noronha-Hostler and N. Yunes, _Mirror neutron stars_ , _Phys. Rev. D_ 106 (2022) 035025 [2103.01965].
  * (136) M. Hippert, H. Tan, J. Noronha-Hostler, N. Yunes, J. Setford and D. Curtin, _Mirror Neutron Stars: How QCD Can Be Used to Study Dark Matter Through Gravitational Waves_ , _Acta Phys. Polon. Supp._ 16 (2023) 1 [2207.13063].
  * (137) A. Mirizzi, G. G. Raffelt and P. D. Serpico, _Photon-axion conversion in intergalactic magnetic fields and cosmological consequences_ , _Lect. Notes Phys._ 741 (2008) 115 [astro-ph/0607415].
  * (138) R. Beck, _Galactic and Extragalactic Magnetic Fields_ , _AIP Conf. Proc._ 1085 (2009) 83 [0810.2923].
  * (139) G. Raffelt and L. Stodolsky, _Mixing of the Photon with Low Mass Particles_ , _Phys. Rev. D_ 37 (1988) 1237.
  * (140) A. Cooray, _Extragalactic Background Light: Measurements and Applications_ , 1602.03512.
  * (141) NuSTAR collaboration, F. A. Harrison et al., _The Nuclear Spectroscopic Telescope Array (NuSTAR) High-Energy X-Ray Mission_ , _Astrophys. J._ 770 (2013) 103 [1301.7307].
  * (142) 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].
  * (143) J. W. Foster, M. Kongsore, C. Dessert, Y. Park, N. L. Rodd, K. Cranmer et al., _Deep Search for Decaying Dark Matter with XMM-Newton Blank-Sky Observations_ , _Phys. Rev. Lett._ 127 (2021) 051101 [2102.02207].
  * (144) L. Bouchet, E. Jourdain, J. P. Roques, A. Strong, R. Diehl, F. Lebrun et al., _INTEGRAL SPI All-Sky View in Soft Gamma Rays: Study of Point Source and Galactic Diffuse Emissions_ , _Astrophys. J._ 679 (2008) 1315 [0801.2086].
  * (145) J. G. Stacy, W. Collmar, A. Strong, V. Schonfelder and A. Carraminana, _Limits on MeV Gamma-Ray Emission from Active Galaxies and Other Unidentified High-Latitude Gamma-Ray Sources Observed with COMPTEL_ , in _International Cosmic Ray Conference_ , vol. 3 of _International Cosmic Ray Conference_ , pp. 1085–1088, Jan., 2008.
  * (146) A. Strong and W. Collmar, _COMPTEL Reloaded: a heritage progect in MeV astronomy_ , _Mem. Soc. Ast. It._ 90 (2019) 297 [1907.07454].
  * (147) A. W. Strong, I. V. Moskalenko and O. Reimer, _Diffuse galactic continuum gamma rays. A Model compatible with EGRET data and cosmic-ray measurements_ , _Astrophys. J._ 613 (2004) 962 [astro-ph/0406254].
  * (148) Fermi-LAT collaboration, M. Ackermann et al., _The spectrum of isotropic diffuse gamma-ray emission between 100 MeV and 820 GeV_ , _Astrophys. J._ 799 (2015) 86 [1410.3696].
  * (149) Fermi-LAT collaboration, G. Principe, _The Fermi Gamma Ray Sky: summary of recent Observations_ , 9, 2022, 2209.03652.
  * (150) LSST collaboration, J. A. Tyson, _Large synoptic survey telescope: Overview_ , _Proc. SPIE Int. Soc. Opt. Eng._ 4836 (2002) 10 [astro-ph/0302102].
  * (151) Y.-Y. Mao et al., _Snowmass2021: Vera C. Rubin Observatory as a Flagship Dark Matter Experiment_ , 2203.07252.
  * (152) J. P. Gardner et al., _The James Webb Space Telescope Mission_ , _Publ. Astron. Soc. Pac_ 135 (2023) 068001 [2304.04869].
  * (153) R. Janish and E. Pinetti, _Hunting Dark Matter Lines in the Infrared Background with the James Webb Space Telescope_ , 2310.15395.
  * (154) S. Roy, C. Blanco, C. Dessert, A. Prabhu and T. Temim, _Sensitivity of JWST to eV-Scale Decaying Axion Dark Matter_ , 2311.04987.
  * (155) XRISM Science Team collaboration, _Science with the X-ray Imaging and Spectroscopy Mission (XRISM)_ , 2003.04962.
  * (156) A. Weltman et al., _Fundamental physics with the Square Kilometre Array_ , _Publ. Astron. Soc. Austral._ 37 (2020) e002 [1810.02680].
  * (157) R. Braun, A. Bonaldi, T. Bourke, E. Keane and J. Wagg, _Anticipated Performance of the Square Kilometre Array – Phase 1 (SKA1)_ , 1912.12699.
  * (158) O. Ghosh, J. Salvado and J. Miralda-Escudé, _Axion Gegenschein: Probing Back-scattering of Astrophysical Radio Sources Induced by Dark Matter_ , 2008.02729.
  * (159) D. Cadamuro and J. Redondo, _Cosmological bounds on pseudo Nambu-Goldstone bosons_ , _JCAP_ 02 (2012) 032 [1110.2895].
  * (160) J. W. Foster, S. J. Witte, M. Lawson, T. Linden, V. Gajjar, C. Weniger et al., _Extraterrestrial Axion Search with the Breakthrough Listen Galactic Center Survey_ , _Phys. Rev. Lett._ 129 (2022) 251102 [2202.08274].
  * (161) ADMX collaboration, N. Du et al., _A Search for Invisible Axion Dark Matter with the Axion Dark Matter Experiment_ , _Phys. Rev. Lett._ 120 (2018) 151301 [1804.05750].
  * (162) ADMX collaboration, C. Boutan et al., _Piezoelectrically Tuned Multimode Cavity Search for Axion Dark Matter_ , _Phys. Rev. Lett._ 121 (2018) 261302 [1901.00920].
  * (163) ADMX collaboration, T. Braine et al., _Extended Search for the Invisible Axion with the Axion Dark Matter Experiment_ , _Phys. Rev. Lett._ 124 (2020) 101303 [1910.08638].
  * (164) N. Crisosto, P. Sikivie, N. S. Sullivan, D. B. Tanner, J. Yang and G. Rybka, _ADMX SLIC: Results from a Superconducting $LC$ Circuit Investigating Cold Axions_, _Phys. Rev. Lett._ 124 (2020) 241101 [1911.05772].
  * (165) ADMX collaboration, C. Bartram et al., _Dark matter axion search using a Josephson Traveling wave parametric amplifier_ , _Rev. Sci. Instrum._ 94 (2023) 044703 [2110.10262].
  * (166) ADMX collaboration, C. Bartram et al., _Search for Invisible Axion Dark Matter in the 3.3–4.2 $\mu$eV Mass Range_, _Phys. Rev. Lett._ 127 (2021) 261803 [2110.06096].
  * (167) S. Lee, S. Ahn, J. Choi, B. R. Ko and Y. K. Semertzidis, _Axion Dark Matter Search around 6.7 $\mu$eV_, _Phys. Rev. Lett._ 124 (2020) 101802 [2001.05102].
  * (168) J. Jeong, S. Youn, S. Bae, J. Kim, T. Seong, J. E. Kim et al., _Search for Invisible Axion Dark Matter with a Multiple-Cell Haloscope_ , _Phys. Rev. Lett._ 125 (2020) 221302 [2008.10141].
  * (169) CAPP collaboration, O. Kwon et al., _First Results from an Axion Haloscope at CAPP around 10.7 $\mu$eV_, _Phys. Rev. Lett._ 126 (2021) 191802 [2012.10764].
  * (170) Y. Lee, B. Yang, H. Yoon, M. Ahn, H. Park, B. Min et al., _Searching for Invisible Axion Dark Matter with an 18 T Magnet Haloscope_ , _Phys. Rev. Lett._ 128 (2022) 241805 [2206.08845].
  * (171) J. Kim et al., _Near-Quantum-Noise Axion Dark Matter Search at CAPP around 9.5 $\mu$eV_, _Phys. Rev. Lett._ 130 (2023) 091602 [2207.13597].
  * (172) A. K. Yi et al., _Axion Dark Matter Search around 4.55 $\mu$eV with Dine-Fischler-Srednicki-Zhitnitskii Sensitivity_, _Phys. Rev. Lett._ 130 (2023) 071002 [2210.10961].
  * (173) B. Yang, H. Yoon, M. Ahn, Y. Lee and J. Yoo, _Extended Axion Dark Matter Search Using the CAPP18T Haloscope_ , _Phys. Rev. Lett._ 131 (2023) 081801 [2308.09077].
  * (174) C. M. Adair et al., _Search for Dark Matter Axions with CAST-CAPP_ , _Nature Commun._ 13 (2022) 6180 [2211.02902].
  * (175) CAST collaboration, A. A. Melcón et al., _First results of the CAST-RADES haloscope search for axions at 34.67 $\mu$eV_, _JHEP_ 21 (2020) 075 [2104.13798].