---
abstract: 'The MAGIC telescopes observed S2 0109+22 in 2015 July during its flaring activity in high energy gamma rays observed by *Fermi*-LAT. We analyse the MAGIC data to characterise the very high energy (VHE) gamma-ray emission of S2 0109+22, which belongs to the subclass of intermediate synchrotron peak (ISP) BL Lac objects. We study the multi-frequency emission in order to investigate the source classification. Finally, we compare the source long-term behaviour to other VHE gamma-ray emitting (TeV) blazars. We performed a temporal and spectral analysis of the data centred around the MAGIC interval of observation (MJD 57225-57231). Long-term radio and optical data have also been investigated using the discrete correlation function. The redshift of the source is estimated through optical host-galaxy imaging and also using the amount of VHE gamma-ray absorption. The quasi-simultaneous multi-frequency spectral energy distribution (SED) is modelled with the conventional one-zone synchrotron self-Compton (SSC) model. MAGIC observations resulted in the detection of the source at a significance level of $5.3\,\sigma$. The VHE gamma-ray emission of S2 0109+22 is variable on a daily time scale. VHE gamma-ray luminosity of the source is lower than the average of TeV BL Lacs. The optical polarization, and long-term optical/radio behaviour of the source are different from the general population of TeV blazars. All these findings agree with the classification of the source as an ISP BL Lac object. We estimate the source redshift as $z = 0.36 \pm 0.07$. The SSC parameters describing the SED are rather typical for blazars.'
bibliography:
- 'ref.bib'
date: 'Accepted XXX. Received YYY; in original form ZZZ'
title: 'The broad-band properties of the intermediate synchrotron peaked BL Lac S2 0109+22 from radio to VHE gamma rays'
---

\[firstpage\]

galaxies: active – galaxies: jets – gamma rays: galaxies – BL Lacertae objects: individual: S2 0109+22

Introduction
============

BL Lac objects dominate the extragalactic very-high-energy (VHE, $E>100$GeV) gamma-ray sky. A relativistic jet shoots from the region of the central super-massive black hole, hosted at the center of BL Lac objects, in the line of sight of the observer. Jets are typically characterized by featureless spectra in the optical band, highly polarized radiation in radio and optical, and variable radiation at all frequencies. The jet emission is non-thermal and described as a continuous spectral energy distribution (SED), spanning from radio to VHE gamma-ray frequencies, and featuring two wide peaks. Synchrotron emission by highly relativistic electrons spiralling in the magnetic field of the jet is used to explain the lower frequency peak. Different scenarios within various models are used to explain the high-frequency peak: external Compton [@1989ApJ...340..162M; @1994ApJS...90..945D; @1994ApJ...421..153S] and synchrotron self-Compton [SSC, @1992MNRAS.258..657C; @1992ApJ...397L...5M] as leptonic models, proton synchrotron emission [@1996SSRv...75..331M; @2000NewA....5..377A; @2001APh....15..121M] and photo-pion production [@2014ApJ...797...89A] as hadronic models. Traditionally, in view of their relative simplicity and agreement with the data, single-zone SSC models have been used to describe BL Lac SEDs [e.g. @2011ApJ...727..129A; @2011ApJ...736..131A]. However, there is growing evidence that these models do not reproduce all the observed features of BL Lac objects , and, in some cases, more complicated models should be considered. BL Lac objects are classified according to the peak frequency of their lower energy peak, $\nu_{syn}$ [@1994MNRAS.268L..51G]: low synchrotron peaked (LSP; $\nu_{syn}<10^{14}$Hz), intermediate synchrotron peaked (ISP; $10^{14}\leq\nu_{syn}<10^{15}$Hz), and high synchrotron peaked (HSP; $\nu_{syn}\geq10^{15}$Hz) [@2010ApJ...716...30A].

S2 0109+22 (also known as GC 0109+224), at coordinates (J2000) $\text{RA = 01h12m05.8s}$ and $\text{DEC=+22d44m39s}$, was first detected as a compact radio source in the 5GHz Survey of the NRAO 43m dish of Green Bank, West Virginia [@1971AJ.....76..980D; @1972AJ.....77..265P]. In 1976, it was optically identified as a stellar object of magnitude 15.5 on the Palomar Sky Survey plates, @1977AJ.....82..776O also measured a strong millimetre emission (1.53Jy at 90GHz)[^1] and defined it as a BL Lac object. Since then it was continuously monitored in radio and optical . performed extensive studies on the radio and optical behaviour and the broad-band SED of this source. It remarkably shows high polarization variability, from 7% to 30% . It is classified as an ISP BL Lac object using different approaches and datasets to calculate the location of its synchrotron peak.

Since the launch of the *Fermi* satellite in 2008, the source has been listed in most of the *Fermi*-LAT catalogues, i.e. 1FGL [@2010ApJS..188..405A]; 2FGL [@2012ApJS..199...31N]; 1FHL [@2013ApJS..209...34A]; and 3FGL [@2015ApJS..218...23A]. However, the source is not listed in the catalogue of sources detected $>50\,\text{GeV}$ by the *Fermi*-LAT [2FHL, @2016ApJS..222....5A]. The source is variable in the high energy (HE: 100 MeV$<E<$100 GeV) gamma-ray band with the variability index equal to 489 and the maximum monthly flux value of $F_{(0.1-100\,\text{GeV})}=(2.14\pm0.17) \times 10^{-7}$phcm$^{-2}$s$^{-1}$ which is reported in February 2011 [3FGL, @2015ApJS..218...23A]. @2008ApJS..175...97H reported a redshift value for the source of $z=0.265$, which was disfavoured by @2016MNRAS.458.2836P using a high signal-to-noise optical spectrum from Gran Telescopio Canarias. Based on this spectrum, $z>0.35$ was measured, assuming the source is hosted by a massive elliptical galaxy typical for this class of sources. VHE gamma-ray observations of this source carried out with MAGIC between 2015 July 22 and 28 (MJD 57225–57231), were triggered when the reported HE gamma-ray daily flux, July 20 (MJD 57223), was about two times higher than the average flux reported in the 3FGL catalogue (private communication with Luigi Pacciani). The MAGIC observations led to the first detection of this source in VHE gamma rays [@2015ATel.7844....1M].

In this paper, we present the multi-frequency observations and data analysis in Section \[sec2\]. A long-term behaviour study, the comparison with other VHE gamma-ray emitting (TeV) blazars, and estimations of the source distance are presented in Section \[sec3\]. Finally, Section \[sec4\] summarizes our results.

