© The Authors 2023

1. Introduction

Blazars are radio-loud active galactic nuclei whose relativistic jets point towards the Earth. Blazars can be classified according to the spectral features in the optical band as BL Lacertae (BL Lacs) objects and flat spectrum radio quasars (FSRQs). While BL Lacs have an (almost) featureless optical spectrum, FSRQs show strong, broad emission lines in the optical band (Urry & Padovani 1995). Blazars are also characterised by high variability over very different timescales. They emit mostly non-thermal radiation at all wavelengths, from radio to γ rays. However, most of the blazars detected in very-high-energy (VHE; ≳100 GeV) γ rays are BL Lacs.

Blazars display a spectral energy distribution (SED) characterised by the presence of two bumps (Ghisellini et al. 2017). The first bump originates from synchrotron radiation by relativistic electrons. The second peak is commonly explained by a leptonic scenario through inverse Compton (IC) scattering of synchrotron photons – synchrotron self-Compton (SSC) scattering – with the same electron population (see e.g. Celotti & Ghisellini 2008; Ghisellini et al. 2010) and/or IC scattering of photons coming from outside the jet in an external Compton process (Dermer & Schlickeiser 1994). Alternatively, different models with a hadronic origin have been proposed to explain the high-energy bump in the SED of blazars (e.g. Mannheim 1993; Cerruti et al. 2015). BL Lacs can be divided into three groups depending on the frequency of the synchrotron peak: low- (ν_peak < 10¹⁴ Hz), intermediate- (10¹⁴ Hz < ν_peak < 10¹⁵ Hz), and high-energy-peaked BL Lacs (HBLs; ν_peak > 10¹⁵ Hz; Padovani & Giommi 1995). Another category of BL Lacs was introduced by Costamante et al. (2001), naming those whose peak is above ν_peak > 10¹⁷ Hz extreme HBLs (EHBLs). These objects display a high X-ray flux with respect to their optical/UV emission. They can also show a high-energy peak shifted to VHE γ-ray frequencies (see for instance the BL Lac 1ES 0229+200 and the sources reported by Acciari et al. 2020).

1ES 0647+250 is a BL Lac object previously catalogued as an HBL (Costamante & Ghisellini 2002; Aleksić et al. 2011). It has an uncertain redshift, with various values reported in the literature. A lower limit of z > 0.3 was first derived from imaging of the source by Falomo & Kotilainen (1999). However, based on deep observations of the host galaxy, Meisner & Romani (2010) derived a value of $z=0.{45}_{-0.10}^{+0.11}$, and Kotilainen et al. (2011) estimated a redshift of z = 0.41 ± 0.06 after using the imaging redshift method from Sbarufatti et al. (2005). Spectral lines have not been detected in its spectrum. These non-detections were used to derive lower limits on the redshift of z > 0.3 by Scarpa et al. (2000) and z > 0.47 by Sbarufatti et al. (2005). An accurate measurement of the redshift is still lacking, though the most recent work by Paiano et al. (2017) set a lower limit of z > 0.29.

It was first reported as a VHE γ-ray emitter by the Major Atmospheric Gamma-ray Imaging Cherenkov (MAGIC) collaboration with a flux above 100 GeV of (3.0 ± 0.7)% Crab Nebula flux Units (C.U.; De Lotto & MAGIC Collaboration 2012). Later on, it was detected by the Very Energetic Radiation Imaging Telescope Array System (VERITAS) with a γ-ray flux of (2.7 ± 0.7)% C.U. above 140 GeV as part of the VERITAS blazar programme carried out between 2010 and 2013 (Dumm 2013).

As for most blazars, this source is bright and variable in all the electromagnetic bands and has been observed as part of many programmes in radio (Piner & Edwards 2014), optical (Kapanadze 2009), and X-rays (Perlman et al. 2005). It is also detected at high energies (HEs, 30 MeV < E < 100 GeV), and it can be found in each of the Fermi Large Area Telescope (LAT; Atwood et al. 2009) source catalogues from 1FGL onwards, including the 4FGL catalogue (Abdollahi et al. 2020). Significant variability has been detected in the optical band, but no intra-night or short-burst variability has been claimed for this source. Furthermore, it has not been possible to come to a firm conclusion about the variability timescales of 1ES 0647+250 in the optical band due to long gaps in the historical light curve (Kapanadze 2009). Nilsson et al. (2018) also detected significant variability for this source using optical data from between 2002 and 2012.

In this paper we perform the first long-term multi-wavelength (MWL) study of 1ES 0647+250. We report on the detection in the VHE γ-ray band by MAGIC in four different epochs (2009−2011, 2014, 2019, and 2020), each of which corresponds to a different state of the source in terms of its optical, X-ray, and VHE flux. The MAGIC observations performed between 2009 and 2011 were triggered by the first studies done on Fermi-LAT data above 10 GeV, which later on led to the first Fermi-LAT catalogue of > 10 GeV sources (1FHL; Ackermann et al. 2013), and identified several VHE γ-ray emitter candidates with only one year of LAT data. The VHE γ-ray observations in 2014 were triggered by an optical flare detected by several optical facilities (Kiehlmann et al. 2014). In December 2019, 1ES 0647+250 showed a historically high X-ray flux (Kapanadze 2019), leading to the detection at VHE γ rays of this blazar (Mirzoyan 2019). Finally, the source displayed its highest state in the VHE γ-ray domain in December 2020, after an X-ray activity comparable to the 2019 flare (Kapanadze 2020).

This paper is structured as follows: In Sect. 2 the data sets used in the analysis are introduced. Based on the MWL data collected in this work, in Sect. 3 we present MWL variability studies for the first time from radio to the VHE γ-ray band. In Sect. 4 the spectral analysis of the X-ray and γ-ray data is performed. A redshift estimation of the source based on the γ-ray spectrum is presented in Sect. 5 and compared with previous measurements. In Sect. 6 we model for the first time the broadband SED of this source for the different observed periods and compare them to one another. In Sect. 7 a discussion and interpretation of the results are presented, and the main results are provided as a conclusion in Sect. 8.

2. Multi-wavelength data

2.1. VHE γ rays: MAGIC telescopes

MAGIC is a stereoscopic system of two 17 m imaging atmospheric Cherenkov telescopes located on the Canary island of La Palma, Spain, at an altitude of ∼2200 m above sea level. They work in an energy range between 50 GeV and tens of TeVs, with a sensitivity above 100 GeV (300 GeV) of about 2% (about 1%) of the Crab Nebula flux after 25 h of observations at zenith angle ZA < 30° (Aleksić et al. 2016). These characteristics make the MAGIC telescopes very well suited for blazar observations in the VHE γ-ray range.

1ES 0647+250 was first observed by MAGIC-I in mono mode in 2008 (Aleksić et al. 2011), with no detection of the source. However, an upper limit of the integral flux above 120 GeV of 1.6 × 10⁻¹¹ cm⁻² s⁻¹ was estimated. For the present work, we use approximately 45 h of stereoscopic good-quality data taken by MAGIC between November 2009 and December 2020: 26.7 h after quality cuts between November 2009 and March 2011. In November 2014, the observations (2.2 h after data quality cuts) were triggered under the target of opportunity programme following an enhanced flux in the optical and HE γ-ray (above 10 GeV) bands, which was measured with the procedure described in Pacciani (2018). A historically high X-ray flux triggered the observations in December 2019 (2.7 h after cuts) and in December 2020 (13.5 h after cuts). The data were analysed using the MAGIC Analysis and Reconstruction Software (MARS; Zanin et al. 2013; Aleksić et al. 2016).

The data analysis was performed by separating the data into four different epochs: the first corresponds to the data taken from 2009 to 2011 (MJD 55131–55620, hereafter epoch E1); while the second (MJD 56986–56987, E2), third (MJD 58819–58821, E3) and fourth (MJD 59198–59208, E4) correspond to the target of opportunity observations in 2014, 2019, and 2020, respectively. Table 1 shows the significances of the detection during the different epochs as estimated following Eq. (17) in Li & Ma (1983).

2.2. HE γ rays: Fermi-LAT

The GeV γ-ray emission from 1ES 0647+250 was characterised with the LAT on board the Fermi Gamma-ray Space Telescope. The Fermi-LAT data presented in this paper were analysed using the standard Fermi analysis software tools (version v11r07p00), and the P8R3_SOURCE_V2 response function. We used events from 0.3 − 300 GeV selected within a 10° radius region of interest (ROI) centred on 1ES 0647+250 and having a zenith distance below 100° to avoid contamination from the Earth’s limb. The usage of events above 0.3 GeV (instead of above 0.1 GeV) is advantageous for sources with hard γ-ray spectra (photon index < 2.0), especially if the source is weak and long integration times are needed for significant detections. The higher minimum energy somewhat reduces the detected number of photons from the source, but this effect is small for hard sources. On the other hand, the angular resolution (68% containment) improves from ∼5° to ∼2° when increasing the energy from 0.1 GeV to 0.3 GeV, which reduces the diffuse backgrounds (which are always softer than photon index 2), thus making the analysis less sensitive to possible contamination from non-accounted (transient) neighbouring sources, and reducing the systematic uncertainties. The diffuse Galactic and isotropic components were modelled with the files gll_iem_v07.fits and iso_P8R3_SOURCE_V2_v1.txt, respectively¹. All point sources in the fourth Fermi-LAT source catalogue (4FGL, Abdollahi et al. 2020) located in the 10° ROI and an additional surrounding 5° wide annulus were included in the model. In the unbinned likelihood fit, the normalisation and spectral parameters of all the sources were fixed to the 4FGL values, with the exception of the seven sources within the ROI identified as variable and with a detection significance larger than 10σ, where the normalisation parameters were allowed to vary. For the three objects located within an angular distance of 5° of 1ES 0647+250 (i.e. 4FGL J0650.6+2055, 4FGL J0653.7+2815, and 4FGL J0709.1+2241), the spectral parameters were also allowed to vary. The normalisation of the diffuse components (Galactic and isotropic) was also allowed to vary in the unbinned likelihood fits. In the 4FGL-DR3 (Abdollahi et al. 2020, 2022), which integrates over 12 years, the log-parabola (LogP) function is preferred to reproduce the spectrum over the power-law (PL) function with a significance of 4.1σ. However, owing to the much shorter timescales used in this study (shorter than 2 years), we decided to parameterise the γ-ray spectral shape of 1ES 0647+250 with a PL, where both the normalisation (flux) and the PL index were kept as free parameters.

