Multi-band optical variability on diverse timescales of the TeV blazar TXS 0506+56, the first cosmic neutrino source

We report the first extensive optical flux and spectral variability study of the TeV blazar TXS 0506+056 on intra-night to long-term timescales using BVRI data collected over 220 nights between January 21, 2017 to April 9, 2022 using 8 optical ground-based telescopes. In our search for intraday variability (IDV), we have employed two statistical analysis techniques, the nested ANOVA test and the power enhanced F-test. We found the source was variable in 8 nights out of 35 in the R-band and in 2 of 14 in the V-band yielding Duty Cycles (DC) of 22.8% and 14.3%, respectively. Clear colour variation in V - R was seen in only 1 out of 14 observing nights, but no IDV was found in the more limited B, I, and B - I data. During our monitoring period the source showed a 1.18 mag variation in the R-band and similar variations are clearly seen at all optical wavelengths. We extracted the optical (BVRI) SEDs of the blazar for 44 nights when observations were carried out in all four of those wavebands. The mean spectral index (\alpha) was determined to be 0.897+-0.171


INTRODUCTION
Blazars are Active Galactic Nuclei (AGN) characterized by highly collimated relativistic jets with half opening angles 5 • closely aligned with the observer's line of sight (e.g., Urry & Padovani 1995).The Doppler enhanced intense non-thermal radiation from this jet dominates the spectral energy distribution (SED) from radio to very high energy (VHE) γ-ray energies.Blazars are usually considered to be comprised of BL Lacertae (BL Lac) objects which show featureless spectra or very weak emission lines (equivalent width EW 5 Å) (Stocke et al. 1991;Marcha et al. 1996), and flat spectrum radio quasars (FSRQs) with prominent emission lines in their composite optical/UV spectra (Blandford & Rees 1978;Ghisellini et al. 1997).Blazars show flux and spectral variability across the entire electromagnetic (EM) spectrum, emit predominantly non-thermal radiation showing strong polarization (> 3%) from radio to optical frequencies, and usually have core dominated radio structures.
⋆ Email: dhiman@aries.res.in† Email: acgupta30@gmail.com The multi-wavelength (MW) SEDs of blazars in log(νFν) vs. log(ν) representation show double-humped structures in which the low energy hump peaks in infrared (IR) through X-ray bands while the high energy hump peaks at γ-ray energies (Fossati et al. 1998).The location of SED peaks are often used to classify blazars into two sub-classes, namely LBLs (low-energy-peaked blazars) and HBLs (high-energy-peaked blazars).In LBLs, the first hump peaks in IR to optical bands and the second in GeV γ-ray energies (Padovani & Giommi 1995).Whereas in HBLs, the first hump peaks in UV to X-ray bands and the second in up to TeV γ-ray energies (Padovani & Giommi 1995).The emission of the lower energy SED hump is due to synchrotron radiation which originates from relativistic electrons in the jet.The high-energy hump can be attributed to two mechanisms.One of these is inverse Compton (IC) scattering of low-energy photons by the same electrons responsible for the synchrotron emission (synchrotron-self Compton, SSC) or external photons (external Compton, EC), collectively known as the leptonic model (e.g., Böttcher 2007).The other mechanism is emission from relativistic protons or muon synchrotron radiation, referred to as the hadronic model (e.g., Mücke et al. 2003).
Flux variability over a wide range of timescales is one of the definitional properties of blazars.On the basis of the timescales over which it is observed, blazar variability can be divided into three somewhat arbitrary categories: microvariability (Miller et al. 1989) or intra-day variability (IDV) (Wagner & Witzel 1995) or intra-night variability (Gopal-Krishna et al. 1993) (occurring on a timescale of a few minutes to less than a day); short-term variability (STV; taking place on a timescale of days to months); and long-term variability (LTV; over a timescale of several months to years or even decades) (Gupta et al. 2004).
Colour-magnitude diagrams (CMDs) for blazars can be analyzed to find any colour trends accompanying brightness changes.Three types of CMD behaviour could be discerned: redder-when-brighter (RWB), bluer-when-brighter (BWB), and achromatic.FSRQs mostly show RWB chromatism because the contribution of the accretion disc to the total emission is significant (Gu et al. 2006;Gaur et al. 2012c).The BWB behaviour seen in many BL Lacs is thought to arise from processes occurring in the relativistic jet, such as particle acceleration and cooling in the framework of the shock-in-jet model (e.g., Marscher & Gear 1985;Kirk et al. 1998).Alternatively, the BWB chromatism could arise from a Doppler factor variation on a convex spectrum (e.g., Villata et al. 2004b;Papadakis et al. 2007).Finally, achromatic behaviour is frequently interpreted as being due to the variations of the Doppler factor, which are most likely explained in a framework involving changes in the viewing angle to the dominant emission region (e.g., Villata et al. 2002;Gu et al. 2006;Gaur et al. 2012c;Agarwal & Gupta 2015;Agarwal et al. 2016).TXS 0506+056 is registered as a blazar in the Texas Survey of Radio Sources catalog (Douglas et al. 1996).The first detection of a high-energy neutrino event from a blazar was reported from TXS 0506+056 on 22 September 2017 by the IceCube collaboration and was coincident in direction and time with a γ−ray flare (IceCube Collaboration et al. 2018a).Prompted by this discovery, an investigation was carried out of 9.5 years of IceCube neutrino observations to search for excess emission at the position of the blazar, and an excess of high-energy neutrino events between September 2014 and March 2015 at energies around 290 TeV at a 3.5 σ level was indeed detected (IceCube Collaboration et al. 2018b).This object is the highest energy γ-ray-emitting blazar among those detected by the Energetic Gamma Ray Experiment Telescope (EGRET) satellite in the γ-ray energy range (30 MeV−30 GeV) (Dingus & Bertsch 2001).It is identified as a jet-dominated in the low-hard state during neutrino flaring in 2014/2015, and so provides evidence for the blazar jet acting as an accelerator of cosmic-ray particles which produce neutrinos (Padovani et al. 2018).TXS 0506+056 was independently detected at high energy γ-rays with the Large Area Telescope (LAT) onboard the Fermi satellite (Tanaka et al. 2017), the MAGIC telescope (Ansoldi et al. 2018), and the AGILE γ-ray telescope (Lucarelli et al. 2019), which strengthens the case for TXS 0506+056 being a very high energy (VHE) γ-rays emitting BL Lac as well as a neutrino emitting source.However, Padovani et al. (2019) claimed that TXS 0506+056 is a masquerading BL Lac, i.e., a FSRQ with hidden broad lines and a standard accretion disc that is outshined by the jet emission.During the intensive follow-up observations, the redshift of TXS 0506+056 was successfully determined to be z = 0.3365 (Paiano et al. 2018).The analysis of single-dish 15 GHz radio flux densities from the Owens Valley Radio Observatory (OVRO) spanning between 2008 and 2018 indicates that the core of TXS 0506+056 is in a highly flaring state coincident with the neutrino event 170922A (Britzen et al. 2019;Kun et al. 2019).In this paper we are reporting an extensive optical variability study of the first neutrino emitting TeV blazar TXS 0506+056 on diverse timescales.This paper is organised as follows: Section 2 provides an overview of the telescopes, photometric observations, and the data reduction procedure.Analysis techniques we used to search for flux variability and correlations between bands are discussed in Section 3. Results of our study are reported in Section 4. A discussion and conclusions are provided in Section 5.

