Implications of photon-ALP oscillations in the extragalactic neutrino source TXS 0506+056 at sub-PeV energies

We consider the photon-ALP oscillation at the source in the presence of blazar jet magnetic field (BJMF). The magnetic field in the jet region can be modeled with poloidal (along the jet axis, $B\propto r^{-2}$) and toroidal (transverse to the jet axis, $B\propto r^{-1}$) components. At distances large enough from the central black hole, the toroidal component dominates over the poloidal component and thus the latter can be neglected. We adopt the toroidal magnetic field strength $B^{jet}(r)$ given by [54, 55]:

where $B^{jet}_{0}$ is the magnetic field strength at the core and $r_{VHE}$ is the distance between the VHE $\gamma$-ray-emitting region and the central black hole. We assume $r_{VHE}\sim R_{VHE}/\theta_{j}$, where $R_{VHE}$ is the blob radius of the VHE emitting region and $\theta_{j}$ is the angle between the jet axis and the line of sight.

Assuming equipartition between the magnetic field and particle energies, the electron density profile $n_{el}(r)$ can be modeled as a power law given by [56]:

where $n_{0}$ is the electron density at $r_{VHE}$. In Ref. [57], a more realistic model was provided that takes into account the fact that the electron distribution is nonthermal in a relativistic active galactic nuclei Jet.

The photon energy $E^{{}^{\prime}}_{\gamma}$ in the jet frame is related to the lab-frame energy $E_{\gamma}$ by $E^{{}^{\prime}}_{\gamma}=E_{\gamma}/\delta_{D}$, where $\delta_{D}=\left[\Gamma_{L}(1-\beta_{j}^{2}cos\theta_{j})\right]^{-1}$ is the Doppler factor with $\Gamma_{L}$ and $\beta_{j}$ being the bulk Lorentz and beta factor, respectively. We assume that at $r>1$ kpc the BJMF strength is negligible. Further details of the BJMF model can be found in Refs. [7, 8].

The intercluster medium magnetic field (ICMF), $B^{ICMF}$ can be modeled as,

with $0.5\leq\eta\leq 1.0$, where $B_{0}^{ICMF}$ and $n_{el}(r_{0})$ are the magnetic field strength and electron density at the cluster center, respectively. The electron density distribution $n_{el}(r)$ at a distance $r$ from the cluster center is

with $r_{core}$ is the core radius and $\zeta=-1$. The typical order of the electron density and the core radius is $\sim\mathcal{O}(10^{-3})$ cm^-3 and $\sim\mathcal{O}(100)$ kpc, respectively.

In a cluster-rich environment, the turbulent magnetic field is of order $\sim$ $\mathcal{O}$(1)$\mu$G [58, 59, 60]. Such a cluster can have a significant effect on the conversion between photons and ALPs [61, 62]. Since there is no evidence that TXS 0506+056 is located in a cluster-rich environment, we do not consider the photon-ALP oscillations in the ICMF model.

The actual strength of the extragalactic magnetic field on the cosmological scale $\sim\mathcal{O}$(1) Mpc, is still unknown, but the currently accepted limit is $\sim\mathcal{O}$(1) nG [63, 64]. In this work, we neglect the effect due to the magnetic field (See Ref. [65] for possible effects of the magnetic field) and consider only the absorption effect due to EBL.

The optical depth of EBL attenuation can be written as [66]

with

and $E_{th}=2\cdot(m_{e}c^{2})^{2}/(E_{\gamma}(1-cos\theta)),$ where $z_{0}$ is the redshift of the source, $H(z)$ is the rate of Hubble expansion, $E_{th}$ is the threshold energy for pair production, $dn(z)/dE^{BG}_{\gamma}$ is the proper number density of the EBL, $\sigma_{\gamma\gamma}(E_{\gamma},z,E^{BG}_{\gamma})$ is the pair-production cross section, $\theta$ is the angle between the projectile and target photons, $E_{\gamma}$ is the projectile photon energy, and $E^{BG}_{\gamma}$ is the target background photon energy. Several EBL models have been proposed in the literature [67, 68, 69, 70, 71, 72, 73], and we consider the model from Ref. [74] in this work.

Finally, we consider the effect of photon-ALP oscillation in the presence of a Galactic magnetic field (GMF). This effect can have both a large-scale regular component and a small-scale random component. Due to the fact that the coherence length is smaller than the oscillation length, we neglect the random component [61] in our analysis and consider only the regular GMF component model given in [75]; the latest model can be found in Refs. [76, 77].

We analyze the Fermi-LAT data for TXS 0506+056 in three phases:

• •

  Neutrino Flare : September 2014 (MJD 56937.81) to March 2015 (MJD 57096.21).

• •

  VHE Flare 1 : 04 September, 2017 (MJD 58000) to 03 November, 2017 (MJD 58060).

• •

  VHE Quiescent : 04 November, 2017 (MJD 58061) to 10 October, 2018 (MJD 58422).

We select the above-mentioned time period of Fermi-LAT (Pass 8) processed data from the Fermi Science Data Center ³³3https://fermi.gsfc.nasa.gov/ssc/data/access/. We choose the spacecraft files and the instrument response functions (IRFs) of P8R3_SOURCE_V2 to match the extracted data. We select the SOURCE event class (evclass=128 and evtype=3) in the energy range of 100 MeV to 300 GeV. The region of interest is selected to be 2^∘ centered on the target source. We consider all of the 4FGL sources around 10^∘ of TXS 0506+056 as background sources together with preprocessed templates of Galactic diffuse emission, gll_iem_v08.fits , and the extragalactic isotropic diffuse emission, iso_P8R3_SOURCE_V2_v2.fits.

