Anisotropic cosmic optical background bound for decaying dark matter in light of the LORRI anomaly

Recently anomalous flux in the cosmic optical background (COB) has been reported by the New Horizon observations. The COB flux is $16.37\pm 1.47\, \rm nW m^{-2} sr^{-1}$, at the LORRI pivot wavelength of $0.608\,\rm \mu m$, which is $\sim 4\sigma$ level above the expected flux from the Hubble Space Telescope (HST) galaxy count. It would be great if this were a hint for the eV scale dark matter decaying into photons. In this paper, we point out that such a decaying dark matter model predicts a substantial amount of anisotropy in the COB flux, which is accurately measured by the HST. The data of the HST excludes the decay rate of the dominant cold dark matter larger than $10^{-24}$-$10^{-23}\,{\rm s}^{-1}$ in the mass range of $5$-$20$eV. As a result, the decaying cold dark matter explaining the COB excess is strongly disfavored by the anisotropy bound. We discuss some loopholes: e.g. warm/hot dark matter or two-step decay of the dark matter, to explain the COB excess.


Introduction
The origin of dark matter (DM) is one of the biggest mysteries in particle theory, astronomy, and cosmology. Many attempts are being made both from the theory side and experimental side. Recently the most precise measurement of the cosmic optical background (COB) was reported by the Long Range Reconnaissance Imager (LORRI) instrument on NASA's New Horizons mission [1,2]. The COB flux is measured to be 16.37 ± 1.47nWm −2 sr −1 , at the LORRI pivot wavelength of 0.608µm [3]. This is about ∼ 4σ level above the expected flux from the Hubble Space Telescope (HST) galaxy counts. It would be exciting if this is a hint of the DM [3,4].
The sub keV (single-component) DM should be either a spin-zero or spin-unity bosonic particle but should not be the fermion due to the Tremaine-Gunn bound [5,6]. In the spin-zero case, the QCD axion that solves the strong CP problem [7][8][9][10] in the hadronic axion window may be a good candidate [11,12]. Such a QCD axion as well as a more generic axion-like particle (ALP) can be produced non-thermally consistent with the cold DM paradigm [13][14][15][16][17]. 1 Alternatively, the hypothesis that the inflaton and DM are unified by a single ALP predicts the mass to be around eV [28,29] (see also Ref. [30]). In those cases, the sub-keV axions naturally decay into two photons with the photon couplings around g φγγ ∼ 10 −11 − 10 −10 GeV. See Refs. [31][32][33][34][35][36][37] for reviews of the axion and ALPs. In the ALP mass range of O(1-10) eV, the intensity of the cosmic optical background light was used to constrain the DM [32,38]. In this paper, we may not regard it as a constraint since the excess in the optical background has been found [3]. In such a situation, we may need to 1 Thermally produced QCD axions produced by pion interactions are too abundant to be consistent with the cosmic microwave background observations unless the axion is lighter than ∼ 0.5 eV [18][19][20][21][22]. However, this bound is significantly relaxed for low enough reheating temperature [23,24]. See also Refs. [25][26][27] for recent discussions about the theoretical estimation of thermal axion abundance.
consider an independent constraint to check whether or how the decaying DM can explain the LORRI anomaly.
In this paper we consider severe constraints on the DM scenario to explain the LORRI anomaly by using the COB anisotropy data [39]. Since the DM density spatially fluctuates in the universe, the photon flux from the DM decay not only contributes to the mean intensity but also to the anisotropy. A similar analysis has been made in Refs. [40,41] for the ALP model to explain the cosmic infrared background (CIB) mean intensity excess observed by the CIBER experiment [42], whose wavelength is longer than the COB measured by LORRI. 2 We apply the same idea to the DM model for the LORRI anomaly and derive constraints on such a scenario from the COB anisotropy measurements.
In this paper, we derive a robust bound from the anisotropic COB for simple cold DM models in the mass range of O(1-10) eV (Fig.2), and show the exclusion limit for the ALP DM in the mass-photon coupling plane (Fig.3). By taking account of the non-linear evolution of the density perturbation we found the bound from the data from HST [45] is so stringent that excludes all the parameter regions for the cold DM explanation of LORRI. We also discuss the possible loopholes and some more exotic DM models for explaining the LORRI anomaly.
2 Isotropic and anisotropic extragalactic background light We introduce a particle χ, which comprises a fraction R (≤ 1) of the total cold DM density, and it is assumed to have a decay mode into two particles including a photon γ: where x is a particle that may or may not be a photon. For simplicity we assume x is massless. We focus on the mass m χ of χ in the range 5 eV m χ 25 eV.
As shown in the Figure later in this paper, the lower bound comes from the optical measurement of a galaxy [46]. We set the upper bound on the mass because of the severe constraint from the reionization history [47,48]. As noted in the Introduction, χ (and x) is unlikely to be a fermion due to the so-called Tremaine-Gunn bound [5,6], m χ 0.5 keV, which is derived from the upper bound on the phase space density in dwarf spheroidal galaxies if it is dominant and if the DM does not have intrinsic multiplicity. A nice candidate may be a spin-zero/two boson which can decay into a pair of photons or a photon plus an exotic vector boson or a spin-one boson that can decay into a photon plus an exotic scalar boson. We assume the decay rate into photon is Γ. Since it comprises a fraction R of the cold DM, the averaged number density is expressed as n χ = ρ DM R(1 + z) 3 /m χ , with ρ DM being the measured present energy density of DM, z being the redshift. As we will see the parameters that are relevant in our analysis are as follows, as long as the χ lifetime is much longer than the age of Universe: • m χ which determines the wavelength of the resulting photon, m χ /(2(1 + z)).
The estimation of the photon flux does not depend on R, q γ and Γ independently, but depends on the combinationΓ. This means that our conclusions will also apply to the case that the dark matter is sub-dominant or/and decaying into a photon and a dark particle. However, it is not easily applied with q γ > 2 since the spectrum of the resulting photon will be quite different.
As a concrete example model, we can consider an ALP as the dominant DM. In this case, with g χγγ being the ALP photon coupling and q γ = 2.

