Estimating the Jet Power from Broadband SED modeling of Mkn 501 for different particle distributions

We consider the broadband spectral energy distribution of the high energy peaked (HBL) blazar Mkn 501 using $\textit{Swift}$-XRT/UVOT, NuSTAR and $\textit{Fermi}$-LAT observations taken between 2013 and 2022. The spectra were fitted with a one-zone leptonic model using synchrotron and synchrotron self-Compton emission from different particle energy distributions such as a broken power-law, log-parabola, as well as distributions expected when the diffusion or the acceleration time scale are energy dependent. The jet power estimated for a broken power-law distribution was $ \sim 10^{47} (10^{44})$ erg s$^{-1}$ for a minimum electron energy $\gamma_\text{min} \sim 10 (10^3)$. However, for electron energy distributions with intrinsic curvature (such as the log-parabola form), the jet power is significantly lower at a few times $ 10^{42}$ erg s$^{-1}$ which is a few percent of the Eddington luminosity of a $10^7$ M$_\odot$ black hole, suggesting that the jet may be powered by accretion processes. We discuss the implications of these results.


INTRODUCTION
Blazars are jetted objects, the subclass of Active Galactic Nuclei (AGN) where the jet is directly pointed towards the observer (Urry & Padovani 1995b).Significant properties include its Spectral Energy Distribution (SEDs) ranging from Radio -TeV gamma rays, non-thermal radiation, rapid variability, radio-loud features, and double hump structure (Urry 1998;Massaro et al. 2004a).Blazars are broadly classified into BL Lacs and Flat Spectrum Radio Quasars (FSRQs) based on their optical emission lines.FSRQs show strong emission lines in the optical band whereas BL Lacs show faint emission lines (Urry & Padovani 1995a;Antonucci 1993).Blazars are also classified into three classes based on the peak of the emission of the synchrotron component.For low-Energy peaked (LBL)/Low -Synchrotron peaked (LSP) BL Lacs, the low-energy hump peaks at the optical/UV band.For High energy peaked (HBL)/high synchrotron peaked (HSP) BL Lacs, the component peaks in the X-ray band and the Intermediate peaked BL Lacs (IBL) have a peak in between the X-rays and optical/UV band (Ghisellini et al. 1997;Fossati et al. 1998).The broadband SED is composed of two peaks with the lower one attributed to synchrotron emission.The radiative processes responsible for the second component in the high-energy regime are still unresolved.In the leptonic model interpretation, the component is assumed to be due to inverse Compton (IC) scattering (Dermer et al. 1992;Dermer & Schlickeiser 1993).The seed photons ★ E-mail: hritwikbora@gmail.com† E-mail: rukaiyakhatoon12@gmail.com ‡ E-mail: rmisra@iucaa.in§ E-mail: rupjyotigogoi@gmail.com for the IC processes may either be the synchrotron photons (i.e. the Synchrotron Self Compton (SSC) process) or they may be external photons from outside the jet (Marscher & Gear 1985;Dermer & Schlickeiser 1993;Khatoon et al. 2022).
One of the main challenges in understanding these systems, is the identification of the energy-generating process which results in these jets having such high powers, as inferred from broadband spectral fitting.For example, Paliya et al. (2017) have fitted the spectra of several blazars (both -ray loud and quiet) and have inferred total jet powers with an average value of few times 10 47 erg s −1 and with some sources exceeding 10 48 erg s −1 .Such jet powers exceed the Eddington luminosity of super-massive black holes and hence they cannot be powered by accretion and instead require alternate energy production processes (e.g., Ghisellini et al. (2014)).In this interpretation, most of the jet power is in the kinetic energy of the protons whose number density is assumed to be the same as that of the non-thermal electrons.The inclusion of electron-positron pairs may decrease the power requirement, but it is unlikely that the pair fraction is substantial due to Compton drag effects and unobserved spectral signatures (e.g.Sikora & Madejski 2000;Celotti & Ghisellini 2008;Sikora et al. 2016).If instead of a leptonic origin for the VHE, one considers a hadronic one, then the jet power required is even higher and can exceed 10 49 erg s −1 (e.g.Zdziarski & Böttcher 2015;Abe et al. 2023).
Since the jet power is dominated by the number of protons, the inferred jet power depends on the number of relativistic electrons.Typically the energy distribution of electrons is assumed to be a broken power-law form and the number of electrons then depends on the lowest lorentz factor of the distribution i.e.  min , which is often assumed to be of the order of unity or tens.Employing large values of  min ∼ 1000, as permitted by the data can lead to a reduction in the jet power.However, the rationale behind having the minimum energy of the electrons to be that high remains unclear.Even then, assuming  min = 1000, for the low flux state of Mkn 501, the jet power can still be as high as ∼ 10 44 erg s −1 (e.g.Abe et al. 2023).
The jet power estimations above have been primarily based on the assumption that the underlying electron energy distribution is a broken power-law (or a smoothened broken power-law).There have been some evidences based on the X-ray and -ray data that the observed spectrum has significant curvature deviating from a power-law (Massaro et al. 2004a;Tanihata et al. 2004;Tramacere, A. et al. 2007).This may indicate that the underlying electron energy distribution has a curvature which may have the shape of a logparabola distribution (Massaro et al. 2004a).Physically motivated models have also been used such as distributions where the high energy cutoff is due to radiative losses (-max models), or where the curvature is due to energy-dependent acceleration (EDA) or energydependent diffusion (EDD) (Sinha et al. 2017;Goswami et al. 2018;Hota et al. 2021;Khatoon et al. 2022).Typically, these different interpretations are often spectrally degenerate, but perhaps may be distinguished by studying how physical are the observed correlations between the best fit spectral parameters for different models (Hota et al. 2021;Khatoon et al. 2022).
Electron energy distributions which are different from a broken power-law one, are expected to provide different estimates of the jet power.The motivation of this work is to estimate these different jet powers with statistical errors, by fitting a set of broadband SED of a bright blazar, with models having different electron distributions.To do this, we consider the HBL source Mkn 501 at a redshift of z = 0.034 which is among the most well studied blazars with VHE emission.It shows extreme HBL characteristics during flaring as reported by Ahnen, M. L. et al. (2018) and Sahu et al. (2020).It was first detected in 1995 by the Whipple observatory (Quinn et al. 1996).It is one of the brightest blazar sources in the catalog after Mkn 421, and has been studied through several multi-wavelength missions.Moreover, Mkn 501 shows significant curvature in the SED as investigated by Tavecchio et al. (2001), which may be described by a log parabola distribution (Massaro et al. 2004b).As reported by Abe et al. (2023), during 15 years of swift observations, the Swift-XRT flux rate was high from 2012 to 2014, with the highest X-ray flaring activity during March-October 2014.A recent multiwavelength study by Furniss et al. (2015) with NuSTAR observations reported hard X-ray variability in the timescale ∼ 7 hours in the range from 7-30 keV which was not previously observed.For this work, we analyzed multi-wavelength observations from Swift-XRT/UVOT, NuSTAR, and Fermi-LAT during the period from 2013 to 2022 covering the high and historically low activity state of Mkn 501, giving ten near simultaneous broadband SED, which we fit using different models.The section 2 describes the observations and the data analysis technique.In section 3, the broadband spectral analysis is reported, whose results are presented in section 4. We summarise and discuss the results in the section 5.In this work, the cosmological constant is considered in accordance with the ΛCDM model with  0 = 71 km s −1 Mpc −1 , Ω Λ = 0.73, Ω  = 0.27 (Komatsu et al. 2011).

