Detection of an intranight optical hard-lag with colour variability in blazar PKS 0735+178

Blazars are a highly variable subclass of active galactic nuclei that have been observed to vary significantly during a single night. This intranight variability remains a debated phenomenon, with various mechanisms proposed to explain the behaviour including jet energy density evolution or system geometric changes. We present the results of an intranight optical monitoring campaign of four blazars: TXS 0506+056, OJ287, PKS 0735+178, and OJ248 using the Carlos S\'anchez Telescope. We detect significant but colourless behaviour in OJ287 and both bluer- and redder-when-brighter colour trends in PKS 0735+178. Additionally, the g band shows a lag of ~10 min with respect to the r,i,z_s bands for PKS 0735+178 on 2023 January 17. This unexpected hard-lag in PKS 0735+178 is not in accordance with the standard synchrotron shock cooling model (which would predict a soft lag) and instead suggests the variability may be a result of changes in the jet's electron energy density distribution, with energy injection from Fermi acceleration processes into a post-shocked medium.


INTRODUCTION
At the centre of most galaxies resides a supermassive black hole (SMBH) and if subject to accretion, this SMBH can be referred to as an active galactic nucleus (AGN).A blazar is an AGN where the relativistic jet emanating from the polar regions of the black hole is oriented along our line of sight within ≲ 15 • (Hovatta et al. 2009).The resulting relativistic beaming produces highly variable emission across the electromagnetic spectrum (Urry 1998).
Blazars can be categorised into two classes based on the features within their optical spectra.BL Lac-types have near featureless optical spectra whereas flat spectrum radio quasars (FSRQs) show emission lines with EW ≥ 5 Å (Stickel et al. 1991).The spectral energy distribution of blazars has a distinct double-hump structure with a lower-energy peak at IR to X-ray frequencies and a higher-energy peak at X-ray to HE -rays (Fossati et al. 1998).The lower-energy peak is attributed to synchrotron emission from the jet, whereas the source of the higher-energy peak is still debated (Ghisellini et al. 1998;Prandini & Ghisellini 2022).This second peak can be explained through leptonic models, where the emission is attributed to lower-energy seed photons being up-scattered by relativistic electrons through inverse-Compton scattering (Böttcher et al. 2013).The source of these seed photons can be from within (synchrotron self-Compton; Maraschi et al. 1992) or from outside (external Compton; ★ E-mail: c.mccall@2017.ljmu.ac.ukDermer & Schlickeiser 1993) the jet.This means relationships between the lower-and higher-energy regimes are likely, as their origin is the same population of electrons.Conversely, hadronic models involving interactions between protons and/or mesons might explain the high energy variability independently of the lower-energy emission (Böttcher et al. 2013).These hadronic models are also important for the production of HE neutrinos (Osorio et al. 2023).
The location of the lower-energy peak can be used to distinguish between different blazar classes and BL Lac-type subclasses.Low synchrotron peaked (LSP) sources have a synchrotron peak < 10 14 Hz, corresponding to IR emission, and can be FSRQs or BL Lacs (low-frequency peaked BL Lacs; LBLs).Only BL Lacs have synchrotron peaks at higher frequencies.High synchrotron peaked (HSP or high-frequency peaked BL Lacs; HBLs) sources have synchrotron peaks > 10 15 Hz in UV or X-ray; those with a synchrotron peak frequency in the optical regime between 10 14 ≤  ≤ 10 15 Hz are classed as intermediate synchrotron peaked (ISP) sources or intermediate-frequency peaked BL Lacs (IBLs) (Urry & Padovani 1995;Abdo et al. 2010).
Blazars can show variability on a variety of time-scales, from years down to minutes, with the latter being referred to as intranight variability.The intranight optical variability (INOV) of blazars has been the focus of numerous campaigns over the past 30 yr (Miller et al. 1989;Sagar et al. 2004;Chand et al. 2021), but is still not a wellunderstood characteristic of their emission.There have been many models proposed to explain the different types of variability such as quasi-periodic oscillations (QPOs) (Jorstad et al. 2022), micro-flares (Bhatta et al. 2015), and gradual flux changes (Subbu Ulaganatha Pandian et al. 2022), observed over hour-long time-scales.These models can predict time lags between different wavebands, colour evolution, and polarization degree and angle changes.Such models include geometric changes related to the jet angle and Doppler factors of emitting blobs (Gopal-Krishna & Wiita 1992), intrinsic changes relating to particle energy distributions (Bachev et al. 2012;Bachev 2015), and extrinsic changes including microlensing from objects outside the blazar system along our line of sight (Paczynski 1996).
Colour evolution in blazars generally takes one of two forms: bluerwhen-brighter (BWB) or redder-when-brighter (RWB).In general, BL Lac-type sources tend to exhibit BWB trends and where this is not the case and a RWB trend is seen, the objects tend to be host-galaxydominated (Negi et al. 2022).FSRQs can be found showing both BWB and RWB trends, but the proportion of sources that exhibit the latter is much greater than among BL Lacs (Zhang et al. 2015;Negi et al. 2022).In any case, the colour changes of blazars can be used to determine the origin of the emission and refine emission mechanisms.The BWB chromatism may arise from a number of different mechanisms.One proposed mechanism is synchrotron cooling of internal shock accelerated electrons, where higher-energy electrons cool faster making bluer light appear more variable than redder light (Kirk et al. 1998).Another is the one-component synchrotron model where an increase in the energy output of the blazar increases the average particle energy and thus the frequency, making the blazar appear bluer-when-brighter (Fiorucci et al. 2004).When in faint states or periods of jet quiescence, redder emission and variability from the accretion disc may become visible leading to RWB behaviour; this is true for both BL Lac and FSRQ-type blazars.A potential explanation of why more FSRQs show general RWB behaviour may come from the flattening of their spectral slope at optical frequencies from the presence of the 'UV bump'.The RWB behaviour, or steepening of the composite thermal and non-thermal spectrum, could be explained if the non-thermal component had a greater contribution towards the total flux during brightening (Gu et al. 2006).
Similarly to colour variability, time lags (or a lack thereof) can help constrain emission processes.Time-lags between different frequencies have been observed in numerous studies (Chatterjee et al. 2008;Gaur et al. 2012;Liodakis et al. 2018), but inter-band lags within the optical waveband are much less commonly observed, with the first proposed detection in 2009 (Wu et al. 2009).The internal shock model (Kirk et al. 1998) predicts time-lags across all synchrotron emission frequencies and therefore within the optical waveband.The rate of synchrotron cooling is frequency dependent, meaning variability at higher-energy precedes that at lower energies (Kirk et al. 1998).
INOV including colour evolution and inter-band time-lags is difficult to detect for a number of reasons, both physical to the system and logistically in terms of observational cadence, and can result in data with irregular or limited time resolution (Böttcher & Dermer 2010).
In this paper, we present the results of a five-night INOV monitoring campaign executed over the period 2023 January 15-19 on four blazars: TXS 0506+056, OJ287, PKS 0735+178 and OJ248.These sources were chosen due to a combination of their historic degrees of INOV and their observability during the campaign.The sample was kept small to ensure the data were sufficient to observe hour-long variability time-scales with high sampling.
In Section 2 we describe the facilities and instruments used and photometric analysis procedures including statistical processes.In Section 3 we describe the individual sources and their historical behaviour.In Sections 4, 5, and 6 we present the results of the correlation analyses including variability indicators, colour behaviour, and time-lags.In Section 7 we discuss the implication of the results from the analysis.

