© ESO 2020

1. Introduction

Blazars are a subclass of active galactic nuclei, with jet axes oriented close to the line of sight of the observer. They are divided into two subclasses, flat-spectrum radio quasars (FSRQs) and BL Lac objects (BL Lacs), which are thought to be intrinsically different. The former show broad emission lines in their optical spectra while BL Lacs have featureless spectra with weak or no emission lines (Stocke et al. 1991; Stickel et al. 1991). The spectral energy distribution (SED) of blazars exhibits a generic two-bump structure: one peak with a maximum in the spectral range from radio to X-rays and a second peak in the interval from X-rays to gamma rays. The radiation is produced in a highly beamed plasma jet and the double-peaked SED is often explained by a single population of relativistic electrons. The low-energy SED bump is thought to arise from synchrotron emission of particles within the magnetic field of the jet. The origin of the high-energy SED bump is less certain. It is commonly attributed to inverse Compton (IC) scattering of low-energy photons (Rees 1967). The low-energy photons can originate externally to the jet (external Compton scattering, Dermer & Schlickeiser 1993) or be produced within the jet via synchrotron radiation (synchrotron self-Compton scattering, SSC, Konigl 1981; Maraschi et al. 1992). As there is no observational evidence for strong external photon fields present in BL Lacs, the main population of seed photons for Compton scattering should originate from the synchrotron emission. In favour of this hypothesis comes from the fact that most of the SEDs of BL Lacs are well described with a simple one-zone SSC model (Bloom & Marscher 1996; Tavecchio et al. 1998; van den Berg et al. 2019). An alternative framework to explain the high-energy emission is the acceleration of hadrons, along with leptons (Mannheim & Biermann 1989). In the following, we focus on leptonic models.

Blazars are classified according to the frequency of the first peak of their SED into low- (LSP, ν_syn <  10¹⁴ Hz), intermediate- (ISP, 10¹⁴ ≤ ν_syn <  10¹⁵ Hz), and high- (HSP, ν_syn ≥ 10¹⁵ Hz) synchrotron-peaked objects (Abdo et al. 2010). Within the very high-energy (VHE; > 100 GeV) gamma-ray-emitting extragalactic objects, the most numerous sources are the HSP BL Lacs. With a large number of multi-wavelength (MWL) campaigns performed since the launch of the Fermi Gamma-Ray Space Telescope (Fermi), there is growing evidence that a one-zone SSC model, typically with a single spherical blob dominating the emission from optical to VHE gamma rays, is too simple to describe the SEDs of these objects (e.g. Ahnen et al. 2017a). Two component models, such as the spine-layer model by Ghisellini et al. (2005), have gained popularity. However, two-component models require a larger number of free parameters (twice as many as in single-zone ones) and therefore they often end up with a large degeneracy for the parameters involved (see e.g. Barres de Almeida et al. 2014). Also, the nature and the location within the jet of these two components is still unclear.

One way to constrain the two-component model is to derive the contribution of the different components from long-term variability. Aleksić et al. (2014) found a common increasing trend in radio and optical light curves of PKS 1424+240 and used this to constrain the contribution of the two components to the optical part of the SED. Lindfors et al. (2016) found a similar increasing/decreasing trend in radio and optical light curves of an additional 12 sources when analysing radio and optical light curves of 32 northern-sky BL Lacs. The authors argued that as the radio variability very closely traces the variability of the core flux in very large baseline interferometry (VLBI) images, the slowly varying optical flux also originates from the core. The fast varying component of the optical flux could instead originate from a distinct, smaller emission region. In this work we selected a subsample of five of the BL Lacs from Lindfors et al. (2016) based on the availability of MWL data (VER J0521+211, PKS 1424+240, 1ES 1727+502, 1ES 1959+650, 1ES 2344+514) with the aim being to place observational constraints on two-component models. Independent from these common trends on long-term light curves, we also use optical polarisation data to disentangle the contribution of the two components and we take into account constraints on jet parameters from the VLBI observations.

The paper is organised as follows: the observations, analysis methods, and wavelength-specific results of our subsample are described in Sect. 2. The observational constraints for SED modelling from VLBI data, MWL light curves, and optical polarisation observations are derived in Sect. 3. In Sect. 4, the SED modelling of all five sources is described. Section 5 includes the discussions of the results of the SED modelling. Finally, in Sect. 6 we present the summary and conclusions of the main results of the paper.

2. Observations, data analysis, and results

The general properties of our sample are listed in Table 1. Power law (PL) and log-parabola (LP) are the two mathematical functions which are employed for our spectral analysis in different bands. They are defined as follows: A simple power law

$$\begin{array}{c}\hfill \frac{\mathrm{d}F}{\mathrm{d}E}(E)=F_0(\frac{E}{E_0}{)}^{-\mathrm{\Gamma }},\end{array}$$(1)

and a log-parabola

$$\begin{array}{c}\hfill \frac{\mathrm{d}F}{\mathrm{d}E}(E)=F_0(\frac{E}{E_0}{)}^{-\mathrm{\Gamma }-\beta ({log}_{10}(E/E_0))},\end{array}$$(2)

where dF/dE is the differential flux as a function of the energy E. Here, F₀, Γ, and β are the flux at the normalisation energy E₀, the spectral index, and the curvature parameter of the spectrum at E₀, respectively.

2.1. Very high-energy gamma rays (MAGIC)

The Major Atmospheric Gamma-ray Imaging Cherenkov experiment (MAGIC, Aleksić et al. 2016) is a system of two 17m diameter telescopes located at the Observatorio del Roque de los Muchachos (ORM), La Palma, Canary islands, Spain. The objects of our sample were observed by MAGIC between 2013 and 2016 as part of different observation campaigns (see Table 2 for a detailed list of included epochs for each source). The data were analysed using the MAGIC Standard Analysis Software (MARS, Moralejo et al. 2009; Zanin et al. 2013) taking into account the instrument performance under different observation conditions (Aleksić et al. 2016; Ahnen et al. 2017b).

We calculated the VHE gamma-ray integral flux of each object and searched for variability at different timescales (from 10 min to 1 week). The constant-flux hypothesis on a one-day timescale is rejected at the 3σ confidence level for 1ES 1727+502, 1ES 1959+650, and 1ES 2344+514 (see below). For PKS 1424+240, no variability was found during the 2014 (MJD 56740-56826) and 2015 (MJD 57045-57187) campaigns. However, the VHE gamma-ray flux of the 2015 campaign was ∼60% of that observed during the 2014 campaign. Therefore, the data are divided into the 2014 and 2015 campaigns. In the case of VER J0521+211, we do not find any significant variability over the four nights of the MAGIC observation in 2013 (Prokoph et al. 2015).

For 1ES 1727+502, there is one night (MJD 57309, 2015 October 14) when the VHE gamma-ray flux was 52% of the average flux. The VHE gamma-ray flux during MJD 57309 was 3.3σ away from the average flux computed using all of the five nights of observation. However, the VHE gamma-ray spectrum could not be computed using the observations on this single night. Therefore, we reproduced the VHE gamma-ray spectrum of this source using all available observations. Exclusion/inclusion of the observation on MJD 57309 did not affect the parameters describing the VHE gamma-ray spectrum.

1ES 1959+650 was in a flaring state during 2016. We selected three different nights based on the level of VHE gamma-ray flux of the source during 2016 and availability of the simultaneous MWL observations at lower energy bands. The highest, intermediate, and lowest VHE gamma-ray flux was observed on 2016 June 14 (MJD 57553), June 8 (MJD 57547), and November 20 (MJD 57711), respectively. No intra-night variability was detected in the data from these selected observations. This is in agreement with the results reported by MAGIC Collaboration (2020a), where a detailed variability analysis was performed on three nights of the highest detected fluxes (including MJD 57553) during the 2016 campaign and intra-night variability (with a timescale of 35 minutes) was found on the nights of 2016 June 13 and 2016 July 1 (MJD 57552 and 57570; see Table 3 in MAGIC Collaboration 2020a).

1ES 2344+514 showed variability on a daily timescale but no shorter variability timescale was detected in the VHE gamma-ray band. MAGIC Collaboration (2020b) performed a detailed analysis on different emission states of this source and found the spectral shape to be similar during different observational epochs despite different levels of VHE gamma-ray flux. The results of the VHE gamma-ray flux study of our sample are summarised in Table 2. The derived variability timescales are further used in Sect. 4.

The VHE gamma-ray spectra are computed for each source and epoch separately in case the source shows variability. The effect of the extragalactic background light (EBL) on VHE gamma-ray spectra was taken into account using the model of Domínguez et al. (2011). Two different models (PL and LP) were then tested. The LP model was preferred over the PL model at 3σ confidence level if the F-test probability value was less than 0.27%. The results of the spectral analysis in the VHE gamma-ray band are summarised in Table 3.

2.2. High-energy gamma rays (Fermi-LAT)

The Large Area Telescope (LAT), based on the pair conversion technique, is the high-energy instrument on-board the Fermi. It has continuously monitored the high-energy (HE, 100 MeV < E <  300 GeV) gamma-ray sky (Atwood et al. 2009) since its launch in 2008. The six-year MJD 56200 (2012 September 4) to 58340 (2018 August 9) light curve for each source was obtained by applying a weekly binning to the events collected by the LAT with an energy higher than 100 MeV over a region of interest of 10° centred on the selected sources. Time intervals coinciding with bright solar flares and gamma-ray burst were excised from the data set as is done in the fourth Fermi-LAT source catalogue (4FGL, Fermi-LAT Collaboration 2020). The data reduction and analysis of the events belonging to the Pass8 source class was performed with the FermiTools software package version 11-07-00 and fermipy (Wood et al. 2017) version 0.17.4. To reduce the Earth limb contamination, a zenith angle cut of 90° was applied to the data. To calculate the weekly flux of the selected sources, a likelihood fit to the data was performed including each source of interest, modelled with a power-law spectrum, the Galactic diffuse-emission model¹ (gll_iem_v07.fits), and isotropic component (iso_P8R3_SOURCE_V2_v1.txt) recommended for the Pass8 Source event class as well as the sources of the Fermi-LAT 4FGL within 15° of the position of the source of interest. The normalisation of both diffuse components in the source model were allowed to vary during the spectral fitting procedure. The normalisation was allowed to vary for the sources located at a distance of less than 2° from the source of interest and with a detection test statistic (TS²) higher than 50 integrated over the full data set. The sources located at distances between 2° and 7° had their normalisation set as a free parameter if their variability index was higher than 18.48³. The spectral index of all the sources with free normalisation was left as a free parameter if the source showed a TS value higher than 25 over an integration time of one week; in all the other cases the indexes were frozen to the value obtained in the overall fit⁴. We apply the correction for the energy dispersion to all sources except for the isotropic background. The HE light curves are shown in Figs. 1–5. The spectrum was obtained from analysing only the data collected over the selected epochs, which were (quasi-)simultaneous to MAGIC data and had sufficient statistics to compute at least two spectral data points per decade in the energy range between 100 MeV and 300 GeV (Table D.1). In all cases, the LP model can describe the spectra of the sources better than the PL model at 4σ confidence level, except for 1ES 1727+502 where the LP model was not statistically preferred over the PL model. These findings are inline with the results reported in the 4FGL catalogue. Moreover, except for PKS 1424+240, which showed a harder spectrum during the 2015 campaign, the spectral parameters were in agreement with those reported in the 4FGL at 3σ confidence level.

Fig. 1. MWL light curves of VER J0521+211 in the range from MJD 56200 (2012 September 30) to 58400 (2018 October 9). From top to bottom panels: radio and VLBI core flux (15 GHz), optical (R-band), optical polarisation degree, electric vector polarisation angle, X-ray flux (2–10 keV), HE gamma-ray photon flux (0.1–300 GeV), and VHE gamma-ray photon flux above the threshold energy given in the panel. Black arrows show the 95% confidence level upper limits. The data, which are marked with vertical lines/area and squares in different bands, are used in the SED modelling.

Fig. 2. Same description as in Fig. 1 for PKS 1424+240.

Fig. 3. Same description as in Fig. 1 for 1ES 1727+502.

Fig. 4. Same description as in Fig. 1 for 1ES 1959+650.

Fig. 5. Same description as in Fig. 1 for 1ES 2344+514.

2.3. X-ray and UV (Swift)