NuSTAR
Nuclear Spectroscopic Telescope Array (NuSTAR) is a space-based X-ray telescope used for nuclear spectroscopy.It was launched by NASA on 13 June 2012 which uses a Wolter Telescope.It consists of two Focal Plane modules FPMA and FPMB (Harrison et al. 2013).Several observations of Mkn 501 were obtained between 2013 to 2022 which are summarized in Table 1.
The software package, "nupipeline" which is integrated with Heasoft-6.28,has been used for the generation of cleaned data files for the analysis."XSELECT V2.4k" was used for plotting the event fits files (cl.evt) in ds9.Source extraction of 30 arcseconds and background of 60 arcseconds had been carried out for all the NuSTAR observations.The "nuproducts" tool was used for the generation of spectrum files (PHA), ancillary files (ARF), and response matrix files (RMF).The "grppha" tool was used to obtain the grouped spectra with all the obtained products to have counts of 30 per bin.

Swift-XRT
We considered observations undertaken by Swift-XRT which are nearly simultaneous with the NuSTAR ones.The detector has an area of 135 cm 2 covering an energy range of (0.2 -10 keV).Observations for the instruments NuSTAR and Swift-XRT are provided by NASA's HEASARC archive 1 .Corresponding to the NuSTAR observation during March 2022, there are two Swift-XRT observations for which the analysis was carried out separately i.e. without summing the two Swift-XRT spectra.The observation details are listed in Table 1.
We have used the software package "xrtpipeline" which is integrated with Heasoft-6.28,for the generation of cleaned data and used "XSELECT V2.4K" for plotting the event fits files (cl.evt) in ds9.For the observations taken during 2013, 2021, and 2022 the Swift-XRT counts were comparatively higher than the ones during 2017 and 2018.Among the 10 observations, 9 observations of Swift-XRT are in the Windowed Timing (WT) mode, which are free from pile-up.For all the non-piled-up observations, a source region of 10 pixels and a background region of 20 pixels are extracted.However, "Obs ID 00035023170" was found to be piled up since it is in the "photon counting" (PC) mode.For this observation, we extracted an annular region with an inner radius of 7 pixels and an outer radius of 18 pixels, and a background region of 36 pixels which was free from source contamination.Using "XSELECT V2.4k", the light curve and spectrum are extracted, and using "xrtmkarf" and "quzcif" tools the ancillary (ARF) and response matrix files (RMF) were generated respectively.Like NuSTAR, Swift-XRT spectrum was grouped such that there were at least 30 counts per bin.

