BROADBAND SPECTRO-TEMPORAL STUDY ON BLAZAR TXS 1700+685

Fermi-LAT is a pair conversion $\gamma$-ray detector in orbit with a field of view (FOV) of about 2.4 sr, which is sensitive to photon energies between 20 MeV - 500 GeV (Atwood et al., 2009). The FOV covers $\sim$ 20% of the sky at any given instant, and it scans the whole sky in every three hours. At 100 MeV, LAT’s single-photon resolution is $<3.5^{\circ}$, but at higher energies ($>1$ GeV), this resolution improves to $<0.6^{\circ}$. In order to identify the potential high energy $\gamma$-ray sources within a given region of interest (ROI), a likelihood analysis is needed to be undertaken using an input model using the gtlike task (Cash, 1979; Mattox et al., 1996). In order to perform our analysis, we have primarily used fermitools¹¹1https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/ package, maintained by the Fermi-LAT collaboration. Usage of fermipy²²2https://fermipy.readthedocs.io/en/latest/quickstart.html was also required for some specific aspects during spectral analysis (Wood et al., 2017).

We collected 0.1–300 GeV Fermi-LAT data within the time-span between MJD 58970 (May 1, 2020) to 59487 (September 30, 2021) for the purpose of our analysis. A circular region of $10^{\circ}$ around the source position was chosen to be the ROI, within which the presence of high energy $\gamma$ ray was searched. As recommended by the Fermi-LAT team in fermitools documentation, the user-specified cut ‘evclass=128’ and ‘evtype=3’ were applied for the purpose of the selection of rows from the input event data and the filter (DATA_QUAL>0)&&(LAT_CONFIG==1) was applied to update the good time interval (GTI) information based on spacecraft specifications. In order to filter any contamination due to the Earth’s limp, a zenith angle cut of $90^{\circ}$ was applied. The latest instrument response function P8R3_SOURCE_V6 was selected for the analysis, and in order to take care of the isotropic and the diffuse background, we used iso_P8R3_SOURCE_V3_v1.txt and gll_psc_v21.xml respectively, available from Fermi Science Support Center (FSSC). The spectral profiles, as well as parameters, were initialized as per the values published in the 4FGL catalog, and beyond 3 degrees of ROI center, all sources were kept frozen. The $\gamma$-ray spectral profile was compared with three different empirical models:

1. 1.

   Power-law profile, defined as

   $$dN(E)/dE=N_{0}(E/E_{0})^{-\alpha},$$

   (1)

   where $E_{0}$ is is the energy scaling factor and $\alpha$ is the spectral index.

2. 2.

   A log-parabola (LP), defined as

   $$dN(E)/dE=N_{0}(E/E_{0})^{-\alpha-\beta\ln{(E/E_{0})}},$$

   (2)

   where $E_{0}$ and $\alpha$ stand as before, and $\beta$ is the curvature index.

3. 3.

   A broken power-law (BPL), defined as

   $$dN(E)/dE=N_{0}(E/E_{\text{break}})^{-\alpha_{i}},$$

   (3)

   with $i=1$ for $E<E_{\text{break}}$ and $i=2$ for $E>E_{\text{break}}$.

TXS 1700+685 was observed seven times between 16th May 2021 and 7th June 2021 by Swift-XRT/UVOT (details of the observations are provided in Table 1) with nonzero exposure. We used FTOOLS task xrtpipeline version 0.13.5 to produce cleaned events. The HEASoft version 6.28 and XSPEC version 12.11.1 have been used for the purpose of our analysis. Cleaned event files corresponding to the Photon Counting mode were considered for analysis and circular regions of radius 20 arcsec centered at the source and slightly away from the source were chosen for the source and the background regions respectively. The source and background spectra were produced using xselect task. The redistribution matrix file (RMF) was obtained from the latest HEASARC calibration database (CALDB) and the ancillary response files (ARF) were produced using the task xrtmkarf. Using grppha, the source event file, background event file, response file, and ancillary response files were tied together and binned in such a way that each bin contains a minimum 10 counts. The spectra were fitted using simple power-law model, with the galactic absorption column density $n_{H}$ kept fixed at $4.0\times 10^{20}~{}\text{cm}^{-2}$.

The UV-optical data were obtained via The Swift Ultraviolet/Optical Telescope (Swift-UVOT, Roming et al. (2005)) which gathers data in six filters: U, V, B, W1, M2, and W2. The source image was extracted from a region of 5 arcsec radius centered at the source and the background image was extracted from a region of 10 arcsec centered away from the source.

