Optical intraday variability of the blazar S5 0716+714 Open Access

ABSTRACT

We present an extensive recent multiband optical photometric observations of the blazar S5 0716+714 carried out over 53 nights with two telescopes in India, two in Bulgaria, one in Serbia, and one in Egypt during 2019 November – 2022 December. We collected 1401, 689, 14726, and 165 photometric image frames in B, V, R, and I bands, respectively. We monitored the blazar quasi-simultaneously during three nights in B, V, R, and I bands; four nights in B, V, and R; two nights in V, R, and I; five nights in B and R; and two nights in V and R bands. We also took 37 nights of data only in R band. Single band data are used to study intraday flux variability and two or more bands quasi-simultaneous observations allow us to search for colour variation in the source. We employ the power-enhanced F-test and the nested ANOVA test to search for genuine flux and colour variations in the light curves of the blazar on intraday time-scales. Out of 12, 11, 53, and 5 nights observations, intraday variations with amplitudes between ∼3 and ∼20 per cent are detected in 9, 8, 31 and 3 nights in B, V, R, and I bands, respectively, corresponding to duty cycles of 75, 73, 58, and 60 per cent. These duty cycles are lower than those typically measured at earlier times. On these time-scales colour variations with both bluer-when-brighter and redder-when-brighter are seen, though nights with no measurable colour variation are also present. We briefly discuss possible explanations for this observed intraday variability.

1 INTRODUCTION

The centre of active galactic nuclei (AGNs) host supermassive black holes (SMBHs) of masses in the range of |${\sim } 10^{6}\!-\!{10}^{10}\, \rm {M}_{\odot }$| (Rees 1984) that accrete matter from their environments and have been observed to emit both thermal and non-thermal radiations. A small fraction of these AGNs emit strongly at radio wavelengths and the ratio of 5 GHz radio to optical B-band emission, called the radio-loudness parameter R, is one of the common traditional ways to characterize these sources as well as a measure of non-thermal emission. If R < 10 then the AGN is called radio-quiet (RQ) and if R ≥ 10 the AGN is known as radio-loud (RL). The latter hosts powerful large-scale bipolar relativistic jets – one of the primary sources of non-thermal emission. A large fraction, ∼ 85–90 per cent belongs to the RQAGNs class and the rest ∼10–15 per cent are RLAGNs (e.g. Kellermann et al. 1989). The orientation of the jet is found to be tightly correlated with the dominance of non-thermal component and broad-band emission. This has led to a new designation scheme for these sources, namely, jetted or non-jetted, in which RQAGNs lack significant emission above hard X-rays while RLAGNs have strong emission up to GeV–TeV gamma-rays (e.g. Padovani 2017).

Blazars are one of the most enigmatic subclasses of RLAGNs and they emit radiation throughout the entire electromagnetic (EM) spectrum from radio to very high energy γ-rays (VHEs; |$\rm E \gtrsim 100$| GeV), with their relativistic plasma jets oriented almost along the line of sight of the observer (e.g. Urry & Padovani 1995). The blazar category is a concatenation of BL Lacertae objects (BL Lacs) and flat spectrum radio quasars (FSRQs). This separation is based on the strength of emission lines in optical-UV bands, measured in terms of equivalent width (EW), with BL Lacs exhibiting a featureless continuum or very weak emission lines (EW |$\le \rm {5}\, \mathring{\rm A}$|⁠) (Stocke et al. 1991; Marcha et al. 1996), while FSRQs show stronger emission lines (Blandford & Rees 1978; Ghisellini et al. 1997). In the complete EM spectrum, the emission from a blazar is predominantly non-thermal primarily due to the Doppler-boosted jet emission. Blazars display double humped spectral energy distributions (SEDs) across the entire EM spectrum (Fossati et al. 1998). Based on the location of the hump peaks, blazars historically have been classified into low-energy peaked blazars (LBLs) and high-energy peaked blazars (HBLs). The first hump peaks in the infrared or optical bands for LBLs and in UV or X-ray for HBLs. The high-energy hump, respectively, peaks in GeV and TeV γ-ray ranges for LBLs and HBLs (Fossati et al. 1998).

Over the last few decades it has been well established that blazars show flux, spectral, and polarization variations in all accessible EM bands on diverse time-scales ranging as short as a few minutes to as long as several years (e.g. Miller, Carini & Goodrich 1989; Urry et al. 1993; Heidt & Wagner 1996; Raiteri et al. 2001, 2006, 2008, 2011; Gupta et al. 2012, 2016, 2017, 2022; Gaur et al. 2014; Hayashida et al. 2015; Bhatta et al. 2016; Itoh et al. 2016; Kiehlmann et al. 2016; Kushwaha 2022, and references therein). Blazar variability on a few minutes to less than a days is variously called as micro-variability (Miller, Carini & Goodrich 1989) or intraday variability (IDV; Wagner & Witzel 1995) or intranight variability (Gopal-Krishna, Sagar & Wiita 1993); variability time-scales from a few days to several weeks are known as short-term variability, while the variability over months to decades is called long-term variability (Gupta et al. 2004). The nature of blazar variability in whole EM spectrum is mostly non-linear, stochastic, and aperiodic (e.g. Kushwaha et al. 2017).