Swift-UVOT
The UV observations for the data sets were taken from Swift-UVOT.The Swift-UVOT covers both the optical and UV bands of the EM spectrum and consists of three optical filters (B, V, U) and three UV filters (uvw1, uvw2, uvm2) (Roming et al. 2005).However, for the observations under consideration, only the UV filters were available.The source and background regions of 8 arcsecs and 16 arcsecs were extracted.The tool "UVOTSOURCE" was used to calculate the magnitude in the AB magnitude system.We have used the relation    (− ) = 3.1 for the magnitude correction, where A  is the galactic extinction and  ( − ) is the interstellar reddening which has been taken to be 0.017 mag (Schlafly & Finkbeiner 2011).
From the magnitude values, the flux was calculated with the help of photometric zero points and conversion factors adopted from Breeveld et al. (2011); Larionov et al. (2016).The flux points and energies of UVOT observations were converted to "XSPEC version 12.11.1"(Arnaud 1996) readable (PHA) format using the tool "ft-flx2xsp".

Fermi-LAT
The Large Area Telescope onboard the Fermi Gamma-ray Space Telescope Fermi-LAT, operates in the energy range from 20 MeV to more than 300 GeV (Atwood et al. 2009).In order to get a significant detection and spectra, we considered the 4 months Fermi-LAT observations keeping the NuSTAR observations in the middle of the total duration.We analyzed the data spanning from February 2013 to May 2022 using the open-source tools fermipy2 .We considered the region of interest (ROI) to be 15 • to carry out the analysis and followed the standard procedure as mentioned in the documentation (Wood et al. 2017).The energy range considered was from 100 MeV to 300 GeV with evclass=128 and evtype=3.The latest instrument function (IRF) "P8R3 SOURCE V3" has been used.To generate XML files, the galactic diffusion "gll iem v07.fits" and the isotropic background model "iso P8R3 SOURCE V3 v1.txt" have been used.Using the latest Fermi-LAT 4FGL catalog (Abdollahi et al. 2020), we have acquired the XML model file encompassing all sources within the ROI.Along with it, the fourth Fermi source catalog providing the spectral models and parameters for the sources located within the ROI are considered.Finally, the flux points and energies of Fermi-LAT observations are converted to XSPEC readable (PHA) format using the tool "ftflx2xsp".

