Multi-wavelength spectral modelling of the candidate neutrino blazar PKS 0735+178

The BL Lac object PKS 0735+178 was in its historic $\gamma$-ray brightness state during December 2021. This period also coincides with the detection of a neutrino event IC211208A, which was localized close to the vicinity of PKS 0735+178. We carried out detailed $\gamma$-ray timing and spectral analysis of the source in three epochs (a) quiescent state ($E_{1}$), (b) moderate activity state ($E_{2}$) and (c) high activity state ($E_{3}$) coincident with the epoch of neutrino detection. During the epoch of neutrino detection ($E_{3}$), we found the largest variability amplitude of 95%. The $\gamma$-ray spectra corresponding to these three epochs are well fit by the power law model and the source is found to show spectral variations with a softer when brighter trend. In the epoch $E_{3}$, we found the shortest flux doubling/halving time of 5.75 hrs. Even though the spectral energy distribution in the moderate activity state and in the high activity state could be modeled by the one-zone leptonic emission model, the spectral energy distribution in the quiescent state required an additional component of radiation over and above the leptonic component. Here we show that a photo-meson process was needed to explain the excess $\gamma$-ray emission in the hundreds of GeV which could not be accounted for by the synchrotron self-Compton process.


INTRODUCTION
Blazars, among the persistent luminous objects in the Universe with luminosity ranging from 10 42 − 10 48 erg/sec are a class of active galactic nuclei (AGN).They are believed to be powered by the accretion of matter onto supermassive black hole with masses of the order of 10 6 − 10 10 M⊙ located at the centre of galaxies (Lynden-Bell 1969;Shakura & Sunyaev 1973).Their radiation is dominated by non-thermal processes from relativistic jets that extend from scales of pc to kpc and are oriented close to the line of sight to the observer.They show rapid and large amplitude flux variations across the electromagnetic spectrum on time scales ranging from minutes to years (Wagner & Witzel 1995;Ulrich et al. 1997).In addition to flux variations, the optical emission from blazars is highly polarized (Angel & Stockman 1980) with significant polarization variations (Andruchow et al. 2005;Rakshit et al. 2017;Pandey et al. 2022).They are also believed to be cosmic accelerators and could be the sources of astrophysical neutrinos (Mannheim 1993a;Plavin et al. 2020;IceCube Collaboration et al. 2018;Giommi & Padovani 2021;Bhatta 2022;Buson et al. 2023).
Blazars are conventionally divided into two sub-classes ⋆ E-mail: athira.bharathan@res.christuniversity.innamely flat spectrum radio quasars (FSRQs) and BL Lacertae type objects (BL Lacs) based on the optical spectral properties.FSRQs show broad emission lines in their optical spectra while BL Lacs have either featureless optical spectra or spectra with weak emission lines with equivalent widths lesser than 5 Å.According to Ghisellini et al. (2011), a more physical distinction between FSRQs and BL Lacs could be based on the luminosity of the broad line region (BLR) with FSRQs having LBLR/L Edd > 5 × 10 −4 .Here, L Edd is the Eddington luminosity given by L Edd = 1.3 × 10 38 (M/M⊙) erg s −1 .The broadband spectral energy distribution (SED) of blazars displays two prominent peaks.The low energy component that peaks between the infrared and X-ray bands is well understood as the synchrotron emission from a relativistic population of electrons in the blazar jet.The origin of the high energy component, peaking at the MeV -GeV energy range (Fossati et al. 1998), is highly debated with models that promote a radiation scenario involving either leptons or hadrons or combination of both (Böttcher et al. 2013;Cerruti 2020).Depending on the location of the low energy synchrotron peak of the SED, blazars are sub-divided into low synchrotron peaked (ν peak < 10 14 Hz), intermediate synchrotron peaked (10 14 Hz < ν peak < 10 15 Hz) and high synchrotron peaked (ν peak > 10 15 Hz) blazars (Abdo et al. 2010).While FSRQs are predominantly low synchrotron peaked blazars, BL Lacs belong to all three classes with a ma-jority of them belonging to the high synchrotron peak category.
In the leptonic emission model of blazar jets, the high energy emission component is interpreted as the inverse Compton scattering of low-energy photons (Abdo et al. 2010).The target low energy photons for inverse Compton emission can be the synchrotron photons, commonly referred as synchrotron self Compton or SSC (Konigl 1981), or the photons external to the jet called the external Compton or EC (Begelman & Sikora 1987).Under this emission scenario, one would expect to see a close correlation between optical and GeV flux variations since the same electron population is responsible for the emission at these energies.In the hadronic emission model, the high energy component is interpreted as a result of hadronic processes.Under this scenario one need not expect to see a correlation between optical and GeV flux variations since different species are responsible for the emission at these energies.However, systematic investigation of the correlation between optical and GeV flux variations have led to varied results with instances of (a) correlation between optical and GeV variations (b) instances of optical flare without a GeV counterpart and (c) GeV flare without the corresponding optical counterpart (Rajput et al. 2019(Rajput et al. , 2020(Rajput et al. , 2021)).Thus, correlated studies of flux variations in the optical and GeV-band could not constrain the high energy emission process unambiguously.
Similarly, the correlation observed between the radio and GeV γ-ray fluxes during various activity states support the region responsible for these emission to be co-spatial.The radio emission is due to synchrotron radiation by the relativistic electrons in the jet, while the GeV emission is due to the inverse Compton process.Interestingly, often radio variations are found to lag the GeV variations (Max-Moerbeck et al. 2014;Rani et al. 2014;Esposito et al. 2015;Kramarenko et al. 2022;Yuan et al. 2023).This lag may be associated with the slow cooling of the radio emitting electrons compared to the γ-ray.
In the hadronic model of the high energy emission from blazar jets, the plausible emission mechanisms can be proton synchrotron (Aharonian 2000) or hadronic cascades (Mannheim 1993b).Accelerated protons on interaction with cold protons or low energy photons in the surrounding medium could lead to the production of neutrinos as well as high energy γ-rays (Prince et al. 2023).Additionally, there are now increased evidence of blazars being extragalactic neutrino sources (Buson et al. 2023;Plavin et al. 2023b).The correlation found between the hard X-ray emission and neutrino emission in blazars (Plavin et al. 2023a) tends to suggest that neutrinos are likely to come from highly beamed blazars strong in X-rays.This led to arguments favoring neutrinos being produced in proton-photon interaction.The association of neutrinos with blazars has also revived the investigation on the hadronic origin of high energy emission from blazar jets.Therefore, a comparative analysis of the SED characteristics of candidate neutrino blazars during the epoch of neutrino detection and quiescent/other flaring epochs not coincident with neutrino detection could provide the needed constraints on the high energy emission process in blazars.
The first blazar found to be associated with the detection of neutrinos by the IceCube collaboration is TXS 0506+056 observed on 22 September 2017 (IceCube Collaboration et al. 2018).The detection of neutrinos was coincident in direction and time with the γ-ray flare from the source.Since then, few blazars have been found to be spatially coincident with the IceCube neutrino events (Paliya et al. 2020;Boettcher et al. 2022;Sahakyan et al. 2023).Generation of broadband SED of a large sample of neutrino blazars in a homogeneous manner and their systematic modeling can provide the key to enhancing our understanding of the high energy component in their broadband SED.We are carrying out such an investigation and in this work, we present our results on the source PKS 0735+178.
The intermediate synchrotron peaked blazar, PKS 0735+178 is one of the brightest BL Lac objects in the sky.It was identified as a BL Lac first by Carswell et al. (1974) and is situated at a redshift of z = 0.45 (Rector & Stocke 2001).However, recently, Falomo et al. (2021) suggest a redshift of z = 0.65.It was detected in the γ-rays by the Energetic Gamma Ray Experiment Telescope onboard the Compton Gamma Ray Observatory (Hartman et al. 1999) and is also detected by the Fermi Gamma Ray Space Telescope, (hereafter Fermi; Abdollahi et al. (2020)).Recently, it was found to be in spatial coincidence with multiple neutrino events by IceCube (IceCube Collaboration 2021), Baikal (Dzhilkibaev et al. 2021), Baksan (Petkov et al. 2021), andKM3NeT (Filippini et al. 2022) detectors.This detection of neutrinos by multiple detectors was coincident with the largest flare ever observed in this source in the optical, UV, soft X-ray, and γ-ray bands.In this work, we carried out a detailed analysis of this source, including the variability, gamma-ray spectral study, and broad-band SED analysis with the help of multiwavelength data.We also modeled the broadband SEDs using simple one-zone leptonic model and the discrepancy of this model in explaining the γ-ray spectrum is inferred with the hadronic photo-meson process.This paper is organized as follows: Multi-wavelength data and reduction is described in Section 2, the analysis, results, and the broadband spectral fitting are discussed in Section 3 and the results of the work are summarized in Section 4. Throughout this paper, we used the cosmological constants ΩM = 0.3, ΩΛ = 0.7, and H0 = 71 kms −1 M pc −1 .

MULTI-WAVELENGTH DATA AND REDUCTION
The observational data used in this work was from the Large Area Telescope (LAT; Atwood et al. (2009)) on board the Fermi Gamma-ray Space Telescope over the period of about 14 years from August 2008 to February 2022.In addition to the γ-ray data, we also used UV, optical, and X-ray data covering the same period from the Swift telescope.

γ-ray
We used all data for PKS 0735+178 collected for the period 2008 August to 2022 February (MJD: 54749−59611) from the Fermi archives.We utilised the fermipy package 1 , a python package that facilitates analysis of data from the LAT with the Fermi Science Tools for the generation of spectral points.We extracted photon-like events categorized as 'evclass=128, evtype=3' with energies 0.1⩽E⩽300 GeV γrays within a circular region of interest (ROI) of 15 • centered on the source.We applied the recent isotropic model, "iso_P8R2_SOURCE_V6_v06" and the Galactic diffuse emission model "gll_iem_v06".For generation of the γ-ray light curves, we considered the source to be detected if the test statistics (TS) > 9, which corresponds to a 3σ detection (Mattox et al. 1996).For bins with TS < 9, we considered the source as undetected.

X-ray
We used data from the Swift X-ray Telescope (XRT) collected from the HEASARC2 archives for X-rays with energies between 0.3 − 10 keV from August 2008 to February 2022.We reduced the data using the default parameter values in accordance with the instructions provided by the instrument team.The source spectra were chosen from a region of radii 60", whereas the background spectra were taken from a region of radii 120".We combined the exposure map using the tool XIMAGE and created the ancillary response files with xrtmkarf.We used an absorbed simple power law model with the Galactic neutral hydrogen column density of 3.74 ×10 20 cm −2 to perform the fitting within XSPEC (Arnaud 1996).

UV and Optical
For the analysis of UV and optical data from August 2008 to February 2022, we utilized the data obtained from the Swift-UV-Optical telescope (UVOT), an instrument onboard the Swift spacecraft.We used the data from the filters V, U, W1 and W2.The central wavelength(FWHM) of these filters are 5468 Å(796 Å), 3465 Å(785 Å), 2600 Å(693 Å) and 1928 Å(657 Å) (Poole et al. 2008).We processed the data from Swift-UVOT using the online tool provided by the telescope's data archive.In order to account for the effects of galactic absorption, we applied corrections to the UV and optical data points during SED analysis.

Multi-wavelength light curve
We utilized the one week binned γ-ray light curve in the 100 MeV to 300 GeV band generated using the procedure outlined in Section 2.1.The multi-wavelength light curves spanning about 14 years that include γ-ray, X-ray, UV,  During this epoch, the source showed the largest flare ever observed in the optical, UV, X-rays and γ-rays.

Long term variability
To quantify the variability of the source on time scales of weeks, we calculated the fractional variability amplitude (Fvar) in different energy bands following Vaughan et al. (2003).We define where S 2 is the variance, x is the mean and σ 2 err is the mean square of the measurement error on the flux points.The uncertainty on Fvar is given by (Vaughan et al. 2003) Here, N represents the number of flux points in the light curve.We calculated the γ-ray variability for E1, E2, E3 and for the total light curve.The results of the variability analysis are given in Table 1.

Short term variability
Blazars are known to show γ-ray flux variations on time scales lesser than an hour.(Foschini et al. 2011;Saito et al. 2013;Pandey & Stalin 2022).Detection of such short time scale variation could enable one to constrain the size and the location of the γ-ray emission region.We, therefore, searched for the presence of flux variations on very short time scales (of the order of hours) in the γ-ray light curve of PKS 0735+178.For this, we identified the epoch of very high γ-ray activity, namely, E3 and generated one day binned light curves.We then calculated the flux doubling/halving time scale using the relation (3) Here F(t) and F(t0) are the fluxes at time t and t0 respectively and τ is the flux doubling/halving time scale.This calculation was done with the condition that the flux difference between epochs t and t0 is greater than 3σ (Foschini et al. 2011).Using the flux doubling time scale, we constrained the size of the γ-ray emitting region as where δ is the Doppler Factor and τ is the flux doubling/having time-scale and c is the speed of light.

Spectral Variability
To investigate spectral variations in the source, during the epochs analyzed in this work, we followed a modelindependent approach of estimating the hardness ratio (HR) and then investigated the dependence of HR on the total flux of the source.For this, we generated one day binned γ-ray light curves in two energy ranges, namely, 0.1−10 GeV and 10−300 GeV.We define HR and the associated errors in HR as The generated light curves in the two energy ranges and the corresponding HR for the three epochs are shown in Fig. 2. From the light curves, it is evident that the hard flux counts are lower compared to the soft flux counts.We carried out a weighted linear least square fit to the HR v/s intensity diagram in Fig. 2. The results of the fit are given in Table 2.We found significant evidence of spectral variation in only one epoch, namely E1, wherein we found a softer when brighter trend.This is at odds with the recent finding of a harder when brighter trend seen in the γ-ray band (Fang et al. 2022).Close to the period of neutrino detection, in the X-ray band the source was found to show a softer when brighter trend, generally seen in the FSRQ category of blazars (Prince et al. 2023).

γ-ray spectra
The slope of the γ-ray spectrum can provide hints on the intrinsic accelerated electron distribution responsible for this emission process.We carried out the analysis of the γ-ray spectra of the three epochs using two spectral models namely the power law (PL) and log parabola (LP) models.The PL model has the functional form as Here, dN(E)/dE is the differential photon number (cm −2 s −1 MeV −1 ), N0 is the normalisation and Γ is the photon index.
The LP model is defined as (Nolan et al. 2012) where, α is photon spectral slope at energy E• and β defines the peak spectral curvature of the SED.We used the Maximum Likelihood estimator gtlike and the likelihood ratio test (Mattox et al. 1996) to check the PL model (null hypothesis) against the LP model (alternative hypothesis).We calculated the curvature of test statistics T Scurve = 2(log LLP -log LP L) following Nolan et al. (2012).We tested the presence of a significant curvature by setting the condition T Scurve > 16. fit are given in Table 3. Significant curvature in the γ-ray spectrum is not found in the three epochs.

Broad band SED analysis
We modeled the generated SEDs of the three epochs using the one zone leptonic model (Sahayanathan et al. 2018).In this model, the emission region was assumed to be a spherical blob of size R filled with non-thermal electrons following a broken power law distribution where, γ is the electron Lorentz factor and, p and q are the low and high energy power-law indices with γ b the Lorentz factor corresponding to the break energy.The emission region is permeated with a tangled magnetic field B and moves down the jet with a bulk Lorentz factor Γ. The broadband SEDs were modeled using synchrotron, SSC, and EC emission mechanisms.XSPEC modeled SEDs of the source PKS 0735+178 for the selected epochs E1, E2 and E3 are shown in Fig. 4. The numerical model is coupled as a local model in X-ray spectral fitting package XSPEC (Arnaud 1996) and the fitting was performed.We found from the fitting, that the γ-ray emission is better explained by the EC process with the target photons from the IR dusty torus at a temperature of 1000 K.The best fit parameters obtained are listed in Table 4.The best fit model along with the observed data are shown in Fig 4 .Though classified as a BL Lac object, the broadband SED of PKS 0735+178 resembles that of a FSRQ with the requirement of EC process to explain the observed high energy emission.This is also similar to the case of the first high energy neutrino blazar TXS 0506+056, which was found not a BL Lac but a FSRQ (Padovani et al. 2019).The requirement of the external target photon field to explain the observed SED in PKS 0735+178 as found here has also been hinted recently (Sahakyan et al. 2023;Acharyya et al. 2023).For all the epochs except E1, the spectral fit satisfies all the γ-ray flux points.The highest γ-ray flux point of E1 corresponding to energy 21 GeV cannot be explained under the leptonic model and we investigated this by considering the additional emission from photo-meson process.In this leptohadronic model, we considered the interaction of relativistic protons with the synchrotron photons.
We considered a particle distribution involving a combination of leptons and hadrons to model the excess high energy emission in E1.The contribution of these protons is subdominant in the lower energy region.We considered a photomeson process in which relativistic protons interact with the synchrotron photons.In this model, we only considered the production of π 0 meson and its decay into two γ photons.We used the model developed by Kelner & Aharonian (2008)to incorporate the photo-meson emission process in the modeling of the SED of E1.The energy distribution of γ-rays produced by the decay of π 0 meson is given by  The epoch E1 is re-modeled by including a photo-meson process and shown in Fig. 5.

SUMMARY
In this work, we carried out γ-ray spectral, timing, and multiband SED analysis of the neutrino blazar PKS 0735+178 based on γ-ray data collected over a period of 14 years.The main motivation is to see how the broadband SED of the blazar during the neutrino detection epoch compares with the SED of the same source generated at other epochs.The results are summarised below • The source was in its historic brightness state in the γray band during E3.During this period, a neutrino event was localized close to the vicinity of PKS 0735+178.During this period, in addition to its brightness state in the γ-ray band, it also flared in the UV, optical and X-ray bands • From variability analysis of the three epochs, we found a maximum Fvar of 95% during the epoch E3.
• On shorter time scales, we found a flux doubling/halving time scale of 5.75 hr, and the size of the γ-ray emission region is estimated to be 8.22×10 15 cm.
• The γ-ray spectra during all three epochs are well fit with the PL model.The photon indices obtained by the PL model fit to all three epochs range between 1.9 and 2. This is similar to that known BL Lacs (Singal 2015), however, harder when compared to FSRQs.
• We found the source to show spectral variations, A moderately softer when brighter trend was noticed in all three epochs.
• The SED modelling of the epochs E2 and E3 leads to the conclusion that the observed γ-ray emission is due to IC scattering of thermal IR photons by relativistic jet electrons.
Even during the epoch of neutrino detection, the SED of PKS 0735+178 can be well modeled under pure leptonic emission scenario.Probably during this epoch the hadronic emission contribution is subdominant to the leptonic process.This is consistent with Sahakyan et al. (2023), where the authors showed that the SED of PKS 0735+178 during the neutrino detection epoch was leptonic dominated.
• The leptonic model alone was unable to reproduce the quiescent state γ-ray spectrum; however, a combination with the photo-meson process suggests significant improvement in the spectral fit.This suggests that the source PKS 0735+178 has a non-negligible hadronic emission contribution during its quiescent state but during the flaring periods E2 and E3 leptonic process dominates over the hadronic process.
With the identification of the association of neutrinos with more blazars, followed by multi band data accumulation and SED modelling, it would be possible to constrain the leptonic and/or hadronic origin of the high energy γ-ray emission in blazars.
The blazar PKS 0735+178 9 and optical from 2008 August to 2022 February (MJD: 54749−59611) are shown in Fig 1.From Fig 1, it is evident that PKS 0735+178 has gone through both quiescent and active phases.During this period, we identified three time intervals.They are denoted as epochs E1, E2, and E3.E1 covers the period MJD 55370−55635 (265 days), when the source was in a quiescent state, E2 covers the period MJD 56025−56150 (125 days), when the source was in an active state and E3 covers the period MJD 59500−59650 (150 days) when the source was in a historically high state coincident with the epoch of neutrino detection.The neutrino event detected on December 8, 2021

Figure 1 .
Figure 1.Multi wavelength light curves of the source PKS0735+178.From the top, the first panel shows the weekly binned γ-ray light curve for the time range MJD 54749−59611; the second panel shows the SWIFT-XRT light curve in photon counting mode, and the third panel shows the Swift UVOT light curves in W1, W2, U, and V-bands.The blue vertical line is the epoch of neutrino detection.The region encompassed by the black vertical lines, red vertical lines and green vertical lines refer to the epochs E 1 , E 2 and E 3 .The inset is the zoomed region around the flaring epoch E 3 .

Figure 2 .
Figure2.Left panels: One week binned gamma-ray light curves for the epochs E 1 (top left), E 2 (middle left) and E 3 (bottom left) in the 0.1 -10 GeV (red) and 10−300 GeV multiplied by a of 20 (green).Right panels: Hardness ratio versus total intensity in the 0.1 -300 Gev band for the Epochs E 1 (top right), E 2 (middle right) and E 3 (bottom right).Here the solid blue lines are the weighted linear least squares fit to the data.

]Figure 3 .
Figure 3. Simple power law (red solid line) and log parabola (green solid line) fits to the γ-ray spectra of PKS 0735+178 during epochs E 1 (top panel), E 2 (middle panel) and E 3 (bottom panel).
and n ph (ϵ) dϵ are the proton and photon number densities in the energy ranges dEp and dϵ respectively.