Observations and data analysis {#sec2}
==============================

In this section, we introduce the instruments and their respective data analysis procedures.

Very high energy gamma rays (MAGIC)
-----------------------------------

MAGIC is a system of two Imaging Atmospheric Cherenkov Telescopes (17m diameter) located in the Canary Island of La Palma ($28.7^{\circ}$N, $17.9^{\circ}$W), at the elevation of 2200m a.s.l. [@2016APh....72...76A]. The use of the stereoscopic technique, combined with large mirror size makes MAGIC one of the most sensitive instruments for VHE gamma-ray astronomy. The corresponding trigger threshold is $\gtrsim 50$GeV [@2016APh....72...76A]. S2 0109+22 is visible from the MAGIC site at zenith angle below 40$^{\circ}$ between mid-July and February.

Triggered by increased activity in HE gamma rays, MAGIC observed S2 0109+22 for 9.63h in 2015 July within a multi-wavelength blazar monitoring program. The observations were performed during 7 consecutive nights from July 22 to July 28 (MJD 57225–57231) with zenith angle range between 11$^\circ$ and 39$^\circ$. The data have been analysed using the MAGIC Standard Analysis Software [MARS, @2009arXiv0907.0943M; @2016APh....72...61A; @2017MNRAS.468.1534A]. Part of the data were affected by clouds, therefore we applied atmospheric transmission correction based on the information obtained with the MAGIC elastic LIDAR [@2015EPJWC..8902003F].

High-energy gamma rays (*Fermi*-LAT)
------------------------------------

The Large Area Telescope (LAT) is the primary instrument on-board the *Fermi Gamma-ray Space Telescope*. Based on the pair-conversion technique, it is designed to investigate the gamma-ray sky in the energy band from 30MeV to &gt;300GeV [@2009ApJ...697.1071A]. In its standard operation mode it surveys the sky, covering it fully every 3h.

The data analysed in this paper were selected from a region of interest around S2 0109+22 with a radius of $15^{\circ}$, in a period lasting around three weeks (MJD 57220–57240) roughly centred on the MAGIC detection peak on MJD 57228 (2015 July 25). The data reduction of the events of the Pass8 source class was performed with the ScienceTools software package version v10r0p5[^2] in the energy range 0.1–300GeV. To reduce Earth limb contamination a zenith angle cut of 90$^{\circ}$ was applied to the data. The un-binned likelihood fit of the data was performed using the suggested Galactic diffuse-emission model and isotropic component [@2016ApJS..223...26A] recommended for Pass8 Source event class[^3].

The normalizations of both diffuse components in the source model were allowed to freely vary during the spectral fitting. The source model also includes the sources of the *Fermi*-LAT third source catalogue [3FGL, @2015ApJS..218...23A] within 25$^{\circ}$ of the source of interest. Spectral indices and fluxes are left to freely vary for sources within 5$^{\circ}$; fluxes are also left to freely vary for sources flagged as ’variable’ in the 3FGL catalogue that lie from 5 to 10$^{\circ}$. The spectral parameters of the sources from 10 to 25$^{\circ}$, were instead fixed to their catalogue value.

To construct the light-curve (LC) with 1-day time bins, only the source of interest (normalization and spectral index) and the diffuse models (normalization) were left free to vary, while the remaining 3FGL sources were fixed to the values obtained for the three week analysis of the region. An upper-limit is shown when the detection significance was $<3\,\sigma$[^4]. The SED was obtained analysing data collected between the 2015 July 22 and 2015 July 28 (MJD 57225–57231), corresponding to the MAGIC observing period.

X-ray and UV (*Swift*)
----------------------

Since 2006, *Neil Gehrels Swift observatory (Swift)* has pointed to the source fifteen times in photon counting mode. Ten of the raw images by the X-ray Telescope [*XRT*, @2004SPIE.5165..201B] on-board the *Swift* satellite, are qualified for analysis[^5]. The multi-epoch (8) event list for the period from 2015 July 21 (MJD 57224.95) to 2015 August 1 (MJD 57235.86) with a total exposure time of $\sim 6.15$h, were downloaded from the publicly available SWIFTXRLOG (*Swift*-XRT Instrument, Log)[^6]. These observations have an average integration time of 2.8ks each. They were processed using the procedure described by , assuming fixed equivalent Galactic hydrogen column density $n_H = 4.24 \times 10^{20}\,\rm cm^{-2}$ reported by . Additionally, *Swift* observed this source two more times in 2006. We analysed those two additional event lists to get a broader view of the source’s X-ray properties.

The Ultraviolet/Optical Telescope (*UVOT*, $[4.9-16.6] \times 10^{5}$GHz) on-board the *Swift* satellite [@2008MNRAS.383..627P], observed the source 15 times during the MAGIC campaign, out of which eight were simultaneous to the *XRT* data taking[^7]. An iterative data calibration procedure was used to calculate the Galactic extinction[^8], the effective frequency, and the flux conversion factor for each filter.

Optical
-------

### Light-curve (KVA, KAIT, and Catalina)

S2 0109+22 was added to the Tuorla blazar monitoring program[^9] when HE activity was reported in 2015 July. The monitoring observations were performed in optical R-band using a 35cm Celestron telescope coupled to the KVA (Kunglinga Vetenskapsakademi) telescope located at La Palma. Data analysis was performed using a semi-automatic pipeline for differential photometry assuming the comparison star magnitudes in . The magnitudes were corrected for Galactic extinction using values from @2011ApJ...737..103S.

In order to study the long-term optical behaviour of S2 0109+22, its optical LC is retrieved from the publicly available online database of 76-cm Katzman Automatic Imaging Telescope (KAIT) at Lick Observatory[^10]. The LC from KAIT is produced through a pipeline that utilizes aperture photometry and performs brightness calibrations using USNO B1.0 catalogue stars in the source field. The long-term optical LC is extended back to 2005 by including available online data from the Catalina Real-Time Transient Survey [@2009ApJ...696..870D]. KAIT and Catalina data are obtained from unfiltered observations, whose effective color is close to the R-band [@2003PASP..115..844L].

### Host galaxy imaging (NOT) {#HGIMAGE}

To investigate the host galaxy of S2 0109+22, we obtained a deep I-band image at the Nordic Optical Telescope (NOT) on 2015 November 11. In total, 26 exposures, each 150 seconds long, were obtained using the ALFOSC[^11] instrument. After subtracting the bias, flat-fielding and fringe map correction, the images were registered using stars in the field and summed. The resulting image has a total exposure time of 3900 seconds with $FWHM \cong 1.14\arcsec$. The comparison stars in were used to calibrate the field.

### Polarization (NOT)

Polarization observations were carried out using the ALFOSC instrument in the standard linear polarization set-up (lambda/2 retarder followed by calcite) in optical *R*-band. Weekly observations were performed from November 2015 to September 2017 within three observing seasons. In order to determine the zero point of the position angle, polarization standards were observed on a monthly basis. The instrumental polarization was measured observing zero-polarization standard stars, and was found to be negligible. Most of the observations were conducted under good sky condition (seeing $\sim 1\arcsec$).

Using aperture (radius of $1.5\arcsec$) photometry, the sky-subtracted target counts were measured for ordinary and extraordinary beams. By using the intensity ratios of two beams and standard formulae in @2007ASPC..364..495L, we calculated normalized Stokes parameters, polarization fraction, and position angle for each observation. Systematic uncertainties are included in our error estimation.

Radio (OVRO and Metsähovi)
--------------------------

S2 0109+22 was observed at 15GHz as part of a high-cadence gamma-ray blazar monitoring program using the Owens Valley Radio Observatory (OVRO) 40 m telescope [@2011ApJS..194...29R]. The observations are calibrated by using a temperature-stable diode noise source to remove receiver gain drifts, and the flux density scale is derived from observations of 3C 286 assuming the value of 3.44Jy at 15.0GHz . The systematic uncertainty of about 5% in the flux density scale is not included in the error bars. Complete details of the reduction and calibration procedure are found in @2011ApJS..194...29R.

The Mets[ä]{}hovi radio telescope, operating at 37GHz, has been observing the source for two decades. We selected radio data obtained after mid-2005 for the long-term study of the source. The instrument and data reduction procedures are described by .

Results {#sec3}
=======

Very high energy gamma rays {#sec3-1}
---------------------------

The VHE gamma-ray signal from the source is estimated after applying energy dependent selection cuts to the signal. Residual background of the observation is measured around a control region [@2017MNRAS.468.1534A]. The distribution of the events is shown in Figure \[theta2\]. In total, there was an excess of ($365.8\pm 69.1$) events in the signal region $\theta^2 < 0.02\,\deg^2$, where $\theta^2$ is the squared angular distance between the reconstructed source position of the events and the nominal position of the expected source. The data taken during MJD 57228 (2015 July 25) contribute $\geq 61\%$ of excess events of the whole sample of data. The source was detected at a significance level of $7.24\,\sigma$ during MJD 57228 (Fig. \[theta2\]).

The LC of the VHE gamma-ray integral flux above 100GeV ($F_{>100\,\text{GeV}}$) is shown in Figure \[fig\_lc\] with the details presented in Table \[tab\_tev\_flux\]. The constant flux hypothesis is disfavoured with $\chi^2/\text{d.o.f.} = 14.5/4$ ($P_{value}=0.005$). The peak flux, detected on MJD 57228 (hereafter flare night), is twice the average flux over the whole period of observation, $F_{>100\,\text{GeV, ave}}=~(4.7 \pm 1.2) \times 10^{-11}$phcm$^{-2}$s$^{-1}$. Real correlation analysis for such a short period around the flare night is beyond the reach with the available data sample shown in Figure \[fig\_lc\]. However, there seems to be an increased flux in X-rays, optical and UV bands around the flare night, which suggests that emission in these bands could originate from a single region.

We compare the integral flux ($F_{>200\,\text{GeV}}$) of S2 0109+22 to that of other TeV BL Lac objects (21 sources) with variable flux in VHE gamma rays presented in the most recent population study by , who studied a time independent correlation between several lower-energy bands and VHE gamma-ray luminosity, and predicted the VHE gamma-ray flux for 182 non-TeV BL Lac objects. The comparison is shown in Figure \[histo2\]. Both the lowest and the largest observed flux of S2 0109+22 during the MAGIC campaign are among the faintest of the population. High and low state VHE gamma-ray predicted energy flux (&gt;200GeV) in are $(4.5\pm1.9)\times 10^{-12}$ and $(9.8\pm2.1)\times 10^{-14}$ergcm$^{-2}$s$^{-1}$, respectively. The largest observed flux over the same energy range, $F_{>200\,\text{GeV}}^{\text{high~obs}}=(4.6\pm1.5)\times 10^{-12}$ergcm$^{-2}$s$^{-1}$, is in good agreement with the predicted flux. The lowest observed flux of this source is $F_{(>200\,\text{GeV})}^{\text{low obs}}=(1.5\pm0.7)\times 10^{-12}$ergcm$^{-2}$s$^{-1}$. The observed VHE gamma-ray flux of the source is fainter than the sample of variable TeV BL Lacs.

Figure \[vhe\_sed\] shows the spectrum of S2 0109+22 in the VHE gamma rays. We assume a simple power-law model, $$\frac{dN}{dE}=F_0 (\frac{E}{E_{\text{dec}}})^{-\Gamma},$$ where $E_{\text{dec}}$ and $F_0$ are the decorrelation energy and differential flux at $E_{\text{dec}}$, and $\Gamma$ is the spectral photon index. The spectral parameters are obtained via forward-folding using Poissonian maximum likelihood procedure described by @2017MNRAS.472.2956A. In order to calculate the intrinsic spectral parameters, the same estimation procedure is used by assuming $z=0.35$ (see Sect. \[sec3-5\]) and Extragalactic Background Light (EBL) absorption model described by @2011MNRAS.410.2556D. The spectral parameters are summarized in Table \[tab\_tev\_spec\] for the flare night and the average spectrum. The fitted model statistics are calculated in the energy range of 65–370GeV and 65–250GeV for average and flare night spectra, where MAGIC detected the source.

High-energy gamma rays {#sec3-2}
----------------------

We have found that there is no significant HE gamma-ray spectral and flux variability on a daily basis during the investigated period (MJD 57220–57240). These results are shown in Figure \[fig\_lc\] (Panels b and c). The HE gamma-ray constant fit flux is $F_{(0.1-300\,\text{GeV})} = (1.4 \pm 0.4) \times 10^{-7}$phcm$^{-2}$s$^{-1}$, which is $\sim\,2$ times higher than the average flux reported in the 3FGL catalogue [@2015ApJS..218...23A] for this source.

To model the HE gamma-ray spectrum of S2 0109+22, a power-law function which uses integrated flux as a free parameter[^12] is used. $$\frac{dN}{dE}=\frac{N(\Gamma +1)E^{\Gamma}}{E^{\Gamma +1}_{\text{max}}-E^{\Gamma +1}_{\text{min}}}$$ where $\Gamma$ is the photon index, $E_{\text{min}}=100\,\text{MeV}$, $E_{\text{max}}=300\,\text{GeV}$, and N is the integral flux between $E_{\text{min}}$ and $E_{\text{max}}$.

We analysed the source in the period MJD 57225–57232, modelling its spectrum with a simple power-law. The likelihood fit obtained a Test Statistic of $TS= 111$. The resulting power law index of the fitted model is $\Gamma=1.81 \pm 0.14$. The spectral index of the investigated period is within the error bars of the one reported in 3FGL. In Figure \[vhe\_sed\], we show the flux values in six logarithmically spaced bins from 100MeV to 300GeV. Upper-limits are shown when the detection significances are lower than $3\,\sigma$.

X-rays {#xray3}
------

The results of our X-ray analysis are shown in Figure \[fig\_lc\] with the details available in Table \[tab\_xray\]. The X-ray flux peaks three nights before the VHE gamma-ray peak. The X-ray spectrum is usually soft (Photon index, $\Gamma_X \geq 2.4$). The constant flux hypothesis is rejected with $>10\,\sigma$ level of confidence. However, only a hint of brighter-harder trend with $2\,\sigma$ level of confidence is present in our data sample. The trend between X-ray spectral index and flux ($F_{0.3-10\,\text{keV}}$) can be described by a linear model (Fig. \[spec\_flux\]) with the test statistics of $\chi ^2/\text{d.o.f.}= 2.97/5$, corresponding to Pearson correlation coefficient of 0.76. Moreover, we tried to fit a log-parabola model to the data obtained on MJD 57228. It reveals that the power law model with an index $\Gamma_X = 2.58 \pm 0.05$ ($\chi ^2/\text{d.o.f.}= 48.2/50$) can describe the spectrum better. The X-ray flux ($F_{0.3-10\,\text{keV}}$) on the flare night was $>6$ times higher than the flux from 2006 observations.

Long-term behaviour
-------------------

Recently two studies of optical and radio behaviour of TeV blazars have been published. studied the long-term optical and radio behaviour of 32 VHE gamma-ray blazars using data from the OVRO and Tuorla blazar monitoring programs. They found correlated flares in half of the sources, and correlated long-term trends in 13 sources. performed a first statistical study of the optical polarization variability of TeV blazars, and found that they are not different from the control sample of non-TeV blazars. S2 0109+22 was not part of those studies. In order to compare its optical and radio behaviour with the sample of VHE gamma-ray blazars, we have performed the same analysis of the long-term optical and radio data and optical-polarization data as done in and .

Moreover, the long-term correlation studies between radio/optical and gamma-ray bands were already performed by @2014MNRAS.445..437M and @2014ApJ...797..137C using similar radio and optical datasets as those presented in this analysis. Therefore, we only attempt to study the long-term radio-optical cross-correlation behaviour of the source together with its optical polarization behaviour.

### Radio-Optical cross-correlation analysis {#3-4-1}

Figure \[fig\_longterm\] illustrates the long-term optical and radio data of S2 0109+22. The coverage is of 12 years in the optical band (R-band and open filters) and at 37GHz, and 10 years at 15GHz.

Following , we calculated the cross-correlation function between the optical and 15GHz LCs using the Discrete Correlation Function [DCF; @1988ApJ...333..646E] with local normalization [LCCF; @1999PASP..111.1347W]. We use temporal binning of 10 days and require that each LCCF bin has at least 10 elements. Following @2014MNRAS.445..437M, the significance of the correlation is estimated using simulated LC. In the simulations, we used a power spectral density index of -1.8 for the radio LC (Max-Moerbeck et al. in prep.), which is slightly smaller than the values between -1.4 and -1.7 reported in @2004MNRAS.348.1379C for the (8 to 37GHz) radio LCs. For the optical, we used a power spectral density index of -1.5 (Nilsson et al. in prep.). While there are several peaks (features) in the LCCF, shown in Figure \[lccf\], none of them reach the $2\,\sigma$ significance level. We also calculated the cross-correlation functions between the optical-37GHz and 37-15GHz. The only significant correlation is between 37-15GHz, with significance $>3\sigma$. The peak is rather broad from -40 to +30 days (Fig. \[lccf\]) and is consistent with zero lag. Typically, for evolving synchrotron self-absorbed components [e.g. @1994ApJ...437...91S; @2014MNRAS.441.1899F], one would expect the higher frequency to lead the lower frequency variations, which is consistent with our finding. However, as stated the peak is rather wide and also consistent with zero time lag. These results may indicate co-spatiality.

The optical-radio correlations of this source have been previously studied by and @2004MNRAS.348.1379C. Both works found several weak peaks in the correlations with lags 190, 400 days , and 190, 789 and 879 days [@2004MNRAS.348.1379C]. In Figure \[lccf\], there is a single ‘feature’ covering all these lags, peaking at $\sim$500 days. This feature is not significant and in general the results of our calculation agree with those by and @2004MNRAS.348.1379C.

We also searched for common long-term trends from the optical and radio data by fitting linear trends to these LCs. No long-term trends were found at these wavelengths.

We then compared the results of the correlation and trend analyses to the results obtained for other TeV blazars in . The sources in which no connection between flaring behaviour nor long term behaviour were found were a minority in that sample and were either very weak sources, or bright sources with clear outbursts like S2 0109+22. These other bright sources in the sample were S5 0716+714, ON 325 and W Com and it was suggested that as there were several $2\,\sigma$ peaks in their correlation function, there might be several time-scales involved, blurring the correlation. However, for S2 0109+22 we do not find any correlation peaks above $2\,\sigma$. This result may indicate that a major fraction of the optical flux in this source is not originating from the same emission region as the radio, or that the radio-optical correlation is more complex than can be probed by the simple cross-correlation function used in this paper.

### Optical Polarization

The optical emission in active galaxies is dominated by synchrotron emission of the their jet, which is intrinsically highly polarized. In an optically thin jet with uniform magnetic field, the polarization fraction can be up to 70% [e.g., @1970ranp.book.....P]. The more typically observed levels of fractional polarization reach a few tens of percent at maximum , which have been taken as evidence for disordered magnetic fields. The linearly polarized emission is described using the Stokes parameters *I* (for total intensity), and *Q* and *U* (for linear polarization). Using the Stokes parameters, the polarization fraction and the electric vector position angle (EVPA) can be defined as $m=(\sqrt{Q^2+U^2})/I$ and EVPA$=1/2\tan^{-1}(U/Q)$. The polarization fraction and EVPA for S2 0109+22 are shown in Figure \[pollc\].

We estimate the long-term polarization variability of S2 0109+22 by using the methods described in where the optical polarization of a sample of TeV and non-TeV-detected BL Lac objects was studied. We calculate the intrinsic mean polarization fraction and its modulation index (standard deviation of the polarization fraction over the mean), by assuming that the polarization fraction follows a Beta distribution, which is confined to values between 0 and 1, similarly as the polarization fraction. A single polarization observation is assumed to follow a Ricean distribution, so that our probability density function is obtained by convolving the Beta and Ricean distributions as follows,

$${\rm PDF}\left(p;\alpha, \beta\right)=\frac{p^{\alpha -1} \left( 1-p\right)^{\beta -1}}{B\left(\alpha, \beta\right)},$$ where $p$ is the polarization fraction and $\alpha$ and $\beta$ determine the shape of the Beta distribution $B\left(\alpha, \beta\right)$. If the parameters $a, \beta$ of this distribution are known, the mean polarization fraction and the intrinsic modulation index are then given by

$$p_\text{int}=\frac{\alpha}{\alpha + \beta}$$ and $$m_\text{int}=\frac{\sqrt{Var}}{p_\text{int}}=\frac{\sqrt{\frac{\alpha\beta}{\left(\alpha + \beta\right)^2 \left(\alpha + \beta
       +1 \right)}}}{\frac{\alpha}{\alpha + \beta}},$$ where $Var$ is the variance of the distribution. Details of the method are described in Appendix A of @2016MNRAS.457.2252B. The intrinsic mean polarization fraction of S2 0109+22 is $0.090^{+0.010}_{-0.008}$, which is higher than the sample mean values of $0.054\pm0.008$ and $0.079\pm0.009$ obtained for the TeV and non-TeV BL Lac objects in . Similarly, the intrinsic modulation index of the polarization fraction $0.54^{+0.08}_{-0.06}$ is higher than the sample mean values for the TeV ($0.29\pm0.03$) and non-TeV ($0.38\pm0.04$) sources.

The polarization angle variability can be quantified by calculating the derivative of the polarization angle variations. First we account for the n$\pi$ ambiguity of the polarization angle by requiring that each subsequent point is within 90$^\circ$ from the previous observation. We obtain a median derivative of 2.4 degrees per day, which translates to 3.3 degrees per day in the source frame when multiplied by $(1+z)$ ($z=0.35$, see Sect. \[sec3-5\]). Comparing this to the histograms in Figure 4 of shows how S2 0109+22 varies more rapidly in polarization angle than the average TeV (mean $1.11\pm0.29$ deg./day) and non-TeV (mean $1.66\pm0.45$ deg./day) sources. This is also seen when we examine the polarization variations in the $Q/I-U/I$-plane (see the inset in Fig. \[pollc\] for the $Q/I-U/I$ plot). As described in a tightly clustered distribution of the points in the $Q/I-U/I$-plane is an indication of a preferred polarization angle. For S2 0109+22 the weighted average of the $Q/I$ and $U/I$ values places the mass center at a distance of 0.039 from the origin, which is smaller than the mean value of $0.050\pm0.008$ for the TeV sources in . However, the spread in the points, quantified as the distance of each point from the mass center, is 0.077, which is much higher than the mean values ($0.021\pm0.003$ for TeV and $0.041\pm0.005$ for non-TeV sources) in . In fact, there is only one non-TeV source with a value higher than we obtain for S2 0109+22.

These results are in good agreement with previous studies and indicate that the optical polarization of S2 0109+22 is more variable both in fractional polarization and position angle than other high-energy BL Lac objects, and that there does not seem to be a preferred polarization angle in the source, at least over our monitoring period. This is not unexpected based on the analysis of which showed that the polarization variability depends more on the position of the synchrotron peak rather than the detection of TeV emission . As shown in Figure \[SED\], in the ISP-type S2 0109+22 the optical emission probes the peak of the synchrotron component, where the variability is expected to be higher [see also @2016MNRAS.463.3365A]. Comparing the obtained intrinsic mean polarization fraction to the values presented by @2016MNRAS.463.3365A, this source seems to be a rather typical ISP-type object. The maximum polarization fraction is over 15%, which is high, but not uncommon for ISP sources, as shown in where about 30% of the ISP objects reach fractional polarization values as high as or higher than 15%. This indicates that the magnetic field order must be fairly high in the emission region.

Redshift estimation {#sec3-5}
-------------------

The lack of emission lines in the optical spectrum of BL Lacs objects makes the determination of the redshift of these sources particularly challenging. An estimation on the distance can be obtained from basic assumption on the host galaxy luminosity . Alternatively, an upper limit on the distance can instead be estimated by studying the deformation induced by the EBL on the VHE gamma-ray spectrum.

### Host galaxy {#sec3-5-1}

We use the deep I-band image (see \[HGIMAGE\]) to search for the host galaxy emission. Two-dimensional surface brightness models were fitted to the light distribution of S2 0109+22 in order to study its host galaxy. Prior to the fitting, the background level was measured and subtracted, removing also a small tilt in the background. Two models were considered: 1) a point source (jet) model and 2) a point source + elliptical galaxy model. Both models had three free parameters, point source x-y position plus flux in the first model, and point source flux, host galaxy flux and host galaxy effective radius in the second model. The first model was used to fix the position of the nucleus, i.e. the second model was fit using the position from the first model to fix the point source and the host galaxy into the same position. Moreover, the ellipticity of the host galaxy was fixed to $\varepsilon = 0$ and the Sersić index to $n = 4$. Both models were convolved with the PSF, determined from two nearby stars, located at 61and 84away from, and with similar peak intensity to S2 0109+22. The fit was performed using pixels within 10.5of the center of S2 0109+22.

We used a Metropolis sampler [e.g. @2017arXiv170404629M] to map a posteriori distribution in three-dimensional parameter space. We employed 10 independent walkers, each completing 30 000 iteration steps and with flat priors. The walkers were initially distributed randomly over a fairly wide range of values, but they all quickly converged towards the same area in the parameter space corresponding to the maximum likelihood. The calculation of likelihood assumed that the pixel values had an uncertainty consisting of four components, each normally distributed: 1) Photon noise, 2) readout noise, 3) error in background determination and 4) error in the PSF model. The background uncertainty was determined by measuring the background around the source in 10 rectangular regions. For the PSF error, we subtracted the PSF from a star close to S2 0109+22 and examined the residuals. The residuals were the strongest near the center of the star, where they amounted to 2% of the local signal.