BROADBAND SPECTRAL ANALYSIS
The source Mkn 501 is well known to have the X-ray emission from synchrotron radiation, however, the hump in the higher energy is attributed mainly to the inverse Compton (IC) process.In our broadband spectral fitting of Mkn 501, we have considered the leptonic model with a single emitting region.We assume a spherical region having a radius R, filled with a tangled magnetic field, B, and isotropic electron distribution, n(), moving with bulk lorentz factor Γ at an angle  with respect to the observer's direction.The seed photons of the IC process are the synchrotron photons i.e. the synchrotron self-Compton (SSC) process.We fit the broadband using a local model in XSPEC as described below.For the spectral fit, we have used the "Tbabs" model (Wilms et al. 2000) to account for the Galactic absorption.The hydrogen column density for the X-ray observations,  H = 1.69 × 10 20 cm −2 was considered and kept fixed, which was obtained in the LAB survey (Kalberla et al. 2005).However, for the UV the  H value was fixed at 0 since the UV flux was dereddened initially.For error calculation, we have adopted Markov Chain Monte Carlo (MCMC) method using the python-based code XSPEC_EMCEE3 .We have used 50 walkers with 5000 iterations and burnt the first 2000.The errors were calculated within a confidence range of 2.EMCEE is an extensible, pure-Python implementation of Goodman & Weare's Affine Invariant Markov chain Monte Carlo (MCMC) Ensemble sampler.
Since the objective is to fit the broadband spectra with different electron energy distributions, (), it is convenient to use a convolution model in XSPEC such that it can be used along with a function that represents the electron distribution.In particular, we follow Hota et al. (2021); Khatoon et al. (2022) to represent the synchrotron flux at an energy  to be, where   is the luminosity distance and V is the volume of the emission region, A = and  () is the synchrotron emissivity function.Instead of using the electron's lorentz factor , the electron energy distribution is represented by ().The transformation is given by  =  √ C, where C = 1.36 × 10 −11  1+ keV with  as the doppler factor, B as the magnetic field, and z being the redshift of the source.Note that  represents  and  2 has dimension of keV.Similarly, the synchrotron self-Comptonization flux is determined as a function of  syn () and ().The expressions used for the flux are given in Sahayanathan et al. (2018) and references therein, with  transformed to .Thus we have a local XSPEC convolution model, sscicon.which can be convolved with any particle distribution using the XSPEC model form (sscicon*n()).The convolution model sscicon has the following parameters, bulk lorentz factor Γ, redshift , magnetic field , viewing angle , and size .The other parameters including normalization , arise from the specific particle distribution form as described below.The model has the provision to have as a parameter, the jet power  j , instead of normalization, .We describe the different particle distributions used in this work below.

Logparabola model
We consider the particle density to be described by the log parabola function, Here, particle spectral index is  at the reference energy   , while  and  are the spectral curvature parameter and the normalization.The reference  2 was fixed at 1 keV during the spectral fit and the other three parameters ,  and norm  were kept free.

Broken Power Law
We considered a Broken power law (BPL) distribution of the particles to fit the spectrum.
Where  break is the break energy, p is the electron spectral index for  <  break and q is the electron spectral index for  >  break and the transformation is given by  =  √ C. We consider the broken power-law particle distribution to exist only between a minimum and maximum  i.e. () is non-zero for  min <  <  max , where

