© The Authors 2025

1. Introduction

Blazars, which are active galactic nuclei (AGN) with one of their jets pointing towards our line of sight, exhibit profuse emission over the entire electromagnetic spectrum. Depending on whether or not they show an optical spectrum that contains features with |EW|> 5 Å in the rest frame (Stickel et al. 1991), they are classified into BL Lacertae (BL Lac) objects or flat-spectrum radio quasars (FSRQs; Urry & Padovani 1995), respectively. Their spectral energy distribution (SED) is known for its double-bump shape (Abdo et al. 2010), owing to the two physical mechanisms typically invoked to explain their emission. There is consensus over the theory that the low-energy bump is produced by synchrotron radiation (see e.g. Königl 1981). The high-energy bump is most commonly interpreted within the leptonic inverse Compton (IC) scenario, either via synchrotron self-Compton (SSC; e.g. Maraschi et al. 1992) or external Compton (EC, see Dermer & Schlickeiser 1993) scattering. However, hadronic processes have gained importance, especially given results suggesting that blazars could be sources of the extragalactic high-energy neutrinos detected by the IceCube Observatory (Aartsen et al. 2018a,b).

The main feature of blazar multi-wavelength (MWL) emission is its remarkable and unpredictable variability in all wavelengths, and on all possible timescales (Wagner & Witzel 1995). Variability on timescales of larger than a few months or years is usually referred to as long-term variability (LTV). Moreover, timescales of the order of several days to weeks are considered as short-term variability (STV). Finally, blazars also display very fast flux variations that take place over timescales of a few hours or even a few minutes, which is known as intraday variability (IDV). These timescales, present in the broadband emission of blazars and thought to be caused by different physical mechanisms owing to their different nature, are typically entangled, which prevents the straightforward interpretation of this observed variability. Therefore, characterising and disentangling the variability signatures of blazars provides a means to understand and interpret the physical origin of the emission observed from these objects.

Here we aim to study the broadband emission of 3C 371, a blazar belonging to the BL Lac object subclass located at a redshift of z = 0.0510 ± 0.0003 (de Grijp et al. 1992). As concluded by Wrobel & Lind (1990), 3C 371 is a BL Lac object observed under a sufficiently large viewing angle that the two radio lobes are visible and do not overlap. Recently, this source has shown remarkable periods of fast IDV in the optical band, which were revealed thanks to observations performed by the TESS satellite over almost 1 year, with a two-minute observing cadence. In parallel, the Whole-Earth Blazar Telescope (WEBT) collaboration monitored the source extensively in the optical and radio between 2018 and 2021. These optical data were analysed in Otero-Santos et al. (2024, hereafter Paper I), where we characterised the IDV and LTV signatures in the emission of 3C 371. Complementing these optical data, several MWL instruments and telescopes observed 3C 371 in this same time window, allowing a broadband interpretation of the emission from radio to γ rays.

In this work, we characterise the broadband variability and SED of 3C 371 using data over the entire spectrum. The paper is structured as follows: In Sect. 2 we describe the MWL datasets and instruments used in the analysis. In Sect. 3 we present the MWL light curves, evaluating the variability and correlations. We present our analysis of the jet features and jet kinematics in Sect. 4, and a spectral analysis of optical-UV, X-ray, and γ-ray bands in Sect. 5. In Sect. 6, we analyse the variability and behaviour of the radio and optical polarised emission, and in Sect. 7 we model the broadband SED and describe our interpretation of the results. Some final remarks and conclusions are provided in Sect. 8.

2. Observations and data reduction

2.1. Gamma-ray observations: Fermi-LAT

We analysed data from the pair-conversion Large Area Telescope (LAT) on board the Fermi satellite. The LAT monitors the γ-ray sky every three hours in the energy range of 20 MeV to ∼2 TeV, and exhibits maximum sensitivity below a few GeV (Atwood et al. 2009). For the present data analysis, a region of interest (ROI) with a radius of 15° was selected, centred at the position of 3C 371. The data between 2018 and 2021 were retrieved from the LAT Data Server¹. These data were analysed by selecting all the Pass8 P8R3 source events between 100 MeV and 300 GeV with version 1.2.23 of the standard software FERMITOOLS. A zenith angle cut of 90° was introduced to reduce the contamination from the limb of the Earth. Moreover, we used the recommended Galactic diffuse emission model and the isotropic component for the event selection performed here².

First, we analysed one year of data with a binned likelihood method to obtain the model that describes the sources in the field. For this analysis, we included not only all the sources within the ROI but also the sources contained in an extra annular region of 10° in radius. The spectral shape and normalisation of the sources within the ROI are left as free parameters in the model, while they were fixed to the catalogue values for those in the annular region. For this, we used the 4FGL-DR2 catalogue (Abdollahi et al. 2020; Ballet et al. 2020). The normalisation of the diffuse components was also left free in the iteration used to obtain the model file. In order to converge to the final model, all the sources with a test statistics (TS) of < 2 were discarded.

To obtain the light curve, we used the resulting model as input for an unbinned likelihood analysis on each 15 day data bin. This bin width was selected due to the relatively faint emission of 3C 371 in the γ-ray band, which led to mostly upper limits on shorter timescales. We assumed a power-law spectral shape for each bin, leaving both the spectral index and normalisation as free parameters. Finally, the parameters of all the point-like sources were fixed, except for five sources located at a distance of < 10° from 3C 371 with TS > 25 and a variability index as reported in the 4FGL-DR2 catalogue of higher than that of 3C 371. The γ-ray flux light curve and the spectral index evolution with time are shown in Fig. 1.

Fig. 1. Fermi-LAT γ-ray light curve of 3C 371. Top: 15 day binned γ-ray flux light curve between 100 MeV and 300 GeV obtained from the Fermi-LAT data. Middle: γ-ray spectral index as a function of time. Bottom: γ-ray spectral index as a function of γ-ray flux.

2.2. X-ray observations: Swift-XRT

The Swift satellite observed 3C 371 in X-rays with its X-Ray Telescope (XRT; see Burrows et al. 2004) between approximately MJD 58682 and MJD 58935. Swift data were reduced with the HEASoft package³ version 6.30.1 with calibration files updated as of November 2022. Observations in pointing mode were processed with the xrtpipeline task, selecting grades 0−12 for the photon counting (PC) mode and grades 0−2 for the windowed timing (WT) mode. In the period of interest, there were 54 PC observations and only two short WT exposures, which we do not consider hereafter. The count rates were low, and so pile-up was not an issue. We extracted the source counts in a circular region with 50 arcsec aperture centred on the source; the background was estimated from a close-by source-free circular region with a radius of 80 arcsec. We used the package Xspec to fit the spectra in the 0.3−10 keV energy range with an absorbed power law. The hydrogen column density was fixed to the Galactic value of N_H = 4.11 × 10²⁰ cm⁻² (Kalberla et al. 2005) with elemental abundances from Wilms et al. (2000).

Because of the low number of counts, in order to have a stability cross-check of the analysis, we used both the Cash statistics and the chi-squared statistics, and in this case we binned the spectra to have at least 20 counts per energy bin. For all the analyses presented in this paper, we used the results of the Cash statistics, because this method is more sensitive than the chi-squared analysis in situations where observations are affected by low number counts, and can recover more flux points (Cash 1979). Figure 2 shows the XRT flux density light curve at 1 keV obtained with both statistics, together with the power-law index Γ, which shows no trend with brightness.

Fig. 2. Swift-XRT X-ray light curve of 3C 371. Top: X-ray flux densities at 1 keV from Swift-XRT obtained with both the Cash (black squares) and the chi-squared (blue triangles) statistics. Middle: Power-law index Γ as a function of time. Bottom: Power-law index Γ as a function of 1 keV flux density.

2.3. Optical-UV observations: Swift-UVOT

The Ultra-Violet and Optical Telescope (UVOT, see Roming et al. 2005) on board the Swift satellite observed 3C 371 during the same period as XRT with three optical filters (v, b, and u) and three UV filters (w1, m2, and w2). A total number of 3, 3, and 10 observations are available for the v, b, and u filters, respectively, in the period considered here. On the other hand, the UV filters have a better coverage, with 28, 30, and 30 observations in the w1, m2, and w2 filters, respectively. The UVOT data were processed using a 5″ aperture radius for the source extraction and an annulus with radii of 15″ and 25″ centred on the source for the background. We derived the source photometry both for each single exposure of the same observation and on the co-added exposures. In the former case, we used the uvotmaghist task, while in the latter case we first checked the image alignment before co-adding the exposures with the task uvotimsum, and finally performed the photometry with the task uvotsource. The resulting light curves in both cases are shown in Fig. 3.

