© The Authors 2024

1. Introduction

Blazars, a subtype of radio-loud active galactic nuclei (AGN), develop relativistic jets that are closely aligned with the line of sight, which leads to a relativistic boosting of their emission. They can be divided into BL Lacertae (BL Lac) objects and flat-spectrum radio quasars (FSRQs) depending on the properties of their optical spectrum (Urry & Padovani 1995). BL Lacs typically present an (almost) featureless optical spectrum, while FSRQs show broad emission lines with an equivalent width |EW|> 5 Å in their rest frame (Stickel et al. 1991). Their broadband emission is mainly non thermal, is associated to the relativistic jet, and can extend from radio frequencies up to γ-ray energies. Their spectral energy distribution (SED) shows a typical double-bump structure (see e.g., Abdo et al. 2010). The low-energy bump is produced by synchrotron radiation of the electrons moving under the influence of the magnetic field of the jet (see e.g., Konigl 1981). The high-energy bump is commonly modeled through inverse Compton (IC) scattering of low-energy photons with the same population of relativistic electrons, that is, within the leptonic interpretation. This can happen with the same low-energy photons of the synchrotron radiation through synchrotron self-Compton scattering (SSC; see Maraschi et al. 1992), or with low-energy photons from outside the jet – if there is an external injection of photons – through external Compton scattering (EC; see Dermer & Schlickeiser 1993). While this is the commonly adopted explanation, hadronic models have also been used in recent years (see e.g., Cerruti et al. 2015).

Variability is a key signature of the blazar emission. This variability is present in their broadband emission, manifesting on all timescales (Fan 2005; Singh & Meintjes 2020). Variations in their emission on timescales of several months or years are a typical feature of the emission of blazars, and this is known as long-term variability (LTV; see e.g., Rajput et al. 2020; Bhatta 2021). Faster variations are also often observed, that is, on timescales of several days or weeks, as short-term variability (STV; see Gupta et al. 2008; Rani et al. 2010). In addition, in some cases blazars have shown remarkably fast variability, on timescales of shorter than one day, down to even a few minutes, this is intraday variability (IDV; e.g., Wagner & Witzel 1995; Villata et al. 2002; Raiteri et al. 2021a,b). Different timescales have been associated with different processes and emission mechanisms. However, variability patterns on different timescales are typically entangled and very complex, resulting in a challenging physical interpretation. Nevertheless, variability has become a major tool for understanding the physics taking place in the relativistic jets of blazars.

3C 371 is a BL Lac object located at a distance of z = 0.0510 ± 0.0003, as estimated from optical spectroscopic observations by de Grijp et al. (1992). Miller (1975) also classified this source as an intermediate case between the BL Lac type and a radio galaxy. 3C 371 is hosted by a bright galaxy with R = 14.22 ± 0.03 (Pursimo et al. 2002; Nilsson et al. 2003), which makes an important emission contribution in the optical and near-infrared bands (Raiteri et al. 2014). It has been reported to show intense variability in its broadband emission on a wide range of timescales: several years in radio (Nieppola et al. 2009) and ultraviolet (UV) wavelengths (Paltani & Courvoisier 1994); changes of more than 1 mag on timescales of years and smaller variations in scales of several days in the optical (see Oke 1967); or long-term variations in X-ray observations performed by EXOSAT (Giommi et al. 1990). It has also been cataloged as variable in the γ-ray regime, as shown by the results presented in the Fermi-LAT Fourth Source Catalog Data Release 3 (4FGL-DR3) catalog, where a variability index of ∼133 is reported (Abdollahi et al. 2020, 2022). Moreover, faster variability has also been observed for 3C 371. For instance, Kraus et al. (2003) report this source as a fast variable radio source, with variability on timescales of ≤2 days which is, similar to the 2.4 day timescales observed by Heidt & Wagner (1996) in the optical R-band. However, interstellar scintillation cannot be ruled out as the origin of such fast variability in the radio band. It has also been reported to show IDV in the optical band by Xilouris et al. (2006), making this object an extremely interesting target for performing variability studies down to the shortest variability timescales.

These interesting variability patterns led to an extensive monitoring of 3C 371 by the Whole Earth Blazar Telescope¹ (WEBT) Collaboration in optical and radio frequencies. A proposal to the Transiting Exoplanet Survey Satellite (TESS) led to observations with an exceptional cadence of 2 min between 2019 and 2020. Taking advantage of these extensive datasets, we have conducted a variability study of 3C 371, evaluating the different characteristic variability timescales in its optical emission between 2018 and 2020, when the source was extensively monitored by the WEBT Collaboration.

The present paper is structured as follows: All the optical datasets are detailed in Sect. 2. We present an evalutaion of the IDV in Sect. 3 that makes use of the TESS data. A variability analysis taking advantage of all the WEBT data collected is presented in Sect. 4. In Sect. 5 we show how we evaluated the optical spectral behavior during the observed period. A study of the variability of the optical polarized emission is presented in Sect. 6 and a discussion of the possible origin of the observed variability in relation to the measured flux distributions on different timescales is presented in Sect. 7. Finally, the main conclusions of this work are summarized in Sect. 8.

2. Observations and data reduction

2.1. TESS photometric data

3C 371 has been a target of the TESS mission with an observing cadence of 2 min for almost 1 year, with data available during 12 TESS observing sectors during cycle 2, detailed in Table 1, leading to a total of 214236 data points. These data, downloaded from the Mikulski Archive for Space Telescopes² (MAST), contain information of the simple aperture photometry electron flux (SAP_FLUX, in electrons per second), that will be used hereafter for the data reduction and analysis presented in this study (Feinstein et al. 2019). They also contain presearch conditioned simple aperture photometry fluxes (PDCSAP_FLUX). However, as discussed by Raiteri et al. (2021a,b), they are not suitable for blazar variability studies as they remove long-term trends in the data that may be genuine variability signatures of the source. SAP_FLUX have a poorer background subtraction. However, owing to the brightness of the source and the good agreement with the WEBT data, we found SAP_FLUX values to be suitable for this study. In order to correct for systematic scaling offsets between different sectors, we scaled each sector by a certain offset, calculated from a direct comparison with the WEBT R-band light curve (see Sect. 2.2). Moreover, we removed outliers and points showing an extremely high scattering (∼0.01% of the points), resulting in a total of 214216 observations included in the light curve. Finally, we also normalized the scaled light curve by the mean SAP_FLUX value. The light curve of the scaled, normalized aperture photometry flux is shown in Fig. 1.

Fig. 1. TESS SAP_FLUX light curve of 3C 371. The offsets between the different sectors were estimated using the observed R-band WEBT light curve, and were then converted into arbitrary flux units and normalized by the mean value. The inset panel corresponds to a zoom onto a 1 h of the light curve to illustrate the cadence and uncertainties of the observations.

2.2. WEBT optical observations

This blazar is also a regularly monitored target by the WEBT Collaboration with several telescopes distributed around the globe. It has been observed in the Johnson–Cousins BVRI bands, with a specially remarkable coverage in the R band. In particular, we focus here in the period between January 2018 and December 2020. The observed optical light curves constructed with the data taken by the WEBT Collaboration are displayed in Fig. 3, and all the telescopes that provided data are listed in Table 2.

Data were processed with standard reduction procedures. The magnitude of 3C 371 was obtained with aperture photometry and calibrated with respect to reference stars in the same field of the source. A photometric sequence including five stars was published by Xilouris et al. (2006), but only in the B and I bands. We then used data acquired at the Teide Observatory in good nights to derive new calibrations in all BVRI bands for the same stars (see Fig. 2). The standard magnitudes obtained by our procedure are reported in Table 3. Values in B and I bands are in good agreement with those published by Xilouris et al. (2006). Magnitude offsets between the various datasets, which were mainly due to different choices of the aperture radius for the source photometry, leading to the inclusion of different contributions from the host galaxy light, were corrected by aligning all datasets to the trend traced by those that adopted a 7.5″ aperture radius. Data points with very large errors and clear outliers were removed. Sequences of largely scattered data taken close in time by the same telescope were binned over a few minutes intervals.

