Cosmic infrared background excess from axion-like particles and implications for multi-messenger observations of blazars

The first measurement of the diffuse background spectrum at 0.8-1.7 $\mu \rm{m}$ from the CIBER experiment has revealed a significant excess of the cosmic infrared background (CIB) radiation compared to the theoretically expected spectrum. We revisit the hypothesis that decays of axionlike particle (ALP) can explain this excess, extending previous analyses to the case of a warm relic population. We show that such a scenario is not excluded by anisotropy measurements nor by stellar cooling arguments. Moreover, we find that the increased extragalactic background light (EBL) does not contradict observations of blazar spectra. Furthermore, the increased EBL attenuates the diffuse TeV gamma-ray flux and alleviates the tension between the detected neutrino and gamma ray fluxes.


I. INTRODUCTION
Recently, the Cosmic Infrared Background Experiment (CIBER) collaboration has claimed the detection of an unexpectedly high flux compared to theoretical expectations in the 0.8-1.7 µm range of wavelengths [1]. This measurement is complementary to other observations in the infrared band like the ones carried by AKARI [2] and IRTS [3]. Even if an astrophysical explanation of the detected excess or systematic errors are not ruled out, it is worthwhile to speculate about a possible flux due to Big Bang relics, such as an axion-like particles (ALP) with mass around 1 eV. ALPs generalize the concept of the axion, introduced to solve the so-called Strong CP problem, which has multifaceted phenomenology [4]. The contribution of ALP decays was examined in Ref. [5]. Here we will revisit the hypothesis, taking into account the detector energy resolution, the possibility of warm dark matter, and the implications of increased EBL for blazar time-variability and multi-messenger observations. While a solid lower bound to the Cosmic Infrared Background (CIB) radiation can be obtained through deep sky galaxy counts [6], the precise shape and intensity of the diffuse, unresolved spectrum in the nearinfrared wavelength range is still unknown. Direct measurements [2,7,8] are difficult because of the large uncertainties caused by zodiacal light. Theoretical models are also subject to uncertainties, which result in different predictions [9][10][11]. Very high energy gamma rays from blazars have been used to set an upper limits on infrared background radiation [12], but the possibility of secondary gamma rays produced by cosmic rays along the line of sight [13][14][15][16][17] undermine these upper bounds.
The uncertainties make it difficult to identify any ad- * kalashev@inr.ac.ru † kusenko@ucla.edu ‡ edovita@mpp.mpg.de ditional contribution to the Extragalactic Background Light (EBL) besides the standard flux due to galaxy emission. Possible enhancements could come from ultraviolet redshifted photons produced by bottom-up astrophysical accelerators, ranging from high redshift galaxies [18] to black holes [19].
In the last few years, searches for indirect probes of portals connecting the Standard Model of particle physics with the dark matter sector have been pursued (see, e.g., [20]). ALPs as a dark matter candidate have recently received great attention due to the non detection of Weakly Interactive Massive Particles [21]. It is, therefore, important to examine the CIB data in light of the ALP hypothesis.
Apart from the increase of the EBL, this hypothesis has an observable impact on the propagation of TeV photons because it implies an enhanced opaqueness through γγ → e + e − processes. A higher level of EBL would help alleviate the tension between the observed neutrino spectrum and the gamma-ray spectrum of blazars, as discussed below. We present a case study in which multi-messenger, multi-wavelength observations can be exploited to obtain new tools to indirectly probe fundamental physics beyond the Standard Model, making use of data from neutrino telescopes (IceCube), gamma-ray satellites (Fermi-LAT) and sounding rockets equipped with infrared cameras (CIBER), extending the already flourishing multi-messenger astronomy tools [22].
The paper is organized as follows. We will review the CIBER data and the particle physics content of an ALP model; we will then tackle the bounds coming from anisotropy observations by the Hubble Space Telescope and CIBER itself. Later we show how the increase of the CIB affects the propagation of ultra-TeV gamma-rays. This brings us to a final discussion and our conclusions.