OBSERVATIONS AND DATA REDUCTION
Optical photometric observations of the TeV blazar TXS 0506+056 were carried out using 8 ground-based telescopes.Two telescopes are located in India: the 1.04 m Sampurnanand Telescope (ST), and the 1.3 m Devasthal Fast Optical Telescope (DFOT), ARIES, Nainital.Both of these telescopes are equipped with CCD detectors and broadband Johnson UBV and Cousins RI filters.The source was observed with alternate observations in the V and R bands for a total of 37 nights between November 7, 2019 and January 31, 2021.One or two B and I image frames were also taken on each night of observations.
Observations of this source with the 60 cm Cassegrain telescope located at the Astronomical Observatory Belogradchik, Bulgaria, were carried out over the course of 40 nights from 10 October 2018 to 17 August 2020, consisting on a single optical frame in B, V , R, and I bands each night.These 40 nights of observations were presented in Bachev et al. (2021).The 2-m Ritchey-Chretien telescope at the National Astronomical Observatory Rozhen, Bulgaria, observed only a single night in V , R, and I bands on 17 August 2020.Observations of the source with the 1.3 m modified RC telescope at Skinakas Observatory, Crete, Greece were taken during six nights, 26 -31 August 2019 in the optical B, V, R, and I bands.
Three additional nights of observations were taken with the 0.6 m Helen Sawyer Hogg (HSH) telescope at CASLEO, Argentina (on loan from the University of Toronto, Canada) in B, V, R, and I bands.Additional V-band observations from AAVSO1 (American Association of Variable Star Observers) were carried out from amateur astronomers' two telescopes in Spain and Italy.
The technical parameters and instrumental details are summarized in Table 1.A total of 220 nights of optical photometric observations of TXS 0506+056 were carried out between 21 January 2017 and 9 April 2022.The AAVSO data are included, with most of these observations being done in the V and R bands.In 35 nights, the observation duration is 1h, which we use to study IDV behaviour of the blazar.R band observations were carried out in each of them but they were performed in the B, V, and I bands in only some of the nights.Considering these 35 nights of observations, we obtained data for 35, 14, 7, and 6 nights, respectively, to look for IDV in the R, V, I, and B bands.During 44 nights we have at least 1 frame in all four optical bands, which are useful for studying this portion of the spectral energy distribution (SED).The observation log is provided in Table 2.
For the preliminary processing of the raw data, we used standard procedures of IRAF 2 software, following the steps described below.For image pre-processing, we generated a master bias frame for each observing night.This master bias was subtracted from all twilight flat frames and all source image frames taken on that night.A master flat was generated for each passband by taking the median of all the sky flat frames and then normalizing the master flat.Next, to remove pixel-to-pixel inhomogeneity, the source image was divided by the normalized master flat.Finally, cosmic ray removal was carried out for all source image frames using the IRAF task cosmicrays.To find the instrumental magnitudes of TXS 0506+056 and its comparison stars, we performed the concentric circular aperture photometry technique with the DAOPHOT 3 II software (Stetson 1987;1992).For aperture photometry, we explored four different concentric aperture radii defined in terms of the Full Width at Half Maximum (FWHM): 1 × FWHM, 2 × FWHM, 3 × FWHM, and 4 × FWHM for every night.After examining the results for these different aperture radii, we found that setting the aperture radii = 2 × FWHM provided the best S/N, so we used this value for our final analysis.
Each night we observed more than three local standard stars on the same field.Then we selected three non-varying standard stars (stars A, C and D) from Figure 4 of Bachev et al. (2021) that are of nearly the same magnitude and colour as the source, to avoid any error occurring from differences in the photon statistics in the differential photometry of the source.Since the magnitudes of TXS 0506+056 and the standard stars were obtained simultaneously under the same air mass and weather conditions, there is no need for correction of atmospheric extinction.Finally, a comparison star (Star C) was used to calibrate the instrumental magnitude of the TeV blazar TXS 0506+056.Data were collected from different telescopes in the form of instrumental magnitudes of the blazar and 2 Image Reduction and Analysis Facility (IRAF) is distributed by the National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.
3 Dominion Astrophysical Observatory Photometry reference stars, so that we can apply the same calibration procedures and analysis to all the data sets.

