Multi-wavelength variability and broadband SED modeling of BL Lac during a bright flaring period MJD 59000-59943

We carried a detailed temporal and spectral study of the BL\,Lac by using the long-term \emph{Fermi}-LAT and \emph{Swift}-XRT/UVOT observations, during the period MJD\,59000-59943. The daily-binned $\gamma$-ray light curve displays a maximum flux of $1.74\pm 0.09\times 10^{-5} \rm ph\,cm^{-2}\,s^{-1}$ on MJD\,59868, which is the highest daily $\gamma$-ray flux observed from BL\,Lac. The $\gamma$-ray variability is characterised by power-spectral-density (PSD), r.m.s-flux relation and flux-distribution study. We find that power-law model fits the PSD with index $\sim 1$, which suggest for long memory process at work. The observed r.m.s.-flux relation exhibits a linear trend, which indicates that the $\gamma$-ray flux distribution follows a log-normal distribution. The skewness/Anderson-Darling test and histogram-fit reject the normality of flux distribution, and instead suggest that the flux distribution is log-normal distribution. The fractional-variability amplitude shows that source is more variable in X-ray band than in optical/UV/$\gamma$-ray bands. In order to obtain an insight into the underlying physical process, we extracted broadband spectra from different time periods of the lightcurve. The broadband spectra are statistically fitted with the convolved one-zone leptonic model with different forms of the particle energy distribution. We found that spectral energy distribution during different flux states can be reproduced well with the synchrotron, synchrotron-self-Compton and external-Compton emissions from a broken power-law electron distribution, ensuring equipartition condition. A comparison between the best fit physical parameters show that the variation in different flux-states are mostly related to increase in the bulk-Lorentz factor and spectral hardening of the particle distribution.


INTRODUCTION
Blazars are extreme class of extragalactic sources with powerful relativistic jet pointing close to the line of sight of observer (Urry & Padovani 1995).The relativistic jet of nonthermal plasma originates near to the supermassive black hole and is launched to kiloparsec/mega parsec scales (see e.g.Blandford et al. (2019)).The close pointing of relativistic jet produces Doppler boosted emission, which results in extreme observational properties.Blazars show variability in flux and spectrum across the electromagnetic band, the timescales of variability can go down to few minutes ( see e.g, Aharonian et al. 2007;Ackermann et al. 2016) to years.The detailed analysis of these observational features often help us to obtain insight into the underlying physics and structure of the emission region.For example, the short timescale variability in flux hints the emission occurs in the ⋆ E-mail: shahzahir4@gmail.cominner region of jet near to the SMBH, such regions are difficult to explore otherwise with the present imaging telescopes.Blazars are broadly divided into two classes: BL Lac objects and Flat Spectrum Radio Quasars (FSRQ).This divide is usually based on the characteristics of their optical spectrum.BL Lac objects have weak or absent emission line features with equivalent width (EW < 5 Å), while FSRQs show prominent emission line features with EW > 5 Å (Urry & Padovani 1995).
The broadband spectral energy distribution (SED) of blazars exhibit a characteristic double peaked spectrum.The low energy component peaking at optical/UV/X-ray energies is due to Doppler boosted synchrotron emission.While the high energy component which peaks at γ-ray energy is explained in different scenarios.The most commonly considered scenario is the inverse-Compton (IC) process, where either the synchrotron photons or external photons or both gets up-scattered to high energies by the relativistic particles in the jet.In case the seed photons for the IC scattering are synchrotron photons, then the process is known as synchrotron self Compton (SSC: Jones et al. 1974;Maraschi et al. 1992;Ghisellini et al. 1993), if the photons are external to jet then the process is known as external Compton (EC: Dermer et al. 1992;Sikora et al. 1994;B lażejowski et al. 2000;Shah et al. 2017).Another possible explanation for the high energy component is through the proton-synchrotron process or proton-photon interactions (Mannheim & Biermann 1992;Mannheim 1993;Mücke & Protheroe 2001;Mücke et al. 2003).These interactions are likely responsible for producing not only γ-ray photons but also high energy neutrinos.In fact, TXS 0506+056 was the first flaring blazar to be associated with the high energy neutrino IC-170922A (IceCube Collaboration et al. 2018).Since the detection of neutrino event from blazars, the hybrid models i.e., leptohadronic models have become more attractive for modeling the broadband SED of blazars (e.g.Ansoldi et al. 2018;Sahakyan 2018;Keivani et al. 2018;Gasparyan et al. 2022).Blazars are further classified based on peak frequency of the synchrotron component (ν p syn ) into low synchrotron peaked blazars (LSP; ν p syn ≤ 10 14 Hz), intermediate synchrotron peaked blazars (ISP; 10 14 ≤ ν p syn ≤ 10 15 Hz), and high synchrotron peaked blazars (HSP; ν p syn ≥ 10 15 Hz) (Abdo et al. 2010) .
BL Lacertae (BL Lac) is a blazar located at a redshift of z ∼ 0.069 (Miller et al. 1978).The source is generally classified as LSP (Nilsson et al. 2018), but has also been identified as ISP (see eg. Ackermann et al. 2011;Hervet et al. 2016).BL Lac falls in the category of TeV detected blazars (Neshpor et al. 2001;MAGIC Collaboration et al. 2019), first VHE emission detected from source was reported in Neshpor et al. (2001).The sources had been reported in enhanced flux states over different energy bands by number of Astronomy telegrams (some recent telegrams are Buson & de Menezes 2021;Blanch et al. 2021;La Mura 2022;Prince 2022).These made it a target of number of multi-wavelength observational studies (see e.g Gaur et al. 2015;Wierzcholska et al. 2015;Abeysekara et al. 2018;MAGIC Collaboration et al. 2019;Prince 2021;Sahakyan & Giommi 2022) and gave us broader understanding of its emission properties.Moreover, a wide range of interpretations have been used to explain the variable emission of the source in different flux states.For example, the low-activity state emission from the source with moderate correlated variability between UV and X-rays bands are explained by the SSC and EC processes (Abdo et al. 2011).A broadband SED modelling performed on the sources in the low and high flux states showed that the two emission regions located at two different sites are required to explain the high energy emission in the source (Prince 2021).These studies suggest that the broadband SEDs of BL Lac object cannot be well described by a one-zone SSC model, and instead complex models involving the multi-zone SSC or EC components are needed.Raiteri et al. (2013) explained the variable emission from the source due to change in orientation of the emitting regions, while Wehrle et al. (2016) used an extended data and explained the flaring emission in the source as turbulent plasma flowing across quasi-stationary shocks.Using the longterm optical photometric and polarimetric data of BL Lac, Blinov & Hagen-Thorn (2009) showed that the observed variability of BL Lac can be explained with a steady component of high degree of polarisation ∼ 40%.Giebels & Degrange (2009) used the X-ray light curve of BL Lac to investigate the nature of its variability.They showed that the X-ray light curve of BL Lac follows the lognormal distribution, infact BL Lac was the first blazar in which lognormality was detected (Giebels & Degrange 2009).The X-ray light curves used in their work were less variable compared to other blazars.Moreover, the amplitude of variability was found to be proportional to the flux level.The lognormal variability and linear r.m.s-flux relation in the X-ray light curves had been observed in various compact systems, such as Seyfert galaxies and X-ray binaries.In these sources, the variability in the emission is a result of fluctuations in the accretion disk (Uttley et al. 2005;McHardy 2008), these fluctuation propagate inwards and produces a multiplicative emission.This implies that the lognormal distribution of flux is being powered by multiplicative processes rather than additive.In case of blazars, the observation of lognormal distribution possibly imply the disc variability is imprinted on the jet, hence provides a possible link between the accretion and jet properties in a blazar.
Number of multiwavelength works have been carried out to understand the possibly physical scenario responsible for the variable emission (Gaur et al. 2015;Wierzcholska et al. 2015;Abeysekara et al. 2018;MAGIC Collaboration et al. 2019;Prince 2021;Sahakyan & Giommi 2022).These studies suggest complex physical mechanisms are required for the emission in low and high flux states.In this work, we conducted a comprehensive multiwavelength study of BL Lac by utilizing more than two and half years of Fermi data, which includes the most recent and brightest γ-ray flare ever detected from the source.The brightest γ-ray flare has not been studied in previous works.We conducted a comprehensive analysis of the γ-ray flux distribution of the BL Lac source, which has not been explored in earlier works.We also developed a convolved one zone leptonic model and incorporated it as a local convolution model in XSPEC.The model provides the statistical broadband fit to the observed SED for any input particle distribution.In this work, using the convolved one-zone leptonic, underlying physical parameters responsible for the variations in different states are statistically constrained in BL Lac for the first time.Moreover, we constrained the underlying particle distribution responsible for the broadband emission using the χ 2 test.The multiwavelength data used in this work is acquired from Fermi-LAT, Swift-XRT and Swift-UVOT.The framework of this paper is as following:-the details of the multiwavelength data and the data analysis procedure are given in section §2.We presents the results of multiwavelength temporal and spectral analysis in section §3.We summarise and discuss the results in section §4.A cosmology with ΩM = 0.3, ΩΛ = 0.7 and H0 = 71kms −1 Mpc −1 is used in this work.