II. FLUX FROM AXION-LIKE PARTICLE DECAY
We are interested in the redshift evolution of the diffuse infrared radiation produced by the decay of a relic axion-like particle to a photon and a hidden photon, a → γ + χ [5,24]. The possibility of having axions with suppressed two-photons coupling has received some attention recently due to the peculiar phenomenology of photophobic axions [25]. The decay is due to the Chern-Simons [26] interaction Lagrangian whereF µν = µνρσ F µν , and the nonrelativistic decay rate for the ALP is found to be where the squared amplitude averaged over final polarization states is |M| 2 = (m 2 a −m 2 χ ); this correctly reduces to the usual axion decay rate when m χ = 0 and one includes a factor of 2 due to the final state involving identical photons [27]. Interestingly, the decay rate depends just on one kinematic quantity in the nonrelativistic approximation, namely, the maximum available energy for the outgoing photon ω max = (m 2 a − m 2 χ )/m a . The degeneracy would be broken if the ALP were non-negligibly relativistic.
The energy intensity (energy flux per unit of energy, time, surface per steradians) is read from where H(z) = H (0) Ω Λ + Ω m (1 + z) 3 is the Hubble function, z is the redshift at which the flux is "observed", z is the redshift at which a decays with a squared amplitude |M| 2 , the momentum at the production point is ω = ω(1 + z )/(1 + z) = ω (0) (1 + z ) (as well as p = p(1+z )/(1+z)), f a (p a ) is the momentum distribution of the ALPs, so that the number density (when there is no decay) is n a = d 3 p a /(2π) 3 f a (p a ). In the following the superscript (0) will indicate comoving quantities. We include the reduction in the number density due to decay with rate Γ over the cosmic time whereas we do not need to account for absorption; the latter is negligible in the wavelength range under study. The only relevant process reducing the flux of a single source is due to Thomson scattering [28]. However, Thomson scattering preserves the energy of the scattering photon. As such, it is irrelevant in the case of diffuse production with no sensible fluctuations in the electron spacial distribution, which we consider in first approximation to be homogeneous.
Equation (3) correctly reduces to Equation (50) of [29], when one takes a Cold Dark Matter (CDM) distribution for the ALP population, f a (p a ) = n (0) 3 , and gets rid of Dirac deltas. Integration over z yields is the Heaviside function, and we have defined the window function W (z , ω ). As expected, the comoving intensity is simply found by multiplying times a (1 + z) −4 factor (one power coming from ω). For z = 0 this agrees with Equation (3) of [5]. In the same paper, the (three-fold) parameter space to explain the CIBER excess has been explored. The maximum available energy must be ω max 10.2 eV to avoid constraints due to reionization and more stringently to the Lyman-alpha forest absorption spectrum.
The lifetime Γ −1 should be roughly of the order of the age of the Universe, and cannot be too small because ALPs can be produced in astrophysical systems, modifying the stellar evolution [5].
On the one hand, the ALP decay to a photon plus a hidden photon avoids the direct detection bounds on the coupling g aχγ , which instead constraint the g aγγ of standard ALPs (decaying to two photons), 1 as well as astrophysical bounds due to Horizontal Branch stars and SN1987a [30]; however, ALPs could still contribute to stellar cooling via plasmon decay γ → a + χ. This has been shown in [5] where only the transverse plasmon decay is considered. The longitudinal excitation of the electromagnetic field produce a longitudinal plasmon population [30]; the latter contributes negligibly to the cooling, because there is no resonant conversion from longitudinal plasmon to pseudoscalars [31]. Finally, there is another parameter which can be varied to fit CIBER data, the ALP number density n (0) a = R ρ DM /m a , where R/m a is a numerical factor and ρ DM is the total DM energy density. 1 In principle there should be also the operator L ⊃ gaγγ 4 aF µνF µν , but it can be technically natural to set gaγγ = 0 assuming a Z2 symmetry of which a and γ are different representations. This also sets to zero the kinetic mixing L ⊃ g kin F µν F χ µν , which would also contribute to stellar cooling.  [1,3]. The total flux (solid black) include the flux from ALP decay and the astrophysical diffuse (dotted black), which we assume to be the upper bound of the band reported in [23], In the following we will also consider a Warm Dark Matter (WDM) scenario, implying an additional fourth tunable parameter, namely the effective temperature. WDM can be produced both thermally or nonthermally [32,33]. In the first case, we suppose that the abundance is given by f a (p a ) = 1/[exp(|p a |/T th (z )) + 1], with the temperature T (0) th of the ALP population fixed to the temperature of relic neutrinos. The aforementioned distribution would be the result of relativistic decoupling, similar to what happens to neutrinos, and the decoupling should happen before e + e − annihilation. This allows us to avoid constraints coming from Cosmic Microwave Background measurements, as well as bounds due to the Big Bang Nucleosynthesis. There are other bounds to the existence of such a thermally produced ALPs, which indeed, being a contribution to radiation, affect temperature and polarization power spectra of the CMB anisotropies, the largescale matter power spectrum, and the Hubble expansion rate [34]. Measurements slightly favor a m a < 1 eV ALPs; however, m a ∼ 2 eV ALPs cannot be confidently excluded [35].
Alternatively, a nonthermally produced dark matter can have a momentum distribution with a strongly model dependent functional dependence [36]. One can assume the distribution to be f a (p a ) = R nth /[exp(|p a |/T nth (z )) + 1], where R nth is again a numerical factor and T (0) nth can be even higher than the CMB temperature. This distribution would be due to the decay of a Big Bang relics population (e.g. the inflaton). Equation (3) can be expressed in terms of special functions with these distribution. The Dirac delta function can be used to get rid of the angular part of the p a , and this would introduce a minimum absolute value of the momentum |p (0),min a |: where Li 2 is the polylogarithm of order 2, T = T depending on whether ω is smaller or bigger than ω max . The photon intensity spectrum due to the astrophysical diffuse, assumed to be [23], plus the ALP decay contribution is shown in Figure 1 for three different models. We plot two CDM scenarios with small and large ω max (model A and B), and a WDM scenario with large T nth to make more evident the ALP kinetic energy effect on the photon spectrum (model C nth ). Considering a thermally produced ALP population (in the following named model C th ), the intensity spectrum (not shown) is indistinguishable from the model A CDM spectrum. Model A and C th however differ strongly for what concerns the intensity anisotropies, as we will see below.