ANALYSIS TECHNIQUES
In this section we briefly explain various analysis techniques we have used to analyse these optical data of the blazar TXS 0506+056 on diverse timescales.Zibecchi et al. 2017;2020).

Power-enhanced F-test
To explore IDV, we used the power-enhanced F -test following the approach of de Diego (2014) and de Diego et al. (2015).In recent studies, this test has been used for finding microvariations in blazars (e.g., Gaur et al. 2015b;Polednikova et al. 2016;Kshama et al. 2017;Pandey et al. 2020a;Dhiman et al. 2023;and references therein).In this test we compare the variance of the source light curve (LC) to the combined variance of all comparably bright comparison stars.In this work, we have three comparison field stars (A, C, and D) (Bachev et al. 2021) from which star C is considered as the reference star, and the remaining two field stars as the comparison stars.For details about our implementation of the powerenhanced F -test see Dhiman et al. (2023).

Nested ANOVA
The nested ANOVA test is an updated ANOVA test that uses several stars as reference stars to generate different DLCs of the blazar.In contrast to power-enhanced F -test, no comparison star is needed in the nested ANOVA test, as all the comparison stars are used as reference stars, so the number of stars in the analysis has increased by one (de Diego et al. 2015;Pandey et al. 2019).The ANOVA test compares the means of dispersion between the groups  V − R colour LCs in magenta, labeled with its observation date and telescope code (from Table 1).
of observations.In our case, we have used three reference stars (A, C, and D) (Bachev et al. 2021) to generate DLCs of the blazar.For details about the nested ANOVA test see Dhiman et al. (2023).
The results of both the statistical tests are given in Table 3, where an LC is conservatively labeled as variable (Var) for IDV only if both the tests found significant variations in it, otherwise it is labeled as non variable (NV), though of course we cannot exclude the presence of weak intrinsic variability even in some of those cases.