Figure \[hostplot\] shows the marginalized posterior distributions of the two host galaxy parameters: the host galaxy flux and effective radius. The parameters are correlated and in addition both correlate strongly with the point source flux. The best-fit (mode of the posteriors) parameters of model no.2 correspond to AGN flux = (6.651$\pm$0.003) mJy, host galaxy flux (0.149$\pm$0.003) mJy and effective radius (1.40$\pm$0.04). The host galaxy flux in the *I*-band optical is I = 18.05 mag.

If we make the assumption that the host galaxy is a passively evolving early type galaxy with absolute magnitude $M_R = -22.8$ [@2005ApJ...635..173S] with $R - I$ = 0.7 and using $A_I$ = 0.057 for the Galactic absorption [@2011ApJ...737..103S], then we obtain $z = 0.36 \pm 0.07$. This value and its error are a result of 1000 trials where we first drew $M_R$ from a Gaussian distribution with average $-22.8\pm0.5$ and then determined the redshift compatible with the observed I-band magnitude taking into account the evolution, K-correction and Galactic absorption.

### Redshift upper limit based on absorption of VHE gamma rays

The absorption of VHE gamma rays through interaction with the EBL increases with source distance and photon energy. Basic assumptions on the intrinsic spectrum can be used to infer a limit on the distance of the blazar [e.g. @2007ApJ...655L..13M; @2010MNRAS.405L..76P]. In order to determine an upper limit for the source distance, we assumed that the intrinsic spectrum is described by a power law or a concave function (i.e. hardness does not increase with energy). The archival data (Fig. \[SED\]) indicates that the spectrum of the source in the HE gamma-ray band is variable. Considering that the source is not located at $z>1$, we assume the hardest possible spectrum for this redshift as an intrinsic power-law index. As a conservative approach we assume a fixed photon-index limit of 1.5 following @2006Natur.440.1018A and . We obtain a 95% confidence level limit to the S2 0109+22 redshift of $z\leq 0.60$. The value is obtained by means of a maximum likelihood fit to the observed event rates vs. the reconstructed energy, modelling the intrinsic spectrum with a power-law function, using the EBL model of @2011MNRAS.410.2556D, and performing a scan in redshift. The limit is obtained, following @2001NIMPA.458..745R, from the resulting profile likelihood vs. redshift, with the intrinsic source parameters, and the background rates vs. reconstructed energy, treated as nuisance parameters. A more conservative limit can be estimated by varying the simulated total light throughput of the instrument by $\pm15\%$. This yields an 95% upper limit on the redshift of $z\leq 0.67$. To estimate the uncertainties caused by EBL model selection, we test eight different EBL models . The results show that the uncertainties due to EBL model selection are negligible compared with the instrumental uncertainties. Finally, in order to verify the assumed intrinsic photon index (1.5), we compare the results with the ones obtained by assuming the photon index in Section \[sec3-2\] ($\Gamma=1.81 \pm 0.14$). The comparison shows that results are consistent with each other.

