Cosmic-ray boosted inelastic dark matter from neutrino-emitting active galactic nuclei

TL;DR

This paper investigates how cosmic rays interacting with dark matter in active galactic nuclei can produce boosted dark matter particles detectable on Earth, providing new constraints on light and inelastic dark matter models.

Contribution

It calculates the flux of boosted dark matter from specific active galactic nuclei considering realistic cosmic-ray distributions and dark sector properties, deriving novel experimental bounds.

Findings

Cosmic-ray boosted dark matter can probe otherwise inaccessible parameter space.

New bounds on light and inelastic dark matter models from AGN sources.

Potential to test dark matter relic abundance scenarios.

Abstract

Cosmic rays may scatter off dark matter particles in active galactic nuclei, where both the densities of cosmic rays and dark matter are expected to be very large. These scatterings could yield a flux of boosted dark matter particles directly detectable on Earth, which enhances the sensitivity of dark matter direct detection and neutrino experiments to light and inelastic dark matter models. Here we calculate the cosmic-ray boosted dark matter flux from the neutrino-emitting active galactic nuclei, NGC 1068 and TXS 0506+056, by considering realistic cosmic-ray distributions, deep inelastic scatterings, and mass splittings in the dark sector. From this we derive novel bounds from these sources on light and/or inelastic dark matter models with Super-K and XENONnT. We find that cosmic-ray boosted dark matter from neutrino-emitting active galactic nuclei can test regions of parameter space favored to reproduce the observed relic abundance of dark matter in the Universe, and that are otherwise experimentally inaccessible.

Tables

Table 1. Table 1: Parameters considered in this work, needed to compute the CR proton spectra and DM distribution on NGC 1068 and TXS 0506+056. For NGC 1068 the corresponding values for the Case C2b are written within [.] if different from Case C2a. The minimum/maximum proton boost factor in the frame of the blob γ min / max , p ′ \gamma^{\prime}_{\min/\max,p} is provided, along with the spectral index of the cosmic rays in the frame of the blob , α p \alpha_{p} . The variable θ LOS \theta_{\rm{LOS}} is the angle between the jet and line-of-sight, here assumed to be 0 in all cases (i.e. either there is no jet, or it is directed at Earth). The radius of cosmic ray acceleration R a ​ c ​ c R_{acc} and the scale radius of the host galaxy r 0 r_{0} are both provided, which are important in calculating the integrated dark matter density through which these cosmic rays pass.

Parameter	TXS 0506+056	NGC 1068: C1	NGC 1068 : C2a [C2b]
$d_L\left(\mathrm{Mpc}\right)$	1750	10	10
$M_{\mathrm{BH}}\left(M_{\odot }\right)$	$5\times {10}^8$	$6\times {10}^6$	$6\times {10}^6$
${\Gamma }_B$	20	1	1
${\theta }_{\mathrm{LOS}}{(}^{\circ })$	0	0	0
${\alpha }_p$	2.0	3.2	3.2
${\gamma }_{\min ,p}^{\prime }$	1.0	${10}^4$	$1.4\times {10}^2$ [$3\times {10}^4$]
${\gamma }_{\max ,p}^{\prime }$	$1.6\times {10}^7$	$3\times {10}^4$	$3\times {10}^5$
$L_p\left(\mathrm{erg}/s\right)$	$8.6\times {10}^{49}$	$7\times {10}^{43}\left(\sim L_X\right)$	$7.8\times {10}^{44}\left(\sim L_{\mathrm{edd}}\right)$ [$1.2\times {10}^{42}$]
$R_{\mathrm{acc}}\left(R_S\right)$	$1.2\times {10}^5$	$30$	$30$
$r_0\left(\mathrm{kpc}\right)$	10	10	10

Equations

$$L_p$$

$$T_p\frac{d{\Gamma }_p}{d{T}_p}\approx \frac{4(1+K)}{3K}{f}_{p\gamma }{T}_{\nu }\frac{d{\Gamma }_{\nu }}{d{T}_{\nu }},$$

$$L_{\mathrm{Edd}}=\frac{4\pi Gc{m}_p}{{\sigma }_T}{M}_{\mathrm{BH}}\simeq 1.3\times 10^{45}\mathrm{erg}{\mathrm{s}}^{-1}\left(\frac{M_{\mathrm{BH}}}{10^7{M}_{\odot }}\right),$$

$$\frac{d{\Gamma }_p}{d{T}_{p}d\Omega }=\frac{1}{4\pi }{c}_p{\left(1+\frac{T_p}{m_p}\right)}^{-{\alpha }_p}\frac{{\beta }_p{\left(1-{\beta }_p{\beta }_B\mu \right)}^{-{\alpha }_p}{\Gamma }_B^{-{\alpha }_p}}{\sqrt{{\left(1-{\beta }_p{\beta }_B\mu \right)}^2-\left(1-{\beta }_p^2\right)\left(1-{\beta }_B^2\right)}},$$

$$L_p=\int d\Omega \int d{T}_p\left(T_p+m_p\right)\frac{d{\Gamma }_p}{d{T}_{p}d\Omega }=c_i{m}_p^2{\Gamma }_B^2{\int }_{{\gamma }_{\min ,p}^{\prime }}^{{\gamma }_{\max ,p}^{\prime }}d{\gamma }_p^{\prime }{\left({\gamma }_p^{\prime }\right)}^{1-{\alpha }_i},$$