Particle distribution with a maximum electron energy
We consider the case where there is acceleration and electrons lose energy through the radiative processes.The evolution of the particles Where the  acc and   are the acceleration and escape time scales of the electrons.  is the cooling term is given by, 8    2 , where B is the magnetic field,   is the Thomson's cross-section, and   is the mass of the electron.Q is the mono-energetic injection of electrons at minimum energy  0 .The solution is given by, (5) where, The particle distribution exists only for  >  min , where  min =  min √ C .

Energy dependent time-scale models:
We also employed models incorporating energy-dependent escape or acceleration timescales to provide a possible interpretation of the curvature in the observed spectrum (Hota et al. 2021;Khatoon et al. 2022).

Energy dependent diffusion model (EDD)
The diffusion takes place within a region containing a tangled magnetic field, causing the escape time scale to depend on the electron's gyration radius.In this context, the escape time scale ( esc ) maybe energy dependent as where,  , corresponds to   , when the electron energy is    2 , and  is the index of the assumed power-law energy dependence.As   cannot be greater than the free streaming value  , which restricts the equation to  <   , where   denotes the energy at which this limit is obtained.If we make the assumption that   is significantly greater than any  of interest and ignore the synchrotron losses, we obtain the following electron energy distribution, where, It is convenient to recast the distribution as follows, with , , and  representing the free parameters.It can be shown that, and that the normalization is given by,

Energy dependent acceleration (EDA) model
We consider the prospect that acceleration time-scale is energydependent and we assume its functional form to be given by, Where,  , corresponds to   , when the electron energy is    2 and  is the index of the energy dependence.Similar to the case of the EDD model, the corresponding electron energy distribution is given by, which can be recast as, where,

𝑅
K is the normalization which is given by,