The estimated redshift ($z=0.36\pm0.07$) and the calculated redshift 95% upper limit ($z\leq 0.67$) in this paper are consistent with the value reported by @2016MNRAS.458.2836P [$z>0.35$]. Therefore, we used $z=0.35$ based on the accuracy of the technique and other uncertainties, to calculate the intrinsic properties of the source.

![\[hostplot\]Marginalized posterior distributions of the host galaxy flux and effective radius *(bottom-right)*. The likelihood distribution of effective radius *(bottom-left)* and host galaxy flux *(top-right)*. The color is proportional to the probability.](bayesplot_final.eps){width="47.00000%"}

Spectral energy distribution
----------------------------

[\[3-6\]]{} In this section, we first present the physical modelling of the SED based on the quasi-simultaneous data described in previous sections. Then, we use a mathematical approach in order to calculate the location of SED peaks in other epochs using archival data (as well as the quasi-simultaneous data near the flare night).

### Broad band SED {#BBSED}

In Figure \[SED\], we plot the broad-band SED of S2 0109+22 using the multi-frequency data described in previous sections. For the SED modelling, the HE and VHE gamma-ray spectra are constructed from MAGIC and *Fermi*-LAT data obtained on MJD 57225-57231. The VHE gamma-ray data are corrected for the EBL absorption effect using the @2011MNRAS.410.2556D model. The VHE gamma-ray spectrum is dominated by the signal from the flare night. However, for the night of the flare, we do not have enough strictly simultaneous data to produce a robust model. As discussed in Section \[sec3-1\], the low significance of the signal outside the flare night did not allow us to construct a low-state SED. For X-ray, UV and optical, we selected the data points which are near the flare night, to avoid averaging a variable source with different distribution of observation times during the MAGIC campaign. The *Swift*-UVOT and *Swift*-XRT data are used to reproduce the UV and X-ray spectra of the source on MJD 57228.41. The optical data point, obtained by the KVA telescope on MJD 57228.22, is corrected for Galactic extinction. The host galaxy contribution to the optical flux is neglected (see Sect. \[sec3-5-1\]). The radio data points were collected on MJD 57227.41 and 57227.08 in 15 and 37GHz respectively, but are not used for SED modelling (see below).