Formalism
In this part we follow Refs. [40] and [41] for calculating the extragalactic background light (EBL) from decaying particle. As done in Ref. [41], we take the mean intensity of the flux detected at the energy ω with an observation bandwidth ∆ω, 3 Here where we have taken the speed of light to be unity. The Hubble parameter at the redshift z is given by 4 where Ω Λ , Ω m and Ω r denote the density parameter of the dark energy, total matter and radiation, respectively. The photon spectrum at the 2 body decay of χ has a delta-function shape 4 with E = (1 + z)E and ω max = m χ /2 in our massless approximation of the decay products.
Since we are focusing on χ as (a part of) the DM, which has a lifetime longer than the age of the Universe, we made the approximation e −Γt → 1 in the χ comoving number. The exponential neglected in Eq. (5) is included in the numerical estimation. The resulting isotropic COB energy flux is In order to explain the LORRI anomaly, we need m χ = 4-20 eV, withΓ ∼ 10 −23 -10 −22 s −1 [4]. The same setup also predicts the anisotropy of the photon flux since the DM density fluctuates in the universe [49,50]. To discuss the anisotropy, let us expand the angulardependent flux with spherical harmonics Y m (Ω), The relevant angular power spectrum is defined as In terms of W we obtain, with r(z) = z 0 dz/H(z) being the comoving distance, and j (kr(z)) the spherical Bessel function. The power spectrum of the matter density fluctuation δ is defined as The power spectrum will be discussed in detail later. When the power spectrum varies slowly with k the Limber approximation is applicable [49,51] 2 π dkk 2 P δ (k; r(z 1 ), r(z 2 ))j (kr(z 1 ))j (kr(z 2 )) Defining z max = ω max /(ω − ∆ω/2) − 1 and z min = ω max /(ω + ∆ω/2) − 1 as the maximum and minimum redshift observed in the anisotropy measurement, we have Note that the integral depends on the observation frequency ω. The observation bandwidth ∆ω depends on the experimental setup. The next task is to evaluate P δ (k, r, r) . To evaluate the power spectrum, we should take into account the non-linear structure formation effect. We include the one and two-halo contributions [52]: and For small (large) distance scales the dominant one is the one-halo (two-halo) contribution. For comparison with the COB anisotropy data discussed later, the one-halo term is dominant for the most region of the relevant multipole moment 10 3 10 6 . Since the estimation is complicated we list the various relevant functions and our strategy as follows.
• dn/dM denotes the comoving number density of halo with the mass of M , where ν = [δ c (z)/σ(M )] 2 , with the critical overdensity δ c (z) and σ(M ) is the variance of the linear density field in spheres containing a mean mass M , We adopt the Sheth-Tormen form for the multiplicity function f (ν): [53] where ν = 0.707 ν, p = 0.3 and A = 0.322, which is fixed from dνf (ν) = 1.
• ρ m is the averaged present (baryonic + dark) matter energy density.
• P (lin) δ (k) is the linear matter density perturbation for which we use the output of the public code Class [54] from the Planck best-fit cosmological parameters [55].
• u M (k) is the Fourier transform of the density profile of each halo [52]. For the Navarro-Frenk-White (NFW) density profile ρ dp (r) = ρ s r 3 s /r(r 2 + r 2 s ) [56], where the sine and cosine integrals are For given halo mass M and redshift z, the concentration parameter, c vir , and the parameters r s , ρ s are obtained as follows. First, let us define the Virial radius R vir as where [57] ∆ vir (z) = (18π 2 + 82y − 39y 2 ) with y =Ω m (z) − 1. Using the Virial radius, the concentration parameter is defined as For the NFW profile, r s = r −2 . The M, z dependence of the concentration function is model-dependent. For instance, in the power-law model [58][59][60] By using this fitting function, we will obtain r s . Finally, we obtain ρ s by the condition