We did not find any publicly available contemporaneous radio observation corresponding to the 2021 flare. However, for reference, we add non-simultaneous radio data released by the Planck mission as a part of their second release of the Planck Catalogue of Compact Sources (PCCS). Even though the primary objective of Planck mission is to measure the cosmic microwave background (CMB) anisotropy, as a by-product of their operation detailed catalogue of compact sources as well as galactic diffusion maps are released, and we use the data between 2009 August 12 and 2010 November 27 for the purpose of our analysis (Planck Collaboration et al., 2014, 2016). The data is processed by the low-frequency instrument (LFI; 30–70 GHz) and high-frequency instrument (HFI; 100–857 GHz) Data Processing Centres (DPCs) and is downloaded via the publicly available SED Builder tool³³3https://tools.ssdc.asi.it/ developed at the Space Science Data Center (SSDC).

In Figure-1, we provide the $\gamma$-ray lightcurve for the source in 0.1-300 GeV corresponding to the time between May 1, 2020 - September 30, 2021. In panels (b) and (c), we also show the same light curve for the lower energy band (LEB, 0.1-1 GeV) and higher energy band (HEB, 1-300 GeV). In Section 3.1, we describe the variability properties and the rise as well as decay times of the prominent flares that we calculated. We provide the correlation features between the LEB and HEB in Section 3.2, which hints at possible delays of onset of these flares in different wavebands and in turn can indicate the possible emission mechanisms.

We observe that the entire domain of observation is interspersed with several flaring episodes, but frequent and short-duration flaring events are most prominent during the $\sim$ 150 day span from February 05, 2021 (MJD 59250) to July 05, 2021 (MJD 59400).

In order to quantify the flux variability, we use the fractional variability amplitude (Vaughan et al., 2003) defined as

where $\sigma$ is the standard deviation of the sample of size N and $\sigma_{i}$ is the statistical uncertainty of each of the data comprising the full sample. The normalized excess variance is defined to be $\sigma_{NXS}={\sigma_{XS}^{2}}/{\mu^{2}}$, and the fractional variability happens to be the square root of this normalized excess variance. The error corresponding to the fractional variability is estimated as

where the meaning of the respective quantities are as explained above (Prince et al., 2018).

During the initial phase of low flaring activity (between MJD = 58970 and MJD = 59250), we observe the normalized excess variance in 0.1-300 GeV to be $0.85$, and the fractional variability to be $0.92\pm 0.07$. However, towards the final later segment of observation (beyond MJD = 59250) where frequent flaring activities are noticed, the fractional variability becomes $1.18\pm 0.70$.

The rise and decay times of the flaring events were calculated by fitting the 150 days long $\gamma$-ray lightcurve with sum of exponential profiles of the form

where $F_{0}$ is the baseline flux, $F_{i}$ is the source flux at time $t_{0}$ representing the approximate amplitude of the flare, $T_{r}$ and $T_{d}$ are respectively the rise and decay times of the flares. In HEB, we observe frequent flaring activity, but the photon intensity in this band is at least one order of magnitude lesser compared to broadband flux (0.1-300 GeV). We examine the rise and decay time, as well as the flux doubling time-scale only for the 150 days time window spanning between MJD = 59250 and MJD = 59400 by fitting the broadband flaring region using a combination of sum of exponential profiles, as shown in Figure-2. Between MJD = 59265 and MJD = 59270, the 0.1-300 GeV flux rises from ($3.61\pm 1.52$)$\times{10^{-8}}~{}\text{erg cm}^{-2}\text{s}^{-1}$ to ($2.31\pm 0.57$)$\times{10^{-7}}~{}\text{erg cm}^{-2}\text{s}^{-1}$, and on MJD = 59275 the flux again drops to ($5.17\pm 2.40$)$\times{10^{-8}}~{}\text{erg cm}^{-2}\text{s}^{-1}$. The computation of the flux doubling time was undertaken using the formula

where $F_{t1,t2}$ are the fluxes at time instances $t_{1}$ and $t_{2}$ respectively and $\tau_{D}$ is the flux doubling time. From this formula, we find the flux doubling time for this flare to be 3.22 days. The rise and decay time obtained from the sum of exponential fits were respectively found to be 2.40 days and 0.55 days respectively. Subsequently, between MJD = 59300 and MJD = 59350, we do not observe either a fully quiescent state or a full-fledged flare, and we call this a ‘plateau’ state after (Prince et al., 2017).