The quasi-simultaneous SED was modelled using a one-zone synchrotron self-Compton model [@2003ApJ...593..667M]. It assumes a spherical, relativistically moving emission region characterized by its radius $R$, magnetic field B and Doppler factor $\delta$. It contains an electron population following a broken power-law distribution with index $p_1$ for $\gamma_\text{min}<\gamma<\gamma_\text{break}$ and $p_2$ for $\gamma_\text{break}<\gamma<\gamma_\text{max}$. The normalization of this electron distribution at $\gamma$=1 is $K$. We use the redshift of $z=0.35$ for the source (see Sect. \[sec3-5\]).

The goodness of the fitted model is judged by a $\chi^2$-test ($\chi^2$/d.o.f = 22.4/16) assuming fixed $\gamma_\text{min}=1.0 \times 10^{3}$ because there is no instrument available to probe the energy range where the influence of $\gamma_\text{min}$ would be significant. Therefore, the curve represents only one possible set of SED parameters. The other parameters used for the model are: $R=5.5 \times 10^{16}$ cm, $B=0.054$ G, $\delta=21.7$, $\gamma_\text{break}=1.2 \times 10^{4}$, $\gamma_\text{max}=4.5 \times 10^{5}$, $p_1=1.94$, $p_2=3.68$ and $K=3.1 \times 10^3$ cm$^{-3}$. The assumed emission region size is compatible with a daily variability time scale. There is no evidence of a shorter variability time scale in the multi wavelength data during the MAGIC campaign. The parameters are rather typical for TeV BL Lac objects [See e.g. @2010MNRAS.401.1570T].