The X-ray Telescope (XRT, Burrows et al. 2004) on-board the Neil Gehrels Swift observatory (Swift) has been observing the sources in the above-mentioned sample since 2004 in both photon-counting (PC) and window-timing (WT) modes. The multi-epoch event lists for the period from 2012 September 30 to 2018 October 9 were downloaded from the publicly available SWIFTXRLOG (Swift-XRT Instrument Log)⁵. Following the standard Swift-XRT analysis procedure described by Evans et al. (2009), the PC observation data were processed using the configuration described by Fallah Ramazani et al. (2017) for blazars. For the data from WT observations, we defined the source region as a box with a length of 40 pixels centred on the source position and aligned to the telescope roll angle. The background region is defined by a box with a length of 40 pixels aligned to the telescope roll angle and 100 pixel away from the centre of the source. For both modes of observation, due to the open issues for analysing the Swift-XRT data⁶, we fitted the spectra of each observation in the 0.3–10 keV energy range assuming all possible combinations of pixel-clipping and point-spread-function together with two mathematical models (i.e. PL and LP), a normalisation energy E₀ = 0.3 keV, and the fixed equivalent Galactic hydrogen column density reported by Kalberla et al. (2005) and listed in Table 1. In total, for each XRT observation, 6 and 16 spectra (for PC and WT modes, respectively) were extracted and the best-fitted model was selected using least χ² and F-test methods. The results of this analysis are partially (only X-ray flux in range of 2–10 keV) presented in Figs. 1–5 for each source. An example of full version of the results is presented in Table A.2, while the complete version of the results is available at the CDS. All sources are variable in the X-ray band in the studied time period.

During the Swift pointings, the UVOT instrument observed the sources in our sample in its optical (V, B, and U) and UV (W1, M2, and W2) photometric bands (Poole et al. 2008; Breeveld et al. 2010). We selected the UVOT data (quasi-)simultaneous to the MAGIC observations and analysed the data using the uvotsource task included in the HEAsoft package (v6.25)⁷. Source counts were extracted from a circular region of 5″ radius centred on the source, while background counts were derived from a circular region of 20″ radius in a nearby source-free region. The contribution of the host-galaxy flux in the UVOT bands is derived from the R-band values reported by Nilsson et al. (2007) and the conversion factors reported by Fukugita et al. (1995). The host-galaxy subtracted (when applicable) UVOT flux densities, corrected for extinction using the E(B − −V) values from Schlafly & Finkbeiner (2011) and the extinction laws from Cardelli et al. (1989), are used in the SED modelling (Sect. 4).

2.4. Optical and radio (Tuorla, OVRO, and MOJAVE)

The optical (R-band) data from MJD 56200 (2012 September 30) to 58320 (2018 July 21) obtained in the framework of the Tuorla blazar monitoring programme⁸ using the 50 cm Searchlight Observatory Network telescope (San Pedro de Atacama, Chile), the 40 cm Searchlight Observatory Network telescope (New Mexico, USA), and the 60 cm telescope at Belogradchik (Bulgaria) in addition to the Kungliga Vetenskapsakademien (KVA) telescope (ORM, La Palma, Canary islands, Spain). Most of the observations were performed with the KVA telescope. The data were analysed and calibrated using the method described by Nilsson et al. (2018), and were corrected for Galactic extinction and host galaxy emission for aperture photometry. The correction coefficients and the aperture used for each individual source are summarised in Table 1. The results of these observations are presented in Figs. 1–5. An example of the results is presented in Table A.1, while the complete version of the results is available at the CDS.

We used the long-term light curves from the Owens Valley Radio Telescope (OVRO, 15 GHz). This programme, the observations, and the analysis methods are described in Richards et al. (2011). In this work we have included the data from the time interval MJD 56200-58320 (2012 September 30–2018 July 21). We note that the light curves cover data from the period between 2015 August 1 and November 24 when the instrument was not working optimally. Therefore, the noise in the data is higher during this period. Moreover, we collected the core flux at 15 GHz using the data from the Monitoring Of Jets in Active galactic nuclei with VLBA Experiments (MOJAVE) programme (Lister et al. 2019).

2.5. Optical polarisation (NOT)

The sources in our sample were monitored with the Nordic Optical Telescope (NOT). The ALFOSC⁹ instrument was used in the standard setup for linear-polarisation observations (λ/2 retarder followed by a calcite plate). The observations were performed two to four times per month in the R-band between 2014 and 2018. The data were analysed as described by Hovatta et al. (2016) and MAGIC Collaboration (2018b) with a semi-automatic pipeline using standard aperture photometry procedures. The data were corrected for Galactic extinction and host-galaxy emission using the values listed in Table 1. The results of these observations are presented in Figs. 1–5.

3. Observational constraints for two-component emission modelling

Two-component models are models where two emission regions are responsible for the observed radiation. There is some evidence showing that there is a correlation between the X-ray and VHE gamma-ray bands in HSP BL Lacs (Acciari et al. 2011; Aleksić et al. 2015). This suggests that the observable emission in these two wavebands originates from the same region. The second component is the one we see dominating in the radio band. In the optical band, we presumably see a mixture of these two components. In this section, we use the radio-to-X-ray data to obtain constraints for these two emission regions to be used in the SED modelling.

As discussed in the introduction, Lindfors et al. (2016) found common trends in the long-term optical and radio variability for all five sources of our sample. They also showed that the brightness of the core in the 15 GHz VLBI images (hereafter VLBI core) closely follows the 15 GHz light curve, as had been previously found at higher frequencies (37 and 43 GHz Savolainen et al. 2002), and therefore suggested that the common slowly varying radio-optical emission region is located at the VLBI core. We therefore collect the observational constraints on the jet parameters from VLBI observations to be used directly in the SED modelling (Sect. 4.1).

On top of the slow variability, the optical band also shows a fast variability, which could originate from a second emission component. For simplicity, we assume that this component is the one dominating the X-ray and VHE gamma-ray emission. In order to constrain the contribution of these two emission regions to the optical flux, we use the long-term light curves and optical polarisation data described in Sect. 2 by implementing two independent procedures described in Sects. 3.2 and 3.3. We also searched for the correlations between different long-term light curves to determine if the same region produces the emission observed at different wavelengths.

3.1. Constraints on jet parameters from VLBI

The arguably strongest observational evidence for two-component models comes from VLBI observations. In polarimetric VLBA observations, Attridge et al. (1999) discovered a clear difference in the polarisation direction of the inner jet and outer layer of FSRQ 1055+018 and similar polarisation structures have been observed in several sources since then (Pushkarev et al. 2005; Gabuzda et al. 2014). Another indication of a spine–sheath structure of the jets is the so-called limb brightening, where the edges of the jet appear brighter than the central spine which has been observed in several radio galaxies and blazars (Giroletti et al. 2004; Nagai et al. 2014; Gabuzda et al. 2014). In particular, such limb brightening has been observed in Mrk 501 (Piner et al. 2009) which is a rather similar source to those in our sample (in terms of VLBI jet properties and synchrotron peak frequency).

The sources in our sample are rather weak in the radio band and therefore potential spine–sheath structures would be impossible to resolve. Their VLBI images all show compact jets in which the core is the brightest component. The core fluxes follow the total intensity variations observed at 15 GHz (see Figs. 1–3, the two other sources had no or only one simultaneous core flux measurement), which agrees with the results found in larger samples suggesting that the radio emitting component is located at the 15 GHz core. We used VLBI data to constrain some of the jet parameters: the apparent speed of the jet, the size of the VLBI core, the jet position angle, and the core polarisation properties. The jet velocities and the size of the core are used directly in the SED modelling. The jet position angle and polarisation properties are only used for comparison with the optical polarisation data in Sect. 3.3. These were collected from Hodge et al. (2018), the authors of which report an uncertainty in the fractional polarisation of approximately 7% of the given values and an electric vector polarisation angle (EVPA) accurate to within 5°, while no error is given for the jet position angle.

The major fraction of the jet parameter constraints are from the MOJAVE programme (Lister & Homan 2005; Lister et al. 2009, 2016, 2019; Hodge et al. 2018), during which observations were performed at 15 GHz. All of the sources in our sample were observed as part of this programme. However, not all of the collected data were obtained between 2013 and 2018, and most of the sources were observed only a few times. Additionally, we collected the results reported by Piner & Edwards (2004, 2018), Piner et al. (2008, 2010), and Tiet et al. (2012).

VER J0521+211 has been observed seven times in the framework of the MOJAVE programme between 2014 and 2018. Lister et al. (2019) reported several moving components in the jet. The fastest component has a maximum jet speed of μ = 774 ± 45 μas yr⁻¹ which corresponds to an apparent projected speed of β_app = (7.72 ± 0.42)c considering z = 0.18. Assuming viewing angles of 3° and 5° gives Doppler factors of δ ∼ 15 and δ ∼ 11, respectively. The median fractional polarisation and the EVPA of the core correspond to 0.5% and ∼200° respectively, and any significant variability is ruled out within the observations. The jet position angle is 250° (Hodge et al. 2018).

PKS 1424+240 was observed three times between 2013 and 2016 in the framework of the MOJAVE programme. Significant motion was detected for two components. The fastest component has a maximum apparent speed of β_app = (2.83 ± 0.89)c. The latter corresponds to Doppler factors of 10 or 7 for viewing angles of 3° and 5°, respectively. The core is polarised with a median fractional polarisation of 2.3% and the EVPA lies in the range 140 − 150°. The jet position angle is 140° (Hodge et al. 2018).

1ES 1727+502 was observed five times between 2013 and 2015. Piner & Edwards (2018) fitted four components to the MOJAVE data of 1ES 1727+502 and they were all consistent with no motion (Lister et al. 2019 reports 0.041 ± 0.043 as maximum jet speed). The polarisation of the core is not significant, and therefore the EVPA cannot be derived. The position angle of the jet is 270° (Hodge et al. 2018).

1ES 1959+650 was dropped from the MOJAVE programme in 2009 as it is too compact and weak at 15 GHz and there is no data from the source between 2013 and 2018. In the earlier data, the source showed a polarisation degree of 2.6–4.5% and its polarisation angle was relatively stable at 144 − 160°. The position angle of the jet is 140 − 160°. The apparent speeds are in agreement with no motion (Piner & Edwards 2004; Piner et al. 2010) which is also confirmed by Hodge et al. (2018) and Lister et al. (2019).

1ES 2344+514 was observed three times between 2013 and 2018. Lister et al. (2016) reported one component with β_app = (0.055 ± 0.036)c and the most recent measurements are in line with this result (β_app = (0.037 ± 0.012)c, Lister et al. 2019). The polarisation of the core is not significant (Hodge et al. 2018), while the jet position angle is 130 − 145° (Piner & Edwards 2004).

Finally, we collected the measured full width at half maximum (FWHM) values of the major axis core region to estimate the size of the VLBI core. For each source we selected a MOJAVE epoch at which the core was resolved and if there were several, we selected the one closest to the epoch used for the SED modelling. As mentioned above, for 1ES 1959+650, the latest MOJAVE observation epoch was in 2009, which is 7 years before our SED data, but the values reported for 2000–2009 were all very similar (Lister et al. 2019). We used the measured FWHM values of the major axis values 0.08, 0.14, 0.09, 0.09 and 0.13 mas as the diameter of the core emission region for VER J0521+211, PKS 1424+240, 1ES 1727+502, 1ES 1959+650, and 1ES 2344+514, respectively.

3.2. Long-term light curves

Lindfors et al. (2016) analysed the long-term radio (15 GHz) and optical (R-band) light curves of the sources studied in this paper using the data from 2008–2013. We repeated the same analysis procedure using data collected between 2013 and 2018 to investigate whether the results obtained by Lindfors et al. (2016) are temporary or not and to use these results for the two-zone SED modelling, in particular to constrain the contribution of the two components in the optical band.

3.2.1. Common emission component at radio and optical frequencies

Following Lindfors et al. (2016), we analysed the long-term radio and optical light curves to separate the slowly varying component from the optical light curves and estimate its minimum contribution to the optical flux.

Briefly, the analysis steps are as follows:

– Test if there are linear correlations between time and flux density in radio and optical light curves. The Spearman r-values for optical and radio light curves are reported in Table 4.