The average flux from the source is $(8.10\pm 2.10)\times 10^{-8}~{}\text{erg cm}^{-2}\text{s}^{-1}$, with the flux peaking at $(1.78\pm 0.17)\times 10^{-7}~{}\text{erg cm}^{-2}\text{s}^{-1}$. Subsequent to MJD = 59350, corresponding to which simultaneous X-ray data is available, we observe twin sharp flaring events within a span of $\sim$ 25 days. The peak flux of these individual flaring events were $(6.05\pm 2.24)\times 10^{-7}~{}\text{erg cm}^{-2}\text{s}^{-1}$ and $(5.04\pm 0.50)\times 10^{-7}~{}\text{erg cm}^{-2}\text{s}^{-1}$ respectively. The rise and decay times for the first twin flare appeared to be 1.6 days and 2.9 days respectively, and corresponding to the second twin flare they turned out to be 0.9 days and 1.3 days respectively. The shortest flux doubling time-scale thus turns out to be 0.9 days, which can be used to constrain the physical size of the emission region as well. From the shortest flux doubling time-scale of 0.9 days, we quantify the variability time-scale to be $t_{\text{var}}=\tau_{D}\times\ln{2}=0.63$ days.

We have carried out a temporal analysis of the $\gamma-$ ray emission from a blazar TXS 1700+685 with the Fermi- LAT data from May 1, 2020 to September 30, 2021. In the low-frequency domain, we have detected a sharp peak in the periodogram at $\sim 0.06~{}\text{days}^{-1}$ which has $\sim 3\sigma$ significance (Figure-3(a)). The significance of the periodogram peak is estimated by the method proposed by (Emmanoulopoulos et al., 2013) and simulating 1000 $\gamma$-ray light curve ⁴⁴4https://github.com/samconnolly/DELightcurveSimulation, as elucidated in (Connolly, 2015). However, we also detect another peak at $\sim 0.62~{}\text{days}^{-1}$ which has marginally $>2\sigma$ detection significance (Figure-3(b)).

A variety of mechanisms have been proposed to explain the background physical behavior in blazer $\gamma$-ray emission. These include periodically varying Doppler beaming from precessing jet (Marscher & Gear, 1985; Camenzind & Krockenberger, 1992; Abraham, 2000; Caproni et al., 2013), processed related to the innermost stable orbit of the accretion disc (Broderick & Loeb, 2006; Pihajoki et al., 2013) and several others.

One possible explanation of the QPO of a blazar involves an accretion disc hot-spot orbiting close to the innermost stable circular orbit of the central SMBH (Gupta et al., 2019). These hotspot emissions are quasi-thermal and could be directly observed in the optical/X-ray band. The optical/X-ray emission could generate the variation in seed photons and after external Compton interaction with jet, it can trigger the flux modulation in the $\gamma-$ray lightcurve (Gupta et al., 2017). However, in this case, the expected period in the optical/X-ray band of QPO depends on the mass and spin of the central SMBH. Furthermore, it would not be the same in the gamma-ray band, which would be Doppler boosted. Due to the lack of optical/X-ray data, we can not comment on optical/X-ray variability.

Another potential explanation of the QPO in Gamma-ray energy is instability in the disk/jet system. Magnetically choked accretion flows (MCAFs) in general-relativistic three-dimensional magnetohydrodynamic simulations are observed to produce quasi-periodic fluctuations in the energy outflow efficiency of the relativistic jet close to the black hole (Tchekhovskoy et al., 2011). The dominant mode of these fluctuations has a period of $\sim 70r_{g}/c$ for the spin $a=0.9375$. For a central mass of $\sim 10^{8}$ M_⊙ has a period of $\sim 1$d when the Schwarzschild radius $r_{s}\approx 3\times 10^{11}$m. Slower-spinning black holes are expected to have longer periods McKinney et al. (2012). It is unclear, whether these oscillations can produce an observable signature in the radio jet or not. However, The different time-scales make it improbable that this mechanism may explain the QPO seen in the Gamma-ray light curve in TXS 1700+685.