Numerical result
So far, we see that a cold DM that explains the isotropic background flux necessarily induces an anisotropic one as well. Here we perform the numerical simulation to check how sizable the anisotropic contribution is. The M integration is performed in the range (10 −6 −10 17 )M for numerical calculation with M being the solar mass. We also take ∆ω = ω as assumed in Ref. [41]. This actually is a conservative choice (see the last paragraph of this section). The resulting angular power spectrum of the COB is shown in Fig.1 in -2 C /(2π) plane. We have taken λ obs = 0.85 µm 5 in the left panel and λ obs = 0.606 µm in the right panel. In each panel prediction from the decaying DM is shown for m χ = 10 eV and 15 eV with fixedΓ = 2 × 10 −23 s −1 . Also shown are the observed data points from the Hubble Space Telescope [45] with the error bars.
By requiring that the one-halo contribution in C (see Eq. (13)) does not exceed the upper error bar of any of the data points for λ obs = 0.85 µm and λ obs = 0.606 µm, we derived the upper bound onΓ as shown in Fig. 2. Here we adopt the power-law model, which predicts O(0.1) smaller C 2 from the previous analysis used in Fig. 1. The bound corresponds to the 95%CL exclusion limit by using a χ 2 distribution, the degrees of freedoms of which are chosen as the number of center values of data points that are smaller than the model predictions  at the correspondingΓ. 6 One can see that even in this case, the region explaining the COB excess by the LORRI is, unfortunately highly in tension with the COB anisotropy bound. We also present the reionization bound [32,38] and the indirect detection bound from the observations of galaxy clusters, VIMOS Abell 2667 and 2390 [61]. 7 In Fig. 3 we translate the constraint onΓ from the COB anisotropy measurement in Fig. 2 into the constraint on ALP-photon coupling g χγγ for the ALP dominant DM, by using Eq. (3). We also show the bound from the Horizontal branch star cooling for the photon coupling. We can see that the bound derived by us is more stringent than the cooling one.
Let us discuss uncertainties of the derived constraint.
• We used the Sheth-Tormen fitting formula for the halo mass function, which should be compared with the result of N-body simulation. For parameter ranges we are interested in, the relevant redshift is z 3.9 for m χ 20 eV for the observation frequency λ = 0.606 µm. For such a redshift, the difference between the Sheth-Tormen fitting formula and the N-body simulation coincide within the accuracy of about factor two [68]. Also we used the formula (25) for c vir . We checked that the formula given in Ref. [69] 6 We also checked that how we define statistics from the data is not very important in deriving the bound since the center values dominate over the error bars (see Fig.1). As we emphasized in the main text, the systematics are more important in deriving the bound, which we expect to strengthen the bound. 7 It is translated from the bound on g χγγ taken from the webpage https://cajohare.github.io/AxionLimits. The region that explains the LORRI anomaly is adopted from Ref. [4] Also shown are the translated bounds from reionization and the optical telescope [32,38,61].
gives almost the same constraint.
• As for the DM density profile, we employed the NFW one [56]. If we adopt other ones like the Burkert profile [70,71] (see also [72]), the resulting COB angular power spectrum may change only slightly at large . We have also checked that the dominant one-halo contribution changes by at most O(1) factor at 10 6 by taking a lower cutoff of M integration in Eq. (15) from 10 −6 M to 10 9 M . Since the actual lower cutoff is expected to be much smaller than this value due to the smallness of the DM free streaming length and Compton wavelength, essentially, it does not lead to any uncertainty unless we assume a non-standard DM production scenario so that the free streaming length becomes extraordinarily long. For clarification, in Fig. 4 we plot the relative contribution to the total one-halo power spectrum from logarithmic • For m χ 15 eV, the LORRI anomaly is explained by the low-frequency tail of the photon spectrum from the DM decay. Thus it predicts a much larger COB flux at a higher frequency. Such a case is highly disfavored by the TeV gamma-ray observation [73][74][75][76]. Although constraints from the TeV-gamma-ray observations are not shown in the Figure, one should note that the LORRI excess itself is slightly in tension with the TeV-gamma ray observations, and the tension becomes more and more serious for heavier axion. • We also comment on the choice of the bandwidth ∆ω. It should be noticed that, according to Eq. (13), C ∝ 1/∆ω roughly, and hence the expected flux would become larger for narrower bandwidth for the case of line photon spectrum from decaying DM. Although we have chosen ∆ω = ω for numerical calculation, the actual bandwidth may be narrower by a factor 2-3 for λ = 0.606 µm [77]. Thus the constraints that we have derived should be regarded as conservative ones in this respect. In other words, the anisotropy measurement of the photon flux has the potential to greatly improve the constraint or the detection probability of decaying DM into a photon with a line spectrum.