Intraday variability amplitude
For each of the variable LCs, we calculated the flux variability amplitude (Amp) in percentage, using the standard equation given by Heidt & Wagner (1996): Here Amax and Amin are the maximum and minimum magnitudes, respectively, in the calibrated LCs of the blazar, while σ is the mean error.The amplitude of variability is also mentioned in the last column of Table 3 for the variable LCs.

Discrete and auto correlation functions
The discrete correlation function (DCF) analysis is used to find possible time-lags and cross-correlations between LCs of different energy bands.This technique was introduced by Edelson & Krolik (1988) and later modified by Hufnagel & Bregman (1992) to produce better error estimates.In general, astronomical LCs are unevenly binned, and, for such LCs this technique is very useful as it can be used for unevenly sampled data.Details about the computation of the DCF we employ here are provided in Pandey et al. (2017) and Dhiman et al. (2023).In computing the DCF, one compares time series from different bands; however, when we correlate a data series with itself, we obtain the auto correlation function (ACF), which consistently exhibits a peak at t = 0.The presence of this prominent peak in a DCF serves as an indicator of the absence of any time lag in the data.The presence of any additional strong ACF peak indicates the presence and value of a variability timescale in the time series data (Rani et al. 2011;Gaur et al. 2015a).

Duty Cycle
The duty cycle (DC) provides a direct estimation of the fraction of time for which a source has shown variability.We have estimated the DC of the blazar TXS 0506+056 by using the standard approach (Romero et al. 1999).For DC calculations, we considered only those LCs which were continuously monitored for 1 hour.For details about DC estimation, see Dhiman et al. (2023).