Owing to the moderate sensitivity of Fermi-LAT for the detection of 1ES 0647+250 on day/week timescales (especially when the source is not flaring), we performed the unbinned likelihood analysis on consecutive 30-day time intervals (not centred on the MAGIC time) to determine the light curve in the energy band 0.3 − 300 GeV, as reported in Fig. 1. The source is detected with a maximum-likelihood test statistic (TS)² above 4 for most of the 30-day time intervals. There are seven 30-day time intervals that yielded a TS below 4, for which we computed 95% confidence level upper limits using a fixed PL index of 1.70, reported in the 4FGL-DR3 catalogue (Abdollahi et al. 2022). The PL spectral index of the source for each time bin is computed using the same procedure as the light curve analysis.

Fig. 1. MWL light curves of 1ES 0647+250. From top to bottom: MAGIC (> 100 GeV, 30-day binning), Fermi-LAT flux and spectral index (300 MeV–300 GeV, 30-day binning), Swift-XRT (0.3−2 keV and 2−10 keV), Swift-XRT spectral index, Swift-UVOT (B, V, U, UVW1, UVM2, and UVW2 filters), optical R band (KVA, Liverpool Telescope, Las Cumbres, and PIRATE telescope), and OVRO (15 GHz). We note that all the optical observations except those from 2020 were performed by KVA. The smaller error bars of the OVRO light curve after 2016 are due to a major upgrade of the instrument. Blue contours correspond to the MAGIC campaigns during the different observed states.

For the flaring episodes in 2014, 2019, and 2020, where we wanted to combine the data with VHE γ-ray spectra from MAGIC, we decided to use a time interval of 12 days for the first two, and 10 days for the third. We also computed the spectrum contemporaneous to the 2-year-long MWL observations in 2009–2011. The spectral results are reported in Sects. 4 and 6.

2.3. X-ray observations: Swift-XRT

The X-Ray Telescope (XRT; Burrows et al. 2004) on the Neil Gehrels Swift Observatory carried out 70 distinct observations of this blazar between May 2010 and December 2020. In particular, the source was observed several times distributed in the different observing campaigns previously defined. Swift-XRT pointed to 1ES 0647+250 a total of 25 times during E1, between May 2010 and March 2011 (MJD 55322–55623). Moreover, the source was also observed during the different flaring states in E2, E3 and E4. It was observed five times in November 2014 (E2, MJD 56981–56987). Another 19 observations were performed during and after the enhanced activity of 2019 (E3), from December 2019 until March 2020 (MJD 58816–58914). Finally, it was targeted eight more times in December 2020 (E4, MJD 59196–59207). An additional 13 more observations were performed, non-simultaneously to those from MAGIC. The Swift-XRT observations were carried out in the windowed timing and photon counting readout modes. The data were processed using the XRTDAS software package (v.3.6.0), which was developed by the Agenzia Spaziale Italiana (ASI) Space Science Data Center (SSDC) and released by HEASARC in the HEASoft package (v.6.28). The data were calibrated and cleaned with standard filtering criteria using the xrtpipeline task and the calibration files available from the Swift-XRT CALDB (version 20200724). For the spectral analysis, events were selected within a circle of 20-pixel (∼46 arcsec) radius, which encloses about 90% of the point spread function, centred at the source position. The background was estimated from nearby circular regions with a radius of 40 pixels. The ancillary response files were generated using the xrtmkarf task applying corrections for point spread function losses and Charge-Coupled Device (CCD) defects using the cumulative exposure maps. Before the spectral fitting, the 0.3−10 keV source energy spectra were binned using the grppha FTOOL to ensure a minimum of 20 counts per bin.

2.4. UV/optical observations: Swift-UVOT

The Ultra-Violet and Optical Telescope (UVOT; Roming et al. 2005) on board the Swift satellite, Swift-UVOT, has performed photometric observations in three optical (U, B and V) and three UV (UVW1, UVM2, and UVW2) filters, for a total number of 70 observations from May 2010 to December 2020. All the UVOT observations are simultaneous to those performed by XRT.

We evaluated aperture photometry for each total exposure applying the official software included in HEASoft package (v6.23), with a final check for attitude stability. We extracted the source counts within the standard circular aperture of 5″ radius, and the background counts from an annular region of inner radius 26″ and 9″ size. We applied the official calibrations (Poole et al. 2008; Breeveld et al. 2011) from the Swift-UVOT CALDB (version 20201026) to convert source counts to fluxes and then, a mean Galactic E(B − V) value of 0.0835 mag (Schlafly & Finkbeiner 2011) and an interstellar extinction curve (Fitzpatrick 1999) were used to obtain νF(ν) values at filter effective frequencies.

2.5. Optical data

Optical monitoring of the source in the R band was also performed by the Tuorla blazar monitoring programme³. For these observations, the 35 cm Kungliga Vetenskapsakademien (KVA) telescope, located in La Palma, was used. The data analysis was performed following the procedure described in Nilsson et al. (2018). This analysis includes the subtraction of the stellar emission from the host galaxy and the correction for Galactic extinction. The monitoring of this source started in December 2002, and continued until December 2019 (MJD 52615–58835).

Optical observations of 1ES 0647+250 were also performed in December 2020 with the 0.4 m robotic telescopes of Las Cumbres Observatory (LCOGT; Brown et al. 2013). During this period, the source was also observed by the robotic 2.0 m Liverpool Telescope (LT) at the Roque de los Muchachos Observatory in La Palma (Steele et al. 2004). These observations were performed with the Infrared-Optical (IO) instrument and its optical imaging component, the IO:O. Furthermore, it was observed during the night of 22 December 2020 by the 43 cm PIRATE (Physics Innovation Robotic Astronomical Telescope Explorer) telescope located at the Teide Observatory, on the Canary island of Tenerife (Holmes et al. 2011).

2.6. Radio observations: OVRO

1ES 0647+250 is also part of the Owens Valley Radio Observatory (OVRO) blazar monitoring programme (Richards et al. 2011)⁴. These observations were conducted with the OVRO 40 m radio telescope, working at a frequency of 15 GHz. The source was monitored by OVRO from January 2008 until December 2020, covering all the MAGIC observing periods (MJD 54476–59199). The data reduction was performed according to the procedure described in Richards et al. (2011). Observations with a signal-to-noise ratio < 3 were treated as non-detections and thus were not included in the MWL light curve and analysis. This resulted in 446 observations after excluding these measurements from the analysis.

3. Multi-wavelength light curve analysis

The MWL light curves of 1ES 0647+250, from VHE γ rays to radio wavelengths, are presented in Fig. 1. All the curves are daily binned except for the MAGIC and Fermi-LAT light curves, for which 30-day binning is used due to the limited ability of these two instruments to detect 1ES 0647+250 in the HE and VHE γ-ray bands when the source is not flaring. In the following subsections, the variability and interband correlations of the MWL data set are evaluated. A description of the light curves, with maximum, mean and minimum flux values for each band, is included in the Appendix A.

3.1. Variability

Emission from blazars is known to be variable across the electromagnetic spectrum. We performed a variability analysis, that is, we searched for significant flux variations and patterns in the data at different timescales on the light curves presented in Fig. 1 by testing the steady-flux hypothesis in the different bands.

The sparse time coverage and the large number of times that the source was not detected significantly with MAGIC prevent one from performing a reliable variability analysis in the VHE γ-ray energy range. In the HE band, this source is catalogued as variable in the 4FGL-DR3 Fermi-LAT catalogue (Abdollahi et al. 2020, 2022), with a variability index of ∼315, estimated as defined in Table 6 of Abdollahi et al. (2022). A fit to a constant function to the Fermi-LAT fluxes reported in Fig. 1 yields a χ²/d.o.f. = 346/143 ≃ 2.4 (p_value = 10⁻²⁰). Therefore, the HE emission from this source is clearly variable on timescales of 30 days, showing an increasing trend in the flux over time. This long-term trend was evaluated by a linear fit with an increasing flux, resulting in a χ²/d.o.f. = 191/142 ≃ 1.3 (p_value = 0.004), which is preferred to a fit to a constant average flux. Due to the low HE flux level and limited sensitivity of the LAT, we cannot detect possible variations in timescales shorter than one month.

In radio to X-rays the source is significantly variable on long-term timescales of the order of several months or years, as is typical for blazars and in line with previous studies for this source (see e.g. Kapanadze 2009, 2019, 2020; Kiehlmann et al. 2014; Nilsson et al. 2018). The results of the variability analysis and the goodness of the constant fit performed for each wave band are shown in Table 2. Moreover, the optical R-band and radio light curves also show the same increasing trend observed in the HE γ-ray light curve. An increasing linear fit is able to describe this steady flux increase over the years. However, since these bands also show variability on shorter timescales, as observed in Fig. 1, the χ²/d.o.f. of both fits is still > 1. Additionally, we also investigated the variability of the data from E1 and E4, as they are the periods with best MWL coverage. The results are also displayed in Table 2. Significant variability in the X-ray band was found during these two epochs, as well as variability in the optical R-band for E1.

Moreover, it is also important to quantify the amount of variability displayed in each band. This can provide useful information about the dynamics of the particle population that dominates the emission in the energy band probed. For this purpose, we used the fractional variability as a measurement of the degree of variability. We followed the prescription of Vaughan et al. (2003), where this parameter is estimated as

$$\begin{array}{c}\hfill F_{\mathrm{var}}=\sqrt{\frac{S^2-⟨{\sigma }_{\mathrm{err}}^2⟩}{{⟨x⟩}^2}},\end{array}$$(1)

where ⟨x⟩ and S² are the mean and the variance of the distribution of measured fluxes, respectively, and $\langle {\sigma }_{\text{err}}^2\rangle$ is the mean square error of the data. The uncertainty associated with the fractional variability is estimated following the prescription of Poutanen et al. (2008), as described by Aleksić et al. (2015a):