Figure 5 .
Figure5.Broad band spectral energy distribution along with the one zone leptonic emission model fits for E 1 with the inclusion of the photo-meson process.The synchrotron emission is represented by the black solid line, the EC emission by the red line, the SSC emission by the green line, and the sum is represented by the yellow solid line.The photo-meson process is shown by the pink dotted line.

Table 1 .
The γ-ray fractional variability amplitude (Fvar) in different epochs in the 100 MeV to 300 GeV band.Here, E 1 corresponds to the quiescent state, E 2 corresponds to the active state and E 3 corresponds to the high state of the source that coincides with the neutrino detection.(MJD 59556), is marked as a blue vertical line in Fig 1 in E3.

Table 2 .
Results of the weighted linear least squares fit to the HR v/s total intensity diagram.Here, the epochs E 1 , E 2 and E 3 correspond to the quiescent, active and high state of the source, R is the linear correlation coefficient and p is the probability for no correlation.

Table 3 .
Results of the power law (PL) and log parabola (LP) model fits to the three epochs E 1 , E 2 and E 3 .Here, E 1 , E 2 and E 3 represent the quiescent, active and high state of the source.Here, Gamma is the photon index, Flux is the γ-ray flux value in units of 10 −12 ph cm −2 s −1 , TS is the test statistics, α is the spectral index, β is the curvature in the spectra and T Scurve signifies the presence of curvature in the spectra.

Table 4 .
Results of the broadband SED analysis carried out for the epochs E 1 , E 2 and E 3 .Here, E 1 represents the quiescent state, E 2 , represents the active state and E 3 represents the high state of the source.The parameters p and q are the low and high energy power law indices of the electron distribution, γ b is the the break energy, R is the size of the emission blob in cm, Γ is the bulk Lorentz factor and B is the magnetic field in Gauss.