Fig. 2. Finding chart obtained from the PanSTARRS-1 (PS1) Image Cutout Service (http://ps1images.stsci.edu/cgi-bin/ps1cutouts) in r band. The field of view is 8 × 8 arcmin. Stars 1−5 represent the photometric sequence of Xilouris et al. (2006) recalibrated by us (see Table 3).

In addition to the data reduction specified above, we have also taken into account the effect of the Galactic extinction. We corrected the observed magnitudes using the Galactic extinction values A_B = 0.123, A_V = 0.093, A_R = 0.074 and A_I = 0.051, extracted from the NASA/IPAC Extragalactic Database³ (NED) and reported by Schlafly & Finkbeiner (2011). Moreover, we also converted the magnitude measurements into flux units in order to construct the optical SEDs presented in Sect. 5. For this, we made use of the zero-mag fluxes provided by Bessell et al. (1998).

Finally, we have also taken into account the emission introduced by the host galaxy. The R-band magnitude of the host galaxy and the radius containing half of the light, r_e, were estimated to be R_host = 14.22 ± 0.03 and r_e = 8.6″ ± 0.1″ by Nilsson et al. (2003). We used a De Vaucouleurs profile (De Vaucouleurs 1948) to model the brightness profile of the host galaxy in order to estimate its flux contribution within the 7.5″ aperture. This was found to correspond to ∼46% of the total light of the host galaxy. We adopted the SWIRE template⁴ (Polletta et al. 2007) for a 13 Gyr old Elliptical Galaxy to infer the host contributions within a 7.5″ aperture radius in the other bands, using the R-band brightness as normalization. These contributions are reported in Table 4 for the different filters, before and after correcting for Galactic extinction.

The WEBT R-band data were also used to convert the TESS counts flux to R-band magnitudes. This calibrated R-band light curve is shown in the third panel of Fig. 3 with the WEBT data from the same temporal period. The calibration was performed as m_TESS = −2.5log(SAP_FLUX) + m₀, where m₀ is an offset calculated by comparing each TESS observing sector with simultaneous R-band WEBT data, accordingly to previous studies by Raiteri et al. (2021a,b). As can be observed in Fig. 3, the calibration leads to an excellent agreement between the WEBT and TESS light curves. The differences in the variability amplitude observed between the different optical filters and the TESS data are related to the spectral response function and the spectral range covered by TESS⁵. This spans from 600 nm to 1000 nm, covering all the spectral range of the I band, most of the R band, but it has only a partial overlap with the V band and almost no overlap in the B band. We note however that the calibration was performed with respect to the optical R band owing to the expected wider variations at this wavelength in comparison to the I band and the fact that is the best sampled among the WEBT light curves (see Raiteri et al. 2021a,b).

Fig. 3. WEBT multiband optical light curves (observed magnitudes) of 3C 371 during the 2018–2020 observing season. Different colors and symbols are used to distinguish the contributing datasets, as specified in Table 2. Vertical lines represent the start of the different TESS sectors during which the source was observed. The R-band calibrated TESS data are represented with gray markers in the third panel.

2.3. WEBT optical polarimetric observations

In addition to the total optical observations in the BVRI bands, the WEBT Collaboration has also performed a monitoring of the polarized emission of 3C 371 in the optical band making use of the telescopes listed in Table 2 with available polarimetric instruments, that is, the Crimean, San Pedro Martir and St. Petersburg Observatories. This monitoring obtained long-term measurements of the polarization degree and the electric vector position angle (EVPA) during the 3 year monitoring period. The EVPA has been corrected from the ±n ⋅ 180° ambiguity following the procedure typically used in the literature (adding/subtracting n ⋅ 180° to minimize the difference between consecutive measurements when Δθ_i = |θ_i + 1 − θ_i|> 90°, see for instance Raiteri et al. 2023a; Otero-Santos et al. 2023a). The observed polarization degree and EVPA trends in time are represented in Fig. 4 with the R-band light curve for comparison.

Fig. 4. Optical polarization of 3C 371 during the 2018–2020 period compared to the jet optical flux. From top to bottom: R-band optical light curve (dereddened and host-galaxy-subtracted flux densities), observed polarization degree, intrinsic polarization degree, corrected EVPA behavior. Different colors and symbols are used to distinguish the contributing datasets, as specified in Table 2.

However, the contribution of the host galaxy to the total emission from blazars has a depolarizing effect of the polarization degree intrinsic to the jet (see e.g., Sosa et al. 2017). Making use of the host galaxy contribution estimation performed in Sect. 2.2, we calculated the intrinsic polarization degree as

$$\begin{array}{c}\hfill P_{\mathrm{jet}}[\%]=\frac{P_{\mathrm{obs}}[\%]}{1-\frac{F_{\mathrm{host}}}{F_{\mathrm{total}}}},\end{array}$$(1)

where P_jet and P_obs correspond to the intrinsic and observed polarization fractions, F_total is the total brightness of the source and F_host is the contribution of the host galaxy. This intrinsic polarization degree is also represented in Fig. 4. We note that no bias correction was applied here to the polarization to account for the difference between the Rice and Gaussian distributions, as discussed by Blinov et al. (2021). Further discussion and analysis of the behavior of the polarization degree and EVPA of 3C 371 during this period is presented in Sect. 6.

3. Intra-day variability

Several blazars in the past have shown signatures of IDV (see e.g., Wagner & Witzel 1995; Villata et al. 2002; Raiteri et al. 2008, 2021a,b; Gupta et al. 2008; Gaur et al. 2015). In particular for 3C 371, Xilouris et al. (2006) have investigated the presence of IDV features in this source. Here, we take advantage of the extended temporal coverage provided by the TESS satellite for 3C 371 with a 2 min cadence to evaluate the possible IDV present in this BL Lac object. For this analysis we use the R-band calibrated light curve after Galactic extinction correction and host galaxy subtraction.

3.1. Variability amplitude and timescales

Several episodes of noticeable IDV in the R-band emission of 3C 371 are observed over the period monitored by TESS in the R band emission. Some remarkable periods of IDV can be observed during different flares occurring on MJD 58741, MJD 58803 or MJD 58916 as displayed in Fig. 1. The total variability amplitude over the roughly 1 year period is ΔR_max = 1.27 mag, varying between magnitudes R_max = 15.51 and R_min = 14.24. We also evaluated the maximum magnitude intraday variations within different time intervals ranging from 24 h down to 10 min. In Fig. 5 we represent the distributions of the maximum magnitude variations for the different time intervals considered here. The peak of the distribution shifts towards higher values of $\mathrm{\Delta }{R}_{\max }^{\mathrm{IDV}}$ for longer time intervals. The maximum noticeable variation that we observe in our entire dataset within a time interval of 10 min is $\mathrm{\Delta }{R}_{\max }^{\mathrm{IDV}}=0.204$ mag, while changes of $\mathrm{\Delta }{R}_{\max }^{\mathrm{IDV}}=0.431$ mag are observed within a 24-h interval. This is also visible in the tails of the distributions shown in Fig. 5, with for example the 6-h histogram reaching values of $\mathrm{\Delta }{R}_{\max }^{\mathrm{IDV}}\sim 0.35$ mag, while the 10 min distribution extends only up to $\mathrm{\Delta }{R}_{\max }^{\mathrm{IDV}}\sim 0.25$ mag. We note that we have taken into account the gaps existent in the light curves due to the end of each observing sector, as well as the gaps introduced in the middle of each sector by the data transfer, by not considering the magnitude difference between the last point of each sector and the first of the next one, and the last point of the first half and first of the second half of each sector.

Fig. 5. Distributions of the maximum IDV magnitude variations observed from TESS data in different time intervals, which are distinguished with different colors and line styles. The vertical black dashed line represents the mean uncertainty of the data.

We can quantify the amplitude of these intraday flux variations relative to the average emission by using the variability amplitude parameter, A_mp(%), following the definition from Romero et al. (1999),

$$\begin{array}{c}\hfill A_{\mathrm{mp}}(\%)=\frac{100}{⟨F⟩}\sqrt{{(F_{\max }-F_{\min })}^2-2{\sigma }^2},\end{array}$$(2)

where ⟨F⟩ is the mean value of the flux in the considered interval, F_max and F_min are the maximum and minimum flux values measured, and σ is the flux uncertainty. The associated uncertainty of A_mp can be estimated as

$$\begin{array}{cc}\hfill \mathrm{\Delta }{A}_{\mathrm{mp}}(\%)=& 100\times \left(\frac{F_{\max }-F_{\min }}{⟨F⟩A_{\mathrm{mp}}}\right)\hfill \\ \hfill & \times \sqrt{{\left(\frac{{\sigma }_{\max }}{⟨F⟩}\right)}^2+{\left(\frac{{\sigma }_{\min }}{⟨F⟩}\right)}^2+{\left(\frac{{\sigma }_{⟨F⟩}}{F_{\max }-F_{\min }}\right)}^2{A}_{\mathrm{mp}}^4}.\hfill \end{array}$$(3)

In this last equation, σ_⟨F⟩ corresponds to the standard deviation of the mean flux value. Considering the complete 2-min cadence TESS light curve, the estimated value for the variability amplitude is A_mp = (123.4 ± 1.1)%. We report in Table 5 the derived values of A_mp for the entire dataset, as well as each sector individually. We observe than on shorter timescales of ∼1 month, i.e. the duration of each sector, the variability is smaller, with variability amplitudes between ∼20% and ∼50% for the least and most variable sectors, respectively.

Moreover, we also estimated A_mp for the flux density variations occurring in time intervals of 10 min, 30 min, 1 h, 3 h, 6 h, 12 h and 1 day. The results are displayed in Table 6 for each time interval considered. We observe that for the time intervals evaluated here, the amplitude of the variations is ≲13.08%, with maximum changes between ∼19% and ∼36%, depending on the width of the interval. Xilouris et al. (2006) report a very low amplitude of the variations (∼1−1.5%) on timescales as short as a few hours. While this value is lower with respect to our estimations, there are some differences in the analyses that are probably leading to this discrepancy. Xilouris et al. (2006) perform a 1.5-h binning that may be smoothing real variations. Their observing cadence is also much lower than that from TESS data. Finally, the flux contribution of the host galaxy, not subtracted in the study by Xilouris et al. (2006) can also be leading to a lower estimation of the variability of the jet.

In order to quantify the fastest variability that we are able to significantly resolve, we follow a similar approach to that from Weaver et al. (2020), where an ANOVA (ANalysis Of VAariance) test was implemented to evaluate the IDV of BL Lacertae. This test has also been used in other variability studies in blazars (see e.g., de Diego 2010; Fraija et al. 2017; Weaver et al. 2019). We divide the dataset in intervals of different width. Then, breaking each interval in subgroups, we evaluate the null hypothesis of every subgroup having the same mean value, which can be translated into a nonvariable nature. When the contrary occurs with a confidence of ∼3σ (p-value ≤ 0.001), we can assume that the source is variable on timescales shorter than the interval considered (see for more details Weaver et al. 2020). We report in Table 6 the number of intervals where we can reject the null hypothesis of non-variability for the different time interval widths considered.

We observe that for most of the intervals corresponding to Δt ≤ 1 h, we barely see any sign of significant variability, with only 24 one-hour intervals fulfilling the aforementioned 3σ variability criterion. On the other hand, for Δt ≥ 6 h we can reject the nonvariable null hypothesis, as at least ∼36% of the intervals show significant variability according to the ANOVA test. Considering the 6-h intervals, we can determine the shortest variability timescale that we are able to resolve in our dataset following the procedure from Weaver et al. (2020, for more details see also Burbidge et al. 1974). This procedure determines, for all pairs of points of the subset of variable 6-h time intervals with S​₂ − S​₁ > 3(σ₁ + σ₂)/2 and S​₂ > S​₁, their characteristic variability timescale as

$$\begin{array}{c}\hfill \tau =\mathrm{\Delta }t/ln(S_2/S_1),\end{array}$$(4)

where S​₁ and S​₂ are the flux values of each point of the pair, σ₁ and σ₂ are their corresponding uncertainties, and Δt = |t₂ − t₁| is the time separation between them. We note that this method is applied on the flux measurements rather than on the magnitudes. Therefore, as in Weaver et al. (2020), we use for this the TESS dereddenned, host galaxy subtracted flux density light curve. In Fig. 6 we represent a histogram of the characteristic variability timescales determined with Eq. (4). While most of the variability shows timescales longer than 1 day, with a mean value of ⟨τ⟩ = 56.5 h and a standard deviation of 24.5 h, we can still observe significant IDV signatures. This is expected, given the long-term coverage of the TESS data, spanning almost 1 year and therefore allowing us to observe variations on much longer timescales than IDV. Here we focus on the significant IDV signatures identified with the ANOVA analysis and timescale quantification commented above. The fastest variability we were able to measure is quite remarkable, with a characteristic timescale of τ_min = 0.47 h.

Fig. 6. Histogram of the minimum variability timescales measured from the variable intervals according to the ANOVA test performed here. The inner panel shows a zoomed-in histogram for timescales of shorter than 6 h.

Apart from this fast variability, and as can be observed in the TESS light curve represented in Fig. 1, several periods of variability are characterized with different variability timescales longer than the minimum variability timescale estimated above. For evaluating the presence of other timescales involved in the IDV observed, we used the Fourier-based technique of the wavelet analysis, which has long been used as a powerful tool for evaluating in the past variability timescales in blazars, owing to its ability to decompose the signal in the two-dimensional time-frequency space, allowing the user to find the different timescales present in a dataset, and contextualize them in time (see Gupta et al. 2009; Zhou et al. 2018; Raiteri et al. 2023a). When these timescales are persistent in time, they reveal the presence of a periodic variability of the data. Hence, this tool has also been widely used for performing periodicity searches in the emission of blazars (e.g., Otero-Santos et al. 2020, 2023b; Peñil et al. 2020; Jorstad et al. 2022; Roy et al. 2022; Li et al. 2023; Tripathi et al. 2024). Moreover, it has been also employed for identifying nonperiodic characteristic variability timescales in blazar light curves (see Raiteri et al. 2023a). Here we use the PYTHON implementation of the wavelet function code developed by Torrence & Compo (1998)⁶.

It is well known that long-term trends can mask the impact of variability on shorter timescales (see Raiteri et al. 2021a,b). Owing to the roughly 1 year time span of the TESS data, we have modelled and subtracted variability signatures on timescales of longer than 1 day, which are further studied using the 3 year data from WEBT (see Sect. 4). Following the procedure employed by Raiteri et al. (2021a,b), we model the long-term variations through a cubic spline interpolation of a binned light curve. This cubic spline interpolation has been widely used in blazar studies to represent the observed LTV in these sources (e.g., Ghisellini et al. 1997; Raiteri et al. 2017). In our case, we perform the interpolation on the daily binned TESS light curve. In Fig. 7 we show the derived wavelet spectrum with associated frequency-space power spectral density (PSD).

Fig. 7. IDV wavelet analysis performed on the dereddened, host-galaxy-subtracted R-band flux light curve of 3C 371. Top: R-band light curve. Middle: time-period power spectrum (left) and time-folded power spectrum (right). White contours in the left panel and the dashed line in the right panel represent variability signatures with a significance > 3σ (> 99.73% confidence level). Bottom: averaged variance between 4 min and 1 day. The horizontal dashed line represents the 3σ (99.73%) confidence level.

We observe both from the time-period spectrum (middle left panel) and the period-folded PSD (middle right panel) represented in Fig. 7 that a large amount of the IDV variability of 3C 371 appears on timescales between 0.5 days and 3−5 days. Nevertheless, several features of faster variations are visible during different periods, most of them coincident with bright, fast developing flares (e.g. on MJD 58916, approximately). A wide variety of intraday characteristic timescales are easily identified at a 3σ confidence level. As an example, we represent the wavelet analysis performed on sectors 18 and 22 individually, where interesting flaring events on very fast timescales can be identified. The results of these analyses are shown in Fig. 8. Again, we see several signatures on timescales ranging from 0.5 days to ∼3 day. In addition, we are able to resolve even faster variations during bright flares occurring in these particular sectors. For instance, during the flare on MJD 58916, the power spectrum reveals significant flux variability down to timescales < 0.25 days (see right panels in Fig. 8). Another example is the flare on MJD 58802 (left panels in Fig. 8), where we can observe one of the fastest signatures of 3C 371 with this analysis, showing a characteristic timescale ≤0.05 days, compatible with the fastest signatures found with the ANOVA approach.

Fig. 8. IDV wavelet analysis performed on the dereddened, host-galaxy-subtracted R-band flux light curves of 3C 371 for sectors 18 (left) and 22 (right). Top: R-band light curves. Bottom: time-period power spectra. White contours represent variability signatures with a significance > 3σ (> 99.73% confidence level). Gray grids represent the cone of influence regions affected by edge effects.

Moreover, also in line with the ANOVA results, we see significant flux variations on timescales longer than 1 day. However, as mentioned above, we go into more detail on longer timescales in Sect. 4. Hence, the period scale of the IDV wavelet analysis is truncated at 8 days. We also note that, due to the gap present in the complete TESS light curve during sector 19, when the source was not in the field of view, an artificial signature of variability appears in the wavelet representations, which should be ignored as it does not represent a real characteristic variability timescale.

Minute and hour timescale variability has been attributed in the past to changes related to turbulence of plasma in the emitting jet region. This fast variability suggests a very compact emitting region. Therefore, the minimum variability timescale and the causality condition can be used to set a constraint on the size of the region responsible for such emission and variability, given by (see e.g., Aharonian et al. 2007)

$$\begin{array}{c}\hfill R\le \frac{{\tau }_{\min }c\delta }{1+z},\end{array}$$(5)

where τ_min is the variability timescale used for this estimation, c is the speed of light, δ represents the Doppler factor of the emitting region and z is the redshift of the source. We use the Doppler factor δ = 9.5 reported by Chen (2018) and the value of τ_min = 0.47 h estimated here. Under these assumptions, the estimated size of the emitting region responsible for the IDV observed from 3C 371 is R < 6.5 × 10¹⁴ cm (i.e. ∼2 × 10⁻³ pc), roughly a factor 2 more constraining than the estimation performed by Carini et al. (1998). This size is smaller than the general emitting region size typically assumed for blazars, as discussed by Raiteri et al. (2023b) for the case of S4 0954+65. Therefore, the observed variability can likely be ascribed to flux fluctuations within a subregion of the relativistic jet.

Moreover, the minimum variability timescale can also be used to make a rough estimation of the mass of the central black hole as proposed in previous works (Xie et al. 2002a,b; Dai et al. 2007; Gupta et al. 2012)

$$\begin{array}{c}\hfill M_{\mathrm{BH}}\le 1.62\times {10}^4\frac{{\tau }_{\min }\delta }{1+z}{M}_{\odot },\end{array}$$(6)

where M_BH and M_⊙ are the masses of the black hole and the Sun, respectively. Using the aforementioned Doppler factor and τ_min values, we derive an estimate of the black hole mass M_BH ≤ 3.4 × 10⁸ M_⊙, which is consistent with the typical mass ranges of 10⁶ − 10¹⁰ M_⊙ assumed for supermassive black holes in AGN (see for instance Kormendy & Richstone 1995; Richstone et al. 1998), and > 10⁸ M_⊙ for the case of radio-loud AGN (Chiaberge & Marconi 2011). This value is also consistent within errors with different estimations reported by other authors (see e.g., Dai et al. 2007; Chen 2018; Pei et al. 2022).

3.2. Power spectral density analysis

We conducted an evaluation of the PSD associated to the data taken by TESS. PSD analyses in blazars have turned to be extended and effective tools for evaluating the type and origin of the variability, characteristic timescales, connection of the variability with the accretion disk, as well as the nature of the underlying noise (see for instance Bhatta & Dhital 2020; Goyal et al. 2022). For these reasons we calculated the PSD corresponding to each TESS sector. Due to the frequency sampling, the PSD is better sampled at low frequencies with respect to high frequencies. To account for this difference, we binned the derived PSDs with an equispaced binning of 0.06 Hz in logarithmic space, leading to a regularly sampled PSD. All the PSDs are represented in the different panels of Fig. 9.

Fig. 9. Binned PSDs of the different TESS sectors. The dashed black lines represent the broken power law spectral shapes used for modeling the PSDs reported in Table 7.

A common way of modeling the PSD of blazars has been through power law functions (𝒫 ∝ ν^−α, see for instance Bhatta & Dhital 2020; Goyal et al. 2022), where the value of the spectral index provides information of the type of noise dominating the variability of the data (α = −2 corresponding to red noise, α = −1 representing flicker noise and α = 0 indicating white noise). In fact, typical values of blazar PSDs tend to be close to α = −2, indicating a dominant red noise on long timescales. Moreover, when the cadence of the observations allow us to study the shortest timescales, a drastic flattening of the slope is often observed. In these cases, a broken power law function is typically used to describe the PSD, accounting for this spectral index change at high frequencies.

Under these considerations, owing to the regular and high cadence of the TESS data, we modeled the PSD of our data using a broken power law function defined as

$$\begin{array}{c}\hfill \mathcal{P}=\{\begin{array}{cc}A{(\nu /{\nu }_{\mathrm{b}})}^{-{\alpha }_1}:\hfill & \nu <{\nu }_{\mathrm{b}}\\ \hfill \\ A{(\nu /{\nu }_{\mathrm{b}})}^{-{\alpha }_2}:\hfill & \nu >{\nu }_{\mathrm{b}}\end{array}\end{array}$$(7)

where α₁ and α₂ are the spectral indices below and above the break frequency ν_b, and A is the amplitude at ν_b. The need for this spectral shape is also clear from visual inspection of the PSDs represented in Fig. 9, where a change of the spectral index is clearly observable at high frequencies. The best-fit parameters derived are reported in Table 7.

All the PSDs derived from the TESS data are well represented by a broken power law function, where the low-frequency power can be modeled with a red-noise power law (α ∼ 2). At frequencies between 5.69 days⁻¹ and 19.88 days⁻¹ (timescales between ∼4.2 h and ∼1.2 h) we observe a break between the low- and high-frequency regime, where white noise becomes dominant, as observed from the spectral index α ∼ 0 above the break. This flattening and white-noise dominance has been observed in the past for IDV PSDs of blazars on the shortest variability timescales (e.g., Raiteri et al. 2021a,b; Pavana Gowtami et al. 2022), and has been interpreted as a characteristic variability timescale, which is found to be within the range of fast variation signatures reported in the previous section. Therefore, IDV and STV at frequencies below the break can be related to a red-noise like, stochastic nature of the synchrotron variability, most likely due to physical processes in the jet itself within small emitting regions, supporting the scenario where these variations come from a jet subregion (Ryan et al. 2019). Above ν_b, white noise overcomes the power variability of red noise. These breaks could be related with electron escape, cooling, or light-crossing timescales (Ryan et al. 2019).

4. Long-term light curves

Apart from the remarkable very fast variability observed from 3C 371 thanks to the extraordinarily high observing cadence of TESS, this source also displays variations in its emission on longer timescales. Using the 3 year optical light curves provided by the WEBT Collaboration presented in Sect. 2, we study the optical LTV observed from this blazar.

4.1. Variability

During the 3 year monitoring of 3C 371 we can observe different emission states. The first year, the variability is mainly dominated by several flares (see Fig. 3). After this period, on MJD 58500 approximately, the source experienced a significant decrease of its emission, leading to a low and much less variable emission state. After this period, a relatively steady increasing trend in the flux was detected. Finally, a state of fairly constant emission on long timescales was reached after MJD 58900.

As for the IDV, we characterized the amplitude of the observed LTV for the different optical bands. The results are reported in Table 8. The maximum variation on long timescales of the R-band flux density is ΔF_max = 12.65 mJy, with maximum and minimum values of F_max = 14.49 mJy and F_min = 1.84 mJy, and A_mp = (264.4 ± 13.9)%. All the optical bands report similar results, with a slightly lower variability in the I band.

We also calculated A_mp within different time intervals following the same procedure as for the IDV and reported in Table 6. The results for the LTV are presented in Table 9 for the optical R band. As expected from the variability amplitude estimated in different IDV timescales and for the complete long-term light curves, A_mp increases for longer time intervals, meaning that typically, LTV shows much larger amplitudes than IDV.

The broadband emission of blazars often involves several physical processes with different characteristic variability timescales entangled. To evaluate these timescales, we used the wavelet method introduced above, which allows us to identify (and locate in time) periods of significant variation. We represent in Fig. 10 the wavelet transform, with its associated PSD, corresponding to the R-band flux density light curve.

Fig. 10. LTV wavelet analysis performed on the dereddened, host-galaxy-subtracted and flux-calibrated R-band light curve of 3C 371. Top: R-band light curve. Middle: time-frequency power spectrum (left) and time-folded power spectrum (right). White contours in the left panel and the dashed line in the right panel represent variability signatures with a significance > 3σ (> 99.73% confidence level). The gray grid represents the cone of influence region affected by edge effects. Bottom: averaged variance between 1 and 100 days. The horizontal dashed line represents the 3σ (99.73%) confidence level.

However, as already introduced in Sect. 3, blazar time series are typically dominated by the longest variability timescales present in the data, masking the fastest variability signatures (see Raiteri et al. 2021a,b). Therefore, the timescales revealed by this analysis refer to the longest variability timescales present in the emission, dominating over faster variations. This is also visible through the location of the observed significant variability, appearing during the first part of the light curve, coinciding with the flaring period, when the amplitude of the variations is larger.

In order to unveil the effect of variability on shorter timescales we performed a detrending of the long-term R-band light curve following again the procedure employed by Raiteri et al. (2021a,b). Owing to the timescales revealed by the wavelet analysis performed on the unbinned, observed light curve (between ∼100 and ∼300 days), we perform the cubic spline interpolation and detrending using bin widths of 100, 80, 60, 50, 40, 30, 20, 15, 10, 7, 6, 5 and 3 days. For bins smaller than 3 days, the average uncertainties of the data are larger than the standard deviation of the detrended light curve. Therefore, we cannot claim statistically significant variations below ∼3 days with this approach.

Figure 11 shows the time-period space power spectrum of the detrended light curves estimated from different bin widths. It becomes clear by visual inspection the masking effect of long timescales, as shorter significant variability signatures are revealed in the wavelet analysis when shorter temporal bins are used for obtaining the detrended light curve. Also, as the binning used for the detrending procedure becomes smaller, the features of the original light curve are better represented and thus, slower variability signatures are removed. This can be seen from the white contours representing the detected periods of significant variation. For the largest bins (100 days and 50 days), most of the variability appears with timescales > 50 days. On the other hand, as we go to smaller bins, the slowest variations are more accurately removed and most of the significant variations appear at timescales < 30 days. After removing the slowest variations (see binnings < 10 days), several variability timescales are identified. The most remarkable ones are three ∼10−15-day signatures corresponding with the three flares occurring on MJD 58250, MJD 58400, and MJD 58500, approximately, as well as several characteristic timescales of ∼2−3 days, mainly concentrated between MJD 58200 and MJD 58600 associated with the flaring period, but also visible during the smaller flare on MJD 58900. On longer timescales, a rather stable timescale appears at ∼120 days, which is representative of the whole flaring period, as well as a signature at ∼50 days associated with the rise of the first flare covered by our observations.

Fig. 11. Wavelet analysis performed on the detrended R-band light curves. Top left: 100 day binning. Top right: 50 day binning. Middle left: 30 day binning. Middle right: 10 day binning. Bottom left: 7 day binning. Bottom right: 5 day binning. The top panels correspond to the time-frequency power spectra of each light curve. White contours represent variability signatures with a significance of > 3σ (> 99.73% confidence level). The gray grid represents the cone of influence region affected by edge effects. The bottom panels represent the averaged variance between 1 day and the corresponding binning. The horizontal dashed line represents the 3σ (99.73%) confidence level.

We observe that the period of greatest variability occurs approximately between MJD 58200 and MJD 58600, coincident with the intense flaring period displayed by 3C 371 (see Fig. 3). Two periods of almost no significant variability are identified. The first one spans from MJD 58600 and MJD 58800, while the second one starts roughly on MJD 59000 and lasts until the end of the observing period. Finally, in line with the red-noise-like PSD observed for the different TESS sectors, both the amplitude of the averaged variance and the power of the variability decrease as we go to shorter timescales.

4.2. PSD and periodicities

Similarly to the PSD analysis carried out with the data from TESS, we evaluated the PSD corresponding to the LTV of 3C 371 in its optical emission. In Fig. 12 we present the PSD of each band. The best-fit parameters are reported in Table 10. Owing to the much better sampling of the optical R band, we were able to obtain a more precise characterization of the PSD in this band and thus, we were able to model it with a broken power-law function as expressed in Eq. (7). On the other hand, due to the more sparse sampling of the IVB bands, we used a simple power law function 𝒫 = Aν^−α to describe their corresponding PSDs.

Fig. 12. Binned PSDs of the different optical bands. The dashed black lines represent the simple or broken power-law spectral shapes used for modeling the PSDs reported in Table 10.

We observe that the R-band PSD shows a behavior similar to that shown by the PSDs corresponding to the TESS data. The variability is represented by a power law with index α₁ = 1.84, close to a red-noise like spectrum, on small frequencies. A break appears at ν_b = 0.012 days⁻¹, where the spectrum changes to a white-noise dominated regime. For the other three bands we use a power law function to model the PSD. Interestingly, these three PSDs show a smaller index, much closer to the index expected from flicker-noise like power spectrum (α ∼ 1). Flicker noise signatures in blazar PSDs have been associated with long-term variability processes with a certain “memory” of the variations taking place in the accretion disk, that is, disk variations have an imprint in the jet variability (see for instance Ryan et al. 2019; Bhatta & Dhital 2020). However, these three bands differ from the behavior of the R-band PSD, with a better sampling and coverage. Therefore, no strong conclusion can be extracted in this regard with the current data.

Moreover, low-frequency breaks in the PSD of blazars have been detected in several occasions (e.g., Ryan et al. 2019). These breaks have been interpreted as characteristic timescales related with the accretion disk (for instance, diffusive timescale in the outer region of the accretion flow or dynamical or thermal timescales driven by turbulence). Therefore, important information can be extracted from these breaks at low frequencies. However, no evident breaks are observed in the PSDs presented here for our data.

Following the analysis of the power spectrum, we also conducted a periodicity search based on the peaks present in the periodograms/PSDs derived. Having a dataset as regularly sampled as possible is of great importance when performing a periodicity analysis, especially with Fourier transform based techniques (see Otero-Santos et al. 2023b). Therefore, we have conducted this search with the R-band data, for the aforementioned reasons. We employed the widely used Lomb-Scargle periodogram (LSP, see Lomb 1976; Scargle 1982) due to its performance on unevenly sampled data (Peñil et al. 2020; Raiteri et al. 2021a,b; Otero-Santos et al. 2020, 2023b). The significance estimation was performed as detailed in Otero-Santos et al. (2023b). In this calculation, trial factors due to the number of independent frequencies sampled have been taken into account following Vaughan (2005). The resulting LSP is shown in Fig. 13.

Fig. 13. LSP of the R-band light curve of 3C 371. Blue and black solid lines correspond to the LSP and the red-noise model used for the significance estimation, respectively. Colored dotted lines represent the different confidence levels between 1σ and 5σ.

We observe that the LSP presents several peaks, associated with the different characteristic variability timescales displayed by this source. However, none of these peaks have enough statistical significance to claim a possible (quasi-)periodic variation of the optical emission of 3C 371, with the highest peak barely reaching the 2σ confidence limit. Therefore, no signs of periodicities are revealed by our analysis. This is in line with the results presented by Kelly et al. (2003) and Rani et al. (2009) in the radio and X-ray bands, respectively.

5. Optical color variations

In order to study and understand the spectral behavior of 3C 371, we analyzed the color variations between the different optical bands over time. For this color analysis, owing to the different sampling between each band, we select the quasi-simultaneous data taken with the different filters within a range of ±1 h. All the color relations are estimated taking the optical R band as the reference because it has the largest number of observations available. This analysis has been performed on the dereddened, host-galaxy corrected, optical flux densities. In this study, we compute the color indices as the ratios between flux density in each band and the flux density in the optical R band, similarly to the approach used in Otero-Santos et al. (2022). The relations between the different bands are represented in Fig. 14, and the fit parameters are reported in Table 11.

Fig. 14. Dereddened, host-galaxy-corrected color indices of 3C 371 for the different optical bands with respect to the optical R-band flux density. The colour scales represent the MJD of the observations. Black dashed-dotted lines correspond to the linear fits describing the color changes. Left: FR/FI color index. Middle: FV/FR color index. Right: FB/FR color index.

As can be seen from the optical color-magnitude representations and the coefficients of the linear fits shown in Table 11, with values relatively close to zero, this source shows a rather achromatic evolution over long timescales. This behavior has been observed in the past for the LTV of other BL Lac objects (e.g., S4 0954+65, S5 0716+714 or BL Lacertae, see Raiteri et al. 2021a,b, 2023a). These authors interpreted this behavior as changes in the Doppler factor introduced to orientation changes of the emitting region under a helical jet model.

The aforementioned approximately achromatic variability is also observed when considering only the flaring period. On the other hand, the opposite is observed during the two periods of low flux variability (∼2 mJy and ∼4 mJy), associated however to a largely chromatic variation of the optical colors. These periods are coincident with the two low emission states of 3C 371 (dotted-dashed lines in Fig. 15, MJD 58650–58800 and MJD > 59000, approximately). Similar behaviors of very large color variation for a very small flux variability have been reported in the past for other blazars such as Mrk 421 during a high state; and 3C 454.3 or CTA 102 during a low state (see Otero-Santos et al. 2022).

Fig. 15. Dereddened, host-galaxy-corrected color index of 3C 371 for the optical B band with respect to the optical R-band flux density. The plot is presented without error bars for clarity. The color scale represents the MJD of the observations. Black dashed-dotted lines indicate the color change of the two low states. Black arrows highlight some examples of short-term RWB variability patterns.

In contrast, short-term color changes display a pronounced chromatic behavior, as can be clearly seen in Fig. 15, where the F_B/F_R color index is represented as an example without error bars for clarity. More remarkably, this short-term color variability displays a variety of behaviors rather than just the commonly claimed bluer-when-righter (BWB) trend for BL Lac objects. Therefore, mixed BWB and redder-when-brighter (RWB) trends occurring on short timescales can also be identified in the optical regime as signatures of fast variability. The existence of a long-term BWB variability in the color variation, accompanied with mixed BWB and RWB patterns on short timescales is an indication that LTV and STV signatures may be produced or driven by different physical processes (Isler et al. 2017). For instance, Negi et al. (2022) interpret these partial shorter trends as transitions between jet-dominated and disk-dominated states. The rather achromatic flaring state could indicate changes of Doppler factor (Raiteri et al. 2021a,b). Finally, long-term BWB trends are mainly expected from variations intrinsic to the synchrotron emission of the jet.

6. Optical polarimetry

We evaluated the behavior of the polarized emission observed from 3C 371 as part of the WEBT monitoring. First, in order to calculate the polarization fraction intrinsic to the jet, we corrected the data from the depolarizing effect introduced by the host galaxy following Eq. (1). After accounting for this effect, we observe an increase of the polarization degree up to a factor ∼3 during the periods with the faintest emission, when the relative contribution of the host galaxy is highest. The intrinsic polarization degree ranges from 0.05% to 23.46%, with an average value of 8.78%. This results in a variability amplitude A_mp = (266.52 ± 104.72)%, consistent with the variability amplitude observed from the different optical bands. We also evaluated the possible correlation between the variability of the optical emission and the polarization degree of the jet. This can be seen in Fig. 16, where the polarization degree is represented with respect to the optical flux. We estimated the correlation coefficient between P and F_R, resulting in a value of ρ = −0.49 and p-value = 10⁻¹⁰, suggesting a mild anticorrelation. Rajput et al. (2022) observe this behavior for three blazars on long timescales, all of them BL Lac objects. Other examples of this are 3C 345 (Otero-Santos et al. 2023a), BL Lacertae (Raiteri et al. 2013, 2023a), and S4 0954+65 (Raiteri et al. 2023b). In fact, Raiteri et al. (2013) explain such behavior as a cause of a lower Lorentz factor in BL Lacs with respect to FSRQs.

Fig. 16. Intrinsic polarization degree after correcting by the depolarizing effect introduced by the host galaxy versus the derredened, host galaxy corrected, optical flux density.

Concerning the behavior of the EVPA, the angle shows an ∼70° alignment at the beginning of the observing period (see Fig. 4). At the end of the first flare, a period of large variation of the angle orientation takes place, between MJD 58200 and MJD 58600 approximately, were a total angle variation of ∼560° occurs, from −37° to 523°. At the end of this interval, the EVPA shows a rather stable orientation at around 430°, which corresponds to an effective EVPA change of 360° (430°  = 70 + 2 × 180°), with the angle going back to its initial orientation.

We studied this interesting feature of the EVPA evolution during the flaring period in detail. Starting on MJD 58230, the source underwent a long EVPA swing that lasted for roughly 400 days. This event showed a total effective rotation of ∼360°, with a remarkable angle variation that ended with the EVPA aligned with the same preferential orientation observed before this event. In fact, BL Lac objects are known for showing relatively stable and preferred orientations of the EVPA (see e.g., Angel & Stockman 1980). Therefore, this could also be the case of 3C 371, with an intrinsic process temporally disturbing the magnetic field and giving place to this feature before adopting again its original orientation. Similar features have been detected in the past for other blazars such as 3C 279, where two of these events were characterized by Kiehlmann et al. (2016), one during a very low and stable emission state and a second one during an optical flare.

Among other recent works, Kiehlmann et al. (2016) have investigated the possibility of this kind of events happening by chance, due to the stochastic variability of the emission, rather than being caused by an intrinsic process taking place in the blazar. In particular, these authors show that long and slow EVPA rotations can appear through random-walk processes, although less smooth on average. In order to assess the possibility of a physical origin for this event, we followed a similar procedure to that from Kiehlmann et al. (2016). This involved defining a smoothness parameter, s, for the potential rotation (see Eq. (4) in Kiehlmann et al. 2016, for more details), where a curve 1 is considered smoother than a curve 2 if s₂ > s₁. With this definition we can assess the smoothness of the observed feature over similar artificially simulated EVPA rotations under a considered model. Therefore, we calculated the smoothness parameter on the daily binned, corrected EVPA values to avoid biasing s with short timescale changes, in which the errors are expected to be larger than the real variation (see Kiehlmann et al. 2016, for further discussion). The derived value of the smoothness parameter is s = 11.6. Then, using artificial EVPA light curves obtained through MC simulations with the approach presented in Sect. 4, we estimated the smoothness parameter for 10⁴ data series with random-walk, red-noise-like variability, with only 0.49% of the simulations showing an EVPA behavior with the same characteristics as the real data and a smaller smoothness parameter. We performed this test for the simulations before and after the ±180° correction to account for the possibility that this event is artificially introduced by an overcorrection of the EVPA.

We note that our MC simulations are based entirely on the information of the sampling and values of the EVPA. More complex simulations are needed to confirm this result more significantly such as the ones performed by Kiehlmann et al. (2016), with information on the Stokes parameters; these allow us to constrain not only the behavior of the simulated EVPA, but also the polarization fraction and variability within different random walk theoretical models. Nevertheless, we are able to exclude at almost a 3σ confidence level a stochastic behavior with respect to our simple stochastic model.

Several models have been proposed for explaining such features. Kiehlmann et al. (2016) comment on the possibility that a similar feature identified in 3C 279 could be produced by moving plasma along a helical magnetic field, encountering some kind of disturbance such as a shock. In fact, helical magnetic fields have successfully reproduced the long-term variation of the polarization degree and EVPA of blazars in the past (see, for instance, Raiteri et al. 2013, 2021b, 2023b).

7. Flux variability distribution and its origin

Finally, we tested the nature of the flux distribution. Despite the common assumption of a Gaussian distribution of blazar fluxes, several authors have claimed that blazars data series are better described by a log-normal distribution (see for instance Sinha et al. 2018; Shah et al. 2018, 2020). This could be indicative of a nonlinear multiplicative origin of the variability, possibly linked to an imprint of the accretion disk in the jet emission, pointing towards a disk-jet connection. Nevertheless, other interpretations have also been considered in the literature. In order to investigate possible differences between the IDV and LTV signatures and origin, we evaluated the flux distribution of each day observed by TESS individually, and compared the results with the flux distribution of the 3 year R-band WEBT light curve.

Most of the evaluated flux distributions for the observations performed by TESS within each 1 day interval show a good agreement with a Gaussian distribution, with values of χ²/d.o.f. ≲ 1. A log-normal distribution was also tested, and was found to be also compatible with the observed PDFs, with typically no significant preference for one distribution over the other. An example of one of the fitted distributions is shown in the bottom panel of Fig. 17. Exceptionally, we observe that 17 out of the total 318 days show a preference for a log-normal PDF shape. Moreover, several days present multimodal distributions, which are not well-represented by either a single Gaussian nor a log-normal distribution. These time intervals are mostly coincident with the periods of fastest variability and with the most intense flares observed by TESS, especially around MJD 58810 between MJD 58850 and MJD 58920, approximately. For these cases, we find that both bi- (or multi)modal normal and log-normal PDFs are compatible with the observed distributions.

Fig. 17. Examples of IDV flux distributions for three 1 day time intervals. Black solid lines represent the histograms of the distributions. Black dotted and blue dashed-dotted lines correspond to the Gaussian and log-normal fits to the data. Top: flux distribution with both a compatible Gaussian and log-normal distribution. Middle: flux distribution with a preferential log-normal distribution. Bottom: multimodal flux distribution.

These two distributions have been associated with different origins of the observed variability. Intrinsic jet processes are known to be additive in nature and thus, leading to Gaussian PDFs (see for instance Sinha et al. 2018). The good performance of a Gaussian distribution in reproducing the IDV flux distributions may indicate that these variations are produced in a compact region of the jet through linear processes, disconnected from the disk. In fact, IDV is often interpreted within the framework of jet processes (e.g., turbulence or magnetic reconnection, see Narayan & Piran 2012), and are difficult to explain with disk-related processes. Therefore, the results found here on the IDV flux distributions are in line with the IDV being produced by intrinsic jet processes.

We have also studied the 3 year LTV flux distribution using the R-band WEBT data, as shown in Fig. 18. A log-normal function represents this PDF better than a normal one. However, we observe a multimodal distribution, clearly caused by the existence of different emission states, i.e. an intense flaring period until MJD 58600 (P1); a first low-flux and relatively stable period between MJD 58600 and MJD 58820 (P2); a period showing a small flare developing between MJD 58820 and MJD 58950 (P3); and a second rather stable and low emission state starting on MJD 58950 (P4), approximately. Therefore, none of these fits results in a good agreement of the observed flux distribution, with χ²/d.o.f. ≫ 1. We employed a combination of three log-normal distributions to model the PDF peaks observed at ∼2 mJy and ∼4 mJy introduced by the low states, plus the largest flux values due to the high emission state. Alternatively, we also used a composite of three Gaussian distributions. A significant improvement on the representation of the multimodal PDF can clearly be seen in the bottom panel of Fig. 18, still with a preference of the log-normal distribution over the Gaussian (${\chi }_{\log -\mathrm{normal}}^2\sim 3.6$ and ${\chi }_{\mathrm{Gaussian}}^2\sim 47.7$).

Fig. 18. Long-term flux distributions for the dereddened, host-galaxy-corrected flux density R-band light curve. Black solid lines represent the histograms of the distributions. Black dotted and blue dashed-dotted lines correspond to the Gaussian and log-normal fits to the data. Top: single log-normal and Gaussian distribution fits. Bottom: multimodal log-normal and Gaussian distribution fits.

Despite the better representation of the LTV flux distribution when using a multilog-normal function with respect to the single log-normal PDF, we observe from the reduced χ² value that the fit still shows discrepancies with the data, especially for the data with F_R > 6 mJy, i.e. during the high emission state. In order to better model the long-term flux distribution, we also evaluated the flux distribution of each period individually (see Fig. 19). We observe a clear preference of a log-normal distribution over a normal distribution for all four periods, with ${\chi }_{\mathrm{Gaussian}}^2>{\chi }_{\log -\mathrm{normal}}^2$. The two low-state periods and the smaller flare (P2, P3 and P4) have a reasonably good fit, with reduced χ² values of 1.6, 1.9 and 2.2, respectively. The intense flaring period P1 shows a slightly higher reduced χ² of 2.9, mainly due to the large number of measurements with a flux density of ∼6.5−7.0 mJy commented above. Nevertheless, the reduced χ² value for the log-normal distribution is lower than the value obtained from the Gaussian PDF, showing a preference for the former distribution.

Fig. 19. Long-term flux distributions for the different emission states. Black solid lines represent the histograms of the distributions. Black dotted and blue dashed-dotted lines correspond to the Gaussian and log-normal fits to the data. Top left: intense flaring period (P1). Top right: first low state period (P2). Bottom left: small flaring period (P3). Bottom right: second low state period (P4).

The prevalence of log-normal distributions has been mostly associated with processes with an origin that can be traced back to the accretion disk (see e.g., Kushwaha et al. 2016; Bhatta & Dhital 2020; Pininti et al. 2023). This kind of distribution suggests that the mechanism leading to the emission and its variability is multiplicative in nature (Shah et al. 2020), which is the case of accretion disk instabilities occurring at different radii in local viscous timescales, which are later propagated to the jet as multiplicative processes (Bhatta & Dhital 2020). In our case, we observe that the long-term flux density follows a multimodal flux distribution, which has been associated in the past with different activity levels or different dominant contributing emitting regions (for instance Pininti et al. 2023). These could be caused by, for instance, changes in the jet orientation (Raiteri et al. 2021a,b) or changes in the disk activity transmitted to the jet (Shah et al. 2018). Here, different emission states can be clearly identified, leading to the observed multimodality. Evaluating the distribution of each of these emission states, we see that the LTV of each one is better represented by a log-normal distribution. Therefore, LTV may be driven by multiplicative processes, possibly produced in the accretion disk, later connected to the jet (Bhatta & Dhital 2020). Nevertheless, there have been several claims of linear processes that could lead to non-Gaussian flux distributions. Biteau & Giebels (2012) state that a “minijets-in-a-jet” model could resemble a log-normal distribution, considering observational uncertainties, invoked however to explain fast variations down to timescales of hours/minutes. Another interpretation proposed by Sinha et al. (2018) involves linear fluctuations in the underlying particle acceleration and/or diffusive escape rates of the electrons, which could produce nonlinear flux perturbations, possibly explaining the observed log-normality.

8. Summary and conclusions

In this work, we present a detailed study of the optical variability of the BL Lac object 3C 371. We investigated its properties on all timescales over which this source presents significant variability, taking advantage of the large amount of data available from the TESS satellite, with 1 year of two-minute cadence observations, and the extensive optical monitoring performed by the WEBT Collaboration. We summarize the main results in the points below:

1. We characterized the behavior of the IDV signatures displayed by the source, observing an increasing amplitude of the IDV with increasing length of the time intervals considered. We detected significant intraday flux variations with timescales ranging from several minutes to a few hours using the wavelet and ANOVA analyses, with the fastest characteristic timescales being τ_min = 0.47 h. This value served to constrain the size of the emitting region responsible for the observed IDV, R ≤ 6.5 × 10¹⁴ cm, and the mass of the central black hole, M_BH ≤ 3.4 × 10⁸ M_⊙. We also identified the periods showing the fastest significant variations through the use of the wavelet technique, with several remarkable, very fast flaring periods.

2. We evaluated the PSDs for the different TESS sectors (τ < 28 days) – which were found to follow a broken power-law shape –, revealing several IDV and STV signatures in accordance with the variability analysis. The PSD on timescales of days was found to be red-noise dominated, with α ∼ 2. Variations below the break frequency ν_b showed a Gaussian behavior, with the break frequencies found between ν_b ∼ 5.69 days⁻¹ and ν_b ∼ 19.88 days⁻¹. This corresponds to minimum variability timescales of between ∼1 and ∼4 h, similar to those derived from the ANOVA and wavelet analyses.

3. The LTV was also characterized, showing a much larger amplitude of the variability than the IDV. As for the IDV, we identify several characteristic variability timescales, in this case ranging from a few days to timescales of > 100 days.

4. We obtained the PSDs of the different optical bands. The R-band data, for which the best sampling was obtained, show a broken power-law shape with a red-noise like spectrum above ν_b, typical of the LTV of blazars (see e.g., Raiteri et al. 2021a,b).

5. The spectral variability analysis carried out on the different optical bands shows a mild BWB behavior over long timescales; it is almost achromatic, as typically observed in the variability of BL Lac objects (e.g., Raiteri et al. 2021a,b). In addition, when evaluating shorter variations, we observe stronger and different color patters such as short BWB and RWB trends, roughly achromatic variability, or periods of large color change with almost no flux variation.

6. We have analyzed the variability of the polarization degree, revealing a variability amplitude comparable to that of the optical emission. A mild anticorrelation between the polarization degree and the flux is observed. We also characterized the variability of the EVPA, showing a slow and large EVPA rotation during the flaring period. We ran simulations to decipher the probability of this feature having a physical origin. We are able to discard a simple random walk origin of the EVPA variation at a 3σ confidence level based on these simulations.

7. Finally, we evaluated the flux distributions over different variability timescales. We find that on timescales of 1 day, the flux is well represented by a Gaussian PDF, which could indicate that the IDV is produced in a compact region of the jet through additive processes (see Shah et al. 2018). On the other hand, a normal function fails to model the long-term flux distribution, which was observed to be better represented by a multilog-normal PDF. This could be pointing towards a disk-jet connection and a variability produced by multiplicative mechanisms (Kushwaha et al. 2016), where the multimodal nature is due to the presence of different activity levels.

With these results we characterize the variability of the optical emission in 3C 371 at different timescales from intraday (IDV) to short-term (STV) and long-term (LTV) scales. From the observed distributions we derive information about its possible origin: The IDV is in agreement with Gaussian additive processes in a compact region of the jet. LTV is consistent with changes of the Doppler factor, possibly due to orientation changes of the jet, as found in previous works. However the STV shows a more stochastic behavior that may be related to intrinsic jet processes.

Acknowledgments

J.O.S. acknowledges financial support through the Severo Ochoa grant CEX2021-001131-S funded by MCIN/AEI/10.13039/501100011033 and through grants PID2019-107847RB-C44 and PID2022-139117NB-C44; and through grant FPI-SO from the Spanish Ministry of Economy and Competitiveness (MINECO) (research project SEV-2015-0548-17-3 and predoctoral contract BES-2017-082171). C.M.R., M.V. and M.I.C. acknowledge financial support from the INAF Fundamental Research Funding Call 2023. J.A.P. acknowledges financial support from the Spanish Ministry of Science and Innovation (MICINN) through the Spanish State Research Agency, under Severo Ochoa Centres of Excellence Programme 2020-2024 (CEX2019-000920-S). C.M.R. acknowledges support from the Fundación Jesús Serra and the Instituto de Astrofísica de Canarias under the Visiting Researcher Programme 2022-2024 agreed between both institutions. Partly based on data collected by the WEBT collaboration, which are stored in the WEBT archive at the Osservatorio Astrofisico di Torino – INAF (http://www.oato.inaf.it/blazars/webt/); for questions regarding their availability, please contact the WEBT President Massimo Villata (massimo.villata@inaf.it). This article is partly based on observations made with the IAC80 operated on the island of Tenerife by the Instituto de Astrofísica de Canarias in the Spanish Observatorio del Teide. Many thanks are due to the IAC support astronomers and telescope operators for supporting the observations at the IAC80 telescope. This article is also based partly on data obtained with the STELLA robotic telescopes in Tenerife, an AIP facility jointly operated by AIP and IAC. M.D.J. thanks the Brigham Young University Department of Physics and Astronomy for continued support of the ongoing extragalactic monitoring program at the West Mountain Observatory. This research was partially supported by the Bulgarian National Science Fund of the Ministry of Education and Science under grant KP-06-H68/4 (2022) and by the Ministry of Education and Science of Bulgaria (support for the Bulgarian National Roadmap for Research Infrastructure). The data presented here were obtained in part with ALFOSC, which is provided by the Instituto de Astrofisica de Andalucia (IAA) under a joint agreement with the University of Copenhagen and NOT. E.B. acknowledges support from DGAPA-PAPIIT GRANT IN119123. This work is partly based upon observations carried out at the Observatorio Astronómico Nacional on the Sierra San Pedro Mártir (OAN-SPM), Baja California, Mexico. G.D., O.V. and M.S. acknowledge support by the Astronomical Station Vidojevica, funding from the Ministry of Science, Technological Development and Innovation of the Republic of Serbia (contract No. 451-03-47/2023-01/200002), by the EC through project BELISSIMA (call FP7-REGPOT-2010-5, No. 256772), the observing and financial grant support from the Institute of Astronomy and Rozhen NAO BAS through the bilateral SANU-BAN joint research project “Gaia astrometry and fast variable astronomical objects”, and support by the SANU project F-187. This work is partly based upon observations carried out at the Hans Haffner Observatory. This paper includes data collected by the TESS mission, which are publicly available from the Mikulski Archive for Space Telescopes (MAST). Python wavelet software provided by Evgeniya Predybaylo based on Torrence & Compo (1998) and is available at URL: http://atoc.colorado.edu/research/wavelets/. The authors thank the referee for the very careful review of the manuscript and the useful comments.

References

All Tables

All Figures

	Fig. 1. TESS SAP_FLUX light curve of 3C 371. The offsets between the different sectors were estimated using the observed R-band WEBT light curve, and were then converted into arbitrary flux units and normalized by the mean value. The inset panel corresponds to a zoom onto a 1 h of the light curve to illustrate the cadence and uncertainties of the observations.
In the text

	Fig. 2. Finding chart obtained from the PanSTARRS-1 (PS1) Image Cutout Service (http://ps1images.stsci.edu/cgi-bin/ps1cutouts) in r band. The field of view is 8 × 8 arcmin. Stars 1−5 represent the photometric sequence of Xilouris et al. (2006) recalibrated by us (see Table 3).
In the text

Fig. 3. WEBT multiband optical light curves (observed magnitudes) of 3C 371 during the 2018–2020 observing season. Different colors and symbols are used to distinguish the contributing datasets, as specified in Table 2. Vertical lines represent the start of the different TESS sectors during which the source was observed. The R-band calibrated TESS data are represented with gray markers in the third panel.

In the text

	Fig. 4. Optical polarization of 3C 371 during the 2018–2020 period compared to the jet optical flux. From top to bottom: R-band optical light curve (dereddened and host-galaxy-subtracted flux densities), observed polarization degree, intrinsic polarization degree, corrected EVPA behavior. Different colors and symbols are used to distinguish the contributing datasets, as specified in Table 2.
In the text

	Fig. 5. Distributions of the maximum IDV magnitude variations observed from TESS data in different time intervals, which are distinguished with different colors and line styles. The vertical black dashed line represents the mean uncertainty of the data.
In the text

	Fig. 6. Histogram of the minimum variability timescales measured from the variable intervals according to the ANOVA test performed here. The inner panel shows a zoomed-in histogram for timescales of shorter than 6 h.
In the text

Fig. 7. IDV wavelet analysis performed on the dereddened, host-galaxy-subtracted R-band flux light curve of 3C 371. Top: R-band light curve. Middle: time-period power spectrum (left) and time-folded power spectrum (right). White contours in the left panel and the dashed line in the right panel represent variability signatures with a significance > 3σ (> 99.73% confidence level). Bottom: averaged variance between 4 min and 1 day. The horizontal dashed line represents the 3σ (99.73%) confidence level.

In the text

	Fig. 8. IDV wavelet analysis performed on the dereddened, host-galaxy-subtracted R-band flux light curves of 3C 371 for sectors 18 (left) and 22 (right). Top: R-band light curves. Bottom: time-period power spectra. White contours represent variability signatures with a significance > 3σ (> 99.73% confidence level). Gray grids represent the cone of influence regions affected by edge effects.
In the text

	Fig. 9. Binned PSDs of the different TESS sectors. The dashed black lines represent the broken power law spectral shapes used for modeling the PSDs reported in Table 7.
In the text

Fig. 10. LTV wavelet analysis performed on the dereddened, host-galaxy-subtracted and flux-calibrated R-band light curve of 3C 371. Top: R-band light curve. Middle: time-frequency power spectrum (left) and time-folded power spectrum (right). White contours in the left panel and the dashed line in the right panel represent variability signatures with a significance > 3σ (> 99.73% confidence level). The gray grid represents the cone of influence region affected by edge effects. Bottom: averaged variance between 1 and 100 days. The horizontal dashed line represents the 3σ (99.73%) confidence level.

In the text

Fig. 11. Wavelet analysis performed on the detrended R-band light curves. Top left: 100 day binning. Top right: 50 day binning. Middle left: 30 day binning. Middle right: 10 day binning. Bottom left: 7 day binning. Bottom right: 5 day binning. The top panels correspond to the time-frequency power spectra of each light curve. White contours represent variability signatures with a significance of > 3σ (> 99.73% confidence level). The gray grid represents the cone of influence region affected by edge effects. The bottom panels represent the averaged variance between 1 day and the corresponding binning. The horizontal dashed line represents the 3σ (99.73%) confidence level.

In the text

	Fig. 12. Binned PSDs of the different optical bands. The dashed black lines represent the simple or broken power-law spectral shapes used for modeling the PSDs reported in Table 10.
In the text

	Fig. 13. LSP of the R-band light curve of 3C 371. Blue and black solid lines correspond to the LSP and the red-noise model used for the significance estimation, respectively. Colored dotted lines represent the different confidence levels between 1σ and 5σ.
In the text

	Fig. 14. Dereddened, host-galaxy-corrected color indices of 3C 371 for the different optical bands with respect to the optical R-band flux density. The colour scales represent the MJD of the observations. Black dashed-dotted lines correspond to the linear fits describing the color changes. Left: FR/FI color index. Middle: FV/FR color index. Right: FB/FR color index.
In the text

	Fig. 15. Dereddened, host-galaxy-corrected color index of 3C 371 for the optical B band with respect to the optical R-band flux density. The plot is presented without error bars for clarity. The color scale represents the MJD of the observations. Black dashed-dotted lines indicate the color change of the two low states. Black arrows highlight some examples of short-term RWB variability patterns.
In the text

	Fig. 16. Intrinsic polarization degree after correcting by the depolarizing effect introduced by the host galaxy versus the derredened, host galaxy corrected, optical flux density.
In the text

Fig. 17. Examples of IDV flux distributions for three 1 day time intervals. Black solid lines represent the histograms of the distributions. Black dotted and blue dashed-dotted lines correspond to the Gaussian and log-normal fits to the data. Top: flux distribution with both a compatible Gaussian and log-normal distribution. Middle: flux distribution with a preferential log-normal distribution. Bottom: multimodal flux distribution.

In the text

	Fig. 18. Long-term flux distributions for the dereddened, host-galaxy-corrected flux density R-band light curve. Black solid lines represent the histograms of the distributions. Black dotted and blue dashed-dotted lines correspond to the Gaussian and log-normal fits to the data. Top: single log-normal and Gaussian distribution fits. Bottom: multimodal log-normal and Gaussian distribution fits.
In the text

Fig. 19. Long-term flux distributions for the different emission states. Black solid lines represent the histograms of the distributions. Black dotted and blue dashed-dotted lines correspond to the Gaussian and log-normal fits to the data. Top left: intense flaring period (P1). Top right: first low state period (P2). Bottom left: small flaring period (P3). Bottom right: second low state period (P4).

In the text