– Fit a polynomial to the radio light curve (see Fig. 6, left panel). The polynomial is determined by adding new orders until the χ²/d.o.f. of the fit no longer improves. The resulting polynomial does not describe all of the radio variability, in particular short flares are not described by this polynomial.

Fig. 6. Steps of the analysis used to determine the contribution of the common emission component at radio and optical frequencies. Left: fitting a polynomial to the radio light curve. Middle: polynomial fit scaled to the average flux of the optical light curve and multiplied with a scaling factor (different colours correspond to different scaling factors, see text). Right: scaled polynomial is subtracted from the optical data (green filled circles with errors). The RMS of the resulting light curve (purple filled circles) is compared with the RMS of the original data. These analysis steps are shown for all sources from bottom to top: VER J0521+211, PKS 1424+240, 1ES 1727+502, 1ES 1959+650 and 1ES 2344+514. In the case of PKS 1424+240 subtracting the polynomial did not decrease the RMS of the optical light curve and therefore the purple dots are under the green symbols and are not well visible in the rightmost panel.

– The polynomial fit is scaled to the average flux of the optical light curve (see Fig. 6, middle panel), and then it is multiplied with a scaling factor assuming values from 0.1, 0.2, 0.3, ...to 1 and the resulting curve is subtracted from the optical data. The root mean square (RMS)¹⁰ of the resulting light curves are calculated and the one that minimises the RMS is selected as the best fit.

– Estimate the fractional contribution of this slowly varying component to the optical flux density by dividing the RMS of the best-fit-subtracted data (RMS1) with the RMS of the original data (RMS2; see Fig. 6, right panel). The contribution of the slowly varying component to the optical variability is consequently 1-RMS1/RMS2 and is given in Table 4. As discussed in Lindfors et al. (2016), this fraction represents the minimum contribution, because in addition to this slow variation, there can be flares in the radio band that are not reproduced by this polynomial. The minimum fraction is then used to guide the two-component modelling (see the following section).

We find that, in all of the five sources, the increasing or decreasing trends in the radio light curves have persisted and are in the same direction as found in Lindfors et al. (2016). This is interesting because the time-spans of the light curves used in these works are different. This means that the increasing or decreasing trends in the radio light curves have persisted for a timescale of ∼10 years.

In the optical, significant trends have persisted in four sources out of five, the exception being PKS 1424+240, where there is no significant rising or decaying trend. Accordingly, the minimum fraction of a slowly varying component (common with radio) contributing to the optical flux is zero for this source. For other sources, the fraction varies from 6% to 52% (see Table 4).

3.2.2. Correlation studies

We calculated the cross-correlation function between three pairs of light curves (radio – optical, radio – X-rays, and optical – X-rays) for each source following the method described by Max-Moerbeck et al. (2014a) and Lindfors et al. (2016). We did not include the HE gamma-ray light curves as the uncertainties of the data points are rather large. Moreover, Liodakis et al. (2018) performed a cross-correlation analysis between the radio/optical and HE gamma-ray bands on a sample of 145 blazars which includes four of the sources in our sample (except 1ES 1727+502). They found only one significant correlation (> 3σ confidence level) between optical and HE gamma-ray bands for VER J0521+211. The VHE gamma-ray light curves were too sparse to be included in the correlation study.

We used the discrete correlation function (DCF; Edelson & Krolik 1988) with local normalisation (LCCF; Welsh 1999). We used temporal binning of 10 days and require that each LCCF bin have at least ten elements. Following Max-Moerbeck et al. (2014b), the significance of the correlation was estimated using simulated light curves. In the simulations, we used power spectral density (PSD) indices of −2.35, −1.55, and −1.7 for the radio light curves of PKS 1424+240, 1ES 1959+650, and the three other sources in our sample, respectively (determined from the radio data: Max-Moerbeck, priv. comm.; Lindfors et al. 2016). For the optical, we used PSD indices reported by Nilsson et al. (2018) for each source except for VER J0521+211 which was not included in their sample. We used the same method as described by Nilsson et al. (2018) and calculated the PSD index of the optical light curve of VER J0521+211 to be −1.6. For the X-ray light curves we used the PSD index value of −1.4 (Aleksić et al. 2015).

The results of the cross-correlation analysis are illustrated in Appendix B. While there are several peaks (or rather wide features) in the LCCFs which reach 2σ significance level, only the radio – optical data sets of two sources (1ES 1727+502 and 1ES 1959+650) show significant correlations (3σ level of confidence). The two significant radio – optical correlations are rather wide (90 and 60 days for 1ES 1727+502 and 1ES 1959+650, respectively) and compatible with zero-days time lag. The LCCF peaks are located at −10 and +40 days for 1ES 1727+502, while they are located at −30 and −20 days for 1ES 1959+650¹¹. Again, PKS 1424+240 is an exception, as there are no 2σ peaks in the radio-optical LCCF, which is in agreement with the results in Sect. 3.2.1 and different from results in Lindfors et al. (2016). In general we find our radio–optical results to be in good agreement with those reported by Lindfors et al. (2016) and Liodakis et al. (2018). Only one radio–X-ray correlation is found for the case of 1ES 1727+502 with the time lag of 680 ± 20 days where the radio flare is leading the X-ray outburst. However, the length of the X-ray light curve is rather short (1200 days) and this correlation could be the artefact of associating the X-ray outburst with one of the previous flares in radio when the X-ray data were not available. It is notable that the optical–X-ray correlation for 1ES 1727+502 shows many features. However, these features are the result of a single dominating outburst in X-rays which results in a time-delay peak with every optical flare, which is one of the limitations of the LCCF method (Emmanoulopoulos et al. 2013).

Correlations are generally used to probe whether or not the emission regions in different energy bands are causally connected. Our results support that at least in the case of 1ES 1727+502 and 1ES 1959+650 the radio and optical emission would partially originate from the same emission region (as the time lag is consistent with zero), which is in line with the result in Sect. 3.2.1.

3.3. Polarisation analysis

The observed optical polarisation of blazars usually contains signatures of two components: an optical polarisation core and a stochastic component (see e.g. Valtaoja et al. 1991; Villforth et al. 2010; Barres de Almeida et al. 2010). Barres de Almeida et al. (2014) made a first attempt to separate the two components and evaluate their relative strengths from the optical polarisation data. We follow this idea, but instead of the iterative fitting applied there, we used a physical model and Bayesian fitting methods.

To do this, we assumed that the R-band flux originates from two components, referred to as the “constant” and “variable” components in the following. We thus have for the Stokes parameters

$$\begin{array}{cc}& I=I_C+I_V\hfill \\ \hfill & Q=Q_C+Q_V\hfill \\ \hfill & U=U_C+U_V,\hfill \end{array}$$(3)

where the subscripts C and V refer to the constant and variable components, respectively. The observed degree of polarisation (PD) and EVPA are then

$$\begin{array}{c}\hfill \mathrm{PD}=\sqrt{{\left(\frac{Q}{I}\right)}^2+{\left(\frac{U}{I}\right)}^2},\end{array}$$(4)

and

$$\begin{array}{c}\hfill \mathrm{EVPA}=\frac{1}{2}{tan}^{-1}\left(\frac{U}{Q}\right),\end{array}$$(5)

where −π ≤ EVPA ≤ π. The constant component was modelled directly by letting I_C, Q_C, and U_C be free parameters, whereas the variable component had nine free parameters (see below and Appendix C). We modelled the variable component as a homogeneous cylindrical emission region in a jet with a helical magnetic field and computed the Stokes parameters using the formulae described by Lyutikov et al. (2005). We assumed that the orientation of the variable component remains constant with respect to the observer, which means that any change in the polarisation of the source must arise from the change of the relative flux ratio between the constant and variable component. This is because in the formulation by Lyutikov et al. (2005) the EVPA of the radiation is always either parallel or perpendicular to the direction of the relativistic outflow.

We describe the parameters of the model, the assumptions we made, and the details of the fitting procedure in Appendix C. In short, the model has 12 free parameters. Most of the parameters in the model cannot be constrained with monochromatic observational data due to a high degree of degeneracy. We fixed 5 of the parameters (of the variable component): index of the electron spectrum, p, to 2.1; the radius of the emitting region, r, to 2.5 × 10¹⁵ cm; the length of the emission region, l, to 5 × 10¹⁵ cm; the magnetic field strength, B₀, to 0.1 Gauss; and the apparent speed β to 0.99. These values are similar to those applied for the SED modelling in the following section (see Sect. 4.1). This model was fitted to the observed R-band polarisation data (Sect. 2.5) in the Q − U plane. One important ingredient of the model is σ, the standard deviation on random variations of Q and U. This parameter adjusts itself according to the predictive power of the other parameters. This parameter is discussed in more detail in Appendix C.

The results of the fitting procedure are reported in Table 5. The errors represent the 68% confidence intervals derived from marginalised distributions. One of the main goals of this fitting procedure was to obtain some constraints on the flux ratio of the two emission components in the optical band to be compared to the ratio derived in Sect. 3.2.1. Unfortunately, this was not achieved in all cases. For instance, in the case of VER J0521+211 the priori range for I_C was from 0 to 3.0 mJy (see Col. 3 in Table 5) and the posteriori averages in the middle of this range with errors that fill the priori completely. The flux ratio I/I_C is best constrained in the case of PKS 1424+240. Therefore, the polarisation study performed here provides limited additional constraints for the SED modelling in this work, but we intend to perform a more detailed study of this method in future work. The results on the flux ratio of the two components in the optical band are compared to those obtained with the decomposition of the long-term light curves in Sect. 5.5. For that purpose, we calculated I/I_C for the SED modelling epochs, i.e. I is the total optical flux in the periods reported in Appendix D. These values are reported in Sect. 5.5.

Finally, we compare the observed optical EVPAs and jet position angles that we derive with our fitting to those from VLBI observations (see Sect. 3.1). If the radio and optical emission originate from the same region, there should be agreement between the optical and VLBI results. BL Lacs objects have a preferred orientation of position angle, i.e. the EVPA is often stable. This feature can be interpreted as the stability of the emission region geometry in the optical band (Angel et al. 1978; Jannuzi et al. 1994; Jorstad et al. 2007) and is also seen in our optical polarisation data (see middle panel in Figs. 1–5). In two cases (1ES 1959+650 and 1ES 2344+514), our jet position angle agrees well with the VLBI angle (∼50% probability of being the same), indicating that the EVPA is parallel to the jet for these sources. The EVPA of the radio core in 1ES 1959+650 is similarly aligned. In the case of VER J0521+211, if we pick the solution with φ₀ = 13° and take into account the 180° ambiguity, good agreement is again achieved with the radio core EVPA = 20°. For PKS 1424+240, similarly good agreement is found if we assume the EVPA to be perpendicular to the jet. For 1ES 1727+502, the agreement is not so clear. Out of two solutions for φ₀, one is too noisy to draw conclusions and the other one cannot be made compatible with the jet radio position angle. As a general conclusion, there appears to be a close correlation between the radio and optical results, which is in line with the results from previous comparisons on TeV BL Lacs Hovatta et al. (2016). These latter authors found a difference between the EVPA and the jet position angle of less than 20° (i.e. the magnetic field is perpendicular to the jet direction) for two-thirds of the sources within a sample of 9 TeV BL Lacs. Given that our errors are approximately 10° and that we can choose from four different angles in the range from 0 to 360°, it is not clear if this agreement is statistically significant in our case, but it is certainly in line with a common origin of radio and optical emission in these sources.

4. SED modelling

4.1. Two-component model

The SEDs are modelled with a two-component model based on Tavecchio et al. (2011) which calculates synchrotron and SSC emission for spherical emission regions and takes also into account synchrotron-self absorption. It is similar to the one used in Aleksić et al. (2014), but the two emission regions are considered to be co-spatial and interacting as in MAGIC Collaboration (2018a); MAGIC Collaboration (2019) to mimic a simple spine–sheath model (see Sect. 3.1). We call the two emission regions “core” and “blob”, with sizes R_core >  R_blob (see Fig. 7). These two regions correspond to the constant and variable components defined in Sect. 3.3, respectively.

Fig. 7. Sketch of the geometrical modelling setup. The two emission regions are located several parsecs from the central black hole (at dcore). The smaller emission region (blob) is embedded in the larger emission region (core) and the interaction of the two emission regions provides additional seed photons for the Compton scattering; see discussion in Tavecchio et al. (2011, Appendix B).