The blazar S5 0716+714¹ (z = 0.31 ± 0.08) (Nilsson et al. 2008) is one of the brightest BL Lacs in optical bands. The discovery of very high energy (VHE) γ-ray emission from S5 0716+714 was reported in MAGIC (Major Atmospheric Gamma Imaging Cherenkov) observations (Anderhub et al. 2009), and it is listed in the catalogue of TeV emitting sources.² Wagner & Witzel (1995) reported that the duty cycle of the object in optical bands is one, which indicates that the source is almost always in an active state displaying variability. S5 0716+714 is therefore among the most well-studied blazars for variability studies across the EM bands on diverse time-scales (e.g. Wagner et al. 1990; Quirrenbach et al. 1991; Wagner et al. 1996; Ghisellini et al. 1997; Raiteri et al. 2003; Tagliaferri et al. 2003; Bach et al. 2005; Pian et al. 2005; Foschini et al. 2006; Villata et al. 2008; Gupta et al. 2012; Bhatta et al. 2013; Rani et al. 2013; Chandra et al. 2015; Wierzcholska & Siejkowski 2016; Sandrinelli et al. 2017; MAGIC Collaboration et al. 2018, and references therein). Over the last three decades, the source has received a great deal of attention in terms of searches for optical IDV (e.g. Heidt & Wagner 1996; Villata et al. 2000; Wu et al. 2005; Montagni et al. 2006; Poon, Fan & Fu 2009; Agarwal et al. 2016; Raiteri et al. 2021, and references therein). Using the major optical outbursts in the source, a possible long-term optical period was estimated to be 3.0 ± 0.3 yr (Gupta et al. 2008). On IDV time-scales, periodic or quasi-periodic oscillations were apparently detected on a few occasions (Gupta, Srivastava & Wiita 2009; Rani et al. 2010). Blazar emission in high-flux states, i.e. pre/post outburst states and outburst states, as well as detected IDVs during those states, are non-thermal Doppler-boosted emission from jets (Blandford & Rees 1978; Marscher & Gear 1985; Hughes, Aller & Aller 1992; Marscher, Gear & Travis 1992). But when the blazar is in very low-flux states, optical IDV may be explained by disturbances or hotspots on the accretion disc surrounding the central SMBH (e.g. Mangalam & Wiita 1993). An alternative interpretation based on low-luminosity AGN behaviour argues that the accretion discs in low-luminosity AGNs are radiatively inefficient, so a weak jet emission could still be responsible for detection of optical IDV in the low-flux state of blazars (Chiaberge et al. 2006; Capetti et al. 2007).

The flux and spectral variabilities seen in blazars on IDV time-scales are some of the most puzzling issues in AGN physics (e.g. Czerny et al. 2008). To address the IDVs in blazars, we started a long-term project in the mid-2000s and have made extensive searches in optical and X-ray bands using various ground and space-based telescopes observations (e.g. Gupta et al. 2008, 2012; Gaur et al. 2010, 2015; Kalita et al. 2015; Agarwal et al. 2016; Pandey, Gupta & Wiita 2017; Aggrawal et al. 2018; Zhang et al. 2019, 2021; Dhiman et al. 2021; Devanand et al. 2022; Pavana Gowtami et al. 2022, and references therein). Here, we present an extended study of optical flux and spectral variabilities of the blazar S5 0716 + 714 on IDV time-scales using observations from two telescopes in India, two in Bulgaria, and one each in Serbia and Egypt. We have carried out these optical photometric observations of the blazar S5 0716+714 during a recent span of ∼3 yr (2019 November – 2022 December) that provide additional useful data on this object. These multiband optical observations also allow us to explore the spectral evolution of the source over a long-term scale during its diverse flux states.

The paper is organized as follows. In the next section, we describe the details of data gathering and the reduction procedure. Section 3 describes the various analysis tests we employ. Results including our inferences and discussion are in Section 4. We end with the summary of this study in Section 5.

2 OBSERVATION AND DATA REDUCTION