Fig. 3. Swift-UVOT optical–UV light curves (observed magnitudes) in the period considered in this paper. Red plus symbols correspond to photometry on single exposures; blue dots correspond to the results of co-adding the exposures of the same observations.

We also evaluated possible calibration effects in the counts-to-flux conversion for each filter. To this end, we calculated the b − v index for all observations, finding that they deviate from the validity range provided by Breeveld et al. (2010), which was originally calculated for Pickles stars and gamma-ray bursts (GRBs). Therefore, the counts-to-flux conversion factors are not suitable for the spectral type of our source. In order to recalibrate the effective wavelengths and conversion factors for each filter, we followed the prescriptions by Raiteri et al. (2010). We constructed an averaged optical-UV SED, which was later convoluted with the effective areas of the UVOT filters to obtain the new effective wavelengths and counts-to-flux conversion factors. The new effective wavelengths were found to be 5431 Å, 4362 Å, 3478 Å, 2627 Å, 2258 Å, and 2092 Å for the v, b, u, w1, m2, and w2 filters, respectively. Moreover, the counts-to-flux conversion rates for these filters were estimated to be 2.602, 1.467, 1.647, 4.392, 8.356, and 5.968, in units of 10⁻¹⁶, respectively.

In order to correct the observations performed by UVOT from Galactic extinction, we calculated the interstellar extinction in the effective wavelengths reported above using the interstellar extinction curve from Cardelli et al. (1989) and the procedure from Raiteri et al. (2010). We obtained the Galactic extinction values of A_v = 0.095, A_b = 0.124, A_u = 0.150, A_w1 = 0.210, A_m2 = 0.261, and A_w2 = 0.248 for the different bands.

Finally, we corrected the data for the contribution of the host galaxy. To estimate its value in each UVOT filter, we adopted the 13 Gyr elliptical galaxy template by Polletta et al. (2007) and calculated the ratio of the galaxy flux in the UVOT filters with respect to that in the BVRI fluxes. Considering the aperture of 5″ used for the UVOT data analysis, we estimate that ∼35.6% of the total host galaxy emission is contributing to the observed flux. The contributions of the host galaxy to the different optical-UV bands in flux density units are reported in Table 1. These contributions were subtracted from the observed flux density values.

2.4. WEBT optical observations

The optical data presented in Paper I were also included as part of the monitoring performed by the WEBT⁴ on 3C 371. These data include photometric data in the Johnson–Cousins BVRI bands and polarimetric data in the R band provided by several observatories and telescopes around the world, as indicated in Table 1 of Paper I. All the data have been corrected for the reddening effect of Galactic extinction and the contribution of the host galaxy to the total flux in each band was subtracted. Further details of the analysis are provided in Paper I.

2.5. Radio observations

Radio observations coordinated as part of the WEBT observing campaigns have also been performed at 4.8 GHz, 6.1 GHz, 8.5 GHz, 24 GHz, and 230 GHz with several radio telescopes, as listed in Table 2. The flux density calibrated radio light curves are presented in Fig. 4 in comparison with the dereddened and host-subtracted optical R-band light curve in flux density units. The data reduction was carried out following a standard procedure (Teraesranta et al. 1998; Giroletti & Righini 2020).

Fig. 4. Radio light curves of 3C 371 compared with the dereddened and host-subtracted optical R-band flux density (top panel). Different colours and symbols are used to distinguish the contributing datasets, as specified in Table 2.

Visual inspection of Fig. 4 shows that, despite the lower observing cadence of the radio data, the radio emission does not follow the same variability pattern as the optical light curve. While the optical R band goes through a period of fast flares and noticeable variability, the radio emission shows a much slower and less intense variability, with an increasing trend that reaches its maximum when the optical emission is at its minimum. This cannot be extended to the 230-GHz band, where the lower observing cadence prevents us from drawing any strong conclusion as to the evolution of the flux at this frequency. We discuss the behaviour and variability of the radio emission further in Sect. 3.

2.6. Radio Very Large Base Array (VLBA) observations

The Monitoring Of Jets in Active galactic nuclei with VLBA Experiments (MOJAVE)⁵ is a long-term AGN monitoring programme performing Very Large Baseline Interferometry (VLBI) 15 GHz radio observations with the goal of studying the kinematics of extragalactic relativistic jets (Lister et al. 2018). The programme has been providing long-term coverage of AGN jets on parsec scales for over 100 sources since 2002.

We present a set of 18 VLBI ultra-high (milliarcsecond) resolution images of the jet in 3C 371 at 15 GHz taken by the MOJAVE collaboration from August 15, 2019, to February 5, 2021. The visibility data, already fully calibrated both in phase and amplitude for the total flux and linear polarisation, were taken from the publicly available MOJAVE archive⁶. We performed no further calibration of the data. Then, following the same procedure as in Escudero Pedrosa et al. (2024a), we processed the calibrated visibility file from every single observing epoch using difmap (see Shepherd 1997) to produce the I, Q, and U Stokes’ parameter images. We assumed no circular polarisation emission to be produced by the source. We produced the linear polarisation flux and angle maps from the Q and U images, respectively. After mapping, the total-flux brightness distribution of the jet for all the observing epochs was modelled using difmap in the visibility plane to describe every two-dimensional total-flux image as a set of circular Gaussian emission components. Each of these components is characterised by a position, a full width at half maximum (FWHM), and an associated brightness. By identifying every one of the modelled Gaussian components over time, we were then able to follow the time evolution of the most prominent emission regions in the jet of 3C 371 (see Sect. 4). The model fit results from this analysis and the estimates of the radio polarisation are included in Appendices A and B, respectively.

3. Multi-wavelength light curves

3.1. Variability

We evaluated the LTV of the MWL light curves, comparing them to the variability estimated for the optical bands presented in Paper I. To this end, we used the amplitude variability parameter (see Romero et al. 1999), which quantifies the amount of variability shown by each dataset as

$$\begin{array}{c}\hfill A_{\mathit{mp}}(\%)=\frac{100}{⟨m⟩}\sqrt{{(m_{\max }-m_{\min })}^2-2{\sigma }^2},\end{array}$$(1)

where ⟨m⟩ is the mean flux density value of the complete dataset, m_max and m_min are the maximum and minimum measurements, and σ is the observational uncertainty. The uncertainty of A_mp can be expressed as

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

In this equation, σ_⟨m⟩ represents the standard deviation of all flux measurements. The results for the different bands are gathered in Table 3 and are represented in Fig. 5. The values obtained for the optical BVRI bands are also included for comparison.

Fig. 5. Amplitude variability for the different bands as a function of frequency. Different markers and colours represent different bands of the spectrum, as indicated in the legend.

We observe that the difference between the maximum and minimum flux is maximum in the optical and γ-ray bands, which also show the maximum values of A_mp. On the other hand, the variability shown by the rest of the bands is clearly lower, with the caveat that the UV and X-ray data have a much shorter time coverage than the other bands. However, and despite this limitation, the variability of the UV and X-ray bands is closer to that of the optical bands than that of the radio data. The v and b UVOT bands have only three measurements during the considered period, which is insufficient to evaluate their variability. The radio bands clearly show variability of a lower amplitude than that observed for the MWL emission, despite the fact that the data cover approximately the same time span as the optical datasets. This is a common result for γ-ray blazars, which usually show the least variability in the radio bands (see e.g. Ahnen et al. 2017, 2018; Acciari 2023).

3.2. Correlations

Understanding how the emissions in the different bands evolve and behave with respect to each other can provide key information on, for instance, the structure of the jet and the location of the emitting region(s) (see e.g. Max-Moerbeck et al. 2014a; Liodakis et al. 2018, 2019). For this reason, we performed a correlation analysis to evaluate the variation and evolution of the emission in the different bands. To this end, we used the z-discrete correlation function (ZDCF) implementation from Alexander (2013). In particular, we employed the PYTHON package pyZDCF⁷ implemented by Jankov et al. (2022), which is based on the original FORTRAN code developed by Alexander (2013). As reference, we took the WEBT R-band light curve with the best sampling among the long-term light curves presented in Sect. 2, performing a cross-correlation function (CCF) using the light curves corresponding to the different bands with respect to the optical R band emission. Owing to the definition of the CCF, positive (negative) time lags mean that the second data series is lagging (leading) the R-band light curve.

