Investigating Possible Correlations between Gamma-Ray and Optical Lightcurves for TeV-Detected Northern Blazars over 8 Years of Observations

Abstract

Blazars are a subclass of active galactic nuclei (AGN) having relativistic jets aligned within a few degrees of our line-of-sight and form the majority of the AGN detected in the TeV regime. The Fermi-Large Area Telescope (LAT) is a pair-conversion telescope, sensitive to photons having energies between 20 MeV and 2 TeV, and is capable of scanning the entire gamma-ray sky every three hours. Despite the remarkable success of the Fermi mission, many questions still remain unanswered, such as the site of gamma-ray production and the emission mechanisms involved. The Asteroid Terrestrial-impact Last Alert System (ATLAS) is a high cadence all sky survey system optimized to be efficient for finding potentially dangerous asteroids, as well as in tracking and searching for highly variable and transient sources, such as AGN. In this study, we investigate possible correlations between the Fermi-LAT observations in the 100 MeV–300 GeV energy band and the ATLAS optical data in the R-band, centered at 679 nm, for a sample of 18 TeV-detected northern blazars over 8 years of observations between 2015 and 2022. Under the assumption that the optical and gamma-ray flares are produced by the same outburst propagating down the jet, the strong correlations found for some sources suggest a single-zone leptonic model of emission.

1. Introduction

Blazars form a subclass of radio-loud active galactic nuclei (AGN) that have jets closely aligned to our line-of-sight, resulting in the emission from these objects being highly Doppler-boosted and making them some of the brightest gamma-ray sources in the extragalactic sky [1]. Blazars are generally characterized by non-thermal, highly-polarized continuum emission, spanning the entire electromagnetic spectrum and characteristically show very fast variability, which has been observed down to timescales of minutes in the gamma-ray regime [2,3,4], as well as in the optical regime [5,6,7].

The spectral energy distribution (SED) of a typical blazar comprises two distinct peaks. Although the first peak, occurring in the radio to the X-ray regime, has been attributed to synchrotron emission from electrons and positrons within the jet, the physical mechanisms responsible for the second peak, produced in the X-ray to gamma-ray regime, is still a matter of debate and two main scenarios have been postulated to explain it. Leptonic models [8,9] attribute the high-energy peak to the inverse Compton (IC) scattering between the energetic leptons in the jet and a field of low energy photons (for example, [10]), either the same photons emitted through synchrotron emission (synchrotron self-Compton (SSC)) or photon populations external to the jet (external inverse Compton (EIC)). On the other hand, hadronic models (for example, [11]) suggest that the second peak may be a result of high energy photons produced in cosmic-ray interactions through the decay of neutral pions (for example, [12]) or proton synchrotron emission (for example, [13]).

Localizing the gamma-ray emission is an indirect process and a variety of different methods have been used previously in the literature. Assuming constant jet geometry, the size of the emission region, r, has been used to infer its distance from the supermassive black hole (SMBH), R, using $\mathrm{r}=\psi \mathrm{R}$. Here $\psi$ is the semi-aperture opening angle of the jet [14,15]. This relation has been used to constrain the emission region to be within a few parsec from the base of the jet. For example, Ref. [16], using ∼2 years of Fermi-LAT observations, constrained the emission to be from within the broad line region (BLR) under the assumption that the full width of the jet is responsible for the emission.

Moreover, the observation of very high energy (VHE) photons ($E_{\gamma }\ge 20$ GeV) suggest the emission originates farther out, at parsec scale distances and from within the molecular torus (MT) region [17,18]. VHE photons are expected to be severely attenuated from interactions with the photons in the BLR and the detection of blazars with ground-based instruments is difficult to explain if the emission were assumed to originate in regions near the central engine. A possible solution which accommodates both the short variability timescales and VHE detection is to abandon the one-zone emission model and invoke the presence of multiple emission regions [19].