The quasi periodic flux variation can be induced because of the shock advancement along the internal helical structure of the jet, or the precession or twisting configuration of the jet (Camenzind & Krockenberger, 1992; Mohan & Mangalam, 2015). There could be several underlying mechanisms behind a precessing jet. The modulation of the jet axis could be induced by the Lense–Thirring precession of the disk (Fragile & Meier, 2009). It can also be triggered if the AGN is part of a binary supermassive black hole system (Ackermann et al., 2015; Sandrinelli et al., 2016b; Sandrinelli et al., 2016a). However, the oscillation period arising out of such a mechanism appears to be $\sim$ few years, which has been detected recently in the case of a few AGN sources (Ackermann et al., 2015; Sandrinelli et al., 2016b; Sandrinelli et al., 2016a). In our case, we found a periodicity about $\sim$ 17 days, which is much lower than the expected value in the case of binary SMBH scenario.

However, in the scenario of the helical motion of the plasma blob along the internal helical structure of the jet Mohan & Mangalam (2015), the variation of the Doppler boosting factor in association with the viewing angle of the emission region can induce a significant flux modulation in the high energy domain. Depending on the viewing angle and the boosting factor, the expected QPO may vary from day to month time scale. Our analysis of the gamma-ray emission from TXS 1700+685 suggests a similar time scale expected from this model. The observed modulation of the flux could be the enhanced emission region moving helically within the curved-jet or curved helical jet itself Sarkar et al. (2021). As we have a limited idea about the emission mechanism of blazar and lack of multiwavelength data of this source, it is possible that the observed modulation in gamma-ray flux is due to some other intrinsic effect or combinations of these effects.

We have performed Auto Correlation Function (ACF) for the $\gamma$-ray light curve (0.3 to 500 GeV) using the discrete correlation function (DCF)⁵⁵5https://github.com/astronomerdamo/pydcf/blob/master/dcf.py (Edelson & Krolik, 1988) and $\zeta$-discrete correlation function ($\zeta$-DCF)⁶⁶6https://www.weizmann.ac.il/particle/tal/research-activities/software (Alexander, 2014) for comparison. Both correlation functions are followed the same pattern. The most correlated peak is found at zero time lag which is expected as here we have done the ACF analysis. The 2nd most significant peak is found $82\pm 5$ days on both sides of the main peak. The correlation coefficient value of this peak is about 50 for both peaks. In estimating the significance of the ACF peak, the same procedure as elucidated in the previous section has been followed. We perform the computation of ACF over 1000 $\gamma$-ray light curves and estimated the 1$\sigma$, 2$\sigma$ and 3$\sigma$ for each time lag. The $\sigma$ levels (black and gray curves) are plotted on the correlation function as shown in Figure-4. We have found 82$\pm$5 days peaks are significant up to 3$\sigma$.

In order to estimate the dominance of high energy photons during the flaring region, we carried out the Fermi-LAT analysis with $2^{\circ}$ of ROI. We choose this threshold because the flux of the source was low, the likelihood analysis could not be performed for most of the time intervals with smaller ROI. We do not perform the analysis with the ULTRACLEAN class of events because of the same reason. In Figure-5, we plot the lightcurve corresponding to the high energy ($>1$ GeV) photons detected within $2^{\circ}$ ROI within the flaring phase. In (Figure-5), we show the number of such photons with $>90\%$ probability of detection at successive time intervals.

We observe the highest number of high energy photons (total 42) on MJD = 59370, and the highest energy of these photons being 6.74 GeV with the detection probability 99.80%. The number of such high-energy photons during the pre-flare plateau region is significantly lower, with the highest number being 10 on MJD = 59270. The highest energy of the detected photons with high detection significance corresponding to the entire span of observation pertaining to our analysis, however, turns out to be 25.41 GeV on MJD = 59275 with a detection probability of 99.97%. Such high energy photons could be produced as a consequence of the external Compton scattering of the photons originating from the broad-line region (BLR), or dusty torus by the jets or synchrotron self-Compton emission.