In order to assess the statistical significance of the results, we applied the prescription from Max-Moerbeck et al. (2014b), which is widely used in correlation studies (see e.g. Max-Moerbeck et al. 2014b; Raiteri et al. 2023; Acciari 2023). This method is based on the simulation of artificial light curves with statistical properties that mimic those from the real datasets as closely as possible. In this case, we produced a large number of simulated light curves (10⁴ simulations) using the light-curve simulation method from Emmanoulopoulos et al. (2013), which makes use of the power spectral distribution (PSD) and probability density function (PDF) of the real data to produce the artificial time series, achieving the best solution in terms of statistical resemblance. In particular, we used the PYTHON implementation DELCgen⁸ developed by Connolly (2015) following the instructions provided by Emmanoulopoulos et al. (2013). Finally, for the estimation of the uncertainty of the measured time lag, we implemented the approach of Peterson et al. (1998, 2004). This formalism is based on an independent Monte Carlo (MC) randomisation and a random subset selection, deriving the time lag of maximum correlation as the centroid r_cen of all the significant correlations (above a 2σ threshold) that also fulfill r > 0.8 R_max, and its uncertainty as the 68% confidence limit of the estimation of r_cen.

With these considerations, we observe that the optical and UV bands show a high degree of correlation, which is compatible with zero time lag, in agreement with the literature. All the ZDCF curves are represented in Fig. 6. The results of the correlation analysis, with the estimated correlation coefficients, time lags, and significance are reported in Table 4. We note that, due to the poor sampling of the optical filters of Swift-UVOT, we only performed the correlation analysis for the three UV light curves.

Fig. 6. ZDCF of the radio, optical, UV, X-ray, and γ-ray light curves with respect to the optical R-band observations. From top to bottom: R − γ rays (left), R − 1 keV (right), R − w2 (left), R − m2 (right), R − w1 (left), R − B (right), R − V (left), R − I (right), R − 24 GHz (left), R − 8.5 GHz (right), and R − 4.8 GHz (centre). Different dotted lines denote the contours corresponding to a statistical significance of 1σ, 2σ, and 3σ.

The X-ray emission seems to show a delay. However, the sparseness of the data and the large uncertainties due to the faintness of the source in the X-ray domain result in a large uncertainty on the lag estimation, making this delay still compatible with zero. Also, the low level of statistical significance stemming from the same factors prevents us from extracting more conclusive results for this band.

We also performed the inter-band correlations between the R-band emission and the different radio frequencies. In particular, we looked for possible correlations with the 4.8 GHz, 8.5 GHz, and 24 GHz light curves. Fewer observations are available for the 6.1 GHz and 230 GHz radio bands, and so no correlations can be obtained at these frequencies. We again observe time lags compatible with zero for the 4.8 GHz, 8.5 GHz, and 24 GHz bands. However, in this case, the radio–optical emission seems to be anti-correlated, with significance ranging between 2.5σ and 2.8σ depending on the frequency. In addition to the anti-correlations observed at zero time lag, the 8.5 GHz and 24 GHz bands also show an increase in the correlation degree at longer delays of ∼400 days, albeit with a rather low significance (∼2σ). This behaviour has typically been observed for the radio emission of blazars in the literature (see e.g. Ramakrishnan et al. 2016). However, the lower coverage of the radio data analysed here prevents us from making conclusions related to long-term behaviour, and the anti-correlation observed could also be an effect of variability on timescales inadequately sampled by our analysis.

Finally, for the γ-ray emission, there is a very faint hint of correlation compatible with zero lag; however, with only a 1.7σ significance. Also, the larger bin of the γ-ray light curve (and therefore lower number of data points) may be affecting the results found here due to a smoothing of the light-curve features. Ramakrishnan et al. (2015) did not find significant evidence of a correlation between the radio and γ-ray emission of 3C 371. We evaluated the 8.5 GHz radio and γ-ray correlation in the time period considered here, and note the same behaviour reported by these authors. The correlation curve between these two bands is shown in Fig. 7.

Fig. 7. ZDCF between the 8.5 GHz radio and the γ-ray-band light curves. Different dotted lines denote the contours corresponding to a statistical significance of 1σ, 2σ, and 3σ.

4. Jet kinematics

We performed a kinematic analysis of 2 cm VLBA images of 3C 371 from the MOJAVE monitoring program. The images were first fully calibrated in phase and amplitude by the MOJAVE collaboration. After downloading the calibrated data, we used difmap to produce a VLBA total flux and polarisation image for each of the 18 observing epochs used in this study. The most prominent features were then fitted to circular Gaussian components using difmap and were visually cross-identified along the 18 epochs. Following the approach presented by Weaver et al. (2022), we aim to characterise the variability timescale of the components by considering that the flux of each knot or component evolves as a decaying exponential function F = F₀e^−t/t_var. This commonly used assumption allows us to estimate the variability timescale as the time at which the flux of a component decays by a factor e. Among all nine identified components, highlighted in Fig. 8, only one of them (B4) was adequately fitted as a ballistic component with an exponential flux decay (see Fig. 9), allowing us to extract its speed and timescale of variability. From these estimates, the Doppler factor and apparent speed could be computed (Jorstad et al. 2005) together with the bulk Lorentz Factor and viewing angle. For this component, we obtain an apparent speed of β_app = 5.9 ± 0.8 in units of the speed of light, a Doppler factor of δ = 6.0 ± 1.1, a bulk Lorentz factor of Γ = 6.0 ± 1.8, and a viewing angle of θ = 9.6 ± 1.6 deg. The uncertainties on these parameters were computed through standard error propagation using the uncertainties in the Gaussian model fitting, as given by difmap, and the uncertainties associated to the fits presented in Fig. 9, following the procedure indicated in Weaver et al. (2022). An independent preliminary kinematics analysis performed by the MOJAVE team also revealed superluminal motions within the jet with β > 4, which is consistent with the analysis performed here (Kovalev & Homan, private communication). These results are also compatible with the Doppler factor estimated by Homan et al. (2021) –based on data obtained up to 2019– by comparing brightness temperatures and proper motions.

Fig. 8. Epoch 2020 Jan. 4 VLBI 2 cm image of 3C 371, showing the total flux density (black contours) and the polarised flux density (colour scale). Black line segments overlaid in the image represent the EVPA orientation. The innermost and outermost contours correspond to flux density levels of 5.5 and 0.0015 Jy, with the rest of the contours logarithmically equispaced. The green ellipse represents the beam size. The image highlights the core A0 and all identified components B1–L0.

Fig. 9. Results of the kinematic analysis for component B4. Top: Distance of B4 from the jet core A0 over time. Bottom: Flux density of B4 over time.

Components B1–B3 were also observed to decay in flux; however, their speed is too low to accurately trace their positions and perform the fit. This is in line with previous studies by MOJAVE on this source, where no superluminal motion of the identified components was observed (see e.g. Lister et al. 2019, 2021). Nevertheless, we note that these studies include all MOJAVE data before the period considered here. The variability of this source was rather low up until 2018, at which point it underwent an unprecedented period of several outbursts and large variations that could be related to the propagation of material along the jet. This propagation of material and acceleration along the jet has been reported for several sources in the past by Homan et al. (2015). Moreover, this brightening and larger variability can also be observed in the VLBI images used here with respect to previous epochs (see Lister et al. 2019), with a clear increase in the core flux contemporaneous to the movement observed in the identified components. In addition, the ejection and motion of material along the jet is often observed in blazars after a strong flaring period such as the one observed in 3C 371 between 2018 and 2019 (see e.g. Liodakis et al. 2020; Lico et al. 2022; Escudero Pedrosa et al. 2024a).

We compare the total 15 GHz flux obtained from the VLBI data analysis shown in Fig. 10 with the radio flux at 8.5 GHz measured by the Svetloe and Medicina radio telescopes. Despite the shorter coverage, the trend followed by the 15 GHz data is correlated with that seen in the other radio bands. We calculate a correlation coefficient between the simultaneous 8.5 GHz and 15 GHz data of ρ = 0.66-pagination (p-value = 0.001).

Fig. 10. Comparison of radio flux at 15 GHz and 8.5 GHz. Top: Total flux at 15 GHz derived from the MOJAVE VLBI images. Bottom: 8.5 GHz radio flux measured by the Svetloe and Medicina radio telescopes.

Finally, the more external components B5–L0 were clearly observed to move faster from the core. However, their flux does not resemble a decaying exponential, most likely due to the fact that they are too dim as reflected by their large uncertainties. Therefore, a reliable kinematic analysis and an estimation of their physical parameters were not possible.

