Study of Intra-Day Flux Distributions of Blazars Using XMM-Newton Satellite

We present a study of the flux distribution of a sample of 15 Intermediate and Low-energy peaked blazars using XMM-Newton observations in a total of 57 epochs on short-term timescales. We characterise the X-ray variability of all of the light curves using excess fractional variability amplitude and found that only 24 light curves in 7 sources are significantly variable. In order to characterise the origin of X-ray variability in these blazars, we fit the flux distributions of all these light curves using Gaussian and lognormal distributions, as any non-Gaussian perturbation could indicate the imprints of fluctuations in the accretion disc, which could be Doppler boosted through the relativistic jets in blazars. However, intra-day variability, as seen in our observations, is difficult to reconcile using disc components as the emissions in such sources are mostly dominated by jets. We used Anderson-Darling (AD) and $\chi^{2}$ tests to fit the histograms. In 11 observations of 4 blazars, namely, ON 231, 3C 273, PKS 0235+164 and PKS 0521-365, both models equally fit the flux distributions. In the rest of the observations, we are unable to model them with any distribution. In two sources, namely, BL Lacertae and S4 0954+650, the lognormal distribution is preferred over the normal distribution, which could arise from non-Gaussian perturbations from relativistic jets or linear Gaussian perturbation in the particle time scale leading to such flux distributions.


Introduction
Blazars are well known for their collimated relativistic jets aligned to the line of sight of the observer [1]. Blazars constitute BL Lacertae objects (BL Lacs) and Flat-Spectrum Radio Quasars (FSRQs). The significant difference between these two classes is their optical emission which is strong for FSRQs but weak/absent for BL Lacs [2]. Emissions from blazars span across the whole electromagnetic spectrum (EM) from radio to high-energy gamma rays. The spectral energy distribution of blazars consists of two broad humps, where the low-energy hump peaks in the optical/X-ray band while the high-energy hump peaks at GeV/TeV energies. Recently, Fan et al. (2016) [3] (see also [4,5]), classified blazars based on the peak frequency of their low-energy component, namely, low-synchrotron frequency peaked blazars (LSP; ν s ≤ 10 14 Hz), intermediate-synchrotron frequency peaked blazars (ISP; 10 14 < ν s < 10 15.3 Hz) and highsynchrotron peaked blazars (HSP; ν s ≥ 10 15.3 Hz). The low-energy component of a SED is ascribed to synchrotron emission by relativistic electrons in the jet (e.g., [6]); whereas the high-energy component is likely due to Inverse Compton (IC) scattering of seed photons off the relativistic electrons. The seed photons are the synchrotron photons in the synchrotron self Compton (SSC) model, whereas the external Compton (EC) scenario foresees that seed photons may come from the accretion disc, broad line region (BLR) or hot corona surrounding XMM-Newton observations for this sample of blazars with an exposure time greater than 10 ks. We also included LBL blazar 3C 273 in our studies. Due to a large number of XMM-Newton observations for blazar 3C 273, we selected 3σ variable observations from Gowtami et al. (2022) [30] including the latest observation of 28 June 2022. XMM-Newton data are obtained from the HEASARC archive 1 . We have used European Photon Imaging Camera (EPIC) pn instrument data. EPIC-pn is most sensitive and less affected by the photon pile-up effects [31]. Science Analysis System (SAS-19.1.0) is used to analyse the EPIC-pn data. SAS task epproc is used to generate the calibrated and concatenated event lists. An X-ray energy range of 0.3-10 keV is considered since data below 0.3 keV are considerably dominated by detector noise components, and data above 10 keV are dominant with background flaring events and hot pixels. Background flaring in the X-ray energy range of 10-12 keV in the observations is checked and avoided from the event list files. A circular region of 40 arcsec radius and circular annulus of 50 (inner radius) to 60 (outer radius) arcsec centred on the source is used to extract the source and background light curve, respectively, for the X-ray energy range of 0.3-10 keV with a time binning of 100 s. Pile-ups in the observations are checked and are not found in any of the observations. SAS task epiclccorr is used to correct for vignetting, PSF variation, bad pixels and to get a background-subtracted light curve.