DATA ANALYSIS
In order to examine the temporal and spectral properties of BL Lac, we first obtained mutiwavelenth data including Optical/UV, X-ray and γ-ray.The data is acquired from Swift-UVOT, Swift-XRT and Fermi-LAT.Details of the observations and the data analysis procedure followed in our works are given below.

Fermi -LAT
Fermi-LAT is a high energy space based telescope with wide field of view ∼ 2.3 Sr.It is one of the two instrument onboard Fermi Gamma-ray Space Telescope (formerly called GLAST), which was launched by NASA in 2008.Fermi-LAT principally operates in scanning mode, it surveys the entire sky in the energy range ∼ 20 MeV-500 GeV every three hours (Atwood et al. 2009).In this work, γ-ray data of BL Lac is retrieved from Fermi-LAT for the time period MJD 59000-59943.The data is converted to science products using the Fermitools (formally Science Tools) with version 2.2.0, hosted on an Anaconda Cloud channel that is maintained by the Fermi Science Support Center (FSSC).We followed the standard analysis procedure described in the Fermi-LAT documentation 1 for the data reduction.The P8R3 events were extracted from the 15 degree ROI centred at the source location.The events having high probability of being photons are included in analysis by using the SOURCE class events as "evclass=128, evtype=3".Since Earth limb is strong source of background γ-rays, we avoided the contamination of γ-rays from the bright Earth limb by using a zenith angle cut of 90 degree to the data.This zenith cut is recommended by the LAT instrument team above 100 MeV.The good time intervals (gti) in which the satellite was working in standard data taking mode are selected by using a filter expression (DAT A−QU AL > 0)&&(LAT−CON F IG == 1).We modelled the Galactic diffuse emission component and the isotropic emission components with gll−iem−v07.fits and iso−P 8R3−CLEAN−V 3−v1.txt,respectively, and the post launch instrument response function used is P 8R3−SOU RCE−V 3.All the sources in the 4FGL catalog within (15+10) degree ROI centred at the source location were included in the XML model file.We initially carried the likelihood analysis for the entire period of data, while keeping the spectral parameters of sources free within 15 degree ROI, and freezing the parameters of sources lying outside the ROI.We allowed to vary the photon index and normalisation of the Galactic diffuse component, and normalisation of isotropic component during the spectral fitting.In the output model file, we freezed the spectral parameters for the background sources with test statistic, TS < 25, and the output model is finally used for the generation of light curve and the spectrum of the source.We considered the detection of the source only if TS > 9 (∼ 3 sigma detection; Mattox et al. 1996).