An important technique to localize the emission region in blazars is the correlation of multi-wavelength lightcurves. In particular, it helps in identifying potential relationships between emission zones and in understanding the dominant emission mechanisms. For instance, a strong correlation between synchrotron produced optical data and IC scattering produced gamma-ray observations (for example, [20,21]) is expected in some leptonic models. The lags and leads, if highly significant, can help to constrain the location of the gamma-ray emission region relative to the optical emission region and also discriminate between SSC and EC models. On the other hand, orphan flares in one waveband having no correlation with others, for example, the 2002 flare of 1ES 1959+650 [22], can be interpreted as evidence of multiple emission zones or as support for hadronic emission models.

In this respect, the Large Area Telescope (LAT) on board the Fermi satellite [23] has been particularly important. This pair-conversion telescope, launched in June 2008, is sensitive to photon energies between 20 MeV and 2 TeV and scans the entire gamma-ray sky every three hours. For example, Ref. [24] presented an investigation of correlations between optical data collected with the robotic 0.76 m Katzman Automatic Imaging Telescope (KAIT) at Lick Observatory and gamma-ray observations with the Fermi-LAT using discrete correlation functions (DCFs). The same method was used by [25,26] to study radio–gamma-ray correlations in blazars.

The Asteroid Terrestrial-impact Last Alert System (ATLAS [27]) is a high cadence all sky survey system comprising four independent units, one on Haleakala (HKO), and one on Mauna Loa (MLO) in the Hawaiian islands in the Northern Hemisphere and one each at the El Sauce Observatory, Chile and the South African Astronomical Observatory in the Southern Hemisphere. It is optimized to be efficient for finding potentially dangerous asteroids, as well as in tracking and searching for highly variable and transient sources, such as AGN, and is capable of discovering more bright, less than 19 mag, supernovae candidates than other ground based surveys.

Blazar observations with ATLAS are in the R-band, centered at 679 nm, having a typical cadence of one data point per two days. In this work, we apply local cross-correlation functions (LCCFs [28]) to investigate possible correlations between ATLAS optical data and gamma-ray observations with the Fermi-LAT for a sample of 18 TeV-detected northern blazars, over 8 years of observations between 2015 and 2022. For brighter sources (for example, Mrk 501 [29]), it is also possible to undertake an analysis of microvariability in outbursts using optical data and this will be investigated in detail in a future publication.

2. Source Selection and Data Reduction

The main goal of this study is to investigate possible correlations between ATLAS optical data and Fermi-LAT gamma-ray observations in blazars. The identification of suitable sources was primarily governed by having sufficient photon statistics to allow for a detailed study of the LCCFs. The sample of sources chosen for this study comprised 18 Fermi-LAT northern blazars detected in the TeV regime1 with ground-based gamma-ray observatories is shown in Table 1.

Optical observations were made with the ATLAS, a high cadence all sky survey system of four 0.5 m telescopes that scan the entire sky periodically. The telescopes are located in Hawaii (Mauna Loa, Hawaii and Haleakala, Maui), in South Africa (Sutherland Observing Station), and Chile (El Sauce Observatory, Rio Hurtado). Between declinations −50 degrees and + 50 degrees, the cadence is one set of observations each day, and in the polar regions once every two days, during observing season and weather permitting. The filter used in these observations is roughly equivalent to the Johnson–Cousins R filter. The transmission curve is centered at 679 nm. The filter is nonstandard, but well calibrated, see [27]. All the data is archived in a database for retrieval.

A forced photometry system is used in which the point spread function (PSF), which represents the distribution of light from a point source on a detector, is obtained for each object based on bright stars, and the function is forced onto the object at the chosen coordinates. Data processing is described in [27,30]. When the data is retrieved from the database, software automatically calculates the magnitude according to the AB system, as well as the flux for each object. A catalog of variable stars produced by this method is shown by [31]. The optical lightcurves for the sample of blazars are shown in Figure 1, Figure 2 and Figure 3.