$$\begin{array}{c}\hfill \mathrm{\Delta }{F}_{\mathrm{var}}=\sqrt{F_{\mathrm{var}}^2-\mathrm{err}({\sigma }_{\mathrm{NXS}}^2)}-F_{\mathrm{var}},\end{array}$$(2)

where err$({\sigma }_{\text{NXS}}^2)$ is the normalised excess variance taken from Vaughan et al. (2003), calculated as

$$\begin{array}{c}\hfill \mathrm{err}({\sigma }_{\mathrm{NXS}}^2)=\sqrt{{\left(\sqrt{\frac{2}{N}}\frac{⟨{\sigma }_{\mathrm{err}}^2⟩}{{⟨x⟩}^2}\right)}^2+{\left(\sqrt{\frac{⟨{\sigma }_{\mathrm{err}}^2⟩}{N}}\frac{2F_{\mathrm{var}}}{⟨x⟩}\right)}^2},\end{array}$$(3)

where N is the number of data points. A deeper discussion on the estimation and caveats of the fractional variability is presented in Aleksić et al. (2015a) and Schleicher et al. (2019), and references therein. The results of the fractional variability for 1ES 0647+250 are presented in Fig. 2.

Fig. 2. Fractional variability of 1ES 0647+250. Filled markers represent MWL observations quasi-simultaneous to the MAGIC observations. Open markers correspond to the Fvar of the complete data sets. Filled markers represent the Fvar using only the simultaneous MWL data with respect to the MAGIC observations.

Owing to the remarkably different temporal coverage of 1ES 0647+250 at 15 GHz, R band, and HE γ rays (where the data span over the entire multi-year period considered in this study), in comparison to the UV, X-rays and VHE γ rays, we decided to apply two strategies to quantify the fractional variability. On the one hand, we used all flux values reported in the light curves from Fig. 1 to compute the variability at radio, R-band, UV, X-ray, and HE γ-ray energies. The results are displayed with open markers in Fig. 2, and show a slight increase in the overall flux variability with increasing energy. Additionally, we computed the fractional variability for all the energy bands sampled, but this time selecting only flux measurements that relate to observations that were performed quasi-simultaneous to those from MAGIC (±0.5 days). In the case of Fermi-LAT, where the flux measurements relate to 30-day bins, we use the GeV flux from the 30-day bin that contains the MAGIC observations. The results obtained with this strategy are displayed with filled markers in Fig. 2. The highest variability occurs at VHE γ-ray energies, although the statistical errors are large due to the relatively large flux measurement errors and the somewhat limited data set collected with MAGIC, biased towards bright flares. On the other hand, the variability at X-rays has small uncertainties, because of the smaller flux measurement errors and the larger data set, and it is clearly larger than the variability at radio, optical/UV, and HE γ rays. Finally, for comparison, we also estimated the fractional variability of the radio/optical/UV/X-ray data sets with a 30-day binning matching the one from the HE γ-ray data, with no significant differences with respect to the values shown in Fig. 2. This indicates that the short-scale variations are less important and the dominant variability corresponds to the long-term variability.

1ES 0647+250 shows lower variability in radio and optical wavelengths than that at X-rays and HEs. Its fractional variability has a double-maximum shape like, for instance, Mrk 421 (Aleksić et al. 2015a; Baloković et al. 2016; MAGIC Collaboration 2021), MAGIC 2001+439 (Aleksić et al. 2014) and 1ES 1959+650 (Kapanadze et al. 2018). The variability increases from its minimum at radio and optical frequencies, reaching a maximum at X-ray wavelengths, followed by a drop at HE γ rays and an increase in the VHE γ-ray regime. This behaviour differs from other sources such as Mrk 501 (Ahnen et al. 2017, 2018) or TXS 0506+056 (Acciari et al. 2022), whose fractional variability progressively increases with the frequency, displaying its maximum variability in the VHE γ-ray domain. This may be a sign of different particle populations, environments and/or processes in the jet; and a higher synchrotron dominance in the jet (Aleksić et al. 2015a). We note, however, that for 1ES 0647+250, the F_var may be biased in γ-ray energies due to the 30-day binning of the Fermi-LAT and MAGIC data. Also, the γ-ray F_var value shows very large statistical uncertainties due to the poor sampling and hence, it is not conclusive. Moreover, the structure of the F_var plot may change with time, indicating that the population or the processes in the jet may also change (see e.g. Furniss et al. 2015, where a double-maximum structure is also reported for Mrk 501). The discrepancy between the R-band observations and the UV and optical Swift filters is well understood and caused by the lower coverage of the filters. The effect of this coverage difference is especially noticeable between the optical and UV Swift filters, causing the minimum F_var to occur in the former band. However, when calculating the F_var using only simultaneous optical and UV observations, one obtains the same F_var value (within statistical uncertainties). Therefore, this indicates that the low F_var in the optical filters is mostly due to the limited time coverage, and that, for equal time coverage, the F_var is slightly larger at optical and UV than in the radio band.

3.2. Correlation

We carried out a correlation analysis between the light curves of the different wave bands available. For this, we made use of the Z-transform discrete correlation function (ZDCF; Alexander 2013). This tool is a modification of the classical discrete correlation function (DCF; Edelson & Krolik 1988) with better performance under uneven sampling conditions. To estimate the significance of the cross-correlations, we followed the procedure described in detail in Max-Moerbeck et al. (2014a). We simulated 3×10⁵ artificial light curves with the same sampling and power spectral density as our data set, as described in Emmanoulopoulos et al. (2013)⁵. The power spectral density slopes derived for each wave band are α_radio = 1.26 ± 0.18, α_R = 1.62 ± 0.13 and α_γ rays = 0.75 ± 0.41, assuming a PL shape. This method has been widely used in the past to compute the significance of auto-correlation and cross-correlation studies (see e.g. Max-Moerbeck et al. 2014a; Lindfors et al. 2016; Otero-Santos et al. 2020). To perform this analysis, we used the radio, optical and HE γ-ray light curves from 2008 to 2019 to avoid introducing the gap present in the optical light curve between 2019 and 2020. We performed the correlation analysis for each pair of light curves.

The results of the correlation analysis are shown in Fig. 3. We found a maximum positive correlation between the optical R-band and the HE γ-ray light curves of r = 0.60 with a significance of ∼3σ with respect to the no-correlation hypothesis, for a time lag of −17 days. Moreover, a long-term correlation (r = 0.67, ∼3σ significance) was also found between the radio and optical light curves, with its maximum degree of correlation observed at a delay of −398 days. Finally, a 4σ long-term correlation (r = 0.50) was also observed between the radio and HE γ-ray wavelengths, with the maximum correlation found at a similar delay as the radio-optical pair, −393 days, meaning that the radio emission is delayed with respect to the optical and γ-ray bands. However, the relatively wide range of time lags for which the correlation remains highly significant (which includes a 3σ correlation at time lag zero for the light curves at radio and γ rays), indicate that the correlation is dominated by the long-term trend, and excludes that the correlation occurs only for a specific time lag. In any case, one can estimate the most representative time lag and its related uncertainty using various strategies. When using the formalism described in Alexander (2013), one obtains that largest ZDCF (ZDCF_max) between the optical R-band and the HE γ-ray light curves occurs at a time lag of 17 ± 30 days; between the radio and optical light curves, the maximum degree of correlation occurs at a delay of −398 ± 80 days; and between the radio and HE γ-ray wavelengths at a delay of −393 ± 40 days. Additionally, we also used the model-independent Monte Carlo flux randomisation and random subset selection method described in Peterson et al. (1998, 2004), obtaining that the centroid of the DCF for correlations above 2σ (DCF_cen) and the related 68% confident limit uncertainties for the optical and γ-ray light curves is −7 ± 105 days; for the radio and optical light curves −380 ± 88 days; and for the radio and γ-ray light curves −332 ± 144 days. The 95% confidence limit uncertainties (using Peterson et al. 2004) are ±174 days and ±219 days for the last two cases. Therefore, we can confirm that the highest degree of correlation between radio and the optical and γ-ray light curves occurs with a time lag: the radio emission is lagging the optical and the γ-ray emission. We note that for light curves with a large time coverage and these dominant long-term trends, the effects of short-term correlations are generally masked by the variations in long timescales (Smith et al. 1993; Lindfors et al. 2016; Raiteri et al. 2021). Thus, these correlations refer to the aforementioned long-term trend that can be seen in the MWL light curves, where the flux increases over the years in radio, HE γ rays, and especially in optical.

Fig. 3. Long-term cross-correlation curves. The ZDCF is represented in black and its 1σ uncertainty by the grey contour. Coloured dotted lines correspond to different significance levels, from 1σ to 4σ. Left: cross-correlation between the optical and HE γ-ray light curves. Middle: cross-correlation between the radio and optical light curves. Right: cross-correlation between the radio and HE γ-ray light curves.

We also investigate whether the light curves showed any correlation at shorter timescales of the order of weeks or months. For this purpose, we performed a detrending of the long-term flux increase for all three of the bands. We followed the procedure described in Lindfors et al. (2016) and MAGIC Collaboration (2020), detrending the light curves in pairs (radio-optical, radio-HE γ rays, and optical-HE γ rays).

First, we fitted the lower-frequency light curve of each pair (e.g. radio light curve for the radio-optical pair) to a polynomial function. The polynomial order was determined by adding orders until the fit that minimises the χ²/d.o.f. value was found. This function describes the long-term variation of the light curve.

The polynomial fit was then scaled so that its variance equalled that of the high-frequency light curve, and the average flux of this curve was added. Then, the polynomial was multiplied by 0.1, 0.2, ..., 1.0 and subtracted from the high-frequency light curve. The best subtraction was determined by calculating the factor that minimises the root mean squared (rms)⁶ of the subtracted light curve. The result of this subtraction is a light curve where the common long-term variation of both the low and high frequency data sets is removed.

Third, the fractional contribution of the subtracted long-term slowly varying component was estimated by dividing the rms of the original data set and the rms of the light curve obtained after subtracting the polynomial function:

