https://orcid.org/0000-0001-5507-7660
ABSTRACT
Variable |$\gamma$|-ray flares up to minute time-scales reflect extreme particle acceleration sites. However, for high-redshift blazars, the detection of such rapid variations remains limited by current telescope sensitivities. Gravitationally lensed blazars serve as powerful tools to probe |$\gamma$|-ray production zones in distant sources, with time delays between lensed signals providing crucial insights into the spatial distribution of emission regions relative to the lens’s mass-weighted centre. We have utilized 15 yr of Fermi–LAT |$\gamma$|-ray data from direction of PKS 1830–211 to understand the origin of flaring high-energy production zone at varying flux states. To efficiently estimate the (lensed) time delay, we used a machine learning-based tool – the Gaussian Process regression algorithm, in addition to – Autocorrelation function and Double power spectrum. We found a consistent time delay across all flaring activity states, indicating a similar location for the |$\gamma$|-ray emission zone, possibly within the radio core. The estimated time delay of approximately 20 d for the five flaring epochs was significantly shorter than previously estimated radio delays. This suggests that the |$\gamma$|-ray emission zone is closer to the central engine, in contrast to the radio emission zone, which is expected to be much farther away. A linear relationship between lag and magnification has been observed in the identified source and echo flares. Our results suggest that the |$\gamma$|-ray emission zone originates from similar regions away from the site of radio dissipation.
1 INTRODUCTION
The accretion of matter onto supermassive black holes powers Active Galactic Nuclei (AGNs). The interplay of magnetic field and rotation either of the black hole (Blandford & Znajek 1977) or of the accretion disc (Blandford & Payne 1982) are believed to generate collimated plasma jets that extend from the central engine to large distances. A subset of these jetted AGNs, known as Blazars, are aligned with our line of sight and exhibit extremely luminous and highly variable emissions across a broad electromagnetic spectrum. Notably, Multifrequency variability of these point-jetted AGNs provides insights into the size of the emission region in the jet (Madejski & Sikora 2016). The detection of extremely rapid variability, on time-scales comparable to the light-crossing time of the black hole, in both high-energy (HE) and very high-energy (VHE) emissions suggests that the emission zone is compact and located close to the black hole, within a few tens of gravitational radii (Aleksić et al. 2011; Ackermann et al. 2016; Shukla et al. 2018; Agarwal et al. 2023). However, HE and VHE |$\gamma$|-ray photons originating near the black hole are highly susceptible to attenuation from photon–photon pair production due to the dense field of ultraviolet (UV) and optical seed photon (Liu & Bai 2006). This suggests that |$\gamma$|-ray emission must occur farther from regions dominated by such external seed photons. None the less, the lack of sufficient seed photons at larger distances from central engine complicates the understanding of where high-energy dissipation occurs.
Variable radio emissions on parsec (pc) to megaparsec (Mpc) scales often correlate with |$\gamma$|-ray emissions, suggesting co-spatial origins within the jet (Marscher et al. 2008; Ghirlanda et al. 2011). However, the absence of rapid variability in radio bands and potential synchrotron self-absorption up to hundreds of GHz limit the detection of radio emissions in smaller structures (Rybicki & Lightman 1979). Observations across radio to X-ray wavelengths indicate radiation regions spanning subparsec to megaparsec scales (Harris & Krawczynski 2006; Tavecchio et al. 2007; Marscher et al. 2008; Fuentes et al. 2023). The limited resolution of high-energy telescopes further complicates identification of |$\gamma$|-ray emission regions, leaving fast variability as the primary probe for high-energy processes. This limitation constrains our understanding of the link between radio and |$\gamma$|-ray variability (Blandford & Levinson 1995; Jorstad et al. 2001). The consistency of radiation sources across energy levels remains uncertain, particularly in high-redshift blazars, where rapid variability detection is hindered by current telescope sensitivities.
Another probe of the |$\gamma$|-ray production regions is gravitational lensing. One of the earliest predictions of Einstein’s general relativity was the deflection of light by the Sun (Einstein 1936). Later, Zwicky (1937a, b) proposed that galaxies, like stars, could also act as gravitational lenses. In such a lensed system, photons from a background galaxy are bent around a foreground lensing galaxy, creating a magnified and distorted image of the background source. Radiation from the same point on the source follows different paths, causing time-variable sources to show similar variability patterns in different lensed images, but with time delays and magnification. These time delays and magnifications depend on the geometry of the source-lens-observer system.
Time delays and magnifications across energy ranges directly reflect the size of the emission region and the distribution of the emission zone around the mass-weighted centre of the lens (Barnacka et al. 2014). Distribution of such time delays provide an alternative approach to understanding the origin of |$\gamma$|-rays in blazars. Continuous observations by telescopes like Fermi enable long-term monitoring of distant sources, providing insights into the evolution of flux variability over time. Among known gravitationally lensed quasars, two have been detected at |$\gamma$|-ray energies: PKS 1830–211 (Abdo et al. 2015) and QSO B0218+357 (Cheung et al. 2014).
PKS 1830–21 was first identified as a gravitationally lensed system by the Very Large Array Radio Telescope, revealing two compact components in the northeast and southwest (Subrahmanyan et al. 1990). The Australian Telescope Compact Array later showed these double components separated by 0.98 arcmin and connected by an elliptical Einstein ring (Jauncey et al. 1991; Nair, Narasimha & Rao 1993). The flat spectrum radio quasar PKS 1830–211 (|$z=2.507$|), which would typically appear as a point source, exhibits a double radio structure indicative of an intervening lens galaxy at redshift |$z=0.89$| (Wiklind & Combes 1996; Winn et al. 2002; Koopmans & de Bruyn 2005). Evidence also suggests a second intervening galaxy at |$z=0.19$|, with H i and OH absorption, though its effect on lensing is expected to be negligible (Nair et al. 1993; Lovell et al. 1996; Winn et al. 2002; Muller et al. 2020).
Radio observations with the Australia Telescope Compact Array at 8.6 GHz measured radio time delays of |$26^{+4}_{-5}$| d (Lovell et al. 1998). Within two years of Fermi/LAT observations, the first gravitational time delay in |$\gamma$|-rays during its quiescent state was measured at |$27.1\pm 0.6$| d (Barnacka, Glicenstein & Moudden 2011). The consistency of time delays in |$\gamma$|-rays and radio bands may suggest a co-spatial origin during low states of |$\gamma$|-ray activity (Barnacka et al. 2014). Later searches in Fermi–LAT during active states revealed shorter time delays of |$23\pm 0.5$| and |$19.7\pm 1.2$| d (Barnacka et al. 2015). An independent study using molecular absorption lines derived a differential time delay of |$24^{+5}_{-4}$| d, with the north-east component leading (Wiklind & Combes 2001).
The search for time delays relies on the length of the light curve. Understanding the origin of |$\gamma$|-ray flares can be enhanced by studying various flaring states at different flux levels. PKS 1830–211 has shown significant activity over the past decade, with multiple flares detected in the Fermi–LAT light curve. We estimated the time delays during such high flaring periods using the autocorrelation function and double power spectrum. We also used a machine-learning technique called Gaussian Process regression to estimate time delays in different flux states. A comprehensive 15-yr search for time delays was conducted, focusing on the active states of the source. The paper is structured as follows: Section 2 describes the data analysis and the tools and techniques used for time delay estimation. Sections 3 and 4 present the results and discussion, respectively. The summary of our results is provided in Section 5.
2 METHODS AND TECHNIQUES
2.1 Data Reduction – Fermi–LAT
Despite being a lensed quasar, PKS 1830–211 appears as a point source due to the limited spatial resolution of the Fermi–LAT telescope. Regardless of this limitation, the Fermi–LAT telescope proves to be useful. It allows us to utilize the combined flux from the two anticipated lensed images, which have now coalesced into a signal exhibiting a distinct time delay. These delays are expected to provide valuable constraints on the emission size of the source.
Fermi–LAT is a pair-conversion |$\gamma$|-ray telescope, sensitive to photons in the energy range of 20 MeV–300 GeV (Atwood et al. 2009). We selected data from 15.5 yr of Fermi observations, covering the period from MJD 54 683 to MJD 60373, within 10|$^{\circ }$| of the location of PKS 1830–211. This source, located approximately |${\sim}5^{\circ }$| from the galactic centre, is prone to galactic contamination. To minimize this, we selected photons with energies within 0.2–300 GeV.
To extract the photon statistics, the standard analysis procedure suggested by the FermiScience Tools and the open-source fermipy package (Wood et al. 2017) was used in the energy range between |$0.2\!-\!300$| GeV, employing the latest instrument response function P8R3_SOURCE_V3. A zenith angle cut of |$90^\circ$|, a GTMKTIME cut of DATA_QUAL> 0 && LAT_CONFIG= = 1, and evtype = 3 were used in the analysis. Only those events highly probable of being photons were considered for further analysis by applying a GTSELECT cut on event class to account for SOURCE class events using evclass = 128.
A source model was prepared by including the source at RA = 278.413 and Dec. = −21.075 and considering all the 4FGL catalogue sources within |$20^{\circ }$| around the region of interest. The source is modelled with log-parabola model, parametrized as
where scale parameter |$E_\mathrm{ b}$| was fixed to 4FGL catalogue value of 645.56 MeV, |$\alpha$| is specral index, |$\beta$| is curvature parameter, and |$N_{o}$| is the Normalization. Spectral parameters for the sources within |$5^\circ$| of the region of interest were allowed to vary. Sources outside |$5^\circ$| were fixed to 4FGL catalogue values. Additionally, the background modelled with the diffuse galactic emission model (gll_iem_v07) and the extragalactic isotropic diffuse emission model for point source analysis (iso_P8R3_SOURCE_V3_v1) was allowed to vary.
We performed a binned likelihood analysis using GTLIKE to evaluate the best-fitting model parameters, including the source’s spectrum and intensity at different epochs. The significance of detection is quantified using the Test Statistic (TS), defined as |$\rm {TS}=-2 \ln (\mathcal {L}_0/\mathcal {L}_1)$|, where |$\mathcal {L}_0$| and |$\mathcal {L}_1$| are the likelihood values without and with the point source at the position of interest, respectively. Only significant epochs with TS > 9, predicted photons > 3, and bins with flux greater than its uncertainty (|$F_{\mathrm{ t}} \gt \sigma _{\mathrm{ t}}$|) are considered for further analysis.
2.2 Analysis tools and techniques
2.2.1 Bayesian Block and HOP algorithm
To enable the detection and characterization of localized variability structures over time, we represent flux points and their uncertainties as step-functions using Bayesian Block (BB; Scargle et al. 2013). Each point of change in the block in the BB representation highlights a 3|$\sigma$| variation from the previous block. The BB output is then processed by the HOP algorithm, which is based on a watershed concept derived from topological data analysis (Eisenstein & Hut 1998). HOP algorithm by itself identifies flaring states or periods of higher flux by clustering data points from neighbouring regions where the flux exceeds a threshold. The combination of BB and HOP identifies a block higher than those before and after it as a peak. This approach then traces down from the peak in both directions, stopping when each subsequent block is lower than the previous one. Here, we use the mean flux as a lower threshold.
This technique divides the light curve into flaring and quiescent epochs, with consecutive connected BBs above the mean flux baseline referred to as a HOP group. The flare identification code developed by Wagner et al. (2021) differentiates between various HOP groups, leading to the identification of nine HOP groups, represented by grey patches for our source in Fig. 1. The maximum time delay between the source and its lensed counterpart is expected to be approximately 70 d, as suggested by Barnacka et al. (2015). Out of the nine HOP groups represented in Fig. 1, some are less than 70 d apart. These close intervals suggest probable pairs of the source and its echo flare, arising from lensing. Therefore, we group HOP groups that are separated by less than 70 d, leading to five flaring states: F1 (MJD 55450–55600), F2 (MJD 56063, 56173), F3 (MJD 58363–58963), F4 (MJD 59063–59153), and F5 (MJD 59683–59943) as shown in Fig. 1. The lag and magnification from these flaring groups are further discussed.
10-d binned light curve of |$\approx$|15.5 yr of Fermi–LAT observation of gravitationally lensed FSRQ PKS 1830–211. The grey region represents the high-flux state and the white region represents the low-flux state. The light curve is divided into flaring epochs identified using HOP groups, marked by grey patches. HOP groups separated by less than 50 d are combined into flaring states, labelled as F1 to F5 (indicated by horizontal lines). The secondary y-axis (right) shows the detection significance as |$\sqrt{\rm {TS}}$|. Periods with TS < 9 are represented by upper limits.
2.2.2 Power spectrum
To identify the intrinsic temporal behaviour of the time series, we evaluate the power spectral density (PSD) of high-energy |$\gamma$|-ray light curves. For a stochastic time series, the power distribution at each frequency typically follows a power-law (PL) (|$P(f)\propto f^{-k}$|) across various wavelengths and time-scales, with an index ranging from approximately 1 to 3 (Nakagawa & Mori 2013; Finke & Becker 2014; Sobolewska et al. 2014). The average slope for the |$\gamma$|-ray PSD of the brightest 22 flat-spectrum radio quasars and 6 BL Lac objects is 1.5 and 1.7, respectively (Abdo et al. 2010). During the quiescent state, blazars typically exhibit temporal variability characterized by pink noise behaviour, with |$\alpha \sim 1$|.
We compute the power-law variability index for the |$1\,$| d and 12-h binned Fermi–LAT light curve (LC) for flares F1–F5 to quantify the temporal variability during the observed period, using the PSRESP method described in Max-Moerbeck et al. (2014), based on Uttley, McHardy & Papadakis (2002). The obtained PSD is fitted with a PL model of the form PSD|$(\nu)\propto \nu ^{-k}$|, where k and |$\nu$| are the spectral index and frequency, respectively. We simulate 1000 LCs with similar flux distribution and statistical variability as the observed LC using Connolly (2015). We have accounted for red noise leakage and aliasing effects as described in Goyal et al. (2022).
2.3 Estimating time delay
Gravitational lensing is often used as a promising tools for determining cosmological distances (Refsdal 1964; Blandford & Narayan 1992; Schechter et al. 1997). Additionally, the time delay and magnification ratio derived at any wavelngth can explicate the location of the emission region relative to the central black hole (Barnacka et al. 2014). Atwood (2007) predicted that LAT could detect delayed emission from bright lensed objects. High-energy observations of blazars exhibit significant variability due to the small emission region. The lensing-induced delay in photon arrival is expected to alter the intrinsic flux pattern of the source. Unlike radio and optical telescopes, which can resolve magnified, multiple images of a lensed source, high-energy observations are often limited by poor spatial resolution. Consequently, the composite flux from source and its echo image appears as a point source to Fermi. This results in a repeated flux pattern spaced by a days and demagnified by a factor of b in the time domain. The total observed flux can be expressed as
Thus, the total flux from the two images is integrated into the combined light curve when observed by high energy telescopes like Fermi. Thus, an added challenge is disentangling the repeated flares imprinted in the combined light curve from the apparent point source. Cheung et al. (2014) attempted to separate these flares and identified a leading and trailing component in the lensed blazar B0218+537.
In this work, we have used three techniques to estimate the lags in data using (1) autocorrelation Function, (2) double power spectrum, and (3) Gaussian process regression. Fig. 1 represents the 5 flaring epochs that have been explored in further work.
2.3.1 Autocorrelation function (ACF)
The autocorrelation function (ACF) is a standard statistical tool for assessing the similarity of a time series with a delayed copy of itself. ACF allows the identification of periodicity or repeated patterns in a signal, making it suitable for estimating lags in data.
For a noise-dominated signal, variable structures are inherently present in power-law noise with an index greater than 1. Barnacka et al. (2015) described the role of ACF in deciphering the lag in noisy lensed signals and found that time delay detection is easier for a source with large variability index (k). Thus, a steep spectrum with spurious peaks improves the chances of confident detection.
The significance of the estimated lags is evaluated using the Monte Carlo simulation described in Section 2.3.4. To make the time series continuous, the epochs of no or less significant observation are interpolated with zeros as in Barnacka et al. (2015). We then performed autocorrelation on 1-d and 12-h binned time series for flares F1 to F5. The results are discussed in Section 3.
2.3.2 Double Power Spectrum (DPS)
A lensed time series is represented in equation (2). Given the long, continuous, and evenly spaced nature of Fermi–LAT light curves, it is feasible to extract the lag from the data using the Fourier transform, as described by Barnacka et al. (2011, 2015) and Barnacka (2013). This method was first used for lag estimation in lensed light curves by Barnacka et al. (2011). The idea is to take a Fourier transform of the first power spectrum of equation (2). The Fourier transform of the first component, |$s(t)$|, is |$\tilde{s}(f)$|, and for the second component, it is |$\tilde{s}(f)e^{-2\pi i f a}$|. Therefore, the observed signal in the frequency domain is transformed into
The first power spectrum (FPS) is the square modulus of Fourier transform of |$S_{\rm {obs}}$|, i.e.
If an intrinsic lag is imprinted onto the signal, it should be encoded into the periodicity of the FPS with a time period inversely proportional to the lag. Therefore, a power spectrum of the FPS is expected to show a large amplitude of the time delay signal. Barnacka et al. (2011) found that DPS is 90 per cent efficient in detecting the encoded lag in the signal, significantly improving over the 10
To correct for the smearing effect caused by the finite length of the signal and sampling in the data, the signal must undergo specific processing steps. We use the method described in Barnacka et al. (2015), based on Brault & White (1971), to extract the time delay present in the signal. This method can accurately extract time delays regardless of whether the light curve is white noise or red noise dominated and eliminates fake time delay peaks expected in red noise signals.
2.3.3 Gaussian Process Regression (GPR)
A Gaussian process (GP) is a random process where any point x in the real domain is a random variable |$f(x)$|, and the joint distribution of a finite number of these variables follows a Gaussian distribution. Mathematically, for a set of inputs |$x_1,x_2$|,...,|$x_n$| with corresponding outputs |$y_1,y_2$|,...,|$y_n$|, wherein y = f(x), the function values f(x) follow a joint Gaussian distribution. GP can be seen as a generalization of the infinite-dimensional multivariate Gaussian distribution. In a finite-dimensional Gaussian distribution, the correlation between random variables is defined using a covariance matrix. For an infinite-dimensional Gaussian distribution, this matrix is replaced by a ‘covariance function’, known as a Kernel (Rasmussen 2004). The mean function of the infinite-dimensional Gaussian is typically set to zero for easier computation, but the mean of the observational data is later added to obtain predictions on the original scale, leveraging the scaling property of Gaussian distributions. Standardization, which involves subtracting the mean and dividing by the standard deviation of the data before fitting the GP, is a common practice.
Kernel selection: The choice of kernel requires some prior information about the data. In our analysis, we incorporate the prior knowledge that the data exhibit a lag effect. This leads to the selection of the following kernel
where l is the length scale, and p is the distance between repetitions. The first multiplicative element in this kernel corresponds to a Gaussian-shaped correlation function, while the second multiplicative element represents a periodic correlation. Fig. 2 depicts the shape of the chosen kernel after multiplication. In this context, the periodicity parameter p effectively functions as the lag.
(Left) Kernel visualization using covariance between each sample location and zeroth point for RBF, Periodic, and |$\rm {RBF}\times \rm {Periodic}$|. (Right) Covariance matrix of the sample space for |$\rm {RBF}\times \rm {Periodic}$| kernel where warmer colours indicate higher correlations.
Hyper-parameter estimation: The likelihood function can serve as an objective function for a non-linear optimization algorithm to obtain the maximum likelihood parameter values. In GPR, this is replaced by the log marginal likelihood, which incorporates both the data fit term and a penalty term to prevent overfitting. The log marginal likelihood consists of three terms added together: The first term, |$-\frac{1}{2} \mathbf {y}^T (\mathbf {K}(\mathbf {X}, \mathbf {X^{\prime }}) + \sigma _n^2 \mathbf {I})^{-1} \mathbf {y}$|, quantifies the quality of the fit. The second term, |$- \frac{1}{2} \log \det (\mathbf {K}(\mathbf {X}, \mathbf {X^{\prime }}) + \sigma _n^2 \mathbf {I})$|, helps avoid overfitting. The last term |$- \frac{n}{2} \log (2\pi)$| is a normalization term to ensure a valid probability distribution, where |$\mathbf {K(X,X^{\prime })}$| is the covariance matrix, |$\mathbf {I}$| is the identity matrix, and n is the number of data points.
We optimize the hyperparameters of the kernel |$\mathbf {K}$|. This study utilizes the scikit-learn GPR module (Pedregosa et al. 2011) for easier computation. We optimize the length scale hyperparameter for various fixed periodicity hyperparameter values and obtain the log marginal likelihood profile across different lag values, as shown in Fig. 8. We select a grid of lag values from 1 to 70 d in steps of 1 d. The upper limit for the lag is chosen based on prior information about the gravitationally lensed source (Zhang et al. 2008). Since the marginal likelihood can exhibit very similar values across different lags, we introduce a metric for better comparison, termed the ‘likelihood metric’. Given that the marginal likelihood can be negative, we multiply it by −1 and subtract the maximum of this value from each marginal likelihood. The optimal lag then corresponds to the maximum value of the likelihood metric. Given the probabilistic nature of the method, the estimated lag is expected to exhibit a distribution centred around the true lag value. The uncertainties in the derived lag are quantified by analysing the spread within this distribution.
2.3.4 Statistical significance
The significance of spurious peaks in ACF and DPS must be assessed to determine whether the observed time delay is a result of chance or represents an intrinsic time delay within the signal.
We simulated |$10^{5}$| light curves using the techniques described in Emmanoulopoulos, McHardy & Papadakis (2013) to generate artificial light curves with similar flux distribution and temporal variability as the observed light curve. This method allows for generating light curves with non-Gaussian distributions, overcoming a limitation of Timmer & König (1995). High-energy |$\gamma$|-ray light curves of blazars typically follow a lognormal flux distribution (Romoli et al. 2018; Bhatta 2021). As a result, the simulated light curves have identical statistical properties to the observed light curve.
To make the simulated light curves as realistic as possible, we included data gaps identical to the observed periods and interpolated them with zeros to ensure similar effects on the ACF and DPS of the simulated signal. We constructed cumulative probability distributions of the derived powers at each time delay. These distributions were then analysed for |$1\sigma$|, |$2\sigma$|, |$3\sigma$|, |$4\sigma$| chances of detection. Any significant (|$\gt 3\sigma$|) powers corresponding to a time delay are considered as intrinsic time delays in the signal.
3 RESULTS
The HE light curve of PKS 1830−211 appears quite complex, with multiple flaring periods over the quiet states. Fig. 1 shows significant variability in the 10-d binned high-energy (200 MeV–300 GeV) flux over time. This variability is evident from the fluctuating flux levels, with some periods displaying higher fractional variability than others (see Table 1). The flaring periods exhibit dominant pink noise behaviour with a PSD power-law index of approximately 1. Notably, a transition from pink to red noise behaviour is observed for the brightest flux state of the source, identified as F3 in this work (Table 1). We use these flaring states to identify dominant emission zones, which should appear as twin pairs of flares separated by a specific time interval for a lensed blazar. The flaring epochs, which are above the mean flux levels, are identified using BBs, indicated by grey patches in Fig. 1. Multiple blocks spaced less than 70 d apart are merged together, resulting in flaring states labelled F1, F2, F3, F4, and F5.
PSD results.
| Flux state|$^{1}$| | Time period|$^{2}$| | |$T_{\mathrm{ obs}}$||$^{3}$| | |$N_{\mathrm{ TS}\gt 9}/N_{\mathrm{ tot}}$||$^{4}$| | |$\Delta T_{\mathrm{ min}}$||$^{5}$| | |$\Delta T_{\mathrm{ max}}$||$^{6}$| | |$T_{\mathrm{ mean}}$||$^{7}$| | |$k \pm k_{\mathrm{ err}}$||$^{8}$| | |$p_{\beta }$||$^{9}$| | |$F_{\mathrm{ var}} \pm \Delta F_{\mathrm{ var}}$||$^{10}$| |
|---|---|---|---|---|---|---|---|---|---|
| [d] | [d] | [d] | [d] | [d] | |||||
| F1 | MJD 55450–55600 | 150 | |$123/150$| | 1 | 5.0 | 1.22 | |$0.91 \pm 0.39$| | 0.85 | 0.66 |$\pm$| 0.02 |
| |$168/360$| | 0.5 | 5.0 | 0.77 | |$0.96 \pm 0.24$| | 0.60 | 0.60 |$\pm$| 0.02 | |||
| F2 | MJD 56063–56173 | 110 | |$101/110$| | 1 | 3.0 | 1.08 | |$0.72 \pm 0.33$| | 0.05 | 0.46 |$\pm$| 0.03 |
| |$153/220$| | 0.5 | 4.0 | 0.71 | |$1.15 \pm 0.24$| | 0.70 | 0.37 |$\pm$| 0.03 | |||
| F3 | MJD 58363–58963 | 600 | |$517/600$| | 1 | 21.0 | 1.18 | |$1.36 \pm 0.24$| | 0.83 | 0.88 |$\pm$| 0.01 |
| |$941/1200$| | 0.5 | 21.0 | 0.64 | |$1.37 \pm 0.22$| | 0.54 | 0.84 |$\pm$| 0.01 | |||
| F4 | MJD 59063–59153 | 90 | |$75/90$| | 1 | 4.0 | 1.19 | |$0.88 \pm 0.47$| | 0.14 | 0.19 |$\pm$| 0.05 |
| |$105/180$| | 0.5 | 5.0 | 0.86 | |$0.50 \pm 0.27$| | 0.34 | 0.12 |$\pm$| 0.08 | |||
| F5 | MJD 59683–59943 | 260 | |$213/260$| | 1 | 22.0 | 1.22 | 1.27 |$\pm$| 0.37 | 0.72 | 0.53 |$\pm$| 0.02 |
| |$355/320$| | 0.5 | 21.5 | 0.73 | |$1.23\pm 0.20$| | 0.08 | 0.45 |$\pm$| 0.02 |
| Flux state|$^{1}$| | Time period|$^{2}$| | |$T_{\mathrm{ obs}}$||$^{3}$| | |$N_{\mathrm{ TS}\gt 9}/N_{\mathrm{ tot}}$||$^{4}$| | |$\Delta T_{\mathrm{ min}}$||$^{5}$| | |$\Delta T_{\mathrm{ max}}$||$^{6}$| | |$T_{\mathrm{ mean}}$||$^{7}$| | |$k \pm k_{\mathrm{ err}}$||$^{8}$| | |$p_{\beta }$||$^{9}$| | |$F_{\mathrm{ var}} \pm \Delta F_{\mathrm{ var}}$||$^{10}$| |
|---|---|---|---|---|---|---|---|---|---|
| [d] | [d] | [d] | [d] | [d] | |||||
| F1 | MJD 55450–55600 | 150 | |$123/150$| | 1 | 5.0 | 1.22 | |$0.91 \pm 0.39$| | 0.85 | 0.66 |$\pm$| 0.02 |
| |$168/360$| | 0.5 | 5.0 | 0.77 | |$0.96 \pm 0.24$| | 0.60 | 0.60 |$\pm$| 0.02 | |||
| F2 | MJD 56063–56173 | 110 | |$101/110$| | 1 | 3.0 | 1.08 | |$0.72 \pm 0.33$| | 0.05 | 0.46 |$\pm$| 0.03 |
| |$153/220$| | 0.5 | 4.0 | 0.71 | |$1.15 \pm 0.24$| | 0.70 | 0.37 |$\pm$| 0.03 | |||
| F3 | MJD 58363–58963 | 600 | |$517/600$| | 1 | 21.0 | 1.18 | |$1.36 \pm 0.24$| | 0.83 | 0.88 |$\pm$| 0.01 |
| |$941/1200$| | 0.5 | 21.0 | 0.64 | |$1.37 \pm 0.22$| | 0.54 | 0.84 |$\pm$| 0.01 | |||
| F4 | MJD 59063–59153 | 90 | |$75/90$| | 1 | 4.0 | 1.19 | |$0.88 \pm 0.47$| | 0.14 | 0.19 |$\pm$| 0.05 |
| |$105/180$| | 0.5 | 5.0 | 0.86 | |$0.50 \pm 0.27$| | 0.34 | 0.12 |$\pm$| 0.08 | |||
| F5 | MJD 59683–59943 | 260 | |$213/260$| | 1 | 22.0 | 1.22 | 1.27 |$\pm$| 0.37 | 0.72 | 0.53 |$\pm$| 0.02 |
| |$355/320$| | 0.5 | 21.5 | 0.73 | |$1.23\pm 0.20$| | 0.08 | 0.45 |$\pm$| 0.02 |
Note. (1) Flux states extracted from use of BB and HOP algorithm (2) Period of the flaring states (3) Total exposure time (4) Fraction of points having TS > 9 (5) Minimum sampling time in observed LC (6) Maximum sampling time in observed LC (7) Mean Sampling time i.e. total observation time over a number of data points in that interval (8) The power-law index for the power-law model of PSD analysis using PSRESP method (9) p-value corresponding to the power-law model. The power-law model is considered a bad fit if |$p_{\beta }\le 0.1$| as the rejection confidence for such model is |$\gt 90~{{\ \rm per\ cent}}$| (10) Fractional variability
PSD results.
| Flux state|$^{1}$| | Time period|$^{2}$| | |$T_{\mathrm{ obs}}$||$^{3}$| | |$N_{\mathrm{ TS}\gt 9}/N_{\mathrm{ tot}}$||$^{4}$| | |$\Delta T_{\mathrm{ min}}$||$^{5}$| | |$\Delta T_{\mathrm{ max}}$||$^{6}$| | |$T_{\mathrm{ mean}}$||$^{7}$| | |$k \pm k_{\mathrm{ err}}$||$^{8}$| | |$p_{\beta }$||$^{9}$| | |$F_{\mathrm{ var}} \pm \Delta F_{\mathrm{ var}}$||$^{10}$| |
|---|---|---|---|---|---|---|---|---|---|
| [d] | [d] | [d] | [d] | [d] | |||||
| F1 | MJD 55450–55600 | 150 | |$123/150$| | 1 | 5.0 | 1.22 | |$0.91 \pm 0.39$| | 0.85 | 0.66 |$\pm$| 0.02 |
| |$168/360$| | 0.5 | 5.0 | 0.77 | |$0.96 \pm 0.24$| | 0.60 | 0.60 |$\pm$| 0.02 | |||
| F2 | MJD 56063–56173 | 110 | |$101/110$| | 1 | 3.0 | 1.08 | |$0.72 \pm 0.33$| | 0.05 | 0.46 |$\pm$| 0.03 |
| |$153/220$| | 0.5 | 4.0 | 0.71 | |$1.15 \pm 0.24$| | 0.70 | 0.37 |$\pm$| 0.03 | |||
| F3 | MJD 58363–58963 | 600 | |$517/600$| | 1 | 21.0 | 1.18 | |$1.36 \pm 0.24$| | 0.83 | 0.88 |$\pm$| 0.01 |
| |$941/1200$| | 0.5 | 21.0 | 0.64 | |$1.37 \pm 0.22$| | 0.54 | 0.84 |$\pm$| 0.01 | |||
| F4 | MJD 59063–59153 | 90 | |$75/90$| | 1 | 4.0 | 1.19 | |$0.88 \pm 0.47$| | 0.14 | 0.19 |$\pm$| 0.05 |
| |$105/180$| | 0.5 | 5.0 | 0.86 | |$0.50 \pm 0.27$| | 0.34 | 0.12 |$\pm$| 0.08 | |||
| F5 | MJD 59683–59943 | 260 | |$213/260$| | 1 | 22.0 | 1.22 | 1.27 |$\pm$| 0.37 | 0.72 | 0.53 |$\pm$| 0.02 |
| |$355/320$| | 0.5 | 21.5 | 0.73 | |$1.23\pm 0.20$| | 0.08 | 0.45 |$\pm$| 0.02 |
| Flux state|$^{1}$| | Time period|$^{2}$| | |$T_{\mathrm{ obs}}$||$^{3}$| | |$N_{\mathrm{ TS}\gt 9}/N_{\mathrm{ tot}}$||$^{4}$| | |$\Delta T_{\mathrm{ min}}$||$^{5}$| | |$\Delta T_{\mathrm{ max}}$||$^{6}$| | |$T_{\mathrm{ mean}}$||$^{7}$| | |$k \pm k_{\mathrm{ err}}$||$^{8}$| | |$p_{\beta }$||$^{9}$| | |$F_{\mathrm{ var}} \pm \Delta F_{\mathrm{ var}}$||$^{10}$| |
|---|---|---|---|---|---|---|---|---|---|
| [d] | [d] | [d] | [d] | [d] | |||||
| F1 | MJD 55450–55600 | 150 | |$123/150$| | 1 | 5.0 | 1.22 | |$0.91 \pm 0.39$| | 0.85 | 0.66 |$\pm$| 0.02 |
| |$168/360$| | 0.5 | 5.0 | 0.77 | |$0.96 \pm 0.24$| | 0.60 | 0.60 |$\pm$| 0.02 | |||
| F2 | MJD 56063–56173 | 110 | |$101/110$| | 1 | 3.0 | 1.08 | |$0.72 \pm 0.33$| | 0.05 | 0.46 |$\pm$| 0.03 |
| |$153/220$| | 0.5 | 4.0 | 0.71 | |$1.15 \pm 0.24$| | 0.70 | 0.37 |$\pm$| 0.03 | |||
| F3 | MJD 58363–58963 | 600 | |$517/600$| | 1 | 21.0 | 1.18 | |$1.36 \pm 0.24$| | 0.83 | 0.88 |$\pm$| 0.01 |
| |$941/1200$| | 0.5 | 21.0 | 0.64 | |$1.37 \pm 0.22$| | 0.54 | 0.84 |$\pm$| 0.01 | |||
| F4 | MJD 59063–59153 | 90 | |$75/90$| | 1 | 4.0 | 1.19 | |$0.88 \pm 0.47$| | 0.14 | 0.19 |$\pm$| 0.05 |
| |$105/180$| | 0.5 | 5.0 | 0.86 | |$0.50 \pm 0.27$| | 0.34 | 0.12 |$\pm$| 0.08 | |||
| F5 | MJD 59683–59943 | 260 | |$213/260$| | 1 | 22.0 | 1.22 | 1.27 |$\pm$| 0.37 | 0.72 | 0.53 |$\pm$| 0.02 |
| |$355/320$| | 0.5 | 21.5 | 0.73 | |$1.23\pm 0.20$| | 0.08 | 0.45 |$\pm$| 0.02 |
Note. (1) Flux states extracted from use of BB and HOP algorithm (2) Period of the flaring states (3) Total exposure time (4) Fraction of points having TS > 9 (5) Minimum sampling time in observed LC (6) Maximum sampling time in observed LC (7) Mean Sampling time i.e. total observation time over a number of data points in that interval (8) The power-law index for the power-law model of PSD analysis using PSRESP method (9) p-value corresponding to the power-law model. The power-law model is considered a bad fit if |$p_{\beta }\le 0.1$| as the rejection confidence for such model is |$\gt 90~{{\ \rm per\ cent}}$| (10) Fractional variability
The flaring periods, except for F4, exhibit similar |$\alpha$| parameters in HE Fermi–LAT spectrum, indicating a consistent physical process within a |$3\sigma$| range. Additionally, the |$\beta$| values for all periods are consistent, suggesting a similar influence of external seed photons in the production of high-energy |$\gamma$|-rays. The spectral parameters for the flaring periods F1-F5 are summarized in Table 2.
Fermi–LAT Spectral parameters for the chosen flaring states and the quiet state. The parameters are derived from the fitted log-parabola model.
| Flare state | |$\alpha \pm \delta \alpha$| | |$\beta \pm \delta \beta$| |
|---|---|---|
| F1 | 2.39 |$\pm$| 0.03 | 0.16 |$\pm$| 0.02 |
| F2 | 2.29 |$\pm$| 0.03 | 0.10 |$\pm$| 0.02 |
| F3 | 2.38 |$\pm$| 0.01 | 0.15 |$\pm$| 0.01 |
| F4 | 2.55 |$\pm$| 0.04 | 0.13 |$\pm$| 0.04 |
| F5 | 2.41 |$\pm$| 0.02 | 0.14 |$\pm$| 0.02 |
| Quiet state | 2.47 |$\pm$| 0.02 | 0.08 |$\pm$| 0.01 |
| Flare state | |$\alpha \pm \delta \alpha$| | |$\beta \pm \delta \beta$| |
|---|---|---|
| F1 | 2.39 |$\pm$| 0.03 | 0.16 |$\pm$| 0.02 |
| F2 | 2.29 |$\pm$| 0.03 | 0.10 |$\pm$| 0.02 |
| F3 | 2.38 |$\pm$| 0.01 | 0.15 |$\pm$| 0.01 |
| F4 | 2.55 |$\pm$| 0.04 | 0.13 |$\pm$| 0.04 |
| F5 | 2.41 |$\pm$| 0.02 | 0.14 |$\pm$| 0.02 |
| Quiet state | 2.47 |$\pm$| 0.02 | 0.08 |$\pm$| 0.01 |
Fermi–LAT Spectral parameters for the chosen flaring states and the quiet state. The parameters are derived from the fitted log-parabola model.
| Flare state | |$\alpha \pm \delta \alpha$| | |$\beta \pm \delta \beta$| |
|---|---|---|
| F1 | 2.39 |$\pm$| 0.03 | 0.16 |$\pm$| 0.02 |
| F2 | 2.29 |$\pm$| 0.03 | 0.10 |$\pm$| 0.02 |
| F3 | 2.38 |$\pm$| 0.01 | 0.15 |$\pm$| 0.01 |
| F4 | 2.55 |$\pm$| 0.04 | 0.13 |$\pm$| 0.04 |
| F5 | 2.41 |$\pm$| 0.02 | 0.14 |$\pm$| 0.02 |
| Quiet state | 2.47 |$\pm$| 0.02 | 0.08 |$\pm$| 0.01 |
| Flare state | |$\alpha \pm \delta \alpha$| | |$\beta \pm \delta \beta$| |
|---|---|---|
| F1 | 2.39 |$\pm$| 0.03 | 0.16 |$\pm$| 0.02 |
| F2 | 2.29 |$\pm$| 0.03 | 0.10 |$\pm$| 0.02 |
| F3 | 2.38 |$\pm$| 0.01 | 0.15 |$\pm$| 0.01 |
| F4 | 2.55 |$\pm$| 0.04 | 0.13 |$\pm$| 0.04 |
| F5 | 2.41 |$\pm$| 0.02 | 0.14 |$\pm$| 0.02 |
| Quiet state | 2.47 |$\pm$| 0.02 | 0.08 |$\pm$| 0.01 |
The time lag between these counterpart flares from the lensed images is estimated using the three methods described in Section 2.3. The maximum time delay between the lensed images is given by |${\sim}6 \left(\frac{z_g}{0.1} \right) \left(2\,h \right)^{-1}$| d, where |$z_g$| is the redshift of the lensing galaxy and h is the Hubble constant in units of 100 |$\rm {km\, s^{-1}\, Mpc^{-1}}$| (Zhang et al. 2008). Using a redshift of |$z_g=0.89$| and |$h=75\, \rm {km\, s^{-1}\, Mpc^{-1}}$|, the maximum time lag is approximately 71 d. This represents the maximum time delay between the mirage images when the source is near the Einstein ring. For large delays, the magnification ratio is expected to be significant. Thus, detecting a large delay is unlikely since the trailing component would be demagnified beyond the sensitivity of Fermi–LAT. Additionally, to explore the full range of potential time delays, the length of the lightcurve must be at least twice the longest time delay.
Detecting the trailing counterpart during a quiescent state of a blazar is challenging due to its expected demagnification, which impedes detection. However, the sensitivity of Fermi–LAT allows for the detection of multiple flaring states (F1 to F5). Dominant pink noise during flaring epochs creates spurious peaks, increasing the chances of detecting both trailing and leading components. The source shows the largest variability for F3 and F1, with respective fractional variability values (F|$_{var}$|; Vaughan et al. 2003) of |$0.88\pm 0.01$| and |$0.66\pm 0.02$| (Table 1). Our objective is to identify time delays from flares with sufficient magnification to detect both leading and trailing components. Whenever possible, pairs of leading and trailing flares, spaced within the identified time lag, are selected to study the spectral properties of source and echo flare and draw the relationship between lag and magnification. To analyse the flare variability properties, we fitted exponential flares to sharp distinct features in the light curve using the form
Here, |$\tau _{\rm {rise}}$| and |$\tau _{\rm {decay}}$| represent the rise and decay time-scales, respectively; |$t_o$| is the peak time, and |$F_0$| represents half of the peak flux of the flare at time |$t_o$|.
3.1 Flare 1 – MJD 55453–55583
The 1-d and 12-h binned light curve of flare F1 is shown in Fig. 3(a). The flaring period spans 150 d, with 18 per cent of the flux points not resulting in significant detection. The flare exhibits a sharp peak from MJD 55 483 to MJD 55 488 and another distinct feature from MJD 55 553 to MJD 55 573 in the |$\gamma$|-ray light curve.
(Top panel) 1-d binned and 12 h binned light curve of the flaring epoch F1 for FSRQ PKS 1830-211, covering the period from MJD 55450 to MJD 55600 in black and red respectively, as marked in Fig. 1. GPR predictions, computed using the parameters with the highest marginal likelihood, overlaid in blue on the 1-day binned data. (Middle panel) ACF on 1-d binned light curve of F1 period (Bottom panel) DPS on 1-d binned light curve.
The light curve exhibits prominent pink noise behaviour with a power-law index of |$k = 0.91 \pm 0.39$|. Monte Carlo simulations were conducted to quantify the significance of the time delay by generating light curves with similar power spectral density indices. The Autocorrelation function did not result in a significant detection, but a feature emerged at |$55 \pm 2$| d (|${\sim}2\sigma$|), likely an artefact of the Fermi light curve due to the spacecraft’s precession period of 53.4 d (Fig. 3b).
The DPS method is more prone to detecting spurious time delays. The DPS on a 1-d binned light curve peaked at |$17 \pm 1.5$| d with above |$3\sigma$| significance as shown in Fig. 3(c). Similarly, GPR analysis on the 1-d binned light curve identified the maximum marginal likelihood metric at |$19.0 \pm 1.5$| d. The corresponding best-fitting GPR light curve is shown in Fig. 8(a). The marginal likelihood peaking at |$9.8 \pm 2.9$| d could represent a lower harmonic of the |$19.0 \pm 1.5$| d delay, which aligns with the periodic nature of the selected kernel. Consistent results have been reported in previous studies, such as a lag of |$19 \pm 1$| d by Abdo et al. (2015) and |$17.9 \pm 7.1$| d by Barnacka et al. (2015) using ACF.
(Top panel) 1-d binned and 12 h binned light curve of the flaring epoch F2 for FSRQ PKS 1830-211, covering the period from MJD 56063 to MJD 56173 in black and red respectively, as marked in Fig. 1. GPR predictions, computed using the parameters with the highest marginal likelihood, overlaid in blue on the 1-day binned data. (Middle panel) ACF on 1-d binned light curve of F2 period (Bottom panel) DPS on 1-d binned light curve.
(Top panel) 1-d binned and 12 h binned light curve of the flaring epoch F3 for FSRQ PKS 1830-211, covering the period from MJD 58363 to MJD 58963 in black and red respectively, as marked in Fig. 1. GPR predictions, computed using the parameters with the highest marginal likelihood, overlaid in blue on the 1-day binned data. (Middle panel) ACF on 1-d binned light curve of F3 period (Bottom panel) DPS on 1-d binned light curve.
(Top panel) 1-d binned and 12 h binned light curve of the flaring epoch F4 for FSRQ PKS 1830-211, covering the period from MJD 59063 to MJD 59153 in black and red respectively, as marked in Fig. 1. GPR predictions, computed using the parameters with the highest marginal likelihood, overlaid in blue on the 1-day binned data. (Middle panel) ACF on 1-d binned light curve of F4 period (Bottom panel) DPS on 1-d binned light curve.
(Top panel) 1-d binned and 12 h binned light curve of the flaring epoch F5 for FSRQ PKS 1830-211, covering the period from MJD 59683 to MJD 59943 in black and red respectively, as marked in Fig. 1. GPR predictions, computed using the parameters with the highest marginal likelihood, overlaid in blue on the 1-day binned data. (Middle panel) ACF on 1-d binned light curve of F5 period (Bottom panel) DPS on 1-d binned light curve.
The likelihood metric for lags is derived using GPR. The bar represents the likelihood value for each lag, with the highest value indicating the most probable lag. The Gaussian fit over the likelihood values represents the mean lag value estimated and its corresponding error bar.
The maximum peak in the light curve at MJD 55 483 to MJD 55 488 and the resulting echo flare (F11 in Fig. 9a) after 20.0|$\pm$|1.1 d can be a proxy for evaluating the magnification ratio, resulting in a ratio of |$\sim$|2.8. Another distinct feature at MJD 55 553 to MJD 55 573 (F12), considering a similar time delay, shows a magnification ratio of |$\sim$|1.6 as in Fig. 9(a).
High flux states of Flare F1, F2, and zoomed section of F5 (MJD 59880–59940) and fitted exponential flare using equation (6). The vertical lines represent the source and echo pair for the lensed flares. The exponential fits with similar colours are considered possible pairs of source and echo flares.
3.2 Flare 2 – MJD 56063–56173
Fig. 4(a) shows the 1-d and 12-h binned light curves. This epoch is above the mean flux level, with periods before and after significantly below the mean. The 110-d-long light curve follows a pink noise power-law time series with an index of |$k=0.72 \pm 0.33$|. Simulated light curves with similar indices were generated to measure the significance of the lag in the signal.
The autocorrelation on F2 in Fig. 4(b) shows two features with significance close to |$2\sigma$|: |$11.0 \pm 2.3$| and |$20.3 \pm 0.5$| d, consistent with Barnacka et al. (2015). An additional lag close to |$2\sigma$| at |$55.7 \pm 2.2$| d appears in both the 1-d and 12-h binned light curves, possibly an artefact of the Fermi telescope’s processing period, similar to Flare F1.
The DPS method on the 1-d binned light curve detected a time delay of |$20 \pm 0.5$| d with more than |$2\sigma$| significance (Fig. 4c). Consistent lags are found using GPR, with increased marginal likelihood metrics at |$13.3 \pm 4.3$| and |$22.1 \pm 2.6$| d, aligning with the ACF results. The likelihood distribution of GPR on Flare F2 is as in Fig. 8(b). The |$13.3 \pm 4.3$| d lag in GPR could be a lower harmonic of the 22.1 d delay in the light curve.
Disentangling the light curve to identify flares and their echo flare image is challenging for Flare F2. The double-peak structure from MJD 56 081 to MJD 56 098 has a demagnified delayed counterpart appearing from MJD 56 101 to MJD 56118, with a demagnification of approximately 1.9 as shown in Fig. 9(b). No echo counterpart to the broad feature from MJD 56 142 to MJD 56 156 is visible within the 20-d period. This absence suggests a much larger demagnification corresponding to longer time delays.
3.3 Flare 3 – MJD 58363–58963
The 1-d and 12-h binned light curve for Flare 3 is shown in Fig. 5(a). This flare represents the brightest flux state of the source, with the flux reaching 14 times its average level. Notably, it is also the longest-lasting flare analysed in this work, spanning a total of 600 d. The flare exhibits the highest power spectral index, with characteristic noise behaviour between pink and red noise types. Multiple flares appear to be superimposed on an underlying envelope, as illustrated in Fig. 5(a).
The ACF, as displayed in Fig. 5(b), reveals two features with significance greater than |$3\sigma$|: one at |$12 \pm 1.8$| d and another at |$21.1 \pm 1.2$| d. For comparison, the DPS method calculates a prominent lag at |$21.0 \pm 0.5$| d (|${\gtrsim}3.5\sigma$|) and |$19.0 \pm 0.5$| d (|$=3\sigma$|). Less significant detections, but still above |${\gtrsim}2.5\sigma$|, are found at |$14 \pm 0.5$| and |$25 \pm 0.5$| d. This suggests multiple lag values imprinted on the flares. The estimated lag values could also be a reflection of the time difference between subsets of flares.
A consistent lag detected by GPR at |$14.0 \pm 2.2$| and |$21.1 \pm 1.2$| d is simultaneously observed by ACF and DPS (Fig. 8c).
Due to the dense overlap of flaring periods, identifying associated leading and trailing counterparts is extremely challenging during such a flaring period. Therefore, estimating magnification by flare identification is not performed for Flare 3.
3.4 Flare 4 – MJD 59063–59153
Fig. 6 shows the 1-d and 12-h binned light curve for Flare 4, which spans 90 d and is the least bright of the five states analysed. The power spectral density for the 1-d binned data yields |$k = 0.88 \pm 0.47$|. Ideally, the sample length should be twice the maximum expected time delay of 70 d, so the 90-d light curves reduce the likelihood of detecting time delays. Additionally, F4 has the least variability in light curve with |$F_{\mathrm{ var},F4}=0.19\pm 0.05$| d. Such small fractional variability highlights the absence of significant variable points in the light curve.
No significant time delay is detected by ACF and DPS in the 1-d binned data. GPR estimates a lag of |$22.4 \pm 2.2$| d as shown in Fig. 8(d). The absence of significant echo flares for the flux rise at MJD 59113–59123 (Fig. 6a) makes it difficult to conclude if a lensed image is detectable within the telescope’s sensitivity.
3.5 Flare 5 – MJD 59683–59943
Fig. 7(a) shows the 1-d and 12-h binned light curve for Flare 5, with multiple peaks visible over the 260-d period. At least 18.1 per cent of the data contains gaps or periods with significantly low detection, which we interpolated with zeros. No significant time lags were found using ACF. However, DPS estimates a |$17 \pm 0.5$| d time lag with more than |$2\sigma$| significance, consistent with the predicted time delay from GPR at |$19.4 \pm 2.7$| ds (See Fig. 8e).
The bright peaks at MJD 59 890 to MJD 59 928 (F51) were fitted with exponential function as shown in Fig. 9(c). However, due to the emergence of multiple overlapping flares, associating flares to identify magnification is challenging.
4 DISCUSSION
The time delay between two lensed images of a quasar, such as PKS 1830–211, is typically evident from systematic changes in flux density in the radio band for the leading and trailing components. However, in the |$\gamma$|-ray band, the resolution of high-energy telescopes constrains the detection of resolved images. As a result, the combined flux evolution from the two lensed images appears merged into a single point source.
The time delay observed from the merged images reflects the distribution of emission sites on the lens plane around the mass-weighted centre of the lens (Barnacka et al. 2014). Time delays, along with magnifications, are fundamental for localizing the emitting region relative to the black hole. Typically, a smaller time delay is associated with areas close to the base of the jet. Emissions originating near the base of the jets are typically reflected as the fast variability in high-energy |$\gamma$|-ray light curves.
In this work, we employed three methods to investigate the observed time lag in the stochastic time series. The resulting lags are listed in Table 3. We devised a novel technique, Gaussian Process Regression, to extract the lag in the signal. From our analysis, we found that the autocorrelation function is unable to efficiently extract the intrinsic time delay present in most signals, partially due to its noise dependence. However, among the five flares studied, autocorrelation successfully detected the intrinsic lag in the time series for Flare 3 with a significance of more than |$3\sigma$|. Flare 3 exhibited noise behaviour between pink and red, and the emergence of multiple flares resulted in significant detection. Similar results were obtained using the Double Power Spectrum and Gaussian Process Regression.
Estimated lags for Flares F1, F2, F3, F4, F5 using the three methods (a) Autocorrelation Function (b) Double Power Spectrum (c) Gaussian Process regression.
| Flare | |$\rm {Lag_{ACF}}$| | |$\rm {Lag_{DPS}}$| | |$\rm {Lag_{GPR}}$| |
|---|---|---|---|
| F1 | – | 17 |$\pm$| 1.5 (|${\gt}3\sigma$|) | 19.0 |$\pm$| 1.5 |
| F2 | 20.3 |$\pm$| 2.3 | 20.0 |$\pm$| 0.5 (|${\sim}2\sigma$|) | 22.1 |$\pm$| 2.6 |
| F3 | 20.5 |$\pm$| 1.0 | 21.0 |$\, \pm \,$| 0.5 (|${\sim}3\sigma$|) | 21.1 |$\pm$| 1.2 |
| F4 | – | 14.0 |$\pm$| 0.5 (|$\lt 2\sigma$|) | 22.4 |$\pm$| 2.2 |
| F5 | – | 17.0 |$\pm$| 0.5 (|${\gt}2\sigma$|) | 19.4 |$\pm$| 2.7 |
| Flare | |$\rm {Lag_{ACF}}$| | |$\rm {Lag_{DPS}}$| | |$\rm {Lag_{GPR}}$| |
|---|---|---|---|
| F1 | – | 17 |$\pm$| 1.5 (|${\gt}3\sigma$|) | 19.0 |$\pm$| 1.5 |
| F2 | 20.3 |$\pm$| 2.3 | 20.0 |$\pm$| 0.5 (|${\sim}2\sigma$|) | 22.1 |$\pm$| 2.6 |
| F3 | 20.5 |$\pm$| 1.0 | 21.0 |$\, \pm \,$| 0.5 (|${\sim}3\sigma$|) | 21.1 |$\pm$| 1.2 |
| F4 | – | 14.0 |$\pm$| 0.5 (|$\lt 2\sigma$|) | 22.4 |$\pm$| 2.2 |
| F5 | – | 17.0 |$\pm$| 0.5 (|${\gt}2\sigma$|) | 19.4 |$\pm$| 2.7 |
Estimated lags for Flares F1, F2, F3, F4, F5 using the three methods (a) Autocorrelation Function (b) Double Power Spectrum (c) Gaussian Process regression.
| Flare | |$\rm {Lag_{ACF}}$| | |$\rm {Lag_{DPS}}$| | |$\rm {Lag_{GPR}}$| |
|---|---|---|---|
| F1 | – | 17 |$\pm$| 1.5 (|${\gt}3\sigma$|) | 19.0 |$\pm$| 1.5 |
| F2 | 20.3 |$\pm$| 2.3 | 20.0 |$\pm$| 0.5 (|${\sim}2\sigma$|) | 22.1 |$\pm$| 2.6 |
| F3 | 20.5 |$\pm$| 1.0 | 21.0 |$\, \pm \,$| 0.5 (|${\sim}3\sigma$|) | 21.1 |$\pm$| 1.2 |
| F4 | – | 14.0 |$\pm$| 0.5 (|$\lt 2\sigma$|) | 22.4 |$\pm$| 2.2 |
| F5 | – | 17.0 |$\pm$| 0.5 (|${\gt}2\sigma$|) | 19.4 |$\pm$| 2.7 |
| Flare | |$\rm {Lag_{ACF}}$| | |$\rm {Lag_{DPS}}$| | |$\rm {Lag_{GPR}}$| |
|---|---|---|---|
| F1 | – | 17 |$\pm$| 1.5 (|${\gt}3\sigma$|) | 19.0 |$\pm$| 1.5 |
| F2 | 20.3 |$\pm$| 2.3 | 20.0 |$\pm$| 0.5 (|${\sim}2\sigma$|) | 22.1 |$\pm$| 2.6 |
| F3 | 20.5 |$\pm$| 1.0 | 21.0 |$\, \pm \,$| 0.5 (|${\sim}3\sigma$|) | 21.1 |$\pm$| 1.2 |
| F4 | – | 14.0 |$\pm$| 0.5 (|$\lt 2\sigma$|) | 22.4 |$\pm$| 2.2 |
| F5 | – | 17.0 |$\pm$| 0.5 (|${\gt}2\sigma$|) | 19.4 |$\pm$| 2.7 |
Our results suggest the presence of a consistent time delay of approximately 20 d during the flaring state of the source, as determined by the three methods used in this work. This indicates a similar orientation of the emitting site around the mass distribution of the lens. Consequently, |$\gamma$|-ray emission consistently occurs in similar regions of the jet across all flaring states. The inferred time delay aligns with the estimated lag reported in (Barnacka et al. 2015) during the high states, indicating that the origin of the |$\gamma$|-rays is likely within the core. The detection of rapid variability, with |$t_{var} \sim 0.38 \pm 0.22$| d, implies an emission region size of |$r_{emm} = c \delta t_{var}/(1+z) = 2.8 \times 10^{15}$| cm at a distance of |$R_{diss} = 2c\Gamma ^2 t_{var} = 0.064$| pc, thereby confirming that the high-energy emission is localized within the core on sub-parsec scales.
Lovell et al. (1996) reported a lag of |$26_{-5}^{+4}$| d using the Australian Telescope Compact Array at 8.6 GHz. Similar time delays were estimated in the |$\gamma$|-ray band during the quiet state of the source (Barnacka et al. 2011). The lag observed during high states in this study suggests a different origin for flaring |$\gamma$|-ray emission. The inconsistency between |$\gamma$|-ray and radio lags indicates different dissipation sites for these emissions, especially during the source’s high state. Radio emissions are typically expected from the outer regions of the parsec-scale jet. A small compact jet leads to synchrotron self-absorption, making radio emission unlikely in the inner parsec scales. Shorter time delays during active states imply that |$\gamma$|-ray dissipation occurs closer to the central engine, whereas radio dissipation occurs farther out in the jet. This has been observed as the absorption of high-energy photons with energies greater than 10 GeV during high states under the influence of BLR photons at sub-parsec scale jet (Agarwal et al. 2024).
The high-energy spectral properties of the source and the echo flare are consistent within 3|$\sigma$| (See Fig. 10). A change in the spectral properties would imply a difference in the influence of soft seed photons on |$\gamma -$|ray photon through |$\gamma -\gamma$| absorption on passing through a more luminous region of the lensing galaxy. A |${\sim}2.8\sigma$| deviation was observed in the spectral index for flare F21. The identical beta parameters for the four possible lensed flares suggest a similar influence of external seed photons from the local jet environment, the EBL, and the intervening galaxy on the high-energy spectrum. Since absorption affects all of the flares in the same way, they originated from similar regions of the jet. Our analysis focused on flares with a clear source and a demagnified lensed echo flare at an average flux level, leading us to select flares F1, F2, and F5. Due to the presence of multiple overlapping flares, identifying individual flares and their echoes for flare F3 was inefficient, likely due to the merging of multiple flares. Exponential fitting on individual flares reveals a linear relationship between lag and magnification. However, further studies are needed to identify more clean flares. This suggests that a smaller emission region is confined close to the base of the jet, while a larger magnification implies a larger emission region located farther out in the jet.
High energy spectral parameter of the flare and its associated echo flare. (left) |$\alpha$| and (right) |$\beta$| for the fitted log-parabola model.
5 SUMMARY
Strong gravitational lensing in |$\gamma$|-ray bright blazars can identify the locations of |$\gamma$|-ray dissipation during both quiescent and active states. The variation in time delays observed during periods of active |$\gamma$|-ray flux suggests different emission regions within the jet compared to those during low flux states (Barnacka et al. 2011). The consistent lag across five flaring states indicates a similar origin for the high-energy |$\gamma$|-ray activity within the radio core. This contrasts with the larger lag observed during quieter |$\gamma$|-ray periods and the consistent time delays from radio observations, which suggest that such emission occurs farther from the central engine than that during flaring periods. Such time delays caused by gravitational lensing of a background source by a foreground object could help constrain the Hubble parameter (Refsdal 1964).
We introduce a novel technique for estimating time delays in long continuous light curves from Fermi–LAT. Detecting these time delays could be crucial in identifying hidden lensed blazars during flaring periods that are not recognized as lensed sources in radio wavelengths. The signatures of such time delays could provide insights into distant blazars and previously unidentified |$\gamma$|-ray sources. Future surveys, including those conducted by the SKA, will likely discover many such lensed quasars. Extensive multiwavelength searches of these systems could offer valuable insights into the origins of radiation and provide a magnified view limited by current telescopes.
ACKNOWLEDGEMENTS
We thank the refree, Dr Nachiketa Chakraborty, for constructive feedback, which has helped improve the manuscript. SA and AS acknowledge Dr Bhargav Vaidya and Prof. Karl Mannheim for useful discussions on the work which helped improve the manuscript. AS acknowledges support for computational facility from DST–SERB grant CRG/2022/009332. This research work has made use of archival data, software, and web tools obtained from NASA’s High Energy Astrophysics Science Archive Research Center (HEASARC) and Fermi gamma-ray telescope Support centre, a service of the Goddard Space Flight Center and the Smithsonian Astrophysical Observatory.
DATA AVAILABILITY
The data will be made available upon reasonable request.