III. ANISOTROPIES
The gravitational clustering of dark matter makes the photon flux produced by the decaying ALP anisotropic. In this section we revisit the calculations as done in [5,37]. We take into account the energy resolution of the detector, following [38]. The average intensity of the flux detected in an energy band centered in ω with width ∆ω is given by assuming a ∆ω = ω flat passband filter for the detector [39]; the fluctuations toward a direction of the skyn can be expanded as spherical harmonics δI(ω, ∆ω,n) =I(ω, ∆ω,n) − I(ω, ∆ω) = l,m a l,m (ω, ∆ω)Y l,m (n) .
Our redshift dependence agrees with the one of Equation (A10) of [38], because we are considering the angular power spectrum of the energy flux (units are energy squared per time, per surface, per steradians and per energy); to compare the results of [38] and ours, Equation (A1) of the same reference shall be multiplied times ν, which gives an additional (1+z) −2 in the final expression. The anisotropy power spectra for lighter (ω max = 1 eV) and heavier (ω max = 8 eV) dark matter are shown in Figure 2, where they are compared with data of CIBER [44] and of the Hubble Space Telescope (HST) [45]. The matter power spectrum has been calculated through CLASS code [46], publicly available at [47]. In the first case, we explored both the CDM and the WDM cases (assuming m a = 2.2 eV for the latter case). The WDM power spectrum has been computed in the adiabatic approximation [48], P δ,WDM = (T WDM /T CDM ) 2 P δ,CDM , where T is the transfer function. The latter relates the primordial and the present day power spectra [49], and is another CLASS output [50]. In all cases, given that we needed to integrate over the redshift, we assumed conservatively a linear evolution for the matter power spectrum, using the non-linear matter power spectrum P δ obtained with CLASS, calculated at redshift z = 0 and evolved backwards Here, D(z) ∝ H(z) ∞ z dz (1 + z )H(z ) −3 is the linear growth factor, to be normalized with D(0) = 1 [42].
As heuristically expected, WDM (light red in the left panel) evades quite easily the constraints due to anisotropy measurements, as understood by showing the model C th anisotropy spectrum (left panel of Figure 2). These become unrestrictive when considering a nonthermally produced hot dark matter with high effective temperature, as in model C nth . Light CDM (model A) is just mildly excluded in our analysis.
For what concerns heavier dark matter (model B), our results are shown in the right panel of Figure 2, where the anisotropy power spectrum is computed both for the 1.6 µm wavelength band (light red) and for the 0.85 µm band (light blue). The 0.85 µm band slightly overshoots the observed data in the relevant wavelength (green points); however, the exclusion is much weaker than what has been found in previous analysis [5], due to averaging over the detector bandwidth.