Analysis Techniques
Excess variance is the measure of intrinsic variability in the light curve and is determined by subtracting the variance due to measurement errors from the total variance of the observed light curve [32]. The fractional rms variability amplitude F var is the square root of the normalised excess variance [33,34]. Histograms are used to visualise the density or frequency distribution of any measured variable. The number of bins affects the shape of the histogram. Therefore, we need an optimal data-based binning for a histogram. The optimal binning algorithm, i.e., the Knuth method [35], is used to obtain the histogram of a count rate distribution. A histogram of a count rate distribution for individual light curves is obtained and fitted with normal and lognormal probability distribution functions. Normal and lognormal distribution functions are defined as follows: where µ and σ are the mean and standard deviation of the distribution, respectively, in units of counts/sec. The Anderson-Darling test (AD) [36,37] is used to quantify the lognormal and normal distribution function fitting to the histogram of the count rate distribution. When the p-value of the AD test is less than or equal to 0.01, we reject the hypothesis of normality/lognormality. The χ 2 goodness-of-fit test is also used to quantify the normal and lognormal distribution fitting for count rate histograms, and we choose the value of reduced χ 2 close to unity as the one representing the data better with null hypothesis probability value lying within the range of 0.75-0.05. If the probability value lies outside this range, then the null hypothesis is rejected. The distribution model is preferred only when both AD and χ 2 tests accept the null hypothesis, where the null hypothesis is that the data follow the specified distribution.