Intraday Flux and Colour Variability
We   reduce the likelihood of detecting any small variations that might be present.No intraday flux variations were detected in either the B or I passbands, during the relatively few nights for which we collected sufficient measurements.The amplitude of the IDV was estimated for the confirmed variable LCs, as shown in the second last column of Table 3.On December 11, 2020, the R band exhibited the lowest detected variability amplitude, which was only 5.45%, while the largest (15.43%) variation was observed in the V band on November 12, 2020.Typically, the blazar variability amplitude is larger at higher frequencies, as was seen on the one of the nights in which both the R and V band showed variability.However, on some occasions the variability amplitude of blazars at lower frequencies has been found to be comparable to, or even larger than, that at higher frequencies (e.g., Ghosh et al. 2000;Gaur et al. 2015b).
To estimate the IDV timescale of the LCs which have shown genuine variation and are listed in Table 3, we used ACF analysis and plotted the results in Fig. 2. IDV timescales are estimated for November 20, 2020 in V and R bands and listed in the last column of Table 3. From Fig. 2, it is seen that for the rest of the variable IDV LCs, either the timescale is longer than the data length, or ACF plots are too noisy to argue for the presence of a timescale.
In our case there are 35 nights during which we could search for IDV (total duration is 69.75 hrs), with observations lasting between 1 − 6 hrs, and we consider only those observing nights that have at least 10 data points.IDV flux variability plots are presented in Figure 1, and the analysis results are reported in Table 3 We studied colour variations of the TeV blazar TXS 0506+056 on IDV timescales with respect both to time and to V -band magnitude (colour magnitude variation).To do so, we used the 14 nights of data on which quasi-simultaneous observations were carried out: V and R-bands in 7 nights; V , R, I bands in 1 night; and B, V , R, I bands in 6 nights.We calculated the V − R colour indices (CIs) for each pair of V and R magnitudes and plotted these V − R CIs with respect to time in the bottom panel of Figure 1.We performed the enhanced F -test and nested ANOVA test on each night's data set and we found V − R colour variation on only 1 night, November 20, 2020.

Short Term Flux Variability
We divided the long term LC of TXS 0506+056 into 6 Short term LCs, ST1, ST2, ST3, ST4, ST5, and ST6, as plotted in Fig. 3 and the variability results are listed in Table 4.The lengths of these short term LCs, along with their brightest and faintest magnitudes, mean magnitude, and amplitude of flux variations the in R-band are also given there.

Long Term Flux Variability
Figure 4 illustrates the long-term (LT) LCs of TXS 0506+056 in the B, V, R, and I bands over the entire monitoring period.

Spectral Variability
The colour-magnitude (CM) relationship serves as a valuable tool for investigating different variability scenarios and gaining a deeper understanding of the origin of blazar emission.In this study, we conducted a search to identify any potential relationship between the source's colour indices (CIs) with brightness in the R-band and with respect to time.We fitted the plots of the optical (B − V ), (B − I), (V − R) and (R − I) CIs with respect to both R-band magnitude and time using straight lines of the form Y = mX + C as shown in Figs. 5 and 6, respectively (Pandey et al. 2020a).The values of the parameters related to the colour-time and colour-magnitude plots are   respectively provided in the accompanying Tables 5 and 6.In our analysis, we observed modestly significant variations (p < 0.51) in the B−V and B−I colours with respect to the R magnitude.On the other hand, the CIs involving the R-band exhibited very weak trends in the same directions.A positive slope defines a positive correlation between CIs and blazar R magnitude, meaning that the source follows a bluer-when-brighter (BWB) trend, while a negative slope defines redder-when-brighter (RWB) trend (e.g., Gupta et al. 2017;and references therein).We find a negative correlation of the (R − I) CI with time while the B − I colour shows a very weak (about 1.5σ) negative slope.The B −V CI is essentially constant while the V −R one shows a very weak positive trend with time.

Spectral Energy Distribution (SED)
Throughout our observation period, we extracted the optical (BV RI) SEDs of the blazar on 44 nights.These SEDs were derived from quasi-simultaneous observations conducted across all four wavebands.For this, we first dereddened the calibrated B, V , R, and I magnitudes by subtracting the Galactic extinction, with A λ having the following values: AB = 0.392 mag, AV = 0.297 mag, AR = 0.235 mag, and AI = 0.163 mag.The values of A λ were taken from the NASA Extragalactic Database (NED)4 .After applying the necessary dereddening and calibration procedures, the magnitudes in each band were converted into corresponding flux densities Fν that had been corrected for extinction.The optical SEDs of TXS 0506+056, in log(ν) versus log(Fν ) representation, are plotted in Figure 8.We measured the source's faintest and brightest fluxes on August 29, 2021 and February 24, 2020 respectively.
To determine the optical spectral indices, we fitted each SED with a first-order polynomial model of the form log(Fν ) = −α log(ν) + C. The obtained fitting results are presented in Table 8.The values of the spectral indices (α) range from 0.791 ± 0.154 to 1.029 ± 0.194 and their mean was 0.897 ± 0.171.This mean value of the spectral index closely matches previously reported results for TXS 0506+056 (Hwang et al. 2021;Paiano et al. 2018).
In Figure 7, the top and bottom panels display the spectral indices of TXS 0506+056 with respect to time and R-band magnitude, respectively.We fitted each panel in Figure 7 with a first-order polynomial to investigate any variations in the spectral index.The corresponding fitting parameter values are provided in Table 7.The optical spectral index demonstrates a hint of a decreasing trend with time and apparently exhibits a weak positive correlation with the R-band magnitude.However, given the slope of the variation (0.053) and the relatively large error on it (0.022) we cannot have confidence in the reality of this trend.

