Probing the Blue Axion with Cosmic Optical Background Anisotropies

A radiative decaying Big Bang relic with a mass $m_a\simeq 5-25 \,\rm eV$, which we dub"blue axion", can be probed with direct and indirect observations of the cosmic optical background (COB). The strongest bounds on blue-axion cold dark matter come from the Hubble Space Telescope (HST) measurements of COB anisotropies at $606$~nm. We suggest that new HST measurements at higher frequencies ($336$~nm and $438$~nm) can improve current constraints on the lifetime up to one order of magnitude, and we show that also thermally produced and hot relic blue axions can be competitively probed by COB anisotropies. We exclude the simple interpretation of the excess in the diffuse COB detected by the Long Range Reconnaissance Imager (LORRI) as photons produced by a decaying hot relic. Finally, we comment on the reach of upcoming line intensity mapping experiments, that could detect blue axions with a lifetime as large as $10^{29}\,\rm s$ or $10^{27}\,\rm s$ for the cold dark matter and the hot relic case, respectively.


I. INTRODUCTION
The extragalactic background light (EBL) is the accumulated radiation emitted in the Universe by all galaxies and active galactic nuclei over the cosmic history, ranging from far-infrared to ultraviolet bands [1][2][3][4][5].A robust lower bound on the cosmic optical background (COB) is obtained through the Hubble Space Telescope (HST) galaxy counts [6,7].However, the shape and intensity of the COB spectrum contributions from diffuse, unresolved sources need to be obtained with different strategiese.g. by directly detecting the COB with a spacecraft, so to avoid the airglow and artificial light affecting ground-based measurements.The Long Range Reconnaissance Imager (LORRI) instrument on NASA's New Horizons mission has reported the COB to be 16.37 ± 1.47 nW m −2 sr −1 at the pivot wavelength of 0.608 µm [8][9][10].This is about ∼ 4σ above the HST galaxy count estimate, suggesting the existence of an unaccounted for EBL component of 8.06 ± 1.92 nW m −2 sr −1 .Such excess could be due to e.g.high redshift galaxies [11] or directcollapse black holes [12].A more mundane explanation of the excess could be related to the modeling of the zodiacal light (ZL), sunlight scattered by interplanetary dust.Zodiacal light largely dominates the sky brightness in the inner solar system.While it can be neglected at the distances from the Sun (51.3 A.U.) where LORRI observations were obtained, a careful estimate of ZL is necessary to evaluate the contribution of diffuse galactic light to the total sky [10].
Since γ-rays (around 100 GeV) are absorbed through electron-positron pair production when scattering on O(10 eV) COB photons, analyses of the observed blazar spectra allow for an indirect measurement of the COB and its redshift evolution [2][3][4][5].This approach has other systematic uncertainties that hinder an easy evaluation of the COB, most notably the knowledge of the injected blazar spectrum, and the possible production of secondary γ-rays [2,[13][14][15].Thus, several effects can potentially undermine the indirect estimates, which are found to be comparable to the COB inferred from galaxy counts [16,17], and in tension with direct measurements.(See however Ref. [18] for an estimate with larger uncertainties.) An alternative to the above mentioned approaches relies on measuring the anisotropies rather than the diffuse intensity of the COB.The foregrounds such as ZL have smooth spatial distributions and a different correlation function compared to the fluctuations generated by hypothetical extragalactic signals, therefore tackling one of the major drawbacks of direct measurements [19][20][21].
Big Bang relics a decaying to photons through processes such as a → γ + γ and a → χ + γ can potentially contribute to the EBL, and crucially to the COB if they have a mass m a ≃ 5 − 25 eV.The authors of Ref. [22] have shown that the excess in the diffuse COB detected by LORRI can be produced by dark matter (DM) in the form of axion-like particles decaying to blue and ultraviolet light, hereafter dubbed "blue axions".However, HST anisotropy measurements at 606 nm already exclude the hypothesis that the excess is due to cold blue axions [23].In the following, we update the cold dark matter (CDM) anisotropy bound from 606 nm measurements with an improved analysis, which includes two additional anisotropy sources (shot noise and foregrounds) and a refined treatment of the detector response, finding good agreement with Ref. [23].We will show that new dedicated HST measurements at wavelengths 336 nm and 438 nm can constrain the lifetime of the blue axion by a further factor 4 to 10, resulting in the best probe to date, as can be seen in Fig. 1.Additional shorter wavelength anisotropy measurements will be targets for forthcoming sub-orbital and space-based measurements [24].
The angular power spectrum depends on the abundance of the relic ρ a and its power spectrum P δ , therefore the same formalism can be applied to different production mechanism scenarios once ρ a and P δ are evaluated for each case.Besides the cold blue axions produced through misalignment mechanism, we will show that 336 nm and 438 nm anisotropy measurements can competitively probe blue axions produced from the annihilation of string-wall networks arising from a symmetry breaking pattern.In this case, they can constitute 100% of the DM but could have a cutoff in the power spectrum, similar to warm dark matter.Blue axions could also be a non-cold dark matter (NCDM) component.Such a component could have either the temperature equal to neutrino temperature ("hot relic") or a temperature set by a freeze-out mechanism in cosmological scenarios.In the latter case, additional degrees of freedom above the Electroweak scale are required to satisfy structure formation bounds on NCDM abundance [29].Nevertheless, we will show that even a hotter component is excluded by anisotropy measurements, further constraining the interpretation of LORRI excess as a decaying Big Bang relic.
We also tackle the problem of alternative particlephysics inspired interpretations of an excess in the diffuse COB.Any such interpretation needs to avoid anisotropy constraints.Moreover, if blue axions decay to two photons (a → γ + γ), strong bounds arise from globular cluster observations, since they would be copiously produced in horizontal branch (HB) stars through Primakoff production [25].An intuitive way to circumvent both anisotropy and stellar cooling bounds is granted by a small hot dark matter component decaying through a dark portal featuring an additional dark vector, a → χ + γ [30].In this case, the power spectrum is suppressed at small scales, and the stellar cooling bound is given by the required agreement between the predicted and observationally inferred core mass at the helium flash of red giants (RGs), since dark sector particles can be produced by plasmon decay γ * → χ + a.
Finally, we explore the possibilities offered by line intensity mapping (LIM), an emerging tool for cosmology with the goal of measuring the integrated emission along the line of sight from spectral lines emitted in the past.A decaying relic would show up as an "interloper line", and applying simple scaling arguments to recent results in the literature we show that LIM can probe blue axion hot relics with a lifetime as large as 10 27 s.
We begin in Sec.II with a derivation of the general anisotropy power spectrum due to a decaying Big Bang Relic.Then in Sec.III we present the possible production mechanisms of the blue axion, and the corresponding bounds and reaches from anisotropy measurements are discussed in Sec.IV.In Sec.V we explore the dark portal scenario.Sec.VI is dedicated to the projected reach of LIM experiments.Finally, Sec.VII is devoted to a summary and discussion.

II. ISOTROPIC AND ANISOTROPIC COSMIC OPTICAL BACKGROUND
We assume that a population of blue axions exists, and for the time being we will be agnostic about its production mechanism and abundance.Blue axions are assumed to have a coupling to two photons, where Because of this interaction Lagrangian, blue axions decay to two photons with a decay rate Following Ref. [30] to compute the contribution to the COB from the decaying blue axions, we define the average energy intensity in units of energy per time per surface per steradians where we also introduced the window function W .If one assumes that the relic decays at rest, the decay spectrum of photons is monochromatic with energy ω max = m a /2 and the energy gets redshifted.The intensity is obtained integrating over the cosmic time and accounting for ζ = 2 photons produced in each decay [30][31][32], where z = ωmax ω − 1.
The Hubble parameter at the redshift z is given by H 4 where H 0 = h 100 km s −1 Mpc −1 is the value of the Hubble parameter today, h = 0.674, and Ω Λ = 0.685, Ω m = 0.315 and Ω r = 5.38×10 −5 (neglecting the contribution of massless neutrinos) denote the density parameter of the dark energy, total matter and radiation, respectively, taken from Ref. [33].Comparing Eq. (3) to Eq. ( 4), the explicit form of the window function W is obtained.We assume that the depletion in the blue axion population due to decay is negligible, and that the decay happens when blue axions are non-relativistic.We see that, to explain the excess detected by LORRI, assuming blue axions to constitute 100% of the dark matter and ω ≃ m a , the rate needs to be Γ a→γγ ≃ 10 −23 − 10 −22 s −1 .
Since blue axions can be clumped, their decay can show up anisotropically in the sky.To discuss the anisotropy, we first need to define the energy intensity as seen in the FIG. 1. Bounds and projected reaches on the lifetime of the blue-axion cold dark matter.We recompute the 606 nm bound (solid blue lines), and find good agreement with the results of Ref. [23].The forecasts of our proposed observations at 438 nm and 336 nm are given with dashed blue lines.Other constraints are also shown: HB stars [25] (green), VIMOS [26] (yellow), FIRAS [27] (red), and γ-ray attenuation [28] (pink) bounds (see Appendix A).The black band identifies the 95% CL excess detected by LORRI [10,22].
detector while pointing at the direction in the sky n, where ω piv is the pivot frequency.We define ϵ(ω) as a normalized throughput function, as we will compare the blue axion anisotropy spectrum to the true anisotropy spectrum, rather than the spectrum as seen in the detector, the Wide Field Camera 3 on board of the Hubble Space Telescope [34]. 1 Therefore, calling T(ω) the throughput of Ref. [35], which is defined as the number of detected counts/s/cm 2 of telescope area relative to the incident flux in photons/s/cm 2 , our normalized throughput function will be where the integral in the denominator is the efficiency as defined in Ref. [34].The fluctuations towards n can be expanded as spherical harmonics, while the relevant angular power spectrum is defined as In terms of the window function we obtain (see also Refs.[30,36]), where r(z) = z 0 dz/H(z) is the comoving distance, and j l (kr(z)) is the spherical Bessel function.The power spectrum is defined as ⟨δ k1 (r(z 1 ))δ k2 (r(z 2 ))⟩ = (2π) 3 δ (3) (k 1 − k 2 )P δ [k 1 , r(z 1 ), r(z 2 )].Since the power spectrum varies slowly with k, we can apply the Limber approximation [37,38], ; r(z 1 ), r(z 2 ) δ (1) [r(z 1 ) − r(z 2 )].(10) We are now able to find the correlation over each multipole moment due to the decay of relics, where the abundance ρ a and the non-linear spatial power spectrum P δ depend on the production mechanism.We expect the constraints on Γ a→γγ to be the most stringent when m a ≃ 2ω piv , with a sharp cutoff at smaller masses.At larger masses, the bounds get weaker, since the number density ρ a /m a becomes smaller and for a given ω piv the integral is dominated by the redshift z ≃ m a /2ω piv .

III. PRODUCTION MECHANISMS
As we have seen in the previous section, the angular power spectrum contribution from a decaying Big Bang relic depends crucially on its abundance ρ a and the nonlinear spatial power spectrum P δ [see Eq. (11)].In this section, we will explore a variety of production mechanisms, which will give rise to different abundances and power spectra.

A. Misalignment mechanism and decay of topological defects
As a first application, we will assume that blue axions constitute 100% of dark matter.For masses in the range we consider (5 − 25 eV), freeze-out is not a viable mechanism to produce the entirety of dark matter.Nevertheless, it is well known that the QCD axion [39][40][41], the pseudo-Nambu-Goldstone boson associated with the spontaneous breaking of the global U (1) PQ Peccei-Quinn symmetry introduced to solve the strong CP problem, is a good candidate for CDM, since it can be produced through the misalignment mechanism [42][43][44].The potential for the field ϕ = |ϕ|e iθ includes the terms where m a ∝ Λ 2 QCD η −1 , and we introduced the bias term ∝ ϵ b [45,46], possibly related to Planck-suppressed operators [47,48], to make the model cosmologically viable if N DW > 1 and U (1) PQ is broken after inflation.Mutatis mutandis, similar considerations can be applied to a pseudoscalar which does not solve the strong CP problem.
Each of the terms dominates at different times as the temperature of the Universe goes down.When the temperature of the Universe is about η ∼ N DW f a (f a being the axion decay constant), with N DW being the number of minima along the orbit of vacua, cosmic strings form and the phase assumes a value θ 0 , the initial misalignment angle.Once the Compton wavelength of the particle a = θη enters the horizon, the field starts oscillating with frequency m a , producing a QCD axion energy density corresponding to CDM.Notice that the mass depends on temperature, m a (T ), due to the nonperturbative effect of QCD.If the U (1) PQ symmetry is broken before or during inflation, the QCD axion abundance depends on θ 0 .If the symmetry is broken after inflation, axions can be produced in the decay of cosmic strings [49,50] (see Ref. [51] and references therein for recent works) and, depending on the UV completion, namely if N DW > 1, in the annihilation of walls bounded by strings (see e.g.[52]).The latter happens thanks to the bias term, which is needed in order for the N DW postinflationary scenario to be cosmologically viable, since otherwise the string-wall network would come to dominate the Universe.
Interestingly, the N DW = 1 post-inflationary scenario predicts a rather small allowed range for the mass of the QCD axion.The misalignment contribution to the abundance depends only on the QCD axion mass (in turn inversely proportional to the U (1) PQ symmetry breaking scale), since θ 0 is averaged over all its possible values.However, a very precise estimate of the contribution from cosmic strings is needed to predict the QCD axion mass as dark matter in the post-inflationary scenario, though (see e.g.Refs.[51,53]).The N DW > 1 scenario is possible only tuning the value of δ [52,54].
In the 5 − 25 eV mass range the QCD axion is largely excluded by several bounds, including the evolution of HB stars in globular clusters, since the coupling to photons is g aγγ = − α 2πfa C αγ , with C αγ = O(1) a model dependent number and α the fine structure constant.A similar bound, limiting the QCD axion mass to m a ≲ 10 −1 eV, arises from the constraints on the neutron electric dipole moment [55,56], the only model-independent bound as other couplings can be smaller at the cost of some fine tuning [57].Notice however that suppressing the coupling to neutrons (that limits the QCD axion mass to be m a ≲ 10 −2 eV [58]) would require even further nontrivial assumptions [59].Intriguingly, a heavier-thanexpected axion (i.e., with the QCD axion band moved to the right) can result from the existence of N degenerate Standard Model (SM) replicas, with the axion being the same particle in all the replicas [60], a generalization of the scenario proposed in Ref. [61] (see also [62]).In this case, the axion mass gets heavier by a factor √ N .Therefore, for the blue axion to be the QCD axion, one would need N = O(1000) copies of the SM.A more appealing scenario with a heavier-than-expected QCD axion might rely on a single mirror world with a QCD ′ dynamical scale much larger than the SM Λ QCD [63].(Other strategies to enlarge the parameter space for the QCD axion can be found in Ref. [60] and references therein.)Since a rather contrived scenario would be needed for the blue axion to be the QCD axion, one could give up the possibility of solving the strong CP problem, and assume the blue axion to be the pseudo-Nambu-Goldstone boson of a global U (1) symmetry whose breaking scale is unrelated to the coupling to photons g aγγ and to the mass m a .Under these assumptions, one can produce the correct DM abundance through misalignment mechanism and cosmic string decay [64].Moreover, if the orbit of vacua after the U (1) symmetry breaking admits multiple minima, the annihilation of walls bounded by strings due to a small bias in the energy density between the true vacuum and the other minima can easily dominate the production [65].In this case, the blue axion could be associated with the production of supermassive black hole seeds [66].Notice that depending on how large the bias term is (i.e., depending on the annihilation temperature T ann of the string-wall network), blue axions would form potentially at temperatures below 1 keV.In such case the matter power spectrum can be affected, and become similar to a warm (or, for even smaller T ann , hot) dark matter power spectrum.Since the string-wall network annihilation happens at T ann ∼ 10 11 keV ϵ b m a /eV [65,66], blue axions from the decay of the string-wall network will have a CDM power spectrum if 10 −22 ≪ ϵ b ≪ 1.While these figures might seem very unnatural, such small biases can arise from Planck-suppressed operators.Thus, one can produce a power spectrum with a cutoff corresponding to the redshift of the string-wall network annihilation, similar to other late-forming dark matter scenarios [67].
To summarize, we identify the misalignment mechanism and, in the post-inflationary scenario, the decay of cosmic strings and domain walls (the latter for a large enough bias term) as the mechanisms to produce the entirety of dark matter in the form of blue axions, ρ a = ρ CDM = Ω CDM ρ c , where Ω CDM = 0.12 h −2 and ρ c = 1.05 × 10 −5 h 2 GeV cm −3 .In this case, the nonlinear power spectrum P δ (z, r, r) is evaluated with the CLASS code [68], publicly available at [69], from redshift z = 0 to z = 12 with steps of 0.1.If the bias term is small, the production is dominated by the annihilation of the string-wall network.In this case, we will assume ρ a = ρ CDM = Ω CDM ρ c , and P δ (z, r, r) equal to the CDM power spectrum, with a cutoff at the co-moving wave-number k T = 7h Mpc −1 , marginally consistent with Lyman-α forest and galaxy clustering data [66,67,70]. 2

B. Alternative scenarios for non-cold Dark Matter
Axions can be copiously produced in the early Universe plasma, and behave as a NCDM component.In this case, the abundance is set by their different interactions, which depend on the considered model [72][73][74][75][76][77].We will focus on the coupling to photons, so that the thermal equilibrium is kept by means of the Primakoff process γQ → aQ, where Q refers to any charged particle.Additional interactions (e.g. to nucleons) can keep the axion in thermal equilibrium at smaller temperatures [56,75,77].
When the interaction rate Γ between axions and the SM particles is much larger than the expansion rate H, axions are in thermal equilibrium, and eventually decouple from the thermal bath when Γ/H < 1.Therefore, assuming the freeze-out to be instantaneous, decoupling happens at the so-called freeze-out temperature T F at which the condition H(T F ) ≃ Γ(T F ) is satisfied.By requiring Γ ≃ H, the freeze-out temperature for the Primakoff process is approximately given by [78] We notice that for g aγγ ≲ 10 −9 GeV −1 , T F is larger than the Electroweak (EW) scale Λ EW , above which the particle content of the plasma is speculative.After decoupling, the axion population established at the freeze-out temperature redshifts until today, with abundance [78,79] and temperature T a given by where g * ,s (T ) is the number of entropy-degrees of freedom [80], T 0 = 0.235 meV is the temperature of the Cosmic Microwave Background (CMB), at which g * ,s (T 0 ) = 3.91.Above the EW scale, the SM provides g * ,s (T F > Λ EW ) = 106.75degrees of freedom.Thus, from Eq. ( 14), axions with mass m a ≃ 155 eV decoupling at T F ≳ Λ EW would account for all the dark matter of the Universe since Ω a h 2 = Ω CDM h 2 = 0.12.This implies that, assuming a large reheating temperature, larger masses are excluded [78].However, the abundance of NCDM axions is constrained to be much smaller than the one of CDM by structure formation observations [29].A NCDM component has a significant free-streaming length, modifying the matter power spectrum on the smallest spatial scales and affecting the observations of the local Universe.In this context, assuming that the temperature of the NCDM component is the same of the standard neutrinos T NCDM = T ν = 0.716 T 0 and combining the prediction of the number of satellites galaxy with the CMB temperature, polarization and lensing measurements from the Planck satellite and the Baryon Acoustic Oscillation (BAO) data, the authors of Ref. [29] constrained the fraction of the NCDM component with respect to the total DM as a function of its mass in the range (10 −5 − 10 5 ) eV [29].
As shown by the solid line in Fig. 2 (obtained interpolating data from Fig. 5 of Ref. [29]), for m a ≃ O(10) eV the 2 σ constraint on for m a = 4 eV and Ω NCDM h 2 = 8.3×10 −3 for m a = 40 eV.Since the cosmological calculations depend on the ratio m NCDM /T NCDM , the bound found in Ref. [29] can be translated for a model where the NCDM axion has mass m a and temperature T a rescaling the mass accordingly, i.e.
This implies that the structure formation bound on f NCDM is a function of m NCDM = m a T ν /T a .Thus, for T a ≲ T ν it is shifted to lower masses (see, e.g., the dashed line in Fig. 2) and vice versa (dotted line).From Eq. ( 14), for m a ∼ O(10) eV, the axion relic abundance is much larger than the structure formation constraints [29].This implies that a freeze-out mechanism is ruled out in the context of the SM.We will consider NCDM axions produced in two alternative scenarios: • a modified freeze-out scenario, in which we assume additional degrees of freedom in order to increase g * ,s , reducing the axion abundance and temperature to satisfy the structure formation bounds; • a hot relic scenario, in which relic axions have the neutrino temperature T ν and their abundance saturates the bounds found in Ref. [29].

Modified Freeze-out
In this scenario, the blue axion decouples at a temperature larger than the EW scale.Since the particle content FIG. 2. 2 σ upper bound on the fraction fNCDM of NCDM as a function of the NCDM axion mass taken from Ref. [29] for Ta = Tν = 0.716T0 (solid line), and rescaled for Ta < Tν (dashed) and Ta > Tν (dotted).
of the plasma at this high temperature is speculative, there could be additional degrees of freedom increasing the value of g * ,s .For instance, in the minimal supersymmetric SM scenario, above the supersymmetry breaking energy scale g * ,s = 228.75[78], about twice as much as the SM prediction.From Eqs. ( 14) and ( 15), we see that a larger g * ,s (T F ) reduces the relic abundance and increases the today NCDM temperature.For an axion with mass m a and temperature given by Eq. ( 15), bounds on the abundance are rescaled following Eq.( 17) and f NCDM is a function of m a and g * ,s , since Thus, from Eq. ( 16), Since the abundance of a thermally produced axion is given by Eq. ( 14), for a fixed mass m a , equating Eq. ( 14) to Eq. ( 19) we find the value of g * ,s satisfying the structure formation constraints.Notice that the abundance does not depend strongly on the details of the decoupling, provided that the associated freeze-out temperature is larger than the EW scale.As an example, for m a = 10 eV, g * ,s (T F ) = 143.479 is required to satisfy the structure formation bound (independently of the exact value of T F > Λ EW ), implying a temperature T a = 0.301 T 0 .The 2 σ upper bound on f NCDM for an axion with T a = 0.301 T 0 is represented by the dashed line in Fig. 2 (the constraint for m a = 10 eV and T a = 0.301 T 0 is f NCDM ≲ 4.3 × 10 −2 , equal to the one obtained for a model with temperature T a = T ν and m a = 23.8eV, represented by the solid line in Fig. 2).For the largest mass considered in this work, m a = 25 eV, the structure formation constraints are satisfied by g * ,s (T F ) = 199.532,corresponding to T a = 0.270 T 0 .As discussed in Ref. [30], for a thermally produced NCDM axion population, the intensity spectrum is indistinguishable from the one given by decaying CDM.Therefore, for each value of the axion mass, we find the value of g * ,s satisfying the structure formation constraints and we evaluate the angular power spectrum using Eq. ( 11) with the abundance ρ a obtained from Eq. ( 14) and the non-linear NCDM power spectrum P δ computed in the adiabatic approximation.In this approximation, if T NCDM and T CDM are the transfer functions of NCDM and CDM, respectively, the non-linear NCDM power spectrum is [81] Both the transfer functions, relating the primordial and the present day power spectra [82], and the nonlinear CDM power spectrum P δ, CDM are obtained using CLASS [83], with parameters m a , Ω a from Eq. ( 14) and T a from Eq. ( 15), computed using the value of g * ,s satisfying the structure formation bound.

Hot relic
Depending on the presence of additional couplings or particles, blue axions could have a large temperature.We define as "hot relic" a particle with the same temperature as neutrinos but with a much larger density, saturating the structure formation bounds, and stay agnostic about its production mechanism.In this scenario, we compute the angular power spectrum C l using Eq. ( 11), with the axion relic density saturating the structure formation constraint represented by the solid line in Fig. 2 [without satisfying Eq. ( 14)] and the non-linear NCDM power spectrum given by Eq. ( 20), where T NCDM , T CDM and P δ,CDM are computed with CLASS with input parameters m a , T a = T ν = 0.716 T 0 and Ω a from Ref. [29] (the solid line in Fig. 2).
We mention here that blue axions could be produced after the decay of heavy dark sector particles, with an energy related to the decaying particle mass.Since the more energetic are the axions, the larger is the effect on structure formation, the relic axion abundance would be more constrained (see the dotted line in Fig. 2), while the non-linear power spectrum P δ would be suppressed at larger scales.We leave the study of this scenario for future work.

Freeze-in
It is possible that axions never reach equilibrium with the thermal bath.Nevertheless, they can still be produced and linger around as dark relics (they "freeze-in").FIG. 3. Angular power spectrum of the COB anisotropy for the DM axion decaying into photons, for λ obs = 606 nm, in the case of CDM (solid lines), freeze-out (dashed), hot relic (dotted) and CDM produced by annihilation of domain walls (thin-solid).We have fixed the axion mass ma = 10 eV (black lines) and ma = 15 eV (red) and the product ρa Γa→γγ = 10 −29 s −1 GeV cm −3 .The HST measurements at 606 nm are shown with error bars [20].
We consider that the axion production starts at the reheating temperature T RH , when the Universe enters its last phase of radiation domination.For large T RH , axions have more time to be produced and their abundance is higher.The lowest value of the reheating temperature compatible with observations is T RH = 5 MeV [84][85][86][87][88][89].As further discussed in Ref. [90], axions with mass m a ≲ T RH and interacting with photons are mainly produced via the Primakoff effect and the relic density is given by [90] Ω a h 2 ≃ 10 −5 m a eV For the purpose of this work, the freeze-in paradigm will not give an axion relic abundance large enough to have observable effects, so we will not discuss it any further.

IV. RESULTS
In this section we compare the resulting angular power spectra in the different scenarios previously discussed with the HST measurements in order to constrain the axion parameter space.In particular, we use the dataset at the shortest available wavelength, namely 606 nm, obtained by the Wide Field Camera 3 and the Advanced Camera for Surveys, which covers 120 square arcminutes in the Great Observatories Origins Deep Survey [20].FIG. 4. Comparison between our bound on CDM using HST measurements at 606 nm (solid blue line) and the bound obtained in Ref. [23] (the dotted blue line).We also show other bounds.The color code is the same as Fig. 1.

A. Bounds and future reaches on cold dark matter
In Fig. 3 we show the angular power spectrum l 2 C l /2π as measured by HST at 606 nm with error bars, and the angular power spectrum from decaying CDM as seen by the detector at 606 nm, at fixed value of ρ a × Γ a→γγ = 10 −29 s −1 GeV cm −3 for m a = 10 eV (black lines) and m a = 15 (red lines).We constrain the blue axion lifetime by requiring that the angular power spectrum C l in Eq. (11) does not exceed the upper error bar of any of the data points.In the case of CDM (solid lines) the quantity l 2 C l is peaked at l ≳ 10 4 and the peak is shifted towards smaller scales as the mass increases.In this context, the most constraining bin is the one at l ≃ 4.8 × 10 5 and we exclude the light blue region delimited by the solid line and labelled as HST 606 nm in Fig. 1, finding good agreement with the results in Ref. [23] for m a ≳ 10 eV.In particular, as shown in Fig. 4, where the bound from Ref. [23] is represented by the dotted blue line, our bound is stronger by ∼ 30% at m a = 10 eV, by ∼ 15% at m a = 15 eV and the discrepancy becomes negligible for even larger masses.For smaller values of the mass, we find a larger discrepancy between the two results, with our bounds stronger by a factor ∼ 4 at m a ≃ 5 eV.This could be due to the several differences in our analysis.Differently from Ref. [23], we get the non-linear power spectrum directly from CLASS.In addition, we characterize the detector response more precisely, by considering the detector ω piv and the throughput (see Appendix B) instead of using an observational bandwidth equal to the observational frequency.We checked that using the same detector response as the one used in Ref. [23], the discrepancy at lower masses is reduced.
Stronger constraints on the axion lifetime could be obtained with HST measurements at shorter wavelengths.FIG. 5. Non-linear spatial power spectrum for CDM (black solid line), freeze-out (dashed) and hot relic (dotted) at redshift z = 1 for ma = 10 eV (black lines) and ma = 15 eV (red).The vertical dashed line at lT = 1.6 × 10 4 corresponds to the cutoff comoving wave-number kT = 7 h Mpc −1 for CDM produced by annihilation of domain walls.In this latter case, the non-linear power spectrum is assumed to be the one for CDM (solid black line), cut at lT .In order to forecast the possible future reaches, we evaluate the angular power spectrum in Eq. ( 11) using ω piv and ϵ(ω) for measurements at 438 nm and 336 nm, respectively (see Appendix B for more details) and we require that it does not exceed the upper error bar of any data points at 606 nm.In this way we probe the light blue regions delimited by dashed lines in Fig. 1, labelled with HST 438 nm and HST 336 nm, respectively.In the case of 336 nm, the bound would be improved by a factor 4-10 in the mass range 10 − 25 eV, leading to the strongest constraints in this region.We stress that our assumptions are conservative, since at smaller wavelengths the measured angular power spectrum is expected to be even smaller, implying even stronger projected reaches.

B. Alternative scenarios
The anisotropy bounds can be relaxed in alternative scenarios.The resulting angular power spectrum C l can be smaller due to the reduction in the abundance ρ a and the suppression of the non-linear power spectrum P δ , depending on the production mechanism.In Fig. 3 we show the resulting angular power spectrum for the cases of freeze-out (dashed) and of a hot relic with neutrino temperature (dotted).The effect of the abundance suppression is blurred out since the angular power spectra are shown at fixed value of ρ a Γ a→γγ (a larger value of Γ a→γγ compensates the suppressed ρ a ).On the other hand, the effect of the non-linear power spectrum suppression is clearly visible, since for NCDM the quantity l 2 C l is suppressed at l ≲ 10 4 .The suppression point is at larger scales for hotter DM and, for a fixed cosmo- FIG. 6. Bounds and future reaches on the lifetime of the blue axion in the case of freeze-out (upper left), hot relic (upper right) and CDM produced by annihilation of domain walls (lower panel).We show also other bounds.The color code is the same as Fig. 1. logical scenario, the suppression starts at smaller scales as the mass increases.Indeed, this reflects the trend of the non-linear power spectrum, as shown in Fig. 5 at a representative value of the redshift z = 1 for m a = 10 eV (black) and m a = 15 eV (red) in the different scenarios that we considered.The most conservative scenario is the hot relic case, since the angular power spectrum is suppressed at larger scales.Again, by requiring that the computed C l must not exceed the error bar of any data points at 606 nm, we compute the bounds for the observations at 606 nm and forecast the sensitivity for future observations.In Fig. 6 we show how the bounds on the axion lifetime are relaxed in the freeze-out (upper-left panel) and hot relic (upper-right panel) scenarios (for a discussion on the other bounds and how they change in the different scenarios see Appendix A).In these scenarios, future measurements at shorter wavelengths would improve the current anisotropy bound by almost one order of magnitude, setting the strongest probe on the blue axion lifetime for 10 eV ≲ m a ≲ 20 eV (freeze-out) and 15 eV ≲ m a ≲ 18 eV (hot relic).
Finally, we show what happens in the case in which axions are produced by the annihilation of a string-wall network.We model this case using a CDM P δ (the black solid line in Fig. 5) featuring a cutoff at the co-moving wave-number k T = 7 h Mpc −1 , represented by the vertical dashed line in Fig. 5.As shown by the thin-solid lines in Fig. 3, in this scenario the resulting angular power spectrum is the same as the CDM case at large scales and it is strongly suppressed at l ≃ 10 4 , due to the cutoff previously discussed.Therefore, the most constraining data point is at larger scales compared to the CDM case (depending on the axion mass) and the constraints are relaxed by a factor ∼ 2, as shown in the lower panel of Fig. 6.However, also in this case, the HST measurements at 606 nm already exclude the CDM interpretation of the LORRI excess. 1 FIG.7. Best fit in the absence of axions (left panel) and in the case of CDM axion with ma = 10 eV (right panel).In both the panels, we show the HST measurements with error bars, the shot noise (dashed black line), the foreground (dotted black), the axion (solid black) contributions and their sum (red line).

C. Model constraints
Our bounds on the blue axion lifetime are conservative since we require that the blue axion decay contribution saturates the COB anisotropy spectra.Moreover, they are robust, since they parametrically depend on the square root of the angular power spectrum.Stronger constraints can be obtained if we model the HST data through three different components, including the axion decay contribution, the shot-noise term and a power-law component accounting for Galactic foregrounds, e.g. the diffuse Galactic light (DGL) [21].In principle, one should consider also the signal from high-z faint galaxies during the epoch of reionization.However, we neglect it since for the HST measurements at 606 nm there is almost no contribution from high-z signal [20,21].At small scales, the power spectrum is dominated by the shot noise, which is scale-independent and with an angular power spectrum given by [21] where A shot is the shot-noise amplitude factor, which is constant for a given observed wavelength.At large scales, the power spectrum is dominated by C f l for foregrounds, which can be described by [21] where A f is the amplitude factor for the foregrounds.For instance, a possible dominant component in C f l could be the DGL, since the DGL term is proportional to ∼ l −3 .Therefore, closely following Ref.[21], the HST data can be fit as with C axion l given by Eq. ( 11).In order to constrain the axion parameter space, we perform a likelihood analysis, by defining the likelihood function L ∝ exp(−χ 2 /2), with χ 2 given by where N d = 13 is the number of data points, C obs l,i is the angular power spectrum from the i-th observational point with error σ i , and C th l,i is the theoretical angular power spectrum evaluated at l corresponding to the i-th data point.We assume the free parameters to follow a flat prior distribution, with A shot ∈ (10 −13 , 10 −11 ) and A f ∈ (10 2 , 10 4 ) [20].In the left panel of Fig. 7 we show the best fit in the absence of axions, with A shot,0 = 3.66× 10 −12 and A f,0 = 1.91 × 10 3 .As shown by the value of the chi squared χ 2 0 = 9.2 over 11 degrees of freedom, the fit well reproduces the data points, with the largest discrepancy for l ∼ 10 5 where data are underestimated.In order to constrain the axion, we marginalize over A shot and A f by minimizing the χ 2 for each value of m a and g aγγ over A shot and A f , χ2 = min If axions are CDM, for certain values of the coupling, their presence can improve the fit.In this case, as shown in the right panel of Fig. 7 for m a = 10 eV, the quantity l 2 C axion l is peaked at l ∼ 10 4 − 10 5 , reducing the discrepancy between data and fit at those scales.For instance, for m a = 10 eV, the best fit is obtained for g aγγ = 1.24 × 10 −11 GeV −1 , A shot = 3.56 × 10 −12 and A f = 1.28 × 10 3 .However, we mention that this should not be considered as a hint for the blue axion.Indeed, axions improve the fit due to the uncertainty related to the FIG. 8. Difference between the bound from HST at 606 nm obtained requiring that C l must not exceed the upper error bar of any data points (conservative, dashed blue line) and the 95% CL constraint from the likelihood analysis (solid blue line) in the case of CDM (upper-left panel), freeze-out (upper right), hot relic (lower left) and CDM produced by annihilation of domain walls (lower-right panel).We show also the other bounds and the excess in LORRI data with the same color code as Fig. 1 standard physics contributions and including them we overfit the data, obtaining χ 2 min = 2.02 over 10 degrees of freedom for m a = 10 eV.Moreover, the cross-correlations between intensity fluctuations at different wavelengths due to dark matter decay is zero (as structures at different redshifts are mainly uncorrelated), contrary to the observations [21].Therefore, more than one axion (in Ref. [21] a continuous distribution of axion masses is assumed) is needed to fit the cross-correlations.
In the case of NCDM, the presence of axions would not improve the fit, since the angular power spectrum is suppressed at larger scales l ≲ 10 4 , as shown in Fig. 3.In this case, the minimum of χ 2 is χ 2 0 = 9.19 for vanishing axion-photon coupling.We define the test statistic where χ2 min is the minimum chi squared, obtained at the value of the axion photon coupling g aγγ = g min (g min ̸ = 0 only in the CDM case).This quantity follows a half-chisquared distribution [91], which allows us to constrain the axion parameter space at 95% Confidence Level (CL) by requiring χ 2 * ≤ 2.7.As shown in the upper-left panel of Fig. 8, in the CDM case this approach is consistent with the conservative one previously discussed and it strengthens bounds on the rate by ∼ 15% at m a = 5 eV and ∼ 25% at m a = 25 eV.In the other cases, this approach leads to factor ∼ 2 stronger constraints for m a ≲ 10 eV, with a lower discrepancy at larger masses, as shown in the remaining panels of Fig. 8. Indeed, when alternative cosmological scenarios are considered, the angular power spectrum is suppressed at larger scales.Since with our conservative approach we set the bound when the angular power spectrum exceeds any of the data points, the most constraining one is at l ≲ 104 , where l2 C l is larger, without taking into account all the other measurements.On the other hand, by performing a likelihood analysis that includes additional contributions, all the data points are important since by definition the χ 2 accounts for the entire data set, as shown in Eq. ( 25), allowing us to set stronger bounds.

V. THE DARK PORTAL
To circumvent stellar cooling bounds while relaxing anisotropy bounds, we will assume the existence of a small hot dark matter component decaying through a dark portal featuring a very light dark photon (m χ ≪ m a ) [30,[92][93][94], with a decay rate Notice that the number of photons produced per unit time by the coupling in Eq. ( 28) is equal to the ones from the coupling in Eq. ( 1) for g aγχ = g aγγ ≡ g. 3 This interaction was assumed in Ref. [30] as a possible interpretation of an excess observed in the cosmic infrared background by CIBER [95], a sounding rocket equipped with infrared cameras.A similar idea has been recently advanced in Ref. [23] to explain LORRI excess, but no dedicated analysis was carried out.We constrain the dark portal model using the same strategy discussed in Sec.IV, i.e. by requiring that the angular power spectrum C l in Eq. ( 11) must not exceed the upper error bar of any of the data points.Since for this model the number of photons produced per unit time is the same obtained from the coupling in Eq. ( 1), the same bounds apply for both models in all the scenarios under consideration.Indeed, in the case of CDM, the production mechanisms are independent of the coupling g.Thus, the angular power spectrum C l is computed using ρ a = ρ CDM and the CDM power spectrum P δ evaluated with the CLASS code, cutting it at the co-moving wave-number k T = 7 h Mpc −1 when the production is dominated by the annihilation of domain walls.In the case of NCDM, the blue axion is thermally produced mainly via pair annihilations in the s channel e + + e − → a + χ, while plasmon decays γ * → a + χ give a negligible contribution in the early Universe plasma.The freeze-out temperature for the pair annihilation is [30] T F ≃ 4.8 × 10 3 GeV 10 −9 GeV −1 g aγχ also in this case larger than the EW scale for g aγχ ≲ 10 −9 GeV −1 .Thus, also for the dark portal model we consider NCDM axions with temperature equal to the neutrino one or set by the modified freeze-out scenario described in Sec.III B 1.For both the scenarios, we are agnostic of the production processes and the angular power spectrum C l is evaluated using ρ a and P δ, NCDM as described in Sec.III B 1 (for the freeze-out) and Sec.III B 2 (for the hot relic case), obtaining the same value of C l for g aγγ = g aγχ = g.In Fig. 9 we show the bounds and projected reaches on the coupling g, valid for both g aγχ and g aγγ .For axions interacting with two photons, the HB bound (horizontal dotted green line) applies, excluding g aγγ ≳ 0.65×10 −10 GeV −1 [25] and ruling out the NCDM interpretation of the LORRI excess in all the mass range we consider.On the other hand, this constraint vanishes for the dark portal since Primakoff emission is not possible.In this model, plasmon decays in stars can copiously produce axion-dark photon pairs, affecting the standard stellar evolution.The plasmon decay rate was explicitly derived in Ref. [30], and found to be identical to the plasmon decay rate to neutrinos with a magnetic dipole moment µ ν , with the substitution µ ν → g aχγ /2. 4  The strongest bounds to date on the neutrino magnetic dipole moment comes from the brightness of the tip of the red-giant branch [96].The best galactic target is ω Centaury, from which µ ν < 1.2 × 10 −12 µ B at 95% CL with µ B = e/2m e the Bohr magneton, much more stringent than previous bounds, e.g.[97].Therefore, we find g aχγ < 7.1 × 10 −10 GeV −1 (the horizontal solid green line).This is the strongest bound on the dark portal in most of the parameter space, and supersedes all previous astrophysical constraints for any mass smaller than O(10 keV) [93,94].Nevertheless, the NCDM dark portal interpretation of the excess is excluded.Future anisotropy measurements at 438 nm and 336 nm will represent the strongest probe on the dark portal model for 8 eV ≲ m a ≲ 24 eV in the case of CDM, 8 eV ≲ m a ≲ 20 eV for the freeze-out scenario and 8 eV ≲ m a ≲ 18 eV for the hot relic case.

VI. LINE INTENSITY MAPPING PROJECTED REACH
We here briefly comment on the promising possibilities of line intensity mapping (LIM) [98,99].Similarly to the anisotropy measurements previously described, rather than identifying galaxies, LIM aims at measuring the integrated emission of spectral lines from both galaxies and the intergalactic medium (IGM) with small low-aperture instruments.The line-of-sight distribution is then reconstructed through the frequency dependence: FIG. 9. Bounds and future reaches on the coupling g of the blue-axion in the case of CDM (upper-left), freeze-out (upper-right), hot relic (lower-left) and late-time produced CDM produced by annihilation of domain walls (lower-right panel).If g = gaγγ the HB bound (the light green region delimited by a dotted line) applies.If g = gaγχ the HB bound vanishes and the RG bound (horizontal solid green line) must be considered.We show also other bounds, valid for both gaγγ and gaγχ.The color code is the same as Fig. 1.
the redshift of a targeted line is obtained comparing the observed frequency with the frequency at rest.A decaying relic would potentially show up in future surveys as an "interloper line".The approach, based on previous ideas of Refs.[21,26], was proposed in Ref. [100], and finally applied to realistic forecasts of decaying DM [101] a → γγ and neutrinos [102] ν i → ν j γ, where ν i is an active massive neutrino with mass m i decaying into a lighter eigenstate ν j with mass m j .The possible observables are the LIM power spectrum and the voxel intensity distribution (VID), i.e. the distribution of observed intensities in each voxel, the three-dimensional equivalent of pixels [101,103].Here we focus on the most powerful projections, which are obtained with VID measurements.
First, in Fig. 10 we reproduce the results of Ref. [102] for the normal hierarchy (the same can be done for the inverted hierarchy).To obtain the dashed curves, we simply rescale the reach of LIM searches for DM decay of Ref. [101] to the abundance of neutrinos.Namely, we find the projected reach on the neutrino decay rate Γ ν as a function of the sum of the neutrino masses m ν , requiring that where n ν = 56 cm −3 [33] is the neutrino number density per flavor, n DM = ρ DM /m a is the DM number density and Γ Bernal DM is the projected reach on DM decay found in Ref. [101] as a function of the DM mass m a .Here, we rescale the DM mass by replacing m a → (m 2 i − m 2 j )/m i since the rest-frame energy of the emitted photons is m a /2 in the case of DM decays, and (m 2 i − m 2 j )/2m i for neutrino decay.Thus, in the case of a transition ν i → ν j , the rescaled projected reach on Γ ν is explicitly given by where m 2 j − m 2 i is a fixed parameter [33] and m i depends on m ν .In the upper panel of Fig. 10, we show the results of our rescaling (the dashed lines) in the case of ν 3 → ν 1 decay, with m 2 3 − m 2 1 = 24.53× 10 −4 eV 2 [33], while in the lower panel we plot the projected reach on neutrino decay for the transition ν 2 → ν 1 , with m 2 2 − m 2 1 = 0.753 × 10 −4 eV 2 [33].In both cases, we find good agreement between our dashed curves and the results of the dedicated analysis of Ref. [102] (the thin-solid lines).This shows that, as for the previously discussed COB anisotropies, the shape of the power spectrum modifies the prediction only negligibly, and one just needs to account for the abundance.Notice, however, that the projected reach for DM decay has been recently computed again by the same authors in an erratum [101].Therefore, we revisit the neutrino decay rate reach (see the thick-solid lines) by rescaling the DM reach of the erratum.The reach is less competitive than initially claimed in Ref. [102], and weakened to the level of future CMB spectral distortion probes [104].Moreover, it is several orders of magnitude larger than astrophysical probes [96], given by the sum of all the possible diagonal and non-diagonal electric and magnetic moments, so the latter are to be interpreted as conservative bounds on the lifetime of a specific neutrino mass eigenstate.On the other hand, the NCDM component can be much more abundant than neutrinos, offering an interesting target to forthcoming LIM searches.
Without further ado, we can obtain the reach for NCDM by simply rescaling the projected reach for CDM by a factor Ω CDM /Ω a , with Ω a depending on the cosmological scenario.Thus, we compute the LIM projections and show them in the m a − Γ plane for CDM (taken from the erratum) in the left panel of Fig. 11 and the most conservative hot relic case (right panel), together with other bounds (solid lines) and the strongest reaches (dashed lines) on the blue axion.From the updated Ref. [101] we see that the Hobby-Eberly Telescope Dark Energy Experiment (HETDEX) [105] can probe CDM blue axions with lifetime 10 29 s, and by a simple rescaling argument we find HETDEX can probe also hot relic blue axions with lifetime 10 27 s.

VII. DISCUSSION AND CONCLUSIONS
In this paper, we have revisited the parameter space of a radiative decaying Big Bang relic with a mass m a ≃ 5−25 eV that we have dubbed "blue axion".We explored the case of the decay to two photons a → γ + γ and the case of decay to a photon and a dark vector, a → χ + γ.
Existing bounds in this part of the parameter space depend on the production mechanism that dictates both FIG.10.Projected 95% CL sensitivity on the neutrino lifetime in the case of normal hierarchy as function of the total neutrino mass from VID measurements for all LIM surveys considered in [102] (thin-solid lines), namely COMAP1 (orange), COMAP2 (red), CCAT-prime (dark yellow) and At-LAST (light blue).The other coloured lines represent the rescaled sensitivities on neutrinos, obtained from the LIM searches for DM decay shown in the published version of Ref. [101] (dashed) and erratum (thick-solid).We show also CMB [104] (dotted black) limits and the astrophysical bound from RGs [96] (solid black).In the upper panel the transition between mass eigenstates ν3 → ν1 is considered, while in the lower panel ν2 → ν1 is assumed.
the abundance and the power spectrum, except for stellar bounds (evolution of HB stars for the two-photon coupling, and of RGs for the dark portal).Therefore, we have considered four different scenarios-namely, cold dark matter produced through the misalignment mechanism (with additional contributions from strings and rapidly decaying domain walls), a late-time produced cold dark matter from the annihilation of a long-lived string-wall network, a thermally produced population with temperature set by a freeze-out mechanism, and a hot relic with the temperature of neutrinos and the maximum abundance allowed by structure formation bounds.
We found that cosmic optical background anisotropies are a powerful probe, as the results we obtain depend mostly on the abundance, and only marginally on the power spectrum.We have proposed the observation of the COB anisotropies with HST at the pivot wavelengths of 438 nm and 336 nm, which has the most competitive discovery reach to date.In the near future, line inten-FIG.11.Projections on the lifetime of the blue axion for LIM VID (dashed gray line) for CDM (left panel) and the hot relic case (right panel).We show also other bounds and reaches with the same color code of Fig. 1. sity mapping will be a powerful probe of the blue axion parameter space, probing lifetimes as large as 10 29 s or 10 27 s assuming the blue axion to be the totality of cold dark matter or a hot relic with the abundance limited by structure formation bounds, respectively.As a byproduct, our analysis also excludes the dark portal model interpretation of the excess recently detected by LORRI on board of the New Horizons mission, lending further credibility to an astrophysical interpretation of the data.However, the experimental directions proposed in our paper will explore presently unprobed coupling strengths.Therefore, blue axions could still be discovered in the nearby future.
where V cl is the volume of the cluster, and d cl is the distance corresponding to a redshift z cl (z cl = 0.233 for A 2667 and z cl = 0.228 for A 2390).
Two galaxy clusters were observed by VIMOS in the period June 27-30, 2003 [106].The signature of decaying axions is an emission line tracing the DM density in the cluster.The absence of such a signal results in a constraint on axion CDM in the range 4.5 eV ≲ m a ≲ 7.7 eV and g aγγ ≲ 3 × 10 −12 GeV −1 , as shown in Fig. 9.
Note that Eq. (A1) scales like ∼ g 2 aγγ Ω a and the bounds in Ref. [26] (see 7 therein) were obtained assuming a NCDM with the relic density Ω NCDM h 2 = m a /130 eV, which is excluded by structure formation bounds [29].Moreover, besides the axionphoton coupling, additional interactions are needed to guarantee this relic density, obtained from Eq. ( 19) by setting g * ,s = 10.To obtain the constraint in the CDM case shown in Fig. 1 we need to multiply the intensity in Eq. (A1) by a factor Ω CDM /Ω NCDM . 5

FIRAS
The decay of a massive particle leads to the injection of photons in the primordial bath producing spectral distortions of the CMB.This effect is more efficient when the energy of the produced photon exceeds the neutral hydrogen excitation energy (13.6 eV).However, any source of distortion is constrained by the fact that the energy spectrum of the CMB is compatible with a perfect blackbody at ∼ 2.7 K [33].
In this context, in Ref. [27], data from the experiments COBE /FIRAS, EDGES and Planck were used to constrain an exotic energy injection in the primordial plasma.
Since the energy injected is proportional to the abundance of decaying DM, the bounds on the rate shown in Fig. 1 can be rescaled to different scenarios by multiplying the squared coupling from Fig. 29 of Ref. [27] (valid for CDM) by a factor Ω a /Ω CDM , with Ω a depending on the considered cosmological scenario.

γ-ray attenuation
High-energy γ rays from far sources are attenuated by scattering on low-energy photons of the EBL and subsequent pair production, γ + γ → e + + e − .One can determine this attenuation as a function of the source redshift and the observed γ-ray energy through joint analyses of the observed blazar spectra.In this context, Ref. [30]  advanced the idea that the blazar γ-ray spectrum is sensitive to the decay of eV-scale axions decaying into photons (see also Ref. [108,109]).In the recent Ref. [28], blue axions are constrained by performing a likelihood analysis, comparing the predicted optical depths due to axion decays with optical depth measurements obtained from the observations of almost 800 blazars [17,110].More precisely, the blue axion decay contribution is compared with the residual optical depths obtained by subtracting from EBL three standard components, namely the emission from galaxies at z < 6, the emission from galaxies at z > 6, and the intra-halo light (IHL)-emission from a faint population of stars tidally removed from galaxies [19,111,112].Since the contribution from galaxies at z < 6 vastly dominates the EBL, in Ref. [28] the most conservative bound, represented by the solid line in Fig. 12, is obtained increasing the uncertainties of the astrophysical EBL from galaxies at z < 6 to saturate the measured EBL for all the considered source redshifts and observed γ-ray energies.
To give some intuition about the nature of these bounds, here we reproduce the γ-ray constraints within at most a factor of O(1) through a simpler, less refined approach, requiring that the optical depth due to axion decay does not exceed the upper error bar of any data point for the optical depth at z s = 2.4, the largest source redshift shown in Ref. [28] (see Fig. 2

therein).
The optical depth is defined in terms of the γ-ray mean free path l as where z s is the source redshift.Let us consider a γ-ray photon observed on Earth with an energy E γ .At redshift z, it has energy ϵ γ = (1 + z)E γ , and can potentially scatter on a low-energy photon of the EBL with energy ϵ.Electron-positron pair production can occur above the EBL photon energy threshold where m e is the electron mass and µ is the cosine of the incidence angle of the two photons.The pair production cross section is [113], Here, Ω a is the axion density parameter and z * ≡ m a c 2 (1 + z)/(2ϵ) − 1 is the redshift of decay that contributes to the photon energy and redshift of interest.For fixed values of the axion mass and coupling, plugging Eq. (A5) into Eq.(A2) we obtain the optical depth τ as a function of the observed γ-ray photon energy.If one assumes that axions constitute all of the CDM, Ω a = Ω CDM , requiring that the computed τ must not exceed the upper error bar of any of the measured optical depths at z s = 2.4 [28], the dotted line in Fig. 12 is obtained, reproducing the γ-ray bound in Ref. [28] within a factor 2 for all masses in the range of interest.
Since the optical depth τ scales as Ω a Γ a , the bounds on the rate shown in Fig. 1 are obtained by multiplying the CDM constraint (taken from Ref. [28]) by a factor Ω a /Ω CDM , with Ω a depending on the considered cosmological scenario.Notice that these bounds are completely independent of the relic power spectrum.We also stress that neutral hydrogen absorption (for axion masses above 20.4 eV) can modify the constraints.Hence, we show this region as dotted in all the relevant plots of the main text.

Horizontal branch stars
The axion-photon interaction turns on the production of light axions in stars through the Primakoff process [114], altering the stellar evolution.The best astrophysical probes of the axion-photon coupling are Horizontal Branch (HB) stars [115][116][117].Axions produced by Primakoff processes reduce the lifetime of HB stars, without affecting the previous red giant (RG) phase because of the large electron degeneracy and the high plasma frequency which prevents an efficient axion production in degenerate stars [115].
Therefore, the axion-photon interaction can be constrained by means of the R parameter, R = N HB /N RGB , which compares the number of stars in the HB (N HB ) and RGB (N RGB ) phases, or equivalently, the duration of these phases.This observable is also sensitive to the helium mass fraction Y , introducing a degeneracy between g aγγ and Y because a lower lifetime due to the presence of axions could be compensated by the increase of the helium abundance.Recent observations obtained a value of the R parameter equal to R ave = 1.39 ± 0.03 [118].By comparing the R parameter computed through stellar simulations including axions with the measured one, the constraint derived on axions is [25,119] g aγγ < 0.65 × 10 −10 GeV −1 (95% CL) . (A7) A similar bound can be obtained by requiring the exotic energy ϵ exotic emitted per unit time and mass to be ϵ exotic ≲ 10 erg g −1 s −1 , to be computed assuming a onezone model for the core of HB stars, T ≃ 8.6 keV and ρ ≃ 10 4 g cm −3 for the temperature and density respectively (see, e.g., [116,120,121]).
Recently, a different observable has been proposed, namely the ratio of asymptotic giant branch (AGB) to HB stars (the R 2 -parameter) [122].This parameter was measured thanks to the HST photometry of 48 globular clusters, obtaining R 2 = 0.117±0.005[123].This observable is known to be quite insensitive to the initial helium mass fraction, in contrast with the R parameter [123].The constraint found in this way improves the previous one to the value g aγγ < 0.47 × 10 −10 GeV −1 .Evidences from asteroseismology might help to improve this modeling [124], leading to a further improvement of the bound to g aγγ < 0.34 × 10 −10 GeV −1 [122].We mention that, despite being more stringent that the one based on the R parameter, the bound from the R 2 parameter is affected by larger uncertainties related to the description of convective phenomena.(WFC3).The WFC3, combining two ultraviolet/visible (UVIS) CCDs with a NIR HgCdTe array, is capable of direct, high-resolution imaging over the entire wavelength range 200 to 1700 nm, with the UVIS channel sensitive to 200 nm ≲ λ ≲ 1000 nm, and the IR channel 800 nm ≲ λ ≲ 1700 nm.In this work, we are interested in observations through the UVIS channel, in which the detectors are two 4096 × 2051 pixel CCDs (namely UVIS 1 and UVIS 2), butted together to yield a 4096 × 4102 light-sensitive array with a ∼ 31 pixel (1.2 arcsec) gap.In order to characterize the anisotropy spectrum as seen in the detector, we consider three wide-band filters from UVIS 2,6 namely 606 nm, 438 nm, and 336 nm, with the pivot-frequency ω pivot in Table I and the throughput T(ω) shown in Fig. 13, defined as the number of detected counts/s/cm 2 of telescope area relative to the incident flux in photons/s/cm 2 .In this way, the normalized throughput function used in Eq. ( 11) is where the integral in the denominator is the UVIS 2 efficiency as defined in Ref.
[34] and reported in Table I.

10 FIG. 12 .
FIG.12.The most conservative γ-ray bound obtained in Ref.[28] (solid pink line) and the one estimated in this Appendix (dotted pink line).