The regions are filled with electrons distributed in Lorentz factor according to a smoothed broken power law (in the following, physical quantities are expressed in the co-moving frame of each individual region):

$$\begin{array}{c}\hfill N(\gamma )=K{\gamma }^{-n_1}{(1+\frac{\gamma }{{\gamma }_b})}^{n_1-n_2},{\gamma }_{\min }<\gamma <{\gamma }_{\max }.\end{array}$$(6)

The distribution has a normalisation K between γ_min and γ_max and slopes n₁ and n₂ below and above the break in the electron distribution, γ_b (Maraschi & Tavecchio 2003). Each of the emission regions has size R, Doppler factor δ, and magnetic field strength B, for which we searched for constraints from observations:

– The sizes of the core emission region were derived from VLBI observations (see Sect. 3.1). The sizes are of the order of 10¹⁷ cm. We note that the derived sizes would suggest variability timescales shorter than what we obtain for the slowly varying component from the data. This means that the origin of the slow variability cannot be the delay caused by a core-size (unlike for the blob, see below) or acceleration/cooling processes that are generally assumed to be the origin of the faster variability, but rather traces for example injection and/or decay phases of the central engine.

– The existence of strong correlation between X-rays and VHE gamma-ray bands indicates that the observable emission in these two wavebands originates from a single emission region. Therefore, the maximum size of the blob emission region was calculated from the VHE gamma-ray or X-ray variability timescale using the causality relation, R <  ct_varδ/(1 + z). The VHE gamma-ray variability timescale for 1ES 1727+502 and 1ES 2344+514 is 24 h, while for 1ES 1959+650 this timescale is 35 min (Sec. 2.1). The X-ray variability timescale for the case of VER J0521+211 and PKS 1424+240 is 24 h (Appendix D).

– The apparent speeds of the jets can be used to derive the Doppler factor of the core, assuming the viewing angle to be known. We did this for VER J0521+211 and PKS 1424+240 assuming viewing angles equal to 3° and 5°. As discussed in Sect. 3.1, three of our sources show subluminal speeds or even no motion, which is common for TeV BL Lacs (Piner & Edwards 2018, and references therein). Therefore, we use the result from Piner & Edwards (2018), who suggest bulk Lorentz factors with values up to 4. We convert this to Doppler factor assuming a jet viewing angle ∼1/Γ and thus δ ∼ Γ.

– The magnetic field strength of the core can be estimated from the VLBI “core shift”-measurements, assuming a conical jet (Blandford & Königl 1979) and equal energy to be carried by particles and the magnetic field, as done in Pushkarev et al. (2012). The median of the magnetic field strength of the core in the sample of 18 BL Lacs of these latter authors is B_core = 0.10 ± 0.01 G. This sample includes 6 TeV BL Lacs: S5 0716+714, OJ 287, BL Lac, OT 081, Mrk 421, and Mrk501. The first four of these objects have a magnetic field strength of B_core ∼ 0.1 G and the last two sources have B_core ∼ 0.4 G. Another way to estimate the magnetic field strength is to consider the cooling timescale of the electrons, which provides a lower limit to magnetic field strength. In high-synchrotron-peaked sources, the observed emission in the hard X-ray band is due to the high-energy tail of the synchrotron emission. Therefore, the variability timescales are directly linked with particle cooling timescales. Bhatta et al. (2018) studied the variability timescale of 13 blazars in hard X-rays. These latter authors reported a hard X-ray variability timescale between ∼5 min and ∼5 h for six TeV BL Lacs. These timescales were calculated using 18 observation epochs. The average of the estimated variability timescales in their work is ∼1 h. We use Eq. (11) described by Bhatta et al. (2018) to calculate the magnetic field strength. We find that for a variability timescale of ∼1 h the magnetic field strength varies between 0.1 and 0.3 G depending on the assumed Doppler factor and redshift. Therefore, we assume a magnetic field strength of between 0.1 and 0.4 G.

As the emission regions are co-spatial, we use the same magnetic field strength for the blob and core component. The magnetic field strength is generally assumed to scale with distance from the central engine as d⁻¹, and so if the blob was closer to the central engine than the core, it would nominally need to have a stronger magnetic field than the core component, of the order of ∼1 G. Tavecchio & Ghisellini (2016) showed that for single-zone models, the magnetic field strengths tend to be significantly lower than the values required for equipartition values and even in two-component models it is very difficult to reproduce the observed SED with the magnetic field strength values of the order of 1 G. There are ways to invoke reduced local magnetic field strengths in jets such as re-connection layers and radial structures of magnetic fields across the jets (see discussion in Nalewajko et al. 2014), but here we are interested in testing whether or not the observations can be modelled with magnetic field strengths obtained from the core shift measurements (see also Sect. 5.1) and without such reduced local magnetic field strengths. However, we note that different effects can change the blob magnetic field (e.g. internal shocks responsible for the emission from the blob and relativistic movement of the blob with respect to the core) and therefore it is not likely that the magnetic field strengths of the two components would in reality be exactly the same.

Finally, we also try to take into account the derived estimations of the relative strengths of the core and blob components in the optical as derived in Sects. 3.2 and 3.3. As the first method only gives a lower limit to the contribution and the second method did not converge in all cases, these constraints are not very strong, but provide clues as to which of these two components dominates the emission in this waveband. We discuss the comparison of the two methods in Sect. 5.5. Based on these assumptions, we leave the following parameters free to vary in the SED fit:

– γ_min, γ_b, and γ_max of the two components: We limited the range for the values to a physically reasonable regime parameter space; that is, γ_min <  10⁴, 10³ <  γ_b <  10⁵, γ_max <  3 × 10⁶, and demand that the values for the core be always lower than those for the blob.

– n_1, blob and n_2, blob: we considered n_1, blob to be always ∼2. Lower spectral index values are traditionally disfavoured, as for lower values the strong radiative losses of the dominant high-energy electrons would lead to a substantial pressure decrease along the jet and prevent the shock from propagating far out (Marscher & Gear 1985). We also assumed n_2, blob − n_1, blob >  0.5.

– n_1, core and n_2, core: we first considered n_1, core to be always ∼2, but this did not reproduce the shape in the radio part of the SED. Radio observations (e.g. Valtaoja et al. 1988; Hughes et al. 1989) suggest that hard spectral indices are common in AGN and there are also theoretical models (e.g. Stern 2003; Virtanen & Vainio 2005) that can produce indices significantly harder than 2, so we decided to consider values n_1, core >  1.6, which seemed to reproduce the shape of the archival data better. We assumed n_2, core − n_1, core >  1.0.

– The electron energy density normalisation factor K: we limited the range for the values to 10² − 10⁴ cm⁻³ and considered only models where K_blob >  K_core.

– Doppler factor of the blob: we limit ourselves to δ_blob <  30.

We selected eight SED data sets based on the availability of (quasi-)simultaneous data and the observed flux variability at VHE gamma rays. The details of the MWL data selection for each data set is presented in Appendix D. With the observational and theoretical constraints listed above, we check if we can find a set of parameters that reproduces these observed SEDs. Figure 8 shows that in all of the eight cases we find a set of parameters (listed in Table 6) which produces a two-component model in good agreement with the (quasi)-simultaneous observational data. There are some “common trends” in these parameters. In all cases, γ_min, blob is high (> 10³) and n_1, blob is hard. The γ_b, blob varies from 3 × 10⁴ to 9.5 × 10⁴, while the γ_max, blob values fall into a range that is one order of magnitude higher. The n_2, blob values also spread over a large range from 2.45 to 3.85. Both for the cores and the blobs γ_min and n₁ values used in all sources are similar. Interestingly, for the cores, in all but one case a power-law electron distribution without a break was used.

Fig. 8. Broadband SED of the source sample during the selected observation epochs/time. The details of the data selection for each SED is presented in Appendix D. The spectral data points in the panels are: archival non-simultaneous data from the ASI Space Science Data Centre, grey open circles; radio data (15 GHz) from OVRO, blue circle; optical (R-band) from Tuorla, red square; optical and UV data from Swift-UVOT, black stars; X-ray data from Swift-XRT, brown diamonds; HE gamma-ray data from Fermi-LAT, red open circles; and de-absorbed VHE gamma-ray data from MAGIC, blue triangles. The SEDs are modelled within the one-zone SSC (green dotted lines) and two-component scenario (black lines). Within the two-component scenario, the violet dash double dotted and purple dashed lines show the emission from the core and the blob, respectively. Moreover, the result of interaction between emissions from the blob and the core are plotted with blue dash-dotted lines.

The applied model is not time dependent, and so all epochs were modelled independently. We only aim to test the model on “snapshot SEDs” and acknowledge that the fast blob would exit the core region at some point. Therefore, in this model setup the observed changes in the SED can be produced for example by exiting of one blob component and entering of a new one. There are two sources for which we had multiple epoch SEDs. For PKS 1424+240, the two different SEDs are characterised by a lower synchrotron flux in 2015 compared to 2014, while the gamma-ray flux was not changing. We modelled this by decreasing the γ_min, blob and K_blob and by softening the electron energy density spectral index n_1, blob. In the case of 1ES 1959+650, the X-ray and the VHE gamma-ray data of the SED changed significantly. In our models, these parts are largely dominated by the blob emission, for which we altered almost all parameters between the different states, but we also had to alter the core parameters to achieve good representation of the observed SEDs. Finally, we note that the co-spatiality of the emission regions means that we have to consider possible gamma-gamma absorption between the core seed photon field and the highest energy photons emitted by the blob. Our calculations showed that the absorption is negligible.

4.2. One-zone model

In previous works, the SEDs of the sources of our sample were all modelled with one-zone SSC models. The data sets used in those works are not always the same as the ones presented here. As the sources are variable, also parameters used to reproduce the SEDs would vary from one epoch to another. Therefore, for comparison purposes, we also modelled the same SEDs using the one-zone SSC model. The model applied here is the one of Maraschi & Tavecchio (2003), which is the same model used as the basis of the two-component model. We kept the same physically motivated range (but not the same parameters) for electron distribution of the emission region as we used for the blob in the two component model. We also applied the same constraint from the variability timescale for the size of the emission region as we applied for the blob component in the two-component model. We then searched for a set of parameters that described the optical to VHE gamma-ray part of the SED well. The resulting parameters are shown in Table 6 and are similar to those derived in previous one-zone modelling for these sources (see Appendix D). The broadband SEDs, including the one-zone SSC models, are illustrated in Fig. 8. We discuss the differences in parameters and appearance of the SEDs in the following section.

5. Discussion

In the first part of the discussion we compare our observationally constrained two-component models to other modellings of the SEDs (Sect. 5.1–5.3). We then discuss the SED classification of our sample sources (Sect. 5.4) and finally compare the flux ratios of the two components with different methods (Sect. 5.5).

5.1. Comparison of the two-component and one-zone models

In Sect. 4, we model the observed SEDs with the two-component model and one-zone SSC model for comparison purposes. One has to be careful when comparing the two models as we did not perform extensive scans of parameter space that could reproduce the SED. However, some general comparison can be done. For the two-component model we selected the parameters in a way that was as observationally motivated as possible, while for the one-zone SSC model, we looked into previous one-zone SSC models and modified the parameters to fit the data in our epoch (VER J0521+211: Archambault et al. 2013; PKS 1424+240: Aleksić et al. 2014, Kang et al. 2016, Cerruti et al. 2017; 1ES 1727+502: Archambault et al. 2015; 1ES 1959+650: Tagliaferri et al. 2003, 2008, Gutierrez et al. 2006, MAGIC Collaboration 2020a; 1ES 2344+514: Aleksić et al. 2013, MAGIC Collaboration 2020b).

The one-zone SSC models describe the observational data from optical to VHE gamma-rays well, while the radio part of the SED is ignored (the radio part originates from another component as small emission regions are optically thick to radio emission). In all cases, the magnetic field strength is smaller than the one in two-component models, which is in line with very low magnetic field strengths (∼0.01–0.001 G) typically used in one-zone SSC models in the literature. Tavecchio & Ghisellini (2016) showed that this is a general result for one-zone SSC models (see also Sect. 5.2). Also the Doppler factors used in one-zone SSC models are in all cases higher than two-component models, again in line with values (20 ≤ δ ≤ 150) reported in the literature for these sources.

