Long-term multi-wavelength study of 1ES 0647+250

The BL Lac object 1ES 0647+250 is one of the few distant $\gamma$-ray emitting blazars detected at very high energies (VHE, $\gtrsim$100 GeV) during a non-flaring state. It was detected with the MAGIC telescopes during its low activity in the years 2009-2011, as well as during three flaring activities in the years 2014, 2019 and 2020, with the highest VHE flux in the latter epoch. An extensive multi-instrument data set was collected within several coordinated observing campaigns throughout these years. We aim to characterise the long-term multi-band flux variability of 1ES 0647+250, as well as its broadband spectral energy distribution (SED) during four distinct activity states selected in four different epochs, in order to constrain the physical parameters of the blazar emission region under certain assumptions. We evaluate the variability and correlation of the emission in the different energy bands with the fractional variability and the Z-transformed Discrete Correlation Function, as well as its spectral evolution in X-rays and $\gamma$ rays. Owing to the controversy in the redshift measurements of 1ES 0647+250 reported in the literature, we also estimate its distance in an indirect manner through the comparison of the GeV and TeV spectra from simultaneous observations with Fermi-LAT and MAGIC during the strongest flaring activity detected to date. Moreover, we interpret the SEDs from the four distinct activity states within the framework of one-component and two-component leptonic models, proposing specific scenarios that are able to reproduce the available multi-instrument data.


Introduction
Blazars are radio-loud active galactic nuclei whose relativistic jets point towards the Earth.Blazars can be classified accord-Corresponding authors; e-mail: contact.magic@mpp.mpg.de.
ing 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 A&A 670, A49 (2023) 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 14 Hz), intermediate-(10 14 Hz < ν peak < 10 15 Hz), and high-energy-peaked BL Lacs (HBLs; ν peak > 10 15 Hz; Padovani & Giommi 1995).Another category of BL Lacs was introduced by Costamante et al. (2001), naming those whose peak is above ν peak > 10 17 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.11  −0.10 , 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 nondetections 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 multiwavelength (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.

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 −11 cm −2 s −1 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 A49, page 2 of 20  1 shows the significances of the detection during the different epochs as estimated following Eq.( 17) in Li & Ma (1983).

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.fitsand iso_P8R3_SOURCE_V2_v1.txt,respectively 1 .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(Abdollahi et al. , 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 1 https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html 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)2 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.
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.

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

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.

Optical data
Optical monitoring of the source in the R band was also performed by the Tuorla blazar monitoring programme 3 .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, andcontinued 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).

Radio observations: OVRO
1ES 0647+250 is also part of the Owens Valley Radio Observatory (OVRO) blazar monitoring programme (Richards et al.   3 http://users.utu.fi/kani/1m/2011) 4 .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 nondetections and thus were not included in the MWL light curve and analysis.This resulted in 446 observations after excluding these measurements from the analysis.

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.

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(Abdollahi et al. , 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 χ 2 /d.o.f.= 346/143 2.4 (p value = 10 −20 ).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 χ 2 /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 longterm 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 2009Kapanadze , 2019Kapanadze , 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 χ 2 /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  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 where x and S 2 are the mean and the variance of the distribution of measured fluxes, respectively, and σ 2 err 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): where err(σ 2 NXS ) is the normalised excess variance taken from Vaughan et al. (2003), calculated as 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.
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(Ahnen et al. , 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 doublemaximum structure is also reported for Mrk 501).The discrepancy between the R-band observations and the UV and optical A49, page 6 of 20 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.

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×105 artificial light curves with the same sampling and power spectral density as our data set, as described in Emmanoulopoulos et al. (2013)  5 .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. (1998Peterson et al. ( , 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.
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 χ 2 /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)6 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: 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  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 χ 2 /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 shortterm correlation analysis of each individual observing epoch.

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 0 • (E/E 0 ) −α ).A LogP fit was also tested (dN/dE = f 0 • (E/E 0 ) −α−β • log(E/E 0 ) ).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 21 cm −2 (HI4PI Collaboration 2016).However, there is no statistical preference for the LogP over the PL, with χ 2 /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 −4 .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).
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 A49, page 8 of 20 Table 3. Spectral parameters of the X-ray spectrum assuming a simple PL shape.

Epoch
Time interval 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 −11 erg cm −2 s −1 for the band between 0.3 keV and 2 keV, and F < 4 × 10 −11 erg cm −2 s −1 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 −3 .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.

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.
A49, page 9 of 20 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 χ 2 /d.o.f.= 126.8/1430.9, which is similar to α Fermi = 1.73 ± 0.02 reported in the 4FGL-DR3 catalogue (Abdollahi et al. 2020(Abdollahi et al. , 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.

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.

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.

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. 2017Petropoulou et al. , 2020a,b;,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 onecomponent 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 twocomponent 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 0 = 67.4km s −1 Mpc −1 , a matter density of Ω m = 0.315 and a dark energy density Ω Λ = 0.685.The models used here are not timedependent.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.

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 lowenergy 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 A49, page 10 of 20 Table 4. Spectral parameters of the HE γ-ray spectra from Fermi-LAT data.

Epoch
Time interval Notes. (* ) For the 2020 VHE γ-ray spectrum, both PL and LogP fit models were used, with a 3σ preference for the latter.These functions are specified in Sect.4.1.
broken PL distribution with the Lorentz factor, described by The distribution has a normalisation K between γ min and γ max and slopes n 1 and n 2 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.

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 synchrotronself absorption.The strength of the magnetic field is typically assumed to scale with the distance from the central engine as d −1 .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 mag-netic 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 19 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).Table 6.SED modelling parameters for one-component SSC and two-component models. (1) (2)  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 16 −4.2× 10 16 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 A49, page 12 of 20 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.

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.

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.

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 θ)] −1 (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 selfabsorption 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 A49, page 13 of 20  3)-( 5) Kinetic power of the magnetic field, electrons and cold protons, respectively (for the core in the case of the two-component model).( 6) Total kinetic power of the jet.( 7)-( 9) Jet power carried by the jet in form of magnetic field, electrons and cold protons, respectively.( 10) Total power carried by the jet.
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. (2018Liodakis et al. ( , 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.

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-whenbrighter 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).How-ever, 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 harderor 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.

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 A49, page 14 of 20 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 onecomponent 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.

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 twocomponent 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 highenergy 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.

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: A49, page 15 of 20 -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.

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 F var of the complete data sets.Filled markers represent the F var using only the simultaneous MWL data with respect to the MAGIC observations.

Fig. 3 .
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: crosscorrelation between the radio and optical light curves.Right: cross-correlation between the radio and HE γ-ray light curves.

Fig. 4 .
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.

Fig. 6 .
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.

Table 1 .
Significance and integrated flux above 100 GeV of the different detections of 1ES 0647+250.

Table 2 .
Goodness of the constant flux hypothesis for every MWL light curve.

Table 5 .
Notes.The TS value is the likelihood test statistic resulting from the fit to the model.Spectral parameters of the VHE γ-ray spectra from MAGIC data.

Table 7 .
Results of the SED modelling for one-component SSC and two-component models.