Optical photometric observations of the blazar S5 0716+714 have been made with six ground-based optical telescopes between 2019 February and 2022 December, using Johnson–Cousin BVRI filters. The largest number of these observations were carried out by two telescopes at the Aryabhatta Research Institute of Observational Sciences (ARIES), Nainital, India, namely the 1.04 m Sampuranand Telescope (ST), and the 1.3 m Devasthal Fast Optical Telescope (DFOT). The next largest number of observations were obtained at the 1.88 m telescope at Kottamia Astronomical Observatory (KAO), Egypt. Additional data were taken by the 2.0m and 60cm telescopes at National Astronomical Observatory (NAO), Rozhen, Bulgaria, and Astronomical Observatory (AO), Belogradchik, Bulgaria, respectively. Finally, a few nights of multiband data were collected at the 1.4 m telescope at Astronomical Station Vidojevica (ASV), Serbia. Detailed information about these telescopes and their CCD (charged coupled device) detectors are provided in Table 1. Table 2 provides the complete observation log.

The raw data from the 1.3 m DFOT was pre-processed using the standard routines of iraf (Image Reduction and Analysis Facility)³ software. Pre-processing includes bias subtraction, flat-fielding, edge trimming, and cosmic-ray removal. The Dominion Astronomical Observatory Photometry (daophot II) software (Stetson 1987, 1992) was then used on the 1.3m DFOT data to obtain the instrumental magnitudes of the blazar S5 0716+714 and the standard stars in the blazar field. Finally, aperture photometry was done to find the instrumental magnitude of the blazar and standard stars in the same field of view. Each image frame was subjected to aperture photometry using four different concentric circular aperture of radii: 1 × full width at half-maximum (FWHM), 2 × FWHM, 3 × FWHM, and 4 × FWHM. However, we found that the aperture radius 2 x FWHM has the best signal-to-noise ratio (S/N; Gupta et al. 2016; Pandey et al. 2019); therefore, we have taken it for our final results. In every observation, at least three local standard stars (Ghisellini et al. 1997; Villata et al. 1998) were also observed in the same frame from which one standard star was used for calibrating the blazar to find the apparent magnitude of the blazar S5 0716+714, and other two standard stars were used to check their mutual non-variability. Since the Blazar S5 0716+714 and standard stars were in the same frame, there was no need to correct for atmospheric extinction.

The data from the 1.04m ST and 1.88m KAO telescopes were reduced using the python language, utilizing various modules: astropy (Astropy Collaboration 2022), photutils (Bradley et al. 2022), ccdproc (Craig et al. 2017), and astrometry.net (Lang et al. 2010). First, the raw data were calibrated for instrumental biases by bias and flat correction using the ccdproc module. Then it was corrected for cosmic-rays using a technique called lacosmic (van Dokkum 2001). Finally, it was astrometrically calibrated using astrometry.net, where the instrumental magnitude of the blazar and comparison stars were measured within an aperture of twice the FWHM using photutils.

The data from the 2m RC NAO and 60-cm AO Bulgarian Telescope were reduced using ESO-MIDAS2 (Gaur et al. 2012a), the European Southern Observatory Munich Image Data Analysis System, which is developed and maintained by European Southern Observatory. The 1.4m ASV Serbian Telescope data were reduced using MAXIMDL⁴ (Pandey et al. 2020). Data from every night has been calibrated using comparison star 5⁵ because it was the closest in magnitude and colour to the blazar.

3 ANALYSIS TECHNIQUES

To search for intraday variability in the blazar S5 0716 + 714, we have used two statistical tests, namely the power-enhanced F-test and nested analysis of variance (ANOVA) test (de Diego 2014; de Diego et al. 2015). By using these methods, we examined the differential light curves (LCs) of the blazar for intraday variations, as they have been argued to be more reliable and powerful than other widely used statistical tests such as the C-test and the standard F-test (de Diego et al. 2015). In both of these tests, we use several standard stars as comparison stars in the blazar to ascertain the presence of even small amplitude variability in the blazar. We have also found the percentage of amplitude change in the magnitude and performed cross-correlation between different bands in cases where we collected sufficient multiband data for the same observing night.

3.1 Power enhanced F-test

In the power-enhanced F-test, we use comparison star 5 as a reference star to find the differential LCs of the blazar S5 0716 + 714 because it was closest to our object in magnitude and colour. We compare the variance of the blazar LC to a combined variance of comparison stars. It is defined as (Pandey et al. 2019)

where |$s^2_{blz}$| is the variance of the differential LCs of the difference of instrumental magnitude between the blazar and reference star and,

is the variance of the combined differential LCs of the instrumental magnitudes of the comparison star and reference star. N_j is the number of data points of the jth comparison star and |$s^2_{j,i}$| is its scaled square deviation, which is defined as

Here, ω_j, m_j,i, and |$\bar{m}_j$| are the scaling factor of the jth comparison star’s differential LC, its differential magnitude, and its mean magnitude, respectively. The averaged square error of the blazar differential LC divided by the averaged square error of the jth comparison star is used as the scaling factor ω_j.