Conclusions and discussion
We derived upper bounds on the DM decay rate from the COB anisotropy in light of the excess in the COB mean intensity observed by LORRI. We found that the parameter region that explains the LORRI anomaly is highly disfavored by the COB anisotropy measurement. In other words, as long as χ plays the important role in the structure formation for the parametrically large scale, the derived bound is robust. Another remark is that the constraint we have derived should be regarded as a conservative one, as noted at the end of Sec. 2.2, due to the conservative choice of ∆ω. The constraint onΓ would become severer by a O(1) factor by precisely taking account of detector properties.
There are several loopholes to relax the COB anisotropy bound while keeping the mean intensity flux to explain the LORRI anomaly.
• Hot/warm χ: Let us suppose that χ has a relatively large velocity dispersion so that it is a hot/warm component rather than the cold one. If dominant, it is in tension with the structure formation and hence R 1 is required. In particular, if χ has a larger velocity than the escape velocity of a DM halo, the magnitude of the anisotropy on small scales is reduced [41]. 8 As a simple scenario, we may consider the case that χ is thermalized in the very early Universe. Assuming that χ is a real scalar, we obtain the corresponding effective number of neutrino species ∆N eff ∼ 0.027(106.75/g dec ) 4/3 , with g dec being the relativistic degrees of freedom at the decoupling of χ production. This scalar contributes to a fraction of the DM as a hot DM component, and the cosmological bound on such a hot relic reads m χ 10 eV for g dec = O(100).
• Two step decay of χ: Cold χ decays to lighter mediator particles φ i , and then the mediator decay into particles including a photon, In particular, we are interested in the case where the mediator has a long enough lifetime so that the decay length of φ i is longer than the typical size of the DM halo and the typical distance between the halos. Then the resulting photons are almost isotropic due to the randomized second decay vertex positions (see Ref. c.f. [82]). 9 In this sense, it may not even originate from the DM but from the decaying dark radiation from reheating [83,84] (see the case that the dark radiation is a decaying ALP [82]. The axion-photon conversion via cosmic magnetic field [85][86][87][88][89][90][91] may also be a candidate to explain the LORRI anomaly, although it requires further study since the resulting photon has certain anisotropy due to the magnetic field distribution.) In this case, one may confirm the evidence of the reheating by precisely determining the dark radiation spectrum [84] in the future line-intensity mapping experiments [92][93][94][95]. 8 Given the 5 σ Hubble tension, the bound on the hot DM may be different, c.f. the inflationary scenarios for explaining unusual primordial density perturbation required from various models for the Hubble tension [78,79] (See also [80,81] for the discussion to obtain n s ∼ 1 in order to alleviate the Hubble tension.) 9 A more exotic possibility may be the multi-component DM who are interacting with each other. In this case, the component that decays into particles, including a photon, may not follow the usual DM distribution, and hence the estimation of the anisotropy spectrum changes much.
However, in the first case, we need to enhance χ coupling to a photon to compensate for the small fraction, R, and hence the bound from the stellar cooling becomes severer. On the other hand, it would be possible to modify the model to alleviate the stellar cooling bound [40,41,83,[96][97][98]. To summarize, the simple decaying DM into 2-particles, including a photon, is difficult to explain the LORRI anomaly, but there might be several loopholes.