The one-zone model does not reproduce the spectrum at the lowest frequencies, since the emission is self-absorbed below the millimetre band. It is generally assumed that this emission is produced in the outer regions of the jet. This is in agreement with the results in Section \[3-4-1\], where no connection between the long-term behaviour of the optical and radio bands was found for this source. Moreover, the location of the SED peaks are roughly estimated to be $\log \nu_{\text{sync}}\simeq 15.4$ and $\log \nu_{\text{IC}}\simeq 23.3$.

### SED peaks {#syncpeak}

In order to determine the peak frequencies of the SED components, we fitted simultaneously two log-parabolic spectra , one for the synchrotron peak and another for the Inverse Compton (IC), to the SED of the source. We try to calculate the location of the SED peaks for two different states. First, we extracted the archival data from the ASI Space Science Data Center[^13]. Since the archival data are non-simultaneous and $\nu_{\text{peak}}$ is known to change with the activity state in blazars [e.g. @2009ApJ...705.1624A], we can expect the fitted $\nu_{\text{peak}}$ to depend on the frequencies covered and on the number of observing epochs. To roughly estimate how much this could affect $\nu_{\text{peak}}$ we constructed 4 different samples from the archival data, one representing a high state, another for a low state and two ‘mixed’ states. The archival data indicate that the source is an intermediate synchrotron peak BL Lac object based on the classifications in @2010ApJ...716...30A with $\log \nu_{\text{sync}}=14.4 \pm 0.1$ and $\log \nu_{\text{IC}}=22.8 \pm 0.2$, which is consistent with the source classification reported by .

