Free Access
Issue
A&A
Volume 645, January 2021
Article Number A137
Number of page(s) 20
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202039301
Published online 26 January 2021

© ESO 2021

1. Introduction

Blazars are a sub-class of active galactic nuclei, which are characterized by intense non-thermal radiation dominating their spectral energy distribution (SED) from radio to very high energy (VHE) passbands. This radiation is strongly beamed, originating in a relativistic jet pointing towards the Earth at a small angle with the line of sight. From their optical spectra, blazars have been classified into two categories, viz. BL Lac objects (having featureless optical spectra) and flat spectrum radio quasars (FSRQs, with broad emission lines in their optical spectra).

Blazars have SEDs showing two peaks: one at low energies (from radio to optical and X-rays) due to synchrotron emission and another one at high energies (X- to γ-rays) due to inverse-Compton scattering of soft photons off the same relativistic electrons, which are responsible for the synchrotron radiation (Königl 1981).

The frequency at which the synchrotron emission peaks is used for further blazar classification. Padovani & Giommi (1995) classified BL Lac objects into low-energy peaked BL Lacs (LBLs) and high-energy peaked BL Lacs (HBLs). This classification was then generalized by Abdo et al. (2010) to include FSRQs: classes become low, intermediate, and high synchrotron peaked blazars. In particular, the low synchrotron peaked blazar class includes LBLs and virtually all FSRQs (Abdo et al. 2010).

Blazars show variability in all passbands, from radio to VHE γ-rays, as well as strong radio and optical polarization and superluminal motions. Since the relativistic jet is supposed to dominate the SED the strong variability could either be caused by internal processes within the jet itself or by accretion disk instabilities triggering changes in it. At radio frequencies the interstellar scintillation may also play a role (e.g. Heeschen et al. 1987).

Variability in blazars can be detected on diverse timescales. We can broadly divide blazar variability into three temporal classes, viz. intra-night variability (INV, or intra-day variability, or microvariability) on timescales from minutes to hours (e.g. Wagner & Witzel 1995; Kinman 1975; Clements et al. 2003; Rector & Perlman 2003; Xie et al. 2004), short-term variability (STV) on timescales from several days to a few months, and long-term variability (LTV, e.g. Gupta et al. 2008; Agarwal et al. 2017) on timescales from several months to many years; in the latter two classes the amplitude of the optical variability can often exceed even 5 mag (e.g. Bachev 2018). Variability studies are a powerful tool to probe the nature of emission processes occurring in blazars, their magnetic field geometry, dominant particle acceleration mechanisms, etc.

Along with flux variability studies, analyses of colour trends accompanying brightness changes have begun to be implemented in the last decade through the colour-magnitude diagrams (CMDs). Three types of CMD behaviour could be discerned: redder-when-brighter (RWB), bluer-when-brighter (BWB), and achromatic. The RWB chromatism is most frequently associated with FSRQs, for which the contribution of the accretion disk to the total emission could be significant. The BWB behaviour is thought to be related to processes associated to 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 (Villata et al. 2004; Papadakis et al. 2007). Finally, the achromatic behaviour is most frequently interpreted as being due to the variations of the Doppler factor, which are most likely explained in the framework of the geometric scenario (e.g. Villata et al. 2002).

PG 1553+113 was discovered by the Palomar-Green survey of UV-excess stellar objects (Green et al. 1986) and was classified as a BL Lac object due to its featureless optical spectrum (Miller & Green 1983) and significant optical variability (Miller et al. 1988). Furthermore, it was classified as an HBL, as its synchrotron emission peak falls in the UV and X-ray frequency range (Falomo & Treves 1990) and was detected at TeV energies (Aharonian et al. 2006). Due to its featureless optical spectrum and its host galaxy remaining undetected, the redshift of PG 1553+113 has remained uncertain (e.g. Treves et al. 2007). The recent detection of its putative galaxy group would set this object’s redshift at z = 0.433 (Torres Zafra et al. 2017; Johnson et al. 2019); we shall thus use this value throughout the paper.

PG 1553+113 is among the few blazars whose flux variability is claimed to be periodic. The first hint for periodicity comes from Ackermann et al. (2015) – they found a (2.18 ± 0.08) years period analysing the Fermi/LAT γ-ray light curve (LC) of the source. This period was further confirmed by Prokhorov & Moraghan (2017), Tavani et al. (2018), Sandrinelli et al. (2018), Covino et al. (2020), and Peñil et al. (2020)1 at γ-rays and by Cutini et al. (2016), Sandrinelli et al. (2018), and Covino et al. (2020) in the R-band. On the other hand, Covino et al. (2019) and Ait Benkhali et al. (2020) expressed some caution about the significance of the year-long period detected at γ-rays. Nilsson et al. (2018) did not find the claimed year-long period in their R-band LC. In addition, Lico et al. (2020) found no signs for clear periodic variability at the radio frequencies.

The caution expressed by Covino et al. (2019) and Ait Benkhali et al. (2020) is based on the low global significance of the periods they found. Sandrinelli et al. (2018) also reported modest global significance of the periods for both high energy and optical LCs of PG 1553+113. According to the authors, however, the presence of one and the same period in both γ-ray and R-band LCs means that this period could be real in spite of the modest individual global significance. Finally, Covino et al. (2020) found that the addition of a periodic signal to the modelling of the PG 1553+113 LCs improves their statistical description; the modelling itself is done by means of Gaussian processes.

In this paper, we present the results from intensive intra-night monitoring (INM) of PG 1553+113 along with the analysis of its historical VR-band LCs. The paper is organized as follows. We describe the observations and data reduction techniques in Sect. 2. In Sect. 3 we report the results of flux and colour variability of the source from intra-night to longer timescales, while, finally, Sect. 4 presents a discussion of our results followed by the conclusions.

2. Observations and data reduction

The data presented here were obtained on 76 nights over a period from January 2016 to August 2019. Observations were carried out in the optical BVRI-bands using nine different telescopes around the globe which are: 2.15 m Jorge Sahade telescope (JS, telescope A) and 60 cm Helen Sawyer Hogg telescope (HSH, telescope B), CASLEO, Argentina; 2.0 m Ritchey-Chrétien (RC, telescope C) and 50/70 cm Schmidt (telescope D) at the Rozhen National Astronomical Observatory (NAO), Bulgaria; 2.01 m RC Himalayan Chandra Telescope (HCT, telescope E) at Indian Astronomical Observatory, Hanle, India; 1.3 m JC Bhattacharya telescope (JCBT; telescope F) at the Vainu Bappu Observatory (VBO), India; 1.0 m RC telescope (telescope G) and 60 cm RC robotic telescope (telescope H) at TUBITAK National Observatory (TUG), Antalya, Turkey; 50 cm Cassegrain telescope (telescope I) located at the Astronomical Observatory of the Jagiellonian University in Krakow, Poland.

The main characteristics and technical details for telescopes A, B, C, D, and E can be found in Table 1 of Agarwal et al. (2019). In addition to these, for the present work, we have used four other optical telescopes (F, G, H, I), which we describe below.

The JCBT telescope located at the VBO, Tamil Nadu, India, is equipped with a 1k × 1k proEM CCD and a 2k × 4k regular CCD (the proEM CCD has a smaller field of view). As described later in this section, to perform differential photometry, suitable pairs of non-variable stars from the same observation frame are necessary. Therefore, we chose the 2k × 4k CCD for our photometric observations, whose central 2k × 2k region was used.

The 1.0 m RC telescope is located at the Bakirlitepe Mountain and is currently operated remotely from TUG in Antalya. The telescope is equipped with 4k × 4k CCD camera working at −90°C. Each observation was taken with 2 × 2 binning and consisted of several exposures through R filter and 1–2 sets of BVRI frames.

The 60 cm RC robotic telescope is also located at TUG. The observations with this telescope are object-based and fully automatic. Observers are allocated a maximum of 3000 s for each day on this telescope. After observations, the bias, dark, and flat corrected frames by the TALON pipeline are sent to observers. The robotic telescope is equipped with a 2k × 2k CCD camera and with 12 standard filters. As long as the weather is suitable, the observation was taken daily with 1 × 1 binning and 1–2 sets of BVRI frames.

In the period between May 7 and June 5, 2018, we performed observations of PG 1553+113 using the 50 cm Cassegrain telescope located at the Astronomical Observatory of the Jagiellonian University in Krakow, Poland. The telescope is equipped with an Apogee Alta U42 CCD and a set of broadband filters (Bessell 1990). The BVRI data were taken with 2 × 2 binning, with exposure times between 60 s and 100 s, depending on the weather conditions and the filter used. This camera accommodates an e2v, back-illuminated chip with 2048 px of 13 μm in size. The calibration images were taken every night by collecting series of bias, dark and sky flat-field frames.

Finally, the 2.0 m telescope of the Rozhen NAO was used in two additional observing configurations compared to that reported in Agarwal et al. (2019): the CCD camera VersArray:1300B attached directly to the RC focus and the two CCD cameras Andor iKon-L attached at the ports of the two-channel focal reducer FoReRo-2 (Jockers et al. 2000; Bonev 2011).

Additional information about the new telescopes and configurations can be found in Table 1. The complete observations log corresponding to the data collected during our monitoring campaign is given in Table 2.

Table 1.

Details about the telescopes and instruments used.

Table 2.

Log of photometric observations for the blazar PG 1553+113.

A detailed description of the data reduction procedure is given in Agarwal & Gupta (2015). To summarize, the data reduction of all CCD frames, including bias correction, flat-fielding, and cosmic-ray removal was performed using IRAF2 software. After correcting the CCD frames, the Dominion Astronomical Observatory Photometry (DAOPHOT II) software (Stetson 1987; Stetson et al. 1992) was used to perform aperture photometry which gave us the instrumental magnitudes of the source and all the stars in the CCD frames.