Results
We analysed all of the XMM-Newton satellite observations of our sample of ISP and LSP blazars for 57 epochs. We estimated the fractional root mean square (rms) variability amplitude F var of all of these observations and found significant (3 sigma level) variations in 7 sources during 24 observational epochs. In order to study the statistical properties of X-ray variability of these light curves, we considered only significantly variable observations in our studies. Our sample of blazars, classifications and F var are tabulated in Table 1. The light curves of blazars are represented in Figure 1. Histograms of count rates of these light curves are obtained and are fitted with lognormal and normal probability distribution functions (PDFs). Histograms and their corresponding distribution functions (i.e., normal and lognormal PDFs) are plotted in Figure 2. The statistical analysis of the AD test yields the AD statistic, p-value, and the χ 2 test yields a reduced χ 2 (χ 2 r ) value and a probability value, which are presented in Table 1. BL Lac is studied by Giebels and Degrange (2009) [19] using the RXTE-PCA observations, and it is the first blazar in which lognormal X-ray variability is detected, and a linear correlation between the excess rms and the average flux was also found. The flux distribution of blazar S5 0716+714 in X-ray with XMM-Newton observations was studied on a short-term timescale by Kushwaha and Pal (2020) [25] and Mohorian et al. (2022) [38]. 3C 273 is studied in X-ray using 16 years of Rossi X-ray Timing Explorer (RXTE) archival data by Khatoon et al. 2020 [39]. The flux distribution of FSRQ 3C 273 is lognormal, while its photon index distribution is Gaussian. This result is consistent with linear Gaussian perturbation in the particle acceleration timescale, which produces lognormal distribution in flux. On the decade-long timescale, BL Lac, 3C273, S5 0716+714 and PKS 0235+164 follow the lognormal flux distribution and display a strong linear rms-flux relationship in optical [40] and in weekly binned gamma-ray [41] observations. The monthly average γ-ray flux distribution favours the lognormal flux distribution, (i.e., [26]) and γ-ray (0.1-300 GeV) observations from Fermi/LAT telescope using maximum likelihood estimation (MLE) methods suggest that the γ-ray flux variability can be characterised by log-stable distributions [42].
The results for individual sources are provided below: BL Lac: This blazar is observed in 5 XMM-Newton observations and found to be variable in only one light curve of 16 May 2008. The AD test shows that the hypothesis of normality is rejected since the p-value for a normal distribution is less than 0.01, whereas the hypothesis of lognormality is accepted. The χ 2 test also gives a reduced χ 2 value of close to 1, and the probability value lies in the range of 0.75 to 0.05 for lognormal PDF fitting as compared to normal PDF fitting. Therefore, a lognormal distribution is a preferred model over a normal one.
ON 231: For observations taken on dates 26 June 2002, 14 June 2008 and 18 June 2008, the AD test rejects the hypothesis for normal as well as lognormal model. The χ 2 test also yields a reduced χ 2 value much greater than 1 for both PDF fittings. Therefore, in these cases, none of the distributions is the preferred model. For the observation from 12 June 2008, the AD test accepts the hypothesis for normal as well as lognormal distributions, and the χ 2 test also yields a reduced χ 2 value close to 1 with a probability value in the range of 0.75 to 0.05 for both models. Therefore, we cannot distinguish between these two models. For the observation from 16 June 2008, the AD test rejects the hypothesis of normality and accepts the lognormality hypothesis. However, it has a reduced χ 2 value not close to 1 with a probability value lying outside the prescribed range. Therefore, in this observation, we can not prefer any distribution.
S5 0716+714: For the observation from 24 September 2007, the AD test rejects the normal hypothesis as well as the lognormal distribution, and the reduced χ 2 value is also not close to 1 for both models. Therefore, we cannot prefer any distribution for this observation. For the observation from 4 April 2004, the AD test rejects the hypothesis of normality but accepts the lognormality hypothesis. However, the corresponding reduced χ 2 value for the lognormal distribution is not close to 1, with a probability value lying outside the specified range. Therefore, in this case, we can not prefer any model over another.
3C 273: For the observation from 12 December 2011, the AD test rejects the hypothesis of normality, whereas it accepts the lognormality hypothesis, but its corresponding reduced χ 2 value is not close to 1, with a probability value lying outside the specified range. Therefore, no distribution is the preferred model. For observations taken on the dates 13 June 2001, 12 January 2007, 9 December 2008 and 6 July 2020, the AD test accepts both the normality and lognormality hypotheses, but their corresponding reduced χ 2 values are not close to 1, with a probability value lying outside the specified range. Therefore, both distributions are not favourable for these observation cases. Similarly, for observations from 14 June 2000, 30 June  2004, 8 December 2007, 20 December 2009, 10 December 2010, 13 July 2015, 26 June 2017 and  28 June 2022, the AD test accepts both the normality and lognormality hypotheses and their corresponding reduced χ 2 value is also close to 1, with a probability value lying in the specified range. Therefore, both the lognormal and normal distributions are favourable models in these observations.   10 and ρ 25 are the linear correlation coefficients for 10 and 25 data points per bin, respectively; dash '-' represents distribution fits well; and 'None' represents that flux distributions equally fit with both distributions and we cannot prefer any model. S4 0954+650: the AD test rejects the normality hypothesis and accepts the lognormality hypothesis. The reduced χ 2 value for the lognormal distribution is close to unity, with a probability value lying in the specified range. Therefore, the lognormal distribution is the preferred model over the normal one.
PKS 0235+164: the AD test accepts both normality and lognormality hypotheses. Their corresponding reduced χ 2 values are close to unity for both models, with a probability value lying in the specified range. Therefore, we cannot prefer any model from this source.
PKS 0521-365: the AD test accepts both normality and lognormality hypotheses. The χ 2 test yields a reduced χ 2 value close to 1, with a probability value lying in the specified range for normal as well as lognormal models. Therefore, in this observation, we cannot prefer any distribution.
In order to quantify the degree of intrinsic variability, we estimated the rms-flux relationship for two different binnings of 15 and 25 data points per bin. We tabulated the results of the linear Pearson correlation coefficient and their corresponding p-values for all sources in Table  1.