In the second step we used the quasi-simultaneous data described in section \[BBSED\]. The locations of the peaks are $\log \nu_{\text{sync}}=15.1 \pm 0.5$ and $\log \nu_{\text{IC}}=23.1 \pm 0.2$, which are consistent with the results obtained from the physical modelling described in Section \[BBSED\]. Table \[synchtab\] shows the summary of the SED peaks using different approaches and datasets. Based on the broad-band SED modelled for this dataset, the X-ray emission is purely synchrotron, which is normal for HSP BL Lac objects [e.g. @2010ApJ...708L.100A]. The historical X-ray observations of 2006 (Table \[tab\_xray\]) show a hard X-ray spectral index ($\Gamma_X=2.06\pm 0.05$) which is in good agreement with the broad-band SED reported by @2004MNRAS.348.1379C and the normal case for LSP and ISP BL Lac objects [e.g. @2017ApJ...844...58P]. Therefore, there is a hint of a transition from intermediate to high synchrotron peak during the MAGIC observation period. The transition is not only in the peak but the whole SED is appearing as a typical X-ray bright HSP SED .

Summary {#sec4}
=======

S2 0109+22 was discovered for the first time in the HE gamma-ray band by the *Fermi*-LAT during the first three months of sky-survey operation in 2008 [@2009ApJ...700..597A]. Previous EGRET upper limits are reported in @2004MNRAS.348.1379C.

In this paper, we reported the first VHE gamma-ray detection of S2 0109+22 by MAGIC in 2015 July. The MAGIC observation was triggered by the source high state in HE gamma rays. During the MAGIC observation campaign, the HE gamma-ray LC does not show variability on a daily time scale, while the constant fit to VHE gamma-ray flux was rejected with $3\,\sigma$ level of confidence.

We performed a long-term and a short-term multi-frequency study of the source, from radio to VHE gamma rays and compare the source to other TeV blazars. The summary of the main outcomes are:

1.  Compared to the sample of 21 known variable TeV BL Lac objects (Fig. \[histo2\]), the observed VHE gamma-ray flux from S2 0109+22 is relatively low. The predicted low state VHE gamma-ray flux by is below the sensitivity of the current generation of IACTs. Moreover, the source was not detected after its flaring activity by MAGIC. Therefore, this source will be a good candidate to be monitored by the future Cherenkov Telescope Array (CTA) in order to characterize its VHE gamma-ray temporal behaviour and its connection to lower energy bands.

2.  The brighter-harder trend is clear in the X-ray band (Table \[tab\_xray\] and Fig. \[spec\_flux\]). Similar behaviour was observed for many TeV BL Lac objects [@2017ApJ...841..123P]. However, such a trend is not present in the VHE gamma-ray (Table \[tab\_tev\_spec\] and Fig. \[vhe\_sed\]) and HE gamma-ray (Fig. \[fig\_lc\], panel b and c) bands, but this could be due to large error bars in these bands. The absence of such a correlation in HE and VHE gamma-ray bands for non-HSP BL Lac objects is widely discussed in the context of the ‘blazar sequence’ [see @2015ApJ...810...14A and references therein].

3.  In the long-term optical and radio LC (Fig. \[fig\_longterm\]), there was no correlation peak between 15GHz and optical flux. This fact suggests that, unlike for many other TeV blazars, the optical and radio emission do not originate from the same region or that the correlation is too complex to be probed by the method found in .

4.  The optical fractional polarization and polarization angle of S2 0109+22 are more variable than found for typical high-energy BL Lac objects .

5.  We used two methods to estimate the redshift of the source. The result of the photometric host-galaxy method is $z = 0.36\pm0.07$, whereas the 95% upper limit estimation based on the absorption of VHE gamma-ray emission, assuming the EBL model described in @2011MNRAS.410.2556D, gives $z \leq 0.67$. The estimated redshifts are in agreement with the one derived by @2016MNRAS.458.2836P.

6.  When comparing the quasi-simultaneous SED presented in this paper with archival data obtained from the ASI Space Science Data Center, there is a hint of intermediate to high synchrotron peak transition. This has been previously suggested for PKS 0301-243 and 1ES 1011+496 in 2008 [@2016MNRAS.459.2286A].

7.  The broad-band SED of S2 0109+22 (Sect. \[BBSED\]) reveals that the parameters of a single-zone SSC model are rather typical for TeV BL Lac objects. Comparing the SED parameters with the ones reported in @2004MNRAS.348.1379C reveals that the magnetic field strength is an order of magnitude weaker. Weaker magnetic field energy density ($U_B=B^2/8\pi$) increases the radiation to magnetic energy ratio ($U_{\text{rad}}/U_{B}=L_{\text{IC}}/L_{\text{sync}}$). Therefore, SSC luminosity component increases to the level above the sensitivity of VHE gamma-ray instruments.

The long-term radio to optical and optical polarization behaviour of the source agree with the classification of the source as an ISP BL Lac object, which are still a minority in the class of TeV blazars. However, there is a hint of type transition as discussed in Section \[syncpeak\] based on the multi epoch comparison of the SED peak locations and X-ray behaviour of the source. In order to precisely characterize the source-type transition behaviour, more simultaneous multi-wavelength observations during different flux states are needed. Such observations can be performed when CTA enables us to detect VHE gamma-ray emission also during the low state of the source. Moreover, considering the increased SSC luminosity, high polarization degree in the optical and high X-ray luminosity of the source make this source an ideal candidate for physical modelling when the X-ray and soft-gamma-ray (MeV) polarization observations become available by instruments such as Imaging X-ray Polarimetry Explorer [IXPE, @2016SPIE.9905E..17W], e-ASTROGAM [@2017ExA....44...25D], and All-sky Medium Energy Gamma-ray Observatory [AMEGO, @2017HEAD...1610313M].

Acknowledgments {#acknowledgments .unnumbered}
===============

Part of this work is based on archival data, software or online services provided by the Space Science Data Center - ASI.

We would like to thank the Instituto de Astrofísica de Canarias for the excellent working conditions at the Observatorio del Roque de los Muchachos in La Palma. The financial support of the German BMBF and MPG, the Italian INFN and INAF, the Swiss National Fund SNF, the ERDF under the Spanish MINECO (FPA2015-69818-P, FPA2012-36668, FPA2015-68378-P, FPA2015-69210-C6-2-R, FPA2015-69210-C6-4-R, FPA2015-69210-C6-6-R, AYA2015-71042-P, AYA2016-76012-C3-1-P, ESP2015-71662-C2-2-P, CSD2009-00064), and the Japanese JSPS and MEXT is gratefully acknowledged. This work was also supported by the Spanish Centro de Excelencia “Severo Ochoa” SEV-2012-0234 and SEV-2015-0548, and Unidad de Excelencia “María de Maeztu” MDM-2014-0369, by the Croatian Science Foundation (HrZZ) Project IP-2016-06-9782 and the University of Rijeka Project 13.12.1.3.02, by the DFG Collaborative Research Centers SFB823/C4 and SFB876/C3, the Polish National Research Centre grant UMO-2016/22/M/ST9/00382 and by the Brazilian MCTIC, CNPq and FAPERJ.

The *Fermi*-LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l’Energie Atomique and the Centre National de la Recherche Scientifique / Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d’Études Spatiales in France. This work performed in part under DOE Contract DEAC02-76SF00515.

The OVRO 40-m monitoring program is supported in part by NASA grants NNX08AW31G, NNX11A043G and NNX14AQ89G, and NSF grants AST-0808050 and AST-1109911.

affiliations {#affiliations .unnumbered}
============

$^{1}$ [Università di Udine, and INFN Trieste, I-33100 Udine, Italy]{}\
$^{2}$ [National Institute for Astrophysics (INAF), I-00136 Rome, Italy]{}\
$^{3}$ [Università di Padova and INFN, I-35131 Padova, Italy]{}\
$^{4}$ [Technische Universität Dortmund, D-44221 Dortmund, Germany]{}\
$^{5}$ [Croatian MAGIC Consortium: University of Rijeka, 51000 Rijeka, University of Split - FESB, 21000 Split, University of Zagreb - FER, 10000 Zagreb, University of Osijek, 31000 Osijek and Rudjer Boskovic Institute, 10000 Zagreb, Croatia.]{}\
$^{6}$ [Saha Institute of Nuclear Physics, HBNI, 1/AF Bidhannagar, Salt Lake, Sector-1, Kolkata 700064, India]{}\
$^{7}$ [Max-Planck-Institut für Physik, D-80805 München, Germany]{}\
$^{8}$ [now at Centro Brasileiro de Pesquisas Físicas (CBPF), 22290-180 URCA, Rio de Janeiro (RJ), Brasil]{}\
$^{9}$ [Unidad de Partículas y Cosmología (UPARCOS), Universidad Complutense, E-28040 Madrid, Spain]{}\
$^{10}$ [Inst. de Astrofísica de Canarias, E-38200 La Laguna, and Universidad de La Laguna, Dpto. Astrofísica, E-38206 La Laguna, Tenerife, Spain]{}\
$^{11}$ [University of Łódź, Department of Astrophysics, PL-90236 Łódź, Poland]{}\
$^{12}$ [Deutsches Elektronen-Synchrotron (DESY), D-15738 Zeuthen, Germany]{}\
$^{13}$ [ETH Zurich, CH-8093 Zurich, Switzerland]{}\
$^{14}$ [Institut de Física d’Altes Energies (IFAE), The Barcelona Institute of Science and Technology (BIST), E-08193 Bellaterra (Barcelona), Spain]{}\
$^{15}$ [Università di Siena and INFN Pisa, I-53100 Siena, Italy]{}\
$^{16}$ [Università di Pisa, and INFN Pisa, I-56126 Pisa, Italy]{}\
$^{17}$ [Universität Würzburg, D-97074 Würzburg, Germany]{}\
$^{18}$ [Finnish MAGIC Consortium: Tuorla Observatory and Finnish Centre of Astronomy with ESO (FINCA), University of Turku, Vaisalantie 20, FI-21500 Piikkiö, Astronomy Division, University of Oulu, FIN-90014 University of Oulu, Finland]{}\
$^{19}$ [Departament de Física, and CERES-IEEC, Universitat Autónoma de Barcelona, E-08193 Bellaterra, Spain]{}\
$^{20}$ [Japanese MAGIC Consortium: ICRR, The University of Tokyo, 277-8582 Chiba, Japan; Department of Physics, Kyoto University, 606-8502 Kyoto, Japan; Tokai University, 259-1292 Kanagawa, Japan; RIKEN, 351-0198 Saitama, Japan]{}\
$^{21}$ [Inst. for Nucl. Research and Nucl. Energy, Bulgarian Academy of Sciences, BG-1784 Sofia, Bulgaria]{}\
$^{22}$ [Universitat de Barcelona, ICC, IEEC-UB, E-08028 Barcelona, Spain]{}\
$^{23}$ [Humboldt University of Berlin, Institut für Physik D-12489 Berlin Germany]{}\
$^{24}$ [also at Dipartimento di Fisica, Università di Trieste, I-34127 Trieste, Italy]{}\
$^{25}$ [also at Port d’Informació Científica (PIC) E-08193 Bellaterra (Barcelona) Spain]{}\
$^{26}$ [also at INAF-Trieste and Dept. of Physics & Astronomy, University of Bologna]{},\
$^{27}$ [Agenzia Spaziale Italiana (ASI) Space Science Data Center, I-00133 Roma, Italy)]{},\
$^{28}$ [Istituto Nazionale di Fisica Nucleare (INFN), Sezione di Perugia, I-06123 Perugia, Italy]{},\
$^{29}$ [Departamento de Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile]{},\
$^{30}$ [Tuorla Observatory, University of Turku, Väisäläntie 20, FI-21500 Piikkiö, Finland]{},\
$^{31}$ [Aalto University Metsähovi Radio Observatory, Metsähovintie 114, 02540 Kylmälä, Finland]{},\
$^{32}$ [Aalto University Department of Electronics and Nanoengineering, P.O. BOX 15500, FI-00076 AALTO, Finland]{},\
$^{33}$ [Tartu Observatory, Observatooriumi 1, 61602 Tõravere, Estonia]{},\
$^{34}$ [Owens Valley Radio Observatory, California Institute of Technology, Pasadena, CA 91125, USA]{},

\[lastpage\]

[^1]: The eleven meter telescope (National Radio Astronomy Observatory) observed the source in 1976.

[^2]: <https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/>

[^3]: <https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html>

[^4]: The detection significance for a given source is approximately equal to the square root of the Test Statistic, for a given source.

[^5]: <https://swift.gsfc.nasa.gov/analysis/threads/gen_thread_attfilter.html>

[^6]: <https://heasarc.gsfc.nasa.gov/W3Browse/swift/swiftxrlog.html>

[^7]: The difference between the number of data points measured by *UVOT* and *XRT* is due to the usage of *XRT* window timing mode, multiple *UVOT* snapshots during *XRT* exposure, and bad quality of *XRT* raw images.

[^8]: Calculated based on the value obtained from @2011ApJ...737..103S

[^9]: <http://users.utu.fi/kani/1m>

[^10]: <http://herculesii.astro.berkeley.edu/kait/agn>

[^11]: \[ALFOSC\]<http://www.not.iac.es/instruments/alfosc>

[^12]: <https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/source_models.html#PowerLaw2>

[^13]: <http://www.asdc.asi.it/>