In this work, we have three comparison stars: C4, C5, and C6. As C5 was closest in magnitude and colour to the blazar, we have taken C5 as the reference star. Since two more comparison stars are left, k = 2. The blazar and all the comparison stars have the same number of observations N, so the number of degrees of freedom in the numerator is (N − 1) and in the denominator is k(N − 1). We have found the F_enh using equation (1) and compared it with a criterion F_c at the confidence level of 99 per cent, i.e. α = 0.01. If F_enh > F_c, then the differential LC is considered as variable (Var); otherwise, it is considered as non-variable (NV).

3.2 Nested ANOVA test

In the nested ANOVA test, we use all the comparison stars as reference stars to find the differential LCs. So unlike power enhanced F-test, which pulled out one reference star, we have one more star to work with in the nested-ANOVA test (de Diego et al. 2015). We use comparison stars C4, C5, and C6 to generate the differential LCs of the blazar. We split these differential LCs so that we have 5 points in each group. The drawback of this technique is that it cannot identify microvariations within each group of observations that are shorter than the span of the grouping. This is unlikely to cause problems because only a few reports of such very short ‘spike’ events have been reported in the literature (Sagar, Gopal-Krishna & Wiita 1996; de Diego et al. 1998; Stalin et al. 2004).

From equation 4 of de Diego et al. (2015), we have estimated the value of |$\rm {SS}_G$| (sum of squares due to groups) and |$\rm {SS}_{O(G)}$| (sum of squares due to nested observations in groups). These sums of squares (SS) divided by the respective degrees of freedom ν, give the mean square errors (MS = SS/ν). We then estimated the F-statistics using the ratio |$F = \rm {MS}_G/\rm {MS}_{O(G)}$| which follows the F-distribution with a degree of freedom of (a − 1) for the numerator and a(b − 1) for the denominator, where a is the number of groups in the night’s observations and b is the number of data points in each group. At a confidence level of 99 per cent, i.e. α = 0.01, if F > F_C then we say that differential LCs are variable (V) otherwise we say it is non-variable (NV). The results of the power enhanced F-test and nested-ANOVA tests are given in Tables 3 and A1.

3.3 Intraday variability amplitude

We compute the IDV amplitude (Amp) in all of the calibrated LCs which were found to be variable by using the relation given by Heidt & Wagner (1996):

Here, A_max & A_min are, respectively, the maximum and minimum calibrated magnitudes of the blazar and σ is the mean error. The variability amplitudes we found are reported in Tables 3 and A1.

3.4 Duty cycle

Variability has been one of the primary defining/characterizing criterion of blazars as well AGNs. It can be quantified in term of duty cycle (DC) which measures the fraction of duration for which the source remains variable/active. We have estimated the DC of the blazar S5 0716+714 by using the standard approach of Romero, Cellone & Combi (1999):

The variable N_i takes the value 1 if IDV is detected but 0 if it is not. The computation of the DC is weighted by the observing time Δt_i for the ith observation, which is normally different for each observation. Here, we use the redshift-corrected observing time, denoted by Δt_i = Δt_i,obs(1 + z)⁻¹.

3.5 Discrete correlation function

The discrete correlation function (DCF) is a statistical tool used to quantify the correlation and identify possible time lags between different data series. It was introduced by Edelson & Krolik (1988) and generalized to provide better estimates of the errors by Hufnagel & Bregman (1992). The DCF is calculated by comparing pairs of data points at different time lags, and is particularly useful for unevenly sampled data, as it does not require interpolation in the temporal domain. To calculate the DCF, the unbinned correlation function (UDCF) is first determined for each pair of data points (x_i, y_j) in the two data series and is given by the equation

where |$\bar{x}$| and |$\bar{y}$| are the mean values of the two discrete data series x_i and y_j, with standard deviations σ_x and σ_y and measurement errors e_x, e_y. The UDCF is then averaged over a range of time delays, and the resulting value is the DCF for that time delay. It is given by the equation

where τ is the centre of the bin of size Δτ and N is the number of data points in the temporal bin width, Δτ. The DCF measures the degree of correlation between the two data series, with values greater than zero indicating a positive correlation, values less than zero indicating a negative correlation, and values equal to zero indicating no correlation. The error is identified using the standard deviations of the number of bins used to calculate the DCF, and it is written as

3.6 Halving/doubling time-scale

Several distinct approaches have been used to quantify different variability time-scales. A standard method for doing this is the halving/doubling time-scale that describes how rapid a nominal flux rise or decay of the source is. Treating flux rise and decay equally, for a pair of data points in the LC, it is given as

where F_t is the flux at some epoch t, t₀ is the reference epoch with flux |$F_{t_0}$|⁠, and τ is the halving/doubling time-scale which further takes the form