Swift-XRT
In this work, the X-ray data is acquired from the Swift-XRT telescope onboard the Neil Gehrels Swift Observatory (Gehrels et al. 2004).During the period MJD 59000-59943, a total of 100 Swift observations of BL Lac source were available.We obtained the Swift-XRT light curve such that each observation ID corresponds to one point in the X-ray light curve.The X-ray data is processed with the HEASoft package (v6.30) and CALDB (v20220803).We used the XRT-DAS software and the Standard xrtpipeline (v0.13.7) to create the cleaned event files.We run the xrtpipeline by giving 1 http://fermi.gsfc.nasa.gov/ssc/data/analysis/ the appropriate inputs to the task following the Swift X-ray data analysis thread page.The source and background regions are chosen by using the xrtgrblc task (Stroh & Falcone 2013).The task selects the source and background regions based on count rate such that a circular source region is selected for count rate ≤ 0.5 counts s −1 and annular source region is selected for count rate > 0.5 counts s −1 .While background region is chosen as annular region in all cases.Further, for the spectral analysis we have used an automated online X-ray analysis tool available at the UK Swift Science Data Centre (Evans et al. 2009).This online tool gives the source, background, and ancillary response files necessary for the spectral analysis.We used the GRPPHA task to rebin the source spectrum such that resultant spectrum has 20 counts per bin.The grouping is needed in order to evaluate the model based on the C-statistic (cstat).The X-ray spectral analysis in the energy range 0.3-10 keV was performed using the XSPEC package (Arnaud 1996) available with the HEASoft.Since the value of Galactic neutral hydrogen column density (nH) had been reported to vary between (1.7 − 2.8) × 10 21 cm −2 (Weaver et al. 2020), therefore we fitted 0.3-10 keV spectra with an absorbed power-law (PL) by choosing the nH value between (1.7 − 2.8) × 10 21 cm −2 , and keeping the normalisation and the spectral index as free parameters.

Swift-UVOT
In addition to the X-ray data, Swift also provides an Optical/UV data via Swift-UVOT telescope (Roming et al. 2005).It observes in optical and UV with the filters v, b, u and w1, m2, and w2 (Poole et al. 2008;Breeveld et al. 2010).We obtained the data of BL Lac from HEASARC Archive and reduced it to the scientific product by using the HEA-Soft package (v6.26.1).The uvotsource task included in the HEASoft package (v6.26.1) were used to process the images.Multiple images in the filters were added by using the uvotimsum tool.The source counts were extracted by using a circular source region of radius 5 ′′ centred at the source location and background region from a nearby source free circular region of radius 10 ′′ .The observed fluxes were dereddened for Galactic extinction using E(B − V) = 0.2821 and RV = AV /E(B − V ) = 3.1 following Fitzpatrick & Massa (2007).We further corrected the flux densities for the contribution of the host galaxy by using the method outlined in Raiteri et al. (2013).This prescribed method recommends flux density values of 2.89, 1.30, 0.36, 0.026, 0.020 and 0.017 mJy for the host galaxy in the v, b, u, uvw1, uvm2 and uvw2 bands, respectively.Moreover, the method suggests to subtract ∼ 50% of the total host flux density to the observed flux density.

RESULTS
BL Lac has shown number of flaring events over the last years across the energy bands, these events have been reported in number of astronomical telegrams (some recent telegrams are Buson & de Menezes 2021;Blanch et al. 2021;La Mura 2022;Prince 2022).Recently, during the renewed activity from BL Lac, the γ-ray flux ∼ 5 × 10 −6 ph cm −2 s −1 had been reported from the source (La Mura 2022).BL Lac has also shown several bright VHE gamma-ray flares (Blanch 2020;Blanch et al. 2021;Cortina & CTA LST Collaboratoin 2021).Motivated by the renewed flaring activity, significant variability and the availability of multiwavelength observations of BL Lac, we carried a detailed multiwavelength study of BL Lac by using the Fermi-LAT and Swift-XRT/UVOT observations.The aim is to understand the temporal and spectral characteristics of source.