In this section, we describe the key features regarding the $\gamma$-ray spectral energy distribution (SED). We attempt to approximate the shape of the SED using different profiles: power-law (PL), log-parabola (LP), and broken power-law (BPL) as defined in Equations 1-3. The choice of these functional forms is motivated by the spectral fitting over equispaced logarithmic energy bins of blazar flares in 100 MeV - 200 GeV range (Ackermann et al., 2010). In order to determine the quality of the corresponding spectral fits, the Log(likelihood) values obtained from the unbinned likelihood analysis have also been examined. In Figure-6, we show the comparison of different profiles. The spectral index corresponding to the power-law fit turns out to be $2.18\pm 0.03$. However, even by visual inspection the presence of spectral curvature is detected, which motivates us to examine the alternative phenomenological models such as logparabola (LP) and broken power law(BPL). The significance of the spectral curvature is determined via $\text{TS}_{\text{curve}}$, which is defined as the $\text{TS}_{\text{curve}}=2(\mathcal{L}(BPL/LP)-\mathcal{L}(PL))$. The spectral curvature is significant only if $\text{TS}_{\text{curve}}>16$ (Nolan et al., 2012), where L represents the likelihood of the function. In both of the cases (LP and BPL) we found this criterion to be satisfied, implying that spectral curvature is statistically significant. In the case of LP, we find $\alpha=2.26\pm 0.05$ and $\beta=0.17\pm 0.04$, with the $TS_{\text{curve}}=18$. On the other hand in the case of BPL fit, we find the spectral index before the break at 1 GeV to be $1.92\pm 0.07$ and after the break to be $2.63\pm 0.11$, with $TS_{\text{curve}}=34$. Thus, The BPL photon index demonstrates the trend of being $<2$ before the break, and after $>2$ after the break. Such a falling spectrum has earlier been observed in the case of several other FSRQs as well, 3C 454.3 (Abdo et al., 2011), 3C 279 (Paliya et al., 2015; Prince, 2020).

These results can also serve to constrain the possible origin of the $\gamma$-ray emission region. If the power-law could approximate the SED well, it would imply that the emission region is outside the broad-line region (BLR). However, the presence of a significant spectral curvature implies that the photons are further reprocessed in the BLR (like photon-photon pair production), which attenuates the $\gamma$-ray flux leading to a spectral break (Liu & Bai, 2006). The presence of the significant cut-off at 1 GeV explains the non-detection of very high energy photons throughout the period of observation, and also implies that the 25.41 GeV photon detected on MJD = 59275 might have originated from outside of the BLR.

Modeling the multi-waveband SED turns out to be a powerful tool in understanding the flaring mechanism associated with the source under consideration. In our work, we consider a one-zone leptonic model of jet acceleration based on the stochastic or shock acceleration mechanism (Ghisellini & Tavecchio, 2009; Tramacere et al., 2011). We use the publicly available jetset implementation for the purpose of our SED modelling (Massaro et al., 2006; Tramacere et al., 2009, 2011). In this model, it is assumed that a spherically symmetric blob of size $R$ filled with relativistic population of electrons is moving through an entangled magnetic field ($B$). The leptons in the jet are accelerated to an ultra-relativistic regime via a stochastic or shock acceleration mechanism inside the jet. The plasma blobs are assumed to be moving at an angle $\theta$ relative to the observer’s line of sight, implying that the emission will be beamed by a factor $\delta=1/[\Gamma(1-\beta\cos{\theta})]$, where $\Gamma$ is the bulk Lorentz factor of the blob. We assume the relativistic population of electrons within the blob is satisfying an empirical broken power-law distribution of the form,

The first hump of the broadband SED is produced by the synchrotron emission as a consequence of the interaction of relativistic electrons with the magnetic field. These synchrotron photons serve as the seed photons for the inverse Compton process by the same population of relativistic electrons named as synchrotron self-Compton (SSC), and this SSC leads to the high energy hump. Further, several other processes which fall in the category of external Compton (EC) process also contribute to the genesis of this hump. The seed photons for the EC emission come from
(i) direct emission from the accretion disk,
(ii) reprocessed emission from the infrared region of the surrounding dusty torus (DT),
(iii) reprocessed emission in the optical/UV domain from the BLR.

We consider the seed photons from BLR, DT as well as the disk for the purpose of our broadband SED modeling. The thermal emission of the accretion disc is approximated by the multi-temperature blackbody radiation, such that the radial profile of temperature is provided as

where $R_{S}$ is the Schwarzschild radius, $\epsilon$ is the efficiency of accretion, $\sigma$ is the Stefan-Boltzmann constant. The outer boundary of the disc is assumed to be $500R_{S}$. In our SED modeling, we let the parameters regarding the electron distribution as well as the magnetic field vary, while the inner and outer extent of the BLR ($1.0\times 10^{18}$ cm and $2.0\times 10^{18}$ cm respectively) as well as the size of the DT ($5.0\times 10^{18}$ cm) are left frozen. The broadband spectral fitting is shown in Figure-7 and the fitted parameters are provided in Table 3.