To build the instrumental differential LCs, we used four steady comparison stars in the frame (Raiteri et al. 2015). We finally selected those two stars that show smaller scatter of the star-star differential LC. As pointed out by Cellone et al. (2007), spurious INV could be detected if the comparison stars differ much in brightness with respect to the target. Therefore, we selected those comparison stars which had magnitude and colour close to that of the target blazar. The dispersion in the star-star LCs was checked for every night and various apertures, starting from a value equal to the full width at half maximum to four times this value. Finally, the aperture with minimum dispersion in the star-star LCs was selected. An aperture radius of 4″ was thus used. The calibration of instrumental magnitudes into the standard system was made using magnitudes listed in Raiteri et al. (2015).

3. Results

3.1. Intra-night variability

The LCs derived during the 28 nights of INM are shown in Fig. A.1. The date of observation and the telescope used are indicated in each plot. We show only LCs having a minimum of 10 data points acquired during the particular night.

Considering the modest number of data points in a single filter, we used the C-test and the F-test, both set at a 99.5% confidence level to assign variability status to the intra-night LCs. The tests are described in detail in Agarwal et al. (2019). Even if the comparison and control stars were appropriately selected, we applied the Γ scaling factor (Howell et al. 1988) to equalize the target and control error distributions, since this procedure has been shown to give the most reliable results (e.g. Zibecchi et al. 2017, 2020).

After computing the parameters of each test with the original choice of stars (C1, F1), we then recalculated them with comparison and control stars interchanged (C2, F2) as a check for any issue in the LC. We have designated a light curve to be variable (Var) if both tests rejected the null hypothesis (i.e. non-variable LC) at a 99.5% confidence level and non-variable (NV) if at least one of the tests failed to reject the null hypothesis. The results from the INV tests are summarized in Table A.1.

3.2. Duty cycle

We found no significant variability during 28 nights of INM. Thus, our INM campaign resulted in zero duty cycle (DC; for details about the DC see Romero et al. 1999). In order to improve this estimate, we searched the literature for other cases of INM of PG 1553+113.

The first INM of PG 1553+113 to our knowledge was reported by Stalin et al. (2005) – they applied the C-test to the R-band LCs (total duration of about 10 h) and found the source to be variable in one night out of two. Another two nights of INM about 13 h long were presented by Osterman et al. (2006). They found no significant INV, but they did not present statistical tests. To quantify the conclusion of Osterman et al. (2006), we used the C-test in the form C = σ/⟨e⟩, where σ is the standard deviation of the source LC and ⟨e⟩ the mean uncertainty of the source photometric measurements. We applied this test to the data listed in Table 4 of Osterman et al. (2006) and found C = 1.52 and C = 1.28 for the respective nights, which means that PG 1553+113 was non-variable. Andruchow et al. (2007, 2011) reported the results from the INM they carried out during four nights in the VR-bands (April 21 to 24, 2007) and five nights in the BR-bands (April 21 to 25, 2009), respectively. Applying the C-test, the authors found no significant variability. Gopal-Krishna et al. (2011) detected INV during three nights out of three using C- and F-tests (a total monitoring duration of 16 h in the R-band). Gaur et al. (2012) found no INV during six nights of monitoring in the BR-bands (C- and F-tests were applied). Gupta et al. (2016) presented the results from the INM during seven nights in the R-band (a total monitoring duration of about 26 h). The authors found the source to be variable in one night, non-variable in another, and probably variable in the remaining nights (F- and χ2 tests were applied). We, however, should point out that in the latter nights the variability amplitude seems to be quite close to the magnitude uncertainties. In such cases the usage of the C-test could be more appropriate (see Zibecchi et al. 2017, 2020). Pandey et al. (2019) monitored PG 1553+113 for eight nights in the VR-bands for 2–4 h each night. Employing enhanced F- and nested ANOVA tests, the authors found the source to vary on intra-night timescales for three nights. Finally, Pasierb et al. (2020) found no INV in the BVR-bands during a single night of monitoring (duration of about 3.7 h, F-test employed).

It is worth mentioning the study of Meng et al. (2018), which is focused mainly on the LTV of a sample of blazars including PG 1553+113. For this source the observations were performed for 28 nights through three intermediate-band filters. The authors did not find INV during each of the observing nights using four statistical tests. In some of the nights, however, the observing time span is insufficient for the INV testing – less than an hour. So, we re-analysed the PG 1553+113 data selecting those nights for which the observing time span is Δt >  2 h and the number of the acquired data points is N >  3. In this way we got a total of eight nights with 2.2 h ≤ Δt ≤ 4.6 h and 7 ≤ N ≤ 11. We then applied the C-test and assumed the source to be variable for the given night if C ≥ 2.576, at least in two of the passbands. We found this condition not met for the selected nights.

The literature data, combined with ours, resulted in a total of 74 nights of PG 1553+113 INM. Despite the large number of nights of INM the blazar PG 1553+113 was found to show INV during only eight nights; during another five nights, the source was classified as probably variable. Because of the lack of information about the total monitoring duration for some of the literature data, we computed the DC simply as NINV/NINM, where NINV is the number of the nights with INV detected and NINM the total number of nights of INM. The so computed DC is 10.8% if the probably variable cases are considered non-variable and 17.6% – otherwise. Our result agrees with the DC of 11.6% reported by Andruchow et al. (2014) for a sample of 6 HBLs.

3.3. Long-term variability

We based the study of the optical behaviour of PG 1553+113 at long-term timescales on our observations along 76 nights from 2016 to 2019. In the course of this multi-wavelength campaign, we acquired a total of ∼3350 data points in BVRI-bands.

In order to remove effects due to dense intra-night sampling during some of the nights, we performed nightly binning of our data. Firstly, we discarded some unreliable data points and then the magnitudes obtained during each given night were weight-averaged. Regarding the uncertainties of the weight-averaged magnitudes, we adopted the following conservative estimate: we took the larger between (i) the uncertainty of the weighted mean (which accounts only for the individual measurement uncertainties) and (ii) the weighted standard deviation about the weighted mean (which accounts mainly for the scatter of the individual data points with respect to the weighted mean). Finally, we took the median of the corresponding MJDs. The data points that give rise to uncertainties of the weight-averaged magnitudes larger than or equal to 0.1 mag were removed and the binning was done again.

Our long-term BVRI-band LCs are presented in Fig. 1. Both statistical tests applied to them revealed that LTV is highly significant at all four photometric passbands.

thumbnail Fig. 1.

Long-term LC for PG 1553+113 covering the entire monitoring duration. Different colours and symbols denote data in different passbands as indicated in the plot. In each passband the magenta solid lines connect the data points obtained after nightly binning of the respective LC.

The most remarkable feature of the LCs is the flare that was detected in all passbands on April 14, 2019 (MJD ∼8588.4) with the R magnitude reaching 13.212, which seems to be the highest level recorded for PG 1553+113 (see also the Erratum3 to Jankowsky & Wagner 2019). On the other hand, the source attained the faintest state on January 12, 2016 (MJD ∼7400.7) with an R magnitude of 14.444. We observed a trend of a larger variability amplitude at higher frequencies on a long-term basis. The LTV characteristics are given in Table 3. We also assembled the historical VR LCs of PG 1553+113 covering a time interval from 2005 to 2019 (Fig. B.1).

Table 3.

Long-term variability characteristics.

3.4. Colour–magnitude diagram

Multi-band monitoring campaigns in the optical BVRI-bands help us in studying spectral changes on various timescales. The study of CMDs could shed some light on the physical processes causing the blazar flux variations. In this context, we firstly considered our BVRI data set and then the historical VR one. Generally, we calculated colour indices by coupling data taken during the same night. The coupling within the night is justified as the DC of the INV of PG 1553+113 is low (see Sect. 3.2).

We plot in Fig. 2 the colour index with the largest frequency base, viz. B − I, against the R magnitude – we choose the R-band following the prescription of Massaro & Trèvese (1996). The presented CMD – from our data only – shows no clear correlation. Instead, one could trace a RWB chromatism when the source is faint and a BWB one when the source is bright. Similar behaviour was reported by Ikejiri et al. (2011) using V − J vs. V CMD; the source was observed in the VJKS-bands from 2008 to 2010. The number of the data points used to build the above mentioned CMDs is relatively small, so, any conclusions based on their analysis should be considered with caution. A denser sampling on B − R vs. R CMD for PG 1553+113 was achieved by Wierzcholska et al. (2015); the source was observed in BR-bands from 2007 to 2012. The authors reported a statistically non-significant RWB trend considering the whole data set and both RWB and BWB trends on shorter time intervals.

thumbnail Fig. 2.

Colour-magnitude diagram for our data only.

Let us analyse the historical LCs constructed by us from the data sets presented in Appendix B. Although the frequency base is not so large, the time span (from 2005 to 2019) and the total amplitude of variations are sufficient to look for possible CMD peculiarities. We plot in the upper panel of Fig. 3 the CMD for the historical data set: the colour indices are plotted against the mean of the V and R magnitudes. This mean could be considered as representative of the intermediate-band magnitude. The overall trend is BWB with the following characteristics of the corresponding linear least-squares fit: a slope of 0.030 ± 0.004; χ df 2 =2.8 $ \chi^2_{\rm df}=2.8 $; a linear Pearson correlation coefficient rCMD = 0.16 with corresponding p-value of 0.002; a standard deviation about the fitted line of 0.037 mag. In the following, we shall consider a CMD trend (i) significant at 99% confidence level if rCMD ≥ 0.5 and the null hypothesis probability is p ≤ 0.01 (e.g. Gupta et al. 2016) and (ii) achromatic if the absolute value of the slope is less than the corresponding fitted uncertainty. Therefore, the so obtained BWB trend is non-significant. In addition, the presented CMD shows a large scatter about the fitted line. This scatter is caused mainly by the data sets we used: they are obtained using various telescope+filter+CCD combinations and the differences in the resulting response curves were not taken into account during the inter-calibration process. In an attempt to minimize this effect we plot the CMD for our, Steward4, and Kanata VR data sets (Fig. 3, lower panel); the other sets are single-band or had to be transformed to the VR-bands, so, their contribution to the scatter is more substantial. Now, the BWB trend becomes clearer and has the following linear fit characteristics: a slope of 0.041 ± 0.004; χ df 2 =3.6 $ \chi^2_{\rm df}=3.6 $; a linear Pearson correlation coefficient rCMD = 0.35 with corresponding p-value less than 10−5; a standard deviation about the fitted line of 0.027 mag. This weak trend, however, is non-significant following the above criterion. The CMD shows three segments at different mean brightness levels in which local RWB trends could be followed. Similar complicated behaviour – an overall BWB chromatism with local RWB trends – could be traced on the B − R vs. R CMD of Raiteri et al. (2015) which has a smaller time span (about five months) but an excellent time coverage.