IV. GAMMA-RAY ATTENUATION AND TIME VARIABILITY
The increased EBL flux has observable impact on the propagation of very high energy E > 0.1 TeV photons due to enhanced rate of e + e − pair production process. This effect may relax the tension between the predicted γ-ray flux and the Fermi LAT measurement of Isotropic Gamma-Ray Background (IGRB) [51] in traditional multi-messenger scenarios of high energy neutrino origin (see e.g. [52][53][54][55][56][57]) and eliminate need of hidden cosmic-ray accelerator [58]. In Fig. 3 we illustrate the effect. We calculate the neutrino and the accompanying γ-ray flux in the minimal pγ production scenario of Ref. [58] with b ν = 25 TeV, assuming low X-ray luminosity AGN evolution of Ref. [62] for the sources and the minimal EBL model [23] with or without the contribution from ALP, for which we use model A. 2 The spectra shown were obtained by solving transport equations for neutrinos and electron-photon cascades with the public numerical code [63]. The effect of the increased EBL is clearly seen on the γ-ray flux above 100 TeV.
At the same time the enhanced Universe opaqueness for γ-rays predicted in the above scenario will only sharpen the well known problem of unexpectedly hard γray spectra detected from the remote blazars. A possible solution to this problem proposed in Ref. [13][14][15][16][17][64][65][66][67][68] is based on natural assumption that the blazars also emit ultra-high energy cosmic rays which contribute to the observed γ-ray flux through secondary electro-magnetic cascades produced in line of sight cosmic ray interactions. The above scenario allows to avoid exponential γ-ray flux suppression with distance from the source.
Another important advantage of the scenario compared to primary γ-ray attenuation is its weak dependence on the EBL [64]. The mechanism proposed implies low strength of the magnetic field in voids B < 10 −14 G, which is consistent with observations [69][70][71], since otherwise the secondary cascade would be deflected. Furthermore due to expected time delays even magnetic field as low as B 10 −17 G may affect time variability of the signal on scale of hours [65], which imposes constrains on the secondary γ-ray flux contribution in case of variable γ-ray sources. To our knowledge, no variability has been observed for distant blazars in the energy range where secondary γ-ray flux is expected to dominate over primary γ rays.
The straightforward way to find if an extra component is needed to fit the observations is to construct the so called de-absorbed spectrum, i.e. the primary spectrum recovered from the observations assuming no extra components. The negative break in the de-absorbed spectrum can be considered as a good indication of extra component presence. By definition the de-absorbed spectrum F de−absorbed = exp(τ (z, E))F observed (15) depends not only on source redshift but also on the EBL model assumed through optical depth τ . We will illustrate this point on high-frequency peaked BL Lac object PG 1553+113, one of the most variable remote sub-TeV γ-ray sources known today. Its γ-ray flaring activity has been detected by H.E.S.S. telescopes during the nights of 2012 April 26 and 27 when the source flux above 0.3 TeV increased by a factor of 3 with evident signs of variability on scale of hours [61]. In Fig. 4 (left) we show the average spectrum of the object measured during the flare by Fermi LAT and H.E.S.S. (as shown in Fig. 3 of Ref. [61]) together with the de-absorbed spectra calculated using EBL model of Ref. [10] with or without extra contribution from ALP decay models A, B and C th(nth) . We use lower limit z > 0.43 [72] as a conservative source redshift estimate. It is now clear from the figure that increased EBL may lead to negative break in the de-absorbed spectrum, which indicates presence of extra component. Let us assume now that the extra component is not The γ-ray and neutrino fluxes expected in a minimal pγ production scenario of Ref. [58] (see details in text). Also shown are the per-flavor IceCube neutrino flux according to [59] (blue error bars) and more resent estimate [60] (green band).
The gamma-ray flux in the absence of ALP decays (dotted line) is decreased in the presence of an additional EBL component (solid line), which alleviates the tension with Fermi LAT IGRB measurements [51].  Fig.3 of Ref. [61]) together with the de-absorbed spectra calculated using EBL model of Ref. [10] with or without extra contribution from ALP. Right: Time delay distribution of secondary γ-rays from PG 1553+113 in optimistic case indicated in the text (the tail is not shown, the median value is indicated on the figure).
as highly variable as we would expect in the case of secondary γ from cosmic rays. Would it contradict observations? To answer this question in a conservative manner, we calculate the maximal expected integral flux of primary γ above 0.3 TeV during the flare phase F var max and the minimal required integral flux of the constant extra component F ext min . We calculate F var max assuming power low injection and maximal initial γ flux consistent with Fermi LAT observations below 30 GeV. F ext min is then calculated simply by subtraction of the primary component from the average observed flux at flare phase. For ALP models A, B and C we get F var max /F ext min integral flux ratio equal to 2.3, 0.36 and 7.6 respectively. From the observation that average integral flare flux above 0.3 TeV is 3 times higher than pre-flare flux we infer 3 = F var + F ext F const + F ext < F var + F ext F ext = F var F ext + 1, (16) where F const is possible contribution of primary photons in pre-flare flux. Now it is obvious that the condition F var /F ext < 0.36 which we have in case of model B contradicts (16), while other models are still in line with the inequality. Finally, we should note that although the average (median) time delays of the secondary γ-rays are expected to exceed days for reasonable magnetic field strengths B > 10 −17 G, the initial pulse can be preserved in part due to underdeveloped cascade particles produced locally within several tens of Mpc. In Fig. 4 (right) we show the result of secondary γ from PG 1553+113 time delay calculations made with Monte Carlo code [73]. Model B EBL was assumed, while the values of initial proton energy E = 40 EeV and turbulent magnetic field strength B = 10 −17 G close to lower limit [74] with correlation length L c = 1 Mpc were chosen to minimize the time delay. The delay scales with magnetic field strength and correlation length roughly as B × L 1/2 c .

V. CONCLUSIONS
In this paper we have explored the possibility that the high EBL spectrum detected by the CIBER collaboration could be due to the decay of an axion-like particle with mass around an electronvolt. Taking into account multi-messenger, multi-wavelength observations, we have shown that a warm dark matter component, produced either thermally or non-thermally, can explain the enhanced EBL detected by the sounding rocket CIBER. The increased level of EBL alleviates the tension between the neutrino flux detected at IceCube and the gammaray flux measured by Fermi, assuming a pγ production scenario. We have shown that the anisotropy measurements do not exclude this solution, and we have studied the effect on the propagation and time time variability of γ rays detected from distant sources, such as the Blazar Lac PG 1553+113. The ALP we consider is not in contradiction with current astrophysical observations, and the concordance of multi-messenger, multi-wavelength data lends credibility to the hypothesis that a decaying particle contributes to the measured excess of infrared background radiation.