For our analysis, we used the BVRI magnitudes and evaluated τ for every possible pair of data points in each of the nights of observation when the source was found to be variable. In this analysis, we employed the selection criterion from Foschini et al. (2011), accepting only those pairs of data points for which the magnitude difference exceeded three times the associated uncertainty. The use of these data point pairs results to a set of estimates of τ and we quote the minimum value. The uncertainties corresponding to each τ estimate were evaluated using standard error propagation techniques. Defining the ratio, |$(m_t-m_{t_0}) /(t-t_0) = q$|⁠, the uncertainties in τ are given as

in terms of the uncertainties in the measurements of the magnitudes, |$\delta _{m_t}$| and |$\delta _{m_{t_0}}$|⁠.

Tables 3 and A1 include the shortest of the halving/doubling timescales and their errors for those nights during which the source was variable.

4 RESULTS AND DISCUSSION

We have carried out multiband optical photometric observations of the blazar S5 0716 + 714 in a total 53 observing nights from 2019 November 10 to 2022 December 22 using 6 optical telescopes listed in Table 1. One or more of these telescopes observed the blazar quasi-simultaneously in BVRI bands in 3 nights, BVR bands in 4 nights, VRI bands in 2 nights, BR bands in 5 nights, VR bands in 2 nights, and only in the R band on 37 nights. In total 17 081 source image frames were observed during the whole observing programme of which 1401, 712, 14662, and 306 image frames were in B, V, R, and I bands, respectively. This has allowed us to obtain extensive information on the IDV variability of this blazar.

4.1 Intraday flux variability

The calibrated optical B-, V-, R-, and I-band magnitudes of the blazar S5 0716 + 714 are plotted with blue, green, red, and black colours, respectively, in Fig. 1 for a sample the of LCs, and in Fig. A1 for all of the LCs. Observation dates and filter names are mentioned in each panel. If LCs from more than one filter are plotted, offsets are applied for clarity and the offest values are mentioned in the panel. To search for the genuine presence of IDV in LCs, we performed statistical tests described in Sections 3.1 and 3.2, and results for a sample of LCs are listed in Table 3 and for all LCs in Table A1. In the status columns of Tables 3 and A1, NV and Var denote non-variable and variable LCs. We estimated amplitude and variability time-scales of variable IDV LCs described in Sections 3.3 and 3.4, respectively, and these are also reported in Tables 3 and A1.

There are 12, 11, 53, and 5 intraday LCs in B, V, R, and I bands, respectively. Genuine IDV is detected in 9/12 in B, 8/11 in V, 31/53 in R, and 3/5 in I-band LCs, respectively. We estimated intraday duty cycles as described in Dhiman et al. (2023), and found them to be about 75, 73, 58, and 60 per cent in B, V, R, and I bands, respectively. Wagner & Witzel (1995) mentioned that the blazar S5 0716 + 714 was highly variable in early optical observations. The source has been extensively studied to search for optical IDV over the last 3 decades and by using all these intraday LCs, the DC was found to be 84, 80, 82, and 85 per cent in B, V, R, and I bands, respectively (Heidt & Wagner 1996; Wagner et al. 1996; Sagar et al. 1999; Villata et al. 2000; Massaro, Montagni & Nesci 2001; Nesci, Massaro & Montagni 2002; Montagni et al. 2006; Stalin et al. 2006; Gupta et al. 2008, 2012; Poon, Fan & Fu 2009; Carini, Walters & Hopper 2011; Wu et al. 2012; Gaur et al. 2012b; Bhatta et al. 2013; Hu et al. 2014; Agarwal et al. 2016; Hong, Xiong & Bai 2017; Xu et al. 2019, and references therein). Our estimated DCs during this recent epoch appear to be somewhat less than the DCs found in this combined sample of cumulative IDV LCs of S4 0716+714 studied over that more extended period. This may indicate either a recent slight decline in activity or merely arise from the differences inherent in these studies. In several of the earlier studies no rigorous test that involves checking putative variability of a blazar has been performed, whereas we have utilized the conservative nested ANOVA and enhanced F-tests with a confidence level of 0.99. Therefore, it is to be expected that our values might be somewhat lower. Other reasons why our DC results might be lower are the smaller number of LCs in the single present work and the somewhat shorter typical length of our nightly measurements.

During the entirety of our observations, IDV amplitudes are found in the range of ∼ 3–20 per cent which is consistent with the previous studies (e.g. Poon, Fan & Fu 2009; Gaur et al. 2012b; Agarwal et al. 2016, and references therein), and IDV time-scales in the range of 0.1–2.8 h. In ∼ 11.4 h of B-band observation on 2020 December 15, a time-scale of 9.5 h is found. In several nights, we made quasi-simultaneous observations in more than one band. On 2019 December 28, 2020 October 22, 2020 December 17, and 2021 March 7, we found the blazar’s variability amplitude is larger at higher frequencies. Such a trend apparently suggests that the source spectrum gets flatter with increasing brightness, and steeper with decreasing brightness (e.g. Massaro et al. 1998; Agarwal & Gupta 2015). But on several other occasions e.g. 21 October 2020; 12 February 2021, 30 November 2021, 28 February 2022, and 13 March 2022, we noticed that the variability amplitude at lower frequencies was comparable or larger than that at higher frequencies, as has also been seen in the past (e.g. Ghosh et al. 2000; Gaur et al. 2015).