5. Spectral variability

5.1. Optical–UV colour variations

In Paper I, we evaluated the optical colours with respect to the optical R band. These colour relations showed a mild –almost achromatic– long-term bluer-when-brighter (BWB) behaviour, with mixed BWB and redder-when-brighter (RWB) trends on shorter timescales. Here, we extend the analysis to the optical–UV colours between the optical R band and the w1, m2, and w2 filters.

In Fig. 11 we present the optical–UV colour relations derived from the data used here. The colour is obtained as the flux ratio between the R band and each UV band, as in Paper I. The temporal coverage provided by Swift-UVOT is much shorter than that from the optical R-band light curve. However, the long-term evolution seems to follow the mild BWB trend observed between optical bands and with a relatively small colour change, as shown by the low values of the linear fit slopes presented in Table 5. Nevertheless, more data in the UV band are needed to strongly claim a clear long-term colour variability. In addition, as observed from the optical colours, on shorter timescales we can observe a RWB pattern, especially in the m2 − R colour relation.

Fig. 11. Dereddened, host-galaxy-corrected colour indices of 3C 371 for the different UV bands with respect to the optical R-band flux density. The colour scales represent the MJD of the observations. The black dash-dotted lines correspond to the linear fits describing the long-term colour changes. Left: Fw1/FR colour index. Middle: Fm2/FR colour index. Right: Fw2/FR colour index.

This behaviour follows the results observed in Paper I for the optical colour variability. A very mildly chromatic long-term behaviour is observed, which can be explained as changes of the Doppler factor, possibly caused by orientation changes in the jet that can be ascribed, for instance, to the presence of a helical jet (see e.g. Raiteri et al. 2021a,b). On the other hand, a variety of colour trends on shorter timescales can be associated to intrinsic jet processes linked to the low-energy synchrotron emission of blazars (Raiteri et al. 2021a,b). Another possibility involves different relative contributions between the synchrotron emission of the jet and that from the disc. In particular, the RWB trend clearly visible between the R and m2 bands may reveal the presence of a faint accretion disc, as this behaviour can be observed during a low flux period (Gu et al. 2006; Gu & Ai 2011; Isler et al. 2017; Otero-Santos et al. 2022). The m2 band is also the one showing the highest colour variation of the three UV wavelengths, as seen from the coefficients of the linear fits, although all three remain compatible within the errors. In such cases, the RWB colour evolution could be explained as a relatively stable disc contribution that peaks in the UV band and becomes visible when the synchrotron continuum is faint enough that it no longer completely hinders this thermal contribution, therefore producing a RWB trend until the synchrotron radiation contribution outshines that from the disc at such frequencies (see e.g. Villata et al. 2006; Isler et al. 2017). However, this scenario is disfavoured given the BL Lac nature of the source; indeed this blazar type typically has inefficient, faint accretion discs that are highly subdominant with respect to the synchrotron emission.

Another possibility is that there are two physical components, a flaring component and a quiescent component, or an underlying flow and a propagating disturbance along the jet. The variability introduced by such a propagating component could manifest on shorter timescales –leading to the flaring period that can be observed in the optical light curve– and would be highly chromatic, while that from the underlying flow is expected to be smaller and slower, and might be driven by Doppler factor changes, producing the observed almost achromatic long-term behaviour. This scenario would be consistent with the propagation of component B4 identified in the VLBI data analysis, as the period during which this component is identified is approximately coincident with the strong RWB trend observed in the colour analysis.

5.2. Optical spectral variability

Apart from the optical–UV colour variation analysis between the different optical bands performed here and in Paper I, we also constructed a set of optical spectra with our entire simultaneous data set in order to test its spectral shape and variability over the analysed period. To this end, we converted the flux densities in milliJanksys (mJy) to units of SED (νF_ν; i.e. erg cm⁻² s⁻¹). We constructed a series of SEDs, binning the optical data taken within an interval of 1 day, allowing us to model the variability of the optical synchrotron emission over timescales of days, producing a total of 82 SEDs. The constructed SEDs are represented in Fig. 12.

Fig. 12. Optical SEDs constructed with quasi-simultaneous BVRI optical data. The colour bar indicates the MJD of the corresponding SED.

The SEDs show that flux varies by almost one order of magnitude. Upon visual inspection, we observe that the frequency of the synchrotron peak is located in the red/near-infrared part of the spectrum, in agreement with the values of ν_sync = 4.9 × 10¹⁴ Hz and ν_sync = 7.7 × 10¹³ Hz reported in the third and fourth catalogues of active galactic nuclei constructed based on data provided by the Fermi-LAT telescope (Ackermann et al. 2015; Ajello et al. 2020). We tested a power-law spectral model in the form of νF_ν = K ⋅ ν^−α to evaluate whether or not there is any relation between the optical spectral index and the flux of the source. We find that the fitted spectral index varies within α_PL = 0.33 − 0.90. Moreover, despite the sometimes large errors in the estimate of the spectral index, we derived a correlation coefficient of ρ = −0.58 and a p-value = 10⁻⁸, indicating a harder-when-brighter evolution of the optical spectrum with a rather small change of the spectral index, which is therefore compatible with the mild BWB trend found in the optical colour analysis. This pattern can be observed in Fig. 13. Owing to the reduced spectral range used for this analysis, we did not consider other, more complex spectral models (e.g. log-parabola function). Xilouris et al. (2006) also reported a hardening of the spectrum with increasing flux variations (and the opposite for decreasing flux changes) for 3C 371 during different nights, which is consistent with the results reported here.

Fig. 13. Optical R-band flux with respect to the estimated spectral index using a power-law model. The colour bar indicates the MJD of the corresponding SED.

5.3. X-ray spectral variability

We also investigated possible X-ray spectral variability using the Swift-XRT data. We fit the X-ray spectral index light curve presented in Fig. 2 to a constant value to decipher whether or not the source displays significant spectral change over the observing period, and evaluated the fit using a χ² test. The resulting value of χ²/d.o.f. = 0.01 suggests there is no significant variation of the X-ray spectral index and that the light curve is fully compatible with being constant over the period analysed. The fitted constant X-ray spectral index is Γ = 1.63 ± 0.19. This is also observable from the relation between the X-ray spectral index and the X-ray flux at an energy of 1 keV, shown in the bottom panel of Fig. 2. Several blazars have shown correlated behaviours between the X-ray flux and spectral index in past studies, with harder-when-brighter trends in some cases (see e.g. Aleksić et al. 2013; Wang et al. 2018; Acciari 2023). However, even if present in this case, the large errors of the X-ray spectral index of 3C 371 prevent us from drawing any strong conclusions in this regard, as the spectral index reported here is compatible with a constant value.

5.4. γ-ray spectral variability

We evaluated the γ-ray spectral variability using the 15 day binned Fermi-LAT light curves. We model the spectrum in each 15 day light curve bin with a power-law function, as shown in Fig. 1. This power-law model is consistent with a constant spectral index of α = 2.31 ± 0.03 (χ²/d.o.f. = 1.02) over the period considered here. This spectral index is in agreement with the value of α_4FGL = 2.32 ± 0.02 reported in the 4FGL-DR4 catalogue after averaging 14 years of LAT data. Therefore, we observe that the γ-ray spectral index is compatible with a constant and rather steep spectrum, and no significant spectral variability in the γ-ray band is detected on timescales of 15 days.

Moreover, we compute the average SED using the 3 years of Fermi-LAT data as specified in Sect. 2.1. Owing to the much larger integration with respect to the 15 day light-curve bins, we test a power law and a log-parabola spectral shape, calculating the possible preference of the log-parabola over the power-law model through a likelihood ratio, as

$$\begin{array}{c}\hfill {\mathrm{TS}}_{\log -\mathrm{parabola}}=-2(ln{\mathcal{L}}_{\mathrm{power}\mathrm{law}}-ln{\mathcal{L}}_{\log -\mathrm{parabola}}),\end{array}$$(3)

where ℒ is the likelihood of the power law and log-parabola in each respective case. The parameters of both spectral models are reported in Table 6. The spectral parameters of both the power-law and the log-parabola models are consistent with those reported in the 4FGL-DR4 catalogue: α_4FGL = 2.32 ± 0.02 for the power law; α_4FGL = 2.27 ± 0.02 and β_4FGL = 0.048 ± 0.010 for the log-parabola. For a three-year integration between 2018 and 2020, we only see a 2.53σ preference for the log-parabola curvature over the simpler power-law model. This result is nevertheless not surprising, as the 4FGL-DR4 catalogue has a 5σ preference for a curved spectrum only after 14 years of integration, which amounts to more than four times the amount of data used in this work.