$$\begin{array}{c}\hfill \mathrm{Fraction}=1-{\mathrm{rms}}_{\mathrm{subtracted}\mathrm{light}\mathrm{curve}}/{\mathrm{rms}}_{\mathrm{original}\mathrm{light}\mathrm{curve}}.\end{array}$$(4)

This procedure quantifies the common slowly variable component between two wave bands. The detrended light curves can be seen in Fig. 4. Moreover, the fitted trends for the radio and optical data sets are included in the Appendix B (see Fig. B.1). We find that the slow varying component between the radio and optical light curves accounts for a fraction of 0.53 of the total variability. This value is much higher than the one reported by Lindfors et al. (2016) for this source. They found that common radio-optical component has a contribution of 0.1. However, the optical light curve used for their analysis corresponds only to the first half of the data set presented in this work. We estimated the fractional contribution for both halves of our optical curve, reproducing the result presented in Lindfors et al. (2016) for the data between 2008 and 2013, while we obtain a value of 0.4 for the second half of the data. This result explains the different values reported by both analyses. As for the OVRO-LAT light curve pair, this component explains a fraction of 0.23 of the flux variation. Finally, the slow component for the optical and Fermi-LAT light curve pair is responsible for a fraction of 0.24 of the total variability.

Fig. 4. Light curve detrending results. Blue dashed lines represent the scaled low-frequency light curve best fit, red dots correspond to the original data set and black dots show the detrended light curve. Left: Fermi-LAT light curve after optical trend subtraction. Middle: optical light curve after radio trend subtraction. Right: Fermi-LAT light curve after radio trend subtraction.

There may still be correlated emission on shorter timescales. In order to search for these short-term correlations, we applied the ZDCF to the detrended light curves. We do not find any correlated emission in these short timescales between radio, optical and HE γ-ray bands, with correlation coefficients < 0.2. In the case of the HE γ-ray band, we do not detect significant long-term variability after detrending, with a χ²/d.o.f. = 177/143 ≃ 1.2 (p_value = 0.028). Moreover, the low coverage of the MWL data and large time bins of the Fermi-LAT light curve due to the low flux of the source do not allow us to perform a detailed short-term correlation analysis of each individual observing epoch.

4. Spectral analysis

4.1. X-ray spectral analysis

An analysis of the X-ray spectra collected by Swift-XRT was carried out. We compared the X-ray spectral behaviour observed in different time intervals during which the target was also observed by MAGIC. Table 3 displays the spectral parameters and flux states during the different periods. We note that during the low state in epoch E1 and the flare from E4, the source showed significant variability in the X-ray spectrum. To account for this variability, we report the spectral parameters for the maximum, mean and minimum flux states. For the enhanced states from E2 and E3, the closest spectra in time to the MAGIC detections are shown.

The different spectra were fitted with a simple PL function (dN/dE = f₀ ⋅ (E/E₀)^−α). A LogP fit was also tested (dN/dE = f₀ ⋅  (E/E₀)^{−α − β  ⋅  log(E/E₀)}). Both models also included a photoelectric absorption with a neutral-hydrogen column density fixed to the Galactic value in the direction of 1ES 0647+250, namely 1.20 × 10²¹ cm⁻² (HI4PI Collaboration 2016). However, there is no statistical preference for the LogP over the PL, with χ²/d.o.f. ∼ 1.0 − 1.5 for both models. Thus, for this analysis we assume the simpler spectral shape defined by the PL. The spectra from E1 show a variation of the spectral index of ∼0.4, varying from values of α_XRT ∼ 2.4 when the source is fainter, up to much harder values of α_XRT ∼ 2.0 when the source is in a higher state (see Fig. 5). The correlation between the spectral index and the X-ray flux was quantified by estimating the Pearson’s linear correlation coefficient for the integral X-ray flux in the 0.3−10 keV band, obtaining a value of r = −0.65 ± 0.15 for this campaign, with a p-value of 3 × 10⁻⁴. This behaviour of harder-when-brighter has been seen in the past for other blazars such as 1ES 2344+514 (Acciari et al. 2011a; Aleksić et al. 2013), Mrk 421 (MAGIC Collaboration 2021) or TXS 1515−273 during bright X-ray flares (Acciari et al. 2021a). This has found to be a typical behaviour for blazars in the X-ray domain (see e.g. Wang et al. 2018).

Fig. 5. PL index αXRT vs X-ray flux for the different epochs. Left: 0.3−2 keV flux. Right: 2−10 keV flux.

The spectral analysis of the 2014 flare (E2) reveals a steeper spectrum compared to those from the non-flaring state. During this period, the X-ray spectral index varies around α_XRT ∼ 2.5, as can be seen from Fig. 5. We also note that we do not see this harder-when-brighter behaviour during this epoch. However, this could also be due to the sparse time coverage during this period, with only four X-ray observations.

The enhanced state observed in E3 shows the highest X-ray flux ever detected for this source (Kapanadze 2019). The spectral index during this flare displays values of α_XRT ∼ 2.4, reaching a harder index of α_XRT ∼ 2.2 during the maximum of the detection. No clear trend is seen between the spectral index and the X-ray flux for this period. However, a visual inspection of the E3 data set (black points in Fig. 5) reveals a different behaviour for the faintest measurements (those with a flux F < 8 × 10⁻¹¹ erg cm⁻² s⁻¹ for the band between 0.3 keV and 2 keV, and F < 4 × 10⁻¹¹ erg cm⁻² s⁻¹ for the 2−10 keV X-ray band, performed after the VHE flare) than for the brightest measurements (those performed during the historically high X-ray activity and roughly simultaneous to the observations performed by MAGIC). The former subset is characterised by a Pearson’s correlation coefficient between the spectral index and the flux for this data set in the 0.3−10 keV band of r = −0.64 ± 0.20 and a p-value of 8 × 10⁻³. The latter subset, however, shows no correlation between the flux and the spectral index. This result may indicate a saturation of the X-ray spectral index for the highest X-ray fluxes during the flare observed by Swift-XRT.

Finally, during the 2020 observing period (E4) the source displayed significant X-ray variability. This variability was also observed during E1, where the source showed variability both in its X-ray flux (as reported in Table 2) and spectral index, with the harder-when-brighter behaviour already reported. The spectral index during E4 ranges from α_XRT ∼ 2.2 to α_XRT ∼ 2.4. However, contrary to the results from E1, the X-ray observations performed during this period do not reveal any correlated evolution of the spectral index and the flux.

When the integrated flux in the 0.3−10 keV band was considered, the same results as those presented above were obtained. This suggests that the spectral index does not vary between the 0.3−2 keV and the 2−10 keV X-ray bands. In summary, two of the epochs (E2 and E4) do not reveal any clear behaviour in their X-ray spectra. On the other hand, a harder-when-brighter behaviour was detected for E1 and E3, in the latter followed by an index saturation for those observations simultaneous to the brightest X-ray observations, and closest to those from MAGIC. This saturation has also been observed for other blazars like Mrk 421 in the past (Acciari et al. 2021b). Moreover, the harder-when-brighter trends have been explained in terms of, for instance, a change in the maximum energy of the electrons responsible for the emission (see e.g. Abeysekara et al. 2017) or a hardening of the electron distribution (Xue et al. 2006). A detailed discussion on the spectral variability studied here can be found in Sect. 7. We also searched for hysteresis loops in the X-ray spectral index evolution during the X-ray flares. However, the low coverage and uncertainties in the spectral index characterisation do not allow us to observe any clear hysteresis-like behaviour.

4.2. HE γ-ray spectral analysis

The HE γ-ray spectrum and SED of 1ES 0647+250 were extracted for all the analysis epochs making use of the Fermi-LAT data. A PL shape was assumed. This fit was performed in the energy range of (0.3−300 GeV). The low energy photons (< 3−5 GeV, below the IC peak) dominate the fit of the spectrum, leading to the observed ‘hard’ indices. Table 4 shows the spectral parameters of the HE band for the different epochs. No significant variability in the slope of the spectrum was detected in the Fermi-LAT observations. Moreover, contrary to the results for the X-ray spectrum, here we do not see any correlation of the spectral index with the flux. The spectral index light curve is compatible with a constant value, α_Fermi = 1.70 ± 0.02 and a value of χ²/d.o.f. = 126.8/143 ≃ 0.9, which is similar to α_Fermi = 1.73 ± 0.02 reported in the 4FGL-DR3 catalogue (Abdollahi et al. 2020, 2022).

For the monitoring performed in E1, all the Fermi-LAT data from November 2009 to March 2011 were used. For the enhanced states of 2014 (E2) and 2019 (E3), a 12-day integration window centred around the MAGIC detection was used due to the low HE γ-ray flux displayed by 1ES 0647+250 in shorter timescales during these epochs, leading to a TS < 25 and mostly upper limits in the spectral points. Finally, for the data from 2020 (E4), Fermi-LAT data simultaneous to the MAGIC observations were used to calculate the spectrum. The spectral parameters of each period are shown in Table 4.

4.3. VHE γ-ray spectral analysis

The spectrum and SED were also obtained for the VHE γ-ray band for the different time periods of the analysis. The spectra from the observations from E1, E2, and E3 were well modelled with PLs. In contrast, a 3σ preference for a log-parabolic shape was observed in the spectrum of E4. The results of the MAGIC spectral analysis are summarised in Table 5. For the rest of the periods, the statistics are not high enough to evaluate a log-parabolic spectral shape. We note that, since the redshift of the source is still under debate, the spectra and SEDs presented in this analysis correspond to the observed spectra, without the correction for extragalactic background light (EBL) absorption.

5. Redshift estimation