In our two-component model, the blob dominates the emission for all sources from X-ray to VHE gamma-rays, while the core dominates the emission in the radio band. There are two important consequences of taking the core component into account for the modelling of the total SED. Firstly, the core component extends always to the optical band, which then constrains the flux of the low-energy part of the blob component to lower values. This in general forces us to use a relatively narrow electron energy range for the blob component with high γ_min values. Secondly, if the two zones are co-spatial, as is the case here, the core component provides seed photons for the Compton scattering, which then relaxes the requirement of very low B and high δ values for the blob component.

As discussed in Sect. 4.1, the magnetic field strength B ∼ 0.1 G is used for both core and blob components. In our model, the blob is either co-spatial to the core or closer to the central engine than the core. Therefore, it is difficult to motivate how the magnetic field of the small blob would be orders of magnitude weaker than that of the core. We emphasise that we have assumed the same B for the two components, but we cannot exclude B values significantly lower than those used for the blob. On the other hand, B values significantly higher (B ∼ 1 G) are easy to imagine; they could be caused by some dynamo effect or simply compression of the field in shocks as the faster blob moves through the slower core, but with these high B-values we could not reproduce the SEDs. However, we demonstrate that with B values close to those derived by Pushkarev et al. (2012) from the VLBI observations (B ∼ 0.1 G), we can actually reproduce the observed SEDs well from radio to VHE gamma rays.

As one of the main advances of the observationally constrained two-component model with respect to the one-zone SSC model is the inclusion of the radio part of the SED in the modelling, it is rather unfortunate that each SED in this part is constrained only with one quasi-simultaneous data point at 15 GHz. There are some archival data at higher radio frequencies from the Planck satellite at 100–217 GHz. As shown in Fig. 8, our two-component modelling also reproduces this part acceptably for two sources (1ES 1727+502 and 1ES 2344+514) while in the case of the three remaining sources (VER J0521+211, PKS 1424+240 and 1ES 1959+650) the model overproduces the flux in this range. As discussed in Sect. 4.1, it seems that to produce the shape of the SED correctly in this range, n_1, core ≲ 1.6 (or inclusion of further components) is required. However, simultaneous data in this band would be required for further investigation of this issue and therefore it is beyond the scope of this paper.

5.2. Equipartition in the two-component model

Tavecchio & Ghisellini (2016) recently showed that typical SED parameters of one-zone models for BL Lacs are far from equipartition, with the kinetic energy of the particles dominating the energy carried by the magnetic field by several orders of magnitude. These latter authors also showed that solutions close to equipartition can be found for two-component models. Therefore, we used equations 5, 9, 16 and A1 from Tavecchio & Ghisellini (2016) to check the ratio between the magnetic-field and accelerated relativistic-electron energy densities $U_{\text{B}}^{\prime }/U_{\text{e}}^{\prime }$ for our models. We note that for the two-component model, this is calculated for the whole system (see Tavecchio & Ghisellini 2016). The values are reported in Table 7. In all cases, the model parameters of the two-component model suggest solutions close to equipartition, with $U_{\text{B}}^{\prime }$ slightly dominating. The only exception to this is found for the VER J0521+211, where $U_{\text{e}}^{\prime }$ is slightly dominating. For one-zone SSC, the solutions are far from equipartition in all cases, with the kinetic energy of the particles dominating the energy carried by the magnetic field by several orders of magnitude. In conclusion, our results are in agreement with the conclusion of Tavecchio & Ghisellini (2016), that is, solutions close to equipartition can be found in two-component models.

We calculated the kinetic energy associated with the electrons, cold protons, and magnetic field of the core (Table 7) using equations 1 to 3 reported by Celotti & Ghisellini (2008). We performed the calculation for one-zone and two-component models. We note that the two ISP sources (VER J0521+211 and PKS 1424+240) have significantly higher kinetic energies than the three HSP sources (1ES 1727+502, 1ES 1959+650 and 1ES 2344+504). Comparing the results of our sample to the sample reported by Celotti & Ghisellini (2008), we find that this is in line with results on larger samples and that the values we estimated for our sample are similar to those obtained for larger BL Lac samples.

5.3. Comparison to other models

Among the sources of the sample, there has been multiple attempts to model the broadband SED of PKS 1424+240 and 1ES 1959+650 in several frameworks different from that of one-zone SSC modelling, which includes the VHE gamma-ray data. In this section, these results are briefly discussed.

PKS 1424+240. Since the discovery of PKS 1424+240 at VHE gamma rays in 2009 (Ong 2009; Teshima 2009) and the first firm lower limit for the redshift of the source (Furniss et al. 2013), it has become an interesting source for many authors, being the most distant TeV BL Lac object so far detected (z = 0.6; Rovero et al. 2016; Paiano et al. 2017). The source has been monitored since 2009 and there have been several attempts to model the broadband SED using different sets of VHE gamma-ray observations. Kang et al. (2016) and Cerruti et al. (2017) used the VHE gamma-ray data obtained with VERITAS during 2009 and 2013 (Archambault et al. 2014). MAGIC observations during 2009-2011 were used by Aleksić et al. (2014) to build the broadband SED of the source.

In the leptonic scenario, external IC models were tested assuming external photons from a broad-line region (Kang et al. 2016) and/or a dusty torus (Kang et al. 2016; Cerruti et al. 2017). While the model reproduces the observed SED acceptably, there is no observational evidence for a broad-line region and/or a dusty torus in the optical spectra of PKS 1424+240. Hadronic models with the gamma-ray component dominated by pure proton-synchrotron emissions did not lead to any reasonable solution. However, in the framework of a lepto-hadronic scenario, good solutions were found assuming that synchrotron emission from secondary particles was responsible for producing the VHE gamma-ray emission and proton-synchrotron emission produces the radiation at lower energies (Cerruti et al. 2017).

Aleksić et al. (2014) used a two-component model for PKS 1424+240. The model is the same as used in our work, but the two emission zones are not assumed to be co-spatial. The differences in the parameters between the present study at the latter-mentioned study are the following: (i) A lower Doppler factor for the blob is used in this work (18–20 vs. 30); and (ii) the magnetic field strength is higher than the one used in Aleksić et al. (2014) by a factor of three. This comparison demonstrates that while in this work we have sought solutions close to equipartition, the two-component model can reproduce the SED also with parameters that are closer to parameters typically needed for one-zone models (low magnetic field strengths and high Doppler factors).

1ES 1959+650. The source was first detected at VHE gamma rays using the Utah Seven Telescope Array in 1998 (Nishiyama 1999). Being one of the bright, nearby objects among the extragalactic VHE gamma-ray emitters, the source has been observed frequently since its discovery. After the detection of an orphan VHE gamma-ray flare during 2002 (Krawczynski et al. 2004), hadronic models were motivated for describing the broadband SED of the source. MAGIC Collaboration (2020a) applied the proton-synchrotron scenario to describe the broadband SED of the source during the 2016 flaring activity (on MJD 57552 and 57553). These latter authors found that this model requires a high value for the magnetic field strength (B = 150 G) and an acceleration efficiency close to the theoretical limit (η_acc = 1). Unlike the pure proton-synchrotron scenario, the lepto-hadronic models can use lower magnetic field strengths to describe the broadband SED of the source during the flaring activity (Reimer et al. 2005; Bottacini et al. 2010; Sahu et al. 2013; MAGIC Collaboration 2020a).

In the leptonic scenario, the external IC model was tested by Aliu et al. (2013). The external Compton component was motivated by the existence of dust in the central environment of 1ES 1959+650 (Fumagalli et al. 2012). Aliu et al. (2013) found that the external-Compton scenario is able to describe the observational data during the low VHE gamma-ray state. Moreover, they found that this setup was able to generate the anti-correlated X-ray variability seen during the 2007-2011 observation campaign. However, this model needed a relatively low magnetic field (B = 0.02 G) and a high Doppler factor (δ = 30) as typical for one-zone models (see Sect. 5.1).

It has been suggested that the simple one-zone SSC model cannot reproduce the multiple flaring activity of 1ES 1959+650 and multiple-zone SSC models are favoured (Krawczynski et al. 2004; Patel et al. 2018). Those models are different from the one discussed here. Krawczynski et al. (2004) did not take into account the interaction between the two emission zones while Patel et al. (2018) did not assume the co-spatiality of the emission regions. Finally, in both of these studies, the strength of the magnetic field was one order of magnitude lower than the value obtained from VLBI observations assuming equipartition.

5.4. Peak frequencies of the SED

All of the sources in our sample are TeV BL Lacs, but these have quite different SED peak frequencies. The observed synchrotron and IC peak frequencies are estimated from the broad-band SED models (local maxima in the model, first peak corresponding to ν_sync and second to ν_IC; in case of ν_sync we consider the peak to be the higher of the two) presented in Sect. 4. The results are summarised in Table 7.

VER J0521+211. The synchrotron and IC peaks are located at ν_sync = 2.40 × 10¹⁴ and ν_IC = 3.65 × 10²⁴ Hz, respectively. The ratio of the luminosity at the IC peak to that at the synchrotron peak (CD, Compton dominance parameter) is equal to CD = 1.71. These values are in line with the current classification of the source as an ISP BL Lac object (Ackermann et al. 2011), while Fermi-LAT Collaboration (2019) classifies the source as HSP (with ν_sync = 1.40 × 10¹⁵, so very close to ISP-HSP borderline value).

PKS 1424+240. The source is reported to be a HSP (Archambault et al. 2014; Fermi-LAT Collaboration 2019). However, most of the observations which led to this conclusion were obtained during the relatively high state of the source. We found that the source is rather an ISP BL Lac object with ν_sync = (4.68 − 4.99) × 10¹⁴ Hz during its quiescent state in 2014 and 2015. This is in line with the location of the synchrotron peak reported by Aleksić et al. (2014, see Fig. 3) using the data set of the 2010 campaign. Such a transition synchrotron peak also occurs in other BL Lacs (MAGIC Collaboration 2018b).

1ES 1727+502. In our SED the synchrotron and IC peaks are located at ν_sync = 1.92 × 10¹⁸ and ν_IC = 7.48 × 10²⁵ Hz, respectively. It is notable that the source was in high state at VHE gamma rays during the 2015 campaign. This is significantly higher than what was derived by Nilsson et al. (2018) and Fermi-LAT Collaboration (2019) for this source (ν_sync = 2.2 × 10¹⁶ Hz and ν_sync = 7.1 × 10¹⁵), respectively, using archival data. It is also visible in Fig. 8 that the ν_sync has clearly shifted to higher frequency during this high state.

1ES 1959+650. As discussed in Sect. 2, all of the SEDs were observed during high states with the ν_sync evolving from 4.02 × 10¹⁷ to 3.58 × 10¹⁸ Hz from the lowest to the highest state. As clearly visible in Fig. 8, in the archival data the ν_sync was clearly lower, Nilsson et al. (2018) and Fermi-LAT Collaboration (2019) estimated 5.0 × 10¹⁶ Hz and ν_sync = 9.0 × 10¹⁵, respectively.

1ES 2344+514. Also for this source, it is clearly visible in Fig. 8 that in the archival data ν_sync was lower than during the campaign modelled in our work. For the archival data, Nilsson et al. (2018) estimated 5.0 × 10¹⁶ Hz and Fermi-LAT Collaboration (2019) estimated ν_sync = 1.6 × 10¹⁶, while from our modelling we obtain ν_sync = 4.61 × 10¹⁸ Hz.

In summary, our sample includes two ISPs and three HSPs, with peak frequencies shifted to higher frequencies during the flaring states studied in this paper. This is typical behaviour for blazars (Petry et al. 2000; H. E. S. S. Collaboration 2013; Ahnen et al. 2016; MAGIC Collaboration 2018b). Even if the sample is small and therefore not conclusive, we note that the SED parameters γ_max and n₂ we used are more similar between the two ISPs and the three HSPs. Also, the two ISPs are the most luminous sources in our sample (see Table 7). The CD parameters of the SEDs are typical for BL Lacs CD∼1.0, Nalewajko & Gupta 2017) and there is no significant differences in terms of CD between the ISPs and HSPs in the sample.

Even if the sample is small – only five sources –, it still represents different SED peak frequencies, and both flaring and quiescent states. This implies that the model applied here should be applicable to a wide range of BL Lacs, with the exception of very fast, very bright VHE gamma-ray flares (MAGIC Collaboration 2019), where acceptable representation of the observed SED was only found with low magnetic field strength values even within the two-component model.