DISCUSSION AND CONCLUSIONS
The thermal radiation in blazars emitted by the accretion disc is typically overwhelmed by the Doppler-boosted nonthermal radiation originating from the relativistic jet.As a result, any observed variability is more likely to be explained by models based on the relativistic jet (e.g., Agarwal & Gupta 2015;and references therein).During periods when a blazar exhibits a very low flux state, the observed variability could potentially be attributed to hotspots or instabilities occurring on the accretion disc (e.g., Chakrabarti & Wiita 1993;Mangalam & Wiita 1993).IDV/STV observed in the optical bands can be attributed to the presence of turbulence near a shock in the jet, as well as other irregularities within the jet flow resulting from variations in the outflow parameters (e.g., Marscher 2014;Calafut & Wiita 2015;and references therein).
The micro-level flux variations observed in blazars LCs on IDV timescales may be attributed to the turbulent plasma flowing at relativistic speeds within the jet (e.g., Marscher 2014;Pollack et al. 2016) or to mini-jets within the jets (Giannios et al. 2009).Different optical IDV behaviours have been observed in two subclasses of blazars, namely LBLs and HBLs.HBLs display relatively less variability in optical bands on IDV time-scales compared to their amplitude of variability in X-rays and γ-rays (e.g., Heidt & Wagner 1996;Gopal-Krishna et al. 2011;Gaur et al. 2012b;Gupta et al. 2016b;and references therein).The presence of strong magnetic fields within the relativistic jet may be responsible for the different microvariability behaviours observed in the optical bands of both LBLs and HBLs (Romero et al. 1999).The idea is that stronger magnetic fields in HBLs can potentially interrupt the formation of small fluctuations caused by Kelvin-Helmholtz instabilities within relativistic jets.These fluctuations typically interact with shocks in the jets, leading to IDV.However, the generation of very rapid variability can be disrupted if the strength of the magnetic field exceeds the critical value Bc (Romero 2005), where ne is the local electron density, me is the rest mass of the electron, and γ is the bulk Lorentz factor of the flow.The bulk Lorentz factor γ, and the Doppler factor δ are given by δ = [γ(1 − βcos θ] −1 , where βc and θ are the velocity of emitting plasma and jet viewing angle, respectively.For TXS 0506+056, θ, δ, and γ are is found to be 8 • -20 • , 2 -9, and ∼ 5, respectively (Kun et al. 2019;Li et al. 2020;Sumida et al. 2022).The local electron density ne can be estimated using broadband multi-wavelength SED of blazars.However, for TXS 0506+056 a broadband multiwavelength SED is not available, so we took the average value of ne for other possible potential neutrino loud TeV blazars (e.g., Mrk 421, Mrk 501, 1ES 1426+428, PKS 2155-304, 1ES 2344+514, etc.) (Neronov & Semikoz 2002).The electron density ne of these blazars are found to be in the range of 0. Changes in colour or spectral index can aid in understanding of the emission mechanisms in blazars.As mentioned earlier, three distinct types of behaviour can be observed in the colour-magnitude diagram (CMD): redder-when-brighter (RWB), bluer-when-brighter (BWB), and achromatic.The BWB trend is commonly observed in BL Lacs, whereas the RWB trend is typically followed by FSRQs (e.g., Gaur et al. 2012a;2015b).Synchrotron models dominated by one component can explain the BWB behaviour if the energy distribution of injected fresh electrons, which cause an increase in flux, is harder (Kirk et al. 1998).Another possible explanation for BWB behaviour involves precession of the jet, causing variations in the Doppler factor that affect the 'convex' spectrum (Villata et al. 2004b).During our observations, we consistently observed that TXS 0506+056 followed the BWB trend.This suggests that the increasing flux can be attributed to jet synchrotron emission, very possibly indicating an enhancement in particle acceleration efficiency (Agarwal et al. 2015).This BWB trend is seen in Figure 7 which illustrates spectral steepening as the magnitude increases.This trend can also be explained by a two component emission picture: one is a stable component (αconstant) while the other is a variable component with a flatter slope (α1).When the variable component dominates over the stable component then chromatic behaviours are exhibited.On shorter timescales, optical variations are predominantly influenced by strong chromatic components, while longer-term optical variations can be attributed to a mildly chromatic component (Villata et al. 2004a).
The spectral and flux variations observed in TXS 0506+056 may offer valuable insights into the evolutionary processes occurring in radio loud AGN (e.g., Barth et al. 2002;Fan 2003).In an accretion disc-based model, the shortest variability time-scale is connected to the time it takes for light to traverse the variable region which is directly related to the BH mass (Bachev et al. 2005).The likelihood of detecting variability in blazars becomes greater as the duration of observations increases.For instance, extending the observation period from 3 to 6 hours resulted in an increase in the probability of identifying IDV in blazars from 64% to 82% (Gupta & Joshi 2005).
In this study, we conducted an analysis of optical photometric data obtained from seven ground-based telescopes, covering the period from October 2018 to August 2021, for the TeV blazar TXS 0506+056.This study represents the first comprehensive investigation of the temporal and spectral behaviour of this source across a wide range of timescales within the optical domain.We examined a total of 35 R-band, 14 V -band, 7 I-band, and 6 B-band IDV LCs using two powerful and robust statistical methods, namely the nested ANOVA test and the power-enhanced F -test.During the period of our observations, we observed significant variability in the R-band on 8 nights and in the V -band on 2 nights.However, no IDV was detected in the B and I bands for which we took many fewer measurements.Additionally, we observed a variation in the (V − R) colour within only one night throughout the entire observation period.The Duty Cycles (DCs) we measured for the R and V bands are 22.8%, and 14.3%, respectively.So, from our IDV analysis, we conclude that optical LCs of TXS 0506+056 are either constant or show nominal variations on IDV time-scales.The blazar TXS 0506+056 did not show large-amplitude variations during our monitoring period.One important caveat is that the duration of our observations ranges from 1 to 6 hours, so it is certainly possible that if more of our nightly observations had covered longer periods we would have seen more frequent IDV.In a statistical study of the optical IDV of blazars, it was found that chances of detection of IDV is 60 -65% if observations are performed for less than 6h, but if the blazar is observed for more than 6h, the chance of detecting IDV is 80 -85% (Gupta & Joshi 2005).
By analysing the ACFs, we were able to identify evidence of variability timescales on November 20, 2020, in both the V and R bands.The variability timescales were determined to be 3.50 hours in the V band and 3.06 hours in the R band reported in Table 3 and shown in Figure 2. To obtain an upper limit for the size (R) of the emission region, we apply the simple causality constraint, where δ is the Doppler factor.The values of δ for TXS 0506+056 have been estimated to be between 2 and 9 using radio very long baseline interferometry (VLBI) / very long baseline array (VLBA) data in the frequency range 8 GHz to 43 GHz (Kun et al. 2019;Li et al. 2020;Sumida et al. 2022).Using that range of Doppler factor values, taking z = 0.3365 (Paiano et al. 2018), and employing the shortest variability timescale we found (tvar = 11.02ks in the R band), and employing Eqn.
(3), we obtain that the size of the emission region is in the range of 4.95 ×10 14 cm -2.23 ×10 15 cm.
To determine the optical spectral index (α), we constructed optical SEDs using quasi-simultaneous observations in the B, V , R, and I bands at different epochs.The analysis revealed a weak positive correlations between the optical spectral index (α) and R-band magnitude.Over the Long-Term Variability (STV) timescale, the weighted mean value of α was determined to be 0.897 ± 0.171.Since TXS 0506+056 is the best case for a blazar producing neutrino fluxes it should remain a priority target for continuing observations in all bands.

Figure 1 .
Figure 1.Nightly optical light curves of the TeV blazar TXS 0506+056, B-band LCs in blue, V -band LCs in green, R-band LCs in red, I-band LCs in black, V − R colour LCs in magenta, labeled with its observation date and telescope code (from Table1).
observed the TeV blazar TXS 0506+056 sufficiently intensely to investigate intraday flux variation for a span of 35 nights from November 8, 2019 to August 31, 2021.Our observations included quasi-simultaneous monitoring of the source in V and R bands on 7 nights, in I, V , and R bands on one night, and in B, V , R, and I bands on 6 nights.On the remaining 21 nights, we observed the source only in the R band.The calibrated B, V, R, and I band IDV LCs of the blazar TXS 0506+056 are shown in the upper panel of each plot in Figure1, while the lower panels show the V − R colours.IDV during some of the nights seem evident by visual inspection of the LCs.To find the presence of such rapid variability, we performed the statistical tests discussed in sections 3.1 and 3.2, and the results of the analyses are presented in

Figure 2 .
Figure 2. ACF plots of TXS 0506+056 light curves (including one colour index).The blue dashed lines indicate 0 lags and the red dashed lines illustrate intraday timescales.

Figure 3 .
Figure 3.The STV LC of TXS 0506+056 in I, R, V , and B bands.

Figure 6 .
Figure 6.Optical colour variability LCs covering the total observation duration of TXS 0506+056

Figure 8 .
Figure 8.The SED of TXS 0506+056 in I, R, V , and B bands.

Table 1 .
Details of telescopes and instruments used.
To obtain the blazar's intraday variability, we have examined two relatively recently developed statisti- (Gaur et al. 2012b9iego et al. 2015)nhanced F -test and the nested analysis of variance (ANOVA) test(de Diego 2014;de Diego et al. 2015).These tests are usually more reliable and powerful than other statistical tests such as the standard C-test(Romero et al. 1999), F -test (de Diego 2010), χ 2 -test(Gaur et al. 2012b), ANOVA test  (de Diego et al. 1998) because these involve several comparison stars (but see also

Table 3
. The intra-day LCs of these 35 nights have observation duration 1h, and are displayed in Figure1.The IDV analysis results of these LCs are reported in Table3.However, note that the errors in V band LCs are roughly twice in comparison to R band, and therefore

Table 4 .
Short term LCs values The plot depicts the nightly averaged magnitudes in the respective bands as a function of time.During our monitoring period the source was detected in the brightest state of R = 13.98 mag on February 24, 2020 while the faintest level detected was R = 15.16 mag on August 29, 2021.The mean magnitudes were 15.62, 15.02, 14.55 and 13.91 in B, V, R, and I bands, respectively.The presence of variability on LT timescales is clearly evident across all optical wavelengths.Using Equation 1, we have estimated very similar variability amplitudes of 117.6%, 137.9%, 117.6%, and 113.8% in B, V, R, and I bands, respectively.

Table 5 .
Colour-Magnitude Dependencies and Colour-Magnitude Correlation Coefficients on LTV Timescales m 1 = slope and c 1 = intercept of CI against R-mag; r 1 = Correlation coefficient; p 1 = null hypothesis probability a

Table 6 .
Colour Variation with Respect to Time on LTV Timescales

Table 7 .
Variation of optical spectral index (α), with respect to R-magnitude and LTV timescale

Table 8 .
Straight-line Fits to Optical SEDs of TeV Blazar TXS 0506+056 MAGICCollaboration et al. 2020).Considering these parameter values, we get B ≈ 0.1 G.As TXS 0506+056 is a TeV blazar, it is expected to have B > Bc that, within this scenario, would inhibit the development of the small scale structures that could yield IDV in the optical LCs.