Carlos Sánchez Telescope -MuSCAT2
Observations were taken with the four-colour simultaneous imager MuSCAT2 (Narita et al. 2019) located on the 1.52 m Carlos Sánchez Telescope (TCS) at the Teide Observatory, Tenerife.The TCS has a Dall-Kirkham configuration with a focal length of  /13.8 and the MuSCAT2 instrument is fixed at the Cassegrain focus.MuSCAT2 achieves four-colour simultaneous observations through the use of three dichroic mirrors to separate the light into four wavelength bands (g (400-550nm), r (550-700nm), i (700-820nm),  s (820-920nm) where the subscript 's' here denotes the 'shorter' waveband range to a traditional z band filter) to be detected by four fast-readout PIXIS CCD cameras.The derived pixel scales of ∼ 0.44 arcsec/pixel across each band correspond to a field of view (FOV) of 7.4 × 7.4 arcsec 2 in all filters (Narita et al. 2019).
A summary of the TCS observations is presented in Table 1.Each source was observed on three separate nights and, in general, observations were interleaved for two sources with a typical observing sequence of 10 frames per source with a 30-second exposure time (the longest that could be executed without autoguiding; imposed by an autoguider failure).

Liverpool Telescope -MOPTOP
Supplementary observations were carried out with the MOPTOP (Shrestha et al. 2020) polarimeter on the fully robotic, 2 m Liverpool Telescope (LT) located at the Observatorio del Roque de los Muchachos, La Palma.The LT has a Ritchey-Chrétien Cassegrain design with MOPTOP fitted at one of the science fold ports.
MOPTOP boasts a dual-beam configuration utilising a continuously rotating half-wave plate and two fast readout, very low noise CMOS cameras.Together, these allow MOPTOP to achieve high sensitivity and time resolution while minimizing systematic errors.MOPTOP has a 7 × 7 arcsec 2 FOV (Shrestha et al. 2020).The data were taken in  (380--520 nm),  (490--570 nm) and  (580--695 nm) filters quasi-simultaneously (observations with different filters taken in succession).
The observations with MOPTOP were taken much less frequently than those of MuSCAT2, aiming for a few observations each night wherever possible.Observations were taken on 2023 January 15-16 before poor weather conditions in La Palma closed the observatory for the remainder of the campaign.