In this work, we analyzed Fermi-LAT photons detected between MJD 57000 and MJD 59945, which corresponds to midnight on 9 December 2014 until midnight on 1 January 2023. Throughout the gamma-ray analysis, we used the Fermi Science Tools version 11-05-032, FERMIPY version 1.0.13 [32], in conjunction with the latest PASS 8 IRFs [33]. We consider the energy range 100 MeV–300 GeV and a region of interest (RoI) with radius 15${}^{\circ }$ centred on each source. A lower limit of 100 MeV is chosen because the PSF of the Fermi-LAT increases at lower energies, making a point source analysis difficult. Furthermore, most of the blazars selected have relatively soft gamma-ray spectra and are not expected to be significantly detected with the Fermi-LAT at energies above 300 GeV. Moreover, we selected only photon events from within a maximum zenith angle of 90${}^{\circ }$ in order to reduce contamination from background photons from the Earth’s limb, produced by the interaction of cosmic-rays with the upper atmosphere.

Sources in the 4FGL-DR3 catalog [34] within a radius of ${20}^{\circ }$ from the spatial position of each source in the sample were included in the initial model for the analyses, with their spectral parameters fixed to their catalog values. This takes into account the gamma-ray emission from sources lying outside the RoI which might yet contribute photons to the data, especially at low energies, due to the size of the point spread function of the Fermi-LAT. The contributions from the isotropic and galactic diffuse backgrounds were modeled using the most recent templates, iso_P8R3_SOURCE_V2_v1.txt and gll_iem_v07.fits, respectively.

The analysis began with an initial automatic optimization of the RoI by iteratively fitting the sources. This ensured all parameters were close to their global likelihood maxima. The spectral normalization of all modelled sources within the RoI were left free as were the normalization factor of both the isotropic and galactic diffuse emission templates. Furthermore, the spectral shape parameters of all sources within 5${}^{\circ }$ of the centre of the RoI were left free to vary whereas those of other sources were fixed to the values reported in the 4FGL-DR3 catalog. A binned likelihood analysis was then performed in order to obtain the spectral parameters best describing the model using a spatial binning of 0.1${}^{\circ }$ pixel${}^{-1}$ and 8 energy bins per decade.

In order to pursue a study of the temporal behaviour of the gamma-ray fluxes, the 15 year Fermi-LAT data were binned monthly with a likelihood routine applied to each bin separately4. During the production of the lightcurve, the spectral parameters of all sources within 5${}^{\circ }$ of the RoI centre were left free for each bin as were the normalization factors of the background emission models. The resulting lightcurves for the sample of blazars are shown in Figure 1, Figure 2 and Figure 3, along with the corresponding uncertainties. Only time intervals having TS ≥ 10 were considered, which roughly equates to a significance of $3\sigma$.

3. Results

Local cross-correlation functions (LCCFs; [28]) were then applied to investigate correlations between the Fermi-LAT and ATLAS lightcurves. Consider two lightcurves having fluxes $a_i$ and $b_j$, corresponding to times $t_{a_i}$ and $t_{b_j}$, the LCCF can then be computed as:

$$\mathrm{LCCF}\left(\tau \right)=\frac{1}{M}\frac{\mathsf{\Sigma }(a_i-{\overline{a}}_{\tau })(b_j-{\overline{b}}_{\tau })}{{\sigma }_{a\tau }{\sigma }_{b\tau }}.$$

(1)

Here, the sums run over M pairs for which $\tau \le t_{a_i}-t_{b_j}<\tau +\mathsf{\Delta }t$ for a chosen timestep $\mathsf{\Delta }t$, and ${\overline{a}}_{\tau }$ and ${\overline{b}}_{\tau }$ are flux averages and ${\sigma }_{a\tau }$ and ${\sigma }_{b\tau }$ are standard deviations over the M pairs, respectively [28]. As adopted in [35], we choose the higher half median of the time separations between consecutive data points in the two lightcurves as the binning of the timelags, $\tau$. Furthermore, the minimum and maximum values of $\tau$ are chosen to be $\pm 0.5$ times the length of the shorter lightcurve [36]. This method is independent of any difference in sampling rates between the lightcurves. LCCFs are intrinsically bound in the interval [−1,1] and [36] found them to be more efficient than the use of Discrete Correlation Functions (DCFs; [37]). The LCCFs obtained for each source are shown in Figure 4 and Figure 5.