We used the joint Fermi-LAT and MAGIC spectra to constrain the redshift of 1ES 0647+250. To perform this estimation, we followed the procedure proposed by Prandini et al. (2010a). This method is based on the assumption that the VHE γ-ray spectrum of a blazar after correcting for the EBL absorption, cannot be harder than the spectrum in the HE range measured by Fermi-LAT. It makes use of the EBL model developed by Franceschini et al. (2008). First, an upper limit of the redshift, z^*, is calculated as the limit value at which the slopes of the HE and VHE γ-ray spectra are equal by de-absorbing the MAGIC spectral points until the spectral indices of the HE and VHE spectra are equal. Then, the empirical formula that relates z^* and the reconstructed value of the redshift, z_rec, with the updated parameters presented in Prandini et al. (2010b), is used to estimate the reconstructed value. We performed this estimation with the Fermi-LAT and MAGIC SED from E4 (reported in Tables 4 and 5) due to their higher flux and thus, higher statistics and lower uncertainties in the flux and spectral index estimation.

This empirical relation applied to our data led to a redshift estimation of z_rec = 0.45 ± 0.05. This derived value of the redshift is in agreement with the current lower limit of z > 0.29 estimated by Paiano et al. (2017) from spectroscopic observations, and the most reliable measurement of the distance by Kotilainen et al. (2011), who reported a value of z = 0.41 ± 0.06 from the detection of the host galaxy. Moreover, the maximum redshift value obtained with this method is z^* = 0.75 ± 0.11, compatible with the estimations of the distance of this blazar cited above.

This method has proven to report accurate values of the distance of several blazars in the past, for instance, MAGIC J2001+435 (Aleksić et al. 2014) or S5 0716+714 (MAGIC Collaboration 2018). However, it has also reported some inconsistent values, for instance PKS 0447−439, with an estimation of z = 0.20 by Prandini et al. (2012), later measured to be z = 0.343 by Muriel et al. (2015). We note that this empirical procedure has different caveats and assumptions when estimating the redshift of a source. A detailed discussion of these caveats can be found in Prandini et al. (2010a,b), where the redshift of several γ-ray emitting blazars was properly estimated.

Alternatively, a second upper limit was derived using a maximum likelihood fit of the joint Fermi-LAT and the MAGIC SEDs from E4, using a concave LogP as the spectral model, as described in Acciari et al. (2019). Using the EBL model from Domínguez et al. (2011), a scan of redshifts was performed, obtaining a likelihood profile from which the redshift can be constrained. A 15% systematic uncertainty in the overall light throughput was taken into account following the studies in Aleksić et al. (2016). Under these considerations, the 95% confidence level upper limit for the redshift z is 0.81. In the following SED modelling we assume the value of z = 0.41 from Kotilainen et al. (2011), which is in agreement with the value we derived from the empirical relation.

6. Broadband SED

The broadband emission of blazars has been successfully described in the past with models based on leptonic emission processes. However, due to the still arguable origin of the high-energy SED peak, hadronic models have also been found successful on several occasions, especially in the scenario of neutrino emission (see for instance Petropoulou et al. 2017, 2020a,b; Kreter et al. 2020). Each model has its strengths and limitations, with the hadronic models being favoured by uncorrelated optical, X-ray and γ-ray variability (see e.g. Dimitrakoudis et al. 2012; Böttcher et al. 2013). Considering the correlated MWL variability observed for this source, we propose here an interpretation based on leptonic models. While one-component models provide an easier solution due to the smaller number of free parameters, they have not always been adequate to reproduce the emission of γ-ray blazars with respect to two-component models (MAGIC Collaboration 2020).

We present the modelling with both one- and two-component leptonic models, comparing their performance and capability to reproduce the observed features. The modelling was performed assuming the cosmological parameters reported by Planck Collaboration VI (2020): a Hubble parameter of H₀ = 67.4 km s⁻¹ Mpc⁻¹, a matter density of Ω_m = 0.315 and a dark energy density Ω_Λ = 0.685. The models used here are not time-dependent. Hence, the different epochs are modelled independently. Given the sparse data sets and the large time separations between the different SEDs, no firm conclusions can be drawn on the temporal evolution of the model parameters.

6.1. One-component model

As an initial approach, a one-component SSC model is used to reproduce the broadband SEDs of 1ES 0647+250 (Tavecchio et al. 1998). This model assumes the existence of a single, homogeneous and spherical emitting region in the jet with size, R, Lorentz factor, Γ, and magnetic field, B. The low-energy bump of the SED is due to synchrotron radiation, while the high-energy bump is modelled through SSC. The population of electrons inside the emitting region is assumed to follow a broken PL distribution with the Lorentz factor, described by

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

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).

The parameters of the one-component models for each epoch are reported in Table 6. The resulting models during each VHE γ-ray detection are shown in Fig. 6. We scanned several combinations of the parameters described above in order to perform the modelling. The agreement between the model and the data is evaluated through visual inspection of the SEDs shown in Fig. 6. The one-component SSC model is able to describe well the observational data from optical to VHE γ rays. However, the radio data are strongly self-absorbed and thus, the emission at radio wavelengths is assumed to originate from a different region (Tavecchio et al. 1998). Therefore, this model is not able to reproduce the radio emission.

