1. Introduction
Blazars are a category of radio-loud active galactic nuclei with relativistic jets facing the Earth. Blazars emit in the entire observable electromagnetic spectrum, displaying variability over flux and spectra with time. These are classified as BL Lacertae (BL Lac) and flat spectrum radio quasars (FSRQs), depending on the absence or presence of the emission lines. FSRQs show strong and broad emission lines with an equivalent width (EW) of greater than 5 Å and BL Lac objects have a weak or featureless spectrum with EW less than 5 Å. Blazars display variability in flux across different timescales ranging from minutes to hours and even years.
The multi-wavelength spectral energy distribution (SED) of blazars in
$\nu$
-f
$\nu$
representation depicts broadband double-humped structures (e.g. Ghisellini et al. Reference Ghisellini1997). The low-energy hump, spanning the radio to X-ray bands, is produced by synchrotron emission from relativistic, non-thermal charged particles. The high-energy peak, extending from X-ray to TeV energies, originates from either leptonic or hadronic processes. In leptonic models, this high-energy emission results from relativistic electrons upscattering photons via the inverse Compton (IC) mechanism. The target photons are either the synchrotron photons themselves (Synchrotron-Self Compton, or SSC) or photons from an external source like the broad-line region, dusty torus, or accretion disc (External Compton, or EC) (e.g. Kirk, Rieger, & Mastichiadis Reference Kirk, Rieger and Mastichiadis1998; Gaur et al. Reference Gaur, Gupta, Lachowicz and Wiita2010). The hadronic model attributes the emission to the proton synchrotron emission and synchrotron and Compton emission from the decay of charged pions or the decay of neutral pions (e.g. Mücke et al. Reference Mücke, Protheroe, Engel, Rachen and Stanev2003). Blazars with peak of the synchrotron hump lying below
$10^{14}$
Hz are called low-energy peaked blazars (LBLs) and those with peaks in between
$10^{14}$
and
$10^{15}$
Hz are called intermediate-energy peaked blazars (IBLs). High-energy peaked blazars (HBLs) are blazars with the synchrotron peak lying above
$10^{15}$
Hz. The SED of blazars in X-ray helps to understand both the flux and spectral variability in the synchrotron and inverse Compton emission. For LBLs, the X-ray emission consists of inverse Compton emission in the leptonic scenario. The IBLs portray both the synchrotron and inverse Compton emission in X-ray and the spectra can be used to estimate the energy at which the upturn from synchrotron to inverse Compton occurs. The SED of HBLs shows the synchrotron part of the emission. Among many Wierzcholska et al. (Reference Wierzcholska and Wagner2016) and Ferrero et al. (Reference Ferrero, Wagner, Emmanoulopoulos and Ostorero2006), worked on the X-ray spectra of blazars to study the upturn from the synchrotron emission to inverse Compton emission.
The models commonly used to describe the X-ray spectra of BL Lac objects include the single power-law (Worrall & Wilkes Reference Worrall and Wilkes1990), broken power-law (Sambruna, Maraschi, & Urry Reference Sambruna, Maraschi and Urry1996), and log-parabola (Inoue et al. Reference Inoue and Takahara1996) models. The log-parabola model effectively describes the X-ray emission from a large fraction of BL Lac objects (Giommi et al. Reference Giommi2021; Donato, Sambruna, & Gliozzi Reference Donato, Sambruna and Gliozzi2005; Tavecchio, Maraschi, & Ghisellini Reference Tavecchio, Maraschi and Ghisellini1998). The double power-law model provides a better representation of the X-ray emission from blazars than the single power-law when both synchrotron and inverse Compton components contribute to the X-ray emission. Both double power-law model and log-parabola model help to identify the energy at which the spectral upturn occurs (Ferrero et al. Reference Ferrero, Wagner, Emmanoulopoulos and Ostorero2006). As noted by Donato et al. (Reference Donato, Sambruna and Gliozzi2005), HBLs typically exhibit an X-ray spectral index between 3 and 1.8, while for other blazars it can be lower than 1.8. The hardness of the photons also varies with flux variability, with Brinkmann et al. (Reference Brinkmann2001) reporting that flux variability is higher at higher energies.
One of the most common blazar among many that has been studied for both synchrotron and inverse Compton emission in X-ray is BL Lacertae (also known as BL Lac but for clarity, we will refer to it by full name). However, BL Lacertae was referred as an LBL by Nilsson et al. (Reference Nilsson2018) when the logarithmic value of synchrotron peak frequency in Hertz was found 13.99. This object displays complex behaviour with changing spectral behaviour (Prince Reference Prince2021). BL Lacertae like other blazars exhibits timing and spectral variability and has been monitored using various multiwavelength campaigns.
1.1. Previous X-ray results of BL Lacertae
Before BL Lacertae flared in X-rays (in the 0.3–10.0 keV band examined in this work) during 2020–2021, the two most widely studied X-ray flares of BL Lacertae were those of 2001 and 2011–2012 (Prince Reference Prince2021). Ravasio et al. (Reference Ravasio, Tagliaferri, Ghisellini, Tavecchio, Böttcher and Sikora2003) examined X-ray observations of BL Lacertae obtained with BeppoSAX and RXTE-PCA in July and October 2000. The BeppoSAX data showed a low-flux state on July 26–27, followed by a high-flux state during October 31 to November 2, 2000. In July, when the source was in a steady state, it exhibited a hard X-ray spectrum. The flux in the 0.7–2.0 keV energy range varied on timescales of a few hours, while the flux in the 2.0–10.0 keV range remained nearly constant. In October 2000, a transition from a soft to hard spectrum was observed in the 10–20 keV energy range.
An extensive multiwavelength analysis of BL Lacertae was conducted by Böttcher et al. (Reference Böttcher2003) for observations spanning from mid-May 2000 to the end of that year. The group collected data in the radio, optical, X-ray, and very high energy (VHE) bands, presenting results on flux and spectral variability as well as broadband SEDs. This study showed that the source was quiescent before September 2000 and then entered a prolonged flaring state. During the high-flux state, the X-ray band exhibited increased variability on timescales of less than a few hours. Evidence was found for the contribution of synchrotron emission up to 10 keV during the flaring state. An examination of the hardness–intensity relation using BeppoSAX data revealed that in the soft X-ray band (0.5–2 keV), the spectrum softened as the flux increased, while in the medium X-ray band (2–4 keV), the usual ‘harder-when-brighter’ behaviour was confirmed. The study also reported spectral softening with increasing flux in RXTE-PCA observations. Additionally, they found that optical variability traced the soft X-ray variability, but with a time lag of 4–5 h.
Earlier XMM-Newton observations of BL Lacertae were analysed by Raiteri et al. (Reference Raiteri2009), who reported that the X-ray spectra exhibited contributions from both synchrotron and inverse Compton emission. Previous studies on BL Lacertae have shown moderate intraday variability in optical multiwavelength observations (e.g. Gaur et al. Reference Gaur, Gupta, Bachev, Strigachev, Semkov, Wiita, Gu and Ibryamov2017). Bhatta, Mohorian, & Bilinsky (Reference Bhatta, Mohorian and Bilinsky2018) analysed X-ray observations of BL Lacertae taken with NuSTAR in December 2012 and found no correlation between hardness and flux during the observation. Wierzcholska et al. (Reference Wierzcholska and Wagner2016) analysed Swift-XRT observations of BL Lacertae and found that both the power-law and log-parabola models provided good fits to the observations. Their results showed different values of the spectral index and curvature depending on whether Galactic absorption was allowed to vary or kept fixed.
The source entered a flaring state in 2020, and several studies have examined observations taken during this period across multiple wavelengths. Sahakyan & Giommi (Reference Sahakyan and Giommi2022) conducted an extensive analysis of BL Lacertae using observations from 2008 to 2020, studying the flares that occurred in 2012 and 2020 and explaining the complex emission behaviour of the source. (Prince Reference Prince2021) investigated the broadband spectral evolution of BL Lacertae during its 2020 flaring period, analysing the October 2020 flare across multiple wavebands. Their studies suggested that BL Lacertae exhibited HBL-like behaviour during this flaring state.
1.2. BL Lacertae enters in a flaring state
BL Lacertae flared in the optical, X-ray, and gamma-ray bands in 2020. Gamma-ray flaring was observed with Fermi-LAT, coinciding with the optical flare (Cheung Reference Cheung2020). MAGIC detected VHE flaring on September 19, 2020 (Blanch Reference Blanch2020), following the flaring in GeV energy ranges reported by Fermi-LAT, which had been detected on August 19, 2020. (D’Ammando Reference D’Ammando2020) reported that the source reached its historical flux maximum in a Swift-XRT observation on October 6, 2020, which was contemporaneous with gamma-ray flaring observed by Fermi-LAT. On October 11–12, 2020, the source was found to be flaring again with NuSTAR, reaching a flux as high as
$3.44 \times 10^{-10}$
erg cm
$^{-2}$
s
$^{-1}$
in the 0.3–10 keV band (D’Ammando Reference D’Ammando2022). The source reached another maximum on January 18, 2021 in the optical band, accompanied by gamma-ray flaring detected by Fermi-LAT (Cutini Reference Cutini2021). The flaring state continued throughout 2021, and until August 2021, BL Lacertae was reported to be flaring in the optical, X-ray, and gamma-ray bands. (Prince Reference Prince2021) reported continued X-ray flaring activity from BL Lacertae, observed with Swift-XRT from July 8 to July 11, 2021 (Prince, Khatoon, & Agarwal Reference Prince, Khatoon and Agarwal2021).
2. Data selection and reduction
The observations for the study of intraday variability and spectral analysis of BL Lacertae have been taken from the XMM-Newton satellite.
Summary of observations taken from EPIC-PN instrument of XMM-Newton. Photon counts is the average photon count of an observation. Here observation is abbreviated as obs.
2.1. XMM-Newton
XMM-Newton satellite was launched in 1999. The instruments onboard XMM-Newton are – three European Photon Imaging Cameras (EPIC), two Reflection Grating Spectrometers (RGS), and one Optical Monitor telescope.
EPIC is an X-ray instrument with three CCDs – one of pn type and other two of MOS type. Each of the CCD has an effective area of approximately 1 500 cm
$^2$
with a large field of view of 30’ diameter. EPIC can observe in the photon energy range of 0.15–15 keV.
The time resolution of the CCDs of EPIC-pn telescope in timing and burst mode is 30
$\mu$
s which is higher than the time resolution of EPIC MOS telescope. Only 40% of the flux is incident on the MOS cameras because MOS cameras are behind the grating of RGS. Thus we have used X-ray observations taken from EPIC-pn instrument which captures the total incident flux.
2.2. Data selection
The archival observations of BL Lacertae taken using the PN detector of EPIC-pn instrument onboard XMM-Newton have been used for the computation of variability. There are five observationsFootnote a of BL lacerate taken from July 2007 till July 2021 which have been downloaded from XMM-Newton Science Archive.Footnote b Here we used the PN observations because of the higher sensitivity and lower photon pile-up effects of the PN detector as compared to the MOS detector of EPIC. All five observations have an observing time of more than 10 ks which makes all of these suitable for time-dependent spectral variability studies. The five observations that have been selected are given in Table 1 along with time and duration of observation.
The observation taken in July 2021 coincides with the source being reported as flaring in X-ray by Prince et al. (Reference Prince, Khatoon and Agarwal2021). The other four EPIC-pn observations are from the non-flaring state of BL Lacertae.
2.3. Data reduction and analysis
The observations were reduced to get spectrum and light curves using the Scientific Analysis System (SAS) software (SAS 21.0.).Footnote c Observation Data Files (ODFs) downloaded from the XMM-Newton data archive consist of the instrument science data, housekeeping and auxiliary files. It also includes the Current Calibration Files (CCFs) which are needed for data calibration. Tasks cifbuild and odfingest are used to obtain the list of calibration files and a summary file of all files in ODF set. All the tasks needed to produce an event list from PN data are compiled in one epchain task. Task evselect was used to create lightcurve and spectrum from a selected region.
Pile-up correction is another important part of analysing the observation. The pile-up effect occurs when more than one photon falls on a pixel and is counted as one photon of added energy of the striking photons. To correct the pile-up effect, the model distribution of a single photon event is compared to a double photon event. The overlap of the two distributions signifies the absence of any significant pile-up. To correct this effect in observation, the pile-up affected region is excised and only an annular region is selected. The task epatplot has been used for this correction.Footnote d This effect was observed only in one observation, 0501660201, which was corrected by using the standard procedure of choosing an annular region of source while for all other observations, a circular region was chosen. Source region selected from observation 0501660201 was an annulur region lying between the concentric circles with radius of 230 and 524 arcsec centred at the source. The source region selected for observations 0501660301 and 0501660401 was a cicular region of radius 656 and 360 arcsec centred at the source. The background region was selected from a circular region on the same CCD outside the source region. For observations taken in timing mode, source and background spectra were extracted from 1D RAWX strips, not 2D regions. Dedicated Timing mode RMFs and ARFs were generated, correctly accounting for the geometry and fast readout. Timing mode inherently mitigates pile-up, and we confirmed it was negligible. Light curves in soft energy bands were carefully inspected, and no significant soft flares were detected requiring additional filtering beyond standard particle background screening. For the two Timing Mode observations (0504370401 and 0891800501), the source region was extracted from detector columns 37–39, centred on the source. The background region was extracted from columns 30–33, outside of the source region.
In the analysis we use the data from the limited energy range, from 0.3 to 10.0 keV, as there is known detector noise below 0.3 keV and the measurements above 10.0 keV are subject to a high proton background caused by the solar activity (Bulbul et al. Reference Bulbul2020; Gaur et al. Reference Gaur, Gupta, Lachowicz and Wiita2010). We monitored such effects by checking the derived light curves in this energy range.
The light curves and spectrum were obtained for further analysis. The time binning for the light curves was chosen on an observation-by-observation basis to optimally resolve the temporal structure of the source variability. The bin time of light curves was obtained by generating multiple light curves of different bin size and inspecting them. The time bin for each observation was chosen depending on the resolution that displayed characteristic features such as rise and decay of flares and the bin duration had to be long enough to ensure sufficient signal to noise ratio. The spectra were binned to minimum 25 counts per bin to ensure significant signal to noise ratio.
3. Methodology – Analysis techniques
Various methods have been used to analyse the flux and spectral variability of X-ray observations of BL Lacertae taken with the PN instrument onboard the XMM-Newton satellite.
3.1. Flux variability analysis
The method of excess variance and fractional variability amplitude (Vaughan et al. Reference Vaughan, Edelson, Warwick and Uttley2003) was used to test the variability of observations as described in Noel et al. (Reference Noel, Gaur, Gupta, Wierzcholska, Ostrowski, Dhiman and Bhatta2022). Fractional root mean square variability amplitude (
$F_\mathrm{var}$
) is mathematically described as:
\begin{equation}F_\mathrm{var}=\sqrt{\frac{S^{2}-{\overline{\sigma^{2}_\mathrm{err}}}}{\overline{x}^2}}\end{equation}
where
$S^2$
is the variance and
$\sigma_{\mathrm{err}, i}$
is the corresponding uncertainty associated with every flux point. The uncertainty of
$F_\mathrm{var}$
is given by
\begin{equation}\mathrm{err}(F_\mathrm{var})=\sqrt{{\left(\sqrt{\frac{1}{2N}}\frac{\overline{\sigma^2_\mathrm{ err}}}{{\overline{x}^2}F_\mathrm{var}}\right)}^2+{\left(\sqrt{\frac{\overline{\sigma^2_\mathrm{ err}}}{N}}\frac{1}{\overline{x}}\right)}^2}\end{equation}
$F_\mathrm{var}$
along with
$err(F_\mathrm{var})$
serves as a quantification of variability. The variable observations were further analysed for time-dependent spectral analysis which is further explained in the following sections.
3.2. Hardness ratio analysis
The hardness ratio (HR) is used to study the spectral states of the source and to compare photon hardness and variability between flaring and non-flaring states. HR analysis is advantageous as it does not depend on predefined models and allows the use of small time bins with low photon counts, making it efficient for studying spectral variability. The variable observations were split into energy ranges of 0.3–2.0 keV (soft) and 2.0–10.0 keV (hard), using the standard 2.0 keV division to separate the light curve into soft and hard bands (Peretz & Behar Reference Peretz and Behar2018; Noel et al. Reference Noel, Gaur, Gupta, Wierzcholska, Ostrowski, Dhiman and Bhatta2022). To calculate the hardness ratio, same method as (Park et al. Reference Park, Kashyap, Siemiginowska, van Dyk, Zezas, Heinke and Wargelin2006) was used:
\begin{equation} HR = \frac{H}{S}\end{equation}
where H is the flux in the hard energy band and S is the flux in the soft energy band in the analysed time bin. To derive and plot the hardness ratio against intensity, lcurve task of the XRONOS program of HEASARC has been used.
3.3. Discrete correlation function
The soft and hard light curves of variable observations were analysed to study the correlation between each other using the method of discrete correlation function (DCF) (Edelson et al. 1988; Noel et al. Reference Noel, Gaur, Gupta, Wierzcholska, Ostrowski, Dhiman and Bhatta2022). Monte Carlo simulation was performed on artifically generated 1 000 light curves to assess the significance of the DCF peaks from which we estimated the 1, 2, and 3 confidence bands. The significance of the observed DCF peak was determined relative to these confidence levels.
To evaluate the time lag, the maxima of DCF was fitted with a Gaussian model (Zhang et al. Reference Zhang, Gupta, Gaur, Wiita, An, Gu, Hu and Xu2019). The Gaussian distribution is given by:
\begin{equation} \mathrm{DCF}(\tau) = A e^{-(\tau- \mu)/2\sigma^2}\end{equation}
where A is the maximum DCF value and
$\sigma$
is the fitted width of the Gaussian function. A limited range of time
$\tau$
is selected for each DCF function such that the DCF maxima lie in this range. The time lag is given by
${\unicode{x03BC}}$
which is the difference between the time at which the DCF maxima lies and the centre of time range
$\tau$
.
3.4. Spectral analysis
The spectral analysis of all five observations is done to understand the complex behaviour of BL Lacertae. For this the spectra, in the energy range of 0.3–10.0 keV, were obtained using evselect, a function to create spectra in the SAS software. The spectral fitting was done using XSPEC v12.12.0 which is a spectral fitting package included in the HEASoft software.
The models used for fitting are power-law (PL), logparabola (LP), broken power-law (BP), double power-law (DP), which is a combination of two power-law, and blackbody models. Each of the models used include galactic absorption.Footnote e The double powerlaw model used here is a phenomenological model describing the spectral shape and not a physical model. A log-parabola model was employed to characterise both the spectral slope and its curvature. We note that while a log-parabola is not the ideal tool for a simple hardness analysis, its parameters provide a robust way to track spectral changes when significant curvature is present and power-law model may fail to adequately fit the spectrum.
In all logparabola models, the pivot energy (pivotE) was fixed at the standard value of 1 keV (Ehlert et al. Reference Ehlert2023; Wierzcholska & Wagner Reference Wierzcholska and Wagner2025). This conventional choice provides a stable reference point for model normalisation and facilitates comparison with other studies (Komossa et al. Reference Komossa2021).
Each of these models was tested using the
$\chi^{2}$
method of fitting and compared to each other for best fit using Akaike Information Criterion (AIC) (Gaur, Mohan, & Pandey Reference Gaur, Mohan and Pandey2021). To ensure that the analysis is not influenced by small sample-bin sizes, the AIC has been adjusted accordingly and is given by AIC +
$\frac{2k(k+1)}{n-k-1}$
where n is the sample size and k is the number of free parameters in a model. This corrected AIC is referred as AICc in the paper.
Galactic absorption (
$n_H$
) has been taken into account for BL Lacertae because interstellar absorption is observed along with a molecular cloud in the line of sight of BL Lacertae (Ravasio et al. Reference Ravasio2002). Previous studies on BL Lacertae have used varied values of
$n_H$
. Ravasio et al. (Reference Ravasio2002) attempted to fix n
$_H$
to the Galactic hydrogen column density measured at 21 cm, which is 2
$ \times 10^{22}$
cm
$^{2}$
. Earlier, Lucas & Liszt (Reference Lucas and Liszt1993) reported the presence of a molecular cloud along the line of sight to BL Lacertae. To account for molecular Galactic absorption, Raiteri et al. (Reference Raiteri2009) froze n
$_H$
at 0.34
$\times 10^{22}$
cm
$^{2}$
. However, the residuals remained high, so this value was not adopted. Similarly, the value of 0.17
$\times 10^{22}$
cm
$^{2}$
, obtained from the HI4PI full-sky H I survey (HI4PI Collaboration et al. 2016), did not fit our data well. Using this HI4PI value resulted in a reduced
$\chi^{2}$
greater than 1, indicating a poor fit to our observations. The value of 0.27
$ \times 10^{22}$
cm
$^{2}$
as used by (D’Ammando Reference D’Ammando2021) was used for analysis. This two-step approach successfully identified the best-fit continuum model for four of the five observations. This model was then used to investigate the absorption by allowing the nH parameter to vary freely in a subsequent fit.
3.5. Interval-wise spectral analysis
Of the two observations, only those exhibiting extreme variability were selected. The light curve of each variable observation was divided into time intervals based on lower and higher photon counts, as shown in Figure 2. The spectrum for each interval was then generated separately. This approach, used by Ferrero et al. (Reference Ferrero, Wagner, Emmanoulopoulos and Ostorero2006) for the IBL source S5 0716+71, is effective for studying spectral variations on intraday timescales.
Light curves of the non-variable observations each with binning of 100 s (a) 0501660201, (b) 0501660301, and (c) 0501660401.
Light curves of the variable observations (a) 0504370401 with bin time of 1 100 s and (b) 0891800501 with bin time of 100 s.
Of the two variable observations, the observation 0891800501 (observed on 09-07-2021) corresponds to a flaring state, while 0504370401 (observed on 16-05-2008) corresponds to a non-flaring state. Thus, this interval-wise spectral analysis is useful for comparing observations at different flux states while also examining how the spectra change with flux on short timescales within an observation.
The broken power-law fitting was used to determine the break energy where the spectrum shows an upturn and the photon indices of the two power-laws change. Assuming the lower-energy power-law arises from synchrotron emission, the synchrotron contribution (Ferrero et al. Reference Ferrero, Wagner, Emmanoulopoulos and Ostorero2006) in each spectrum was estimated by calculating the flux below the break energy and expressing it as a percentage of the total flux in the 0.3–10.0 keV energy range.
4. Results
4.1. Flux variability
The results of
$F_\mathrm{var}$
are given in Table 1 along with the observation details. Out of 5 total observations of BL Lacertae taken with XMM-Newton until August 2021, only two observations, 0504370401 (date 16-05-2008) and 0891800501 (date 09-07-2021) were found to be extremely variable. The observations that we find non-variable here were also analysed by Raiteri et al. (Reference Raiteri2009) and they found the same non-variable behaviour. The observation ID 0504370401 (date 16-05-2008) is the observation with highest variability out of five with
$F_\mathrm{var}$
value of 19.16%. The observations 0501660201 (date 10-07-2007) and 0501660301 (date 05-12-2007) have extremely low levels of variability with
$F_\mathrm{var}$
value lesser or almost equal to 3 times the
$\sigma_\mathrm{err}$
. The observation 05016601401 (08-01-2008) has an undefined
$F_\mathrm{var}$
and therefore variability values are not listed in Table 1. This can occur if the source shows very low intrinsic variability during the observation or due to low photon count statistics. However, this observation has a count rate of
$3.42 \pm 0.05$
counts/s, indicating that the photon statistics are sufficient, and the undefined
$F_\mathrm{var}$
is likely due to a lack of significant variability during this observation.
The highest variable observation also happens to be the observation with the longest duration, that is, 133.9 ks. The other variable observation that was taken on 9 July 2021, which is 0891800501, coincides with the flaring period of BL Lacertae which started in 2020 and lasted over 2021. The light curves of the non-variable and the variable observations are shown in Figures 1 and 2, respectively. From the visual inspection as well, it is somewhat clear that only two observations shown in Figure 2 display variability. Also from the light curves, it is clear that the observation taken in 2021 which coincides with the flaring stage of the source displays the highest flux count out of the five observations. Later we see a comparison between the spectra of the source in different variable regions of the variable observation. Also, the analysis of hardness ratio and lag in emission is studied of flaring (0891800501; 09-07-2021) and non-flaring observation (0504370401; 16-05-2008).
4.2. Characteristic time scale
The timescale of variability is an important parameter that can help us to derive some of the jet parameters such as the size of the emission region. We estimated the variability timescale for all the XMM-Newton EPIC-PN observations and the formula used is
\begin{equation} F(t_2) = F(t_1) \times 2^{(t_2 - t_1)/T_d}\end{equation}
Here
$F(t_1$
) and
$F(t_2$
) are the observed fluxes at time
$t_1$
and
$t_2$
, respectively.
$T_d$
is the flux doubling/halving or variability timescale. From the variability estimate, we see that the observations 0504370401 (date 16-05-2008) and 0891800501 (date 09-07-2021) are more variable and hence can probe the smaller timescale. The variability timescale estimated for these two observations is 1.28 ks or 21.33 min and 4.35 ks or 1.2 h, respectively. A minute and hour scale of variability in X-ray is common in blazar which puts tight constraints on the size of the emission region.
BL Lacertae is a complex source that sometimes behaves as an LBL type source and sometimes like an HBL type (Prince Reference Prince2021) which suggests the synchrotron emission can produce the X-ray emission and hence the X-ray decay timescale can be linked with the radiation cooling or synchrotron timescale (Homan et al. Reference Homan, Wijnands, van der Klis, Belloni, van Paradijs, Klein-Wolt, Fender and Méndez2001; Belloni, Psaltis, & van der Klis Reference Belloni, Psaltis and van der Klis2002; Uttley et al. 2014).
The synchrotron cooling timescale can be defined as (Rybicki & Lightman Reference Rybicki and Lightman1979):
\begin{equation} t_\mathrm{cool} \simeq 7.74\times10^8 \frac{(1+z)}{\delta}B^{-2}\gamma^{-1} ~\mathrm{s}.\end{equation}
where B is the strength of the magnetic field in Gauss and
$t_{\mathrm{cool}}$
is the synchrotron cooling timescale in seconds. Following Rybicki & Lightman (Reference Rybicki and Lightman1979), we can also derive the characteristic frequency of the electron population, responsible for the synchrotron emission at the SED peak,
\begin{equation} \nu_{\mathrm{ch}, e} = 4.2\times10^6 \frac{\delta}{(1+z)}B \gamma^{2} ~\mathrm{Hz}.\end{equation}
Using the above two equations, we can eliminate the
$\gamma$
since it changes with different states and can derive a single equation given as:
\begin{equation} B^3\delta \simeq 2.5 (1+z)(\nu_{\mathrm{ch}, e}/10^{18})^{-1} \tau_d^{-2}.\end{equation}
where
$\tau_d$
is in ks. Using the above equation, we derive the magnetic field strength for the Doppler factor (obtained from Acciari et al. Reference Acciari2021; Prince Reference Prince2021),
$\delta$
= 20, and a variability time scale of 1.28 ks, as 0.4 G. For a range of Doppler factor between 10 and 30, magnetic field strength takes values between 0.5 and 0.35 G.
Spectral fitting of each observation. Each observation is corrected for galactic absorption model Tbabs of value 2.7
$\times$
$10^{21}\,{cm}^{-2}$
. Here PL stands for power-law, LP for logparabola, BP for broken power-law and DP for double power-law.
$\Gamma_{1}$
is the photon index of PL and first PL, and the first photon index of BP. In the case of LP, this becomes
$\alpha$
.
$\Gamma_{2}$
is the second photon index of BP, and photon index of second PL in DP. In the case of LP, there is
$\beta$
which gives the curvature of the model.
$E_\mathrm{break}$
is the breaking energy of the broken power-law. Flux is in the unit of 10
$^{-12}\,\mathrm{ergs}\,{cm}^{-2}$
s
$^{-1}$
in the energy range of 0.3–10.0 keV. The values where the error bars are high are mentioned as higher in the table.
${}^*$
$\Gamma_2$
was fixed at
$-2.06$
because it was not constrained by the fit for observation 0501660201.
4.3. Spectral variability
The spectral fitting of each observation with different spectral models is shown in Table 2. We notice the following spectral behaviour in the XMM-Newton observations of BL Lacertae:
-
• In the first observation, 0501660201 (observed on 10-07-2007), the
$\chi^{2}$
fitting is below the expected value of 1. The AICc value is smallest for powerlaw followed by logparabola model. These two models appear to fit the observation 0501660201 quite well compared to broken power-law and double power-law which also show a huge error against
$\chi^{2}$
fitting. The broken powerlaw model has the lowest AICc value compared to all the other models fitted to 0501660401 (observed on 08-01-2008) and for the four models fitted to 0501660301 (observed on 05-12-2007). -
• The variable observation 0504370401 (observed on 16-05-2008) is well-described by broken powerlaw models as explained by the lowest AICc value. The AICc values of models fitted to the variable observation 0891800501 (observed on 09-07-2021) shows double powerlaw model to be the best model describing the X-ray emission. The reduced
$\chi^{2}$
of the fitted broken power-law model suggests that it is also a good model. Thus, we pick broken powerlaw model as a good choice to evaluate the upturn energy in SED at which synchrotron hump gradually moves to inverse Compton hump. -
• Each of the observations were also well-fitted by logaparbola model, described by a negative curvature, with
$\chi^{2}$
value within acceptable limit. The negative curvature obtained from this model pinpoints the existence of the upturn from synchrotron to inverse Compton emission in the SED of BL Lacertae in 0.3–10.0 keV energy range. The logparabola photon index given by alpha is higher in the flaring observation than in other observations giving an indication of the observation going in a soft state during the flare in 2021. The value of photon index for 0891800501 (date 09-07-2021) is above 2.5 while for all other observations it is below 2.5. To confirm the softer when brighter behaviour, the spectral parameter of the logparabola model was plotted against flux in Figure 3. Softer when brighter behaviour is evident in BL Lacertae. An increased negative curvature (
$\beta = -0.77$
) has been observed in the brightest observation which is consistent with the softer when brighter behaviour and indicating a steep high-energy tail from inverse Compton contribution in BL Lacertae’s X-ray spectrum. The modelling with double power-law model also helped identify the presence of both synchrotron and inverse Compton emission processes playing a role in the emission in X-ray energy range of 0.3–10.0 keV (Wierzcholska & Wagner Reference Wierzcholska and Wagner2025). -
• All observations also showed good fit with blackbody model and have acceptable reduced
$\chi^{2}$
values of closer to 1 (see Table 3). Addition of blackbody model was useful for getting the value of temperature kT of the disc which has been used to understand its change with varying flux. A positive trend of kT is observed with flux and was found highest for the observation 0891800501 observed to be flaring in 2021 (see Figure 8). -
• The best fit models were obtained with the galactic absorption nH kept frozen at 2.7
$\times$
$10^{21}\,{cm}^{-2}$
which is in accordance with the previous work of D’Ammando Reference D’Ammando2021 and takes into account the intervening molecular cloud in the line of sight of BL Lacertae. To check for any additional absorption or deviation from the frozen absorption value, the best fitted broken powerlaw model was used. The parameters of broken powerlaw model for each absorption were kept frozen to the best fit parameters as obtained in Table 2 and nH was left to vary. Galactic absorption nH was found to be 2.71
$\times$
$10^{21}\,{cm}^{-2}$
for the three observations – 0501660201, 0501660301, and 0501660401. Galactic absorption nH was found 2.70
$\times$
$10^{21}\,{cm}^{-2}$
and 2.68
$\times$
$10^{21}\,{cm}^{-2}$
for observations 0504370401 (observed on 16-05-2008) and 0891800501 (observed on 09-07-2021), respectively.
Photon index from logparabola model plotted against flux in the total energy band.
Spectral fitting of the five observations with blackbody model along with Galactic absorption Tbabs fixed at 2.7
$\times$
$10^{21}$
cm
$^{-2}$
and logparabola model. For a good fitting, a power-law model had to be added to the two variable observations: 0504370401 and 0891800501.
Observation 0504370401 (a) light curves in soft (0.3–2.0 keV) and hard band (2.0–10.0 keV) with 1 100 s time binned, (b) hardness ratio against time, (c) hardness ratio against intensity in which the colourbar represents the timescale in seconds from beginning of observation to end of it, and (d) discrete correlation function estimated between hard and soft energy range.
4.4. Soft and hard bands – variability and correlation
The two variable observations were analysed for change of hardness ratio with time and flux. The correlation between the soft and hard bands was tested using DCF. Below are the results for the two variable observations:
0504370401, In this non-flaring variable observation (observed on 16-05-2008), there is no trend of hardness ratio observed with time (Figure 4b). The emission gets softer with increasing flux in this observation which is evident in Figure 4(c). The DCF peak reaches a significance of 1.5 above the 1 confidence band but does not exceed the 2 or 3 thresholds derived from our Monte Carlo simulations. This indicates that the apparent correlation is not statistically significant and may be consistent with noise.
Observation 0891800501 (a) light curves in soft (0.3–2.0 keV) and hard band (2.0–10.0 keV), (b) hardness ratio against time, (c) hardness ratio against intensity, and (d) discrete correlation function estimated between hard and soft energy range.
0891800501 In this flaring variable observation (observation date 09-07-2021), softer when brighter behaviour is observed in the hardness ratio versus flux plot, given in Figure 5(c). The light curves in the soft and hard band also do not show any obvious correlation except that after 14 000 s of the start time of observation, the two light curves show a slight increase in flux count. While this nature of softer when brighter is consistent with Swift-XRT observations of BL Lacertae during its October 2020 flare (Prince Reference Prince2021), this is opposite to what was observed in the source during its flare in 2012 (Sahakyan & Giommi Reference Sahakyan and Giommi2022). There is no trend in hardness ratio with respect to time. The DCF plot, Figure 5(d), is randomly distributed and did not show any strong correlation between the two wavebands. DCF peak shows a marginal significance of 1.3 above the 1 confidence band but does not exceed the 2 or 3 levels in Monte Carlo simulations performed, indicating that the correlation is not statistically significant.
4.5. Interval-wise spectral variability
The two variable observations- 0504370401 (date 16-05-2008) and 0891800501 (date 09-07-2021) were split into shorter time intervals. The intervals of each observation are shown in their light curve in Figure 2. The observation 0504370401 is of longer time duration and its light curve is showing many variable states. Hence the observation was split into five different sections. Observation 0891800501 is a shorter observation with 19.6 ks. This observation was split into two sections. The first interval of this observation shows no visible variability in the light curve, whereas the second interval exhibits an increasing flux. Spectral fitting was performed separately for each interval.
From the previous section, we are aware of the concave structure of the x-ray spectra of BL Lacertae from the fitting of the logparabola to all the XMM-Newton observations of BL Lacertae. Also, the broken power-law was a good fit for all observations. Thus, we proceeded with fitting each interval with the broken power-law model which gives the value of break energy (
$E_{b}$
).
$E_{b}$
parameter of the broken power-law model gives the energy at which the SED transitions from the first hump to the second hump. Thus, the part before the break energy is from the contribution of synchrotron emission, and after the break energy lies the inverse Compton component. The spectral fitting of broken power-law fitted to the spectra of each section of observation is given in Table 4. Here it is visible that there is variability in spectral fitting in each section of the two variable observations. The break energy is changing. The flux for each section below and above break energy is given in Table 5. The fitted model to the observations are in online supplementary material.
Break energy from broken power-law is plotted against flux in the total energy band for each section of the two variable observations in Figure 7. In Figure 6, we plotted a variation of synchrotron contribution against flux and break energy obtained from broken power-law fitting against total flux respectively for each section of the variable observations. Synchrotron contribution is higher in the section with a higher flux value. Both the plots (Figures 6 and 7) show how the synchrotron contribution is high in observation with higher flux for this source. Thus, it can be said that synchrotron emission played a role in the X-ray flaring of BL Lacertae in 2021.
Spectral fitting of different sections of variable light curves with a broken power-law model with Galactic absorption Tbabs fixed at 2.7
$\times$
$10^{21}$
cm
$^{-2}$
. Obs ID is observation ID.
$\Gamma_1$
and
$\Gamma_2$
are the two-photon index of broken power-law. Break energy (
$E_b$
) is in the unit of keV. FLux is in the unit of
$10^{-12}$
ergs cm
$^{-2}$
s
$^{-1}$
.
$\chi^{2}$
is the reduced
$\chi^{2}$
value with dof as degrees of freedom.
Synchrotron contribution for each interval from the broken power-law model with frozen Tbabs.
Contribution from synchrotron emission plotted against flux over the entire energy range for each section of the two variable observations.
Break energy from broken power-law plotted against flux over the entire energy range for each section of the two variable observations.
5. Discussion
BL Lacertae exhibits complex behaviour, displaying variability ranging from low to extremely high levels. In our set of observations, two out of five (0504370401, taken on May 16, 2008, and 0891800501, taken on July 9, 2008) show high variability, while the remaining three display low variability. One of the two variable observation has a very long exposure time. This appears in accordance with the previous finding of a moderate correlation between exposure time and fractional variability by Noel et al. (Reference Noel, Gaur, Gupta, Wierzcholska, Ostrowski, Dhiman and Bhatta2022) in the blazar Markarian 421.
BL Lacertae has been observed to exhibit LBL characteristics (Nilsson et al. Reference Nilsson2018), IBL characteristics (Ackermann et al. Reference Ackermann2011; Hervet, Boisson, & Sol Reference Hervet, Boisson and Sol2016), and HBL characteristics (Prince Reference Prince2021) at different times. The negative curvature of the log-parabola model and the soft-to-hard spectral break of the broken power-law model both indicate a concave spectrum, suggesting the presence of two distinct emission components or a spectral minimum in the 0.3–10.0 keV band. The broken power-law model shows a softer first photon index compared to the harder second one, with the break energy ranging from 1.81 to 2.58 keV across the observations. A spectral break in this range implies that the synchrotron peak frequency is located at or below these energies. This behaviour is characteristic of BL Lacertae acting as an Intermediate-peaked BL Lac (IBL) object during these epochs.
This source has been well studied in the past, with several works aiming to estimate the exact location of the synchrotron peak (Raiteri et al. Reference Raiteri2009; Nilsson et al. Reference Nilsson2018; Prince Reference Prince2021). Raiteri et al. (Reference Raiteri2009) compiled the broadband spectrum using simultaneous observations across the wavebands (from radio to X-ray) and showed that the synchrotron peak lies in the near-IR band. They also noted that the X-ray spectrum is highly variable, changing from extremely steep to extremely hard during different flux states.
The Galactic absorption obtained from the analysis is significantly higher than the survey value of galactic absorption obtained in HI4PI survey (HI4PI Collaboration et al. 2016). The excess absorption suggests the presence of transient absorbing material – a molecular cloud along the line of sight to the emission region. This molecular cloud has been extensively studied earlier using
${}^{12}$
CO observations Bania et al. (Reference Bania, Marscher and Barvainis1991) and Liszt & Lucas (Reference Liszt and Lucas1998). This value is in compliance with the work of D’Ammando (Reference D’Ammando2021) and much lower than Raiteri et al. (Reference Raiteri2009).
(a) Variation of blackbody temperature kT given in keV with respect to flux. (b) kT plotted against time duration of observation. (c) Variation of kT with the date of observation.
The addition of a blackbody component was statistically required to model a prominent soft X-ray excess in the spectrum. A similar feature was noted by Raiteri et al. (Reference Raiteri2009), who argued that its variability pointed to a dynamic origin. Since the non-thermal jet emission was bright during the observation 0891800501 (observed on 09-07-2021), this soft excess likely represents an additional emission component rather than the unveiling of a faint, underlying disc. Future multi-wavelength campaigns incorporating optical/UV monitoring and X-ray polarimetry will be crucial to distinguish between an accretion-based or jet-based origin for this excess.
In this study, we found that as the source transitions from a low-flux state to a high-flux state, the photon index changes from 2.2 to 2.9, indicating softer-when-brighter behaviour as shown in Figure 3. Similar results were also reported by Prince (Reference Prince2021) using Swift-XRT spectra. While BL Lacertae has previously displayed the ‘harder-when-brighter’ behaviour typical of blazars, it showed a ‘softer-when-brighter’ trend in the two variable observations examined here, consistent with the findings of Prince (Reference Prince2021) during the 2020–2021 flare of BL Lacertae.
The contribution of synchrotron emission in the 0.2–10.0 keV energy range appears to increase with flux (Figure 6). Additionally, the correlation between the soft and hard bands showed no significant correlation in the variable observations. The contribution of the synchrotron component in the 0.3–10.0 keV range increased during the flaring state, indicating that enhanced synchrotron emission is one of the drivers of the source’s flaring activity. Since the soft and hard sub-energy bands show no correlation, the synchrotron-driven flaring may be attributed to factors unrelated to single emission zone such as interaction with external medium, geometric effects such as precession, or multiple emission regions.
BL Lacertae was observed to soften during the flaring state. This could be explained by the extension of synchrotron emission to higher energies, with the inverse Compton component remaining constant. This behaviour was also verified by the correlation of break energy with flux. A similar observation was made for S5 0716+71 by Giommi et al. (Reference Giommi1999), Ferrero et al. (Reference Ferrero, Wagner, Emmanoulopoulos and Ostorero2006), and Tagliaferri et al. (Reference Tagliaferri2003).
IBLs are complex sources exhibiting variable spectral behaviours over time. Synchrotron emission appears to be a dominant mechanism driving the emission during the flaring in 0.3–10.0 keV in both BL Lacertae and S5 0716+71 (Ferrero et al. Reference Ferrero, Wagner, Emmanoulopoulos and Ostorero2006). The concave structure observed in the log-parabola model fitted to BL Lacertae was also detected in the intermediate BL Lac object S5 0716+71 (Ferrero et al. Reference Ferrero, Wagner, Emmanoulopoulos and Ostorero2006). This feature was attributed to the transition from the tail of the synchrotron bump into the lower-energy region of the inverse Compton peak is similar to the log-parabola model fitted to BL Lacertae. The IBL S5 0716+71 also exhibited strong variability on timescales of hours. BL Lacertae exhibits similar behaviour to other IBLs, such as S5 0716+71 and ON 231. However, it is also important to recognise its demonstration of both LBL and HBL behaviour, as highlighted by previous studies. Only a detailed multiwavelength study of the SED of BL Lacertae conducted over an extended period of at least a decade (given its flaring pattern every 10 yr) can fully explain its complex behaviour and determine the exact classification of this blazar.
Acknowledgement
The project is co-financed by the Polish National Agency for Academic Exchange.
Appendix
Top: Spectra fitted with broken power-law model along with the residuals. Bottom: SED plot for Section I of observation 0504370401.
Same as for Figure A1 for Section II of observation 0504370401.
Same as for Figure A1 for Section III of observation 0504370401.
Same as for Figure A1 for Section IV of observation 0504370401.
Same as for Figure A1 for Section I of observation 0891800501.
Same as for Figure A1 for Section II of observation 0891800501.
Same as for Figure A1 for entire observation 0891800501.
Same as for Figure A1 for Section V of observation 0504370401.
Same as for Figure A1 for entire observation 0504370401.
Open access