The likelihood analysis, spectral energy distribution, and light curve are obtained by using the PYTHON-based Fermipy package ⁴⁴4https://fermipy.readthedocs.io/en/latest/index.html [82].

We calculate the photon-ALP oscillation probability for the three phases of TXS0506+056 mentioned above following Sec. IV.

Extending this survival probability to sub-PeV energies, we can estimate the residual photon flux of the observed neutrino’s counterpart. For this calculation, we use the ALP-$\gamma$ oscillation using the CAST upper limits on the parameters, namely $m_{a}=1$ neV and $g_{a\gamma}=6.6\times 10^{-11}$ GeV^-1.

Assuming that IC170922-A resulted from proton-proton interactions, as calculated in Ref. [83], the average counterpart $\gamma$ rays will follow [84], $E^{2}_{\gamma}\cdot\frac{dN_{\gamma}}{dE_{\gamma}}=\frac{2}{3}E^{2}_{\nu}\cdot\frac{dN_{\nu}}{dE_{\nu}}\,,$where $E_{\gamma}\approx 2E_{\nu}$ is the energy of the photons produced from $\pi^{0}$ decay. Subsequently, these VHE photons attenuate by interacting with the synchrotron and synchrotron self-Compton photons, resulting from relativistic electrons inside the blob. The fraction of VHE photons that can escape from the blob is estimated as in Ref. [85],

where the optical depth due to the interaction of high-energy photons of energy $\epsilon_{\gamma}^{\prime}$ in the comoving frame is

$n^{\prime}_{k}(\epsilon^{\prime}_{k})$ is the number density of the ambient photons of energy $\epsilon^{\prime}_{k}$ (in $m_{e}c^{2}$) in the comoving frame,

and $\sigma_{\gamma\gamma}$ is the pair-production cross section [86]. $R^{{}^{\prime}}_{blob}$ is the radius of the blob, and $D^{2}_{L}$ is the luminosity distance of the source. $\Gamma_{k}$ is the bulk Lorentz factor and $F_{k}(\epsilon_{k})$ is the photon flux in the observer frame.

Hence, the intrinsic $\gamma$-ray spectrum can be obtained by counterpart gamma rays as in Sec. VI.2 and this is shown by the grey vertical hatched and dot-dashed curves in the top panel of Fig. 4. This flux generally gets exhausted due to interaction with the EBL and CMB. On multiplying the intrinsic flux with the photon survival probability under photon-ALP oscillations for the corresponding phase, we can have a surviving fraction of the flux. The $\gamma$ rays due to the ALP effect at sub-PeV energies is shown by the red cross hatched and dashed curves. Interestingly, the surviving photons, considering the ALP-$\gamma$ oscillations in VHE Flare 2014, are above the differential sensitivity to a Crab-like point gamma-ray source for 1 year of exposure by LHAASO [87]..

In this section, we calculate the diffuse flux at sub-PeV energies as an effect of ALP-$\gamma$ oscillation in sources like TXS 0506+056. The authors of Refs. [88, 38] emphasized that TXS 0506+056 cannot be considered a blazar of BL Lac type, but rather an intrinsically FSRQ, and all FSRQs are of the low-energy (synchrotron) peaked (LBL). Hence, we calculate the diffuse ALP $\gamma$ flux for the source luminosity function (LF) evolution like FSRQs using

where $dN/d\Gamma$ is the intrinsic photon index distribution which is assumed to be a Gaussian, $d^{2}V/dzd\Omega$ is the comoving volume element per unit redshift per unit solid angle, $dF_{\gamma}^{int}/dE$ is the intrinsic photon flux, here taken as calculated in Sec. VI.2 for TXS 0506+056, $\tau^{alp}_{\gamma\gamma}$ is the opacity under the ALP scenario considering the CAST upper limit on the parameters $m_{a}$ – $g_{a\gamma}$, and $\rho(L_{\gamma},z)$ is the gamma-ray luminosity function (GLF).

We consider here the luminosity-dependent density evolution of the GLF with parametrization given by:

with

where $z_{c}(L_{\gamma})=z_{c}^{*}\left(\frac{L_{\gamma}}{10^{48}}\right)^{\alpha}\,$. We collect the best-fit parameters $A,L_{c},z_{c}^{*},\alpha,\eta 1,\eta 2,\delta_{1}$, and $\delta_{2}$ for three blazar classes, namely, FSRQ [89], HSP, and LISP sources [90]. The limits of integration are $\Gamma_{min}=1.2$, $\Gamma_{max}=3.0$, $z_{min}>0$, and $z_{max}=2$. For FSRQs, we also calculate the diffuse flux using the integration limits of Ref. [89] as a comparison.

Figure 5 shows the diffuse $\gamma$ ray flux expected from photon-ALP conversion for the three classes. The data points of diffuse gamma-ray emission from the Galactic plane recorded by the LHAASO-KM2A [91] and Tibet AS-$\gamma$ [92] experiments are also shown. We show each class with (i) $\L_{\gamma}\in[10^{40},10^{52}]$ erg/sec, and (ii) $\L_{\gamma}\in[10^{42},10^{52}]$ erg/sec.

Dedicated surveys of the extragalactic diffuse gamma-ray flux by observatories like LHAASO, Tibet AS-$\gamma$, and the upcoming Cerenkov Telescope Array [93] will be able to constrain the ALP parameters considering photon-ALP oscillations at sub-PeV energies as proposed here.