5.5. Contribution of the two components in the optical band

One of the main goals of this paper was to constrain the contribution of the two components from the long-term data and use that as an input in the SED modelling to limit the parameter space for two-component models even further. This was not entirely successful. Instead, we decided to compare the resulting ratio of the two components given by the long-term light-curve decomposition method (Sect. 3.2), the polarisation study (Sect. 3.3), and the SED modelling (Sec. t4.1). For the SED modelling, the ratio is calculated as the ratio between the core optical flux and the observed optical flux. In the SED modelling, the constant component usually dominates in the optical band, 1ES 1959+650 being an exception with a rather even ratio. The comparison between the three ratios is presented in Table 8.

For four sources (PKS 1424+240, 1ES 1727+502, 1ES 1959+650 and 1ES 2344+514) the lower limit derived using the first method matches the range given by the polarisation study, while for VER J0521+211 it is larger. As can be seen from Table 5, the flux of the constant component is not well constrained for this source; the posteriori (Col. 3) just fills the priori (Col. 2). The low ratio in Table 8 is then the result of the priori (i.e. the minimum optical flux in the whole period), being much lower than the optical flux during the flare in 2013 for which the SED is constructed (see Appendix D). In all cases, the ratios derived from the SEDs are larger than or comparable to the lower limits derived with the first method (see below). Only for two sources (1ES 1727+502 and 1ES 1959+650) is the ratio derived from the SED modelling within the range derived from the optical polarisation study. For all other sources, the ratio derived from the SED modelling is always larger than the one from the optical polarisation.

The first method gives the overall minimum fraction throughout the five-year period of the slowly variable component’s contribution to the variability.

It is important to understand that the three derived fractions are not directly comparable. The first method gives the overall minimum fraction throughout the five-year period of the slowly variable component contribution to the variability. As discussed in Lindfors et al. (2016) it is not possible to constrain the contribution at a given time to the snapshot SED, while that is exactly what is calculated from the polarisation study and the SED modelling. Moreover, the SED modelling obviously does not give an independent result, as we indeed tried to use the results from the other two methods to guide us on the relative strengths of the two components in the optical band. Also, as discussed in Sect. 3.3, for the polarisation study the posteriori was very wide in several cases, which of course results in a very large range for the ratio. Unfortunately, for the source PKS 1424+240, for which the polarisation study was performing best, we failed to derive constraints from the radio–optical long-term study. We intend to further investigate these methods with a larger sample of sources in future work.

6. Summary and conclusions

Here we present our study of the broadband SED of five VHE gamma-ray emitting BL Lacs. We studied the broadband SEDs of these sources during eight observation epochs based on the variability of the VHE gamma-ray flux of the objects. Traditionally, in view of their relative simplicity and agreement with the data, single-zone SSC models have been used to describe TeV BL Lac SEDs (Tagliaferri et al. 2003, 2008; Gutierrez et al. 2006; Archambault et al. 2013, 2015; Aleksić et al. 2014; Kang et al. 2016; Cerruti et al. 2017). However, in almost all cases the radio parts of the SEDs are excluded from the one-zone SSC model based on the argument that the radio emission originates from an outer region. Moreover, the application of a one-zone SSC model requires a low magnetic field strength and a high Doppler factor in order to reproduce the IC part of the SEDs. This results in jets that are far from equipartition (Tavecchio & Ghisellini 2016). The application of multiple-component SSC models has therefore been suggested by many authors (e.g. Krawczynski et al. 2004; Aleksić et al. 2014; Tavecchio & Ghisellini 2016; Patel et al. 2018). We used a two-component model with two interacting, co-spatial emission regions (core and blob) in order to mimic the so-called spine-layer SSC model described by Tavecchio & Ghisellini (2016).

We collected MWL data of the sample in the time-span between 2012 September 30 and 2018 October 9. We tried to constrain the contribution of each emission zone to the total optical flux employing two different approaches. In the first approach we used the method described in Lindfors et al. (2016) to see if the observed trends in radio and optical light curves were temporary and calculated the minimum fraction of the optical (R-band) flux that originates from the same region as the radio emission at 15 GHz. In the second approach, we used optical polarisation data to distinguish between constant and variable components. Moreover, we searched for the existence and time lag of the flaring activity between the radio, optical, and X-ray flux of the sample following the method described by Max-Moerbeck et al. (2014a) and Lindfors et al. (2016).

In order to reduce the parameter space of our two-component model, we used observational constraints from VLBI data. In particular, the magnetic field strength, the size, and the Doppler factor of the core emission region are derived from the VLBI data using the values obtained from the MOJAVE programme (Piner & Edwards 2004, 2018; Piner et al. 2008, 2010; Tiet et al. 2012; Pushkarev et al. 2012; Lister et al. 2019).

We find parameter sets to describe the broadband SEDs of the sample during the selected epochs. The results of our study can be summarised as follows:

• All of the sources of the sample show variability in at least one of the studied bands. Intra-night variability on the nights of 2016 June 13 and 2016 July 1 was detected at VHE gamma rays during the 2016 observation campaign of 1ES 1959+650 (see Table 3 in MAGIC Collaboration 2020a for details). Two of the sources (1ES 1727+502 and 1ES 2344+514) show significant flux variability on the timescale of one day at VHE gamma rays. The other two sources, VER J0521+211 and PKS 1424+240, show daily variability in the X-ray and optical bands, respectively. These variability timescale constraints were applied to the SED modelling.

• The uncertainty of the data points in the HE gamma-ray light curves of the sources was rather large and prevented us from cross-correlation studies. The HE gamma-ray spectral parameters (for the selected epochs) are generally similar to those reported by Acero et al. (2015) and Fermi-LAT Collaboration (2020).

• We find that the significant increasing and/or decreasing trends in the optical band have persisted in four sources out of five (except for PKS 1424+240). Moreover, the increasing and/or decreasing trends in the radio light curves have persisted for all of the sources in the sample. Taking into account the similar findings reported by Lindfors et al. (2016) for the time-span of 5 years before the starting date of our light curves, the radio light-curve trends have persisted for ∼10 years.

• Among the 15 cross-correlation pairs, we only find two significant correlations between radio and optical for 1ES 1727+502 and 1ES 1959+650 with a time lag compatible with zero days, suggesting that in these two sources the emission in these two bands originates from the same region.

• Optical and VLBI core polarisation angles seem to also align in our sources, which is in agreement with the earlier studies (e.g. Hovatta et al. 2016) and with a common origin for optical and radio emission as assumed in our two-component SED modelling.

• The polarisation analysis provided limited additional constraints to the contribution of the two emission components in the optical band. We attempted to model the optical polarisation with a simple “variable” and “constant” component model, but the limitation of the available data lead to a large range on the calculated flux of the constant component. With the exception of VER J0521+211, the lower limit derived by the long-term radio–optical analysis method matches the range given by the polarisation study.

• In all cases, the ratios between the two components in optical derived from the SEDs are larger than or comparable with the lower limits derived with the long-term radio–optical analysis method.

We presented the first systematic attempt to model the broadband SEDs of BL Lacs considering all the observational constraints provided by radio and optical observations: VLBI, long-term variability, and polarisation. The modelled SEDs are associated with ISP and HSP sources, in low and flaring states, and in all of the cases we find a model that is in good agreement with the observed SED. This implies that the model should be applicable to a large fraction of the BL Lacs. Moreover, within the selected two-component scenario, where the emission regions are co-spatial and located at the VLBI core, which is several parsecs away from the central engine, it is possible to reproduce the SEDs with magnetic field strengths and Doppler factors that are well in agreement with the values derived directly from the VLBI observations. This demonstrates that the usually neglected radio component does not have to originate from a region far away from the region dominating the emission in the X-rays and VHE gamma-ray bands and when this component is properly taken into account, we can reproduce the observed SEDs with parameters where the whole system is close to equipartition.

The two-component model presented in this work can be further validated by future observations. Even if the data presented in this work are comprehensive, they remain sparse. The MWL data set and therefore the model can be improved by new MWL simultaneous observations from the radio to the VHE gamma-ray band. In particular, simultaneous radio, millimetre, infrared, and optical data are needed to better constrain the contribution of the emission regions in lower energy bands. Moreover, the energy band where the IC component starts dominating the synchrotron component (hard X-ray to MeV) has to be studied in more detail. Furthermore, better estimations of the IC emission in the energy range between several GeV and 100 GeV are needed to accurately constrain the IC peak frequency. Finally, the connection between radio, optical, and X-ray polarisation should be studied to understand the topography of the magnetic field. In order to achieve these purposes, simultaneous data including radio-millimetre, mid-infrared (JWST, Euclid), X-ray (NuSTAR, IXPE, Spektr-RG), MeV (e-Astrogam or AMEGO), and HE and VHE gamma-ray (CTA) are of interest.

Acknowledgments

We would like to thank the anonymous reviewer whose comments have helped us improving this manuscript. VFR and EL were supported by Academy of Finland projects 317636 and 320045. We would like to thank the Instituto de Astrofísica de Canarias for the excellent working conditions at the Observatorio del Roque de los Muchachos in La Palma. The financial support of the German BMBF and MPG; the Italian INFN and INAF; the Swiss National Fund SNF; the ERDF under the Spanish MINECO (FPA2017-87859-P, FPA2017-85668-P, FPA2017-82729-C6-2-R, FPA2017-82729-C6-6-R, FPA2017-82729-C6-5-R, AYA2015-71042-P, AYA2016-76012-C3-1-P, ESP2017-87055-C2-2-P, FPA2017-90566-REDC); the Indian Department of Atomic Energy; the Japanese ICRR, the University of Tokyo, JSPS, and MEXT; the Bulgarian Ministry of Education and Science, National RI Roadmap Project DO1-268/16.12.2019 and the Academy of Finland grant nr. 320045 is gratefully acknowledged. This work was also supported by the Spanish Centro de Excelencia “Severo Ochoa” SEV-2016-0588 and SEV-2015-0548, the Unidad de Excelencia “María de Maeztu” MDM-2014-0369 and the “la Caixa” Foundation (fellowship LCF/BQ/PI18/11630012), by the Croatian Science Foundation (HrZZ) Project IP-2016-06-9782 and the University of Rijeka Project 13.12.1.3.02, by the DFG Collaborative Research Centers SFB823/C4 and SFB876/C3, the Polish National Research Centre grant UMO-2016/22/M/ST9/00382 and by the Brazilian MCTIC, CNPq and FAPERJ. The Fermi-LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l’Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d’Études Spatiales in France. This work performed in part under DOE Contract DE-AC02-76SF00515. Part of this work is based on observations made with the Nordic Optical Telescope, operated by the Nordic Optical Telescope Scientific Association at the Observatorio del Roque de los Muchachos, La Palma, Spain, of the Instituto de Astrofisica de Canarias. We acknowledge the support from the Bulgarian NSF through grants DN 18-1220 13/2017, DN 18-10/2017, KP-06-H28/3 (2018) and KP-06-PN38/1 (2019). This research has made use of data from the OVRO 40-m monitoring program which is supported in part by NASA grants NNX08AW31G, NNX11A043G, and NNX14AQ89G and NSF grants AST-0808050 and AST-1109911. This research has made use of data from the MOJAVE database that is maintained by the MOJAVE team (Lister et al. 2019). Part of this work is based on archival data, software, or online services provided by the Space Science Data Centre – ASI. This research has made use of data and/or software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC and the High Energy Astrophysics Division of the Smithsonian Astrophysical Observatory.

References

Appendix A: Examples of online data

Examples of data available at the CDS are presented in Tables A.1 and A.2 following the analysis procedures described in Sects. 2.4 and 2.3, respectively.

Appendix B: Light-curve correlations

In this section the results of the cross-correlation analysis (Sect. 3.2.2) are shown in Figs. B.1–B.5 for each source.

Fig. B.1. Results of the LCCF study of VER J0521+211: (top) between radio (15 GHz) and optical (R-band); (middle) between radio (15 GHz) and X-ray (2–10 keV); and (bottom) between optical and X-ray (2–10 keV). We show 1σ (green solid line), 2σ (blue dashed line), and 3σ (red doted line) significance limits. Positive significant lags show that the flare at lower frequency is preceding the other band.