The gravitational light-bending effect of the emitted radiation of quasars by the intervening foreground sources has been well recognized in the literature (Torres et al., 2003). For a gravitationally lensed blazar source, the $\gamma$-ray lightcurve is expected to show some time-delay, as has already been found in the case of PKS 1830-211 (Barnacka et al., 2011) and B0218+357 (Cheung et al., 2014). The resolution of Fermi-LAT is not sufficient to resolve the lensed images located within 1 arc second spatial resolution (Pramesh Rao & Subrahmanyan, 1988), but the time-lag between different photon paths detected from the self-correlation of the $\gamma$-ray lightcurves could be a plausible way to detect such lensing effects (Abdo et al., 2015; Abhir et al., 2021a).

In our case, we have detected a strong self-correlation ($>3\sigma$ detection significance) of the 0.1-300 GeV $\gamma$-ray lightcurve at a time-delay of 82 days (Figure-4). We conjecture the presence of such a strong signature of a time delay to be a consequence of the microlensing effect of the $\gamma$-ray photons originating from different emission regions. However, conclusive pieces of evidence can only be gathered from detailed radio monitoring of this source, since it would then be possible to infer the presence and impact of such lensing effects using flux ratios (Lovell et al., 1996). Also, since the emission region corresponding to radio and $\gamma$-ray could be different, emission size induced chromaticity would be an important indicator to probe such effects.

In the broadband SED, each segment indicates some constraints on the allowed parameter space. In the case of blazars, some important constraints are:

1. 1.

   The synchrotron emission component suffers a sharp decline below $10^{12}$ Hz because of the synchrotron self-absorption mechanism, in which the photons interact with the ambient electrons and lose energy. Since the radio data is above the predicted curve and demonstrates only the emission from the extended segment from the jet, radio data provides only a weak constraint on the low energy synchrotron spectrum.

2. 2.

   The optical/UV data is primarily contributed by the synchrotron and the optical emission from the disk, this segment of the SED provides constraints on the magnetic field ($B$), particle density ($N$) and the first spectral index ($\alpha_{1}$).

3. 3.

   The high energy peak contributed by the X-ray/$\gamma$-ray is contributed by the synchrotron self-Compton (SSC) and external Compton (EC) processes. This segment is determined by the second spectral slope ($\alpha_{2}$), the physical distances of BLR and DT ($R_{\text{BLR}}$ and $R_{\text{DT}}$).

4. 4.

   Overall emission is also a function of the bulk Lorentz factor of the jet ($\Gamma$), viewing angle ($\theta$), Doppler beaming factor ($\delta$).