thumbnail Fig. 3.

Colour-magnitude diagram for the historical data set (upper panel, 373 data points plotted) and for our, Steward, and Kanata data sets (lower panel, 274 data points plotted). The linear least-squares fits are overplotted. For both panels the median uncertainties of the magnitudes and colour indices are 0.018 mag and 0.036 mag, respectively.

Furthermore, we plot in Fig. 4 the temporal variations of the V − R colour index compared to the scaled R-band LC. The comparison suggests no clear correlation.

thumbnail Fig. 4.

Temporal evolution of the V − R colour index (black dots). The scaled R-band LC is overplotted for comparison (magenta line, the flux increases along +y axis). The horizontal dashed line marks the weighted mean V − R value of (0.339 ± 0.002) mag with a weighted standard deviation of 0.030 mag about the weighted mean. The insert shows the 2019 flare in details – the maximum of the spectral hardness (the minimum of V − R) precedes the flare by about 10 days (see Sect. 3.5 for details).

To study the colour index behaviour on small timescales, we divided the LCs into separate segments. Considering Fig. B.1 one sees 5 distinct minima that divide the LCs into 6 flares; in addition, the 6th flare can be divided into a pre-flare (around MJD of 8000) and a flare5. Eventually, we divided the LCs into seven segments and built a CMD for each of them. The most significant trend (rCMD = 0.49 and p = 0.0007) was found for the brightest flare – this flare will be analysed in details in Sect. 3.5. Regarding the other segments, all of them show either achromatism or non-significant chromatism. In particular, the pre-flare could be considered achromatic – the absolute value of the CMD slope (viz. −0.007 ± 0.011) is less than its uncertainty.

The inspection of Fig. B.1 suggests at least two variable components – a long-term one, responsible for the quasi-periodic behaviour of the source, and a short-term one producing the flux variations onto the top of the long-term component, that is, F(tot)(t) = F(long)(t)+F(short)(t), where F means flux. Therefore, to study the CMD behaviour on small scales, we have to remove the possible impact on their CMDs of the long-term variations. We shall assume that the long-term component is achromatic on the base of the above analysis.

To separate the components presumed, we used the approach of Villata et al. (2002). As a first step, we binned the VR-band LCs over a time interval of 10 days. Then, we transformed magnitudes into fluxes6 and fitted a cubic spline to the binned LCs. The mean flux ratio of splines over the time, F V ( spl ) ( t ) / F R ( spl ) ( t ) t $ \langle F^{\,\rm (spl)}_V(t)/F^{\,\rm (spl)}_R(t)\rangle_t $, is 0.89 ± 0.06 (s.d.), which corresponds to a two-point spectral index α VR ( spl ) = 0.72 ± 0.42 $ \alpha_{VR}^{\mathrm{(spl)}}=0.72 \pm 0.42 $ – these are the flux ratio and spectral index of the long-term component (or the base level). The spline fits are very similar, so, we averaged them for further analysis. After that we subtracted from the VR-band LCs the mean spline, shifted in such a way, that the condition [F(t)−e(t)] − F(spl)(t) > 0 is fulfilled for each MJD and both passbands; here e is the uncertainty of the flux level F. In this way, we corrected the LCs for the variable long-term component. The spline fits and the adopted base levels are shown in Fig. 5.

thumbnail Fig. 5.

Historical VR-band LCs (dots). The cubic splines, fitted through the 10-day binned LCs, are overplotted (V-band – green line; R-band – red line) along with the shifted mean spline representing the base levels for the VR-bands (black lines).

The next assumption, based solely on the achromaticity, is that the base level variability is a result of the Doppler factor change. We recall that the Doppler factor changes if the viewing angle varies. To account for this relativistic boosting we divided each flux data point by the ratio, C(t), between the spline value at that time and the spline global minimum value, C ( t ) = F ( spl ) ( t ) / F min ( spl ) $ C(t) = F^{\mathrm{(spl)}}(t)/F^{\mathrm{(spl)}}_{\mathrm{min}} $.

The resulting “corrected” LCs are presented in Fig. 6. The variation amplitudes are comparable over the entire period covered by the historical data set. This should mean that we see the variations of the short-term component alone. The corresponding CMD, however, still does not show any clear chromatism (Fig. 7). The characteristics of the linear fit presented in Fig. 7 are: a slope of 0.120 ± 0.007; χ df 2 =5.3 $ \chi^2_{\rm df}=5.3 $; a linear Pearson correlation coefficient rCMD = 0.22 with corresponding p-value of 0.003; a standard deviation about the fitted line of 0.131 mJy.

thumbnail Fig. 6.

Historical VR-band LCs corrected for the long-term variations.

thumbnail Fig. 7.

Flux ratio FV/FR plotted against (FV + FR)/2. Shown are our, Steward, and Kanata data sets corrected for the long-term variations. The best fit linear model to the data is overplotted. The median uncertainties of the fluxes and flux ratios are 0.082 mJy and 0.087 mJy, respectively.

The scatter about the fitted line now is larger compared to the pre-corrected CMD scatter of 0.022 mJy. This could be attributed mainly to the unknown true base level and shape. Regarding the shape, we implicitly assumed that it could be represented by the cubic spline fitted to the binned LCs. Thus, by correcting as explained above, we introduce uncertainties, firstly, by the base level subtraction itself and, secondly, by the Doppler boosting correction based on this level. The other source of uncertainties is the quality of the cubic spline interpolation, which depends on the LC sampling – both in terms of the total number of data points and the degree of their clustering. We, however, can consider this source of less importance regarding the CMD scatter. Despite the large scatter, now the CMD shows no significant sub-structure compared to the pre-corrected CMD (see Fig. 3).

Regarding the short time span segments, we divided the historical LCs into, the situation does not change. All segments but the pre-flare and flare ones show either achromatism or non-significant chromatism. Now, the pre-flare shows non-significant BWB chromatism (after removing 3 deviant data points) and the flare BWB trend gets more significant – rCMD = 0.63 and p <  10−5.

The above analysis shows that on both long and short timescales, the variability of PG 1553+113 has no significant chromatism (except the 2019 flare). Another notable result is that the corrected LCs show signs for quasi-periodicity (see Sect. 3.6).

Given the achromaticity alone, we cannot distinguish whether it reflects variations caused by geometric effects or by processes intrinsic to the jet. However, the observed periodicities, discussed in Sect. 3.6, could be considered in support of the geometric origin of the PG 1553+113 variability (e.g. Sobacchi et al. 2017).

3.5. The 2019 flare

3.5.1. Colour hysteresis and time lag