As done in [24], the centroid lags and uncertainties are derived from weighted least-square Gaussian fits to the LCCF points at the location of the most significant correlation peak. Although the peaks of the Gaussian fits give a first order determination of the uncertainty, these do not account for the effects of correlated red-noise between the datasets [38]. The significances of the correlations are obtained by performing Monte Carlo simulations to produce 1000 artificial lightcurves matching the probability distribution function (PDF) and power spectral density (PSD) of each observation using the method outlined in [39]5. The 68%, 95%, and 99% confidence intervals obtained are also shown in Figure 4, in blue from lighter to darker shades.

The centroid lags and significances obtained for each source are summarised in Table 1. Throughout this paper, a positive time delay, $\mathsf{\Delta }t>0$, is defined as corresponding to the optical emission leading the gamma-ray emission. Wider LCCF peaks may indicate a range of characteristic timescales in the correlated response or limitations in the instruments [24]. The presence of a peak, significant at a level of above 3$\sigma$, implies strong optical and gamma-ray correlations. This is seen for W Comae, S3 1227+25, 3C 66A, and 1ES 1215+303, and suggests a single-zone model of emission. Under the assumption that the optical and gamma-ray flares are produced by the same outburst propagating down the jet, both positive and negative time-lags on timescales of days to tens of days are predicted in SSC and EIC models. A significant peak consistent with zero indicates an absence of time-lag and is seen for S3 1227+25 and 3C 66A.

The results obtained in this investigation are broadly in agreement with those seen in [20], where the optical and gamma-ray correlations of 121 blazars were found to be significant at a level of above $68\%$ (∼$1\sigma$). The majority of the corresponding time-lags were within $\pm 20$ days of zero-lag and this was interpreted as evidence of a common origin for the flares in the two bands, implying leptonic processes dominate blazar gamma-rays, while not excluding the possibility of hadronic contribution to the emission. However, it should be noted that although [20] found that the time delays, in general, did not exceed 50 days at $2\sigma$ level and was less at $3\sigma$ level, in this study, we find significant delays for 11 of the 18 considered objects, often exceeding 100 days.

Furthermore, out of the nine objects with a correlation of above $99\%$, five sources, namely, 1ES 0414+009, 1ES 1959+650, 1ES 1440+122, 1ES 1011+496, 1ES 1218+304, and PKS 1424+240, are found have a delay of more than 100 days. In all of these cases a negative time-lag is found, suggesting that the gamma-ray emission lags the optical emission. This long delay could imply a more complex emission mechanism than the simple one-zone leptonic SSC model found for the remaining sources; for example, one with multiple emission regions with the optical emission originating upstream of the gamma-ray emission (this was also seen for gamma-ray and radio correlations in [26]). However, continuous and simultaneous monitoring over a longer time period is required to draw any stronger conclusions.

Moreover, while similar conclusions of leptonic single-zone models where the low- and high-energy emission comes from the same population of electrons were also found in [21] from a study of gamma-ray and optical correlations for a sample of 1180 blazars, it should also be noted that they find evidence of “orphan” flares, gamma-ray flares with no optical counterpart or optical flares having no gamma-ray counterpart, in some of the sources. As seen in Figure 1, Figure 2 and Figure 3, a few sources in this investigation, for example, 1ES 0502+675 (MJD∼59750) and S3 1227+25 (MJD∼58600), also show “orphan” flares. The origin of these flares is uncertain but they have been interpreted as support for hadronic emission models (for example, [40]), evidence of multiple emission zones (for example, [19]) or a result of contamination in the optical regime due to accretion disk emission (for example, [41]).

Table 1. The final sample of Fermi-LAT northern blazars detected in the TeV regime, selected for this study along with their classification, right ascensions (RA), and declinations (DEC) in degrees, redshifts (z; shown as “–” for sources not having a verified redshift measurement). The final two columns list the times corresponding to the peaks of the Gaussian fit along with the associated uncertainties and their significance in percentile derived from Monte Carlo simulations. Listed below are the TeV detection references for each source in the sample.

Source	Class	RA	Dec	z	Timelag	Significance
		[deg]	[deg]		[days]	[%]