Because of the weak constraint provided by the radio data and given the fact that in the absence of simultaneous radio observation we are using only historical data for referencing purposes, the radio segment of the SED is not well approximated by our model, which is common in Blazar SEDs ( (e.g., D’Ammando et al. (2012)). The fitted spectrum is provided in (Figure-5) and the fitted parameters are listed in Table 3. The lower energy hump is mostly contributed by the synchrotron component, and the contribution from the disk emission is found to be sub-dominant. Reharding the high energy hump, we observe the X-ray segment is well modeled with the SSC component, while the $\gamma$-ray segment is produced by both SSC and EC components. The generation of keV photons in synchrotron emission requires the presence of a strong magnetic field or the injection of a large number of high-energy electron populations. The equipartition value that we have observed in our case implies that the bulk of the power is carried away by the electron population, and the magnetic field carries only a tiny fraction of power. This is reflected from the low magnetic field ($\sim 0.12$ G) we have obtained from our spectral fit. Our estimation of the magnetic field is in agreement with the previous estimation for which the equipartition value is similarly on the higher side (e.g. in PKS 1424-418; Abhir et al. (2021b) for which the estimated magnetic field turns out to be $\sim 0.2$ G), and is at variance with low equipartition value estimation corresponding to stronger magnetic fields (Buson et al., 2014). The bulk Lorentz factor ($\Gamma$) and the Doppler factor ($\delta$) are also significantly lower compared to what was obtained in Buson et al. (2014).

The highest-energy photon obtained during the flaring phase was found to be $6.67$ GeV. This is below the threshold for which the photons become opaque (more than $\sim$ 20 GeV) and are reprocessed by the BLR region. Because of this reason, we included the photons from DT, BLR as well as the disk as seeds for the EC mechanism. We observe the EC contribution by the BLR region to be dominant in the $\gamma$-ray segment, and the EC from the DT to be subdominant. The EC component corresponding to the disk photons seems to be negligible (turns out to be 4-5 orders of magnitude smaller compared to the other two components).

We have estimated the approximate power carried by different components of the jet, as well as the entire jet by the following equation

Using the fitted parameters, we find the jet power contributed by the electrons and protons to be $1.61\times 10^{43}~{}\text{erg s}^{-1}$ and $1.06\times 10^{42}~{}\text{erg s}^{-1}$ respectively. The power deposited by the magnetic field component turns out to be $3.02\times 10^{42}~{}\text{erg s}^{-1}$, which is subdominant compared to the leptonic component. The total jet power turns out to be $2.03\times 10^{43}~{}\text{erg s}^{-1}$.

The flux distribution study can provide insight into the underlying mechanisms producing the variations in blazars. The flux variability in blazars has been observed to be stochastic in nature. Typically, the flux distribution is Gaussian in the case of linear stochastic process, while a lognormal flux distribution (i.e. a Gaussian distribution in the logarithm of flux) is expected for a non-linear stochastic process. Several blazars show lognormal flux distributions using multi-wavelength lightcurves, and at different time-scales (Giebels & Degrange, 2009; H.E.S.S. Collaboration et al., 2010; Tluczykont et al., 2010; Kushwaha et al., 2016; Sinha et al., 2018; Khatoon et al., 2020; Scargle, 2020; Prince et al., 2021). The advent of high sensitive Fermi-LAT (100 MeV - 300 GeV) enables us to characterize the flux distribution study in $\gamma$-ray energy. Several blazars seen by Fermi-LAT, show lognormal behavior in their long-term $\gamma$-ray lightcurves (Ackermann et al., 2015; Shah et al., 2018; Romoli et al., 2018; Bhatta & Dhital, 2020). The lognormal variability in blazar is usually explained by assuming the jet emission arising from non-linear multiplicative processes originated in the accretion disc. However, minute-scale variability as observed in many blazars may not be produced by the fluctuations in the disc, and that favors the variability to produce in the jet itself. The “mini jets-in-a-jet" statistical model provides a possible explanation where the addition of emission from a large collection of mini jets can explain the lognormal flux distribution (Biteau & Giebels, 2012). Furthermore, Sinha et al., 2018 showed that the lognormal flux distribution can be obtained when the underlying particle acceleration process is time-dependent. In such cases, a Gaussian perturbation in the particle acceleration time-scale can give rise to a lognormal flux distribution and linear flux-rms correlations at high energies.

In this work, we presented the flux distribution study for the source TXS 1700+685, with three days binned gamma-ray flux lightcurve. To select the statistically significant lightcurve, we considered the flux points for which flux is greater than 2-$\sigma$ level, such that ${\frac{F}{\Delta{F}}}$ > 2. We produced the normalized histogram of the logarithm of flux and fit it with the Gaussian and lognormal probability density functions (PDFs, Shah et al. (2018)) (Figure-8). It is found that the lognormal function fits the distribution significantly better with reduced chi-square (${\chi^{2}}_{red}){\approx}$ 0.91 for 11 degrees of freedom (dof), than the Gaussian function (${\chi^{2}}_{red}{\approx}$ 2.29 for 11 dof). The linear dependence of the mean flux $<F>$ with its excess root mean square (rms) variability, is an essential feature for the lognormal process. We estimated the excess variance (Vaughan et al., 2003), $\sigma_{XS}=\sqrt{\sigma^{2}-\overline{\sigma_{err}^{2}}}$ with $\sigma^{2}$ and $\overline{\sigma_{err}^{2}}$ represent the variance and the average of the square of the measurement errors, using 50d binned data. The rms-flux plot is shown in Figure-9, fitted with a linear function of slope 0.57$\pm$0.1. The Spearman’s rank correlation coefficient ($r_{s}$) and the correlation probability (p) are obtained between $<F>$ and $\sigma_{XS}$, which results $r_{s}$ = 0.78 and p = $4.5\times 10^{-3}$, suggesting linear rms-flux relation. The obtained lognormal feature in the flux distribution and the linearity between the excess variance to the flux implies the variation in flux is lognormal.