$$c_p=\frac{L_p}{m_p^2{\Gamma }_B^2}\times \{\begin{array}{ll}\left(2-{\alpha }_p\right)/\left[{\left({\gamma }_{\max ,p}^{\prime }\right)}^{2-{\alpha }_p}-{\left({\gamma }_{\min ,p}^{\prime }\right)}^{2-{\alpha }_p}\right]& \text{ if }{\alpha }_p\ne 2;\\ 1/\log \left({\gamma }_{\max ,p}^{\prime }/{\gamma }_{\min ,p}^{\prime }\right)& \text{ if }{\alpha }_p=2,\end{array}$$

$${\rho }_{\mathrm{sp}}(r)={\rho }_R{g}_{\gamma }(r)(\frac{R_{\mathrm{sp}}}{r}{)}^{{\gamma }_{\mathrm{sp}}},$$

$${\int }_{4R_{\mathrm{S}}}^{R_0}4\pi r^2{\rho }_{\mathrm{sp}}(r)\mathrm{d}r=\Delta M_{\mathrm{BH}}.$$

$${\rho }_0={\left[\frac{\Delta M_{\mathrm{BH}}}{4\pi {\left({\alpha }_{\gamma }{r}_0{\left(M_{\mathrm{BH}}/r_0^3\right)}^{\frac{1}{3-\gamma }}\right)}^{{\gamma }_{\mathrm{sp}}-\gamma }{r}_0^{\gamma }{\int }_{4R_S}^{R_0}{g}_{\gamma }(r)r^{2-{\gamma }_{\mathrm{sp}}}dr}\right]}^{4-\gamma }$$

$${\rho }_0\simeq {\left[\frac{\left(3-{\gamma }_{\mathrm{sp}}\right)\Delta M_{\mathrm{BH}}}{4\pi {\left({\alpha }_{\gamma }{r}_0{\left(M_{\mathrm{BH}}/r_0^3\right)}^{\frac{1}{3-\gamma }}\right)}^{{\gamma }_{\mathrm{sp}}-\gamma }{r}_0^{\gamma }\left(R_0^{3-{\gamma }_{\mathrm{sp}}}-R_{\min }^{3-{\gamma }_{\mathrm{sp}}}\right)}\right]}^{4-\gamma }.$$

$$\frac{M_{\mathrm{BH}}}{10^8{\mathrm{M}}_{\odot }}\sim 0.7{\left(\frac{M_{\chi }}{10^{12}{\mathrm{M}}_{\odot }}\right)}^{4/3}.$$

$$\rho (r)=\{\begin{array}{ll}0& r\le 4R_{\mathrm{S}},\\ \frac{{\rho }_{\mathrm{sp}}(r){\rho }_{\mathrm{sat}}}{{\rho }_{\mathrm{sp}}(r)+{\rho }_{\mathrm{sat}}}& 4R_{\mathrm{S}}\le r\le R_{\mathrm{sp}},\\ {\rho }_0{\left(\frac{r}{r_0}\right)}^{-\gamma }{\left(1+\frac{r}{r_0}\right)}^{-(3-\gamma )}& r\ge R_{\mathrm{sp}}\end{array}.$$

$${\Sigma }_{\chi ,\mathrm{spike}}={\rho }_{\mathrm{sp}}\left(R_{\mathrm{acc}}\right)\frac{R_{\mathrm{acc}}}{g_{\gamma }\left(R_{\mathrm{acc}}\right)}{\int }_1^{R_{\mathrm{sp}}/R_{\mathrm{acc}}}dy{g}_{\gamma }\left(R_{\mathrm{acc}}y\right)y^{-{\gamma }_{\mathrm{sp}}},$$

$${\Sigma }_{\chi ,\mathrm{spike}}\simeq {\int }_{R_{\mathrm{acc}}}^{R_{\mathrm{sp}}}dr{\rho }_{\mathrm{sp}}(R_{\mathrm{acc}}){\left(\frac{r}{R_{\mathrm{acc}}}\right)}^{-{\gamma }_{\mathrm{sp}}}\simeq \frac{{\rho }_{\mathrm{sp}}(R_{\mathrm{acc}})R_{\mathrm{acc}}}{({\gamma }_{\mathrm{sp}}-1)}\left[1-{\left(\frac{R_{\mathrm{sp}}}{R_{\mathrm{acc}}}\right)}^{1-{\gamma }_{\mathrm{sp}}}\right].$$

$${\Sigma }_{\chi ,\mathrm{host}}={\rho }_0{r}_0\left[\ln \left(1+\frac{r_0}{R_{\mathrm{sp}}}\right)-\frac{1}{1+R_{\mathrm{sp}}/r_0}\right]\simeq {\rho }_0{r}_0[\log \left(\frac{r_0}{R_{\mathrm{sp}}}\right)-1],$$

$$\mathcal{L}\supset -\frac{1}{4}{F}_{\mu \nu }^{\prime }{F}^{\prime \mu \nu }+\frac{m_{Z^{\prime }}^2}{2}{Z}_{\mu }^{\prime }{Z}^{\prime \mu }+i{g}_{\chi }{Z}^{\prime \mu }{\bar{\chi }}_1{\gamma }_{\mu }{\chi }_2+(h.c)+i\sum _f{g}_f{Z}^{\prime \mu }\bar{f}{\gamma }_{\mu }f,$$

$$\frac{d{\sigma }_i}{d{T}_{{\chi }_2}}=\frac{{\sigma }_{\chi -i}}{4}\frac{m_{Z^{\prime }}^4}{(m_{Z^{\prime }}^2+q^2{)}^2}\frac{m_{\chi }\left[\right(s-(m_{\chi }^2+m_i^2+\delta m_{\chi }){)}^2+m_{\chi }{T}_{{\chi }_2}(q^2-2s)]}{2{\mu }_{\chi -i}^2\lambda (m_{\chi }^2,m_i^2,s)}{F}_i^2(q^2),$$

$${\sigma }_{\chi -i}=4\frac{g_i^2{g}_{\chi }^2{\mu }_{\chi -i}^2}{\pi m_{Z^{\prime }}^4},$$

$${\sigma }_{\chi -i}=\frac{64\pi {\mu }_{\chi -i}^2\alpha {\alpha }_D{\epsilon }^2}{m_{Z^{\prime }}^4},$$

$$F_p\left(q^2\right)=\frac{1}{{\left(1+q^2/{\Lambda }^2\right)}^2}$$

$$\cos \theta (T_p,T_{{\chi }_2})=\frac{-{\delta }^2+2m_p(\delta +T_{{\chi }_2})+2m_{\chi }{T}_{{\chi }_2}+2\delta T_p+2T_p{T}_{{\chi }_2}}{2\sqrt{T_p(2m_p+T_p)}\sqrt{T_{{\chi }_2}(2\delta +2m_{\chi }+T_{{\chi }_2})}}.$$

$${\chi }_2\to \{\begin{array}{l}{\chi }_1+Z^{\prime }\mathrm{if}\delta >m_{Z^{\prime }}\\ {\chi }_1+e+\bar{e}\mathrm{if}\delta <m_{Z^{\prime }}\&\delta >2m_e\\ {\chi }_1+\gamma \mathrm{if}\delta <m_{Z^{\prime }}\&\delta <2m_e\end{array}.$$

$$\cos {\theta }_p=\cos {\theta }_s\cos {\theta }_d+\sin {\theta }_s\sin {\theta }_d\cos {\varphi }_s.$$

$$E_{{\chi }_1}=\gamma E_{{\chi }_1,r}+\gamma \beta p_r\cos {\theta }_r,$$

$$\cos {\theta }_d=\frac{\gamma \beta E_{\chi 1,r}+\gamma p_r\cos {\theta }_r}{\sqrt{E_{{\chi }_1}^2-m_{\chi }^2}}.$$

$$\cos {\theta }_r=\frac{E_{{\chi }_1}-\gamma E_{{\chi }_1,r}}{\gamma \beta p_r}.$$

$$\frac{d{\Phi }_{{\chi }_1}}{d{E}_{{\chi }_1}}=\frac{{\Sigma }_{{\chi }_1}}{2\pi m_{\chi }{d}^2}\int d\cos {\theta }_{r}d{\varphi }_{s}d{T}_{{\chi }_2}d{T}_{p}d{E}_{{\chi }_1,r}\frac{d{\Gamma }_p}{d{T}_{p}d\Omega }({\theta }_p)\frac{d\sigma }{d{T}_{{\chi }_2}}\frac{1}{{\Gamma }_{\chi }}\frac{d^2{\Gamma }_{\chi }}{d\cos {\theta }_{r}d{E}_{{\chi }_1,r}}\delta (E_{{\chi }_1}-E_{{\chi }_1}(T_{{\chi }_2},{\theta }_r,E_{{\chi }_1,r})).$$

$$\begin{array}{rl}\frac{d{\Phi }_{{\chi }_1}}{d{E}_{{\chi }_1}}{|}_{iso}=& \frac{{\Sigma }_{{\chi }_1}}{4\pi m_{\chi }{d}^2}\int d{\varphi }_{s}d{T}_{{\chi }_2}d{T}_{p}d{E}_{{\chi }_1,r}\frac{d{\Gamma }_p}{d{T}_{p}d\Omega }({\theta }_p)\frac{d\sigma }{d{T}_{{\chi }_2}}\frac{1}{{\Gamma }_{\chi }}\frac{d{\Gamma }_{{\chi }_1}}{d{E}_{{\chi }_1,r}}\frac{1}{\gamma \beta p_r}\\ & \Theta (E_{{\chi }_1}+\gamma \beta p_r-\gamma E_{{\chi }_1,r})\Theta (\gamma E_{{\chi }_1,r}+\gamma \beta p_r-E_{{\chi }_1}).\end{array}$$

$$\frac{d{\Gamma }_{\chi }}{d{E}_{{\chi }_1,r}}{|}_{{\chi }_2\to {\chi }_1+Z^{\prime }/\gamma }={\Gamma }_{\chi }\delta (\frac{(m_{\chi }+\delta {)}^2+m_{\chi }^2-m_{Z^{\prime }/\gamma }^2}{2(m_{\chi }+\delta )}-E_{{\chi }_1,r}).$$

$$\frac{d{\Gamma }_{\chi }}{d{E}_{{\chi }_1,r}}{|}_{{\chi }_2\to {\chi }_1+e+\bar{e}}=-2(m_{\chi }+\delta )\frac{d{\Gamma }_{\chi }}{d{m}_{e\bar{e}}^2},$$

Full text

††preprint: FERMILAB-PUB-25-0620-T

I Introduction

Evidence for gravitational effects of non-luminous matter at different cosmological scales has been established through multiple observations [23, 34]. A well-motivated paradigm consists of dark matter (DM) being composed of new particles with weak or feeble interactions with the Standard Model (SM). A plethora of experiments have been devoted to search for DM particles from the Galactic halo via their elastic scatterings off nuclei and/or electrons at Earth-based detectors [61, 22, 45, 46]. As of today, no conclusive signal has been found, which allows us to set constraints of DM particles with masses from the MeV to the TeV scale [1, 11, 12, 109].

Constraints on spin-independent interactions in the mass range from 1 GeV to 1 TeV restricts several DM scenarios able to address the electroweak hierarchy problem while yielding the observed relic abundance of the Universe via freeze-out [75]. The constraints on MeV scale DM are still orders of magnitude weaker than for GeV-TeV scale DM, but several thermal and nonthermal production scenarios are already probed with some direct detection experiments, especially in light of recent results from the DAMIC-M and PANDAX-4T experiments [12, 80, 26, 33].

There is an important exception which is poorly tested by direct detection searches. In some DM models, the inelastic scattering channel can naturally dominate over the elastic one. A canonical example is the vector current of a Majorana particle, which is forbidden for the elastic case, but not for the inelastic one. This situation has a close analogue in neutrino physics. For Majorana neutrinos, the diagonal magnetic moment operator vanishes identically by CPT invariance [96, 77, 58], so the elastic channel $\nu_{i}\to\nu_{i}$ is absent. Only off-diagonal magnetic moments are allowed, which mediate inelastic processes of the form $\nu_{i}\to\nu_{j}$ with $i\neq j$. A similar mechanism appears in the Minimal Supersymmetric Standard Model, where the lightest supersymmetric particle can be dominantly Higgsino and much lighter than the other supersymmetric states. In that case, the elastic scattering channel is suppressed by the large masses of the supersymmetric partners, and the inelastic scattering channel induced by the electroweak gauge interactions can dominate [93].

Beyond these examples, several realizations of inelastic DM with mass splittings ranging from $\sim$ eV to $\sim$ MeV have been discussed in the literature, showing that simplified models can account for the DM of the Universe, and that a variety of phenomenological probes can constrain these models, e.g. in [66, 13, 101, 17, 85, 30, 18, 93, 92, 54, 31, 68, 21, 32, 55, 36, 65, 20]. Such inelastic models are interesting, as the scattering is suboptimal for direct detection probes.

It has been discussed in a variety of works that models of light and/or inelastic DM can be probed when considering scatterings of cosmic rays(CRs) off DM in the Milky Way, which would yield a boosted flux of DM particles on Earth that could be detected with DM and neutrino experiments, see e.g [27, 44, 28, 103, 83, 47, 19], and the scattering could also lead to a flux of secondary gamma-rays and neutrinos [64, 19]. These works were only able to probe relatively large scattering cross sections of DM off protons and electrons, coinciding with regions of parameter space that are mostly ruled out by several other probes. However, it was demonstrated in [105] that a larger boosted DM flux on Earth could instead arise from some blazars, due to the large CR proton flux and expected DM density in these environments.

In this work, we calculate the flux of CR boosted DM reaching the Earth from NGC 1068 (a Type-II Seyfert galaxy) and TXS 0506+056 (a blazar), two well known multi-messenger sources and CR accelerators. We derive upper limits on a simplified model of DM from a non-observation of an excess of events induced by DM interactions with nuclei and electrons at Earth-based experiment Super-Kamiokande (Super-K) [10] (and XENONnT [16] in Appendix E). We discuss the complementarity of our constraints with those obtained from other searches, demonstrating that active galactic nuclei (AGN) allow us to improve over previous bounds in a large region of parameter space. Previous works considering scatterings of DM particles off CRs in AGN focused either on spectral distortions on the observed gamma-ray fluxes [62, 14], CR electron scatterings off DM [108, 67], the implications of CR-DM interactions for multi-messenger observations including high-energy neutrinos [69, 65, 38, 40, 106, 39], or studied the boosted DM flux considering only the multi-messenger blazar TXS 0506+056 [105, 63, 24, 74, 38] or considering additional blazars that do not present evidence for CR proton acceleration [39, 39]. These works generally present predictions without addressing CR model uncertainties, and did not consider inelastic DM scenarios. In this work we address all these tasks in detail for TXS 0506+056, further computing the boosted DM flux from the steady neutrino-emitting source NGC 1068.

The paper is organized as follows. In Section II we describe our modeled CR flux from NGC 1068 and TXS 0506+056, discussing the associated uncertainties. In Section III we describe the DM density profile and column density in NGC 1068 and TXS 0506+056. In Section IV we compute the CR boosted DM flux from NGC 1068 and TXS 0506+056, considering deep inelastic scattering (DIS) and mass splittings in the dark sector. We also comment on the attenuation of a DM flux as it passes through Earth. In Section V we discuss the scattering signatures of this flux at Super-K, and derive upper limits on the DM parameters. We further confront these limits with complementary constraints and motivated regions of parameter space able to reproduce the relic abundance of the Universe. Finally, in Section VI, we present our conclusions.

II Cosmic-ray flux in active galactic nuclei

Recent high-energy neutrino observations from TXS 0506+056 [5, 4, 7] and NGC 1068 [6], together with the EM counterparts, which are produced via $pp$, $p\gamma$ interactions, plus leptonic processes such as inverse compton scattering or synchroton radiation, allow us to infer the energy budget of CRs. The so-called leptohadronic models have been able to explain the neutrino and EM spectra simultaneously under certain conditions, and the accelerating region of CRs has been constrained to some extent [29, 91]. Below we discuss the details of the CR proton spectra in these sources.

II.1 TXS 0506+056

TXS 0506+056 is a blazar (an AGN with a relativistic jet pointing towards the observer) at redshift $z\approx 0.33$ ($d_{L}=1750$ Mpc) from which a neutrino event (IceCube-170922A) of energy $0.1-1$ PeV was reported [4] at the level of $3\sigma$ while there was an ongoing gamma-ray flare. The neutrino observation was followed up by various electromagnetic (EM) observations in the high-energy and very high-energy gamma-rays, hard and soft X-rays, ultraviolet, optical, and near-infrared wavelengths [78] (see references therein). In addition to this, an archival search by the IceCube Collaboration at the location of TXS 0506+056 lead to significant association of another $13\pm 5$ neutrino events [5]. The presence of such high energy signatures both in the EM and neutrino channels make TXS 0506+056 an ideal site of particle acceleration [90].

The spectral energy distribution of TXS 0506+056 from optical to gamma-rays can be explained by accelerated primary electrons that undergo inverse Compton and synchrotron processes. In fact, to simultaneously explain the high-energy neutrino and gamma-ray observations such leptonic models are amongst the only viable ones, where the radiation from accelerated protons are hidden due to EM cascades in the source and $\gamma\gamma$-attenuation. The neutrino and proton luminosities are thus strongly constrained by X-ray observations, else the high energy photons produced from photohadronic interactions would violate the X-ray data.

The differential CR proton luminosity is related to the total luminosity by

[TABLE]

where $T_{p}$ is the CR proton energy kinetic in the source frame, $\Gamma_{p}$ is the production rate of protons, $L_{p}$ is the normalization, and $T_{p}^{\rm min}$ and $T_{p}^{\rm max}$ give the minimum and maximum CR kinetic energies respectively. In the case with no jet (that is $\Gamma_{B}=1$), then $T_{p}^{\min/\max}=m_{p}\gamma^{\prime}_{\min/\max,p}$, where $\gamma^{\prime}_{\min/\max,p}$ is the minimum/maximum boost factor in the frame of the blob. In the case with a jet, $\gamma^{\prime}_{\min/\max,p}$ are still used as cutoffs in the blob frame, and the velocity distribution is boosted to obtain $T^{\min/\max}_{p}$ in the lab frame. In Fig. 1 we show the CR proton spectra in terms of $T_{p}^{2}\frac{d\Gamma_{p}}{dT_{p}}$. We show the results from the IceCube Collaboration [5] in blue. For TXS 0506+056, based on IceCube observations, the CR protons are accelerated from 0.2 PeV to 2 PeV. The best fit line along with the upper and lower bounds are obtained by converting the neutrino observations assuming $T_{p}\approx 20T_{\nu}$. The differential CR proton spectrum is approximately related to the neutrino spectrum by [90, 91]

[TABLE]

where we have $K=1(2)$ $p\gamma$ ($pp$) processes, $f_{p\gamma}$ is defined as the efficiency of pion production, that is, the fraction of CR protons producing pions through the photomeson process. Observations of gamma-rays between 10 - 100 GeV imply the inefficiency of neutrino production in the same region. Thus setting the optical depth to be less than 1 at 100 GeV implies $f_{p\gamma}\sim 10^{-4}$ [90] (see Sec. 3 there). The Eddington luminosity assuming a black hole mass $M_{\rm BH}\sim 5\times 10^{8}M_{\odot}$ is shown with a black dotted line111The Eddington luminosity ($L_{\rm Edd}$) is obtained by balancing the outward radiation pressure and the inward gravitational force. In hydrostatic equilibrium

$L_{\rm Edd}=\frac{4\pi Gc\,m_{p}}{\sigma_{T}}\,M_{\rm BH}\simeq 1.3\times 10^{45}\ {\rm erg\,s^{-1}}\left(\frac{M_{\rm BH}}{10^{7}M_{\odot}}\right),$

(3)

where $M_{\rm BH}$ is the mass of the black hole and $\sigma_{T}=6.65\times 10^{-25}\ {\rm cm}^{2}$ is the Thomson cross section. . In this case it is evident that $L_{p}$ for single zone models we consider in this work violate $L_{\rm Edd}$ which is reasonable given that the absolute jet luminosity for blazars can exceed the accretion luminosity (generally a fraction of the Eddington luminosity) [56]. The CR protons are isotropically accelerated within a blob, and the blob itself moves in the jet of TXS 0506+056 with speed $\beta_{B}$ with respect the central black hole. The corresponding Lorentz boost factor is given by $\Gamma_{B}\equiv\left(1-\beta_{B}^{2}\right)^{-1/2}$. For TXS 0506+056 we assume $\Gamma_{B}\approx 20$.

We follow the model ”LMBB2b” from [78] to simultaneously explain the neutrino and EM observations, where the CR proton power-law index $\alpha_{p}=2$. The Lorentz factor $\Gamma_{B}\approx 20$ which results in an isotropic equivalent CR proton luminosity (in the observer frame) $L_{p}=\Gamma_{B}^{4}L_{p}^{\prime}=8.6\times 10^{49}\ \rm erg\ s^{-1}$. Based on the model, the neutrino emission peak occurs at 10 PeV which corresponds to $\gamma_{p,\rm max}^{\prime}=1.6\times 10^{7}$, where primed coordinates refer to boosts in the frame of the blob. As for the low energy cut we assume $\gamma_{p,\rm min}^{\prime}=1$. The comoving blob radius is given by $R_{B}^{\prime}\approx\Gamma_{B}ct_{\rm var}/(1+z)\sim 300R_{s}$, where the Schwarzschild radius $R_{s}=2GM_{\rm BH}/c^{2}$, $t_{\rm var}$ is the variability timescale of the source. Thus, the acceleration radius $R_{\rm acc}\sim\Gamma_{B}R_{B}/{\rm tan}\ \theta_{j}$, where $\theta_{j}\sim 1/\Gamma_{B}$ which gives $R_{\rm acc}\sim\Gamma_{B}^{3}ct_{\rm var}/(1+z)\approx 1.8\times 10^{19}$ cm ($\sim 6$ pc). Based on the X-ray and gamma-ray observations the variability timescale is $t_{\rm var}\lesssim 10^{5}$ s. The CR proton spectra is shown in Fig. 1 (left) as a solid blue line.

II.2 NGC 1068

The IceCube Collaboration also reported an excess of $79^{+22}_{-20}$ neutrinos at $4.2\sigma$ significance associated with NGC 1068 - a Seyfert Type-II galaxy at a distance of $\sim 10$ Mpc [6]. The galaxy consists of a central black hole with mass $M_{\rm BH}=6\times 10^{6}M_{\odot}$ along with an accretion disk. The production of high-energy neutrinos and consequently gamma-rays is mostly attributed to hadronic processes - hadronuclear ($pp$) and photomeson ($p\gamma$) interactions, resulting in the production of mesons. The observed isotropic equivalent neutrino luminosity is $L_{\nu}\simeq 1.4\times 10^{42}\ \rm erg\ s^{-1}$ for neutrino energy $E_{\nu}$ between 1.5 TeV to 15 TeV. The isotropic equivalent gamma-ray luminosity between 100 MeV and 100 GeV is $7.7\times 10^{40}\ \rm erg\ s^{-1}$. Such high-energy emissions naturally require sites of efficient particle acceleration. The various mechanisms for such particle acceleration include magnetically-powered corona models via turbulence and magnetic reconnection, shocks, and reconnection-driven flows [89, 91, 51, 37, 50].

The IceCube neutrino observations from [6] (see Fig. 4 there) are converted into corresponding CR proton energies by assuming $T_{p}\approx 20T_{\nu}$ as before. The differential CR proton spectrum is approximately related to the neutrino spectrum by Eq. (2). For a conservative estimate we take $K=1$ corresponding to the $p\gamma$ scenario, where the system is calorimetric such that $f_{p\gamma}\sim 1$ [89, 91, 37]. This gives the best fit range (blue shared region) in the right plot of Fig. 1, while the light blue band between the dashed blue lines cover the power law CR proton spectra that are consistent with the neutrino observations at $95\%$ C.L. The best-fit spectral index for the CR protons is given by $\alpha_{p}=3.2$ (see Eq. 1). Given this index, the acceleration radius ($R_{\rm acc}$) is $\sim 30\ R_{s}$. This is consistent with the requirements of the minimum CR luminosity [91, 37, 89] in case of both the $pp$ and $p\gamma$ scenarios. Note that the high-energy neutrino production efficiency is inversely proportional to the acceleration radius $R_{\rm acc}$, thus a smaller radius leads to more efficient neutrino production with a smaller value of $L_{p}$. We explore two cases corresponding to the CR spectra for NGC 1068 which are distinguished by the normalization of the CR luminosity ($L_{p}$), the minimum ($T_{p}^{\rm min}$), and maximum ($T_{p}^{\rm max}$) CR proton kinetic energy.

• •

Case 1 (C1): X-ray observations for the total coronal luminosity above $2$ keV is given by $L_{X}=7^{+7}_{-4}\times 10^{43}{\rm erg\ s}^{-1}$. We adopt this as a representative case where we normalize the CR luminosity (see Eq. 1) to $L_{X}$. Note that this is still conservative given that the contributions from the optically thick and geometrically thin disk is such that the bolometric luminosity is $L_{\rm bol}\simeq 4.8\times 10^{44}{\rm erg\ s}^{-1}$. The power law index ($\alpha_{p}$) is consistent with the IceCube best fit. We choose $T_{p}^{\min}=10$ TeV and $T_{p}^{\rm max}=30$ TeV. Given that the IceCube neutrino flux has a lower limit of $E_{\nu}\sim 1$ TeV (implying $T_{p}\sim 20$ TeV), our chosen value of $T_{p}^{\rm min}$ is reasonable. On the other hand, the choice of $T_{p}^{\rm max}$ is conservative since a larger value would imply more neutrinos with higher energies leading to luminosities greater than $L_{X}$. This case is shown as the solid red line in Fig. 1.

• •

Case 2a (C2a): For this case, we choose $T_{p}^{\rm max}\sim 300$ TeV according to the upper limit from the observed IceCube neutrino flux. The lower limit is chosen such that the normalization $L_{p}\approx L_{\rm edd}$, where the Eddington luminosity is given by $L_{\rm edd}\simeq 7.7\times 10^{44}{\rm erg\ s}^{-1}$ and shown as the black dotted line in Fig. 1. This results in $T_{p}^{\rm min}\sim 100$ GeV. The spectral index for the CR protons match the IceCube best fit value of $3.2$. This can in-principle mimic the $pp$ scenario for neutrino production. We show this case using a purple dashed line in Fig. 1. Note that choosing a harder spectrum (resulting from magnetic reconnection and/or stochastic acceleration) like $\alpha_{p}=2$, would in turn lead to a larger value of $T_{p}^{\rm min}$ since $L_{p}\lesssim L_{\rm edd}$.

• •

Case 2b (C2b): We also consider a small variation of the above C2a case where instead of extrapolating the IceCube best fit to lower energies, we just use the IceCube results for the best fit (solid blue line) between $E_{p}^{\rm min}\sim 30$ TeV and $T_{p}^{\rm max}\sim 300$ TeV (same as C2a). The normalization for this case is $L_{p}\simeq 1.2\times 10^{42}{\rm erg\ s}^{-1}$. This serves as the most conservative scenario.

We discuss the results corresponding to C1 in the main text, while the cases C2a and C2b are discussed in Appendix B. The two cases described above complement each other and attempts to fold in the model uncertainties associated with the astrophysical modeling of NGC 1068, including, sites and mechanisms of particle acceleration, neutrino production processes and the resulting spectra. The two cases thus help explore the non-trivial effects on the constraints as will be evident in the following sections.

II.3 Spectral Modeling

The spectral distribution of CRs on Earth’s frame corresponds to a power-law function with energy. We follow the formalism from [105, 63], such that the CR flux reads

[TABLE]

If we want $\frac{d\Gamma_{p}}{dT_{p}}$ as used in Eq. 1, then we integrate over solid angle $\Omega$. In the above equation the subscript $p$ is for protons with mass $m_{p}\simeq 0.938$ GeV, $\alpha_{p}$ is the spectral index of the CR proton flux, $T_{p}$ and $\beta_{p}=$ $\left[1-m_{p}^{2}/\left(T_{p}+m_{p}\right)^{2}\right]^{1/2}$ are respectively the kinetic energy and velocity of the particle, $\mu$ is the cosine of the angle between the particle’s direction of motion and the jet, and $c_{i}$ is a normalization constant proportional to the source luminosity $L_{i}$, via

[TABLE]

thus

[TABLE]

with $\gamma^{\prime}$ the particle Lorentz factor in the frame of the blob. In Table 1, we show the relevant parameters considered in this paper for TXS 0506+056 and NGC 1068. Since NGC 1068 does not show evidence for a jet, we can set $\Gamma_{B}=1,\,(\beta_{B}=0)$ and use the same equations. We use these parameters to show the results in Fig. 1. We note that at high energies, our spectrum follows a power-law, while for low energies at TXS 0506+056, the flux becomes suppressed. This is due to the jet boost factor: even protons which are low energy in the frame of the blob are boosted to relativistic energies. This effect becomes strong at $T_{p}\lesssim\Gamma_{B}m_{p}.$

III Dark matter distribution in Active Galactic Nuclei

Now we turn towards discussing the density profile of DM particles at the source and the corresponding column density that CRs traverse before escaping the source, while boosting DM on its way. It has long been known that adiabatically-growing black holes are expected to form a spike of DM particles in their vicinity [98, 99, 60]. An initial DM profile of the form $\rho(r)=\rho_{0}(r/r_{0})^{-\gamma}$ evolves into

[TABLE]

where $R_{\rm sp}=\alpha_{\gamma}r_{0}(M_{\rm BH}/(\rho_{0}r_{0}^{3}))^{\frac{1}{3-\gamma}}$ is the size of the spike, with $\alpha_{\gamma}\simeq 0.293\gamma^{4/9}$ for $\gamma\ll 1$, and numerical values for other values of $\gamma$ provided in [60]. The steepness of the profile is parameterized by222The steepness of the spike is robustly predicted from dimensional arguments. Let’s assume a model consisting entirely of circular orbits. Further, let’s assume that the initial density cusp is $\rho\sim r^{-\gamma}$ before the addition of the black hole and $\rho\sim r^{-\gamma_{\mathrm{sp}}}$ after the black hole has accreted most of its mass. Conservation of mass implies that $\rho_{i}r_{i}^{2}dr_{i}=\rho_{f}r_{f}^{2}dr_{f}\Longrightarrow r_{i}^{3-\gamma}\sim r_{f}^{3-\gamma_{\mathrm{sp}}}$. Conservation of angular momentum implies that $r_{i}M_{i}(r)=r_{f}M_{f}(r)\simeq r_{f}M_{\rm BH}\Longrightarrow r_{i}^{4-\gamma}\sim r_{f}$. Combining these two results, one gets $\gamma_{\mathrm{sp}}=\frac{9-2\gamma}{4-\gamma}$. $\gamma_{\rm sp}=\frac{9-2\gamma}{4-\gamma}$, and $r_{0}$ denotes the scale radius of the host galaxy. Furthermore, $g_{\gamma}(r)$ is a function which can be approximated for $0<\gamma<2$ by $g_{\gamma}(r)\simeq(1-\frac{2R_{\rm S}}{r})$, with $R_{\rm S}$ being the Schwarzschild radius, while $\rho_{R}=\rho_{0}\,(R_{\rm sp}/r_{0})^{-\gamma}$.

We consider that the initial DM distribution follows a NFW profile [95, 94] with $\gamma=1$, resulting in a spike with $\gamma_{\rm sp}=7/3$ and $\alpha_{\gamma}=0.122$ [60] (although we note that if the dark matter distribution near the galactic center is a core as suggested in [100], this would only mildly change the spike index to 2.25). The masses of the central SMBHs of the two AGN considered in this work are given in Table 1. The normalization factor ($\rho_{0}$) can be obtained from the uncertainty in the black hole mass [62, 82], in such a way that the profile is compatible with both the total mass of the galaxy and the mass enclosed within the radius of influence of the $\mathrm{BH}$, of order $10^{5}R_{\mathrm{S}}$. We follow in this paper the criteria from [62]. The DM spike mass within the region that is relevant for the determination of the BH mass, typically $R_{0}=10^{5}R_{\mathrm{S}}$, must be smaller than the uncertainty on the $\mathrm{BH}$ mass $\Delta M_{\mathrm{BH}}$. For NGC 1068, we adopt $M^{\rm NGC}_{\rm BH}=(1-2)\times 10^{7}M_{\odot}$ [107], and for TXS 0506+056, we take $M^{\rm TXS}_{\rm BH}=(3-10)\times 10^{8}M_{\odot}$ [97]. The normalization constant $\rho_{0}$ is thus obtained by solving the following equation

[TABLE]

Factorizing $\rho_{0}$ in this expression we find a general expression

[TABLE]

where we have used the fact that $\frac{\gamma-3}{\gamma_{sp}-3}=4-\gamma$ to simplify the overall exponent. This expression still contains a non-trivial integral over $r$. It has previously been assumed in the literature that $g_{\gamma}(r)\sim 1$ (only true for $r\gg R_{S}$), and considered that the mass is dominated by the contribution from $r\gg R_{\mathrm{S}}$, i.e., typically $r>R_{\min}=\mathcal{O}\left(100R_{\mathrm{S}}\right)$ [62, 81]. One can then obtain a simpler expression

[TABLE]

This criteria, applied to $\gamma=1$, yields masses of the DM halo below as expected from universal relations between supermassive black hole (SMBH)-galaxy bulge masses [42, 48]:

[TABLE]

If the DM particles self-annihilate, the maximal DM density in the inner regions of the spike is saturated to $\rho_{\text{sat}}=m_{\chi}/(\langle\sigma v\rangle t_{\mathrm{BH}})$, where $\langle\sigma v\rangle$ is the velocity averaged DM self-annihilation cross section, and $t_{\rm BH}$ is the SMBH age. We assume in this work that (co-)annihilations are not efficient in depleting the DM spike and $\rho_{\rm sat}\gg\rho_{\rm sp}$. Furthermore, the DM spike extends to a maximal radius $R_{\rm sp}$, beyond which the DM distribution follows the initial NFW profile. The DM density profile therefore reads [60, 81, 82]

[TABLE]

From the DM distribution we can find the column density of DM particles at these sources analytically. For the accelerating region of CRs, $R_{\rm acc}$, we take $R^{\rm NGC}_{\rm acc}=30R_{S}^{\rm NGC}$ [91] and $R^{\rm TXS}_{\rm acc}=1.2\times 10^{5}R_{S}^{\rm TXS}$ [78] (see Table 1), consistent with inferred values from multi-messenger observations from these sources , and with our CR modeling descriptions in Sec. II. The exact column density in the spike reads

[TABLE]

with $y=r/R_{\mathrm{acc}}$. Dropping the $g_{\gamma}(r)$-dependence, one gets the approximation [49]

[TABLE]

Further, the contribution to $\Sigma_{\chi}$ from the passage through the halo of the host galaxy is

[TABLE]

where the last approximation assumes $r_{0}\gg R_{\rm sp}$. Under this prescription, we find that the column density of DM particles in the spike (with $\gamma=1$) of TXS 0506+056 is given by $\Sigma^{\rm TXS}_{\chi}=1.9\times 10^{28}$ GeV cm*-2*, and for NGC 1068, we find $\Sigma^{\rm NGC}_{\chi}=1.7\times 10^{33}$ GeV cm*-2*. We will use these fiducial values along the paper. The corresponding boosted fluxes would change linearly with the column density of DM traversed at the source.

IV Cosmic-ray boosted dark matter flux

We consider a coupling between the SM and the dark sector of the form

[TABLE]

with $F^{\prime\mu\nu}=\partial^{\mu}Z^{\prime\nu}-\partial^{\nu}Z^{\prime\mu}$ and the sum in the last term being over SM fermions. In what follows we will examine our constraints on two classes of models. In the first case, we adopt the assumption that the $Z^{\prime}$ couples equally to protons, neutrons, and electrons, $g_{e}=g_{p}=g_{n}=3g_{q}$, and set other SM couplings to zero. In the second case, we will examine the dark photon scenario in which $g_{e}=-g_{p}$, with all other SM couplings set to zero. Given that the neutron coupling plays a minor role in the phenomenology we consider, both model predictions can be done simultaneously. The dark sector contains a vector mediator $Z^{\prime}$, a stable Dirac fermion DM particle $\chi_{1}$ of mass $m_{\chi}$, and an excited Dirac fermion $\chi_{2}$ of mass $m_{\chi}+\delta$. The DM ($\chi_{1}$) scattering cross section with CRs $i=p,e$ ($\chi_{1}+i\rightarrow\chi_{2}+i$) is given by

[TABLE]

where $T_{\chi_{2}}$ is the upscattered DM kinetic energy, the momentum transfer is $q^{2}=2m_{\chi}T_{\chi_{2}}-\delta^{2}$, the Mandelstam variable $\mathrm{s}=m_{i}^{2}+m_{\chi}^{2}+2\left(m_{i}+T_{i}\right)m_{\chi}$, $\lambda(a,b,c)=a^{2}+b^{2}+c^{2}-2ab-2ac-2bc$. The non-relativistic DM-proton and DM-electron scattering cross section is given by

[TABLE]

where $g_{i}$ denotes the gauge coupling of the mediator $Z^{\prime}$ to protons or electrons, and $g_{\chi}$ denotes the coupling of the mediator to the DM. The reduced mass of the DM-proton/electron system is $\mu_{\chi-i}$. In models with a vector portal between the DM and the SM sectors, the gauge coupling $g_{i}$ is typically proportional to the kinetic mixing $\epsilon$ between the vector boson and the SM photon [70]. The non-relativistic cross section from Eq. 18 can then be written as

[TABLE]

where $\alpha_{D}=g_{\chi}^{2}/4\pi$ is the dark electromagnetic structure constant. The from factor $F_{i}\left(q^{2}\right)$ can be for either the electron, quark, or proton. For electrons and quarks, it equals 1, while for protons it reads

[TABLE]

where $\Lambda\simeq 770$ MeV [15]. We note that while this form factor is simple, any corrections will enter at $\mathcal{O}(q^{2}/m_{p}^{2})$ or higher [39]. The form factor above already suppresses scattering at these momentum transfers, so any corrections will not change the overall phenomenology, at most altering behavior by $\mathcal{O}(1)$ in a region where our constraints are weak. For now, we will consider only elastic proton scattering, although we will comment on DIS later on (resonances are ignored entirely, which is a mildly conservative approximation, since either elastic scattering or DIS dominate over most energy ranges considered in this work).

We know that the initial scattering angle is uniquely determined by the proton and $\chi_{2}$ kinetic energies

[TABLE]

IV.1 Decays

We now introduce DM decays, and assume the decay length is small compared to intergalactic distances, so all interactions can be assumed to happen at the AGN. We will consider the following two-body and three-body decay processes

[TABLE]

In order to account for such decay processes, let us consider a proton which is propagating at angle $\theta_{p}$ relative to the jet (where the jet iself is pointed at Earth). The proton upscatters a $\chi_{2}$ particle at scattering angle $\theta_{s}$ and azimuthal angle $\phi_{s}$. The $\chi_{2}$ then decays into a $\chi_{1}$ which propagates at angle $\theta_{d}$ relative to the parent particle. For the $\chi_{1}$ particle to reach Earth, the geometry requires

[TABLE]

As indicated in Eq. 21 the scattering angle $\theta_{s}$ is uniquely determined by $T_{p}$ and $T_{\chi_{2}}$. When considering decays, it is helpful to work in the rest frame of the $\chi_{2}$ particle, so we can define the angle, energy, and momentum of $\chi_{1}$ in this frame to be $\theta_{r}$, $E_{\chi_{1},r}$, $p_{r}$ respectively. Letting $\gamma=1+T_{\chi_{2}}/(m_{\chi}+\delta)$ and $\beta=\sqrt{1-1/\gamma^{2}}$, we can find quantities in the lab frame as

[TABLE]

and

[TABLE]

Inverting these expressions we find

[TABLE]

By considering the double differential decay rate of $\chi_{2}$, $\frac{d^{2}\Gamma_{\chi}}{d\cos\theta_{r}dE_{\chi_{1},r}}$, with total decay rate $\Gamma_{\chi}$, we may now compute the flux as

[TABLE]

For all of decays considered, we assume that they are isotropic in the rest frame. We can then integrate out $\cos\theta_{r}$

[TABLE]

As this assumption of isotropic decays continues for the rest of the text, we will drop the “iso” subscript. For two-body decays, we find

[TABLE]

In the case of three-body decays, we define the invariant mass square of the electron-positron pair as $m_{e\bar{e}}^{2}=(m_{\chi}+\delta)^{2}+m_{\chi}^{2}-2(m_{\chi}+\delta)E_{\chi_{1},r}$. The differential rate is then

[TABLE]

where

[TABLE]

Following the described prescription, we show a set of simulated CR boosted DM fluxes for various combinations of free parameters in the left plot of Fig. 2. The red (blue) lines indicate the flux from TXS 0506+056 (NGC 1068). We fix the DM and mediator mass and select gauge couplings to give a non-relativistic p-$\chi_{1}$ cross section (Eq. 18) of $10^{-30}\mathrm{cm^{2}}$. We find that the boosted DM fluxes peak at low DM energies, as expected given the falling power law spectra of CRs in AGN and the fact that scattering is suppressed when the momentum transfer is greater than the mediator mass. For the same choice of DM parameters, we notice that the boosted fluxes from NGC 1068 are larger than those from TXS 0506+056 at energies below $\sim 20$ GeV, while the order flips for higher energies. This behavior, along with the low-energy features on the boosted fluxes from TXS 0506+056, arise from its highly boosted jet.

IV.2 Deep-Inelastic Scattering (DIS)

We can also consider the contribution from DIS. We will parametrize this in terms of the Björken parameter $x$. We must now use a double differential cross section $\frac{d^{2}\sigma}{dxdT_{\chi_{2}}}$. However, as we are working with high energy scatterings, we can use the approximation

[TABLE]

where $f_{i}(x,q^{2})$ is the parton distribution obtained from [35] and we sum over up, down, anti-up, and anti-down quarks. For $\frac{d\sigma_{i}}{dT}$, we use Eq 17 with a form factor set to 1 and neglect the quark mass. The reason that we can take $s_{q}=xs$ where $s_{q}$ is the Mandelstam variable for the DM-quark interaction and $s$ is the variable for the overall DM-proton interaction is as follows. In the rest frame of the DM, the proton is highly boosted, and thus the four momentum of the proton $p_{p}^{\mu}$ and of the underlying quark $p_{q}^{\mu}$ are related by $p_{q}^{\mu}\simeq xp_{p}^{\mu}+\mathcal{O}(m_{p}/T_{p})$. Thus, when this kinetic energy is far higher than the masses of the particles involved, we know $s_{q}=(p^{\mu}_{q}+p^{\mu}_{\chi})^{2}\simeq xp_{p}^{\mu}p_{\chi,\mu}\simeq xs$. We can equivalently think of this as the quark having kinetic energy $T_{q}=\frac{1}{2m_{\chi}}\big(s_{q}-(m_{q}+m_{\chi})^{2}\big)$ which allows us to compute the scattering angle via Eq. 21. In this approximation, we have

[TABLE]

where the sum is over up, down, anti-up, and anti-down quarks. We use Eq. 17 with the form factor set to 1 for the cross section, find the scattering angles from the quark kinematics, and use parton distribution functions $f_{i}(x,q^{2})$ from [71, 35]. We consider that the coupling to quarks is 1/3 of the coupling to protons 333This is a conservative choice when extrapolating constraints to a dark photon model, as in that model up and down quarks have a coupling of 2/3 and -1/3 of the proton respectively.. Also, we only consider DIS for $q>1.295$ GeV, which is the lowest momentum transfer at which ManeParse includes parton distributions [35]. We also simulate the boosted DM fluxes at TXS 0506+056 and NGC 1068 including DIS as described. The results are shown in the right plot of Fig. 2. We note that the flux from NGC 1068 stays above the TXS 0506+056 flux when DIS is considered, as there is no longer a form-factor suppression at high energies. Above a certain energy threshold when the momentum transfer can be above 1.295 GeV, the DIS contribution starts to contribute and the boosted DM fluxes are enhanced by orders of magnitude with respect to the case without DIS. The DM energy threshold for DIS to occur depends on the DM mass and mass splitting, thus in Fig. 2 the fluxes present “bumps” at different energies depending on the labeled assumptions. The effects of DIS become more important at Super-K than at experiments with lower energy thresholds such as XENONnT. The spectral features induced by DIS may also help to differentiate DM from other cosmogenic and neutrino-induced backgrounds.

IV.3 Attenuation from Earth

The DM-SM cross sections considered here can be large enough that the Earth is no longer optically thin to the DM flux. In this case, we must now consider the effects of attenuation. To do so, knowing the location of both the AGN source and detector, we can find the optical depth $\tau$ as

[TABLE]

where $i$ sums over the different elements, $\ell$ is a line-of-sight integral, and $n_{i}$ is the number density of each element [43, 86, 84]. The cross section includes coherent scattering off nuclei, incoherent scattering on nucleons and electrons, and DIS (when considered at the source). We treat the detectors as 1 km underground in a spherically symmetric Earth. We note that the line of sight changes as the Earth rotates, and the cross section is an energy dependent quantity, so we cannot describe our problem with a single $\tau$.

To properly compute attenuation, we need to understand the interior of the Earth. For the density, we use the Preliminary Reference Earth Model [43]. For the elemental compositions, we find the core/mantle composition from [86] and the crust composition from [84]. The results are summarized in Table 2.

To set a “ceiling” on our bounds, we consider a test energy

[TABLE]

where $\Delta E=1$ GeV (0.1 GeV) for Super-K (XENONnT) is chosen as a characteristic energy for detecting scattering signals, and the second value assures that upscattering is possible for DM-proton scattering. We set the “ceiling” to be the minimum value of the couplings for which $\tau(E_{\rm test})\geq 1$ for all points during Earth’s rotation. For such strong attenuation, the flux at our detectors would be too small and we would be unable to set constraints.

For less severe attenuation, we must update the flux to account for attenuation, which at the detector becomes

[TABLE]

where the brackets signify averaging over the Earth’s rotation and the flux on the right-hand side is the flux without attenuation as obtained via Eq. 28 or Eq. 33 depending on if DIS is considered.

V Direct detection of cosmic-ray boosted DM at Super-K

In the previous section, we have calculated the flux of $\chi_{1}$ particles reaching the Earth. The scattering rate at the detector protons, electrons or nuclei (denoted as $i=p,e,A$), is given by

[TABLE]

where $N_{i}$ is the number of target protons or electrons in the detector, and $T_{i}$ denotes their recoil energy of the proton or electron, and $E_{1}$ is the incoming energy of the DM particle. The scattering cross section for $\chi_{1}$ with a particle $i$ is

[TABLE]

where as before the form factor is 1 for electrons and quarks. Blazar-boosted DM has been studied at Super-K previously [63]. In it, they study electron scattering in the 22.5 kton water Cherenkov detector, breaking their study into 3 bins based on recoil energy: Bin 1 [0.1 GeV, 1.33 GeV], Bin 2 [1.33 GeV, 20 GeV], Bin 3 [20 GeV, 1 TeV]. Furthermore, they make angular cuts of $24^{\circ}$, $7^{\circ}$, and $5^{\circ}$ respectively for the 3 bins. We show the electron scattering rates induced by CR boosted DM at Super-K applying the angular cuts in Fig 3. We notice that the scattering rate is enhanced at low energies, where the boosted fluxes are larger. The aforementioned angular cuts were chosen to contain most outgoing electrons from elastic scattering. As the kinematics for inelastic DM are different, we also consider the possible exclusions without the angular cuts, setting our 95% confidence level where the number of DM-induced scatters $N_{\rm BSM}=2\sqrt{N_{\rm Bkg}}$ (with backgrounds obtained from [76]). This allows us to set our limits at 20, 4.5, and 3 (150, 51, and 6) events with (without) angular cuts over 2500 days for Bins 1, 2, and 3 respectively. We assume that the detector has perfect efficiency, energy resolution, and angular resolution so that all events appear in the true angular and energy bins. Constraints can also be placed through interactions with protons at the Super-K detector. Super-K has also analyzed proton data for CR boosted DM, looking for events with momentum above 1.2 GeV and below 3 GeV [10]. In their sample, they predicted $111.7\pm 10.6(\mathrm{stat.})\pm 30.7(\mathrm{sys.})$ (a relative $27\%$ systematic uncertainty). As 126 events were observed, this means we can set a 90$\%$ confidence exclusion for any model producing 80 or more elastic proton-DM induced events. To do this, we consider just scattering off hydrogen atoms, ignoring oxygen.

We also consider DIS in the detector. Such events would produce multiple Cherenkov rings, so we can set a limit by using data in Super-K Run I-IV where 5755.9 multi-ring events are predicted compared to 5770 observed. Considering only statistical uncertainty, we would set $90\%$ confidence constraints at 114 DM-DIS events. We do not consider the energy distribution of DIS events but rather compute the entire rate as

[TABLE]

where $N_{p}$ and $N_{n}$ are the number of protons and neutrons in the detector and

[TABLE]

where the sum is again over up, down, anti-up and anti-down quarks. We can relate $\frac{d\sigma(s)}{dt}$ to Eq. 38 via $\frac{d\sigma(s)}{dt}=\frac{1}{2m_{i}}\frac{d\sigma}{dT_{i}}$ and using the substitutions $T_{i}\rightarrow-t/2m_{i}$ and $E_{1}\rightarrow\frac{1}{2m_{i}}(s-m_{\chi}^{2}-m_{i}^{2})$ where $m_{i}$ is the quark mass. As before, we take the coupling of quarks to $Z^{\prime}$ to be 1/3 that of proton.

We will derive constraints on various combinations of the free BSM parameters for a proper assessment of the strength of this phenomenological probe of DM. In Fig. 4, we show constraints on light DM with mass splittings comparable to the DM mass, a scenario widely studied in the literature since it can be probed at collider and beam dump experiments, and where the thermal relic abundance has not been completely probed. We show constraints from electron scattering at Super-K, both with and without considering DIS at the AGN. We also include constraints from DIS scattering at Super-K when the same is also considered at the AGN. Constraints from elastic proton recoils are weaker than the attenuation “ceiling” and are therefore not shown. As it can be appreciated in the figure, our bounds from CR boosted DM at NGC 1068 and TXS 0506+056 can improve over previous bounds from beam dump experiments and CR cooling in AGN, further testing the thermal freeze-out target for a wide range of inelastic DM masses. Concretely, Super-K is able to probe these models for DM masses $m_{\chi}\lesssim 0.1$ GeV ($\delta=0.4$mχ, $m_{Z^{\prime}}=10m_{\chi}$) and $m_{\chi}\lesssim 0.4$ GeV ($\delta=0.8_{\chi}$, $m_{Z^{\prime}}=3m_{\chi}$). It can also be appreciated in the plots that the inclusion of DIS enhances the sensitivity with respect to elastic scatterings at large DM masses. This is expected, as larger mass splittings require larger momentum transfers, where DIS has the higher cross section.

We can relax the previous assumptions on the mass splitting and relation between mediator mass and DM mass, in order to explore the strength of our constraints in other regions of the parameter space. We show our constraints from Super-K on the non-relativistic DM-proton scattering cross section $\sigma_{\chi-p}$ (see Eq. 18) versus DM mass $m_{\chi}$, for various values of the mass splitting $\delta$ ($=0.01,0.1,1,10\ \rm GeV$). We show results from $\chi_{1}$-electron scatterings at Super-K in Fig. 5. Complementary (weaker) results from proton recoils at Super-K are shown in the Appendix D. The left plots in Fig. 5 correspond to limits from TXS 0506+056, and the right plots correspond to NGC 1068. Moreover, we also show limits for various values of the mediator mass $m_{Z^{\prime}}=0.05,5,500$ GeV.

We note from the plots that the mass splitting significantly impacts the strength of the constraints at high-DM masses, where larger mass splittings generically lead to less stringent constraints, since the scattering process is kinematically suppressed with respect to smaller values of the mass splitting. We further notice that the constraints are more stringent for heavy mediators than for light mediators. This can be understood since the boosted DM fluxes are suppressed by the propagator in the DM-proton scattering cross section as $\sim 1/(m_{Z^{\prime}}^{2}+2m_{\chi}T_{\chi_{2}}-\delta^{2})^{2}$. The inverse suppression with the kinetic energy of the boosted DM $T_{\chi_{2}}$ manifests for lighter mediator masses over a wider range of kinetic energies than for heavy mediators, which leads to somewhat smaller boosted DM fluxes for heavy mediators for the energies relevant of Super-K, see Fig 2.

In the limit $\delta\rightarrow 0$, our constraints can be compared to complementary bounds from direct detection, CR cooling in AGN, and cosmological probes. These bounds are shown as a shaded red region in the figure. In particular, for $m_{\chi}\lesssim 7\times 10^{-4}$ GeV the leading bound arises from measurements of the number of relativistic species ($N_{\rm eff}$) at the time of Big Bang Nucleosynthesis (BBN) [41, 57, 79]. For $7\times 10^{-4}\,\mathrm{GeV}\lesssim m_{\chi}\lesssim 3\times 10^{-3}\,\mathrm{GeV}$, the leading bounds arise from CR cooling in NGC 1068 [69, 65, 87]. In the range $3\times 10^{-3}\,\mathrm{GeV}\lesssim m_{\chi}\lesssim 0.2\,\mathrm{GeV}$, the leading bounds arise from a search on direct detection of galactic CR boosted DM at Super-K [10]. Finally, for $m_{\chi}\gtrsim 0.2$ GeV, traditional direct detection searches of galactic halo DM sets the most stringent bounds [8, 1, 25].

The constraints from CR boosted DM from AGN can overcome previous bounds in the range $10^{-3}$ GeV $\lesssim m_{\chi}\lesssim 1$ GeV in the limit $\delta\rightarrow 0$ and $m_{Z^{\prime}}\gtrsim 0.05$ GeV, further testing large mass splittings. In the presence of sufficiently large mass splittings ($\delta\gtrsim$ 400 keV), direct detection constraints do not apply. In this regime, inelastic DM is better probed in beam dumps and colliders [73, 88, 53]. In Fig. 5, the shaded gray region represents the combined exclusion limits from several experiments. For $10^{-3}~\mathrm{GeV}\lesssim m_{\chi}\lesssim 10^{-2}~\mathrm{GeV}$, the strongest bound is set by NA64 [88], where the $Z^{\prime}$ is produced via dark bremsstrahlung and undergoes semi-visible decays, often with displaced visible products. In the range $10^{-2}~\mathrm{GeV}\lesssim m_{\chi}\lesssim 0.2~\mathrm{GeV}$, comparable limits arise from LSND [72], CHARM [59], and NuCal [104], where the $Z^{\prime}$ can be produced not only through dark bremsstrahlung but also via meson decays such as $\pi^{0}\to\gamma Z^{\prime}$, and the visible decay products of the heavier dark-sector state are typically displaced. For $m_{\chi}\gtrsim 0.2~\mathrm{GeV}$, the most stringent constraints come from collider searches, where the $Z^{\prime}$ is produced in $e^{+}e^{-}$ or $pp$ collisions and decays semi-visibly. These include BaBar, LEP, and LHC results [73]: BaBar and LEP bounds stem from mono-photon plus missing-energy searches and invisible $Z$-width measurements, while LHC limits are derived from mono-jet/photon signatures and Drell–Yan–like processes, with sensitivity to both prompt and displaced visible decays.

We can also compute the relic abundance for non-zero mass splittings in our set-up. In the common limits $m_{\chi_{1}}\gg 2m_{e}$, $m_{Z^{\prime}}>2m_{\chi_{1}}$ and $\delta<m_{\chi_{1}}$, the relic abundance is predominantly set by annihilations into electron-positron pairs, and the thermally-averaged cross section reads [gonzález2021cosmology]

[TABLE]

We confront our limits with the region of the parameter space that yields the DM relic abundance in the Universe via freeze-out ($\langle\sigma v\rangle\simeq 1.5\times 10^{-36}\mathrm{~cm}^{2}$), for two different choices of the mass splitting ($\delta=0.01$ GeV in dashed purple, $\delta=10$ GeV in solid purple). We note however that the displayed thermal benchmarks are approximate and their validity is restricted to the regime $\delta\lesssim m_{\chi_{1}}$.

It can be noticed in Fig. 5 that such freeze-out motivated region is fully tested for some mediator masses, while for others (e.g $m_{Z^{\prime}}=5$ GeV) only the region around the resonance remains viable. Additionally to showing the benchmark values of the DM relic abundance, we also shade in black color the region of parameter space that is excluded by perturbativity, $g_{p}g_{\chi}>16\pi^{2}$. We note that this consideration is not relevant for light mediators, but becomes important for heavy mediators. For instance, for $m_{Z^{\prime}}=500$ GeV, we find that our limits from CR boosted DM probe regions of parameter space which are already excluded by perturbativity.

The constraints from CR boosted DM in AGN will be further strengthened with the future Hyper-Kamiokande experiment [9]. We estimate the sensitivity enhancement to be mild though, within one order of magnitude. A brief discussion can be found in the Appendix A.

VI Conclusions

We have derived novel constraints on the DM-proton and DM-electron scattering cross sections from the consideration that CRs in NGC 1068 and TXS 0506+056 can scatter off the surrounding DM in the vicinity of their supermassive black holes (SMBHs), yielding a flux of boosted DM particles on Earth. This scenario remains viable even when there is a mass-splitting in the dark sector. High-energy neutrino and gamma-ray observations indicate that CRs are accelerated in the vicinity of the central black hole in these sources, with large luminosities. This, together with the expectation that the DM density is expected to be very large in these environments, induces a potentially measurable flux on Earth, compensating the dilution with extragalactic distances.

In particular, we have derived constraints from Super-K (and XENONnT in Appendix E). We considered a realistic choice for the CR acceleration radius in these sources and computed the CR flux using leptohadronic models informed by observations. Importantly, we considered a concrete model of DM-proton and DM-electron interactions consisting of a fermionic DM candidate and a vector mediator. In addition, we also allow for a mass splitting in the dark sector. The latter possibility introduces a novel complexity in the analysis and rich new phenomenology, further allowing to evade the stringent bounds from direct detection experiments on galactic gravitationally-bound DM. We included decays from the up-scattered heavier DM state into the lightest and a photon or an electron-positron pair, which affects non-trivially the corresponding boosted DM fluxes. Finally, we include deep-inelastic scattering (DIS) at the source and at the detector, which has significant implications on the ensuing limits, especially at Super-K.

We report progress on several fronts. First, we find that the boosted DM flux from NGC 1068 generically overcomes that from TXS 0506+056 when DIS is included, due to the expected larger DM density at the source, as inferred from its measured black hole masses. In addition, the CR luminosities that have been considered in previous work at TXS 0506+056 may only apply during flaring periods, much shorter than the running exposures at XENONnT and Super-K. NGC 1068 is expected to accelerate CRs steadily, in contrast. Second, we find limits on inelastic light DM scenarios that overcome complementary bounds from collider and beam dump experiments, and allow us to test uncharted parameter space of thermal DM. Third, we find constraints in the non-relativistic DM-proton scattering cross section that can overcome previous bounds from direct detection, galactic CR boosted DM, CR cooling in AGN and BBN/CMB in the mass range from $m_{\rm DM}\sim 1$ MeV to $m_{\rm DM}\sim 1$ GeV, probing thermal DM models. When considering a mass splitting among the two DM states, the constraints from direct detection shut off, and our phenomenological probes previously uncharted parameter space over a range of DM masses and mass splitting spanning several orders of magnitude.

A caveat of our work (and previous works on CR boosted DM) concerns the expected DM self-scatterings at the source, which we neglected, but that could affect non-trivially the energy distribution of the boosted DM fluxes on Earth. We will address the impact of self-scatterings in depth in upcoming work. Another important caveat is related to the DM density profile and corresponding DM column density around the SMBH in AGNs. Even if DM spikes are not formed and sustained in time, an NFW profile would still yield a large column density, although $\sim$ 3 to 4 orders of magnitude weaker than it is customarily assuming a DM spike [69, 65].

Future observations of high-energy neutrinos at IceCube from other AGN may allow us to set more stringent constraints by accumulating further statistics and reducing the uncertainties on the inferred CR spectra at these sources. Independent observational methods to infer the DM distribution in these environments are also crucial [102], to reduce the degeneracy with the DM-CR scattering cross section. These points, together with the upcoming improvements in direct detection experiments (like XLZD or Darkside-20k [2, 3]) and Hyper-Kamiokande [9], will allow us to dive into thermal light/and or inelastic DM models further, hopefully leading us to even stronger constraints or better – a detection.

Acknowledgments

A.G., G.H., and I.M.S. are supported by the U.S. Department of Energy under the award number DE-SC0020262. M. M. and K. M. are supported by NSF Grant No. AST- 2108466. A.G. is partially supported by the World Premier International Research Center Initiative (WPI), MEXT, Japan. A.G. is grateful to QUP for hospitality during his visit. M. M. also acknowledges support from the Institute for Gravitation and the Cosmos (IGC) Postdoctoral Fellowship and the FermiForward Discovery Group, LLC under Contract No. 89243024CSC000002 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics. The work of K.M. is also supported by the NSF Grants, No. AST-2108467 and No. AST-2308021, and KAKENHI No. 20H05852.

This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Workforce Development for Teachers and Scientists, Office of Science Graduate Student Research (SCGSR) program. The SCGSR program is administered by the Oak Ridge Institute for Science and Education (ORISE) for the DOE. ORISE is managed by ORAU under contract number DESC0014664. All opinions expressed in this paper are the author’s and do not necessarily reflect the policies and views of DOE, ORAU, or ORISE. A.G. is a recipient of the SCGSR Award.

Appendix A Scaling

In this work, we have considered AGN with properties as specified in Table 1. Better future observations or models may alter our understanding of these values. Changes to $\Gamma_{B}$, $\theta_{\mathrm{LOS}}$, $\alpha_{p}$, and $\gamma^{\prime}_{\max/\min,p}$ lead to non-trivial changes in the resulting event rate. Alterations to the other parameters, however, can be captured in a simple scaling relation. In this appendix, we will derive such a scaling.

We know that our rate will follow Eq. 37 (repeated below)

[TABLE]

From Eq. 28, we know that the flux scales as

[TABLE]

where $\Sigma_{\chi}$ itself depends upon the black hole mass, acceleration region, and extent of the DM spike. We know that $\sigma\sim g_{\chi}^{2}g_{i}^{2}$. Thus, our number of events at a detector $N_{ev}$ will scale as

[TABLE]

where $\mathcal{E}$ is the detector exposure. Therefore, if a number of events $N_{des}$ is required by an experiment to set an exclusion, limits on the coupling go as

[TABLE]

Note, this assumes that the attenuation through the Earth is roughly the same. This is valid as long as the coupling does not exceed the ceiling value as calculated in Sec. IV.3.

As a quick exercise, let us estimate how our exclusions will change with data from Hyper-Kamiokande which is expected to have 10 times the exposure of Super-K. For energy/angular bins dominated by background, the desired number of events will scale as $\mathcal{E}^{1/2}$, while for background-free bins, the desired number of events is fixed. To that point we can estimate the future exclusions $(g_{\chi}g_{i})_{HK}$ to relate to our current exclusions $(g_{\chi}g_{i})_{SK}$ by

[TABLE]

Appendix B Different Cases of NGC 1068

As indicated in Fig. 1 and Table 1, there are multiple plausible cases for the proton spectrum at NGC 1068. Although our main text figures only show the results from Case 1 (C1), here we explore the impact from the other cases.

We see that for the case of $\delta=0.4m_{DM}$ and $m_{Z^{\prime}}=10m_{DM}$, Case 2A sets the strongest constraints, while Case 2B sets the weakest. Therefore, for these moderate mass splittings, we can treat Cases 2A and 2B as a range of possible constraints. For the larger mass splittings $\delta=0.8m_{DM}$ and electron scattering, Case 1 becomes the strongest for $m_{DM}\gtrsim 0.3$ GeV. This can be understood as large mass splittings requiring higher energies, and Case 1 has the highest luminosity as large energies. This effect is less extreme for proton scattering at Super-K, and Case 2A continues to set the strongest bounds.

Appendix C Upper limits from $\chi_{2}$ scatterings on Earth

In the main text, we assumed that all upscattered $\chi_{2}$ decay long before reaching Earth. While this is a reasonable assumption considering the mass-splittings and distances in this work, it can still be elucidating to consider the impact of $\chi_{2}$ at Earth assuming no decays 444In principle, one could consider a case where only a fraction of the DM decays to the ground state, but in this work we only consider the two extreme cases.

First, we must find the flux of $\chi_{2}$ at Earth. This is relatively easy considering elastic scattering ($p+\chi_{1}\rightarrow p+\chi_{2}$). Since Earth is on-axis for the jet, we know that if the DM scatters at angle $\theta$, the original proton must have been oriented at angle $\theta$ to the jet (NGC 1068 does not have a jet, but its emission is isotropic, so the production angle is unimportant). Therefore, we can find our flux as

[TABLE]

where $\cos\theta(T_{p},T_{\chi_{2}})$ is given by Eq. 21.

When DIS is included, the contribution to the $\chi_{2}$ flux is

[TABLE]

Then, the differential rate of elastic scattering can be determined as

[TABLE]

with the cross section for $\chi_{2}$ scattering with point-particle $i$ is

[TABLE]

Form factors are included for protons and nuclei as needed.

For DIS, we can compute

[TABLE]

where once again $\frac{d\sigma(s)}{dt}=\frac{1}{2m_{i}}\frac{d\sigma}{dT_{i}}$ and $T_{i}\rightarrow\frac{-t}{2m_{i}}$. We let $E_{2}\rightarrow\frac{1}{2m_{i}}\big(s-(m_{\chi}+\delta)^{2}-m_{i}^{2}\big)$. We show the resulting constraints in Fig. 7.

Appendix D Upper limits from proton recoils at Super-K

As mentioned in the main text, we also derived upper limits on the DM-proton scattering cross section at Super-K from proton recoil data. This scenario, where the DM is boosted by CR protons at the source, and leaves signatures by proton recoils at the detector, allows us to connect our results with models where the DM is leptophobic, or couples predominantly to quarks. The upper limits are shown in Fig. 8, for TXS 0506+056 and NGC 1068 CR model C1 as an example. These limits are weaker than those obtained from electron recoils, but can still be leading in some regions of parameter space. Importantly, if the inelastic DM relic abundance is set by annihilations into protons, we find that thermal freeze-out scenarios are severely constrained.

Appendix E Upper limits from nuclear recoils at XENONnT

We also derive limits from nuclear recoils at the XENONnT experiment [16]. Just as we did for Super-K, we can use Eq. 37 to calculate the differential recoil rate at XENONnT. Then, in order to find the total number of events in the energy window of the experiment, we calculate

[TABLE]

where the minimum and maximum recoil energies of the xenon nuclei A are given by $T_{A}^{\rm min}=3.3$ keV and $T_{A}^{\rm max}=60.5$ keV. $\epsilon(T_{A})$ denotes the efficiency of the experiment, which we take from [16]. $\mathcal{E}$ refers to the exposure of the experiment, which is 1.09 tonne $\times$ yr. Finally, in order to derive an upper limit on the gauge couplings, mass splittings or non-relativistic cross section, we find the 90% C.L upper limit on the number of signal events from a Poissonian likelihood as

[TABLE]

where the number of observed events at XENONnT is $N_{\rm obs}=152$, and the expected number of background events is $N_{\rm bck}=149$. The 90$\%$ C.L upper limit on $N_{\rm sig}$ can then be obtained from the test statistic

[TABLE]

getting $N_{\rm sig}^{90\%\rm C.L}<24.17$. Thus, we impose that the total rate induced by certain combinations of the parameters that we aim to constrain, confer Eq. 52, shall not exceed this upper limit. Under this prescription, we find the limits from XENONnT to be significantly weaker than those from Super-K.

Our results are shown in Fig. 9. Depending on the values of the mass splitting, mediator mass, and $\chi_{1}$ mass, the flux may survive the attenuation in the atmosphere and reach the detector in some instances, but in most regions of parameter space the DM does not reach the detector (see Sec. IV.3 for more details).