4.2 Intraday cross-correlated flux variability

Out of 53 nights of total observations, there are 16 nights in which quasi-simultaneous observations were carried out in 2–4 optical bands. These observations provide us with an excellent opportunity to search for time lags between bands on IDV timescales when the blazar shows genuine variability in two or more bands. Of those 16 nights, S5 0716 + 714 exhibited IDV on 12 nights. We have adopted the DCF analysis technique described in Section 3.4 to find the cross-correlations and to estimate any time lag, if present between bands. We consider the genuine detection of lag if it is >3 times of the bin size of DCF analysis. Using this criterion, we found that in all the 12 nights with variations detected in with multiple bands, the DCF peaks at about zero lag, as shown in Fig. 2 for the sample of V- and R-band variable LCs. All the DCF plots are provided in Fig. A2.

The strong correlated variability with zero time lag suggests that the emission in different optical bands on each of the observing nights are produced by the same physical processes and arise from the same region. The result is expected as these optical frequencies fall in a very narrow window of the complete EM spectrum that belong to the same synchrotron spectral component which powers the optical-UV emission. Similar results were found earlier for a few other extensively studied blazars (e.g. Poon, Fan & Fu 2009; Bachev et al. 2011; Wu et al. 2012; Agarwal & Gupta 2015; Agarwal et al. 2016; Bachev et al. 2017, 2023; Dhiman et al. 2023, and references therein). We note that there have been a few occasions when a time lag in the range of 6–13 min was reported between two optical bands in the blazar S5 0716 + 714 (Qian, Tao & Fan 2000; Villata et al. 2000; Stalin et al. 2006).

4.3 Intraday colour variability

We have studied the colour variation with respect to R-band for all the 12 nights that showed multiband IDV, as R-band data were obtained on all of them. First, we found differences in the calibrated magnitudes of two bands, e.g. (B − I), (V − I), (B − R), and (V − R), then plotted them against R band. Since the number of data points is not the same in all the bands, we have binned the data for the bands for which we have a larger number of data points. To study these behaviours we have fitted the trends with a straight line, y = mx + c where y is a Colour Index and x is R band. The fitting is done using the python language module curvefit. The results are given in Table 4 and plots are displayed in Fig. 3.

There, respectively, are 2, 3, 9, and 8 colour–magnitude plots in (B − I), (V − I), (B − R), and (V − R) with respect to R band. Positive and negative slopes m indicate that the trend is bluer-when-brighter (BWB) and redder-when-brighter (RWB), respectively. We only claim the presence of a genuine colour variation when the slope |$m \gt \rm {3}\sigma$|⁠, the correlation coefficient |r| > 0.5, and the null hypothesis probability p < 0.01 (i.e. 99 per cent) (Gaur et al. 2012b; Dhiman et al. 2023). In a total of five observing nights, we detected both BWB and RWB trends in colour–magnitude plots of the blazar S5 0716+714 on IDV timescales (see Fig. 3 and Table 4). Similar BWB and RWB results were previously reported in earlier observations of the source (Hong, Xiong & Bai 2017; Wang, Xiong & Bai 2019). A BWB trend is commonly seen in BL Lacs, but occasionally a RWB trend is also seen (e.g. Gaur et al. 2012b; Dhiman et al. 2023, and references therein). When a BWB trend is seen in optical bands it indicates that essentially all that emission is synchrotron from the jet, but the simplest explanation for a RWB trend is that the accretion disc is also contributing to the optical emission (Villata et al. 2006; Raiteri et al. 2007; Gaur et al. 2012b).

We estimated the average spectral index, <α_BR>, for all the nights during which quasi-simultaneous B and R band optical observations were carried out. We applied the following expression given by Wierzcholska et al. (2015) for average spectral index α_BR estimates for individual nights

where ν_B = 6.874 × 10¹⁴ Hz and ν_R = 4.679 × 10¹⁴ Hz represent the effective frequencies of the B and R bands, respectively (Bessell 2005). We display the spectral indices as a function of time in Fig. 4. We find no appreciable long-term change in the spectral index over time.

4.4 Relation between average flux and IDV time-scale