Moreover, we compute the γ-ray spectrum for the highest and lowest emission during the period coincident with Swift observations, which is the time interval with complete MWL coverage used in Sect. 7 for a broadband interpretation of the emission. For these models, we used time windows of 60 and 90 days for the high and low states, respectively. These integration intervals were chosen in order to obtain a well-characterised SED, as shorter time intervals led to spectra and SEDs with only upper limits. The results for both spectral models are again reported in Table 6. Owing to the short time intervals considered, the use of a power-law shape instead of the log-parabola is justified, as no significant curvature is expected. This is proven by the spectral analysis performed here, that shows a hint of a curved spectrum observed after integrating 3 years of data.

6. Polarised emission

This source shows a rather low radio polarisation degree at 15 GHz, with values varying between 0.15% and 1.91%, as shown in Fig. 14. As done in Paper I for the optical emission, we evaluated the correlation between the 15 GHz flux density and the 15 GHz polarisation degree. We find that the radio polarisation degree is anti-correlated with the radio flux density, as reflected by the correlation coefficient of ρ = −0.52 (p-value = 0.016). Therefore, we observe the same behaviour as identified for the optical emission and observed in the past for other sources (for instance BL Lacertae or S4 0954+65; see Raiteri et al. 2013, 2023, respectively), where it was interpreted as a cause of the lower Lorentz factor typically measured for BL Lac objects with respect to FSRQs (Raiteri et al. 2013).

Fig. 14. Radio and optical polarisation of 3C 371. From top to bottom: 15 GHz flux density, polarisation degree at 15 GHz, EVPA at 15 GHz, intrinsic optical R-band polarisation degree, and optical R-band EVPA.

Regarding the EVPA at radio frequencies, during the first epochs covered by the VLBI data there is a certain scatter of the EVPA values between 85° and 34°, contrary to the stability shown by the EVPA in the optical R band in the same period. However, the last seven epochs show a rather constant evolution of the EVPA, with a value of ∼80°. This value is compatible with the polarisation angle orientation of 430° in the optical R band after the slow-rotating feature reported in Paper I, considering the difference introduced by the complete 360° rotation (430°  = 70 + 2 × 180°), meaning that both the radio and optical EVPA end up aligned in approximately the same direction after the studied period. This delay in the stabilisation of the EVPA in the radio band could also indicate a delay of the features in radio frequencies with respect to other frequencies, as mentioned above, meaning that the emission at optical and radio frequencies comes from different regions of the jet, as already suggested by the results obtained from the correlations. Unfortunately, there is no information available on the radio polarisation during the slow rotation observed in the optical polarised emission with which to evaluate the behaviour in the 15 GHz band during this event.

7. SED modelling

Finally, we interpret the MWL emission of 3C 371 within the context of the emission processes occurring in blazar jets, providing information on the physical parameters describing the jet and the particle populations responsible for this emission. In order to compare and interpret the observed variability, we construct the SED at the lowest and highest states in which we have complete MWL coverage of the source, that is, during the period in which Swift-XRT and Swift-UVOT data are available. The high-state SED is produced with the MWL data closest in time to MJD 58921, where the maxima of the I- and R-band fluxes during this period are observed. On the other hand, the low-state SED is based on data obtained as closely as possible in time to MJD 58726, following the same criteria for the lowest flux observed. The γ-ray spectrum of each state is obtained by integrating a total of 60 days around the defined MJD for the case of the high state, and 90 days in the case of the low state, as already mentioned in Sect. 5.4. We corrected the γ-ray spectra for the absorption introduced by the extragalactic background light (EBL). Even if this effect is small in the energy range of Fermi-LAT and especially at such a low redshift, we consider the EBL absorption of the γ-ray SED using the EBL model from Saldana-Lopez et al. (2021) in order to model the broadband SED in the most accurate and correct way.

The modelling was performed using the JetSeT⁹ software (Tramacere et al. 2009, 2011; Tramacere 2020). We interpret broadband emission within a one-zone (blob) leptonic scenario, where the low-energy emission is produced by synchrotron radiation and the high-energy emission is mainly the product of SSC scattering of low-energy photons, with an electron population described by a power-law function with an exponential cut off with the Lorentz factor, as

$$\begin{array}{c}\hfill n(\gamma )=N{\gamma }^{-p}exp\frac{\gamma }{{\gamma }_{\mathrm{cut}}},{\gamma }_{\min }\le \gamma \le {\gamma }_{\max },\end{array}$$(4)

where γ_min, γ_max, and γ_cut are the minimum, maximum, and cut-off values of the Lorentz factor, p is the spectral index of the electron distribution and N is the normalisation. Moreover, the emission depends on a set of parameters including the size of the emitting region R, the magnetic field B, the bulk Lorentz factor Γ, and the viewing angle θ of the jet.

Apart from the emission of the jet, we also consider the presence of a dusty torus (DT), a broad line region (BLR), and an accretion disc that can contribute in the frequency range from infrared to UV. We used the measurement of the Hβ broad emission line luminosity of 3 × 10⁴⁰ erg s⁻¹ reported by Torrealba et al. (2012), and considering that the luminosity of this emission line corresponds to approximately 22/556 of the entire BLR (see Francis et al. 1991; Ghisellini & Tavecchio 2015), we estimate that the BLR has a total luminosity of L_BLR ∼ 8 × 10⁴¹ erg s⁻¹. Using the approximation that relates disc and BLR luminosity as L_disc ≃ 10 × L_BLR (Ghisellini & Tavecchio 2015), the accretion disc luminosity can be approximated as L_disc ∼ 8 × 10⁴² erg s⁻¹. We fix the accretion efficiency of the disc and the covering factors of the DT and BLR to standard assumed values of η = 0.08, τ_DT = 0.1, and τ_BLR = 0.1 (see e.g. Smith et al. 1981; Ghisellini & Tavecchio 2015). The radii of the BLR and DT are fixed with respect to the disc luminosity using the relations available from Kaspi et al. (2007) and Cleary et al. (2007), respectively, as

$$\begin{array}{cc}\hfill R_{\mathrm{BLR},\mathrm{in}}& ={10}^{17}{\left(\frac{L_{\mathrm{disc}}}{{10}^{45}}\right)}^{0.5},\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill R_{\mathrm{BLR},\mathrm{out}}& =1.1\times R_{\mathrm{BLR},\mathrm{in}},\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill R_{\mathrm{DT}}& =2.5\times {10}^{18}{\left(\frac{L_{\mathrm{disc}}}{{10}^{45}}\right)}^{0.5}.\hfill \end{array}$$(7)

The presence of the BLR, DT, and disc also introduces the possibility of having an injection of low-energy photons from an external photon field to the jet that can ultimately give rise to significant high-energy emission via EC scattering. Therefore, in our models, we also consider this process when interpreting the high-energy peak of the SEDs. Moreover, as these components are more visible when the synchrotron emission is fainter due to the lower dominance of the synchrotron radiation over their thermal contribution, we first model the SED of the low state, and use the derived values of the DT and disc temperatures, T_DT and T_disc, as fixed parameters in the modelling of the high-state SED.

We also consider a second component designed to model the extended jet radio emission, that is, the emission below the synchrotron self-absorption (SSA) frequency of the blob. We model this emission as a power law with an index α_radio, starting from ${\nu }_{\mathrm{SSA}}^{\mathrm{radio}}$ and with an exponential cut-off at ${\nu }_{\mathrm{cut}}^{\mathrm{radio}}$. Below the ${\nu }_{\mathrm{SSA}}^{\mathrm{radio}}$, the spectral index is set to 5/2. This component mimics an extended jet, where the window between ${\nu }_{\mathrm{SSA}}^{\mathrm{radio}}$ and ${\nu }_{\mathrm{cut}}^{\mathrm{radio}}$ is produced by the envelope of the optically thin synchrotron emission from different layers of cooled particles, where the age of the particles increases with the distance from the ‘blazar region’, and the window below ${\nu }_{\mathrm{SSA}}^{\mathrm{radio}}$ marks the region dominated by the SSA emission from these layers. We stress that this model is purely phenomenological, but nevertheless the inferred value of ${\alpha }_{\mathrm{SSA}}^{\mathrm{radio}}$ could provide interesting information about the magnetic field topology of the jet, if we assume a conical geometry for example. Nevertheless, this is beyond the scope of the present analysis and we aim to investigate this topic in a follow-up analysis. As seen from the correlation analysis, the emission at radio wavelengths does not show a significant correlation at zero lag as in the rest of the bands, either indicating a possible long delay or a lack of correlation with the rest of the MWL data. This can be explained, for instance, in terms of a strong self-absorption at radio wavelengths, with this radiation being emitted from a different region of the jet at a larger distance from the central black hole, where the radio opacity is lower (see e.g. Max-Moerbeck et al. 2014a; Acciari 2023). Another alternative is that the radio emission is observed after the emitting region has propagated along the jet and has adiabatically expanded (Tramacere et al. 2022).