RESULTS
We performed the fitting of the observations listed in Table 1 using the various particle distributions discussed in the preceding section.We fixed the bulk lorentz factor, Γ = 15, and the viewing angle,  = 0, and allowed for the other parameters to vary.Note that for synchrotron and synchrotron-self Compton processes considered here, the relevant parameter is the doppler factor  = 1/(Γ(1 − cos) which for  = 0,  = 2Γ.
For observations made during 2013 (listed as S 2 and S 3 in Table 1), none of the models considered provided a reasonable fit.This was because the Swift-XRT spectra seemed to have a high energy excess, incompatible with the NuSTAR spectra.We were unable to find the reason for this and diferred a more detailed investigation of this to a later work.Thus, for these two observations, the Swift-XRT spectra were omitted.Using the different particle distributions, we achieved reasonable fits for all ten broadband spectra, and the corresponding best-fit parameters are listed in Table 2.A representative spectrum for the observation during MJD 57870 ( listed as S 4 in Table 1) is shown in Figure 1, with the different best-fit model spectra overlayed.For some of the observations, especially the latter ones during 2022, the broken power-law distribution provides a better fit than the others.We note that the majority of the  2 contribution originates from the Xray spectra and in general for all the particle distribution considered, broadband SED is reasonably represented.
The estimated jet power differs significantly for the different particle distributions as shown in Figure 2 where the jet power is plotted against MJD.The models that have intrinsic curvature in the particle energy distributions (i.e.log-parabola, EDD and EDA) require nearly two orders of magnitude less jet power than the broken power-law and gamma-max distributions.For the latter two distributions, the jet power depends sensitively on the minimum energy of the electrons,  min .Figure 3 shows the variation of  2 (top panel) and jet power as a function of  min for the broken power-law distribution for a representative data set.As expected the jet power decreases with increasing  min for large values, the  2 increases because the UV flux is no longer consistent.For a value as large as  min ∼ 10 3 , the jet power for the broken power-law model is ∼ 10 44 erg s −1 , which is still significantly larger than the power estimated for the log-parabola model in Figure 2.
In the above analysis, it was assumed that the bulk lorentz factor Γ = 15.For a representative data set, the top panels of Figure 4 show the variation of  2 as a function of Γ for the broken powerlaw (left panel) and log-parabola (right panel) models.The bottom panels show the dependence of jet power with Γ.As can be seen, Γ is not constrained since the  2 does not change with increasing Γ.The jet power increases with Γ but only marginally, a factor of two increase when Γ increases from 10 to 45.Thus the broad estimates of the jet powers presented here with Γ = 15 are not sensitive to that assumption.

SUMMARY AND CONCLUSIONS
We collated ten broadband spectral energy distributions for Mkn 501, spanning 2013 to 2022, using data from Swift-UVOT, Swift-XRT, NuSTAR, and Fermi-LAT.We fitted the broad band spectra using a single zone leptonic model taking into account synchrotron and synchrotron-self Compton emission from a non-thermal electron distribution.Different electron energy distributions were considered such as a broken power-law, power-law with a maximum energy due to radiative cooling (-max model), log-parabola, energy-dependent diffusion time-scale (EDD) and energy-dependent acceleration timescale (EDA).We determined that all these models adequately explained the broadband spectra and may be considered spectrally degenerate.
We estimated the total jet power with errors for the different observations and for different particle energy distributions.We found as reported before, that the jet power for broken power-law and -max models, was as high as ∼ 10 47 erg s −1 for a minimum lorentz factor  min = 10 and reduced to ∼ 10 44 erg s −1 for  min = 10 3 , which was the maximum value allowed for spectral fitting.On the other hand, for the other energy particle distributions with intrinsic curvature, the jet power turned out to be ∼ 10 43 erg s −1 , independent of  min .We show that these estimates for jet power are nearly independent of the chosen value of the bulk lorentz factor Γ.
To illustrate the reason for the difference in jet power among various particle distributions, we plotted Figure 5 that shows the particle distributions  2 () as a function of  for the broken power-law, log-parabola, EDA, and EDD models for observation S 3 .The peak of the function indicates the energies at which most of the particle energy dominates.In the case of the broken power-law distribution, the peak occurs at  min and is significantly higher than for the other distributions, resulting in a higher jet power.It's noteworthy that the X-ray emission at ∼ 1 keV arises from electrons with  values of ∼ 1.6 × 10 6 , 7.5 × 10 4 , 1.25 × 10 5 , and 1.25 × 10 5 for the broken power-law, log-parabola, EDA, and EDD models, respectively.
The comparatively lower jet power estimated for the intrinsically curved particle energy distributions (especially the log-parabola one) of ∼ 10 43 erg s −1 makes it about a few percent of the Eddington luminosity for a 10 7 M ⊙ blackhole,  Edd ∼ 10 45 erg s −1 .This has an important ramification that such a jet may be directly powered by the accretion process, which in turn may provide a framework to understand the jet-disc connections inferred from the correlation between disc and jet parameters in non-beamed systems (Cao & Jiang 1999;Ho & Peng 2001;Ghisellini et al. 2014;Sbarrato et al. 2016).Furthermore, the results will have a bearing on galaxy evolution studies which invoke AGN feedback through a powerful jet (Fabian 2012;Cielo et al. 2018).A low jet power also indicates a smaller reservoir of energy in the jet, which may be depleted due to shocks and radiative losses in a shorter time scale.This provides an incentive to make a detailed study of the time evolution of such systems that have intrinsically curved particle distributions and hence lower jet power.
The intrinsically curved particle distributions assumed in this work, are empirical in the sense that the curvature is due to an assumed energy dependence of the diffusion or acceleration timescales.Clearly, a more detailed study considering the micro-physics of the processes needs to be undertaken to validate the assumed dependence.Finally, the analysis needs to be extended to different blazars, especially other kinds of blazars, such as FSRQs, to get a more comprehensive picture.
edges the financial support of the NWU PDRF Fund NW.1G01487, the South African Research Chairs Initiative (grant No. 64789) of the Department of Science and Innovation and the National Research Foundation of South Africa.H.B. and R.G. would like to acknowledge IUCAA for their support and hospitality through their associateship program.

Table 1 .
Summary of NuSTAR and Swift-XRT/UVOT observations of the source Mkn 501.