Fig. 6. Broadband SEDs of 1ES 0647+250 for the different observing campaigns described with one-component and two-component SSC emission models. Top left: 2009–2011. Bottom left: 2014. Top right: 2019. Bottom right: 2020. Filled blue dots correspond to the MWL simultaneous data of each period. Empty dots represent the archival data extracted from the SSDC database (https://www.ssdc.asi.it). The dotted lines correspond to the fitted one-component SSC model. Black lines represent the total two-component SSC model. The dashed and dotted-dashed lines correspond to the blob and core components of the two-component model, respectively.

6.2. Two-component model

Alternatively, we also model the SEDs in different epochs with a two-component model based on Tavecchio et al. (2011). This model calculates synchrotron and SSC emission for spherical emission regions while taking into account synchrotron-self absorption. The strength of the magnetic field is typically assumed to scale with the distance from the central engine as d⁻¹. If the two components are separate, the one responsible for the X-ray and VHE γ-ray emission is closer to the central engine than the one responsible for radio and optical emission. Therefore, the former needs to have a stronger magnetic field than the latter component, of the order of ∼1 G. Tavecchio & Ghisellini (2016) show that the magnetic field strengths tend to be significantly lower than the values required for equipartition values in one-component models. Moreover, in two-component models it is difficult to reproduce the observed SED with the magnetic field strength values of the order of 1 G. Re-connection layers and radial structures of magnetic fields across the jet are possible ways to invoke reduced local magnetic field strengths (see discussion in Nalewajko et al. 2014). Therefore, similar to the approach used in MAGIC Collaboration (2020), we assumed two co-spatial and interacting emission regions to mimic a simple spine-sheath model. The two spherical emission regions are called ‘core’ and ‘blob’, with sizes R_core > R_blob. The regions are filled with electrons distributed in Lorentz factor according to a smoothed broken PL (see Eq. (5)) and the physical quantities are expressed in the co-moving frame of each individual region. Each of the emission regions has a Lorentz factor, Γ, size, R, and magnetic field strength, B. The following constraints were employed to reduce the number of free parameters for this model.

First, the measured full-width-half-maximum values of the major axis of the very-long-baseline-interferometry core can be used to calculate the upper limit of the size of the core emission region. The measured full-width-half-maximum value of the major axis is 1.88″ (Piner & Edwards 2014), which corresponds to R_core ≤ 3.3 × 10¹⁹ cm, assuming a flat universe and z ≃ 0.41.

Second, the size of the blob can be constrained from the shortest variability timescales observed as R ≤ δcτ/(1 + z), where δ is the Doppler factor and τ the variability timescale observed.

Third, the bulk Lorentz factor of the core region is limited to 4, which is common for TeV blazars (Piner & Edwards 2018). The bulk Lorentz factor is then converted to Doppler factor assuming a jet viewing angle ∼1/Γ and thus δ ∼ Γ.

Finally, the magnetic field strength of the core can be estimated from the very-long-baseline-interferometry ‘core shift’ measurements (Pushkarev et al. 2012) or considering the cooling timescale of the electrons from the variability timescale in the X-ray band (Bhatta et al. 2018; Acciari et al. 2021a). Such observations are not available for the source. Therefore, we follow the same assumption employed in MAGIC Collaboration (2020; i.e. 0.1 ≤ B ≤ 0.4 G and similar for core and blob).

The parameters of both the core and blob for the different two-component models of 1ES 0647+250 are reported in Table 6, and the models are displayed in Fig. 6. We can only set constraints on the size of the blob through the X-ray variability during 2020, with variability on timescales as short as 1 day. Considering the minimum and maximum Lorentz factors derived for the blobs, this leads to a blob size R ≤ 3.1 × 10¹⁶ − 4.2 × 10¹⁶ cm. For previous epochs, the variability timescales are longer and thus, the blob size can also adopt higher values. As for the one-component model, the fit is evaluated by visual inspection after scanning the parameters describing the broadband emission. All the two-component models are able to satisfactorily reproduce the broadband emission of our source. The comparison and interpretation of both models with the rest of the results from this study are discussed in Sect. 7.

7. Discussion

This manuscript reports the first detailed study of the broadband emission from radio to VHE γ rays of the blazar 1ES 0647+250. The study uses a data set that spans from 2009 to 2020, which allows the multi-band variability and correlations to be evaluated over timescales of years. For this, along with the observations performed by the MAGIC telescopes, we make use of radio data from OVRO, several optical telescopes (KVA, LT, LCOGT, and PIRATE), observations performed by Swift and its instruments UVOT in the optical and UV regime, and XRT in the X-rays, and HE γ-ray data from Fermi-LAT. In the following subsections we discuss the main results obtained and the implications of the observations reported in previous sections.

7.1. Variability

The variability analysis carried out for this source reveals that its broadband emission is clearly variable during the monitored period. The maximum of F_var appears at X-ray and VHE γ-ray wavelengths, as shown in Fig. 2. The minimum timescale detected in these bands with significant variability (> 5σ; see Table 2) is 1 day in X-rays, only during the flaring state of E4. Moreover, when using a 30-day binning for all bands, no significant difference is observed in the structure of F_var. This is in line with the fact that the variability of this blazar is dominated by the long-term variations of the contribution from a common component, as derived by the correlation analysis.

The fact that F_var has its maximum in X-rays has been related in the past with a higher dominance of the synchrotron emission (Aleksić et al. 2015a). This structure of F_var can reveal fundamental differences related to the particle populations and processes producing the broadband emission in blazars. Aleksić et al. (2015a) find a similar F_var structure in Mrk 421 to the one shown by 1ES 0647+250. One can compare its behaviour with that reported for Mrk 501 by Ahnen et al. (2017) or Aleksić et al. (2015b), where F_var increases with the frequency, reaching its maximum at the VHE γ-ray band. In the framework of the typical one-component model, the X-ray emission is mainly generated by the high-energy electrons, contrary to the VHE γ-ray emission, which is due to a combination of low-energy and high-energy electrons in Thomson and Klein-Nishina regimes, respectively (Abdo et al. 2011). Thus, a higher F_var in the X-ray domain, as for the case of 1ES 0647+250, may be an indication of higher variability of the high-energy electron population, while the combined low- and high-energy electron distribution may dominate during the γ-ray flares, which leads to a high F_var in VHE γ rays.

7.2. Correlation and contribution of the two components

Several studies in the past have found long-term correlations between the optical and radio emission of different sets of blazars with the optical emission leading the radio counterpart by a few hundred days (see for instance Hufnagel & Bregman 1992; Tornikoski et al. 1994; Hanski et al. 2002; Ramakrishnan et al. 2016; Acciari et al. 2021b). Here, we detect a correlation between the radio and optical emission of 1ES 0647+250 with the optical contribution leading the radio one by 398 days, and between the radio and HE γ-ray emission, with the γ rays leading the radio by 393 days. These correlations are detected at a level of ≳3σ, and they are in line with the studies mentioned above. Additionally, a correlation between the optical and γ-ray emission was found with a time lag compatible with zero. Due to the long-term nature of the flux variability that leads to the correlations reported here, multi-year data sets are needed to further increase the statistical significance of these results.

Hufnagel & Bregman (1992) interpreted these correlations as physically related regions, however on different timescales. This has also been stated by other authors (see e.g. Zhang et al. 2017), and the time lags interpreted as different cooling times between the electrons responsible for the radio and optical emissions (Bai & Lee 2003). Another plausible explanation is that the radio emission comes from an outer region, but is triggered by the same physical mechanisms as the optical and γ-ray emission. Max-Moerbeck et al. (2014b) explains this behaviour as due to the different opacities for the radio and γ-ray wavelengths to become observable. Under this scenario, we can estimate the distance between the radio and γ-ray emitting regions using Eq. (1) from Max-Moerbeck et al. (2014b). For this estimation, we use the value of the bulk Lorentz factor Γ = 4 obtained from the core of the two-component SED modelling. The Doppler factor, δ, is then obtained from the Lorentz factor assuming the approximation of δ ∼ Γ for a viewing angle θ ∼ 1/Γ. The redshift value of z = 0.41 ± 0.06 from Kotilainen et al. (2011) was used. We obtain for a time lag of −393 ± 40 days a distance d = 3.6 ± 0.4 pc. We also estimated the distances for changes on Lorentz factors by a factor of 2 (Γ = 2 and Γ = 8) as a conservative comparison of the derived value of d for different values of Γ, obtaining distances of d = 0.8 ± 0.1 pc and d = 14.9 ± 1.7 pc, respectively. Under this interpretation, the radio and MWL emission would come from physically separated regions. This estimation can also be performed deriving the Doppler factor as δ = [Γ ⋅ (1 − β cos θ)]⁻¹ (see e.g. Liodakis et al. 2017) instead of the aforementioned small angle approximation. We include this estimation in the Appendix C as a comparison. Another plausible interpretation in the scenario of a one-component model is introduced by Tramacere et al. (2022), where the delay between the radio and MWL emissions would be caused by an adiabatic expansion of the emitting region during its propagation along the jet. This expansion leads to a shift of the synchrotron self-absorption frequency to values comparable to or lower than the frequency of the radio emission.

We also investigated the correlations at shorter timescales, finding no significant correlated emission between the radio, optical and HE γ-ray bands. For this, we performed a detrending of the light curves based on the assumption that the emission is due to a combination of a common and an independent component, following the prescription of Lindfors et al. (2016). We were able to estimate the contribution of this common component to the emission for each pair of light curves. For the radio-optical, radio-γ-ray and optical-γ-ray pairs the values found were 0.53, 0.23 and 0.24, respectively. These values are compatible with those reported by Lindfors et al. (2016) for this source.

The existence of a common component for the different bands, and the fact that the emission is uncorrelated after subtracting this contribution, may indicate that the long-term variations (timescales of years) are driven by the same mechanism and they come from the common emitting region. In this scenario, the two-component model would be favoured, as the emission would not come from physically separated emitting regions. On the other hand, the short-term variations (timescales of days to months) would be due to the components that do not have common emission. This would be in agreement with the results presented by Ramakrishnan et al. (2016), who state that long-term variations and strong flares are correlated between the radio and optical bands in blazars. Additionally, marginal evidence of correlation between the radio and γ-ray emission is also reported by these authors. Finally, Ramakrishnan et al. (2016) also report strong optical-γ-ray correlation for a set of blazars, compatible with the results found for 1ES 0647+250. We note however that this blazar sample is mainly composed of FSRQs, indicating that these correlations may be present in both low- and high-peaking blazars. Additionally, Liodakis et al. (2018, 2019) find significant long-term correlations between the radio, optical and γ-ray bands for more than 100 sources on a sample of 178 blazars, including several HBLs. The optical and γ-ray emission typically shows time lags compatible with zero, while the radio emission is usually delayed by a few hundred days. Moreover, they also find that variability on shorter timescales is not necessarily correlated, with the detection of a large fraction of optical and γ-ray orphan flares.

7.3. Spectral analysis results

The spectral analysis carried out on the X-ray data of 1ES 0647+250 reveals different behaviours for the different observing epochs. During the low state, a clear harder-when-brighter trend between the X-ray flux and the spectral index was observed. X-ray brightening in HBLs has been related in the past to spectral hardening rather than to an overall flux enhancement (Giommi et al. 2021), leading to the commonly observed harder-when-brighter evolution for BL Lac objects (e.g. Pian et al. 1998; Acciari et al. 2011a; Aleksić et al. 2013; Kapanadze et al. 2014; Baloković et al. 2016; Wang et al. 2018; MAGIC Collaboration 2021). Several conjectures have been proposed to explain such behaviour. A possible explanation could be a hardening or softening of the electron distribution responsible for the synchrotron emission, leading to a hardening (softening) trend with the increasing (decreasing) flux (Xue et al. 2006). Additionally, this effect could also be due to an increase in the maximum electron energy or a shift of the synchrotron peak towards higher frequencies (Abeysekara et al. 2017). However, this shift of the peak is not evident looking at the different broadband SEDs.

Moreover, no hardening or softening of the X-ray spectrum was observed during the flares registered. However, the faintest observations from E3, performed after the MAGIC VHE detection, display again the hardening of the spectrum with the flux increase. On the other hand, the brightest observations and closest to the flare show a rather constant index, which could indicate the presence of a spectral index saturation during this flare. This behaviour is also reported in Acciari et al. (2021b) for the nearby blazar Mrk 421, which shows a harder-when-brighter behaviour for a large range of X-ray fluxes, but it shows a sort of saturation for the very high (and very low) X-ray fluxes (Acciari et al. 2021b). In this case, the saturation is present in both the X-ray and VHE bands.

We also searched for hysteresis processes in the X-ray evolution of the observed flares. These processes have been detected in the past during flares for several blazars, both in X-rays and γ rays (see for instance Nandikotkur et al. 2007; Abeysekara et al. 2017; Wang et al. 2018). These processes provide unique information regarding the acceleration, cooling and energy loss timescales (Böttcher & Chiang 2002), or possible lags between hard and soft X-rays (Wang et al. 2018). Nonetheless, no evidence of hysteresis processes was found in the X-ray spectral analysis of 1ES 0647+250.

Concerning the γ-ray band, we do not detect any harder- or softer-when-brighter trend in the long-term variability of 1ES 0647+250. Despite the fact that the harder-when-brighter trend has been detected in the past also in the HE and VHE γ-ray domains for some blazars (see e.g. Acciari et al. 2011b; Abeysekara et al. 2017; MAGIC Collaboration 2021), it has not been systematically observed. Authors like Nandikotkur et al. (2007) or Abdo et al. (2010) have not found clear correlations between the hardness index and the γ-ray flux of HBLs. However, we note that the results here may be biased by the wide binning of the Fermi-LAT and MAGIC data.

7.4. Comparison of one-component and two-component models

A visual inspection reveals that both models are able to reproduce the MWL emission of this blazar. Both the one-component and two-component models explained the broadband SEDs with a magnetic field of B = 0.16 G, except for the one-component model from E3. This one needs a slightly higher magnetic field (B = 0.18 G) due to the high dominance of the synchrotron emission with respect to the SSC scattering after the highest X-ray flux detected for this source. In this case, a magnetic field of B = 0.16 G is unable to well reproduce the broadband emission, leading to large differences between the model and the data regardless of the values of the other parameters. These values are in agreement with the typical magnetic fields derived for BL Lacs through a one-component SED modelling by Ghisellini et al. (2010), with B ∼ 0.05 − 1 G. Moreover, the Lorentz factor of the one-component model is similar to or lower than that of the blob of the two-component, while the core shows a rather small Γ. Ghisellini et al. (2010) infer Γ ∼ 10 − 20 for a large sample of BL Lacs, compatible with those used in this work to reproduce the broadband emission. The Lorentz factors and the parameters describing the distribution of the electron population γ_b and γ_max are also compatible with those derived by these authors.

The one-component model reproduces the broadband emission from the optical regime up to VHE γ rays. However, it cannot explain the radio emission with its emitting region. This is in line with the results obtained through the correlation analysis, where the radio was found to be delayed with respect to the optical and γ-ray emission. Since the radio photons suffer from strong absorption in the inner regions of the jet, this emission is most likely generated in the outer part of the jet. Thus, the one-component model does not reproduce the radio emission since this emission is not co-spatial.

On the other hand, the two-component can naturally reproduce the emission of 1ES 0647+250 from radio to VHE γ rays with co-spatial blob and core regions. In this framework, the core provides seed photons for the IC scattering occurring in the blob. The core dominates the emission from radio to optical wavelengths, while the blob is the main contribution to the emission from X-rays to VHE γ rays. This is in line with the results found in MAGIC Collaboration (2020). This leads to higher γ_min values for the blob, as shown in Table 6.

The use of a two-component model that reproduces the radio is consistent with the long-term slow trend displayed by 1ES 0647+250. We note that the derived sizes would suggest variability timescales shorter than those obtained for the slowly varying component from the data. This means that the slow variability cannot be caused by core-size or acceleration and cooling processes, which are generally assumed to be the origin of the faster variability; rather, it traces, for example, injection or decay phases of the central engine.

7.5. Equipartition and jet power

We calculated the contributions of the magnetic and electron energy densities (see Table 6). For the one-component model, Tavecchio & Ghisellini (2016) show that, for BL Lacs, the magnetic energy density is typically one to two orders of magnitude lower than the electron energy density. Here, we found that, except for the model of the 2020 (E4) flare, the parameters suggest the existence of equipartition within a factor of a few. This is in line with the results presented by Nievas Rosillo et al. (2022), where one-component models are used to reproduce the broadband emission of a sample of HBLs and EHBLs. Tavecchio & Ghisellini (2016) also note that by using two-component models it is possible to reproduce the observed SEDs of BL Lacs assuming equipartition between both energy densities. However, other studies based on two-component models (e.g. Acciari et al. 2021a, for TXS 1515–273) show that, while these models can lead to solutions close to equipartition during low activity states (as reported by Tavecchio & Ghisellini 2016), this does not seem to always be possible, especially during flares.

In the case of two-component models we found that for the core component, the magnetic energy density is dominant over the electron energy density in all cases by an order of magnitude. This ratio for the blob components is similar to that from the one-component model. However, the ratio of the magnetic and electron energy densities is closer to equipartition. This is expected since the blob component is much denser and filled with more-energetic particles than the core. We calculated the equipartition value for the system of two interacting regions using Eqs. (5), (9), (16), and (A.1) from Tavecchio & Ghisellini (2016). The ratios of energy carried out by the magnetic field to that by electrons are 38, 0.3, 1.1, and 0.3 for the E1, E2, E3, and E4 data sets. This is in line with the VHE γ-ray state of the source during these epochs. During E1, the source was in a low state when comparing its X-ray and VHE γ-ray emission with other periods. Therefore, the contribution of energy carried by the electrons is far less than the one carried by the magnetic field. This is reflected in the equipartition values of the whole system as well. It is also evident that during this period, a larger amount of emission in the VHE γ-ray band is coming from the interaction between the components. Such a contribution is smaller for the other periods where the VHE γ-ray emission is dominated by the emission from the blob.

We also estimated the kinetic energy carried by protons, electrons and the magnetic field of the emitting region. We assumed the simple solution where there is one proton per injected electron, and we made use of Eqs. (1)–(3) of Celotti & Ghisellini (2008). The results are shown in Table 7, Cols. (3)–(5). For the one-component model, we find that L_e is about one order of magnitude higher than L_B, as expected from the results reported by Celotti & Ghisellini (2008) for a much larger population of blazars. This comes from the fact that for BL Lacs, the high-energy emission is of the same order as the synchrotron contribution, and SSC scattering is the only radiative process involved. For the case of the core in the two-component models, we found that the kinetic energy carried by protons, electrons and the magnetic field is higher than the one obtained in the one-component model.

Finally, we also estimated the power carried by the jet using Eq. (6) from Ghisellini et al. (2009). This estimation was made again under the assumption of one proton per electron and thus, P_e ∼ P_p. The estimated values are shown in Table 7 along with the total power carried by the jet. The results derived from the one-component model suggest that the power carried by electrons (and protons) tends to be higher than that of the magnetic field. In the context of SSC models, this is a signature characteristic of BL Lacs, where almost all the power of the jet is being used to produce the observed radiation. The values obtained here are also compatible with those estimated for a large sample of BL Lacs, as reported in Fig. 4 from Ghisellini et al. (2010). The values derived from the two-component model are also compatible with those reported by these authors. However, in this case, we observed a higher power carried by the magnetic field, in contrast with the lower value for the one-component model.

8. Conclusions

We present, for the first time, a detailed characterisation and interpretation of the broadband emission of the BL Lac object 1ES 0647+250 between the years 2009 and 2020. The main results of this study are summarised as follows:

• The MWL emission is clearly variable on long timescales, with an increasing flux in radio, optical, and γ-ray wavelengths. Such behaviour has been seen for other blazars (e.g. 1ES 1215+303; Valverde et al. 2020), where the flux increase over year timescales is compatible with that expected from variations in the conditions of the accretion disc. While this interpretation is applicable to the 11-year data set from 1ES 0647+250, more data would be necessary to fully confirm or rule out this hypothesis.

• A long-term correlation with no delay is measured between the optical and γ-ray emission at a confidence level of ∼3σ. Moreover, the radio is correlated with the optical and the γ-ray bands (at a statistical significance of 3σ and 4σ) with time lags of 393 ± 40 days and 398 ± 80 days, respectively. This time delay is compatible with the radio being emitted from a distinct region of the jet, at a distance of d = 3.6 ± 0.4 pc, assuming Γ = 4, δ ∼ Γ, and θ = 1/Γ.

• A harder-when-brighter behaviour in the X-ray spectra is observed during the low state as well as for the flare from 2019 (E3), followed in the latter case by a spectral index saturation. No clear relation is observed in the γ-ray domain.

• We estimate the redshift of this object through the comparison of its simultaneous GeV and TeV spectra during a flaring activity (using the method described in Prandini et al. 2010a), obtaining a value of z = 0.45 ± 0.05, which is in agreement with some tentative measurements reported in the literature.

• The broadband SED is characterised and interpreted within one- and two-component leptonic scenarios for four distinct epochs, namely the low activity in 2009−2011 (E1) and three flaring activities in the years 2014 (E2), 2019 (E3), and 2020 (E4). All the models use a magnetic field of about B = 0.16 G, and the energy loss of the electron population is dominated by synchrotron emission (needed to match the high X-ray flux and spectra collected during the different epochs). The model parameters required by the two theoretical scenarios to describe the broadband SEDs are similar to those from typical BL Lacs reported in Ghisellini et al. (2010) and Tavecchio & Ghisellini (2016) for a large sample of objects.

This paper describes a new and comprehensive MWL analysis of the GeV-emitting blazar 1ES 0647+250, one of the few distant gamma-ray blazars detected at VHEs, thoroughly characterising its long-term evolution over the last decade and evaluating its broadband SED for the first time.

Acknowledgments

We want to thank the anonymous referee for her/his useful comments and discussion. 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, MPG and HGF; the Italian INFN and INAF; the Swiss National Fund SNF; the grants PID2019-104114RB-C31, PID2019-104114RB-C32, PID2019-104114RB-C33, PID2019-105510GB-C31, PID2019-107847RB-C41, PID2019-107847RB-C42, PID2019-107847RB-C44, PID2019-107988GB-C22 funded by MCIN/AEI/10.13039/501100011033; 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-400/18.12.2020 and the Academy of Finland grant nr. 320045 is gratefully acknowledged. This work was also been supported by Centros de Excelencia “Severo Ochoa” y Unidades “María de Maeztu” program of the MCIN/AEI/ 10.13039/501100011033 (SEV-2016-0588, SEV-2017-0709, CEX2019-000920-S, CEX2019-000918-M, MDM-2015-0509-18-2) and by the CERCA institution of the Generalitat de Catalunya; by the Croatian Science Foundation (HrZZ) Project IP-2016-06-9782 and the University of Rijeka Project uniri-prirod-18-48; by the DFG Collaborative Research Centers SFB1491 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’Etudes Spatiales in France. This work performed in part under DOE Contract DE-AC02-76SF00515. The authors acknowledge the use of public data from the Swift data archive. This work makes use of observations from the KVA telescope. This work makes use of observations from the Liverpool Telescope, operated on the island of La Palma by Liverpool John Moores University in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias with financial support from the UK Science and Technology Facilities Council. This work makes use of observations from the Las Cumbres Observatory global telescope network. This article is based on observations made in the Observatorios de Canarias del IAC with the telescope PIRATE operated on the island of Tenerife by the Open University in the Observatorio del Teide. This research has made use of data from the OVRO 40-m monitoring program (Richards et al. 2011), supported by private funding from the California Insitute of Technology and the Max Planck Institute for Radio Astronomy, and by NASA grants NNX08AW31G, NNX11A043G, and NNX14AQ89G and NSF grants AST-0808050 and AST- 1109911. Part of this work is based on archival data, software, or online services provided by the Space Science Data Centre – ASI. D.P. acknowledges support from the Deutsche Forschungs gemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311. T.H. was supported by the Academy of Finland projects 317383, 320085, 322535, and 345899. S.K. acknowledges support from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme under grant agreement No. 771282. W.M. gratefully acknowledges support by the ANID BASAL projects ACE210002 and FB210003, and FONDECYT 11190853. R.R. acknowledges support from ANID BASAL projects ACE210002 and FB210003, and ANID Fondecyt 1181620. Author contribution: J. A. Acosta Pulido: Las Cumbres, PIRATE and Liverpool telescope data reduction; D. Dorner: project management, coordinating MAGIC observations, theoretical interpretation, coordination of MWL data analysis, drafting the manuscript; V. Fallah Ramazani: modelling, optical data analysis, theoretical interpretation, drafting the manuscript; T. Hovatta: OVRO data analysis; S. Kiehlmann: OVRO data analysis; C. Leto: Swift data analysis; D. Morcuende: VHE γ-ray data analysis, redshift estimation, drafting the manuscript; J. Otero-Santos: VHE γ-ray data analysis, redshift estimation, modelling, variability and cross-correlation analysis, drafting the manuscript; D. Paneque: coordinating MAGIC observations, theoretical interpretation, Fermi-LAT data analysis, coordination of MWL data analysis, drafting the manuscript; M. Perri: Swift data analysis; E. Prandini: redshift estimation; A. C. S. Readhead: OVRO data analysis; F. Verrecchia: Swift data analysis. The rest of the authors have contributed in one or several of the following ways: design, construction, maintenance and operation of the instrument(s) used to acquire the data; preparation and/or evaluation of the observation proposals; data acquisition, processing, calibration and/or reduction; production of analysis tools and/or related Monte Carlo simulations; discussion and approval of the contents of the draft.

References

Appendix A: MWL light curve description

The top panel of Fig. 1 shows the 30-day binned MAGIC light curve. The average integral flux above 100 GeV of the complete data set was estimated to be (2.96 ± 0.22) × 10⁻¹¹ cm⁻² s⁻¹, which corresponds to (6.2 ± 0.5)% C.U. During its non-flaring activity, the source showed a VHE γ-ray flux above 100 GeV of (0.97 ± 0.24) × 10⁻¹¹ cm⁻² s⁻¹, equivalent to (2.0 ± 0.5)% C.U. In epoch E2, 1ES 0647+250 experienced an increase of its γ-ray emission during its first high activity state, with an integral flux of (1.62 ± 0.78) × 10⁻¹¹ cm⁻² s⁻¹, corresponding to (3.4 ± 1.6)% of the Crab Nebula flux. Moreover in epoch E3, it displayed another bright state during the enhanced X-ray activity, with a flux of (3.82 ± 0.88) × 10⁻¹¹ cm⁻² s⁻¹ (equal to (8.0 ± 1.8)% C.U.) above 100 GeV. The highest VHE γ-ray emission of this blazar was observed in E4, with a γ-ray flux above 100 GeV of (7.10 ± 0.45) × 10⁻¹¹ cm⁻² s⁻¹, which corresponds to (15.0 ± 1.0)% C.U.

Regarding the Fermi-LAT light curve shown in the second panel of Fig. 1, a 30-day binning was also used. The average flux of 1ES 0647+250 between 300 MeV and 300 GeV was estimated to be (0.79 ± 0.03) × 10⁻⁸ cm⁻² s⁻¹. An increasing trend in the flux can be identified in the light curve, reaching a maximum flux of (2.76 ± 0.45) × 10⁻⁸ cm⁻² s⁻¹ at MJD 58115 (2017), a factor of 3.5 times higher than the average flux. During the VHE γ-ray non-flaring state observed during E1, this blazar also showed its lowest HE γ-ray emission, with an average flux of (0.43 ± 0.06) × 10⁻⁸ cm⁻² s⁻¹. Concerning the different enhanced states, the 30-day binning of the light curve does not allow for the performance of a comparison of the simultaneous emission between the HE and VHE γ-ray bands.

As for the Swift-XRT X-ray light curve, the flux in the energy band between 0.3 keV and 2 keV ranges from a value of (1.47 ± 0.06) × 10⁻¹¹ erg cm⁻² s⁻¹ during its minimum up to a historical maximum emission of (12.30 ± 0.25) × 10⁻¹¹ erg cm⁻² s⁻¹ coincident with the high state detected in E3. The average flux derived from all observations was estimated to be (2.98 ± 0.01) × 10⁻¹¹ erg cm⁻² s⁻¹. During the low state observed in epoch E1 the average X-ray flux is (2.22 ± 0.02) × 10⁻¹¹ erg cm⁻² s⁻¹. Moreover, the quasi-simultaneous XRT observations to the 2014 flare show an emission of (4.92 ± 0.11) × 10⁻¹¹ erg cm⁻² s⁻¹. The same behaviour is seen in the 2-10 keV energy band. Finally, the average X-ray emission during the flare from E4 was estimated to be (6.67 ± 0.06) × 10⁻¹¹ erg cm⁻² s⁻¹.

Swift-UVOT observed the source with almost the same temporal coverage as the XRT instrument. This instrument reported the flux of the source over the different periods in three UV and three optical filters. Here, we report the UW1 filter, since it is the one with the highest coverage out of the three UV filters of the instrument. The emission in the UV band ranges from a minimum flux of (0.65 ± 0.04) mJy during the first campaign up to a maximum of (1.42 ± 0.07) mJy coincident with the 2019 flare. The average emission during the non-flaring state monitored during E1 was (0.74 ± 0.01) mJy. The flux increased in E2 up to an average of (1.17 ± 0.03) mJy, a similar value to the flare from E3, (1.12 ± 0.02) mJy. The highest average emission was detected during the flare from E4, with an average flux of (1.20 ± 0.02) mJy.

The R-band optical light curve shows the same long-term increasing trend observed in the HE γ-ray emission of 1ES 0647+250. The optical flux ranges from (0.94 ± 0.04) mJy to (4.41 ± 0.07) mJy. Between 2009 and 2020, the average optical emission of 1ES 0647+250 was (1.88 ± 0.01) mJy. During the VHE γ-ray low state of epoch E1, the average optical flux was (1.76 ± 0.01) mJy. The emission of this blazar increased over the years, reaching a flux in the R band of (2.57 ± 0.05) mJy during the flaring activity of epoch E2. The optical emission reached its maximum in February 2019, and faded down to a flux density of (2.66 ± 0.05) mJy during the high state detected by MAGIC in E3. Finally, the activity during epoch E4 was similar to the previous period, with an average flux density of (2.72 ± 0.01) mJy.

The activity in radio wavelengths also displays the same increase seen in HE γ-ray and optical wavelengths. The flux detected by OVRO goes from a quiescent state of (0.043 ± 0.007) Jy up to (0.110 ± 0.002) Jy during its highest emission. The average flux at 15 GHz during the 11-year time span was estimated to be (0.090 ± 0.001) Jy. The mean radio emission during the VHE γ-ray low activity level shows a flux of (0.060 ± 0.001) Jy. During the flare observed in epoch E2, the most simultaneous OVRO observation to the MAGIC detection reports a flux of (0.106 ± 0.012) Jy. This corresponds to a factor of ∼1.8 compared to the average emission during the low state. As for the high state of E3, the measured flux at 15 GHz was (0.101 ± 0.001) Jy. Finally, the observation performed by OVRO in epoch E4 reported a flux density of (0.097 ± 0.002) Jy.

Appendix B: Radio and optical long-term trends

Here we show the long-term trends fitted for the radio and optical R-band light curves of 1ES 0647+250 with the method from Lindfors et al. (2016).

Fig. B.1. Estimated long-term trends for the 15 GHz radio and R-band optical light curves with the method from Lindfors et al. (2016).

Appendix C: Estimation of the distance between emitting regions

The distance between two physically separated emitting regions can be derived as described in Sect. 7. Instead of using the small angle approximation (θ ∼ 1/Γ and therefore δ ∼ Γ), the Doppler factor can be calculated as

$$\begin{array}{c}\hfill \delta =\frac{1}{\mathrm{\Gamma }(1-\beta cos\theta )},\end{array}$$(C.1)

where previous information on the viewing angle is needed. In our case there is no estimation on the angle. Therefore, for this calculation we assume a typical value of 1°. This leads to an estimated distance of d = 1.5 ± 0.2 pc, d = 7.1 ± 0.8 pc, and d = 29.1 ± 3.4 pc for a bulk Lorentz factor of 2, 4, and 8, respectively, following the same approach used in Sect. 7. Therefore, using the Doppler factor resulting of Eq. C.1, with θ = 1° yields a distance larger by approximately a factor of 2 with respect to the θ ∼ 1/Γ and δ ∼ Γ approximations.

All Tables

All Figures

Fig. 1. MWL light curves of 1ES 0647+250. From top to bottom: MAGIC (> 100 GeV, 30-day binning), Fermi-LAT flux and spectral index (300 MeV–300 GeV, 30-day binning), Swift-XRT (0.3−2 keV and 2−10 keV), Swift-XRT spectral index, Swift-UVOT (B, V, U, UVW1, UVM2, and UVW2 filters), optical R band (KVA, Liverpool Telescope, Las Cumbres, and PIRATE telescope), and OVRO (15 GHz). We note that all the optical observations except those from 2020 were performed by KVA. The smaller error bars of the OVRO light curve after 2016 are due to a major upgrade of the instrument. Blue contours correspond to the MAGIC campaigns during the different observed states.

In the text

	Fig. 2. Fractional variability of 1ES 0647+250. Filled markers represent MWL observations quasi-simultaneous to the MAGIC observations. Open markers correspond to the Fvar of the complete data sets. Filled markers represent the Fvar using only the simultaneous MWL data with respect to the MAGIC observations.
In the text

Fig. 3. Long-term cross-correlation curves. The ZDCF is represented in black and its 1σ uncertainty by the grey contour. Coloured dotted lines correspond to different significance levels, from 1σ to 4σ. Left: cross-correlation between the optical and HE γ-ray light curves. Middle: cross-correlation between the radio and optical light curves. Right: cross-correlation between the radio and HE γ-ray light curves.

In the text

	Fig. 4. Light curve detrending results. Blue dashed lines represent the scaled low-frequency light curve best fit, red dots correspond to the original data set and black dots show the detrended light curve. Left: Fermi-LAT light curve after optical trend subtraction. Middle: optical light curve after radio trend subtraction. Right: Fermi-LAT light curve after radio trend subtraction.
In the text

	Fig. 5. PL index αXRT vs X-ray flux for the different epochs. Left: 0.3−2 keV flux. Right: 2−10 keV flux.
In the text

Fig. 6. Broadband SEDs of 1ES 0647+250 for the different observing campaigns described with one-component and two-component SSC emission models. Top left: 2009–2011. Bottom left: 2014. Top right: 2019. Bottom right: 2020. Filled blue dots correspond to the MWL simultaneous data of each period. Empty dots represent the archival data extracted from the SSDC database (https://www.ssdc.asi.it). The dotted lines correspond to the fitted one-component SSC model. Black lines represent the total two-component SSC model. The dashed and dotted-dashed lines correspond to the blob and core components of the two-component model, respectively.

In the text

	Fig. B.1. Estimated long-term trends for the 15 GHz radio and R-band optical light curves with the method from Lindfors et al. (2016).
In the text