In addition, and in order to simplify and constrain the model and the number of free parameters, we relate the values of γ_max and γ_cut through their ratio, allowing values of γ_cut/γ_max of between 0.001 and 0.5. The size of the emitting region R is constrained through the opening angle of the jet θ_open, assuming a value of 5°, as

$$\begin{array}{c}\hfill R=R_H\times tan{\theta }_{\mathrm{open}},\end{array}$$(8)

where R_H is the distance of the region from the central black hole, and is fixed to 10¹⁸ cm. We also restrict the range of variability of the viewing angle and Lorentz factor of the jet, θ and Γ, to values around those derived from the VLBI analysis. For the radio component, we fix the value of ${\alpha }_{\mathrm{SSA}}^{\mathrm{radio}}$ from the two radio data points with the lowest frequency to the value of ${\alpha }_{\mathrm{SSA}}^{\mathrm{radio}}=-0.35$, and we fix ${\nu }_{\mathrm{SSA}}^{\mathrm{radio}}=1\times {10}^6$ Hz, that is, without loss of generality, below the lowest-frequency radio data point. These values were crosschecked with those obtained from the archival data available at the Space Science Data Center (SSDC)¹⁰ database and were also found to be relatively constant around the aforementioned values using more than two decades of data. Therefore, we can reliably assume these values to be valid for both models. We then leave the cut frequency ν_cut and νF_ν, p as free parameters. Finally, we assume one cold proton per relativistic electron in the jet. We then performed a model fitting using the JeTSet ModelMinimizer module with an iminuit minimisation (Dembinski et al. 2020). The fit is performed considering the components and constraints mentioned above, adding a systematic error of 20% to the radio, optical, and UV data of the high- and low-state SEDs to avoid biasing the resulting models due to the low statistical uncertainties of these data with respect to the X-ray and γ-ray spectral points. The parameters describing each model are reported in Table 7. The two resulting models are shown in Fig. 15. A detailed discussion on the reliability of the models and the uncertainties of the parameters is included in Appendix C.

Fig. 15. Broadband SED models for the high- (left) and low-brightness states (right) of 3C 371. The top panels contain the broadband SEDs, with the contribution of each component and the best-fit model being represented with different markers, lines, and colours, as indicated in the legend. The bottom panels represent the residuals of the fitted models. We note that the EC scattering contribution via low-energy photons from the BLR was also considered. However, it is even fainter than that from the DT, and therefore does not fall within the scale of the figures and is not included.

The broadband emission is dominated by the non-thermal synchrotron continuum emission, with very faint contributions of the dusty torus and the accretion disc. The accretion disc makes a very minor contribution to the UV emission during the low state, which may support the slight increase in flux in the w1 and m2 bands when the synchrotron emission is faintest. However, we note that this small bump may also be biased by the fact that the data are not strictly simultaneous due to the poorer time coverage of Swift-UVOT. Nevertheless, as expected from BL Lac objects, this contribution is very small even when the synchrotron emission is at its lowest (see e.g. Wang et al. 2002; Plotkin et al. 2012). Moreover, also as expected from BL Lac blazars, the contribution of the EC from the torus, BLR, and disc photon fields is also negligible, as observed in Fig. 15 (e.g. Gao et al. 2011). In fact, the BLR EC contribution is even below the limits of the figure scale, and is therefore not represented in the plots. As a comparison, we included in Appendix D the SED models corresponding to a SSC model, ignoring the EC contribution. These characteristics indicate that this source is indeed a BL Lac object, observed in this case under a rather high viewing angle, as derived from both the VLBI analysis (θ = 9.6 ± 1.6 deg) and the SED models (θ ∼ 10.5 − 11 deg), and as reported in the past by Wrobel & Lind (1990). The values derived for B and Γ are consistent with those reported by Ghisellini et al. (2010) for a larger sample of blazars, combining both BL Lacs and FSRQs.

The variability between the high- and low-emission states can be explained in terms of changes to the distribution of particles –mainly through a hardening of the electron distribution during the high state– as well as a slight increase in the magnetic field B and the bulk Lorentz factor Γ. This is consistent with the short-timescale variations being caused by intrinsic jet processes –such as changes of the electron population– and LTV being attributed to changes of the Lorentz factor. These changes in the electron population, as discussed above based on the colour analysis results, can be due to an injection of plasma with a harder energy distribution, as indicated by the different colour trend and its coincidence with the detection of B4 in the VLBI analysis. Regarding the Doppler factor changes responsible for the LTV, as previously discussed, there may be more than one possible cause. A plausible mechanism could be related to orientation changes of the jet or the emitting region, either based on, for instance, helical (Villata & Raiteri 1999; Ostorero et al. 2004; Raiteri et al. 2021a,b) or precessing jets (see Escudero Pedrosa et al. 2024b). These scenarios have been confirmed in a handful of sources (e.g. Raiteri et al. 2017; Cui et al. 2023; Fuentes et al. 2023). However, Doppler factor variations based on intrinsic Lorentz factor changes of the underlying flow or the component responsible for the LTV could also explain the observed behaviour.

For each model, we also estimated the different contributions to the energy density and jet power, as reported in Table 8. The estimated magnetic and electron energy densities show a dominance of the latter by a factor of ∼10. This is in line with the results reported by Tavecchio & Ghisellini (2016), where one-zone SSC models find solutions that are relatively far from an equipartition state between the electron and magnetic energy densities, being in favour of the electron contribution by a factor of 10−100. The derived contributions to the jet power in forms of electrons (P_e), protons (P_p), Poynting flux (P_B), and radiation (P_rad) are also consistent with the range of values derived by Ghisellini et al. (2010).

We compared the derived disc luminosity and the power carried by the jet in the form of radiation, P_rad, with the relation obtained by Ghisellini et al. (2014) for these two quantities (see Fig. 1 from Ghisellini et al. 2014). The pair of values obtained for the two models presented here are consistent with the relation presented by these latter authors within errors and the 2−3σ confidence levels of this relation. The values of 3C 371 can be placed within the diagram from Ghisellini et al. (2014) in the same region as other BL Lac objects with similar faint accretion discs, again supporting the BL Lac nature of this source.

8. Summary and conclusions

Following the characterisation of the optical emission and variability presented in Paper I, we analysed the MWL emission of 3C 371. To this end, we used the exhaustive monitoring coordinated by the WEBT collaboration in optical and radio wavelengths, as well as UV, X-ray, and γ-ray data provided by the Swift-UVOT/XRT and Fermi-LAT telescopes. In addition, we studied the jet kinematics in the considered period through the analysis of VLBI MOJAVE radio data.

1. The γ-ray band shows variability of the same amplitude as that of the optical emission. Moreover, the UV and X-ray variability is slightly lower in amplitude than that seen in the optical, but with a much shorter temporal coverage coincident with a relatively faint emitting period. Finally, the radio band presents the least variability among all bands, as commonly observed in γ-ray-emitting blazars (see e.g. Ahnen et al. 2017, 2018).

2. We find strong correlations with a zero time lag between the optical and UV bands. The X-ray and γ-ray emission are much more weakly correlated, but with a large uncertainty on the time delays, which are still compatible with zero, and a lower statistical significance. Moreover, the radio emission shows an anti-correlated behaviour at zero time lag, which could however be an effect of the sampling and coverage of the data.

3. We studied the jet kinematics using the MOJAVE VLBI images during the period considered here. Several components are identified, with component B4 showing significant superluminal motion along the jet that allows us to constrain the kinematic parameters, resulting in an estimation of the apparent speed β_app, Doppler factor δ, Lorentz factor Γ, and viewing angle θ. The innermost components do not show significant velocities, in agreement with previous MOJAVE analyses (see Lister et al. 2019, 2021). The most external ones on the other hand show significant movement, although they are insufficiently bright for extraction of the kinematic parameters.