Fig. B.2. Same description as in Fig. B.1 for PKS 1424+240.

Fig. B.3. Same description as in Fig. B.1 for 1ES 1727+502.

Fig. B.4. Same description as in Fig. B.1 for 1ES 1959+650.

Fig. B.5. Same description as in Fig. B.1 for 1ES 2344+514.

Appendix C: Details of the polarisation modelling

Here we provide a detailed description of the model parameters and fitting procedure we used in Sect. 3.3. As discussed there, we modelled the variable polarisation component as a homogeneous cylindrical emission region in a jet with a helical magnetic field and computed the Stokes parameters using the formulae in Lyutikov et al. (2005).

The relativistic outflow with speed β is passing through this cylindrical region, with length l and radius r whose symmetry axis is parallel to the velocity vector $\overline{v}=vz$. This vector forms an angle θ with the line of sight to the observer and a sky projected angle φ₀ with the vector pointing to the North in the sky. The cylinder has a uniform electron density K_e and a power-law energy distribution of electrons with a slope p. The magnetic field in the cylinder has two components, one parallel to $\overline{v}$, and another one perpendicular to it, i.e.

$$\begin{array}{c}\hfill \overline{B\prime }=(B{\prime }_x,B{\prime }_y,B{\prime }_z)=B_0[sin\psi \prime (-sin\varphi ,0,cos\varphi )+cos\psi \prime (0,0,1)],\end{array}$$(C.1)

where ψ′ is the magnetic field pitch angle, ϕ is the azimuthal angle in the plane perpendicular to $\overline{v}$ (z-axis), and we follow the notation in Lyutikov et al. (2005) by marking quantities B and ψ in the co-moving plane as primed. The model for the variable component now has 9 parameters: β, θ, φ₀, l, r, K_e, p, B₀, and ψ′. The total number of parameters to model the variable and constant components is 12; the constant component has three free parameters: I_C, Q_C, and U_C.

Most of the parameters in the model cannot be constrained with monochromatic observational data due to a high degree of degeneracy in the model. Since the variable (blob) emission region is assumed to be homogeneous and unresolved, Eq. (2) of Lyutikov et al. (2005) can be reformulated as follow and the Stokes parameter I becomes

$$\begin{array}{c}\hfill I=\mathrm{const}\times \frac{p+7/3}{p+1}\frac{{\delta }^{2+(p-1)/2}}{D^2{(1+z)}^{2+(p-1)/2}}\frac{K_e\pi r^{2}l}{sin\theta }\int |B\prime sin\psi \prime |\mathrm{d}S,\end{array}$$(C.2)

where D is the luminosity distance of the source and the integration is over the volume of the cylinder. There are infinite combinations of β, θ, p, K_e, r, l and B₀ which produce the same I (and Q and U), so there is no way to constrain these parameters.

To deal with the degeneracy, the fit was made in the Q − U plane following these assumptions and the same procedure: we fixed p to 2.1, r to 2.5 × 10¹⁵ cm, l to 5 × 10¹⁵ cm, B₀ to 0.1 Gauss, and β to 0.99. These values are similar to those applied for the SED modelling (Sect. 4.1). At each iteration, we then determine the K_e that equals the model I to the observed I value for each data point. We use this value of K_e to compute Q_V and U_V (see Eq. (2) in Lyutikov et al. 2005). We thus assumed that changes in the I_V were due to changes in K_e. This also means that we effectively fit two linearly polarised components with constant Q and U, but with I of the other (variable) one changing. It is clear from examining the data that this simple model can fit only major trends in the data. There is a lot of fast variability (see Figs. 1–5 panels 2,3 and 4 from the top), which is probably caused by random turbulence and therefore cannot be explained by a deterministic model. We treat this fast flickering as pure noise by adding one more parameter to the model: σ, which is added in quadrature to the errors of the observed Stokes parameters Q and U when we compute the likelihood.

There are therefore seven free parameters in our final model: the constant-component Stokes parameters I_C, Q_C, and U_C, viewing angle θ, magnetic field pitch angle ψ′, jet position angle φ₀, and standard deviation of the turbulence, σ. The fit was made using a Monte Carlo Markov chain (MCMC) ensemble sampler (Goodman & Weare 2010). The posteriori probability was first sampled with 21 walkers using 10 000 steps and best-fit values for the parameters were then obtained from the marginalised distributions. We found that the fit stabilised quite well after 50 iterations. We therefore treated the first 50 iterations as the “burn-in” phase, and discarded them. Since the model always “predicts” I correctly, the posteriori was computed from the Q and U data only. The priors were set up in the following way: I_C: flat from zero to minimum observed flux (given in Col. 3 of Table 5), Q_C and U_C: Gaussian with μ = 0.0 and σ = 0.3 mJy, θ: flat 0 to 10°, ψ′: flat 0 to 90°, φ₀: flat 0 to 180° and σ: Gaussian with μ = 0.1 mJy and σ = 0.1 mJy.

The results of the fits are illustrated in Figs. C.1–C.5 and summarised in Table 5. The error bars give the 68% confidence intervals derived from marginalised distributions. Figure C.5 gives an example of typical posteriori distributions. This figure illustrates some common trends, which we will now discuss.

Fig. C.1. Results of polarisation analysis for VER J0521+211. Panel a: observed optical R-band flux (blue circles), variable component model (red squares), and model constant component (brown diamonds). Panel b: observed (blue circles) and modelled (red stars) optical polarisation degree. Panel c: observed (blue circle) and modelled (red triangles) electric vector polarisation angle. Panels d–f: posteriori distributions of the polarisation fitting. The colour scale gives the number of visits by the sampler in each cell. In addition, panel f shows the observations in the Q − U plane (orange) and evolution of the model (light green).

Fig. C.2. Same description as in Fig. C.1 but for PKS 1424+240.

Fig. C.3. Same description as in Fig. C.1 but for 1ES 1727+502.

Fig. C.4. Same description as in Fig. C.1 but for 1ES 1959+650. The last season was excluded from the polarisation fitting (see text).

Fig. C.5. Same description as in Fig. C.1 but for 1ES 2344+514.

As discussed above, the EVPA of the variable component is always parallel or perpendicular to the projected direction of the jet in the sky. At low pitch angles (ψ′≲56°) the EVPA is perpendicular and at higher angles (ψ′> 56°) parallel to the jet (see Fig. 3 in Lyutikov & Kravchenko 2017). The degree of polarisation is at maximum when ψ′ = 0° or 90° and approaches zero when ${\psi }^{\prime }={\psi }_{\text{crit}}^{\prime }~56°$. This explains why our fitted ψ′ are so similar among our targets; values close to 56° are the only way to produce the observed low degrees of polarisation. The exact value of ${\psi }_{\text{crit}}^{\prime }$ depends on the viewing angle, which explains why the posteriori is confined to such a narrow strip (see upper left panel of Fig. C.5). However, the narrowness of this strip does not mean that we have put strict constraints on ψ′. We have simply forced it to this range by the assumption of a perfectly organised magnetic field. Had we introduced random variations into the magnetic field, the integrated degree of polarisation would have been at the observed level over a wider range of ψ′.

The posteriori for the jet angle φ₀ has typically two branches separated by 90°. In the upper left panel of Fig. C.4 this shows up as two horizontal bands, one of which is much weaker than the other. The weaker branch is visible also in the upper left panel as a small extension towards the lower left corner. This is again a product of the possibility to have two EVPA orientations with respect to the jet axis. Table 5 shows a comparison between our findings from this optical polarisation analysis and the VLBI results listed in Sect. 3.1. The comparison is not straightforward due to the ±180° ambiguity of the EVPA and the 0 or 90° relative orientation if the jet and the EVPA allowed by our model.

We also note that our model fails to fit the EVPA change in 1ES 1959+650 during the last observing season (Fig. C.4). The EVPA changes by ∼45°, which cannot be done by simply changing the relative contributions of the two components, whose properties are at the same time consistent with the earlier periods. In our formulation, the only way is to change φ₀, that is, the orientation of the jet. Therefore, we excluded the last season of data from the polarisation fitting.

Finally, with our present model, the evolution of Q and U is always along a straight line in the Q, U-plane. The observed evolution is much more complex, but we have simplified our treatment by introducing the parameter σ, the standard deviation of the random turbulent variations of Q and U. This parameter adjusts itself during the fit according to how well the rest of the parameters fit the data. If the predictive power of the latter is weak, σ increases to accommodate the increased residuals. This happens automatically because increasing σ increases the likelihood up to a point, after which the likelihood begins to decrease with increasing sigma. As a result, σ adjusts itself to a value that is roughly equal to the value the standard deviation of the residuals. In the case of Fig. C.5, σ = 0.018 mJy, and so adding the turbulent component would mean adding a Gaussian random variable with this σ to every point in the evolution of the model in Fig. C.5, which would scatter the model points similarly to the observed points. As discussed in Sect. 3.3, the addition of parameter σ makes the model very smooth and therefore it does not reproduce the observed complex variability very well. The MCMC loop finds a good fit simply by increasing the errors, implying that the model we selected for the variable component is probably too simple. If our model had more explanatory power, σ would tend to lower values and the model would follow the data points more closely. However, the approach introduced in this work, that is, the Bayesian fitting of physical model to optical polarisation data, will be further investigated in future work.

Appendix D: Details of the SED modelling

Table D.1 summarises the detailed description of the selection procedure used for MWL observational data sets used for SED modelling (Sect. 4). Here we discuss the simultaneity of the modelled data and some details of the model parameters for each source. In general for each case of SED modelling, the radio data point at 15 GHz was selected based on the shortest time lag between observation performed by OVRO and other instruments. All of the UV and optical data points were corrected Galactic extinction and the contribution of host-galaxy flux.

D.1. VER J0521+211

The VHE gamma-ray spectrum was constructed by combining all of MAGIC observations between MJD 56580.18 and 56627.95 (4.7 h distributed over four nights). Fermi-LAT data obtained between MJD 56580 and 56627 were used for building the HE gamma-ray spectrum. The selected observations performed by Swift-XRT and Swift-UVOT instruments were carried out on MJD 56625.30, which is the most simultaneous observation to one of the four MAGIC observation windows (MJD 56625.04-56625.12). Similarly, we used optical R-band data on MJD 56625.15 obtained by the KVA telescope. For this source, the contribution of the host-galaxy flux in the UV and optical bands was neglected based on the uncertainty of the redshift and the reported redshift lower limit. In this period, the source was clearly in high state in the optical, X-ray, and VHE gamma-ray band compared to the previous observations Archambault et al. (2013).

Panel a in Fig. 8 illustrates the broadband SED of VER J0521+211. The two-component model can reproduce the observed quasi-simultaneous data using the parameters sets for the two emission regions reported in Table 6. The core Doppler factor was calculated from the apparent motion seen in the VLBI data assuming jet viewing angle of 5° and a redshift of 0.18. The size of the blob emission region was set to be smaller than the shortest variability timescale in the data set (24 h detected in X-rays band) as the sampling of the light curves would limit our capability to detect shorter timescales of variability. However, due to degeneracy between different parameters, we could achieve equally good representation of the data by increasing R and decreasing K.

The comparison of the SED parameters of the blob in Table 6 with those obtained from single-zone SSC model tested by Archambault et al. (2013), shows that the main differences are: the size of emission region (4.0 × 10¹⁷ cm for one-zone model vs. 1.3 × 10¹⁶ cm for two-component model); magnetic field strength (0.0025 vs. 0.1 G) which affect the jet equipartition; the Doppler factor (30 vs. 12); and the maximum electron Lorentz factor (2.0 × 10⁶ for one-zone model vs. 4.0 × 10⁵ for two-component model). These differences are in line with the general trends discussed in Sect. 5.1.

D.2. PKS 1424+240

As discussed in Sect. 2.1, we attempted to model the SED of PKS 1424+240 using the data obtained from observations during the two campaigns. Neither of the campaigns were performed during particularly high flux states in lower (optical, X-rays) bands.