Temporal Study
We first analysed the Fermi-LAT data of BL Lac acquired during the period MJD 59000-59943.We obtained a daily binned γ-ray light curve by integrating the photons over the energy range 100 MeV-500 GeV.Specifically, the daily binned differential flux of BL Lac was modelled by power-law (PL) model where N0 is the prefactor, Γ is the spectral index and E0 is the scale energy.During analysis, N0 and Γ were kept free in the fitting process, while E0 is chosen as 870 MeV.We consider the source detection if TS obtained from the maximumlikelihood analysis exceeds 9, which corresponds to approximately a 3 sigma detection level (Nolan et al. 2012).In all the time bins, we have ensured the convergence of the likelihood fit.The resultant daily binned γ-ray light curve acquired is shown in Figure 1, all points displayed on the light curve exhibit a 3 sigma detection.As shown in Figure 1, the source started a major γ-ray activity around MJD 59800 and continued high activity for more than two and half years.During the active period, a maximum daily averaged γ-ray flux of (1.74±0.09)×10−5 ph cm −2 s −1 is observed on MJD 59868.5, with a corresponding spectral index of 1.97 ± 0.04 and a TS value of 8008.This is the highest ever γ-ray flux detected from BL Lac and is factor of ∼ 300 larger than the base flux F0 ∼ 4.96×10 −8 ph cm −2 s −1 , which is the minimum flux detected during this period.The daily binned γ-ray light curve shows that the two maximum flux values are outliers.To check the accuracy of these flux values, we perform further analysis by obtaining the 6 hour binned γ-ray light curves around these flux values, which are depicted by the shaded strips in the Figure 1.We noted that the 6 hour binned light curves (shown in the inset plots) exhibited a similar pattern, with peak flux values of (2.83 ± 0.12) × 10 −5 ph cm −2 s −1 and (2.92 ± 0.25) × 10 −5 ph cm −2 s −1 observed on MJD 59331.12 and MJD 59868.37,respectively.Moreover, the γ-ray light curve shows large flux variations with many low and high flaring components.We calculated the rise and decay times of these components by using the sum of exponential (SOE) function where F b is the base line flux and Here Fp,i is peak flare amplitude at time tp,i; τr,i and τ d,i are rise and decay times of the respective flare component.
The fitted SOE profile along with the daily binned γ-ray light curve points are shown of Figure 2. We used more than 25 exponentials in the SOE function, which resulted in χ 2 /dof of 8754/816.We noted that adding more exponential in the SOE does not improve the fit statistics significantly.The best fit parameters of the components for which the peak flux is greater than 5 × 10 −6 ph cm −2 s −1 are given in Table 1.Using rise/decay time scales, we determined the profile shapes of these components by calculating the parameter, ζ = τ d −τr τ d +τr , such that the component is symmetric if |ζ| < 0.3, moderately asymmetric if 0.3 < |ζ| < 0.7, and asymmetric if 0.7 < |ζ| < 1.We noted that five components are moderately asymmetric, one component is asymmetric and six components are symmetric.
The daily binned γ-ray light curve is obtained by fitting the integrated spectrum with PL model.In Figure 3, we plotted the index as function of daily binned γ-ray flux, the blue solid circles represent the individual flux-index points.We also sorted the flux and index values in the increasing order of flux and then average them over the bin of 13 points.In the Figure 3, red diamond points represent the weighted average index values as function of weighted average flux values.The Figure 3 suggest that source exhibits a mild harder when brighter trend, a usual trend observed in blazars (e.g Britto et al. 2016;Shah et al. 2019).We used the Spearman rank correlation method to calculate the correlation coefficient and null-hypothesis probability between index and flux values.The returned values of the correlation coefficient −0.37 and null-hypothesis probability 4.52 × 10 −31 further confirms a mild anti-correlation between the two quantities.
In order to obtain an estimate of shortest flux doubling time scale, we scanned the daily binned γ-ray light curve with the equation where F(t0) and F(t) are the values of flux at time t0 and t respectively and τ is characteristic doubling time scale.On applying the condition that the significance of difference in flux at t and t0 is ≥ 3 sigma (Foschini et al. 2011), the daily binned γ-ray light curve resulted in shortest flux doubling timescale, tvar = 0.40 d during the time interval 59022-59023.