1ES 0502+675 ${}^{\left[a\right]}$	HBL	76.9	67.6	0.340	−9.75 ± 3.52	≥68
W Comae ${}^{\left[b\right]}$	IBL	185.4	28.2	0.102	74.78 ± 6.85	≥99
VER J0521+211 ${}^{\left[c\right]}$	IBL	80.4	21.2	–	−236.07 ± 6.65	≥95
S3 1227+25 ${}^{\left[d\right]}$	IBL	187.6	25.3	–	20.71 ± 7.10	≥99
1ES 0414+009 ${}^{\left[e\right]}$	HBL	64.2	1.1	0.287	−746.48 ± 5.34	≥99
1ES 0806+524 ${}^{\left[f\right]}$	HBL	122.5	52.3	0.138	626.02 ± 6.65	≥95
RXJ0648.7+1516 ${}^{\left[g\right]}$	HBL	102.2	15.3	0.179	466.59 ± 3.85	≥95
3C 66A ${}^{\left[h\right]}$	IBL	35.7	43.0	0.340	2.61 ± 2.67	≥99
1ES 1215+303 ${}^{\left[i\right]}$	HBL	184.5	30.1	0.131	−56.10 ± 5.97	≥99
RGB J0710+591 ${}^{\left[j\right]}$	EHBL	107.6	59.1	0.125	−130.00 ± 28.90	≥68
1ES 1959+650 ${}^{\left[k\right]}$	HBL	299.9	65.2	0.048	−683.44 ± 9.26	≥99
1ES 0229+200 ${}^{\left[l\right]}$	EHBL	38.2	20.3	0.139	331.57 ± 19.98	≥68
1ES 1440+122 ${}^{\left[m\right]}$	HBL	220.7	12.0	–	−272.09 ± 2.03	≥99
Mrk 421 ${}^{\left[n\right]}$	HBL	166.1	38.2	0.031	28.68 ± 7.08	≥68
1ES 1218+304 ${}^{\left[o\right]}$	EHBL	185.3	30.2	0.182	−565.78 ± 2.09	≥99
1ES 1011+496 ${}^{\left[p\right]}$	HBL	153.8	49.4	0.212	−122.25 ± 8.26	≥99
Mrk 501 ${}^{\left[q\right]}$	HBL	253.5	39.8	0.034	−574.50 ± 9.65	≥68
PKS 1424+240 ${}^{\left[r\right]}$	HBL	216.8	23.8	–	−581.28 ± 1.78	≥99

TeV detection references: [a] VERITAS [42]; [b] VERITAS [43]; [c] VERITAS [44]; [d] VERITAS [45]; [e] H.E.S.S. [46]; [f] VERITAS [47]; [g] VERITAS [48]; [h] VERITAS [49]; [i] MAGIC [50]; [j] VERITAS [51]; [k] Utah [52]; [l] H.E.S.S. [53]; [m] VERITAS [54]; [n] Whipple [55]; [o] MAGIC [56]; [p] MAGIC [57]; [q] Whipple [58] and [r] MAGIC [59].

Table 1. The final sample of Fermi-LAT northern blazars detected in the TeV regime, selected for this study along with their classification, right ascensions (RA), and declinations (DEC) in degrees, redshifts (z; shown as “–” for sources not having a verified redshift measurement). The final two columns list the times corresponding to the peaks of the Gaussian fit along with the associated uncertainties and their significance in percentile derived from Monte Carlo simulations. Listed below are the TeV detection references for each source in the sample.

Source	Class	RA	Dec	z	Timelag	Significance
		[deg]	[deg]		[days]	[%]