Data set of the observation campaign during 2014. The VHE gamma-ray spectrum was constructed by combining all of the MAGIC observations in the time-span from MJD 56740.06 and 56825.99. A similar time window was used to compute the HE gamma-ray spectrum from the data obtained with Fermi-LAT. We searched for variability of spectral parameters (F₀, Γ, and β) in HE gamma-ray band in the selected time window. These parameters did not show any significant (at 3σ confidence level) variability over a one-week timescale. In the presented data, there is a mismatch between the HE and VHE gamma-ray spectra at energies between 40 and 60 GeV. This mismatch is mainly due to non-simultaneity of observations. However, considering the systematic uncertainty of both instruments in that range, the mismatch is negligible. At MJD 56800.93, the VHE gamma-ray flux was the most consistent measurement to the average flux of the entire campaign. Therefore, we selected Swift-XRT (MJD 56801.25), Swift-UVOT (MJD 56801.26) and optical (MJD 56800.96) observations which were quasi-simultaneous to MJD 56800.93. Panel b in Fig. 8 shows the broadband SED of PKS 1424+240 (compiled from observations of the 2014 campaign). The parameter sets for the two emission regions can reproduce the observed quasi-simultaneous observational data (see Table 6). The core Doppler factor is calculated from apparent speed observed in VLBI using a jet viewing angle of 3°. The size of the blob emission region is compatible with a variability timescale seen in the optical and X-rays (1 and 2 days respectively).

Data set of the observation campaign during 2015. The VHE gamma-ray spectrum was constructed by combining all of the MAGIC observations performed between MJD 57045.05 and 57186.06. Similar time span was used to compute the HE gamma-ray spectrum from the data obtained with Fermi-LAT. Following the procedure used for observation campaign during 2014, the spectral parameters did not show any significant (at 3σ confidence level) variability at weekly timescale. We selected Swift-XRT (MJD 57135.26), Swift-UVOT (MJD 57135.26) and optical (MJD 57135.11) observations which were simultaneous to one of the MAGIC observation window. Panel c in Fig. 8 presents the broadband SED of PKS 1424+240 during the 2015 campaign. The observed quasi-simultaneous observational data are reproduced by the sets of parameters summarised in Table 6. For the core, we use parameters similar to those in 2014, with minor changes to reproduce the overall lower state of the synchrotron part of the SED. For the blob, n₂ is softer than in 2014 to reproduce the lower state in X-rays. Also, the size of the blob emission region is slightly smaller than in 2014, even if no variability was detected during this campaign, again to reproduce the lower X-ray state.

SED modelling results in optical band. In contrast with the result in Sect. 3.2.1, the emission from the core dominates the total flux in the optical band (F_core/F_total = 0.93). We tried to find sets of parameters in which the optical emission would be dominated by the blob component. Only solutions we found that could produce the HE and VHE gamma-ray part of the SED, had very low magnetic field strength and was far from equipartition, while the set of parameters that we present here give $U_B^{\prime }/U_e^{\prime }=2.39$. We emphasis that the derived value for F_core/F_total in Sect. 3.2.1 is a minimum value, so there is no contradiction between these results. The optical polarisation method (Sect. 3.3) suggests that there are two components contributing to the optical flux, which further highlights the general conclusion in Sect. 5.5 that these two methods are complementary.

D.3. 1ES 1727+502

The VHE gamma-ray spectrum was computed by combining all of the MAGIC observations performed between MJD 57306.83 and 57327.83. The HE gamma-ray spectrum was built using the data obtained with Fermi-LAT in the time span of MJD 57263 and 57353. The Swift-XRT and Swift-UVOT observations were selected from the night when the VHE gamma-ray flux was consistent with average flux, i.e. on MJD 57307.99. The optical data point obtained from the observation performed with KVA telescope on MJD 57307.83.

Panel d of Fig. 8 shows the broadband SED of 1ES 1727+502 during flaring activity in 2015. The observed quasi-simultaneous data can be reproduced in the framework of a two-component model using the sets of parameters reported in Table 6. The size of blob emission region is compatible with variability timescale of 6.3 h which is in agreement with our observational constrain of 24 h detected in VHE gamma-ray band.

The comparison of the SED parameters obtained with a one-zone SSC model used by Archambault et al. (2015) with those for the blob emission region (Table 6) shows similar differences as seen in VER J0521+211. The size of the emission region is (4.3 − 7.4 × 10¹⁷ cm for one-zone model vs. 7.1 × 10¹⁵ cm for two-component model), magnetic field strength (0.0003–0.0006 vs. 0.1 G), the Doppler factor (30 vs. 11), and the maximum electron Lorentz factor (5.5 − 7.0 × 10⁶ for one-zone model vs. 1.3 × 10⁶ for two-component model).

D.4. 1ES 1959+650

Following the discussion in Sect. 2.1, three data sets were used to build the SEDs of this source in different states of VHE gamma-ray flux.

Low state. The VHE gamma-ray spectrum was calculated using the 1.2 h of the MAGIC data from observation starting at MJD 57711.82. The time-span between MJD 57711.43 and 57712.35 was used for building the HE gamma-ray spectrum adopting Fermi-LAT data. For X-ray, UV and optical, the observations which were performed on MJD 57711.58 (Swift-XRT), 57711.59 (Swift-UVOT) and 57711.82 (KVA) were used.

Intermediate state. The VHE gamma-ray spectrum was constructed using 1.4 h of MAGIC data from observation starting at MJD 57547.13. The Fermi-LAT data obtained between MJD 57545.16 and 57550.19 were used for building the HE gamma-ray spectrum. For X-ray and UV, the observations which were performed on MJD 57547.13 (Swift-XRT and Swift-UVOT) were used. The optical and UV data from (MJD 57547.14) were used in SED modelling.

High state. The VHE gamma-ray spectrum was computed using 2.2 h of MAGIC data from observations that started at MJD 57553.06. The Fermi-LAT data obtained between MJD 57552.00 and 57554.00 were used for calculating the HE gamma-ray spectrum. For X-ray and UV, the observations which were performed at MJD 57553.10 (Swift-XRT and Swift-UVOT) were used. The optical data point from (MJD 57553.13) was used in SED modelling.

It should be noted that even the “low state” SED is somewhat above the archival data from previous “low-state” campaigns, as it is clearly visible in Fig. 8. The comparison of one-zone SSC models is discussed in detail by MAGIC Collaboration (2020a), but the differences follow the general trend we discuss in Sect. 5.1.

D.5. 1ES 2344+514

The VHE gamma-ray spectrum was calculated using 0.5 h of MAGIC data from observations that started at MJD 57612.06 in order to have quasi-simultaneous MWL coverage. Due to the source faintness in the HE gamma-ray band, the Fermi-LAT data obtained between MJD 57520.00 and 57704.00 were used for building the spectrum. For X-ray and UV, the observations which were performed on MJD 57613.52 (Swift-XRT and Swift-UVOT) are used. The optical data (R-band) from (MJD 57611.02) was used in SED modelling. Figure 8h, demonstrates the broadband SED of 1ES 2344+514 during the 2016. As can be seen the source is relatively bright compared to archival observations in X-rays, HE and VHE gamma-rays during this period.

The one-zone SSC can describe the low-state SED of the source (Aleksić et al. 2013) using the data sets of the observation campaign during 2008. However, the radio data is not included in that modelling, and is assumed to origin from different component further out. We find that most of the parameters of the one-zone model reported by Aleksić et al. (2013), are in good agreement with the blob parameters in this work. It is notable that the one-zone model used the Doppler factor (δ = 20) and magnetic field strength (B = 0.07 G), while these parameters for the blob emission in this work are more physically realistic (δ = 6 and B = 0.1 G). Moreover, similar quasi-simultaneous data set is used by MAGIC Collaboration (2020b) and the set of parameters describing the broadband SED have low magnetic field strength (B = 0.02 G), emission region size (R = 1 × 10¹⁶ cm), and high break and maximum electron Lorentz factor (γ_b = 1.8 × 10⁶ and γ_max = 8.0 × 10⁶).

All Tables

All Figures

Fig. 1. MWL light curves of VER J0521+211 in the range from MJD 56200 (2012 September 30) to 58400 (2018 October 9). From top to bottom panels: radio and VLBI core flux (15 GHz), optical (R-band), optical polarisation degree, electric vector polarisation angle, X-ray flux (2–10 keV), HE gamma-ray photon flux (0.1–300 GeV), and VHE gamma-ray photon flux above the threshold energy given in the panel. Black arrows show the 95% confidence level upper limits. The data, which are marked with vertical lines/area and squares in different bands, are used in the SED modelling.

In the text

	Fig. 2. Same description as in Fig. 1 for PKS 1424+240.
In the text

	Fig. 3. Same description as in Fig. 1 for 1ES 1727+502.
In the text

	Fig. 4. Same description as in Fig. 1 for 1ES 1959+650.
In the text

	Fig. 5. Same description as in Fig. 1 for 1ES 2344+514.
In the text

Fig. 6. Steps of the analysis used to determine the contribution of the common emission component at radio and optical frequencies. Left: fitting a polynomial to the radio light curve. Middle: polynomial fit scaled to the average flux of the optical light curve and multiplied with a scaling factor (different colours correspond to different scaling factors, see text). Right: scaled polynomial is subtracted from the optical data (green filled circles with errors). The RMS of the resulting light curve (purple filled circles) is compared with the RMS of the original data. These analysis steps are shown for all sources from bottom to top: VER J0521+211, PKS 1424+240, 1ES 1727+502, 1ES 1959+650 and 1ES 2344+514. In the case of PKS 1424+240 subtracting the polynomial did not decrease the RMS of the optical light curve and therefore the purple dots are under the green symbols and are not well visible in the rightmost panel.

In the text

	Fig. 7. Sketch of the geometrical modelling setup. The two emission regions are located several parsecs from the central black hole (at dcore). The smaller emission region (blob) is embedded in the larger emission region (core) and the interaction of the two emission regions provides additional seed photons for the Compton scattering; see discussion in Tavecchio et al. (2011, Appendix B).
In the text

Fig. 8. Broadband SED of the source sample during the selected observation epochs/time. The details of the data selection for each SED is presented in Appendix D. The spectral data points in the panels are: archival non-simultaneous data from the ASI Space Science Data Centre, grey open circles; radio data (15 GHz) from OVRO, blue circle; optical (R-band) from Tuorla, red square; optical and UV data from Swift-UVOT, black stars; X-ray data from Swift-XRT, brown diamonds; HE gamma-ray data from Fermi-LAT, red open circles; and de-absorbed VHE gamma-ray data from MAGIC, blue triangles. The SEDs are modelled within the one-zone SSC (green dotted lines) and two-component scenario (black lines). Within the two-component scenario, the violet dash double dotted and purple dashed lines show the emission from the core and the blob, respectively. Moreover, the result of interaction between emissions from the blob and the core are plotted with blue dash-dotted lines.

In the text

Fig. B.1. Results of the LCCF study of VER J0521+211: (top) between radio (15 GHz) and optical (R-band); (middle) between radio (15 GHz) and X-ray (2–10 keV); and (bottom) between optical and X-ray (2–10 keV). We show 1σ (green solid line), 2σ (blue dashed line), and 3σ (red doted line) significance limits. Positive significant lags show that the flare at lower frequency is preceding the other band.

In the text

	Fig. B.2. Same description as in Fig. B.1 for PKS 1424+240.
In the text

	Fig. B.3. Same description as in Fig. B.1 for 1ES 1727+502.
In the text

	Fig. B.4. Same description as in Fig. B.1 for 1ES 1959+650.
In the text

	Fig. B.5. Same description as in Fig. B.1 for 1ES 2344+514.
In the text

Fig. C.1. Results of polarisation analysis for VER J0521+211. Panel a: observed optical R-band flux (blue circles), variable component model (red squares), and model constant component (brown diamonds). Panel b: observed (blue circles) and modelled (red stars) optical polarisation degree. Panel c: observed (blue circle) and modelled (red triangles) electric vector polarisation angle. Panels d–f: posteriori distributions of the polarisation fitting. The colour scale gives the number of visits by the sampler in each cell. In addition, panel f shows the observations in the Q − U plane (orange) and evolution of the model (light green).

In the text

	Fig. C.2. Same description as in Fig. C.1 but for PKS 1424+240.
In the text

	Fig. C.3. Same description as in Fig. C.1 but for 1ES 1727+502.
In the text

	Fig. C.4. Same description as in Fig. C.1 but for 1ES 1959+650. The last season was excluded from the polarisation fitting (see text).
In the text

	Fig. C.5. Same description as in Fig. C.1 but for 1ES 2344+514.
In the text