Fourier spectral analysis
Fermi-LAT provides longterm data of BL Lac object covering several years of observation.As shown in Figure 1, the γ-ray light curve of BL Lac shows unpredictable aperiodic variability with flux values changing by several order in magnitude.Such random variability is often result of a stochastic processes, instead of deterministic process.We used PSD to characterise the variability.PSD is one of the common tools used to examine the light curve variability, it gives an estimate of variability power as a function of temporal frequency.In order to obtain an insight into the physical process which causes large variability in the γ-ray light curve, we obtained PSD of source by splitting daily binned γ-ray light curve into equal segments, and calculated periodogram in each segment.We finally averaged them into final periodogram.The averaged power spectrum is normalised by using the fractional r.m.s normalisation (see   4 and best fit values of N and αp are obtained as 14.94 ± 2.76 and 1.18 ± 0.06, respectively.The index ∼ 1 of power spectra suggest for flicker-noise type process.Similar results have also been reported by (Sobolewska et al. 2014;Bhatta & Dhital 2020) in the γ-ray light curves of few other blazars.The flicker noise has the property of maintaining shape over several orders of frequencies up to

Flux distribution
Flux distribution study of the light curves of astrophysical systems is another tool used to probe the nature of underlying physical processes responsible for the variability.For example observation of normal flux distribution would indicate an additive processes, while a lognormal distribution would imply a multiplicative processes.Mostly the distribution of observed flux of compact black hole systems follow a lognormal distribution.The linear r.m.s-flux relation suggests that the flux distribution of the γ-ray light curve is lognormal.To investigate this, we characterise the γ-ray flux distribution of BL Lac by performing the skewness, AD and histogram fitting test.The skewness and AD test rejects the normality of flux distribution.Skewness value of the flux distribution is obtained as 2.55±0.16,and the AD test yielded a statistic value of 44.13, which is much larger than the critical value (CV) of 0.78 defined at a 5% significance level.However, the AD test suggests the lognormality of the flux distribution.
The statistic value of 0.75 for the log of the flux distribution is smaller than CV of 0.78 at 5% significance level, which suggest the null hypothesis that the flux distribution is lognormal can not be rejected.We further checked the PDF of flux distribution by constructing the normalised histogram of logarithm of flux.The histograms were constructed such that each bin contains equal number of points while varying the bin width.The normalised histogram points are plotted in Figure 6.The resulting histogram in log-scale is fitted by and where µ l and σ l are the mean and standard deviation of logarithm flux distribution, µg and σg are the mean and standard deviation of flux distribution.Equation 4 results in lognormal fit, while Equation 5 results in normal fit.The normalised histogram and best lognormal/normal fit to the histogram are shown in the left and right panel of Figure 6.The reduced-χ 2 values obtained from fitting the flux distribution with a normal and a lognormal PDFs were 3.33 and 1.01, respectively.This result further suggests that the flux distribution is more accurately described by a lognormal distribution.Observation of lognormal distribution in the γ-ray light curve suggest the underlying physical processes should be multiplicative in nature.It also suggest for two flux states in the source.

Multiwavelength Light Curve
In order to unveil the behavior of BL Lac at optical, UV and X-ray band, Swift-XRT/UVOT had carried a total of 100 observations of the source during the time period MJD 59000-59943.The X-ray and optical/UV data corresponding to these observation were analysed by following the analysis procedure described in section 2. The acquired X-ray and optical/UV light curves of BL Lac are shown in the second and bottom panels, respectively of the multiwavelength light curve plot (see Figure 7).Each point in the X-ray and optical/UV light curves corresponds to individual observations.The top panel shows the daily binned γ-ray light curve with flux points obtained by integrating over the energy range 0.1-500 GeV.The multiplot shows a simultaneous flux variations in different energy bands.We quantified variability in each energy band by calculating fractional variability amplitude using the equation (Vaughan et al. 2003) where S 2 is variance and F is mean of flux points in the light curve, and σ 2 err is mean of square of the measurement errors.The uncertainty on Fvar is calculated using the equation (Vaughan et al. 2003) where N is the number of points in the light curve.The values of Fvar acquired in the considered energy bands are plotted against the energy in Figure 8.The plot shows that the variability amplitude increases with energy from optical to γ-ray band, the variability amplitude is highest in the γ-ray band.This trend is can be manifestation of cooling of relativistic electron such that higher energy electrons responsible for the γ-ray emission cool faster than the low energy electrons, which results in faster variability in the high energy emission.

Broadband Spectral Analysis
As shown in the multi-wavelength plot (Figure 7), the flux variation in different energy bands are correlated.We employed the Z-transformed discrete correlation function (ZDCF, Alexander 1997)  (the uncertainties are at 1 sigma confidence level) between the γ-ray light curve and the X-ray, B, U, V, W1, M2, W2 light curves, respectively suggest a positive correlation between the emissions in different bands.Specifically, the acquired values indicate that the γ-ray emission does not exhibits a significant time lag with the X-ray and optical/UV emission.In order to identify the underlying particle energy distribution and physical parameters (like magnetic field, size of emission region, bulk Lorentz factor etc) responsible for the simultaneous flux variations, we examined the broadband spectral characteristics of BL Lac by choosing different time intervals from the multiwavelength light curve.Due to significant variability exhibited by BL Lac source on short timescales, we utilized Bayesian analysis to divide the γ-ray light curve into segments, with each segment assumed to exhibit a steady behavior in terms of the underlying physical parameters.This allows us to estimate the most probable values of the param- eters for each segment.We then selected the Bayesian segments for the broadband spectral analysis based on source activity and simultaneous observations available in γ-ray, Xray and UV/optical bands.We have shown these time segments by the vertical stripes in the multiplot (Figure 7) and identified them as S-1 MJD 59079-59084, S-2 MJD 59237-592425, S-3 MJD 59428-59434, S-4 MJD 59894-59897 flux states, Quiescent States: QS-1 MJD 59443-59457 and QS-2 MJD 59760-59774.By comparing the results from different segments, one can identify changes in the physical parameters that govern the flux variations.
We modeled the γ-ray spectrum in the considered flux states with LP model, defined in Equation 1 and PL model, dN dE = N0(E/E0) −Γ , where N0 is the prefactor, Γ denotes the spectral index and E0 represents the scale energy, which is fixed at 856 MeV.The fitting parameters obtained are summarised in Table 2.We determined the statistical significance of the curvature in the γ-ray spectrum by using the re-  As shown in Table 2, significant curvature (TScurve > 16) is observed in S-2, S-3 and QS-1 states.The γ-ray spectral points for the broadband SED modelling are obtained by dividing the total energy (1-500 GeV) into 8 energy bins equally spaced in log scale.For the S-2, S-3, and QS-1 states, the source spectrum is fitted using a LP model, while for the S-1, S-5, and QS-2 states, a PL model was employed.During the spectral fit in each energy bin, the spectral parameters of BL Lac were kept free, while the parameters of other sources in the ROI were frozen to the best fit values acquired in the energy range 0.1-500 GeV.The X-ray spectrum in each flux state is obtained by using an online automated products generator tool (Evans et al. 2009).We binned the acquired X-ray spectrum by using the GRPPHA task, such that their are 20 counts per bin.In case of Swift-UVOT, the images of the observation IDs falling in particular flux state are combined using the UVOTIMSUM task, and finally flux values are obtained from the combined image.The broadband SED points obtained in the S-1, S-2, S-3, S-4, QS-1 and QS-2 are shown in Figures 9, 10 and 11.
We consider one-zone leptonic model in order to model broadband SED in the selected flux states.In this model, we assume emission arises from a spherical blob of radius R filled with relativistic electron distribution, n(γ).The blob moves down the jet with bulk Lorentz factor Γ at a small angle θ with respect to line of sight of observer.Relativistic motion at small angle with respect to observer amplifies the blazar emission and this amplification is determined by a beaming factor δ = 1/Γ(1 − β cos θ), β is the velocity (in units of c) of blob.Further, we assume that the variability is governed by light crossing time scales so that size of emission is obtained through expression R ∼ δtvar/(1 + z).The relativistic electrons in presence of magnetic field, B and target photon field emit radiations through synchrotron and IC processes.In our model, we assume seed photons for IC process are synchrotron photons from jet itself, so that emission is through SSC process.We expressed the electron Lorentz factor γ in terms of new variable ξ such that ξ = γ √ C, where C = 1.36×10 −11 δB/(1+z).Following Begelman et al. (1984), synchrotron flux at energy ϵ can be obtained using the equation where, dL is luminosity distance, V is volume of emission region, A = C , ξmin and ξmax correspond to the minimum and maximum energy of electron, and f(x) is the synchrotron emmisivity function (Rybicki & Lightman 1986).The SSC flux received by the observer at energy ϵ can be obtained using the equation where, ϵi is incident photon energy, Similarly the EC flux recieved by the observer can be obtained using the equation where ph is target photon energy density and here y = 1 − √ Cϵs ξmec 2 .We solved Equations 8, 9 and 10 numerically and the resultant numerical code is incorporated as a local convolution model in XSPEC in order to perform a statistical fitting of broadband SEDs.The convolution code allows us to model the broadband spectrum for any particle energy distribution, n(ξ).In a single XPEC spectral fitting iteration, the synchrotron and IC ( (SSC/EC) process are solved simultaneously and same set of parameters (e.g B, Γ, R etc) are utilized in these process in a single itteration.In our convolution code, XSPEC 'energy' variable is interpreted as ξ = γC.Three cases of particle distribution viz., broken power law (BPL), LP, and the physical model namely energy dependent acceleration (EDA; Hota et al. 2021) were considered in our analysis.We first examined the underlying particle distribution by modelling broadband SED of S-1 flux state.Our results show that the convolved BPL electron distribution undergoing synchrotron and SSC loses provides better fit (χ 2 /dof ∼ 4.62/14) to the broadband SED than the LP and EDA model with χ 2 /dof of 30.45/14 and 32.83/14 respectively.This suggest that the underlying particle distribution responsible for broadband emission in BL Lac is preferably a BPL model.Hence, we used BPL form of electron distribution to fit the broad band SED in the considered flux states.A systematic of 10% is added evenly over the entire data to account for the additional uncertainties in the model.Using the convolved SED model involving synchrotron and SSC processes with BPL electron distribution, the observed broadband spectrum is determined mainly by 10 parameters viz.ξ b , ξmin, ξmax, p, q, Γ, B, R, θ and norm N. The code also allows us to fit the SED with jet power (Pjet) as one of the parameter, however in this case, N must be a fix parameter.We carried fitting with p, q, Γ, and B as free parameters, while other parameter were freezed to typical values required by the observed broadband spectrum.The reason for freezing the parameters is limited information available in the Optical/UV, X-ray and γ-ray bands.Moreover, we used Tbabs model to account for the absorption in the X-ray spectrum.We noted that synchrotron and SSC emission provided a reasonable fit to all the flux states with χ 2 /dof values as 4.62/14, 14.96/14, 12.05/15, 8.16/14, 12.06/14, and 6.18/14 for S-1, S-2, S-3, S-4, QS-1 and QS-2 states, respectively.However, the equipartition values obtained as 2252, 216, 16100, 9060, 2096, 3864 in the S-1, S-2, Table 2.The parameters obtained by fitting the integrated γ-ray spectrum of S-1, S-2, S-3, S-4, QS-1 and QS-2 states of BL Lac with the PL and LP model.Col. 1: flux state; 2: time period of flux state; 3: fitted model; 4: integrated flux in units of 10 −6 ph cm −2 s −1 ; 5: PL index or index defined at reference energy; 6: curvature parameter; 7: test statistics; 8: -log(likelihood); 9: significance of curvature.S-3, S-4, QS-1 and QS-2 flux states are much larger than unity.The results indicate additional process responsible for the high energy emission.The BL Lac source is wellknown for exhibiting emission line features originating from the BLR (Vermeulen et al. 1995;Corbett et al. 1996).These results indicate that the BLR region may play significant role for the high energy emission.Therefore, we investigated the broad band emission from BL Lac by considering the EC emission alongside synchrotron and SSC emissions.We assumed the seed photons for EC scattering are BLR photons.Interestingly, the emission lines detected from the source mainly consists of Hα emission lines (Vermeulen et al. 1995;Corbett et al. 1996), hence for numerical stability, we approximate the BLR emission as a blackbody with a temperature 42000 K (equivalent to the temperature corresponding to Lyman-alpha line emission at 2.5 × 10 15 Hz).Using this model, the observed broad-band spectrum can now be reproduced by introducing two additional parameters viz.

State
target photon temperature T and fraction of the photons which undergo Compton process f.The limited observation information and large number of model parameters introduces a degeneracy in model parameters.The optimum set of parameters would be those for which Pjet is minimum.To check for minimum Pjet, we varied ξmin, which corresponds to variation in γmin through the relation γmin = ξmin/ √ C. The Pjet and the χ 2 /dof obtained for different values of γmin are given in the Table 3.In all the flux states, Pjet shows decreasing trend as γmin increases, and the reduced-χ 2 obtained is reasonable good up to ξmin ∼ 10 −3 .Therefore, keeping the minimum jet power in mind, the value of ξmin is constrained between 10 −3 − 10 −4 in all the flux states in the final SED fit.The resultant best-fit model SED along with observed points are shown in Figures 9, 10, and 11 and corresponding best fit parameters are given in Table 4.

SUMMARY AND DISCUSSION
The Fermi light curve of BL Lac during period MJD 59000-59943 shows the source was in active state for long time.
During this period, γ-ray light curve revealed presence of multiple flaring components, a maximum daily averaged γray flux of (1.74 ± 0.09) × 10 −5 ph cm −2 s −1 is observed on MJD 59868.5 .This is the highest one day binned γ-ray flux detected from source.A shortest flux doubling timescale of tvar = 0.40 d is observed during the time interval 59022-59023.The γ-ray light curve require a series of exponentials to reproduce the profile shape.We noted that among dominant components, five components are moderately asymmetric, one component is asymmetric and six components are symmetric.The symmetry in the flare profile can be due to the light travel time effects, while the asymmetry in the flare profile could be attributed to strengthening and weakening of acceleration process.A slow rise in the asymmetric flare possibly indicates acceleration of particles to higher energy, while fast decay may be associated with the rapid energy loss of high energy particles.Moreover, a usual harder when brighter trend is observed in the γ-ray light curve.In case of BL Lac, the γ-ray spectrum lies near the peak of IC component, therefore the spectral hardening during the flaring indicates the shift in Compton peak towards high energy.
The spectral hardening and shift in SED peak energy of IC component during high flux states had been observed in 3C 279 by Shah et al. (2019).
We calculated PSD of γ-ray light curve in order to obtain an insight into the physical processes causing large variability in BL Lac object.We found that the PSD is a power law with an index ∼ 1, which suggest for a flicker noise type process.Similar results have also been reported by Sobolewska et al. (2014); Bhatta & Dhital (2020).In their work, PSD of γ-ray light curves of sample of blazars were found to be consistent with a power law where majority of the indexes were near to 1.0.The flicker noise is halfway between random walk (index = 2) and white noise (index = 0), it has the property of maintaining shape over several orders of frequencies up to arbitrarily low values.Therefore, observation of such feature in the light curve implies long memory process is at work.In other words, it implies that the shorter and longer timescale variations are coupled together or equivalently the underlying processes should be multiplicative.In blazars, the jet emission can possess memory of the events occurring at the accretion disk, especially the disk modulations and thus indicates for disk-jet connection.
We also characterised the variability of γ-ray light curve of BL Lac object by obtaining a correlation between the flux and r.m.s.The r.m.s-flux plot shows a linear trend, which has been noted in several other astrophysical systems like black hole binaries and Active Galaxies (Uttley & McHardy 2001;Vaughan et al. 2003;Uttley et al. 2005).The linear r.m.s-flux relation implies that short and long variability time-scales are coupled together in multiplicative way (Uttley et al. 2005).It also rules out shot-noise models where the different time-scales of variability are combined additively.In addition to r.m.s-flux relation, we checked the flux distribution of the γ-ray light curve with the skewness, AD test and histogram fitting.All these tests reject the normality of the flux distribution, instead suggest that the flux distribution is lognormal.Observation of lognormal distribution imply that the underlying emission process responsible for the variability is multiplicative rather than additive.These features are mostly believed to be result of perturbations in the accretion disk (Uttley et al. 2005;McHardy 2008;Lyubarskii Table 4. Best-fit parameters obtained by fitting the local convolution SED model involving SSC nad EC processes to S-1, S-2, S-3, S-4, QS-1 and QS-2 states.Row:-1: particle index before the break energy, 2: particle index after the break energy, 3: bulk Lorentz factor of the emission region, 4: magnetic field in units of 10 −3 G, 5: ξ break parameter, which represent break energy γ break through the relation γ break = ξ break / √ C, where C = 1.36 × 10 −11 δB/(1 + z), 6: equipartition parameter value 7: logarithmic jet power in units of erg s −1 , 8: χ 2 /degrees of freedom, and 9: Galactic neutral hydrogen column density (n H ) in units of cm −2 .The subscript and superscript values on parameters are lower and upper values of model parameters respectively obtained through spectral fitting.−− implies that the upper or lower bound value on the parameter is not constrained.For each of the flux states, size of emission region and viewing angle are chosen as 10 17 cm and 0.1 degree, respectively, ξ min values are chosen within the range of 10 −4 − 10 −3 , and ξmax is selected between 1 − 4. In the fluctuating accretion model, the perturbation in the mass accretion rate can propagate inwards and accure in multiplicative way in the inner regions of the disc, thereby producing the multiplicative emission.In case of blazar, the emission mainly comes from the jet, hence possible realisation for the observation of lognormal behavior could be that the disk fluctuations are imprinted on the jet emission (Giebels & Degrange 2009;McHardy 2010;Shah et al. 2018).
On the other hand, the minute time scale variation observed in the γ-ray light curves (Gaidos et al. 1996;Aharonian et al. 2007) implies that the jet emission should be independent of accretion disc fluctuations (Narayan & Piran 2012).In such cases, the lognormal distribution in flux can be explained by the linear Gaussian perturbations in particle acceleration time-scales (Sinha et al. 2018).Moreover, the fluctuations in the escape time scales of the electrons would produce flux distribution shapes other than Gaussian or lognormal.In addition, Biteau & Giebels (2012) have shown that the additive shot noise-model can also produce a lognormal flux distribution, such as the Doppler boosting of emission from a large number of randomly oriented mini-jets results in a flux distribution with features similar to that of a lognormal distribution.
We examined the broadband spectral characteristics of BL Lac by choosing the time intervals for which the simultaneous observations are available in γ-ray, X-ray and UV/optical bands.The convolved one zone leptonic model suggests that the underlying particle energy distribution responsible for the broadband emission is more likely to be a BPL distribution.Therefore, the broad band SEDs in different flux states are statistically modelled by using the BPL electron energy distribution which undergo synchrotron and SSC loss.The statistical fit is carried by keeping p, q, Γ, and B as free parameters, while other parameter are kept fixed to the typical values required by the particular flux state.
We showed that the optimal Pjet in different flux states (see Table 3) are obtained by choosing the ξmin within the range 10 −4 to 10 −3 .Under the conditions of equipartition and minimum jet power, the best fit parameters (shown in Table 4) implies that the increase in flux from low to high flux state is associated with the increase in Γ. Additionaly, the particle spectral indices becomes harder in the high flux state.
The SED of BL Lac has been the subject of multiple modeling attempts in the past.During the era of the EGRET satellite, it had been observed that SED modeling of the high flux states above 100 MeV requires the inclusion of external seed photons for Compton scattering.For example, Böttcher & Bloom (2000) showed that SSC and EC emission is required to yield an acceptable fit to the broadband spectrum.They showed that in BL Lac, unlike other BL Lac objects, the broad emission line region plays an important role for the high energy emission.Using the broadband emission model of Ravasio et al. (2002), Albert et al. (2007) showed that one zone SSC model can explain the broad band emission upto VHE energies during relatively weak γ−ray emission, while the strong γ−ray emission during the 1997 flare requires SSC as well as EC emission components for the broadband emission.Additionally, Abdo et al. (2011) showed that the SED may be described by a single zone or two zone SSC model, but a hybrid SSC plus EC model is preferred based on the observed variability.In this work, we constrained the underlying particle distribution using the χ 2 test.We noted that synchrotron and SSC emission provided a reasonable fit to all the flux states considered in our analysis.However, the equipartition values obtained are much larger than unity.Therefore we modelled the broad band emission from BL Lac by considering the EC emission alongside synchrotron and SSC emissions.These three emission processes ensure the equipartition between the magnetic energy density and the particle energy density.The best fit broken power law indices obtained show that index after break energy are steeper than what would be expected from the cooling break.The results thus rules out the broken power law spectrum as originating from a radiative cooling.The exact explanation for this steep spectrum remains unclear.It could potentially result from multiple acceleration processes or energy dependence of diffusion time scales.For example, a steep electron spectrum may possibly be interms of an arbitrary energy dependence of the diffusion coefficient.Zirakashvili & Aharonian (2007) studied the energy spectra of shock accelerated electrons and their associated radiation by considering an arbitrary energy-dependence of the diffusion coefficient.The authors show that in case of Bohm diffusion, the spectral cutoff takes on a sub-exponential form/steeper cutoff at high energies.Within this framework, the electron spectrum exhibits a significantly faster decay rate than the broken power law observed in our work.Alternatively, Sahayanathan (2008) proposed a two-zone model where broken power law injection into the cooling region undergoing synchrotron losses.introduces an additional break in the electron spectrum with indices p + 1 and p + 2 where p is the index of spectrum before the break energy γ b .Nonetheless, the presence of this additional break in the spectrum within our work can not be confirmed due to the energy gap between the observed optical/UV and X-ray spectra.
The MLC reveal that the BL Lac object exhibits a strong correlated variability at optical/UV, X-rays and γray bands.Analysis, particularly using the ZDCF, shows a positive correlation among the emissions in different bands without any significant time lag between them.This implies that single emission region and same electron population are responsible for the emission in different energy bands during a particular flux state.Moreover, the light curves show large amplitude of variability in X-ray flux compared to the variability in optical/UV and γ-ray band.A similar result was reported by Prince (2021).This trend may be attributed to the shape of the broad band SED.The broadband SED of BL Lac source during high flux states (see Figures 9 and 10) shows that the X-ray spectrum lie after the break energy.This implies that the X-ray emission are due to the high energy electrons, while the γ − ray and optical/UV emission are comparatively due to low energy electrons.Since the high energy electrons cool faster, theoretically one expects larger amplitude variability in flux at X-ray band compared to the γ-ray and optical/UV band.

Figure 1 .
Figure 1.One-day binned integrated Fermi-LAT γ-ray light curve of BL Lac obtained during the period MJD 59000-59943.The shaded strips represent the time periods for which the 6 hour binned γ-ray light curves are obtained.The inset plots display the 6 hour binned γ-ray light curves.

Figure 2 .
Figure 2. Daily binned γ-ray light curve of BL Lac fitted with the SOE function defined in Equation 2 .

Figure 3 .
Figure 3. Daily binned γ-ray index plotted as function of flux.Red diamond points represent the weighted average index values as function of weighted average flux values.

Figure 4 .
Figure 4. Power spectrum of one day binned γ-ray light curve.The fractional rms-squared normalization, also known as rms normalization is used in the plot.

Figure 5 .
Figure 5.The variation of average absolute r.m.s values as a function of one day mean γ-ray flux (in units of ph cm −2 s −1 ).Solid line represent the best linear fit.

Figure 6 .
Figure 6.Left panel: Normalised histogram of BL LAC fitted with the Gaussian PDF.Right panel: Normalized histogram of BL LAC fitted with the lognormal PDF.The flux are in units of ph cm −2 s −1 .

Figure 7 .
Figure 7. Multiwavelength light curves (MLC) of BL Lac obtained by using Fermi-LAT and Swift XRT and UVOT observations.The observations spanned a period from MJD 59000 to 59943.Top panel is daily binned γ-ray light curve, second, third and fourth panel are the X-ray, UV and optical light curves.The colored vertical stripes indicate the regions where broadband spectral modeling is performed.

Figure 8 .
Figure 8.The fractional variability amplitude (Fvar) of BL Lac plotted as function of Energy

Figure 9 .
Figure 9. Broadband SED of BL Lac obtained during the flux state S-1 (left panel)and S-2 (right panel).The flux points are represented by open diamond (Swift-UVOT), filled circles (Swift-XRT), and filled squares (Fermi-LAT).The solid red curve represents the combined best fit synchrotron, SSC and EC spectrum.

Figure 10 .Figure 11 .
Figure 10.SED of BL Lac obtained during the flux state S-3 (left panel) and S-4 (right panel).The labelling are same as that of Figure 9.

Table 1 .
Rise and fall time of dominant components of the light curve obtained by fitting the SOE (Equation2) .Col. 1: peak time (MJD); 2: peak flux (in units 10 −6 ph cm −2 s −1 ); 3 and 4: rise time and decay time of the components (in days); 5: asymmetry parameter.

Table 3 .
The table shows variation of the Jet power with the γ min obtained by using the local convolution SED model.Col: 1. Flux state, 2. ξ min parameter, 3. γ min , 4. Jet power and 5. χ 2 /dof.