During April 2019, PG 1553+113 reached its historical brightness maximum showing a prominent flare. The flare shows some sub-structure – there is a flux increase (sub-flares?) around MJD 8550 and MJD 8610 (the latter is best traced in the R-band). The LCs sampling, however, is not good enough to derive further conclusions about these sub-structures. The detailed inspection of the flux and colour (or spectral) index behaviour during the flare revealed a lag of the flux changes behind the spectral ones (see Fig. 4). The spectrum got its hardest value (data point #8 in Fig. 8) before the flux maximum was reached (data point #19) and then continuously softened till the end of the flare. The corresponding CMD shows a clockwise spectral hysteresis loop – the spectrum hardens (i.e. gets flatter ≡ shows BWB chromatism) as the flux rises and softens (i.e. gets steeper ≡ shows RWB chromatism) as the flux declines (Fig. 8). This kind of hysteresis is typically observed at X-rays in HBLs (e.g. Takahashi et al. 1996) and signals to the existence of a time lag, 𝒯: the flux variations at a high frequency (V-band) lead those at a low frequency (R-band; this is the so-called “soft” lag7). The detected hysteresis loop could also be traced on the CMD built using only our and Steward data sets as well as on the CMD built using the historical VR data sets after their correction for the long-term variations. So, we shall consider this hysteresis as real rather than an artefact from the heterogeneous data used.

thumbnail Fig. 8.

Colour-magnitude diagram for the 2019 flare. The corresponding spectral index: α = [(V − R)−(AV − AR)−0.186]/0.07, is also indicated. Here AV and AR stand for the Galactic extinction in the VR-bands, respectively (see Sect. 3.4). The clockwise spectral hysteresis can be clearly seen – the arrows are used to guide the eye about the hysteresis loop direction. The colours and numbers are sequential in time: violet (#1)→ red (#44).

To search for the expected time lag we utilized both the discrete correlation function (DCF, Edelson & Krolik 1988) and the z-transformed DCF (ZDCF, Alexander 1997, 2014). We cross-correlated the V- and R-band LCs transformed into fluxes.

The DCF is commonly used to search for correlation in two unevenly sampled time series. We did not take into account the measurement uncertainties in the construction of the DCF following White & Peterson (1994). In our particular case the discrete correlations were binned using a time lag bin of width Δ𝒯 = 1.5 days and a Gaussian weighting scheme was applied, so that the higher importance is assigned to the unbinned values closer to the bin centre. In our notation the positive lag means that the variations in the V-band lead those in the R-band. The time lag was determined using both the Gaussian fitting and the centroiding methods; we took into account only those lags, for which the correlation coefficient is larger than 0.5. We found 𝒯gauss = (5.76 ± 1.48) days using the Gaussian fit and 𝒯cent = 3.71 days using the weighted centroid. The DCF is shown in Fig. 9. For more reliable estimation of the time lag and its uncertainty, we applied the flux randomization/random subset selection method (FR/RSS, Peterson et al. 1998) based on Monte Carlo simulations. The data points counted more than once during the RSS process were rejected. For each FR/RSS realization the time lag was found by means of the Gaussian fit. We ran a total of 1000 cycles and the resulted time lag values were used to build the cross-correlation peaks distribution (CCPD, Maoz & Netzer 1989). The CCPD is shown in Fig. 10 and the derived time lag is T CCPD = 4 . 86 2.62 + 2.50 $ \mathcal{T}_{\mathrm{CCPD}}=4.86^{+2.50}_{-2.62} $ days. The quoted time lag value is the 50th percentile (the median) of the CCPD, while the 16th and 85th percentiles serve as the 1σ time lag uncertainties.

thumbnail Fig. 9.

Discrete correlation function. The fitted Gaussian is not shown for the sake of clarity.

thumbnail Fig. 10.

Cross-correlation peak distribution for DCF.

The application of ZDCF gave similar results (Fig. 11). The corresponding estimates of the time lag are 𝒯gauss = (3.71 ± 0.88) days and 𝒯cent = 3.85 days. It is worth noting the better agreement between these estimates compared to the DCF ones. So, we detected a significant positive time lag in agreement with the clockwise spectral hysteresis.

thumbnail Fig. 11.

Same as in Fig. 9, but for ZDCF.

3.5.2. Flare fit

In order to quantify the flare, we fitted it with an exponential law (Valtaoja et al. 1999) plus a linear base level:

F ( t ) = { a + b t + exp [ + ( t t max ) / T ris ] if t t max a + b t + exp [ ( t t max ) / T dec ] if t t max , $$ F(t) = {\left\{ \begin{array}{ll} a+b\,t+\exp [+(t-t_{\rm max})/{\mathcal{T} }_{\rm ris}]&\quad \mathrm{if\,} t \le t_{\rm max} \\ a+b\,t+\exp [-(t-t_{\rm max})/{\mathcal{T} }_{\rm dec}]&\quad \mathrm{if\,} t \ge t_{\rm max} \end{array}\right.}, $$

where (a, b) are the parameters of the local linear base level, tmax the epoch of the flare maximum and (𝒯ris, 𝒯dec) the flare rise and decay e-folding timescales, respectively. The data points belonging to the possible sub-flare at MJD ∼ 8550 and the data point at the R-band maximum were excluded from the fit. The fitted maximum position along with the rise and decay e-folding timescales are listed in Table 4. The resulting fits are shown in Fig. 12.

thumbnail Fig. 12.

Approximation of the 2019 flare with an exponential law plus a local linear base level. The possible sub-flares could be discerned around MJD = 8550 and MJD = 8610.

Table 4.

Flare best fit parameters – observer-frame.

The fits confirmed the time lag found by the cross-correlation techniques – the flare V-band LC leads the R-band one by 𝒯fit = (2.6 ± 0.2) days. This lag equals to the cross-correlation one to within the uncertainties quoted.

3.6. Search for periodicity

We addressed the PG 1553+113 periodicity issue using our historical VR-band LCs and applying three different methods. The periodicity presence is obvious if we consider the spline fits to the historical LCs (see Fig. 5); we recall that the spline fits represent the long-term variability component. To quantify the periodicity, we fitted a sine function of the form F(t) ∝ sin (2πt/P + Φ) to the corresponding spline; here P and Φ stand for the period and phase, respectively. The resulting fits are shown in Fig. 13. We found periods of (807 ± 20) days and (812 ± 30) days considering the V- and R-band LCs, respectively. The corresponding mean differences between the maxima of the splines and of the fitted sine waves are (9 ± 20) days and (17 ± 30) days. The quoted standard errors of the mean were used as the uncertainties of the corresponding periods.

thumbnail Fig. 13.

Sine function fits to the VR splines: V-band – lower green curve with black dots overplotted, R-band – upper red curve with black dots overplotted, sine fit – black curves of constant amplitude.

Next, we applied the classical Lomb-Scargle periodogram (LSP) analysis, which was proposed by Lomb (1976) and extended by Scargle (1982). The LSP has been generalized for more practical use by Press & Rybicki (1989). More details are outlined in VanderPlas (2018) and references therein. The statistical significance of the periods, found by means of LSP and REDFIT (see below) methods, is evaluated locally (e.g. Ait Benkhali et al. 2020).

The periodograms are shown in Figs. 14 and 15. The maxima of the most significant peaks correspond to periods of (803 ± 70) days and (801 ± 62) days, for the V- and R-band LCs, respectively. The quoted uncertainties represent the half-width at the half-power maximum derived by means of Gaussian fit; we ignore the possible asymmetry of the fitted peaks. These periods agree to within the uncertainties with the ones obtained by applying the sine fits to the splines. Because of this, we shall consider these periods as real. There are, however, a number of additional, less significant peaks – their existence is most probably caused either by the departure of the shape of the folded LCs from a sinusoid (see the inserts of Figs. 14 and 15) or by aliasing due to the finite monitoring time span and irregular sampling. So, we shall not consider these peaks further.

thumbnail Fig. 14.

Lomb-Scargle periodogram for the historical V-band LC. The dotted line marks the level corresponding to a maximum peak false alarm probability of 1% under the assumption of white noise. The insert shows the LC folded with the most significant period; the red line denotes the fitted sine wave.

thumbnail Fig. 15.

Same as in Fig. 14, but for the R-band.

Inspecting Fig. 14 one can see that the amplitude of the power spectrum between the major peaks, that is, the noise-induced power spectrum decreases with decreasing period – this is a signature of the presence of red noise; generally, the red noise arises from stochastic processes and is frequency-dependent. To test if the power spectrum peaks are significant against the red noise background, we utilized the improved translation of the REDFIT programme (Schulz & Mudelsee 2002) run within the R environment8. In this programme an auto-regressive process of first-order (AR1) is used to approximate the red noise. Before actual computations, we checked whatever the AR1 model is appropriate to characterize our historical VR-band LCs; the non-parametric runs test indicates that both the VR spectra are consistent with the AR1 model.

We ran REDFIT using the oversampling factor of 4, two overlapping by 50% segments, and Welch spectral window. The results are shown in Fig. 16. The peaks exceeding the 99% significance level correspond to periods of (772 ± 85) days and (817 ± 107) days, for the V- and R-band, respectively. The periods agree to within the uncertainties with the above estimates; the uncertainties are larger owing to the broader peaks.

thumbnail Fig. 16.

Output from the REDFIT programme. The black line shows the bias-corrected power spectrum, the blue dashed line marks the theoretical red noise spectrum and the red line marks the 99% local significance level derived by means of Monte Carlo simulations.

There is another peak above the 99% significance level corresponding to a period of (248 ± 12) days (Fig. 16, R-band9). A hint for quasi-periodicity on such timescales could be found inspecting Fig. 6, especially the R-band LC. A similar period was found by Sandrinelli et al. (2018): (250 ± 60) days. In addition, the most significant rest-frame period, reported by Nilsson et al. (2018) – 174 days – is in agreement with the above periods, transformed to the rest-frame.

To check further the existence of a possible period in the range 200–250 days, we apply the LSP method and the REDFIT programme to the LCs corrected for the long-term variations. The results from the LSP application are presented in Figs. 17 and 18. The most significant peaks correspond to periods of (208 ± 6) days and (212 ± 6) days, for the VR-band LCs, respectively. The REDFIT found the following periods: (211 ± 10) days and (210 ± 8) days, for the VR-band LCs, respectively (Fig. 19). The non-parametric runs test, however, indicated that the spectrum of the corrected V-band LC is inconsistent with the AR1 model, that is, non-AR1 components are present; the REDFIT assumes that the noise in a time series can be approximated by an AR1 process.

thumbnail Fig. 17.

Same as in Fig. 14, but for the historical V-band LC corrected for the long-term variations.

thumbnail Fig. 18.

Same as in Fig. 17, but for the R-band.

thumbnail Fig. 19.

Same as in Fig. 16, but for the historical VR-band LCs corrected for the long-term variations.

In Table 5 we summarize the periods found. The local significance of the periods – except those, obtained by means of the sine wave fit – is greater than 2.5σ. The similar local significance of the same year-long period was found by Sandrinelli et al. (2018) using the AR1 model.

Table 5.

Summary of the periods found in the historical VR-band LCs of PG 1553+113.

3.7. Spectral energy distribution

To further understand spectral variations in the source, we studied its multi-band optical SEDs. In Fig. 20, we have generated optical SEDs of PG 1553+113 considering all dates when we have quasi-simultaneous observations in at least three optical passbands. To display SED variations during our observing campaign, we selected different states of the source as follows: an outburst state (March 13 and April 4, 2019), intermediate state (June 11, 2016, May 10, 2018, and February 28, 2019) and low state (January 12 and August 11, 2016). These seven SEDs display an increasing trend, which could hint towards the fact that the synchrotron peak is located above 1015 Hz. This is in accordance with the HBL classification. Therefore, based upon our optical study, we can conclude that the source behaves more like an HBL rather than an LBL. However, these optical SEDs cover a narrow frequency range, therefore further multi-wavelength observations will offer a more powerful diagnosis of the true class of the source.

thumbnail Fig. 20.

Optical SEDs of PG 1553+113 for selected brightness levels (see text).

The SEDs do not show any significant shape change as the mean flux level varies. This lack of variations in the SED shapes allows us to build the mean SED over the time interval covered by our monitoring campaign using the weighted mean magnitudes listed in Table 3. To get the mean spectral index, we fitted a straight line of the form log(Fν) = − α log(ν)+const to the mean SED. We got α = 0.89 ± 0.06 with χ df 2 =1.1 $ \chi^2_{\rm df}=1.1 $, a result consistent with those of Falomo et al. (1994) and Pandey et al. (2019).

4. Discussion and conclusions

In this paper, we have presented the results from the optical monitoring of the blazar PG 1553+113 from 2016 to 2019 on intra-night and long-term timescales. We have also constructed and analysed the historical 2005–2019 VR-band LCs of the blazar.

In the course of our monitoring campaign, we recorded in April 2019 the brightest state of PG 1553+113 for the period from 2005 to 2019, covered by the historical LC (see Figs. 1 and B.1): R ≃ 13.2 mag. In general, the variations of the PG 1553+113 flux on both long and short timescales are achromatic, except the 2019 flare. The spectral index that characterizes the mean SED of the blazar over the period from 2005 to 2019 is α = 0.89 ± 0.06. We shall discuss below some constraints on the jet parameters imposed by our observing results.

4.1. Intra-night variability

From the analysis of our and literature INM data during a total of 74 nights, we found a low value of the PG 1553+113 DC: 10.8% if the probably variable cases are considered non-variable and 17.6% – otherwise. This result agrees with the low DC reported for HBLs (e.g. Andruchow et al. 2014). The low value of the DC imposes some constraints on the physical jet parameters.

Let us assume that the INV in blazars is associated with interactions of shocks with small scale structures – like density inhomogeneities or eddies – in the otherwise steady jet flow. Then the low DC value of HBLs could be explained with the strong magnetic field present in this sub-class of blazars (e.g. Sambruna et al. 1996). Romero (1995) found that the axial magnetic field halts the build-up of Kelvin-Helmholtz instabilities in the sub-parsec to parsec scale jets when the field value exceeds the critical value defined as:

B c = [ 4 π n e m e c 2 ( Γ 2 1 ) ] 1 / 2 Γ 1 , $$ \begin{aligned} {\mathcal{B} }_{\rm c} = \left[4 \pi n_{\rm e} m_{\rm e} c^2 \left(\Gamma ^{\,2} - 1\right)\right]^{1/2}\Gamma ^{\,-1}, \end{aligned} $$(1)

where ne is the local electron density, me the rest mass of the electron, and Γ the bulk Lorentz factor of the flow. Therefore, we have that, in case of HBLs, there will be no instabilities when ℬ > ℬc, which in turn reduces the incidence of rapid variability in their optical LCs.

4.2. The 2019 flare

The 2019 flare marks the brightest state of PG 1553+113 for the period from 2005 to 2019. The analysis of the flare VR-band LCs revealed a clockwise hysteresis loop and a soft lag, which means that the synchrotron cooling of accelerated electrons is the dominant emission mechanism of the flare (Kirk et al. 1998).

The flare has a slight asymmetry, with the rising part being steeper. This asymmetry, however, is questionable because of the putative sub-flare at MJD ∼8610, which is not well resolved in order to be accounted for in the proper way. So, for further analysis, we shall assume that the 2019 flare is symmetric with the decay timescale equal to the rise one. This means that the cooling time is smaller than or comparable to the light crossing time, tcool ≲ ℛ/c.

4.2.1. Radius of the emitting region

Under the assumption that the injection time is less than the light crossing time, the rising part of the flare LC could be used to constrain the radius of the emitting region: ℛ ≲ c𝒯risδ/(1 + z), where δ is the Doppler factor and z the redshift. Using a rise timescale of 24.9 days, which is the weighted mean over the V- and R-band rise timescales, we got a radius of ℛ ≲ 4.5 × 1016 δ cm.

Regarding the Doppler factor of PG 1553+113, there exist two sets of values. The first set is based on the multi-wavelength SED modelling – δ varies from 20 (Tavecchio et al. 2010) to 40 (Aleksić et al. 2015). The second, single-value set is based on the Very Long Baseline Array observations – δ is equal to 1.4 (Lico et al. 2020). This discrepancy is typical for TeV HBLs and reflects the so called “Doppler crisis” (Tavecchio 2006). The “Doppler crisis” can be resolved if one assumes the presence of the velocity gradient either in the axial or in the transverse direction (for details see, e.g. Piner et al. 2018, and references therein). This assumption could reduce the Doppler factor that characterizes the one-zone synchrotron self-Compton models.

The typical radius of the emitting region, used in the SED modelling, is ℛ ≲ 1017 cm (e.g. Aleksić et al. 2015; Banerjee et al. 2019). Therefore, the Doppler factor must be δ ≲ 3 given our estimate of the weighted mean rise timescale. The Doppler factor of 20–40, as derived from the SED modelling, would result in a too large radius, ℛ ≃ (9−18) × 1017 cm.

4.2.2. Magnetic field strength and electron energies in the emitting region

Assuming tcool ≲ ℛ/c, we could estimate the lower limit of the magnetic field strength in the emitting region using the rise timescales listed in Table 4 (Ghisellini et al. 1997):

B δ 1 / 3 0.67 ν 15 1 / 3 T ris 2 / 3 ( 1 + U S / U B ) 2 / 3 ( 1 + z ) 1 / 3 , $$ \begin{aligned} {\mathcal{B} }\delta ^{1/3}\gtrsim 0.67\,\nu _{15}^{-1/3}\,{\mathcal{T} }_{\rm ris}^{\,-2/3}\,(1+U_{\rm S}/U_{\rm B})^{\,-2/3}\,(1+z)^{1/3}, \end{aligned} $$(2)

where ℬ is in units of Gauss, 𝒯ris is in units of days, ν15 is the observed frequency (in units of 1015 Hz), US the radiation energy density of the synchrotron emission, and UB the magnetic energy density; the energy densities are measured in the comoving frame and we assume that US = UB. We got ℬδ1/3 ≳ 0.071 G for the V-band and ℬδ1/3 ≳ 0.072 G for the R-band. So, if we use Doppler factors of 3 and 20, then the respective lower limits are 0.05 G and 0.03 G and they are the same for both VR-bands.

An independent way to estimate the magnetic field strength is to use the soft time lag and the following expression (e.g. Papadakis et al. 2003):

B δ 1 / 3 300 ( 1 + z ν R ) 1 / 3 [ 1 ( ν R / ν V ) 1 / 2 T ] 2 / 3 , $$ \begin{aligned} {\mathcal{B} }\delta ^{1/3}\simeq 300\,\left(\frac{1+z}{\nu _R}\right)^{1/3}\,\left[\frac{1-(\nu _R/\nu _V)^{1/2}}{{\mathcal{T} }}\right]^{2/3}, \end{aligned} $$(3)

where ℬ is in units of Gauss, νV and νR are the frequencies of the VR-bands (in units of 1017 Hz) and 𝒯 the soft time lag (in units of seconds). Here it is assumed that the lag resulted from synchrotron cooling of high-energy electrons to lower energies (Chiappetti et al. 1999). We got ℬδ1/3 ≃ 0.066 G if 𝒯 = 𝒯CCPD and ℬδ1/3 ≃ 0.100 G if 𝒯 = 𝒯fit. If we use Doppler factors of 3 and 20 then the respective values are 0.05 G and 0.02 G using 𝒯CCPD and 0.07 G and 0.04 G using 𝒯fit. The magnetic field strengths estimated in both ways agree well to each other as well as with the values obtained by Tavecchio et al. (2010) for TeV blazars by means of SED modelling.

Finally, we obtained the lower limit of the energy, γobs (in units of mec2), of the electrons emitting at the observed frequencies corresponding to VR-bands (Ghisellini et al. 1997):

γ obs δ 1 / 3 2 × 10 4 ν 15 2 / 3 [ T ris ( 1 + U S / U B ) ( 1 + z ) ] 1 / 3 . $$ \begin{aligned} \gamma _{\rm obs}\delta ^{1/3}\lesssim 2\times 10^4\,\nu _{15}^{2/3}\,[{\mathcal{T} }_{\rm ris}\,(1+U_{\rm S}/U_{\rm B})\,(1+z)]^{1/3}. \end{aligned} $$(4)

We got γobsδ1/3 ≲ 5.5 × 104 for the V-band and γobsδ1/3 ≲ 5.0 × 104 for the R-band. So, if we use Doppler factors of 3 and 20 then the respective upper limits are 3.8 × 104 and 2.0 × 104 for the V-band and 3.5 × 104 and 1.8 × 104 for the R-band.

4.3. Periodicity analysis

Our periodicity analysis of the historical VR-band LCs using three different methods gives a median (over the methods and VR-bands) period of [2.21 ± 0.04 (s.d.)] years, which confirms the previous findings. In addition, our results increase the number of the wavelengths at which the year-long period is detected. This fact, coupled with the findings of Sandrinelli et al. (2018) and Covino et al. (2020) strongly suggest that this period is real feature in the PG 1553+113 variability pattern.

The year-long periodicity in the PG 1553+113 emission could be explained by the jet precession in a SMBH binary system (Sobacchi et al. 2017; Caproni et al. 2017; Tavani et al. 2018). The precession will change the jet viewing angle and, hence, the corresponding Doppler factor. In the case of a continuous jet, the observed flux is related to the intrinsic one via F(obs) = δ(α + 2)F. If the viewing angle changes because of the jet precession then δ = δ(t) and, hence, F = F(t). We already assumed that the mean spline, derived in Sect. 3.4, represents the long-term variability component of geometric origin (jet precession). So, to account for the flux changes of the mean spline the Doppler factor has to be changed up to 1.3 times relative to its value corresponding to the minimum spline flux.

We also detected a possible secondary period of ∼210 days using the LCs corrected for the long-term variations. As long as the flux variations of the corrected LCs are achromatic, we should rely on geometric causes of this periodicity, for example, the emitting region follows a spiral path within the jet or the jet itself has a helical geometry. The jet wobbling in PG 1553+113, recently found by Lico et al. (2020) could contribute to some extent to the periodic STV. It is worth noting that Bhatta et al. (2016) found a secondary period of ∼400 days for the blazar OJ 287 in addition to the main one of ∼12 years.

The update and further analysis of the historical LCs are of substantial importance in order to confirm or reject the putative secondary period in PG 1553+113. The observed periodicity (if we assume it is real) coupled with the observed achromatic flux variation support the geometric origin of the STV and LTV of PG 1553+113.


1

The authors presented the most comprehensive search for periodicities at γ-rays using ten different methods along with the corresponding techniques for significance estimation of the periods found.

2

IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the National Science Foundation.

4

In what follows, the Steward Observatory and the Kanata telescope will be referred to as Steward and Kanata, respectively, for short.

5

This flare marks the brightest state of PG 1553+113 recorded in our historical LC.

6

The magnitude-to-flux transformation was done as follows. Firstly, the calibrated BVRI magnitudes of the blazar were corrected for the Galactic extinction using the values taken from the NASA/IPAC Extragalactic Database: AB = 0.188 mag, AV = 0.142 mag, AR = 0.113 mag, and AI = 0.078 mag. The extinction corrected magnitudes were then converted into fluxes using the zero points of Bessell et al. (1998).

7

This notation comes from X-ray passbands and means the hard X-rays lead the soft X-rays. The opposite is known as a “hard” lag.

9

The peak could be traced onto the V-band REDFIT periodogram as well, but at lower significance.

Acknowledgments

We thank the anonymous referee whose comments improved the paper. BM and LSM are supported by the Bulgarian National Science Fund of the Ministry of Education and Science under grant DN 18/13-2017. SZ acknowledges NCN grant No. 1028/29/B/ST9/01793. We thank TUBITAK National Observatory for partial support in using T60 and T100 telescopes with project number 1505 and 1486, respectively. AO was supported by the Scientific Research Project Coordination Unit of Ataturk University, Project ID 8418. AR acknowledges the Research Associate Fellowship with Order No. 03(1428)/18/EMR-II under the Council of Scientific and Industrial Research (CSIR). Data from the Steward Observatory spectropolarimetric monitoring project were used. This program is supported by Fermi Guest Investigator grants NNX08AW56G, NNX09AU10G, NNX12AO93G, and NNX15AU81G. Based on observations obtained with the Samuel Oschin 48-inch Telescope at the Palomar Observatory as part of the Zwicky Transient Facility project. ZTF is supported by the National Science Foundation under Grant No. AST-1440341 and collaboration including Caltech, IPAC, the Weizmann Institute for Science, the Oskar Klein Center at Stockholm University, the University of Maryland, the University of Washington, Deutsches Elektronen-Synchrotron and Humboldt University, Los Alamos National Laboratories, the TANGO Consortium of Taiwan, the University of Wisconsin at Milwaukee, and Lawrence Berkeley National Laboratories. Operations are conducted by COO, IPAC, and UW. The iPTF project is a scientific collaboration between Caltech; Los Alamos National Laboratory; the University of Wisconsin, Milwaukee; the Oskar Klein Centre in Sweden; the Weizmann Institute of Science in Israel; the TANGO Program of the University System of Taiwan; and the Kavli Institute for the Physics and Mathematics of the Universe in Japan. The CSS survey is funded by the National Aeronautics and Space Administration under Grant No. NNG05GF22G issued through the Science Mission Directorate Near-Earth Objects Observations Program. The CRTS survey is supported by the U.S. National Science Foundation under grants AST-0909182. Based on data acquired at Complejo Astronómico El Leoncito, operated under agreement between the Consejo Nacional de Investigaciones Científicas y Técnicas de la República Argentina and the National Universities of La Plata, Córdoba and San Juan (proposals JS-2019A-10, JS-2019A-16, HSH-2018B-03, and Staff time). This research has made use of the NASA/IPAC Extragalactic Database (NED), which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology.

References

  1. Abdo, A. A., Ackermann, M., Agudo, I., et al. 2010, ApJ, 716, 30 [NASA ADS] [CrossRef] [Google Scholar]
  2. Ackermann, M., Ajello, M., Albert, A., et al. 2015, ApJ, 813, L41 [NASA ADS] [CrossRef] [Google Scholar]
  3. Agarwal, A., & Gupta, A. C. 2015, MNRAS, 450, 541 [NASA ADS] [CrossRef] [Google Scholar]
  4. Agarwal, A., Mohan, P., Gupta, A. C., et al. 2017, MNRAS, 469, 813 [NASA ADS] [CrossRef] [Google Scholar]
  5. Agarwal, A., Cellone, S. A., Andruchow, I., et al. 2019, MNRAS, 488, 4093 [CrossRef] [Google Scholar]
  6. Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006, A&A, 448, L19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Ait Benkhali, F., Hofmann, W., Rieger, F. M., & Chakraborty, N. 2020, A&A, 634, A120 [CrossRef] [EDP Sciences] [Google Scholar]
  8. Aleksić, J., Ansoldi, S., Antonelli, L. A., et al. 2015, MNRAS, 450, 4399 [NASA ADS] [CrossRef] [Google Scholar]
  9. Alexander, T. 1997, in Is AGN Variability Correlated with Other AGN Properties? ZDCF Analysis of SmallSamples of Sparse Light Curves, eds. D. Maoz, A. Sternberg, & E. M. Leibowitz, Astrophys. Space Sci. Libr., 218, 163 [Google Scholar]
  10. Alexander, T. 2014, ZDCF: Z-Transformed Discrete Correlation Function [Google Scholar]
  11. Andruchow, I., Cellone, S. A., & Romero, G. E. 2007, Boletín de la Asociación Argentina de Astronomía, 50, 299 [Google Scholar]
  12. Andruchow, I., Combi, J. A., Muñoz-Arjonilla, A. J., et al. 2011, A&A, 531, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Andruchow, I., Cellone, S. A., & Romero, G. E. 2014, Rev. Mex. Astron. Astrofís. Conf. Ser., 44, 95 [Google Scholar]
  14. Bachev, R. 2018, Bulgarian Astron. J., 28, 22 [Google Scholar]
  15. Banerjee, B., Joshi, M., Majumdar, P., et al. 2019, MNRAS, 487, 845 [NASA ADS] [CrossRef] [Google Scholar]
  16. Bessell, M. S. 1990, PASP, 102, 1181 [NASA ADS] [CrossRef] [Google Scholar]
  17. Bessell, M. S., Castelli, F., & Plez, B. 1998, A&A, 333, 231 [NASA ADS] [Google Scholar]
  18. Bhatta, G., Zola, S., Stawarz, Ł., et al. 2016, ApJ, 832, 47 [NASA ADS] [CrossRef] [Google Scholar]
  19. Bonev, T. 2011, Gaia follow-up network for the solar system objects : Gaia FUN-SSO workshop proceedings, 89 [Google Scholar]
  20. Caproni, A., Abraham, Z., Motter, J. C., & Monteiro, H. 2017, ApJ, 851, L39 [NASA ADS] [CrossRef] [Google Scholar]
  21. Cellone, S. A., Romero, G. E., & Araudo, A. T. 2007, MNRAS, 374, 357 [NASA ADS] [CrossRef] [Google Scholar]
  22. Chiappetti, L., Maraschi, L., Tavecchio, F., et al. 1999, ApJ, 521, 552 [NASA ADS] [CrossRef] [Google Scholar]
  23. Clements, S. D., Jenks, A., & Torres, Y. 2003, AJ, 126, 37 [NASA ADS] [CrossRef] [Google Scholar]
  24. Covino, S., Sandrinelli, A., & Treves, A. 2019, MNRAS, 482, 1270 [NASA ADS] [CrossRef] [Google Scholar]
  25. Covino, S., Landoni, M., Sandrinelli, A., & Treves, A. 2020, ApJ, 895, 122 [CrossRef] [Google Scholar]
  26. Cutini, S., Ciprini, S., Stamerra, A., Thompson, D. J., & Perri, M. 2016, Active Galactic Nuclei 12: A Multi-Messenger Perspective (AGN12), 58 [Google Scholar]
  27. Drake, A. J., Djorgovski, S. G., Mahabal, A., et al. 2009, ApJ, 696, 870 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  28. Edelson, R. A., & Krolik, J. H. 1988, ApJ, 333, 646 [NASA ADS] [CrossRef] [Google Scholar]
  29. Falomo, R., & Treves, A. 1990, PASP, 102, 1120 [NASA ADS] [CrossRef] [Google Scholar]
  30. Falomo, R., Scarpa, R., & Bersanelli, M. 1994, ApJS, 93, 125 [NASA ADS] [CrossRef] [Google Scholar]
  31. Filippenko, A. V., Li, W. D., Treffers, R. R., & Modjaz, M. 2001, in The Lick Observatory Supernova Search with the Katzman Automatic Imaging Telescope, eds. B. Paczynski, W. P. Chen, & C. Lemme, ASP Conf. Ser., 246, 121 [Google Scholar]
  32. Gaur, H., Gupta, A. C., Strigachev, A., et al. 2012, MNRAS, 425, 3002 [NASA ADS] [CrossRef] [Google Scholar]
  33. Ghisellini, G., Villata, M., Raiteri, C. M., et al. 1997, A&A, 327, 61 [Google Scholar]
  34. Gopal-Krishna, Goyal, A., Joshi, S., et al. 2011, MNRAS, 416, 101 [NASA ADS] [Google Scholar]
  35. Green, R. F., Schmidt, M., & Liebert, J. 1986, ApJS, 61, 305 [NASA ADS] [CrossRef] [Google Scholar]
  36. Gupta, A. C., Cha, S.-M., Lee, S., et al. 2008, AJ, 136, 2359 [NASA ADS] [CrossRef] [Google Scholar]
  37. Gupta, A. C., Agarwal, A., Bhagwan, J., et al. 2016, MNRAS, 458, 1127 [CrossRef] [Google Scholar]
  38. Heeschen, D. S., Krichbaum, T., Schalinski, C. J., & Witzel, A. 1987, AJ, 94, 1493 [NASA ADS] [CrossRef] [Google Scholar]
  39. Howell, S. B., Mitchell, K. J., & Warnock, A. I. 1988, AJ, 95, 247 [NASA ADS] [CrossRef] [Google Scholar]
  40. Ikejiri, Y., Uemura, M., Sasada, M., et al. 2011, PASJ, 63, 639 [NASA ADS] [Google Scholar]
  41. Itoh, R., Nalewajko, K., Fukazawa, Y., et al. 2016, ApJ, 833, 77 [NASA ADS] [CrossRef] [Google Scholar]
  42. Jankowsky, F., & Wagner, S. 2019, ATel., 12631, 1 [Google Scholar]
  43. Jayasinghe, T., Stanek, K. Z., Kochanek, C. S., et al. 2019, MNRAS, 485, 961 [NASA ADS] [CrossRef] [Google Scholar]
  44. Jockers, K., Credner, T., Bonev, T., et al. 2000, Kinematika i Fizika Nebesnykh Tel Supplement, 3, 13 [NASA ADS] [Google Scholar]
  45. Johnson, S. D., Mulchaey, J. S., Chen, H.-W., et al. 2019, ApJ, 884, L31 [CrossRef] [Google Scholar]
  46. Kinman, T. D. 1975, in Variable Stars and Stellar Evolution, eds. V. E. Sherwood, & L. Plaut, IAU Symp., 67 [Google Scholar]
  47. Kirk, J. G., Rieger, F. M., & Mastichiadis, A. 1998, A&A, 333, 452 [NASA ADS] [Google Scholar]
  48. Kochanek, C. S., Shappee, B. J., Stanek, K. Z., et al. 2017, PASP, 129, 104502 [Google Scholar]
  49. Königl, A. 1981, ApJ, 243, 700 [Google Scholar]
  50. Law, N. M., Kulkarni, S. R., Dekany, R. G., et al. 2009, PASP, 121, 1395 [NASA ADS] [CrossRef] [Google Scholar]
  51. Li, W., Filippenko, A. V., Chornock, R., & Jha, S. 2003, PASP, 115, 844 [NASA ADS] [CrossRef] [Google Scholar]
  52. Lico, R., Liu, J., Giroletti, M., et al. 2020, A&A, 634, A87 [CrossRef] [EDP Sciences] [Google Scholar]
  53. Lomb, N. R. 1976, Ap&SS, 39, 447 [Google Scholar]
  54. Maoz, D., & Netzer, H. 1989, MNRAS, 236, 21 [NASA ADS] [Google Scholar]
  55. Marscher, A. P., & Gear, W. K. 1985, ApJ, 298, 114 [NASA ADS] [CrossRef] [Google Scholar]
  56. Masci, F. J., Laher, R. R., Rusholme, B., et al. 2019, PASP, 131, 018003 [NASA ADS] [CrossRef] [Google Scholar]
  57. Massaro, E., & Trèvese, D. 1996, A&A, 312, 810 [Google Scholar]
  58. Meng, N., Zhang, X., Wu, J., Ma, J., & Zhou, X. 2018, ApJS, 237, 30 [CrossRef] [Google Scholar]
  59. Miller, H. R., & Green, R. F. 1983, BAAS, 15, 957 [Google Scholar]
  60. Miller, H. R., Carini, M. T., Gaston, B. J., & Hutter, D. J. 1988, ESA SP, 2, 303 [Google Scholar]
  61. Monet, D. G., Levine, S. E., Canzian, B., et al. 2003, AJ, 125, 984 [NASA ADS] [CrossRef] [Google Scholar]
  62. Nilsson, K., Lindfors, E., Takalo, L. O., et al. 2018, A&A, 620, A185 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Osterman, M. A., Miller, H. R., Campbell, A. M., et al. 2006, AJ, 132, 873 [NASA ADS] [CrossRef] [Google Scholar]
  64. Padovani, P., & Giommi, P. 1995, ApJ, 444, 567 [NASA ADS] [CrossRef] [Google Scholar]
  65. Pandey, A., Gupta, A. C., Wiita, P. J., & Tiwari, S. N. 2019, ApJ, 871, 192 [CrossRef] [Google Scholar]
  66. Papadakis, I. E., Boumis, P., Samaritakis, V., & Papamastorakis, J. 2003, A&A, 397, 565 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  67. Papadakis, I. E., Villata, M., & Raiteri, C. M. 2007, A&A, 470, 857 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  68. Pasierb, M., Goyal, A., Ostrowski, M., et al. 2020, MNRAS, 492, 1295 [CrossRef] [Google Scholar]
  69. Peñil, P., Domínguez, A., Buson, S., et al. 2020, ApJ, 896, 134 [CrossRef] [Google Scholar]
  70. Peterson, B. M., Wanders, I., Horne, K., et al. 1998, PASP, 110, 660 [NASA ADS] [CrossRef] [Google Scholar]
  71. Piner, B. G., & Edwards, P. G. 2018, in Fourteenth Marcel Grossmann Meeting – MG14, eds. M. Bianchi, R. T. Jansen, & R. Ruffini, 3074 [Google Scholar]
  72. Press, W. H., & Rybicki, G. B. 1989, ApJ, 338, 277 [NASA ADS] [CrossRef] [Google Scholar]
  73. Prokhorov, D. A., & Moraghan, A. 2017, MNRAS, 471, 3036 [NASA ADS] [CrossRef] [Google Scholar]
  74. Raiteri, C. M., Stamerra, A., Villata, M., et al. 2015, MNRAS, 454, 353 [NASA ADS] [CrossRef] [Google Scholar]
  75. Rau, A., Kulkarni, S. R., Law, N. M., et al. 2009, PASP, 121, 1334 [NASA ADS] [CrossRef] [Google Scholar]
  76. Rector, T. A., & Perlman, E. S. 2003, AJ, 126, 47 [NASA ADS] [CrossRef] [Google Scholar]
  77. Romero, G. E. 1995, Ap&SS, 234, 49 [NASA ADS] [CrossRef] [Google Scholar]
  78. Romero, G. E., Cellone, S. A., & Combi, J. A. 1999, A&A, 135, 477 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  79. Sambruna, R. M., Maraschi, L., & Urry, C. M. 1996, ApJ, 463, 444 [NASA ADS] [CrossRef] [Google Scholar]
  80. Sandrinelli, A., Covino, S., Treves, A., et al. 2018, A&A, 615, A118 [CrossRef] [EDP Sciences] [Google Scholar]
  81. Scargle, J. D. 1982, ApJ, 263, 835 [Google Scholar]
  82. Schulz, M., & Mudelsee, M. 2002, Comput. Geosci., 28, 421 [NASA ADS] [CrossRef] [Google Scholar]
  83. Shappee, B. J., Prieto, J. L., Grupe, D., et al. 2014, ApJ, 788, 48 [NASA ADS] [CrossRef] [Google Scholar]
  84. Smith, P. S., Montiel, E., Rightley, S., et al. 2009, in 2009 Fermi Symposium, eConf Proceedings C0911022 [arXiv:0912.3621] [Google Scholar]
  85. Sobacchi, E., Sormani, M. C., & Stamerra, A. 2017, MNRAS, 465, 161 [NASA ADS] [CrossRef] [Google Scholar]
  86. Stalin, C. S., Gupta, A. C., Gopal-Krishna, Wiita, P. J., & Sagar, R. 2005, MNRAS, 356, 607 [NASA ADS] [CrossRef] [Google Scholar]
  87. Stetson, P. B. 1987, PASP, 99, 191 [Google Scholar]
  88. Stetson, P. B. 1992, in Astronomical Data Analysis Software and Systems I, eds. D. M. Worrall, C. Biemesderfer, & J. Barnes, ASP Conf. Ser., 25, 297 [Google Scholar]
  89. Takahashi, T., Tashiro, M., Madejski, G., et al. 1996, ApJ, 470, L89 [NASA ADS] [CrossRef] [Google Scholar]
  90. Takalo, L. O., Nilsson, K., Lindfors, E., et al. 2008, in Am. Inst. Phys. Conf. Ser., eds. F. A. Aharonian, W. Hofmann, F. Rieger, 1085, 705 [Google Scholar]
  91. Tavani, M., Cavaliere, A., Munar-Adrover, P., & Argan, A. 2018, ApJ, 854, 11 [NASA ADS] [CrossRef] [Google Scholar]
  92. Tavecchio, F. 2006, The Tenth Marcel Grossmann Meeting. On Recent Developmentsin Theoretical and Experimental General Relativity, Gravitation and Relativistic Field Theories, 512 [CrossRef] [Google Scholar]
  93. Tavecchio, F., Ghisellini, G., Ghirlanda, G., Foschini, L., & Maraschi, L. 2010, MNRAS, 401, 1570 [Google Scholar]
  94. Torres Zafra, J., Cellone, S. A., Andruchow, I., Buzzoni, A., & Portilla Barbosa, J. G. 2017, Rev. Mex. Astron. Astrofis. Conf. Ser., 49, 138 [Google Scholar]
  95. Treves, A., Falomo, R., & Uslenghi, M. 2007, A&A, 473, L17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  96. Valtaoja, E., Lähteenmäki, A., Teräsranta, H., & Lainela, M. 1999, ApJS, 120, 95 [NASA ADS] [CrossRef] [Google Scholar]
  97. VanderPlas, J. T. 2018, ApJS, 236, 16 [NASA ADS] [CrossRef] [Google Scholar]
  98. Villata, M., Raiteri, C. M., Kurtanidze, O. M., et al. 2002, A&A, 390, 407 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  99. Villata, M., Raiteri, C. M., Kurtanidze, O. M., et al. 2004, A&A, 421, 103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  100. Wagner, S. J., & Witzel, A. 1995, ARA&A, 33, 163 [NASA ADS] [CrossRef] [Google Scholar]
  101. White, R. J., & Peterson, B. M. 1994, PASP, 106, 879 [NASA ADS] [CrossRef] [Google Scholar]
  102. Wierzcholska, A., Ostrowski, M., Stawarz, Ł., Wagner, S., & Hauser, M. 2015, A&A, 573, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  103. Xie, G. Z., Zhou, S. B., Li, K. H., et al. 2004, MNRAS, 348, 831 [NASA ADS] [CrossRef] [Google Scholar]
  104. Zibecchi, L., Andruchow, I., Cellone, S. A., et al. 2017, MNRAS, 467, 340 [NASA ADS] [Google Scholar]
  105. Zibecchi, L., Andruchow, I., Cellone, S. A., & Carpintero, D. D. 2020, MNRAS, 498, 3013 [CrossRef] [Google Scholar]

Appendix A: Results from INM

We present the results from the INV tests in Table A.1. The individual intra-night LCs, derived during the 28 nights of INM, are shown in Fig. A.1.

thumbnail Fig. A.1.

Intra-night LCs for PG 1553+113: blue denotes B-band, green – V-band, red – R-band, and black – I-band. In each plot the x-axis is the duration of the INM in hours. Date of observation and the telescope used are indicated in each plot.

thumbnail Fig. A.1.

continued.

thumbnail Fig. A.1.

continued.

Table A.1.

Results of INV observations of PG 1553+113.

Appendix B: Historical light curves

In order to study in detail the variability of PG 1553+113 on long-term timescales, we constructed its historical LC. We used only major public data sets, so, this LC will not include all available photometric data on the blazar.

The first data set we considered was the one obtained at the Steward Observatory10 of the University of Arizona. The observations were performed from 2012 to 2018 through VR filters and were described in Smith et al. (2009). Comparing Steward photometry to ours, we found 7 V and 9 R data points that have been obtained either in one and the same night or in two consecutive nights – in all cases, the magnitude differences were less than the corresponding uncertainties. Therefore, there is no need the Steward photometry to be adjusted to ours and we merged the two data sets into a basic set to which the other archival data will be matched. After each successful matching the basic data set is updated accordingly.

The next major VR data set comes from Itoh et al. (2016). The observations were performed from 2008 to 2014 using the Kanata telescope. This data set was adjusted and then added to the basic one. We have only one common observing night, so we were forced to use an interpolation onto the overlapping parts of the LCs to do the inter-calibration; this procedure was favoured by the fact that the source did not show significant variability during the period of interest.

The V-band LC was successively supplemented with the following data sets:

  1. Catalina Surveys (Drake et al. 2009, The Catalina Surveys Data Release 2). The observations have been obtained without a filter from 2005 to 2014 and the photometry has been transformed by the Catalina pipeline to V-band11. The inter-calibration with the updated basic data set was performed on the base of 33 observing nights in common.

  2. All-Sky Automated Survey for Supernovae12 (ASAS-SN, Shappee et al. 2014; Kochanek et al. 2017; Jayasinghe et al. 2019). We matched ASAS-SN photometry (taken from 2013 to 2018) to the updated basic data set using 21 nights in common.

  3. The Zwicky Transient Facility13 (ZTF, Masci et al. 2019). The data time span is from 2018 to 2019. Before adjustment to the updated basic data set, the ZTF gr magnitudes were transformed to the Johnson-Cousins VR ones. Then the matching was done using 12 nights in common.

The following data sets were added to the R-band LC:

  1. Katzman Automatic Imaging Telescope14 (KAIT, Filippenko et al. 2001). The KAIT unfiltered magnitudes roughly correspond to the R-band (Li et al. 2003, we used the results since March 2019 from a new pipeline). The matching was done on the base of 34 nights in common. We also added to the historical LC the KAIT results obtained with the old pipeline before March 2019 that are missing among the results obtained with the new pipeline. The old photometry has been calibrated with respect to the “R2” magnitudes of USNO-B1 catalogue (Monet et al. 2003). In this case the matching was done on the base of 15 nights in common. The KAIT photometry needs no adjustment before its joining to the basic data set. The time span of KAIT data is from 2010 to 2019.

  2. Tuorla Blazar Monitoring Program15 (Takalo et al. 2008). The observations were performed from 2005 to 2012 in the R-band (Nilsson et al. 2018) and were adjusted to the updated basic data set on the base of 35 nights in common.

  3. Intermediate Palomar Transient Factory16 (iPTF Law et al. 2009; Rau et al. 2009). The observing time span is from 2009 to 2010 and the matching was done using four nights in common.

  4. ZTF (see above). The matching was done using 18 nights in common.

After the different data sets were put together accounting for the offsets existing among them, nightly binning (and cleaning of the deviant data points) was performed as described in Sect. 3.3. The so assembled historical VR LCs for PG 1553+113 are presented in Fig. B.1. The variability amplitude of the historical VR-band LCs is A ≃ 1.4 mag.

thumbnail Fig. B.1.

Long-term VR-band LCs for PG 1553+113. Different colours and symbols denote data from different observatories and telescopes as indicated in the plot. In particular, “others” denote the data collected by the nine telescopes mentioned in Sect. 2 The data plotted are not nightly binned.

All Tables

Table 1.

Details about the telescopes and instruments used.

Table 2.

Log of photometric observations for the blazar PG 1553+113.

Table 3.

Long-term variability characteristics.

Table 4.

Flare best fit parameters – observer-frame.

Table 5.

Summary of the periods found in the historical VR-band LCs of PG 1553+113.

Table A.1.

Results of INV observations of PG 1553+113.

All Figures

thumbnail Fig. 1.

Long-term LC for PG 1553+113 covering the entire monitoring duration. Different colours and symbols denote data in different passbands as indicated in the plot. In each passband the magenta solid lines connect the data points obtained after nightly binning of the respective LC.

In the text
thumbnail Fig. 2.

Colour-magnitude diagram for our data only.

In the text
thumbnail Fig. 3.

Colour-magnitude diagram for the historical data set (upper panel, 373 data points plotted) and for our, Steward, and Kanata data sets (lower panel, 274 data points plotted). The linear least-squares fits are overplotted. For both panels the median uncertainties of the magnitudes and colour indices are 0.018 mag and 0.036 mag, respectively.

In the text
thumbnail Fig. 4.

Temporal evolution of the V − R colour index (black dots). The scaled R-band LC is overplotted for comparison (magenta line, the flux increases along +y axis). The horizontal dashed line marks the weighted mean V − R value of (0.339 ± 0.002) mag with a weighted standard deviation of 0.030 mag about the weighted mean. The insert shows the 2019 flare in details – the maximum of the spectral hardness (the minimum of V − R) precedes the flare by about 10 days (see Sect. 3.5 for details).

In the text
thumbnail Fig. 5.

Historical VR-band LCs (dots). The cubic splines, fitted through the 10-day binned LCs, are overplotted (V-band – green line; R-band – red line) along with the shifted mean spline representing the base levels for the VR-bands (black lines).

In the text
thumbnail Fig. 6.

Historical VR-band LCs corrected for the long-term variations.

In the text
thumbnail Fig. 7.

Flux ratio FV/FR plotted against (FV + FR)/2. Shown are our, Steward, and Kanata data sets corrected for the long-term variations. The best fit linear model to the data is overplotted. The median uncertainties of the fluxes and flux ratios are 0.082 mJy and 0.087 mJy, respectively.

In the text
thumbnail Fig. 8.

Colour-magnitude diagram for the 2019 flare. The corresponding spectral index: α = [(V − R)−(AV − AR)−0.186]/0.07, is also indicated. Here AV and AR stand for the Galactic extinction in the VR-bands, respectively (see Sect. 3.4). The clockwise spectral hysteresis can be clearly seen – the arrows are used to guide the eye about the hysteresis loop direction. The colours and numbers are sequential in time: violet (#1)→ red (#44).

In the text
thumbnail Fig. 9.

Discrete correlation function. The fitted Gaussian is not shown for the sake of clarity.

In the text
thumbnail Fig. 10.

Cross-correlation peak distribution for DCF.

In the text
thumbnail Fig. 11.

Same as in Fig. 9, but for ZDCF.

In the text
thumbnail Fig. 12.

Approximation of the 2019 flare with an exponential law plus a local linear base level. The possible sub-flares could be discerned around MJD = 8550 and MJD = 8610.

In the text
thumbnail Fig. 13.

Sine function fits to the VR splines: V-band – lower green curve with black dots overplotted, R-band – upper red curve with black dots overplotted, sine fit – black curves of constant amplitude.

In the text
thumbnail Fig. 14.

Lomb-Scargle periodogram for the historical V-band LC. The dotted line marks the level corresponding to a maximum peak false alarm probability of 1% under the assumption of white noise. The insert shows the LC folded with the most significant period; the red line denotes the fitted sine wave.

In the text
thumbnail Fig. 15.

Same as in Fig. 14, but for the R-band.

In the text
thumbnail Fig. 16.

Output from the REDFIT programme. The black line shows the bias-corrected power spectrum, the blue dashed line marks the theoretical red noise spectrum and the red line marks the 99% local significance level derived by means of Monte Carlo simulations.

In the text
thumbnail Fig. 17.

Same as in Fig. 14, but for the historical V-band LC corrected for the long-term variations.

In the text
thumbnail Fig. 18.

Same as in Fig. 17, but for the R-band.

In the text
thumbnail Fig. 19.

Same as in Fig. 16, but for the historical VR-band LCs corrected for the long-term variations.

In the text
thumbnail Fig. 20.

Optical SEDs of PG 1553+113 for selected brightness levels (see text).

In the text
thumbnail Fig. A.1.

Intra-night LCs for PG 1553+113: blue denotes B-band, green – V-band, red – R-band, and black – I-band. In each plot the x-axis is the duration of the INM in hours. Date of observation and the telescope used are indicated in each plot.

In the text
thumbnail Fig. B.1.

Long-term VR-band LCs for PG 1553+113. Different colours and symbols denote data from different observatories and telescopes as indicated in the plot. In particular, “others” denote the data collected by the nine telescopes mentioned in Sect. 2 The data plotted are not nightly binned.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.