Our sample of 31 genuinely variable IDV LCs of S5 0716 + 704 measured in the optical R band are obtained from observations carried out over more than 3 yr (2019 November 10 to 2022 December 22). Since this blazar has apparently shown long-term quasi-periodicity on the timescale of 3.0 ± 0.3 yr (Raiteri et al. 2003; Gupta et al. 2008), it has probably been observed in all of its normal flux states (e.g. low, intermediate, and high) during our observations. In Fig. 5, we plot the average nightly flux verses IDV time-scale. We see that there is no consistent relationship between the blazar’s IDV timescale and average nightly flux.

A shock travelling down the inhomogeneous medium in the jet is thought to be the cause of the largest variations in the optical emission from blazars and may dominate the overall flux in intermediate and high brightness states (e.g. Marscher & Gear 1985; Marscher, Gear & Travis 1992; Wagner & Witzel 1995). In the lowest brightness states of blazars, particularly FSRQs, the detected optical IDV may arise from instabilities in the accretion disc (e.g. Mangalam & Wiita 1993). But for a BL Lac such as S5 0716+714, observable disc contributions are very unlikely. Small increases in the Doppler factor, δ, of the relativistic jet can make for both substantial observed increased fluxes (∼∝δ³) and a reduction in the measured timescales, as δ⁻¹. Hence, a correlation between higher average fluxes and shorter IDV time-scales is expected if the IDV arises solely from a single region which undergoes rapid changes in the bulk jet velocity and/or viewing angle. The absence of such an obvious correlation in the average flux and IDV time-scales suggests that more than one zone is involved, at least most of the time. A similar result was found in this blazar in another optical study based on data taken between 1999 November 26 and 2003 March 23 (Gupta, Srivastava & Wiita 2009).

4.5 Possible origins of IDV

Flux variability on diverse time-scales is one of the basic features of blazars and has been one of the characteristics used to classify them as peculiar sources when they were discovered. The inability to spatially resolve the hypothesized AGN structure except perhaps for some of the very nearest AGNs (e.g. M87, Event Horizon Telescope Collaboration et al. 2019) has made variability an indispensable tool with which to infer spatial scales. Significant IDV variability implies a highly compact emitting region within these extra-galactic objects.

Since the Doppler-boosted non-thermal radiation from the relativistic jet typically outweighs the thermal radiation from the accretion disc in blazars, models based on the relativistic jet are most likely to account for the variability on any observable time scale. The detected optical IDV in the blazar S5 0716 + 714 is of intrinsic origin. The dominant fundamental emission model for intrinsic variability over longer timescales is shocks propagating through the jet (Blandford & Königl 1979; Marscher & Gear 1985). For such models to produce intrinsic IDV, a relativistic shock propagating down a jet while interacting with flow irregularities (Marscher, Gear & Travis 1992) or relativistic shocks changing directions due to jet precession (Nesci et al. 2005) are possible. Non-axisymmetric bubbles are carried outward in relativistic magnetized jets in a different emission scenario (Camenzind & Krockenberger 1992). The IDV detected in the low-flux state may be explained by hotspots on, or instabilities in, the accretion discs (Chakrabarti & Wiita 1993; Mangalam & Wiita 1993).

The various IDV behaviours seen in the optical bands of blazars may be caused by the presence of high magnetic fields within the relativistic jet (Romero, Cellone & Combi 1999). The strength of the magnetic fields present in the jet of HBLs has been suggested as the cause of the variations in IDV behaviour. An axial magnetic field (B) greater than a critical value (B_c) may stop the jet base’s Kelvin–Helmholtz instabilities and bends from occurring, which would otherwise result in IDV. This suggests that the magnetic field in the jet of S5 0716 + 714 is weaker than B_c which is given by Romero (1995) as

where n_e, m_e, and γ are the local electron density, the rest mass of electron, and the bulk Lorentz factor of the flow, respectively. The Doppler factor is given by δ = [γ(1 − βcos θ)]⁻¹, where βc is the velocity of emitting plasma, and θ is the jet viewing angle. For S5 0716 + 714, |$\theta , \ \delta, \ \rm {and} \ \gamma$| are reported to be 3.0 ± 0.4 deg, 15.6 ± 4.0, and 14.0 ± 3.7, respectively (Jorstad et al. 2017). The electron density n_e of S5 0716 + 714 is found in the range of 10²–10⁴ electrons cm⁻³ (Liao et al. 2014; MAGIC Collaboration et al. 2018). Considering those values of these parameters, and within the framework of this model for damping of jet instabilities, we obtain estimates for B in the range of 0.07–0.70 G. Using multiwavelength SEDs of S5 0716 + 714, B has been argued to be in the range of 0.01–0.58 G (Liao et al. 2014; MAGIC Collaboration et al. 2018). We note that our estimated range of B values is consistent with those earlier estimates.