4. The optical–UV colour evolution shows clear mixed BWB and RWB trends on short timescales that can be associated with intrinsic jet processes (Raiteri et al. 2021a,b). On longer timescales, there appears to be a very mild BWB behaviour similar to that observed between the optical bands. However, the coverage of the UV observations provided by UVOT is much poorer in comparison to that in the optical R band.

5. The constructed optical SEDs from the quasi-simultaneous BVRI observations show a variability of approximately one order of magnitude in flux, and present a mild harder-when-brighter change of the spectral index after considering a power-law function to describe the spectrum. The X-ray SEDs and γ-ray spectral indices were found to be consistent with constant values of Γ = 1.63 ± 0.19 and α = 2.31 ± 0.03, respectively. Moreover, we calculated the average γ-ray spectrum over the entire period, which shows a ∼2.5σ preference for a log-parabola spectral shape over a power law.

6. We characterised the polarised emission at 15 GHz. Similarly to the optical polarisation, the radio polarisation degree is anti-correlated with the radio flux density. Moreover, the radio EVPA shows a rather erratic behaviour during the first studied epochs, stabilising in the latest ones with the same preferred direction as the optical EVPA.

7. The broadband SED of 3C 371 is successfully modelled within a leptonic scenario, where the bulk of the γ-ray emission is produced via SSC, with the EC having a negligible impact due to the faintness of the dusty torus and accretion disc photon fields. The variability of the high- and low-state SEDs is explained mainly through a hardening of the electron population distribution.

8. The derived parameters through the VLBI and broadband SED analyses and the energetics obtained for the relativistic jet are in agreement with 3C 371 being a BL Lac object observed under a high (∼10 deg) viewing angle.

This thorough characterisation of 3C 371 allows us to interpret its broadband variability and emission, and provides relevant information on the physical parameters describing the jet and confirms its misaligned BL Lac-type nature (Wrobel & Lind 1990), which has been debated in the past, with other authors proposing a radio galaxy classification (Miller 1975).

Data availability

The complete version of Table A.1 can be consulted in Zenodo: https://doi.org/10.5281/zenodo.14281786

Acknowledgments

We thank the anonymous referee for the thorough revision of the manuscript and the comments provided to improve the reported results. We thank Dr. Yuri Kovalev and Dr. Daniel Homan for the fruitful discussion about the VLBI data analysis and the preliminary independent MOJAVE kinematics results. J.O.-S., J.E.P. and I.A. acknowledge 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. J.O.-S. also acknowledges financial support from the projetct ref. AST22_00001_9 with founding from the European Union – NextGenerationEU, the Ministerio de Ciencia, Innovación y Universidades, Plan de Recuperaci’on, Transformación y Resiliencia, the Consejería de Universidad, Investigación e Innovación from the Junta de Andalucía and the Consejo Superior de Investigaciones Científicas; as well as from INFN Cap. U.1.01.01.01.009. C.M.R., M.V. and M.I.C. acknowledge financial support from the INAF Fundamental Research Funding Call 2023. 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. 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). 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). We acknowledge support by Bulgarian National Science Fund under grant DN18-10/2017 and Bulgarian National Roadmap for Research Infrastructure Project D01-326/04.12.2023 of the Ministry of Education and Science of the Republic of Bulgaria. 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-66/2024-03/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. The Medicina and Noto radio telescopes are funded by the Ministry of University and Research (MUR) and are operated as National Facility by the National Institute for Astrophysics (INAF). The Sardinia Radio Telescope is funded by the Ministry of University and Research (MUR), Italian Space Agency (ASI), and the Autonomous Region of Sardinia (RAS) and is operated as National Facility by the National Institute for Astrophysics (INAF). We thank the SRT operators Elise Egron, Delphine Perrodin and Mauro Pili for the help with the observations The Submillimeter Array is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics and is funded by the Smithsonian Institution and the Academia Sinica. We recognise that Maunakea is a culturally important site for the indigenous Hawaiian people; we are privileged to study the cosmos from its summit. This research has made use of data from the MOJAVE database that is maintained by the MOJAVE team (Lister et al. 2018).

References

Appendix A: Model fitting of the MOJAVE data

Here we include the modelfit results from the VLBI data analysis of all 18 epochs included in this work, from August 15, 2019 to February 5, 2021. Figure A.1 shows the VLBI images of all epochs, with the identification of the different components reported in the manuscript.

Fig. A.1. VLBI 2 cm images of 3C 371 for each observing epoch included in the analysis, showing the total flux density (black contours) and the polarised flux density (colour scale). Innermost and outermost contours correspond to flux density levels of 5.5 and 0.0015 Jy, with the rest of the contours logarithmically equispaced. The green ellipse at the top left corner (next to the first images) represents the beam size. The images highlight the positions of the core A0 and all identified components B1-L0. The date corresponding to each image is written on top of each plot.

Moreover, in Table A.1 below we compile the modelfit results obtained from our analysis of the VLBI data during all epochs, for all identified components. This table includes the information of the flux density and the position of each identified component for the different epochs included here. Detailed information on the products of the analysis can be found in the difmap documentation¹¹.

Appendix B: Radio polarisation data analysis

In this appendix we include the products of the polarisation estimation in the 15-GHz radio band based on the MOJAVE VLBI data for the 18 epochs analysed here. All the estimates were found to be consistent with those reported by the MOJAVE team in the dedicated webpage¹² to the observations of 3C 371. Following the typical VLBI radio data analysis approach, we assume a 5% systematic uncertainty to account for amplitude calibration errors when estimating the total flux and the Stokes parameters (Jorstad et al. 2017). The latter is translated into a 7% error of the polarised flux. For the error of the polarisation angle we assume a conservative value of Δθ_sys = 5 deg of systematic uncertainty, consistent with that obtained assuming the 5%-10% typically used for VLBI radio images (see for instance Jorstad et al. 2017; Kim et al. 2019; Weaver et al. 2022)¹³. This systematic error is later added in quadrature to obtain the total EVPA uncertainty as $\mathrm{\Delta }\theta ={(\mathrm{\Delta }{\theta }_{\mathit{sys}}^2+{(\mathrm{\Delta }P/P)}^2)}^{1/2}$.

Appendix C: Uncertainty evaluation of the SED parameters

In order to evaluate the statistical uncertainties and confidence bands of the models, we have performed a Markov chain Monte Carlo (MCMC) approach implemented in JetSeT making use of the emcee package (Foreman-Mackey et al. 2013). The best fit obtained from the SED modelling for each parameter is used for the MCMC procedure as the prior value. Then the asymmetric errors and the ∼3σ (5-95%) confidence level contours of the fit of each parameter are calculated through the MCMC approach with 500 iterations using 50 burn in steps and a number of walkers equal to 10 times the number of free parameters in each case, for a total of 100 and 150 walkers for the high and low state SEDs, respectively. The corresponding high- and low-state models with the 3σ significance contours as derived from the MCMC method are displayed in Fig. C.1. Moreover, the uncertainties of the derived parameters from both models are reported in corner plots resulting from the MCMC simulations, shown in Figs. C.2 and C.3.

We note that the uncertainties reported in Figs. C.2 and C.3 are representative of the set of parameters used in these two specific models. In fact, as a caveat, we note that typically these SED models are based on a set of a larger number of parameters, typically rather unconstrained. This leads to a degeneracy on the solution, with possibly different combinations of parameters leading to a statistically adequate solution in terms of goodness of the fit. Therefore, these uncertainties are not representative of the errors of other parameter combinations. In addition, we note that not all statistically possible solutions may have a physically meaningful interpretation. Therefore, an evaluation of the best fit parameters after the modelling is needed to validate whether the combination of parameters is physically suitable for describing the fitted SEDs and the broadband emission of the source.

Fig. C.1. Broadband SED models of 3C 371. Left: High-state SED. Right: Low-state SED. The spectral points of each model are represented with the blue markers. Downward pointing arrows represent upper limits. The black solid line corresponds to the best fit and the grey contours represent the 3σ confidence region of the fit. The bottom panel of each SED shows the residuals of the model.

Fig. C.2. Corner plot summary of the parameter distributions estimated through the MCMC procedure for the high state SED model. The panels on the diagonal show the distribution of each parameter, indicating the best fit value and the derived uncertainties with the black solid and dashed lines, respectively. The best fit value of each parameter and its corresponding uncertainties are also indicated on top of each histogram. The panels below the diagonal show the 2-dimensional projections of the histograms for each pair of parameters.