Data Reduction
The data were calibrated and analysed using the astropy Python package and standard differential photometry techniques.MuSCAT2 frames required bias/dark subtraction and flat fielding and this was performed using calibration frames taken before/after observations each night.The data also required a WCS fit which was performed using the Astrometry.netAPI (Lang et al. 2010).MOPTOP reduction is carried out using an automated data reduction pipeline running at the telescope which provides bias/dark subtracted, flat-fielded, and WCS fitted frames (Smith et al. 2016).
The MuSCAT2 photometric data were calibrated by calculating the weighted average zero point in each frame using estimates from five in-frame calibration stars with known  s magnitudes from SDSS, Pan-STARRS and APASS catalogues (Abdurro'uf et al. 2022;Flewelling et al. 2020;Henden et al. 2018).This zero point was then used to calibrate the data for the source.
MOPTOP photometric data were calibrated using an in-frame reference star with known   magnitudes.Polarimetric data were calibrated using observations of zero-polarized standard stars to calculate instrumental polarization values ( 0 and  0 ) and polarized standard stars to find the instrumental position angle values.
Before applying statistical methods and generating light curves, the data were sigma-clipped using the Python package astropy as a form of quality control.In this process, data are removed if they are above or below three times the standard deviation from the median value.This is an iterative process that produces updated median and standard deviation values until no further values are removed.Applying this to the flux values at each epoch, disproportionate variability between filters caused by erroneous pixels or poor sky conditions was removed.The final amount of data removed equates to less than 10 per cent (and in most cases less than 5 per cent) per source per epoch.
We note that no correction for host galaxy contribution to the source magnitudes or colours has been applied.In general blazars of subclasses LSP and ISP outshine their host galaxies by several orders of magnitude, so significant host galaxy emission is only observed when the AGN is in a very low state in combination with deep photometric imaging (Nilsson, K. et al. 2012;Gaur 2014;Olguín-Iglesias et al. 2016).None of our sources were in such a state.

TXS 0506+056
TXS 0506+056 is situated at a redshift of 0.337 (Paiano et al. 2018) and until recently was classified as an LSP BL Lac (LBL) object (Fan et al. 2014).It has been proposed that TXS 0506+056 may belong to a subclass of blazars named masquerading BL Lacs, or blue FSRQs (Ghisellini et al. 2012;Padovani et al. 2019;Lewis et al. 2021), where the Doppler-boosted synchrotron radiation in the relativistic jet is bright enough to outshine the broad line region (Rodrigues et al. 2019;Rajagopal et al. 2020).This would make the FSRQ appear as a BL Lac object due to the apparent lack of emission lines.
In 2017 September, TXS 0506+056 was found to be in a consistent location with the IceCube neutrino event EHE 170922A (Kopper & Blaufuss 2017) to within a 3 significance level (IceCube et al. 2018) while in a state of heightened -ray activity (Tanaka et al. 2017).Since then, the object has been the subject of various studies over different frequencies and time-scales (see Keivani et al. 2018;Bachev et al. 2021;Acciari et al. 2022).
We observed TXS 0506+056 over three nights, totalling 7.45 h in the  s filters.We obtained complementary photo-polarimetric data with MOPTOP on the LT in the   bands to assess the polarimetric state of the source during the observing campaign.We also took Fermi -ray data from the Light Curve Repository (LCR) (Abdollahi et al. 2023).The Fermi data showed TXS 0506+056 to be in a ray state a little higher than its median flux level over all time, at 6.63 ± 6.05 × 10 −8 0.1-100 GeV ph cm −2 s −1 (median level: 8.30 ± 0.23×10 −8 0.1-100 GeV ph cm −2 s −1 ).Its optical linear polarization degree varied between roughly 9 and 16 per cent and the polarization angle varied between approximately 150 • and 170 • across all   bands.

OJ287
OJ287 is located at a redshift of 0.306 (Sitko & Junkkarinen 1985) and is a well-known LSP BL Lac-type object (LBL) (Nilsson et al. 2018).OJ287 is one of the best binary supermassive black hole candidates with a double-peaked outburst period of roughly 12 yr (Sillanpaa et al. 1988(Sillanpaa et al. , 1996)).The observed optical outbursts date back over 130 yr, corresponding to the interaction of the secondary black hole with the primary's accretion disc.A second, longer period of approximately 60 yrs has also been reported (Valtonen et al. 2006) which is thought to arise from the orbital precession inducing a precession into the accretion disc (Katz 1997;Sundelius et al. 1997) and causing a subsequent wobble in the jet angle.
We observed OJ287 over three nights, totalling 10.86 h in the  s filters.Additional MOPTOP photo-polarimetric data showed the optical linear polarization degree varied between roughly 10 and 25 per cent and its polarization angle varied between approximately 160 • and 175 • across all   bands.Fermi -ray data from the Light Curve Repository (LCR) (Abdollahi et al. 2023) showed OJ287 to be in a -ray state a little lower than its median flux level over all time, at 3.77 ± 2.93 × 10 −8 0.1-100 GeV ph cm −2 s −1 (median level: 5.94 ± 0.31 × 10 −8 0.1-100 GeV ph cm −2 s −1 ).

PKS 0735+178
PKS 0735+178 is an ISP BL Lac-type object with a disputed redshift as a result of its featureless optical spectrum.A lower limit of  ≥ 0.424 was first proposed by Carswell et al. (1974) after the detection of a strong absorption feature, and has since been refined using the detection of the host galaxy using deep  band imaging when the central engine was in a faint state ( = 0.45 ± 0.06 (Nilsson, K. et al. 2012)), and the surrounding galaxies ( ∼ 0.65 (Stickel et al. 1993;Falomo et al. 2021)).Sahakyan et al. (2023) compare the multiwavelength data to that of TXS 0506+056, another blazar neutrino candidate and utilise conditions set out in Padovani et al. (2019) to conclude that like TXS 0506+056, PKS 0735+178 may also be a masquerading BL Lac as a result of a radio power  1.4 GHz > 10 26 W Hz −1 and a -ray Eddington ratio   /  ≳ 0.1.

TEMPORAL VARIABILITY
The  s light curves for each source on a given night are shown in Figs B1, B2, B3 and B4 for TXS 0506+056, OJ287, PKS 0735+178 and OJ248 respectively.The source name and night of observation are given above each plot.
In order to quantify variability in the light curves and to disentangle intrinsic variability from noise, we employ several statistical tests to determine the variability likelihood.Specifically, we calculate the variability amplitude and fractional variability and perform chisquared and enhanced F-test analyses.These tests are detailed in Sections 4.1, 4.2, 4.3, and 4.4.

Variability Amplitude
Variability amplitude is defined in Heidt & Wagner (1996) as where  max and  min are the maximum and minimum observed values, and ⟨ err ⟩ is the median measurement error.The percentage variability amplitude, VA per is given in Romero et al. (1999) as where ⟨⟩ is the average observed value.Its error, ΔVA per , is given in Singh et al. (2018) as where  err,max and  err,min are the errors on the maximum and minimum values, respectively.The variability amplitude aims to provide a quantification of the absolute range of variability of a given source by simply looking at the range of magnitudes outside of the scatter from the measurement errors.The values obtained via the variability amplitude calculations can be seen in column four in Table 2.We find variability amplitudes ranging from 0.167 per cent up to 1.456 per cent across all sources, dates and filters.The ratio between the error and percentage variability amplitude shows that in 12/48 cases the percentage variability amplitudes are associated with large errors (where the ratio is greater than 3).We note that nine of these cases are attributed to the source OJ248; likely due to it being the faintest of the sample.

Fractional Variability
Fractional variability, described fully in Schleicher et al. (2019), is a method of quantifying variability intensity while accounting for measurement uncertainties.It differs from the variability amplitude by considering the variability relative to the mean brightness level.It is defined in Edelson et al. (2002) as where  2 is the variance of the data set, ⟨ 2 err ⟩ is the median square error, and ⟨⟩ is the median value.It can also be given as the square root of the normalised excess variance ( 2 NXS ).Its associated error is given in Poutanen et al. (2008) as where Δ 2   is the error on the normalised excess variance.This error is given in Vaughan et al. (2003) as where N is the number of data points in the sample.It follows that if the variance is less than the average square error,  2 < ⟨ 2 err ⟩, a real value cannot be computed and will be denoted as < 0, indicating detection of insignificant variability.Where sources had  var > Δ var , this test is deemed to show that an object has shown significant variability.
The fractional variability values are shown in column five in Table 2.We find 12/48 instances across all sources, dates and filters where significant levels of variability have been detected.These detections correspond to OJ287 on 2023 January 15, and PKS 0735+178 on 2023 January 15 and 17 across all filters, with the most significant detections corresponding to PKS 0735+178 on 2023 January 17.

Chi-squared
Chi-squared ( 2 ), as used in Zeng et al. (2017), is given by where   and  err, are the individual values and errors respectively within the data set, and ⟨⟩ is the median value.Its critical value was determined at the 99.9 per cent confidence level ( = 0.001) with the degrees of freedom being equal to the number of data points.Where the value is greater than the critical value, significant variability has been detected. 2 is a useful metric as it quantifies the levels of variability about the median values.Incorporating the critical value allows us to determine the significance of the value.
The  2 values together with critical values are shown in column six in Table 2.With the  2 test, we detect significant variability in 27/48 instances across all sources, dates and filters.In most cases, non-detection is consistent across filters per source per date.We find consistent variability detections across all filters for TXS 0506+056 on 2023 January 15, OJ287 on 2023 January 15 and 18, and PKS 0735+178 on 2023 January 15 and 17.

Enhanced F-test
The enhanced F-test ( enh ) is given in de Diego (2014) and aims to quantify the variability of a target while incorporating the variability of multiple reference stars (Subbu Ulaganatha Pandian et al. 2022).It is given as where  2 blazar is the variance of the blazar and  2 star is the combined variance of the comparison stars. 2 star is defined as where   is the number of observations of the  th comparison star, k is the number of comparison stars and  , is the scaled square deviation. , is given as where  , is the magnitudes of each comparison star, ⟨  ⟩ is the mean magnitude of the comparison star and   is the scaling factor to account for the different SNRs of the comparison stars.  is given as where ⟨ 2 blazar ⟩ and ⟨ 2   ⟩ are the average square errors of the blazar and comparison star magnitudes respectively.Its critical value was determined at the 99.9 per cent confidence level ( = 0.001) with the degrees of freedom for the blazars being one less than the number of observations, and the degrees of freedom for the comparison stars being the sum of one less than the number of observations for each comparison star.
Finally, the  enh values with critical values are shown in column seven in Table 2.They show significant detections of variability in 15/48 instances across all sources, dates and filters.12 of these detections correspond to observations in all filters of OJ287 on 2023 January 15, and PKS 0735+178 on 2023 January 15 and 17.

Temporal Variability Summary
We summarise the results of each test in Table 2 and determine epochs of variability, we look at the results of the fractional variability,  2 , and enhanced F-test analyses and determine variable epochs if all three tests are passed.If one or fewer tests were met, we considered there to be no evidence for intranight variability in these sources/epochs.Overall all TXS 0506+056 and OJ248 epochs display no significant levels of intranight variability.There is possible variability from TXS 0506+056 on 2023 January 16.Given that this possible variability only occurs in the  band, we consider this likely a result of scatter in the data.OJ287 shows one epoch of significant variability on 2023 January 15.PKS 0735+178 shows two epochs of significant variability on 2023 January 15 and 17 and one of possible significant variability on 2023 January 16.The latter consists of very few observations due to poor observing conditions so, based on the light curve, it is possible that this is a false detection.We show the light curves with statistically significant variability for OJ287 and PKS 0735+178 in Figs 1, 2, and 3 and will discuss these variable epochs further in later sections.

COLOUR VARIABILITY
We test for colour variability by investigating how the  −  s colour changes with respect to  band magnitude and with time.We use the Spearman rank correlation coefficients to quantify the level of monotonic variability observed in each source on a given night.We use  = 0.05 for the significance coefficient, p, implying a 95 per cent confidence interval and assign the strength of the correlation, c, by the ranges specified in Table 3.We also utilise the enhanced Ftest to account for the variability of the reference stars as previously described.
The Spearman rank correlation coefficients and enhanced F-test values for each set of  −  s vs  data are shown in Table 4 (full table for all epochs available in Table B1 with the corresponding plots in Fig. B5 in the appendix).There are two epochs that show significant colour variability during observations; PKS 0735+178 on the nights of 2023 January 15 and 17.The former shows a positive correlation with a strength of 0.33, indicating a weak correlation.The positive nature of this correlation implies as the source gets brighter, it also gets redder in colour.Conversely, the latter date shows a negative correlation with a strength of 0.35, again indicating a weak correlation, although the negative nature this time indicates as the source gets brighter it also gets bluer.From left to right, the source name, date of observation and filter is given.VA represents the variability amplitude with error as described in section 4.1,  var is the fractional variability with error as described in section 4.2,  2 is the Chi-squared value given with the critical value as described in section 4.3 and  enh is the enhanced F-test value with the critical value as described in section 4.4.The final column describes whether or not the source was deemed variable in a particular filter on a given night.Of the three variability tests ( var ,  2 and  enh ), if all three showed variability the source was deemed variable (V), two meant possibly variable (PV), and one or none meant likely not variable (NV).The magnitude of c shows whether the correlation is positive or negative.
To confirm these results, we also obtained correlation statistics on the slope of the optical spectral energy distribution (SED) vs the  band magnitude.The slope was obtained by fitting a line through the , ,  s band magnitudes at each epoch.This analysis confirmed the same sources and epochs to show significant variability as the colour analysis (see Appendix A for details).

TIME-LAG ANALYSIS
We test for the possibility of a time lag between  s bands on the nights where sources show statistically significant variability.This would be indicative of a shock or any energy density evolution within the jet, and allow us to rule out geometric variability processes like Doppler factor evolution of spiralling emitting regions.The variability must occur over time-scales less than the duration of the observations (minima and maxima within the lightcurves), which allows us to match up light curve features between bands and test for intra-band lags.Only one source and night fit these criteria, PKS 0735+178 on 2023 January 17.To perform the lag analysis, we utilise the Discrete Correlation Function (DCF) which provides an estimate for the time lag between two unevenly sampled time series without the need for interpolation, while accounting for the effects of correlated errors (Edelson & Krolik 1988).It is defined by where (  ,   ) are the observations, (⟨⟩, ⟨⟩) are the mean value from each light curve, (  ,   ) are the standard deviation of each light curve, and (⟨Δ⟩, ⟨Δ⟩) are the median error values (Liodakis et al. 2018).To find the DCF value associated with each time shift, , we average over the number of (  ,   ) pairs, , where  − Δ 2 < Δ   <  + Δ 2 or in this case, the mean    value What also sets the DCF apart from other correlation methods is that a standard error on  () can be given by assuming the individual UDCF ij values within a bin are uncorrelated.We investigate the possibility of a lag within ± 60 min.While analysing the data using the DCF, its limitations in accounting for regularly unevenly sampled data became apparent.The data consists of an observing sequence over ∼ 5 min before a ∼ 10 min break whilst observing a second target.When performing the DCF, this periodically resulted in a large decrease in the number of overlapping bins, zero in some instances, within Δ 2 .This is seen in the correlation curves (Fig. 4) as periodic peaks and drops in the coefficient values.Fig. 4 shows the results of the DCF on the data from PKS 0735+178 on the night of 2023 January 17 on each  s light curve with respect to the  light curve.In this configuration, a positive lag implies  leading the other bands and a negative lag implies  lagging the other bands.The solid curve shows a Gaussian fit to the DCF correlation values, calculated to offset the structure induced by the periodic number of overlapping bins.The dotted line shows the peak of the Gaussian curve, and therefore the lag value.It shows a significant non-zero lag in each  s light curve with respect to . Between the three bands, the lags are all consistent, with a mean value of −6.94±1.43min.The uncertainty of 1.43 min is the average cadence of the observations and was chosen as the larger value of average cadence and error on the Gaussian peak.
In order to check the significance of the induced correlation curve structure, and to mitigate the scatter in the light curves, we also calculated the DCF after fitting a curve to the data.We fit each light curve using the GaussianProcessRegressor module from scikitlearn in Python (Pedregosa et al. 2011) using the Rational Quadratic kernel.Calculating the DCF on this fitted curve and following the same steps as outlined previously, we obtain the results shown in Fig. 5.We keep the same uncertainties (1.43 min) to reflect the original data cadence.The results for  and  are consistent with the values obtained previously, but the lags obtained in  and  s are significantly larger at −11.68 ± 1.43 min and −15.48 ± 1.43 min, respectively.

DISCUSSION
Blazar intranight variability is thought to arise from geometric changes within the blazar jet; such as the Doppler factor variability of an emitting region travelling in a helical motion in the jet, from the evolution of an emitting region through the jet or from the acceleration/cooling of particles.Additionally, it is entirely possible for the observed behaviour to be a combination of multiple emitting regions or different processes occurring simultaneously.
The mechanism behind Doppler factor variability involves an emitting region, or 'blob', of density inhomogeneity travelling helically along the jet.This causes quasi-periodic oscillations (QPOs) in the light curve resulting from the apparent changing Doppler factor and subsequent bulk Lorentz factor (Camenzind & Krockenberger 1992;Mohan & Mangalam 2015;Bachev et al. 2023).On intranight timescales, this behaviour would present across the optical regime as multiple brightness peaks, depending on the number of blobs, where individual peaks would be observed with no colour changes or timelags (Papadakis et al. 2004;Bachev 2015).If the origin of the variability was many emitting blobs, each with differing SEDs, then one might expect the emission of different blobs to dominate at different times and subsequently cause rapid colour changes in addition to the brightness changes (Bachev et al. 2023).This variability, however, is a relativistic effect rather than any change in the emission output of the source.
Changes in the intrinsic luminosity of the source on intranight time-scales can be attributed to processes such as shocks or magnetic reconnection in the jet.These processes involve a uniform injection of fresh, more energetic electrons which evolve as a function of their energy distribution, where harder electrons cool faster (Urry et al. 1997).This may produce intra-and inter-band time-lags, which can determine cooling times and constrain the homogeneous synchrotron model (Tavecchio et al. 1998).An evolving energy distribution may also produce colour variability (Papadakis et al. 2004).Additionally, emission at optical frequencies can trace slightly different parts of the SED depending on the location of the synchrotron peak.For LSP sources (three of our sources), optical frequencies trace the falling region of the synchrotron peak which means redder frequencies map higher-energy emission and may produce faster-evolving variability, causing colour variability and time-lags between wavebands.Conversely, for HSP blazars, optical frequencies trace the rising part of the SED so one would expect the bluer frequencies to evolve faster (Subbu Ulaganatha Pandian et al. 2022).
In our work, we found that TXS 0506+056 and OJ248 showed no evidence of variability in the epochs studied.OJ248 is the faintest object in our sample and would have benefitted from longer exposure times for better signal-to-noise had the autoguider on the TCS been available.TXS 0506+056 showed significant, weak, colour variability on 2023 January 15, which may be due to the scatter in the data.
OJ287 showed evidence of significant flux variability on the night of 2023 January 15, but no significant changes in colour.There are no significant short-time-scale features in the light curve, and the observed variability consists of a gradual decrease in the brightness over the ∼ 6 hours of observing.
PKS 0735+178 displayed significant variability on two out of three nights, including significant colour correlations showing both redderwhen-brighter and bluer-when-brighter behaviour.Additionally, on the night when BWB colour variability was observed, a hard-lag of order 10 min was detected.
If the hard-lags observed in PKS 0735+178 are caused by the evolution of the electron energy distribution, different shock-in-jet processes can be examined to explain the variability.When the acceleration time-scale during the shock is much less than the post-shock cooling time-scale, ie  acc ≪  cool , soft lags are expected, where the lower energy emission (red) lag behind the higher energy emission (blue).Conversely, when the acceleration time-scale is comparable to the cooling time-scale, ie. acc ≈  cool , hard-lags are expected, where the lower energy emission precedes the higher energy emission (Zhang et al. 2002).
In order to achieve a hard-lag, and produce comparable acceleration and cooling time-scales, an energy injection is required to accelerate electrons after the shock has passed, rather than allowing the shocked particles to cool, which results in soft lags (Mastichiadis & Moraitis 2008).Injecting energy into the post-shocked medium can be achieved using second-order Fermi acceleration processes.Kalita et al. (2023) describe how turbulent magnetic fields built behind a shock travelling through an inhomogeneous medium can produce these processes, resulting in acceleration of the post-shock particles via magnetic reconnection.In this scenario, energy is released to the surrounding particles through the interaction of magnetic field lines with opposite polarity.
While we cannot make a firm conclusion as to the nature of the detected INOV in PKS 0735+178, the detection of a hard-lag favours changes to the electron energy distributions and the internal shock model over any geometric changes.

CONCLUSION
We performed simultaneous , , ,  s photometric observations using MuSCAT2 on the Carlos Sánchez Telescope to study the intranight optical variability of four -ray bright blazars.Our analysis consisted of employing several statistical methods to test for the presence of variability on time-scales of a few hours.Additionally, the DCF was used to test for intra-band time lags between bands , , and  s with respect to band .We found: • TXS 0506+056 and OJ248 showed no evidence for intranight variability on any night.
• OJ287 showed evidence for intranight variability on 2023 January 15.The nature of this variability was a gradual change, around 0.1 magnitudes over 7 hours, and was not accompanied by any significant changes in colour.
• PKS 0735+178 showed evidence for intranight variability on two occasions along with changes in colour; showing both a redderwhen-brighter and a bluer-when-brighter correlation on different dates.
• PKS 0735+178 showed a time lag where the  band lags the , ,  s bands by around 10 min.This suggests the variability may arise from changes in the electron energy-density distribution.
Further observations of blazars during all activity states at high cadences can confirm whether intra-band hard-lags across optical frequencies are a more common feature than previously thought.This would provide strong evidence for changes in the jet's energy density as the cause for INOV in blazars.

Figure 4 .
Figure 4. DCF coefficients testing for a lag on the data from PKS 0735+178 on the night of 2023 January 17.The coefficients (blue, green, red, and purple points for filters ,  , ,  s respectively) are fitted with a Gaussian (black line) to find the peak.This peak value (vertical dotted line) is shown in the legend with an uncertainty.

Figure 5 .
Figure 5.As Fig. 4 but with the fitted data from PKS 0735+178 on the night of 2023 January 17.

Figure B5 .
Figure B5.Colour-magnitude ( vs - s magnitudes) diagrams for each of the four blazars (different rows) on different nights (different columns) as indicated above each plot.Also present above each plot are the corresponding Spearman rank correlation coefficients and significance values.

Table 1 .
List of the blazars used in this analysis including the source Right Ascension (), Declination ( ), type, redshift (), TCS observation date, hours observed (), and the number of observations ( ) in ,  , ,  s bands.2015) including data at radio, NIR, optical, X-ray and -ray frequencies.A large multiwavelength flare was observed in late 2012 by WEBT, Swift (D'Ammando & Orienti 2012) and Fermi (Orienti & D'Ammando 2012) as well as an additional radio outburst in late 2010 and optical-NIR flare in early 2007.
. A long-term multiwavelength analysis of this source was performed by the GASP-WEBT Collaboration from 2006-2013 (Carnerero et al.

Table 3 .
Light curves for OJ287 on the night of 2023 January 15.Panels correspond to ,  , ,  s data separately, from top to bottom.Correlation strengths for Spearman rank correlation coefficients.

Table 4 .
As Fig.1, but for PKS 0735+178 on the night of 2023 January 17.Colour variability statistics for variable sources on a given night. and  refer to the Spearman rank correlation coefficients (significance and strength respectively),  enh is the enhanced F-test value with the critical value as described in section 4.4, and the final column describes whether or not the colour of the source was deemed variable on the given night.If  < 0.05 and  enh >  crit the source was deemed variable (V), otherwise not variable (NV).

Table B1 .
SED gradient using , , and  z band magnitudes against  band magnitude for each of the four blazars (different rows) on different nights (different columns) as indicated above each plot.Spearman rank correlation coefficients and significance values are also shown above each plot.Light curves of TXS 0506+056 on the nights of 2023 January 15, 2023 January 16, and 2023 January 18.The panels of each of the three plots correspond to , , , and  s filters, from top to bottom.Figure B2.As Fig.B1, but for OJ287 on the nights of 2023 January 15, 2023 January 18, and 2023 January 19.As Fig.B1, but for PKS 0735+178 on the nights of 2023 January 15, 2023 January 16, and 2023 January 17.As Fig.B1, but for OJ248 on the nights of 2023 January 2023 January 17, and 2023 January 19.Colour variability statistics for each source on a given night. and  refer to the Spearman rank correlation coefficients (significance and strength respectively),  enh is the enhanced F-test value with the critical value as described in section 4.4, and the final column describes whether or not the colour of the source was deemed variable on the given night.If  < 0.05 and  enh >  crit the source was deemed variable (V), otherwise not variable (NV).