Discussion and Conclusions
In this work, we studied the statistical properties of the X-ray variability of a sample of 15 blazars, consisting of ISPs and LSPs observed using the XMM-Newton satellite. We analysed a total of 57 observation epochs for our sample to search for significant variability (at the 3σ level) by estimating the fractional variability amplitude F var and found 24 light curves to be variable for a total of 7 sources. In order to study the cause of variability, we model histograms of the count rate of all of these variable light curves using normal and lognormal distributions. The lognormality of flux distributions in blazars, if found, is very important as they imply that the X-ray emissions originating from jets have imprints of the accretion disk or its environment, as found in the literature for many X-ray binaries and Seyfert galaxies, where disk contribution is more compared to jets [17,[43][44][45]. Such behaviour was inferred for the famous, nearby HSPs Mrk 421 and Mrk 501 at a Very-High-Energy (VHE) band by [46]. The Third LAT AGN Catalog (3LAC; [47]) also showed a lognormal distribution for blazars. Kushwaha et al. (2016) [24] studied an FSRQ PKS 1510-089 using multi-wavelength observations and found that the flux distribution of optical and γ rays follow two distinctive lognormal profiles, whereas the X-ray flux distribution follows a single lognormal distribution. Sinha et al. (2016) [21] also studied HSP blazar Mrk 421 and found that the radio to VHE frequencies follows a lognormal distribution. Shah et al. (2018) [26] also analysed γ-ray data of 50 blazars and found that their flux distribution follows a lognormal distribution. However, all of the above studies are performed on short-or long-term observations.
Minute or intra-day variability, as seen in most of the blazars, are difficult to explain using the disc component [48], and strongly favour an origin from the relativistic jets. In the statistical studies of IDV light curves of blazars, if lognormality is found to be a better model for a normal distribution, it could arise from a few additive processes under specific scenarios. Many models have been proposed to explain the IDV of blazars, which involves the shock-in-a-jet model [49,50], needle-in-a-jet [51], jets-in-a-jet model [28,52] and turbulent extreme multi-zone (TEMZ) model (e.g., [53]). Biteau and Giebels (2012) [27] studied the properties of blazar light curves using a mini jets-in-a-jet model and found that flux from a single region produces a power law histogram and a lognormal distribution in the case of many randomly oriented mini jets. It has been found that the histograms having skewed/lognormal profiles followed a linear rms-flux relationship. Sinha et al. (2018) [23] also provided an alternative explanation for the lognormal behaviour of blazar light curves and proposed that the non-Gaussian signatures in blazar variability could arise from linear fluctuations in the underlying particle acceleration and/or the diffusive escape rate of the emitting electrons. Such small Gaussian perturbations propagating in blazar jets could produce non-linear flux distributions and can explain the lognormal behaviour. Kushwaha and Pal (2020) [25] studied blazars during high variability phases on IDV timescales and found the flux distributions show a normal profile compared to lognormal ones.
In our studies, we found the lognormal distribution to be a preferred model over the normal one for sources BL Lacertae and S4 0954+650 in two observational epochs. For these observations, we found a linear rms-flux relationship with a Pearson correlation coefficient of 0.61 and a p-value of 0.2 for BL Lacertae and 0.63 with a p-value of 0.56 for S4 0954+650. In 11 other observations from 4 sources, namely, ON 231, 3C 273, PKS 0235+164 and PKS 0521-365, we found that the flux distributions equally fit with both distributions, and we can not prefer any model. In the rest of the observations, we are unable to model them with normal or lognormal distributions.
If the variability is triggered due to fluctuations in the accretion disk and the jet is powered by the central source channelling the energy into non-thermal particles, then the observed flux distributions may be an imprint of the accretion disk onto the jet (e.g., [19]). All of the above scenarios might be plausible on the long-term timescale variability of blazars. The lack of correlation between the mass of the Supermassive Black Hole (SMBH) and the variability timescale in blazars [24] suggests that the jet can modify the imprint of accretion disc fluctuation, at least on IDV timescales. Based on observations and simulations, blazar jets are found to be highly magnetised systems and, therefore, magnetic reconnection is a potential mechanism for X-ray variability. Turbulent magnetic reconnection in the coronal region around the accretion disk can explain the emission from compact accreting sources but fails for blazars where the jet is closely pointing to our line of sight. A magnetic reconnection-based model, mini jets-in-a-jet, can not be favoured for these blazars as we did not find a significant linear relationship between rms and flux. Therefore, the lognormality found in our observations could be attributed to either due to non-Gaussian fluctuations in blazar jets or through linear fluctuations of the underlying particle acceleration [39] or the diffusive escape rate of the emitting electrons [23]. They showed that the perturbation in the acceleration timescale produces a Gaussian distribution in the index that results in a lognormal distribution of the flux. However, the perturbation in the particle cooling rate produces neither a Gaussian nor a lognormal flux distribution, which is also found in many of our observations. Similarly, in the shock acceleration scenario, a small perturbation in the acceleration timescale can lead to the variability of the particle number density, which is a linear combination of Gaussian and lognormal processes. The relative weight of these processes determines the dominant shape of the flux distribution. If the variability arises due to variations in the number density of the accelerated particles, then the flux distribution could appear as both Gaussian and lognormal, which is also seen in our observations. However, we need more observations and a bigger sample to characterise the statistical properties of the flux distribution of blazars.