Fig. C.3. Corner plot summary of the parameter distributions estimated through the MCMC procedure for the low state SED model. The panels on the diagonal show the distribution of each parameter, indicating the best fit value and the derived uncertainties with the black solid and dashed lines, respectively. The best fit value of each parameter and its corresponding uncertainties are also indicated on top of each histogram. The panels below the diagonal show the 2-dimensional projections of the histograms for each pair of parameters.

Appendix D: Comparison with a leptonic SSC model

In this section we include as a comparison the models obtained assuming a simpler leptonic SSC model, without including the contributions of the low-energy BLR and DT photons to the observed γ-ray emission via EC scattering. We also ignore here the thermal emission from the DT and AD contributing to the low-energy IR and optical emission, respectively. Therefore, we consider only the non-thermal emission of the jet. These models are shown in Fig. D.1.

Fig. D.1. Broadband SED models for the high (left) and low (right) brightness states of 3C 371 considering only γ-ray emission produced through SSC scattering. The top panels contain the broadband SEDs, with the contribution of each component and the best fit model being represented with different markers, lines and colours as indicated in the legend. The bottom panels represent the residuals of the fitted models.

Comparing these models with those presented in Fig. 15 it is clear that the EC scattering contribution to the observed γ-ray emission is negligible in both emission states. Moreover, focusing on the IR-optical emission, we can also observe that during the high emission state, the DT and AD have also no impact on the total emission, being completely subdominant and opaqued by the much brighter synchrotron radiation from the relativistic jet. This is also true for the DT contribution during the low emission state.

For the case of the AD in the latter state however, we see that there is a small difference with respect to the model including this component in the fit. The low-state model presented in Fig. 15 shows a slight excess in the UV bands associated with the AD. Despite still being subdominant and one order of magnitude below the synchrotron emission, according to our estimations based on the optical spectra of 3C 371 (see Sect. 7), during the low state the AD contributes with roughly 10% of the emission at UV frequencies, leading to this small excess. This becomes even clearer when comparing the residuals of both models. While for the model including the EC components, DT and AD, all residuals are well within an interval of ±2, when considering only the SSC emission the UV residuals are clearly higher, especially in the UVOT w1 and m2 bands, where the emission from the AD peaks. This serves as additional confirmation of the presence of a faint AD in 3C 371, which becomes relevant when the synchrotron emission is faintest, as previously discussed.

All Tables

All Figures

	Fig. 1. Fermi-LAT γ-ray light curve of 3C 371. Top: 15 day binned γ-ray flux light curve between 100 MeV and 300 GeV obtained from the Fermi-LAT data. Middle: γ-ray spectral index as a function of time. Bottom: γ-ray spectral index as a function of γ-ray flux.
In the text

	Fig. 2. Swift-XRT X-ray light curve of 3C 371. Top: X-ray flux densities at 1 keV from Swift-XRT obtained with both the Cash (black squares) and the chi-squared (blue triangles) statistics. Middle: Power-law index Γ as a function of time. Bottom: Power-law index Γ as a function of 1 keV flux density.
In the text

	Fig. 3. Swift-UVOT optical–UV light curves (observed magnitudes) in the period considered in this paper. Red plus symbols correspond to photometry on single exposures; blue dots correspond to the results of co-adding the exposures of the same observations.
In the text

	Fig. 4. Radio light curves of 3C 371 compared with the dereddened and host-subtracted optical R-band flux density (top panel). Different colours and symbols are used to distinguish the contributing datasets, as specified in Table 2.
In the text

	Fig. 5. Amplitude variability for the different bands as a function of frequency. Different markers and colours represent different bands of the spectrum, as indicated in the legend.
In the text

Fig. 6. ZDCF of the radio, optical, UV, X-ray, and γ-ray light curves with respect to the optical R-band observations. From top to bottom: R − γ rays (left), R − 1 keV (right), R − w2 (left), R − m2 (right), R − w1 (left), R − B (right), R − V (left), R − I (right), R − 24 GHz (left), R − 8.5 GHz (right), and R − 4.8 GHz (centre). Different dotted lines denote the contours corresponding to a statistical significance of 1σ, 2σ, and 3σ.

In the text

	Fig. 7. ZDCF between the 8.5 GHz radio and the γ-ray-band light curves. Different dotted lines denote the contours corresponding to a statistical significance of 1σ, 2σ, and 3σ.
In the text

Fig. 8. Epoch 2020 Jan. 4 VLBI 2 cm image of 3C 371, showing the total flux density (black contours) and the polarised flux density (colour scale). Black line segments overlaid in the image represent the EVPA orientation. The innermost and outermost contours correspond to flux density levels of 5.5 and 0.0015 Jy, with the rest of the contours logarithmically equispaced. The green ellipse represents the beam size. The image highlights the core A0 and all identified components B1–L0.

In the text

	Fig. 9. Results of the kinematic analysis for component B4. Top: Distance of B4 from the jet core A0 over time. Bottom: Flux density of B4 over time.
In the text

	Fig. 10. Comparison of radio flux at 15 GHz and 8.5 GHz. Top: Total flux at 15 GHz derived from the MOJAVE VLBI images. Bottom: 8.5 GHz radio flux measured by the Svetloe and Medicina radio telescopes.
In the text

	Fig. 11. Dereddened, host-galaxy-corrected colour indices of 3C 371 for the different UV bands with respect to the optical R-band flux density. The colour scales represent the MJD of the observations. The black dash-dotted lines correspond to the linear fits describing the long-term colour changes. Left: Fw1/FR colour index. Middle: Fm2/FR colour index. Right: Fw2/FR colour index.
In the text

	Fig. 12. Optical SEDs constructed with quasi-simultaneous BVRI optical data. The colour bar indicates the MJD of the corresponding SED.
In the text

	Fig. 13. Optical R-band flux with respect to the estimated spectral index using a power-law model. The colour bar indicates the MJD of the corresponding SED.
In the text

	Fig. 14. Radio and optical polarisation of 3C 371. From top to bottom: 15 GHz flux density, polarisation degree at 15 GHz, EVPA at 15 GHz, intrinsic optical R-band polarisation degree, and optical R-band EVPA.
In the text

Fig. 15. Broadband SED models for the high- (left) and low-brightness states (right) of 3C 371. The top panels contain the broadband SEDs, with the contribution of each component and the best-fit model being represented with different markers, lines, and colours, as indicated in the legend. The bottom panels represent the residuals of the fitted models. We note that the EC scattering contribution via low-energy photons from the BLR was also considered. However, it is even fainter than that from the DT, and therefore does not fall within the scale of the figures and is not included.

In the text

Fig. A.1. VLBI 2 cm images of 3C 371 for each observing epoch included in the analysis, showing the total flux density (black contours) and the polarised flux density (colour scale). Innermost and outermost contours correspond to flux density levels of 5.5 and 0.0015 Jy, with the rest of the contours logarithmically equispaced. The green ellipse at the top left corner (next to the first images) represents the beam size. The images highlight the positions of the core A0 and all identified components B1-L0. The date corresponding to each image is written on top of each plot.

In the text

	Fig. C.1. Broadband SED models of 3C 371. Left: High-state SED. Right: Low-state SED. The spectral points of each model are represented with the blue markers. Downward pointing arrows represent upper limits. The black solid line corresponds to the best fit and the grey contours represent the 3σ confidence region of the fit. The bottom panel of each SED shows the residuals of the model.
In the text

Fig. C.2. Corner plot summary of the parameter distributions estimated through the MCMC procedure for the high state SED model. The panels on the diagonal show the distribution of each parameter, indicating the best fit value and the derived uncertainties with the black solid and dashed lines, respectively. The best fit value of each parameter and its corresponding uncertainties are also indicated on top of each histogram. The panels below the diagonal show the 2-dimensional projections of the histograms for each pair of parameters.

In the text

Fig. C.3. Corner plot summary of the parameter distributions estimated through the MCMC procedure for the low state SED model. The panels on the diagonal show the distribution of each parameter, indicating the best fit value and the derived uncertainties with the black solid and dashed lines, respectively. The best fit value of each parameter and its corresponding uncertainties are also indicated on top of each histogram. The panels below the diagonal show the 2-dimensional projections of the histograms for each pair of parameters.

In the text

Fig. D.1. Broadband SED models for the high (left) and low (right) brightness states of 3C 371 considering only γ-ray emission produced through SSC scattering. The top panels contain the broadband SEDs, with the contribution of each component and the best fit model being represented with different markers, lines and colours as indicated in the legend. The bottom panels represent the residuals of the fitted models.

In the text