The presence of turbulence in the jet will generate stochastic variations in the synchrotron emission. A shock encountering a turbulent region in the jet (e.g. Marscher, Gear & Travis 1992) is one of the possible explanations for various features in the IDV LCs of the blazar S5 0716 + 714. A shock would accelerate particles in each cell, which then cool by synchrotron emission, increasing the flux. There will be variations in the density, size, and magnetic field direction of each individual turbulent cell. After the shock travels through the turbulent cell, the radiation intensity decreases.

One other way to produce a substantial fast change is through the production of a minijet, where a small portion of the jet is accelerated to ultrarelativistic velocities, most likely through magnetic reconnection (e.g. Giannios, Uzdensky & Begelman 2009), though this mechanism may be more important for making γ-ray flares than IDV. Turbulence in the jet in the vicinity of a shock can provide the fast changes in Doppler factors for a small region that could yield IDV as well as flux variations on longer timescales (Marscher 2014; Calafut & Wiita 2015; Pollack, Pauls & Wiita 2016). More recently, it has been shown that turbulence in magnetized relativistic jets can be generated through the development of kink instabilities (e.g. Acharya, Borse & Vaidya 2021; Kadowaki et al. 2021; Medina-Torrejón et al. 2021). The instabilities and turbulence then drive magnetic reconnection and particle acceleration. The turbulence in both kinetic and magnetic energies eventually approach three-dimensional Kolmogorov spectra that produce fast flux variations (Kadowaki et al. 2021).

5 SUMMARY

We have presented an extensive multiband optical observations of the blazar S5 0716 + 714 taken from 2019 November to 2022 December using six telescopes in India, Bulgaria, Serbia, and Egypt. A summary of our results are as follows:

• S5 0716 + 714 showed frequent and significant IDV in optical fluxes. The duty cycles in B, V, R, and I bands are 75, 73, 58, and 60 per cent, respectively. The maximum variability amplitude is found to be ∼20 per cent.

• The cross-correlated variability in different bands on same night of observation show that the variations are simultaneous within the limits of the cadence of the measurements.

• Colour variations are present rather frequently and both BWB and RWB trends are seen in the colour–magnitude plots. On three observing nights, the LCs could be divided into two or three segments, such that the different segments show different colour variability trends.

• We found there is no clear trend between average nightly flux and the shortest IDV time-scale measured that night. This indicates that more than one emission region may normally be producing these fast fluctuations, as expected in turbulent relativistic jet models.

ACKNOWLEDGEMENTS

We thank the anonymous reviewer for useful comments which helped us to improve the manuscript. TT would like to acknowledge financial support from the Department of Science and Technology, Government of India, through INSPIRE-fellowship grant no. DST/INSPIRE Fellowship/2019/IF190034. ACG is partially supported by a Chinese Academy of Sciences (CAS) President’s International Fellowship Initiative (PIFI) (grant no. 2016VMB073). NRIAG team acknowledges support from the Egyptian Science, Technology & Innovation Funding Authority (STDF) under grant number 45779. This research was partially supported by the Bulgarian National Science Fund of the Ministry of Education and Science under grants KP-06-H38/4 (2019), KP-06-KITAJ/2 (2020), and KP-06-H68/4 (2022). Financial support from the Bulgarian Academy of Sciences (Bilateral grant agreement between BAS and SANU) is gratefully acknowledged. GD and OV acknowledge support by the Astronomical station Vidojevica, funding from the Ministry of Science, Technological Development and Innovation of the Republic of Serbia (contract no. 451-03-47/2023-01/200002), by the EC through project BELISSIMA (call FP7-REGPOT-2010-5, no. 265772), the observing and financial grant support from the Institute of Astronomy and Rozhen NAO BAS through the bilateral SANU-BAN joint research project GAIA ASTROMETRY AND FAST VARIABLE ASTRONOMICAL OBJECTS, and support by the SANU project F-187. PK acknowledges financial support from the Department of Science and Technology (DST), Government of India, through the DST-INSPIRE faculty grant (DST/INSPIRE/04/2020/002586). HG acknowledges financial support from the Department of Science and Technology (DST), Government of India, through INSPIRE faculty award IFA17-PH197 at ARIES, Nainital, India. MFG acknowledges support from the National Science Foundation of China (grant no. 11873073), Shanghai Pilot Program for Basic Research Chinese Academy of Science, Shanghai Branch (JCYJ-SHFY2021-013), the National SKA Program of China (grant no. 2022SKA0120102), the Original Innovation Program of the Chinese Academy of Sciences (E085021002), and the science research grants from the China Manned Space Project with No. CMSCSST-2021-A06. ZZ is thankful for support from the National Natural Science Foundation of China (grant no. 12233005).

DATA AVAILABILITY

The data in this article will be shared after one year of the publication at the reasonable request to the first author.

Footnotes

References

APPENDIX: RESULTS OF IDV