1ES 0502+675 ${}^{\left[a\right]}$	HBL	76.9	67.6	0.340	−9.75 ± 3.52	≥68
W Comae ${}^{\left[b\right]}$	IBL	185.4	28.2	0.102	74.78 ± 6.85	≥99
VER J0521+211 ${}^{\left[c\right]}$	IBL	80.4	21.2	–	−236.07 ± 6.65	≥95
S3 1227+25 ${}^{\left[d\right]}$	IBL	187.6	25.3	–	20.71 ± 7.10	≥99
1ES 0414+009 ${}^{\left[e\right]}$	HBL	64.2	1.1	0.287	−746.48 ± 5.34	≥99
1ES 0806+524 ${}^{\left[f\right]}$	HBL	122.5	52.3	0.138	626.02 ± 6.65	≥95
RXJ0648.7+1516 ${}^{\left[g\right]}$	HBL	102.2	15.3	0.179	466.59 ± 3.85	≥95
3C 66A ${}^{\left[h\right]}$	IBL	35.7	43.0	0.340	2.61 ± 2.67	≥99
1ES 1215+303 ${}^{\left[i\right]}$	HBL	184.5	30.1	0.131	−56.10 ± 5.97	≥99
RGB J0710+591 ${}^{\left[j\right]}$	EHBL	107.6	59.1	0.125	−130.00 ± 28.90	≥68
1ES 1959+650 ${}^{\left[k\right]}$	HBL	299.9	65.2	0.048	−683.44 ± 9.26	≥99
1ES 0229+200 ${}^{\left[l\right]}$	EHBL	38.2	20.3	0.139	331.57 ± 19.98	≥68
1ES 1440+122 ${}^{\left[m\right]}$	HBL	220.7	12.0	–	−272.09 ± 2.03	≥99
Mrk 421 ${}^{\left[n\right]}$	HBL	166.1	38.2	0.031	28.68 ± 7.08	≥68
1ES 1218+304 ${}^{\left[o\right]}$	EHBL	185.3	30.2	0.182	−565.78 ± 2.09	≥99
1ES 1011+496 ${}^{\left[p\right]}$	HBL	153.8	49.4	0.212	−122.25 ± 8.26	≥99
Mrk 501 ${}^{\left[q\right]}$	HBL	253.5	39.8	0.034	−574.50 ± 9.65	≥68
PKS 1424+240 ${}^{\left[r\right]}$	HBL	216.8	23.8	–	−581.28 ± 1.78	≥99

TeV detection references: [a] VERITAS [42]; [b] VERITAS [43]; [c] VERITAS [44]; [d] VERITAS [45]; [e] H.E.S.S. [46]; [f] VERITAS [47]; [g] VERITAS [48]; [h] VERITAS [49]; [i] MAGIC [50]; [j] VERITAS [51]; [k] Utah [52]; [l] H.E.S.S. [53]; [m] VERITAS [54]; [n] Whipple [55]; [o] MAGIC [56]; [p] MAGIC [57]; [q] Whipple [58] and [r] MAGIC [59].

4. Conclusions

In this work, we present an investigation into possible correlations between ATLAS optical data and gamma-ray observations with the Fermi-LAT for a sample of 18 TeV-detected northern blazars, over 8 years of observations between 2015 and 2022. Overall, we find the optical and gamma-ray emission to be highly correlated in our sample, with varied time delays, ranging from timescales of days to even years for some sources. With the exception of one source, 1ES 1218+304, all the correlations are found to correspond to a significance of $68\%$ (∼$1\sigma$) and for nine sources the correlations are found to correspond to significance of $99\%$ (∼$3\sigma$).

The observed strong correlations support leptonic models of IC scattering gamma-ray emission in which the seed photons are scattered to higher energy in the relativistic jet by the same electrons responsible for synchrotron emission. It should be noted that the significance of the correlations and the corresponding time delays do not allow us to make strong conclusions on whether the seed photons are dominated by SSC or EIC radiation. However, the lack of clear trend towards lag or lead in our BL Lac sample agrees with the results presented in [24] and can be interpreted as evidence towards SSC being the dominant mechanism in BL Lac sources, as opposed to FSRQs which, in general, show the gamma-ray leading the optical emission, interpreted as evidence towards EIC emission being the dominant mechanism.

In conclusion, the gamma-ray optical correlation in BL Lac sources appears complex, as also seen in other variability studies (for example [21,24]). Multi-wavelength studies at high cadence over many years is needed to probe the emission mechanisms further and this will be possible with further coverage with ATLAS and other optical instruments in conjunction with the continued successful operation of the Fermi-LAT. Finally, we aim to perform a comprehensive study of the sample further investigating micro-variability in both energy bands in the future.