Astronomy & Astrophysics manuscript no. main˙02˙07˙2018 c© ESO 2018 July 3, 2018 Multi-wavelength characterization of the blazar S5 0716+714 during an unprecedented outburst phase MAGIC Collaboration: M. L. Ahnen1, S. Ansoldi2,20, L. A. Antonelli3, C. Arcaro4,27, D. Baack5, A. Babić6, B. Banerjee7, P. Bangale8, U. Barres de Almeida8,9, J. A. Barrio10, J. Becerra González11, W. Bednarek12, E. Bernardini4,13,23,27, R. Ch. Berse5, A. Berti2,24, W. Bhattacharyya13, A. Biland1, O. Blanch14, G. Bonnoli15, R. Carosi15, A. Carosi3, G. Ceribella8, A. Chatterjee7, S. M. Colak14, P. Colin8, E. Colombo11, J. L. Contreras10, J. Cortina14, S. Covino3, P. Cumani14, P. Da Vela15, F. Dazzi3, A. De Angelis4,27, B. De Lotto2, M. Delfino14,25, J. Delgado14, F. Di Pierro4,27, A. Domı́nguez10, D. Dominis Prester6, D. Dorner16, M. Doro4,27, S. Einecke5, D. Elsaesser5, V. Fallah Ramazani17, A. Fernández-Barral4,14,27, D. Fidalgo10, M. V. Fonseca10, L. Font18, C. Fruck8, D. Galindo19, S. Gallozzi3, R. J. Garcı́a López11, M. Garczarczyk13, M. Gaug18, P. Giammaria3, N. Godinović6, D. Gora13, D. Guberman14, D. Hadasch20, A. Hahn8, T. Hassan14, M. Hayashida20, J. Herrera11, J. Hose8, D. Hrupec6, K. Ishio8, Y. Konno20, H. Kubo20, J. Kushida20, D. Kuveždić6, D. Lelas6, E. Lindfors17?, S. Lombardi3, F. Longo2,24, M. López10, C. Maggio18, P. Majumdar7, M. Makariev21, G. Maneva21, M. Manganaro11 ?, K. Mannheim16, L. Maraschi3, M. Mariotti4,27, M. Martı́nez14, S. Masuda20, D. Mazin8,20, K. Mielke5, M. Minev21, J. M. Miranda15, R. Mirzoyan8, A. Moralejo14, V. Moreno18, E. Moretti8, T. Nagayoshi20, V. Neustroev17, A. Niedzwiecki12, M. Nievas Rosillo10, C. Nigro13, K. Nilsson17, D. Ninci14, K. Nishijima20, K. Noda14, L. Nogués14, S. Paiano4,27, J. Palacio14, D. Paneque8, R. Paoletti15, J. M. Paredes19, G. Pedaletti13?, M. Peresano2, M. Persic2,26, P. G. Prada Moroni22, E. Prandini4,27, I. Puljak6, J. R. Garcia8, I. Reichardt4,27, W. Rhode5, M. Ribó19, J. Rico14, C. Righi3, A. Rugliancich15, T. Saito20, K. Satalecka13, T. Schweizer8, J. Sitarek12,20, I. Šnidarić6, D. Sobczynska12, A. Stamerra3, M. Strzys8, T. Surić6, M. Takahashi20, L. Takalo17, F. Tavecchio3, P. Temnikov21, T. Terzić6, M. Teshima8,20, N. Torres-Albà19, A. Treves2, S. Tsujimoto20, G. Vanzo11, M. Vazquez Acosta11, I. Vovk8, J. E. Ward14, M. Will8, D. Zarić6; from Fermi-LAT Collaboration: D. Bastieri27, D. Gasparrini28, B. Lott29, B. Rani30?, D. J. Thompson30; MWL Collaborators: I. Agudo31, E. Angelakis32, G. A. Borman33, C. Casadio32,31, T. S. Grishina34, M. Gurwell35, T. Hovatta36,37,38, R. Itoh39, E. Järvelä36,37, H. Jermak47, S. Jorstad34,40, E. N. Kopatskaya34, A. Kraus32, T. P. Krichbaum32, N. P. M. Kuin41, A. Lähteenmäki36,37, V. M. Larionov34,42, L. V. Larionova34, A. Y. Lien30, G. Madejski43, A. Marscher40, I. Myserlis32, W. Max-Moerbeck32, S. N. Molina31, D. A. Morozova34, K. Nalewajko44, T. J. Pearson45, V. Ramakrishnan36, A. C. S. Readhead45, R. A. Reeves46, S. S. Savchenko34, I. A. Steele47, M. Tornikoski36, Yu. V. Troitskaya34, I. Troitsky34, A. A. Vasilyev34, and J. Anton Zensus32 (Affiliations can be found after the references) ABSTRACT Context. The BL Lac object S5 0716+714, a highly variable blazar, underwent an impressive outburst in January 2015 (Phase A), followed by minor activity in February (Phase B). The MAGIC observations were triggered by the optical flux observed in Phase A, corresponding to the brightest ever reported state of the source in the R-band. Aims. The comprehensive dataset collected is investigated in order to shed light on the mechanism of the broadband emission. Methods. Multi-wavelength light curves have been studied together with the broadband Spectral Energy Distributions (SEDs). The data set collected spans from radio (Effelsberg, OVRO, Metsähovi, VLBI, CARMA, IRAM, SMA), UV (Swift-UVOT), optical photometry and polarimetry (Tuorla, Steward, RINGO3, KANATA, AZT-8+ST7, Perkins, LX-200), X-ray (Swift-XRT and NuSTAR), high-energy (HE, 0.1 GeV ¡ E ¡ 100 GeV) with Fermi-LAT to the very-high-energy (VHE, E¿100 GeV) with MAGIC. Results. The flaring state of Phase A was detected in all the energy bands, providing for the first time a multi-wavelength sample of simultaneous data from the radio band to the VHE. In the constructed SED the Swift-XRT+NuSTAR data constrain the transition between the synchrotron and inverse Compton components very accurately, while the second peak is constrained from 0.1 GeV to 600 GeV by Fermi+MAGIC data. The broadband SED cannot be described with a one-zone synchrotron self-Compton model as it severely underestimates the optical flux in order to reproduce the X-ray to γ-ray data. Instead we use a two-zone model. The electric vector position angle (EVPA) shows an unprecedented fast rotation. An estimation of the redshift of the source by combined HE and VHE data provides a value of z = 0.31 ± 0.02stats ± 0.05sys, confirming the literature value. Conclusions. The data show the VHE emission originating in the entrance and exit of a superluminal knot in and out a recollimation shock in the inner jet. A shock-shock interaction in the jet seems responsible for the observed flares and EVPA swing. This scenario is also consistent with the SED modelling. Key words. galaxies: active – BL Lacertae objects: individual: S5 0716+714 – galaxies:jets – gamma-rays:galaxies 1 ar X iv :1 80 7. 00 41 3v 1 [ as tr o- ph .H E ] 1 J ul 2 01 8 1. Introduction The blazar S5 0716+714 is a BL Lac object characterized by extreme variability in almost all energy bands. Because of the featureless optical continuum (Paiano et al., 2017) it is hard to estimate its redshift. Nilsson et al. (2008) claimed a value of z = 0.31 ± 0.08 based on the photometric detection of the host galaxy. Detection of intervening Lyα systems in the ultra-violet spectrum of the source confirms the earlier estimates with a red- shift value z < 0.32 (95% confidence) (Danforth et al., 2013). The familiar shape of the Spectral Energy Distribution (SED) of blazars, in which the two bumps are identified as synchrotron and inverse Compton (IC) respectively, is used for their classification. In the case of S5 0716+714 the first peak of the SED is located between 1014 and 1015 Hz, leading to a classification as intermediate synchrotron peaked blazar or IBL (Intermediate-peaked BL Lac object) (Giommi et al., 1999; Ackermann et al., 2011). Because of its remarkable variability, S5 0716+714 has been the subject of many optical monitoring campaigns (Wagner et al., 1996; Montagni et al., 2006; Rani et al., 2011, 2013, 2015). Several authors carried out flux variability studies (e.g. Quirrenbach et al., 1991; Wagner et al., 1996) and morpho- logical/kinematic studies at radio frequencies (Antonucci et al., 1986; Witzel et al., 1988; Jorstad et al., 2001; Bach et al., 2005; Rastorgueva et al., 2011). The observed intraday variability at radio wavelengths is likely to be a mixture of intrinsic and exter- nal (due to interstellar scintillation) mechanisms. Wagner et al. (1996) reported a significant correlation between optical – radio flux variations at day-to-day timescales. Rani et al. (2010) re- ported the detection of a ∼ 15-min quasi-periodic oscillations (QPOs) at optical frequencies, so far the shortest ones ever ob- served in any blazar. The study of the broadband flux variability is a difficult task because of the complexity of the flaring activity, which can vary rapidly on a timescale of a few hours to days, while on a timescale of ∼ 1 year a broader and slower variability trend has been reported (Raiteri et al., 2003; Rani et al., 2013). As seen by Very Long Baseline Interferometry (VLBI) studies (Bach et al., 2005; Britzen et al., 2009), the source presents a core-dominated jet pointing towards the north. Very Large Array observations show a halo-like jet misaligned with it by ∼90◦ on kiloparsec scales (Antonucci et al., 1986): Britzen et al. (2009) suggest that there is an apparent stationarity of jet components relative to the core. However, recently apparent speeds of ∼ 40 c have been reported (Rastorgueva et al., 2011; Larionov et al., 2013; Lister et al., 2013; Rani et al., 2015). (Britzen et al., 2009; Rastorgueva et al., 2011; Rani et al., 2014) observed non-radial motion and wiggling component trajectories in the inner milliarcsecond jet region. Observations by BeppoSAX (Tagliaferri et al., 2003) and XMM-Newton (Foschini et al., 2006) provide evidence for a con- cave X-ray spectrum in the 0.1 – 10 keV band, a signature of the presence of both the steep tail of the synchrotron emission and the rising part of IC spectrum. EGRET on board the Compton Gamma-ray Observatory (CGRO) detected high-energy (HE, 0.1 GeV ¡ E ¡ 100 GeV) γ-ray emission from S5 0716+714 several times from 1991 to 1996 (Lin et al., 1995; Hartman et al., 1999). Two strong γ- ray flares in September and October 2007 were detected in the source with AGILE (Chen et al., 2008) in the HE range. These ? Corresponding authors: M. Manganaro, email: manganaro@iac.es; B. Rani, NPP fellow, email:bindu.rani@nasa.gov; G. Pedaletti, email: giovanna.pedaletti@desy.de; E. Lindfors, email: elilin@utu.fi authors also carried out SED modeling of the source with two synchrotron self-Compton (SSC) emitting components, repre- sentative of slowly and rapidly variable components, respec- tively. S5 0716+714 was first detected in the very-high-energy (VHE, E¿100 GeV) range by MAGIC with 5.8σ significance level in November 2007 and then in April 2008 (Anderhub et al., 2009) during an optical flare. At that time MAGIC was work- ing with a single telescope and the energy threshold of the tele- scope for such an high zenith range (47◦ < zd < 55◦) was 400 GeV. The analysis of multi-wavelength data suggested a corre- lation between the VHE γ-ray and optical emission. A struc- tured jet model, composed of a fast spine surrounded by a slower moving layer (Ghisellini et al., 2005; Tavecchio & Ghisellini, 2009) better described the data compared to a simple one-zone SSC model. This source is also among the bright blazars in the Fermi-LAT (Large Area Telescope) Bright AGN Sample (LBAS) (Abdo et al., 2010) and in the Fermi third catalog (Acero et al., 2015) it is among the ones with the highest variability in- dex. The combined GeV – TeV spectrum of the source might display an absorption-like feature in the 10 – 100 GeV energy range (Senturk et al., 2013). Apart from the γ-ray/optical connection, which is to be ex- pected in one-zone, single population leptonic models, another interesting feature of the broadband activity was reported in Rani et al. (2013), in which the authors found the γ-ray/optical flux variations leading the radio variability by ∼ 65 days. An orphan X-ray flare was detected in 2009 by Swift-XRT but no VHE ob- servations were available to check if there was a counterpart in the VHE range. The behavior of the source in the past years can be contextualized in the scenario of a shockwave propagating along a helical path in the blazar’s jet (Marscher et al., 2008; Larionov et al., 2013). Recently, more attention was given to the electric vector po- sition angle swings: Rani et al. (2014) found a significant cor- relation between the inner jet outflow orientation and γ-ray flux variations, showing how the morphology of the inner jet has a strong connection with the γ-ray flares. Chandra et al. (2015a) studied the rapid variation in the degree of polarization (PD) and in the polarization angle (PA) within the Helical Magnetic Field Model (HMFM) (Zhang et al., 2014), explaining such features as most likely due to reconnections in the emission region of the jet. In the present work, we report the results of a multi- wavelength (MWL) campaign organised to follow an unprece- dented outburst phase of the blazar S5 0716+714 during January 2015. The source was detected at its historic highest brightness at optical and IR bands. On January 11, 2015, (MJD 57033) the NIR photometry reported an increase of its flux by a fac- tor of 2.5 in the NIR band in a rather short lapse of 12 days (Carrasco et al., 2015; Chandra et al., 2015b). During the night of 18 January 2015 (MJD 57040), the source was detected at its historic high brightness, with R band magnitude ∼11.71 (Bachev et al., 2015). The TeV observations triggered to follow the excep- tionally high optical state detected the source at energies above 150 GeV (Mirzoyan et al., 2015) and went on until the flaring activity faded away. The paper is structured as follows. In Section 2 and subsec- tions we describe the various instruments involved in the obser- vations and their data analysis. Then in Section 3 we present the multi-wavelength light curves and discuss the peculiarities of the source activity in the different bands. Section 4 is devoted to the jet analysis of the object by VLBI. In Section 5 the spectral fit- ting of NuSTAR and Swift-XRT data are shown. In Section 6 the M. L. Ahnen et al.: MWL characterization of S5 0716+714 during an unprecedented outburst VHE spectra obtained by MAGIC and the redshift estimation using MAGIC and Fermi-LAT simultaneous data are presented, and in Section 7 the broadband SED for the two phases con- sidered is discussed together with the modeling of the source. Results and Conclusions are given in Section 8. 2. Observations and analysis 2.1. VHE γ-ray observations MAGIC is a stereoscopic system consisting of two 17 m diam- eter Imaging Atmospheric Cherenkov Telescopes located at the Observatorio del Roque de los Muchachos, on the Canary Island of La Palma. The current sensitivity for medium-zenith obser- vations (30◦ < zd < 45◦) above 210 GeV is 0.76 ± 0.04 % of the Crab Nebula’s flux in 50 h (Aleksić et al., 2016). On 19 January 2015 (MJD 57041), triggered by the high opti- cal state and by high-energy photons detected by Fermi-LAT, MAGIC started to observe S5 0716+714. On 23 January 2015 (MJD 57045) the significance of the signal was found to be 6.4 σ, and it reached a maximum value of 13.2 σ on 26 January (MJD 57048): the flux increased from (4.1±1.1)×10−11cm−2 s−1 to (8.9 ± 1.1) × 10−11cm−2 s−1 above 150 GeV, which is the highest level ever detected in the VHE band for this source. The next activity of S5 0716+714 in the VHE range was de- tected by MAGIC on 13 February (MJD 57066), this time last- ing four days only, up to the 16 February (MJD 57069). In the present work the two multi-wavelength periods of observations that include MAGIC data are indicated as Phase A, from 18 to 27 January 2015 (MJD 57040 to MJD 57050), and Phase B from 12 to 17 February 2015 (MJD 57065 to MJD 57070). We collected 17.74 hours of data in the zenith angle range of 40◦ < zd < 50◦ and the analysis was performed using the standard MAGIC anal- ysis framework MARS (Zanin et al., 2013; Aleksić et al., 2016). After the applied quality cuts, the survived events amount to 17.46 hours in total. A statistical significance of 18.9 σ was found for the full sample after cuts. The significance of the signal was calculated as in Eq. 17 in Li & Ma (1983). For the MAGIC SED, the systematic uncertainty on the flux normalization was estimated to be 11%, and on the spectral slope ± 0.15. The anal- ysis energy threshold is ∼ 125 GeV, measured as the peak of the Monte Carlo (MC) energy distribution for a source with the spectral shape of S5 0716+714. A full description of the MAGIC systematic uncertainties can be found in Aleksić et al. (2016). 2.2. HE γ-ray observations The HE γ-ray (0.1–300 GeV) observations for a time period between 1 November 2014 (MJD 56962) and 31 July 2015 (MJD 57234) were obtained in a survey mode by the Fermi- LAT (Atwood et al. , 2009). The LAT data were analyzed using the standard ScienceTools1 (software version v10.01.01, pass8) with instrument response function P8R2 SOURCE V6. Photons in the source event class were selected for the analysis. We an- alyzed a region of interest (ROI) of 10◦ in radius centered at the position of S5 0716+714 using a maximum-likelihood algo- rithm (Mattox et al., 1996). We included all 54 sources of the 3FGL catalog (Acero et al., 2015) within 20◦ of the position of S5 0716+714 in the unbinned likelihood analysis. Model param- eters for sources within 5◦ of the ROI are kept free; we kept the model parameters for the rest fixed to their catalogue values. 1 https://fermi.gsfc.nasa.gov/ssc/data/analysis/ To investigate the source variability at E ¿100 MeV, we gen- erated the daily-binned photon flux and index curves using un- binned likelihood analysis. The daily binned data were com- puted by modeling the spectra by a power-law model (PWL, N(E) = N0 E−Γ, N0 : prefactor, and Γ : power law index). We ex- amined the spectral behavior over the whole energy range with a PWL model fitting over equally spaced logarithmic energy bins with Γ kept constant and equal to the value fit over the whole range. Even if in the Fermi third catalog (Acero et al., 2015) the source spectrum is described by a log parabola model, for the present dedicated analysis the PWL model better fits the data. The upper limits for the light curve are shown as grey triangles in the second panel from top in Fig. 1 and they were calculated for test statistics ¡9. 2.3. X-ray and optical/UV observations 2.3.1. NuSTAR The exceptional flare from the object triggered a Target of Opportunity observation of the object by NuSTAR. NuSTAR, a NASA Small Explorer satellite sensitive in the hard X-ray band features two multilayer-coated telescopes, focusing the reflected X-rays on the pixellated CdZnTe focal plane modules with the half-power diameter of an image of a point source of ∼ 1′. It provides a bandpass of 3-79 keV with spectral resolution of ∼ 1 keV. For more details, see Harrison et al. (2013). After screening for the South Atlantic Anomaly passages and Earth occultation, the pointing (with OBSID 90002003002) re- sulted in roughly 14.9 ks of net observing time on 24 January 2015 (MJD 57046). After processing the the raw data with the NuSTAR Data Analysis Software (NuSTARDAS) package v.1.3.1 (via the script nupipeline), the source data were ex- tracted from a region of 45′′ radius, centered on the centroid of X-ray emission, while the background was extracted from a 1.5′ radius region roughly 5′ SW of the source location. Spectra were binned in order to have at least 30 counts per rebinned channel. We considered the spectral channels corresponding nominally to the 3-60 keV energy range, where the source was robustly de- tected. The mean net (background-subtracted) count rates were 0.174±0.003 and 0.165±0.003 cts s−1, respectively, for the mod- ules FPMA and FPMB. We found no variability of the source as a function of time within the NuSTAR observation, and we summed the data into one spectral file for each focal plane mod- ule. 2.3.2. Swift-XRT The multi epochs (35) event-list obtained by the X-ray Telescope (XRT) (Burrows et al. , 2004), on-board the Neil Gehrels Swift satellite in the period 1 January 2015 (MJD 57023.2) to 28 February 2015 (MJD 57081.2) with total exposure time of ∼16.81 hours were downloaded from publicly available SWIFTXRLOG (Swift XRT Instrument Log). They were pro- cessed using the procedure described in Ramazani et al. (2017). All these observations have been performed in photon count- ing (PC) mode, with an average integration time of 1.7 ks each. The equivalent Galactic hydrogen column density is fixed to the value of nH = 3.11 × 1020cm−2 (Kalberla et al., 2005). We per- formed spectral fits to all 35 epochs using a simple power law model, with Galactic absorption. This model provides a good fit, and the 0.3 -10 keV spectral index is 2.72 ± 0.1. We discuss the Swift observations more extensively in Section 3. We note specifically that one of the pointings was simultaneous with the 3 M. L. Ahnen et al.: MWL characterization of S5 0716+714 during an unprecedented outburst NuSTAR observation (at MJD 57045–24 January 2015) and we use that observation for joint XRT - NuSTAR spectral fitting in Section 5. 2.3.3. Swift-UVOT Photometric observations by the Swift Ultra-Violet and Optical Telescope (UVOT) instrument were made in the three UV (uvw2, uvm2, and uvw1) and three optical (u, b, and v) filters in both imaging and event mode. During the February 2015 ob- servations, good coverage in the UV and u bands exist and in particular uvm2 data in event mode were obtained in order to resolve short timescale variability in the UV. The UVOT data reduction used the Heasarc Heasoft version 6.16 and Swift CALDB (September 2013). The event mode data were typically 800 s long exposures and were binned into shorter time slices, converted into images, then aspect corrected. This analysis allowed the identification of a few observations for which the pointing drift was too large, and those were ex- cluded from further analysis. Data taken in image mode were similarly validated. The magnitudes were determined from the images using the UVOTMAGHIST program using the standard calibration (Poole et al., 2008; Breeveld et al., 2011). The de- tails of the event mode data processing are as follows: GTI ex- tensions were created in the event file for the desired time in- tervals; UVOTATTJUMPCORR was run to improve the attitude file; COORDINATOR and UVOTSCREEN were used to correct the event file, which was then processed using UVOTIMAGE so the individual data could be inspected. 2.4. Optical observations and polarimetry data The Tuorla Blazar monitoring program2 collects blazar opti- cal light curves in the R band from several observatories. The present work shows in particular data from the 1.03 m telescope at Tuorla Observatory, Finland and the 35 cm telescope at the KVA observatory on La Palma, Canary Islands, Spain. The data are analyzed with a semi-automatic pipeline using standard pro- cedures (Nilsson et al., 2017) The Boston University (BU) group uses the 1.83 m Perkins Telescope at Lowell Observatory (Flagstaff, AZ) to carry out optical observations of a sample of γ-ray blazars, including S5 0716+714. The telescope is equipped with the PRISM cam- era operating in photometric (UBVRI) and polarimetric modes. The details of the observations and data reduction can be found in Jorstad et al. (2010). The optical polarimetric data presented here are comple- mented by data from the RINGO3 polarimeter on the 2.0 m fully robotic Liverpool Telescope on the Roque the los Muchachos Observatory (La Palma, Canary Islands, Spain). The RINGO3 data presented here were obtained as part of a monitor- ing program of bright optical blazars (Jermak et al., 2016). S5 0716+714 is one of the targets regularly monitored by this program, with a time cadence of ∼ 3 days. The RINGO3 po- larimeter acquires polarimetric measurements in three different passbands recorded in the so called ”Red”, ”Green”, and ”Blue” cameras3. Here we only present the data from the ”Green” cam- era (the one with the closest wavelength passband to R-band). The RINGO3 data were reduced following the procedure ex- plained in Steele et al. (2017). 2 http://users.utu.fi/kani/1m 3 See the RINGO3 specifications at: http://telescope.livjm.ac.uk/TelInst/Inst/RINGO3/ R-band photometry and polarimetry observations of S5 0716+714 were performed using the HONIR (Hiroshima Optical and Near-InfraRed camera) instrument installed on the 1.5 m Kanata telescope located at the Higashi-Hiroshima Observatory, Japan (Akitaya et al., 2014). A sequence of pho- topolarimetric observations consisted of successive exposures at 4 position angles of a half-wave plate: 0, 45, 22.5 and 67.5 deg. The data were reduced under the standard procedure of CCD photometry. The aperture photometry was performed using the APPHOT package in PYRAF4, and the differential photometry with a comparison star taken in the same frame of S5 0716+714. The comparison star is located at R.A. = 07:21:53.44 and Decl. = +71:20:36.0 (J2000), and its magnitude is R = 14.032 (UCAC-4 Catalog). The polarization angle is defined as usual (measured from north to east), based on calibrations with polarized stars, HD183143 and HD204827 (Schulz & Lenzen, 1983). The systematic error caused by instrumental polarization was smaller than 2 deg using the polarized stars. St. Petersburg University optical photometric and polarimet- ric data are from the 70 cm AZT-8 telescope of the Crimean Astrophysical Observatory5 and the 40 cm telescope LX-200 in St. Petersburg, both equipped with nearly identical imaging photometers-polarimeters. Polarimetric observations were per- formed using two Savart plates rotated by 45◦ relative to each other (see Larionov et al., 2008b). Instrumental polarization was found via stars located near the object under the assumption that their radiation is unpolarized. This is indicated also by the low level of extinction in the direction of S5 0716+714 (AV = 0.085 mag; AR = 0.067 mag; Schlafly & Finkbeiner 2011). 2.5. Radio observations The 230 GHz (1.33 mm) flux density data were obtained at the Submillimeter Array (SMA) near the summit of Mauna Kea (Hawaii). S5 0716+714 is included in an ongoing monitoring program at the SMA to determine the flux densities of compact extragalactic radio sources that can be used as complex gain cal- ibrators at mm wavelengths (Gurwell et al. , 2007). Potential calibrators are from time to time observed for 3 to 5 minutes. Data from this program are updated regularly and are available at the SMA website6, while the present analysis was a dedicated one. The IRAM 30 m millimeter Radiotelescope provided 230 GHz (1.3 mm) and 86 GHz (3.5 mm) data that were ob- tained as part of the POLAMI7 (Polarimetric AGN Monitoring at Millimeter Wavelengths) program, see Agudo et al. (2018a,b) and Thum et al. (2018). The POLAMI data of S5 0716+714 were acquired and reduced as described in detail in Agudo et al. (2018a). The CARMA data were taken with the eight 3.5 m an- tennas as part of the Monitoring of γ-ray Active galactic nu- clei with Radio, Millimeter and Optical Telescopes (MARMOT) project8. We used 7.5 GHz of bandwidth with a center frequency of 94 GHz. The integration time on S5 0716+714 was 5 minutes for each observation, which yields a typical rms of 10−110 mJy. Absolute flux density calibration was done using nearby obser- vations of Mars. The observational errors are dominated by the absolute calibration uncertainty, assumed to be 10%. All data were processed using the Multichannel Image Reconstruction 4 http://www.stsci.edu/institute/software hardware/pyraf/ 5 http://craocrimea.ru/ru 6 http://sma1.sma.hawaii.edu/callist/callist.html 7 http://polami.iaa.es 8 http://www.astro.caltech.edu/marmot/ 4 http://telescope.livjm.ac.uk/TelInst/Inst/RINGO3/ M. L. Ahnen et al.: MWL characterization of S5 0716+714 during an unprecedented outburst Image Analysis and Display (MIRIAD; Sault, Teuben & Wright, 1995). We have analyzed Very Long Baseline Array (VLBA) data obtained for S5 0716+714, within the VLBA-BU-Blazar pro- gram9, which are contemporaneous with the high energy event in January and February 2015. The data include total and polarized intensity images at 43 GHz at 9 epochs from November 2014 to August 2015. They were reduced using the Astronomical Image Processing System (AIPS)10 and Difmap11 software packages, in the general manner described by Jorstad et al. (2017, 2005). The total intensity images were modeled by components with circular Gaussian brightness distributions. This allows us to de- termine the minimum number of components that provides the best fit between the data and model at each epoch, as well as the following parameters of components: flux density, S , distance from the core, r, position angle with respect to the core, Θ, and size of the component, a (FWHM of the Gaussian). These pa- rameters are given in Table A.1 of Appendix A. The 37 GHz observations were made with the Metsähovi ra- dio telescope. The measurements were made with a 1 GHz-band dual beam receiver centered at 36.8 GHz. The observations are ON–ON observations, alternating the source and the sky in each feed horn. A typical integration time to obtain one flux density data point is between 1200 and 1800 s. Data points with a signal- to-noise ratio ¡ 4 are handled as non-detections. The flux density scale is set by observations of DR 21. A detailed description of the data reduction and analysis is given in Teräsranta et al. (1998). The error estimate in the flux density includes the con- tribution from the measurement rms and the uncertainty of the absolute calibration. S5 0716+714 was observed at 15 GHz as part of a high- cadence γ-ray blazar monitoring program using the Owens Valley Radio Observatory (OVRO) 40 m telescope (Richards et al., 2011). The OVRO 40 m uses off-axis dual-beam optics and a cryogenic pseudo-correlation receiver with a 15.0 GHz center frequency and 3 GHz bandwidth. The source is alter- nated between the two beams in an ON-ON fashion to remove atmospheric and ground contamination. The fast gain variations are corrected using a 180 degree phase switch. Calibration is achieved using a temperature-stable diode noise source to re- move receiver gain drifts, and the flux density scale is derived assuming the value of 3.44 Jy at 15.0 GHz in Baars et al. (1997). The systematic uncertainty of about 5% in the flux density scale is not included in the error bars. Complete details of the re- duction and calibration procedure are found in Richards et al. (2011). Observations at 2.6, 4.8, 10, and 15 GHz radio bands were conducted using the Effelsberg 100 m radio telescope12. Measurements for the target source and for the calibrator sources were made quasi-simultaneously using the cross-scan method slewing over the source position, in azimuth and elevation di- rection in order to gain the desired sensitivity. Subsequently, at- mospheric opacity correction, pointing off-set correction, gain correction, and sensitivity correction were applied to the data. Details of the observations and data reduction are referred to Angelakis et al. (2015). The sources 3C 286, NGC 7027 and 3C 84 have been used as common calibrators for the instruments listed in this section. 9 www.bu.edu/blazars/VLBAproject.html 10 http://www.aips.nrao.edu/index.shtml 11 ftp://ftp.eso.org/scisoft/scisoft4/sources/difmap/difmap.html 12 http://www.mpifr-bonn.mpg.de/en/effelsberg 3. Multi-wavelength light curves In Fig. 1 we present the multi-wavelength data collected during the course of the campaign. A summary of the most important dates can be found in Table 1. The top panel shows the daily- binned MAGIC light curve: in the VHE band, the no-variability hypothesis has been discarded, since the fit for a constant flux resulted in a χ2/n.d. f . = 42/7 for Phase A and χ2/n.d. f . = 10/3 for Phase B data. A Gaussian function better fits the flare shape, providing χ2/n.d. f . = 15/5 and χ2/n.d. f . = 1.8/1 for the two phases respectively. Phase A peaks on 25 January, MJD 57047.3 ±0.4 (in Fig. 1 corresponding to the P2 vertical dashed line) with a flux of (5.9 ± 0.5) × 10−11cm−2 s−1; the standard deviation of the fit is σ = (2.8 ± 0.5) days. Phase B peaks on 14 February 2015 (MJD 57067.9 ±0.2) with a corresponding flux of (5.3 ± 0.7) × 10−11cm−2 s−1 indicated by the vertical dashed line P4. The standard deviation of the fit is σ = (1.2±0.3) days. Intra-day variability (in shorter time scales with respect to the daily-binned light curves) was not detected with MAGIC: the light curve was fit at different time intervals down to 5 minutes, but a constant fit was found to be consistent with the data up to the daily scale, where variability is significant. The Fermi-LAT daily-binned light curve is shown in the sec- ond panel from top: the first peak visible in the curve, marked with the vertical dashed line P1, is the precursor of the whole flaring activity. After P1, (MJD 57038.5, 16 January 2017), which triggered VHE observations, two other peaks are visible in Phase A, and the maximum flux reached in HE is (8.8±1.4)× 10−7cm−2 s−1, 4 times the average flux in HE for the source from the Fermi-LAT 3FGL Catalog (Acero et al., 2015). The photon index for Fermi-LAT observations stays very close to the average Fermi-LAT index of the source, indicated in the corresponding panel of Fig. 1 with the horizontal dashed blue line. The average integral photon X-ray flux (0.3–10 keV) re- ported by Swift-XRT in Fig. 1, fourth panel from top, is (1.87 ± 0.46) × 10−11erg cm−2 s−1. The X-ray flux peaking at MJD 57047 (25 January 2015) with F(0.3−10keV) = (3.95±0.12)× 10−11erg cm−2 s−1 which is a factor of ∼ 2 higher than the av- erage flux of the analysed period. The X-ray spectral index dur- ing the analysed period varies between 2.03 ≤ ΓX ≤ 2.56. The softest spectral index was obtained on the night of the highest X-ray and VHE γ-ray flux (MJD 57047–25 January 2015). The X-ray spectra start hardening smoothly afterwards for the next 14 consecutive nights. The X-ray spectra on the nights of VHE γ-ray flares of Phase A and B can be well described by a power law with spectral index of ΓX,MJD 57047 = 2.56 ± 0.06 (reduced χ2/n.d.f. = 0.738/29 ) and ΓX,MJD 57067 = 2.33 ± 0.06 (reduced χ2/n.d.f. = 0.871/26) respectively. From the Tuorla optical monitoring, on 18 January 2015 (MJD 57040) the magnitude in the R-band reached a value of ∼ 11.5, higher than the magnitude value of the source in the quiescent state (Rmag ∼ 13.5). Later, the source faded down to Rmag = 12.18 ± 0.03 on 20 January 2015 (MJD 57042), flared up again to Rmag = 11.77 ± 0.02 on 23 January (MJD 57045), and then demonstrated rapid flux variations above Rmag = 11.9 (Spiridonova et al., 2015). In general, the optical light curves and the Swift-UVOT light curves show the same trend, a double peaked shape with the second peak coincident with the dashed vertical line P2 in Fig. 1. P2 is in fact identifying a peak not only in the VHE light curve, but also in the X-ray, optical, and UV bands. AZT-8+ST7 light curve can be fit by a Gaussian peak- ing on MJD 57047 (25 January 2015) at a flux of (67.49 ± 0.01) mJy. The standard deviation of the fit is σ = 3.366± 0.001 days. The Swift-UVOT light curve can be fit by a Gaussian peaking on 5 M. L. Ahnen et al.: MWL characterization of S5 0716+714 during an unprecedented outburst 57010 57020 57030 57040 57050 57060 57070 ] -1 s -2 cm -1 1 [1 0 F 0 5 10 P1 P2 P3 MAGIC >150 GeV 57010 57020 57030 57040 57050 57060 57070 ] -1 s -2 cm -6 [1 0 F 0 0.5 1 6− 10× Fermi-LAT 0.1-100 GeV 57010 57020 57030 57040 57050 57060 57070 in de x F er m i-L A T 1.5 2 2.5 3 0.1-100 GeV 57010 57020 57030 57040 57050 57060 57070 ] -1 s -2 er g cm -1 2 [1 0 F 10 20 30 40 12−10× SWIFT-XRT 0.3-10 keV 57010 57020 57030 57040 57050 57060 57070 [m Jy ] (R -b an d) ν F 20 40 60 80 Tuorla Perkins AZT-8+ST7 LX-200 Kanata Tuorla Perkins AZT-8+ST7 LX-200 Kanata Tuorla Perkins AZT-8+ST7 LX-200 Kanata Tuorla Perkins AZT-8+ST7 LX-200 Kanata 57010 57020 57030 57040 57050 57060 57070 [m Jy ] (U -V -B -b an ds ) ν F 0 20 40 60 UVOT B UVOT U UVOT W1 UVOT W2 UVOT M2 57010 57020 57030 57040 57050 57060 57070 [J y](M ill im et re ) ν F 2 4 6 SMA 230GHz IRAM 230 GHz CARMA 90 GHz IRAM 86GHz hovi 37 GHzaMets 57010 57020 57030 57040 57050 57060 57070 [J y](R ad io ) ν F 1.5 2 2.5 3 3.5 OVRO 15 GHz Effelsberg 15GHz Effelsberg 10GHz Effelsberg 4.8GHz 57010 57020 57030 57040 57050 57060 57070 P ol ar iz at io n [% ] 0 10 20 Kanata Perkins AZT-8+ST7 LX-200 Steward RINGO3 57010 57020 57030 57040 57050 57060 57070 ]° E V P A [ 0 200 400 Kanata Perkins AZT-8+ST7 LX-200 Steward RINGO3 MJD 57010 57020 57030 57040 57050 57060 57070 31-12-2014 14-01-2015 28-01-2015 11-02-2015 Phase A Phase B Fig. 1: Multi-wavelength flux and index curves of S5 0716+714 during the period from MJD 57010 to MJD 57080 (19 December 2014 to 27 February 2015). The shadowed areas indicate Phase A (from 18 to 27 January 2015 – MJD 57040 to MJD 57050) and Phase B (from 12 to 17 February 2015 – MJD 57065 to MJD 57070) high states in the VHE range and the corresponding activity in the other bands. P1, P2 and P3 (vertical dashed lines) indicate peaks in the HE and VHE emission. MJD 57047 (25 January 2015) at a flux of (50.2± 0.1) mJy. The standard deviation of the fit is σ = 3.36 ± 0.02 days. The activity in the radio band shown in Fig. 1 in Phase A and Phase B is moderate compared to the other energy bands but 6 M. L. Ahnen et al.: MWL characterization of S5 0716+714 during an unprecedented outburst Table 1: Summary of important dates MJD Calendar date Description 57038.5 16 January 2015 P1: first peak of the HE emission→ trigger VHE observations 57040 18 January 2015 start of Phase A 57044/45 22/23 January 2015 1day EVPA rotation of ∼360◦ 57047.3 ± 0.53 25 January 2015 P2: Gaussian fit peak of the VHE emission in PHASE A 57050 28 January 2015 end of Phase A 57050 ± 3 28 January 2015 K14b passage through A1 57056 03 February 2015 R4: Gaussian fit peak of radio emission in the intermediate phase 57065 12 February 2015 start of Phase B 57067.8 ± 0.23 14 February 2015 P3: Gaussian fit peak of the VHE emission in Phase B 57070 16 February 2015 end of Phase B 57092 11 March 2015 R5: Gaussian fit peak of radio emission Energy [GeV] 14−10 10−10 6−10 2−10 210 610 ] -1 s -2 d N /d E [T eV c m 2 E 15−10 14−10 13−10 12−10 11−10 10−10 9−10 va r F 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Phase AvarF Phase BvarF MWL SED 26 Jan Fig. 2: Fractional variability (Fvar) as a function of the energy for Phase A and Phase B. Vertical bars denote 1σ uncertainties. In gray a Spectral Energy Distribution of the source is overlayed (a snapshot of 26 January 2015), to make easier to associate the value of the fractional variability to the corresponding energy band. does not describe a simply quiescent radio state. In the past the source has gone through many high states in the radio band, for instance the one described in Rani et al. (2013), Fig. 2, named R6. During the R6 flare, the highest flux density reported was ∼ 10 Jy while in the present one, the highest level of radio flux density at the same frequency (230 GHz) is only 5 Jy. We study in more detail the possible delay between radio and optical/γ-ray bands in Sec. 3.1. As reported in Chandra et al. (2015a), the Phase A flare presents a double peaked shape in the HE γ-ray and optical bands. The feature is particularly evident in the R-band. The γ- ray light curve has the first sub-peak (in Fig. 1 corresponding to the P1 vertical dashed line) located immediately before Phase A at MJD 57038.5 (16 January 2015). That indicates the optical/γ- ray emission as possible precursors of the VHE activity, whose peak starts to rise after MJD 57040 (18 January 2015) as indi- cated by the Gaussian fit of the VHE light curve in the top panel of Fig. 1. The Phase B flare is very different, being clearly visible in the VHE and X-ray band only. All the other bands are in a quies- cent level, perhaps reproducing the conditions of the X-ray flare in December 2009 in Rani et al. (2013), where VHE data were not available. In that case, the X-ray emission was described by both synchrotron and inverse Compton mechanisms in a single- zone, one-population leptonic model. The fractional variability Fvar has been calculated using equation 10 in Vaughan et al. (2003): Fvar = √ S 2 − σ2 err x2 , (1) which represents the normalized excess variance. S stands for the standard deviation and σ2 err the mean square error of the flux measurements, while x̄ indicates the average flux. The uncer- tainty of Fvar is given by Eq.(2) in Aleksić et al. (2015), after Poutanen et al. (2008). Fvar was calculated for all the light curves shown in Fig. 1 and the results are plotted in Fig. 2 for both Phase A (full black dots) and Phase B (red open circles). To make a direct comparison of the variability determined for the various energy bands, we computed Fvar using only the multi-instrument observations strictly simultaneous to those per- formed by MAGIC. The overall behaviour of the fractional variability shows a rising tendency with increasing energy, at least up to the X- ray frequency. Since Fvar is highly sensitive to the sampling of 7 M. L. Ahnen et al.: MWL characterization of S5 0716+714 during an unprecedented outburst 56960 56980 57000 57020 57040 57060 57080 57100 57120 57140 ] -1 s -2 cm -6 [1 0 -r ay ) γ F ( 0 0.5 1 6−10× P1 P2 R4 P3 R5 Fermi-LAT 0.1-100 GeV Fermi-LAT 3FGL 0.1-100 GeV 56960 56980 57000 57020 57040 57060 57080 57100 57120 57140 [m Jy ] (R -b an d) ν F 20 40 60 80 Tuorla AZT-8+ST7 LX-200 Kanata Tuorla AZT-8+ST7 LX-200 Kanata Tuorla AZT-8+ST7 LX-200 Kanata 56960 56980 57000 57020 57040 57060 57080 57100 57120 57140 [m Jy ] (U -V -B -b an ds ) ν F 0 20 40 60 UVOT B UVOT U UVOT W1 UVOT W2 UVOT M2 56960 56980 57000 57020 57040 57060 57080 57100 57120 57140 [J y](M ill im et re ) ν F 2 4 6 SMA 230GHz IRAM 230 GHz CARMA 90 GHz IRAM 86GHz hovi 37 GHzaMets 56960569805700057020570405706057080571005712057140 [J y](R ad io ) ν F 1.5 2 2.5 3 3.5 OVRO 15 GHz Effelsberg 15GHz Effelsberg 10GHz Effelsberg 4.8GHz MJD 56960 56980 57000 57020 57040 57060 57080 57100 57120 57140 31-10-2014 31-12-2014 02-03-2015 02-05-2015 A B Fig. 3: HE γ-ray, optical and radio activity of S5 0716+714 during the period from MJD 56950 to MJD 57150 (20 October 2014 to 08 May 2015). The three shadowed areas indicate (from the left side) Phase A (from 18 to 27 January 2015 – MJD 57040 to MJD 57050), an intermediate phase (from 28 January to 11 February 2015 – MJD 57050 to MJD 57064), and Phase B (from 12 to 17 February 2015 – MJD 57065 to MJD 57070). Vertical dashed lines mark important dates, as shown in Table 1. the observed data, in the HE and VHE bands where the sam- pling is poorer, the results are affected by very large error bars, making impossible to confirm a general trend of Fvar for the whole energy spectrum. The highest fractional variabilities in Phase A occur in the VHE γ-ray band and in the X-ray band, with MAGIC (Fvar=0.37 ± 0.08), and Swift-XRT (Fvar=0.28 ± 0.010). In Phase B the highest fractional variability is again in the VHE range with MAGIC (Fvar=0.23 ± 0.13), not far from the one obtained by UVOT telescope (Fvar=0.18 ± 0.8×10−3). 3.1. Analysis of radio activity in a larger time-frame When MAGIC detected S5 0716+714 for the first time in 2008 (Anderhub et al., 2009), the radio band was in a quiescent level. Here we see an increased activity in the low radio frequencies, 8 M. L. Ahnen et al.: MWL characterization of S5 0716+714 during an unprecedented outburst especially in the intermediate period between Phase A and B. This activity could be just an effect of a previous smaller flare in high energy delayed by months as it seems typical for this source when considering longer periods of observation (Rani et al., 2013, 2014). The present work includes one month of data, from the be- ginning of Phase A to the end of Phase B (MJD 57040 to MJD 57070). To have a better understanding of the radio be- havior of the source we gathered γ-ray, optical and radio data in Fig. 3, for a longer time-period of 8 months centered on 28 January 2015 (MJD 57050). VHE emission was not detected outside of Phase A and Phase B, but the data from the avail- able instruments in the 8-months time window make it possi- ble to better investigate the radio response and to compare the present dataset with the scenarios described in other previous multi-wavelength flares. In Fig. 3 Phase A and Phase B are still defined by grey shad- owed areas, while the intermediate phase is filled in light red. From Fig. 3, second panel from the bottom, the radio activity in the intermediate zone between Phase A and Phase B (from MJD 57051 to MJD 57065–29 January to 12 February 2015) could be fit by a Gaussian shape (χ2/n.d. f . = 34.9/14 and stan- dard deviation of σ = 9.7±0.3 days), with a maximum flux den- sity of 4.42±0.07 Jy, corresponding to MJD 57056 (03 February 2015). In the bottom panel of Fig. 3, which shows the radio activity in lower frequencies from the Effelsberg telescope, a peak could be identified by a Gaussian fit of the data (χ2/n.d. f . = 28.7/47 and standard deviation of σ = 14.5 ± 0.5 days), for a maxi- mum flux density at 15 GHz of 3.00 ± 0.02 Jy, corresponding to MJD 57057 (04 February 2015). If we consider the dashed vertical line R4 as the position of the radio peak (using the value of MJD 57056–03 February 2015 retrieved from the Metsähovi data which has a smaller error), we can see a delay with respect to the γ-ray/optical peak P2 of ∼ 9 days, which is considerably smaller than the delays of ∼ 65 days or more found in Rani et al. (2014). Another hint of delayed activity, indicated by the dashed vertical line labeled as R5, can be seen in Metsähovi data, with a maximum flux density of 3.48 ± 0.06 Jy, corre- sponding to MJD 57092 (11 March 2015), fit by a Gaussian with χ2/n.d. f . = 84.3/13 and standard deviation of σ = 17.3 ± 0.9 days. In the latter case, the delay from the P2 flare would be longer, ∼ 45 days. Based on data from April 2007 to January 2011, a consid- erable delay of the radio, ∼ 65 days from the optical/γ flare, was found in Rani et al. (2013). In Rani et al. (2014) using a dataset from August 2008 to September 2013, the highest peak in radio flux occurred ∼ 82 ± 32 days after the γ-ray one. The delay in the radio emission from the optical/γ ones is support- ing a scenario in which the γ-ray emission is produced upstream of the core while the radio one has its origin in a shock in the jet, first appearing and evolving in the innermost, ultracompact VLBI core region and subsequently moving downstream the jet at parsec scales with apparent superluminal speeds. Similar re- sults, on a larger sample of blazars, are presented in Fuhrmann et al. (2014). A longer-term study centered on the flaring activity reported here could be interesting for future investigations but is beyond the scope of the present work. 3.2. Electric Vector Position Angle swing An important feature of Phase A is determined by the very fast change in the electric vector position angle (EVPA) happening over the night MJD 57044 (22 January 2015), 4 days after the first peak P1 in the optical band and 2 days before the MAGIC peak P2. This particular feature can be seen in the bottom panel of Fig. 1, as well as in Fig. 4, where the feature is zoomed in the time range MJD 57043-MJD 57048 (21-26 January 2015). MJD 57044.5 57045 57045.5 57046 57046.5 ]° E V P A [ 100− 0 100 200 300 400 500 Kanata AZT-8+ST7 Steward RINGO3 Fig. 4: Zoom of the Electric Vector Position Angle rotation of S5 0716+714 during the period from MJD 57043 to MJD 57048 (21 to 26 January 2015) in Phase A. The dataset we presented put together EVPA data coming from many different instruments, as shown in the bottom panel of Fig. 1. All the EVPA data collected are in agreement and were treated as in Larionov et al. (2013): in particular, to solve the ±180◦ ambiguity, we have added/subtracted 180◦ each time that the subsequent value of EVPA was > 90◦ less/more than the preceding one. In Marscher et al. (2008), a similar behaviour of the EVPA was reported for BL Lacertae: a radio to γ-ray outburst, accom- panied by a rotation of the EVPA, was observed as the conse- quence of a bright feature moving in the jet. In that case, the EVPA rotation was slower: it rotated steadily by about 240◦ over a five-day interval before settling at a value of 195◦. In 2008, when the source was observed in the VHE range for the first time (Anderhub et al., 2009), the simultaneous op- tical outburst was accompanied by a ∼ 360◦ PA rotation of the electric vector, as reported in Larionov et al. (2008a). The rota- tion happened with an approximate rate of 60◦ per day, so slower than in the present case. That rotation was interpreted as a con- sequence of the propagation of polarized emission from a knot spiraling down the jet. Chandra et al. (2015a) investigated the PA rotation swing from Phase A using the data from the Steward Observatory, in the frame of the HMFM (Zhang et al., 2014), suggesting that the fast rotation was due to reconnections in the emission region in the jet. 4. Jets evolution study with very-long-baseline interferometry We analyzed the structure of the jet with very-long-baseline interferometry. Details of the analysis have been presented in Section 2.5. As shown in Fig. 5, the core of the jet is the brightest knot, designated as A0, which is the southern-most feature of the jet, and assumed to be stationary. The position of the core across 9 M. L. Ahnen et al.: MWL characterization of S5 0716+714 during an unprecedented outburst Fig. 5: A sequence of total (contours) and polarized (color scale) intensity images of S5 0716+714 at 43 GHz, convolved with a beam of 0.24×0.15 mas2 at PA=-10◦. The global total intensity peak is 2655 mJy/beam and the global polarized intensity peak is 107 mJy/beam; black line segments within each image show the direction of polarization; the black horizontal line indicates the position of the core, A0. (a) Separation of knots A1 (red), K14a (yellow), K14b (blue), and K15 (light blue) from the core A0 (black dotted line); the yellow and blue line segments on the position of A0 indicate the 1σ uncertainty of the ejection times of K14a and K14b, respec- tively. (b) Light curves of the core A0 (black), stationary feature A1 (red), and moving knots K14a (yellow) and K14b (blue), flux densities of K14b are shifted by 0.1∼Jy for clarity; vertical blue and red lines indicate time of passage of K14b through A0 and A1, respectively. Fig. 6: VLBA-BU-BLAZAR analysis of S5 0716+714. epochs is shown by the black horizontal line. In addition to the core, we have identified a stationary feature, A1, and two moving knots, K14a and K14b. Fig. 6a shows the evolution of the posi- tions of the knots. Note that a stationary feature at a position sim- ilar to that of A1 was reported previously by Rani et al. (2015) and Jorstad et al. (2017) in data obtained from 2007 to 2013. The moving knots have apparent speeds of (10.7 ± 0.8) c and (9.7 ± 1.8) c for K14a and K14b, respectively, although K14a decelerates at ∼0.5 mas from the core. The direction of motion Φ is ΦK14a=(25.4±2.2) deg and ΦK14b=(43.3±5.6) deg respec- tively. Extrapolation of the motion of K14a and K14b back to the VLBI core suggests that K14a and K14b passed through the core on MJD 56880 ± 22 and MJD 56971 ± 30, respectively. Knot K14b is of special interest with respect to the high en- ergy event in January 2015 (MJD 57040-57050): according to its proper motion of (0.51 ± 0.09) mas/yr, K14b passed through A1 on MJD 57050 ± 30. This coincides with the high γ-ray state and TeV detection of S5 0716+714. Moreover, Fig. 6b shows that A1 had an elevated flux density at the epoch closest to MJD 57050 when K14b should be passing through A1 according to its proper motion. In addition a change of the position angle of A1 from ∼80◦ to ∼50◦ is detected around the time of the expected pas- sage (see Table A.1). The latter angle is close to the PA of K14b, Θ=(47.5 ± 4.6) deg. This implies an interaction between super- luminal knot K14b and stationary feature A1, and supports the hypothesis that A1 is a recollimation shock similar to that ob- served, e.g., in BL Lacertae (Marscher et al., 2008; Cohen et al. , 2014). In addition, the average size of A1 is (0.049 ± 0.020) mas, which implies that K14b needs (35 ± 13) days to pass through the stationary feature. The latter agrees very well with the dura- tion of 34 days of the elevated γ-ray flux in the Fermi light curve of S5 0716+714, from MJD 57032 to 57066 (10 January 2015 to 13 February 2015). In this scenario the TeV detections can be associated with the entrance and exit of the superluminal knot in and out of the recollimation shock. 10 M. L. Ahnen et al.: MWL characterization of S5 0716+714 during an unprecedented outburst 5. Spectral fitting of NuSTAR and Swift-XRT data For NuSTAR data, we performed the spectral fitting with XSPEC v12.8.2, with the standard instrumental response matrices and effective area files derived using the ftool nuproducts. We fit the data for both NuSTAR detectors simultaneously, allowing an offset of the normalization factor for module FPMB with respect to module FPMA. Regardless of the adopted models, the nor- malization offset was less than 5%. First, we adopted a simple power-law model modified by the effects of the Galactic absorp- tion, corresponding to a column of 3.11 × 1020 cm−2 (Kalberla et al., 2005). The fit returns the power-law index of 1.93 ± 0.04, but the residuals show that the NuSTAR spectrum is more con- cave (i.e., the spectrum gets flatter towards higher energies) than a simple power-law model would imply. In addition, this index is significantly harder than that inferred from the Swift-XRT data alone (which shows the index of ∼ 2.75 ), which also suggests a more complex spectral model. Since the NuSTAR and Swift data were nearly strictly con- temporaneous (with a significant overlap) we fit the NuSTAR and Swift-XRT data simultaneously, but allowing for the nor- malizations of Swift and NuSTAR to fit independently. We at- tempted two more complex models (both with absorption fixed at the Galactic value as above). First, we considered a broken power law, with steeper low-energy and harder high-energy in- dices. This is similar to the model considered by Wierzcholska & Siejkowski (2016). The low- and high-energy indices are, re- spectively, 2.52 ± 0.07 and 1.81 ± 0.08, the break energy is at 5.2+0.7 −0.5 keV, and χ2 is 351 for 328 PHA channels. We note that the break energy in our fit is somewhat lower than that deter- mined by Wierzcholska & Siejkowski (2016), but this is likely due to a different choice of bandpass, size of the source extrac- tion region, and precise location of the region of the detector from which the background was subtracted. We also attempted a double power-law representation of the data, also modified by Galactic absorption as above: here, the resulting spectrum is a sum of two power-law models, and is probably more physically motivated than a broken power law. The fit returns χ2 = 352 for 328 PHA bins with a low-energy index of 2.62 ± 0.16 and a high-energy index of 1.41 ± 0.22. Since this model can represent a physically sensible superposi- tion of two separate components, we express a preference for the two-power-law spectral form. With this model, the 2 − 10 keV flux is (9.7 ± 0.7) × 10−12 erg cm−2 s−1. We note here that the most reasonable interpretation of such a 2-component spectral shape is that we witness a contribution of two separate compo- nents, namely the “tail” of the low-energy component (presum- ably produced by the synchrotron process) and the onset of the high energy component (presumably due to the IC process). We plot the unfolded spectrum of the Swift XRT and NuSTAR data observed on 24 January 2015 (MJD 57046) and fit to the two- power-law model in Fig. 7. 6. VHE Differential Energy Spectrum and EBL deabsorbtion The VHE γ-rays from distant blazars can interact with the optical-UV photons from the extragalactic background light (EBL)(Gould & Schreder, 1967; Stecker et al., 1992) via pair production, resulting in an attenuation of the intrinsic VHE spec- trum. Finite resolution of the instrument will also modify the in- trinsic spectrum. Unfolding techniques are adopted in the MARS code to unfold the observed spectrum from the instrument re- sponse. The differential spectrum of S5 0716+714 is shown in Fig. 7: Unfolded X-ray spectrum of S50716+714 derived from simultaneous fitting of the contemporaneous Swift-XRT and NuSTAR data obtained on 24 January 2015 (MJD 57046). The adopted model is a sum of two power laws. Swift-XRT data are plotted in green, while the two NuSTAR modules are plotted in red and black, respectively. The “valley” between the two main broad-band spectral peaks is in the X-ray band. Fig. 8 for a simple unfolding considering instrumental response only (hereafter observed spectrum). An unfolding including also de-absorption from EBL with the Domı́nguez et al. (2011) model (hereafter the intrinsic spectrum) was also performed, and pa- rameters of the observed and intrinsic spectra are reported in Table 2. The EBL imprint on the γ-ray spectra from distant blazars could be used to constrain the EBL density, under some assumptions on the intrinsic spectrum of the source (see e.g. Ackermann et al., 2012; Abramowski et al., 2013). The differ- Energy [GeV] 200 300 400 500 600 700 ] -1 s -2 c m -1 /d E [T eV φd 13−10 12−10 11−10 10−10 9−10 8−10 7−10 MAGIC Phase A, with 11.95 hours total MAGIC Phase B, with 5.49 hours total Fig. 8: Unfolded observed differential energetic spectra by MAGIC for Phase A (black full dots) and Phase B (red full squares). Parameters for the spectra (including the ones regard- ing the intrinsic EBL deabsorbed spectra with Domı́nguez et al. (2011) model) are reported in Table 2. ential VHE spectra, observed as well as EBL-corrected using the model of Domı́nguez et al. (2011), can be described by a power- law : 11 M. L. Ahnen et al.: MWL characterization of S5 0716+714 during an unprecedented outburst Table 2: VHE spectrum parameters for a PWL fit f0(cm−2 s−1 TeV−1) Γ χ2/n.d.f. P E (GeV) -Phase A- observed (2.75 ± 0.25stat ± 0.30sys) × 10−10 4.08 ± 0.22stat ± 0.15sys 0.11/3 0.99 127.2-659.1 intrinsic (4.85 ± 0.45stat ± 0.53sys) × 10−10 2.73 ± 0.28stat ± 0.15sys 1.68/3 0.64 127.2-659.1 - Phase B - observed (7.94 ± 1.27stat ± 0.87sys) × 10−10 4.64 ± 0.49stat ± 0.15sys 0.17/2 0.92 127.2-474.3 intrinsic (1.06 ± 0.17stat ± 0.11sys) × 10−10 3.65 ± 0.53stat ± 0.15sys 0.23/2 0.89 127.2-474.3 dF dE = f0 ( E 150 GeV )−Γ , (2) where the normalization constant f0, the spectral index Γ , the goodness of the fit (χ2/n.d.f. and probability P), the energy range of the fit E are indicated in Table 2 for Phase A and Phase B data respectively. 6.1. Redshift estimation The simultaneous spectra from MAGIC and Fermi-LAT were used to estimate the redshift of the source. We apply the method presented in Prandini et al. (2010, 2011) based on the assumption that the slope of the VHE spectrum corrected for EBL absorption should not be harder than the one measured by Fermi-LAT at lower energies. The redshift at which the two slopes match, z*, is then used as upper limit estimate of the source distance if there is no spectroscopic redshift available. If we apply the method to the data presented here assuming the Franceschini et al. (2008) EBL model, a 2 σ upper limit on the redshift of 0.598 is found. The empirical formula proposed in Prandini et al. (2011) ap- plied to this data gives as most probable value for the redshift zrec = 0.31 ± 0.02stat ± 0.05sys, where the first error is related to the statistical errors of Fermi-LAT and MAGIC slopes, while the second error is the error of the method itself, as estimated in Prandini et al. (2011). This value is in agreement with the ones given in literature by Nilsson et al. (2008); Danforth et al. (2013). The value of z = 0.31 ± 0.08 found in Nilsson et al. (2008) was based on the photometric detection of the host galaxy, while the z < 0.322 (95% confidence) result reported by Danforth et al. (2013) was obtained by detection of Ly-α systems in the ultra-violet spectrum of the source. For the SED modelling (see next section), we used the redshift value of 0.26 as in Anderhub et al. (2009). This value is within the errorbars of the redshift determined here as well as within other observations (see the Introduction). 7. Broadband Spectral Energy Distribution The multi-wavelength SEDs for Phase A and B respectively are presented in Fig. 9. Archival data from ASDC (ASI Science Data Center)13 are shown in grey. When the source was detected for the first time in the VHE range (Anderhub et al., 2009) by MAGIC, the only available simultaneous multi-wavelength data were coming from KVA (optical) and Swift (X-rays) and there were no constraints on the second bump beyond the MAGIC data. Nevertheless the very soft X-ray spectra belonging to the synchrotron component, combined with the high VHE γ-ray flux 13 http://www.asdc.asi.it/ challenged the simple one-zone SSC model as it would require a very high flux of γ-rays around 10 GeV, higher than has been ob- served from the source with Fermi-LAT or its precursor EGRET. This condition for one-zone persists also for the new data, but now we actually have simultaneous data in this energy range from Fermi-LAT and we find that we cannot describe the ob- served broadband SED with the one-zone SSC model during this flaring period. While a one-zone SSC model can match the observed Swift- XRT+NuSTAR spectrum as a transition between synchrotron and IC components, and simultaneously the γ-ray data from Fermi-LAT and MAGIC, it tends to under-reproduce the ob- served optical flux. This has been independently verified using the BLAZAR code (Moderski et al., 2003). Based on the multiwavelength data in Section 3 and VLBA data in Section 4, we use two-zone model to describe the SED in Phase A and Phase B. We use two blobs close to each other to represent a situation where a superluminal knot (blob 1) is inter- acting a recollimation shock region (blob 2). The Phase A SED represents a snapshot of a time when the knot enters the rec- ollimation shock region and the Phase B SED a time when the knot has exited the recollimation region. We model the two blobs with the framework similar to one presented for flat-spectrum radio quasar PKS 1222+216 in Tavecchio et al. (2011), modi- fied for the case of no external seed photons as in Aleksić et al. (2014) for PKS 1424+240. Unlike the case of PKS 1424+240, in our case the two emission regions are very close to each other and they provide seed photons for inverse Compton scattering for each other. The model blobs have broken power-law electron spectra (γmin, γb and γmax are the minimum, break and maximum Lorentz factors respectively; n1 and n2 are the low and high en- ergy slope of the smoothed power law electron energy distribu- tion), magnetic field B, normalization of the electron distribution K, radius of the emission region R and Doppler factor δ. We also use the observed variability behaviour as a guide on how the parameters change between Phase A and Phase B. As discussed in Section 3, Phase A and Phase B have different variability behaviours; while Phase A consists of a flare in all bands, in Phase B the activity is constrained to the X-ray and VHE γ-ray bands. To limit the number of free parameters, we fix the larger component (which is representing the recollima- tion shock) to have most of the parameters the same in Phases A and B. We only change K to a lower value to represent the gen- eral lower state in optical and GeV γ-rays of Phase B. We then find parameters for the smaller emission region to describe the observed SEDs in Phase A and B separately. In both panels of Fig. 9 the ”blob1” component is represented by the red dashed line. The ”blob2” component of the model is represented by the blue dash-dotted line in both panels, and the emission resulting from the interaction of two components is reported with a green line. 12 M. L. Ahnen et al.: MWL characterization of S5 0716+714 during an unprecedented outburst Energy [GeV] 12−10 8−10 4−10 1 410 ] -1 s -2 d N /d E [T eV c m 2 E 14−10 13−10 12−10 11−10 10−10 9−10 MWL SED Phase A Energy [GeV] 12−10 8−10 4−10 1 410 ] -1 s -2 d N /d E [T eV c m 2 E 14−10 13−10 12−10 11−10 10−10 9−10 MWL SED Phase B Fig. 9: MWL Spectral Energy distributions for Phase A and Phase B. Archival data form ASDC are shown in grey. The two compo- nents (blobs representing a moving emission feature and a recollimation shock, see text) are shown with blue and red dashed lines. The green line is the emission that is a result of interaction between these two blobs and the black solid line the sum of these three components. The red full circles represent the intrinsic (EBL deabsorbed according to (Domı́nguez et al., 2011)) MAGIC SED used in the model. For data taken in the radio and optical band the error bars are smaller than the size of the marker. Table 3: Input parameters for the emission models of S5 0716+714 γmin γb γmax n1 n2 B (G) K R (cm) δ z -Phase A- ”blob1” 100 1.3 × 104 2 × 105 1.95 3.5 0.1 3.5 × 103 2.6 × 1016 25 0.26 ”blob2” 300 1 × 104 1.5 × 105 2.32 4.6 0.12 1.6 × 104 6.5 × 1016 25 0.26 - Phase B - ”blob1” 4 × 103 8 × 103 2 × 105 1.9 3.2 0.09 9 × 103 1.35 × 1016 25 0.26 ”blob2” 300 1 × 104 1.5 × 105 2.32 4.6 0.12 1 × 104 6.5 × 1016 25 0.26 We report a set of parameters we found to give a reasonable description (but see below) of the observed SED in Table 3 for the both phases. The set of parameters we present is not unique, but the parameters used are within the range typically found for TeV blazars and also for the ones found for PKS 1424+240 using the similar modelling setup. Even if this simple two-zone model provides a better rep- resentation of the observed data from radio to VHE γ-rays re- spect to a one-zone model, the model line is slightly lower than the γ-ray fluxes in the range of 10-100 GeV even though within the systematic uncertainties of the data. The data in this energy range suggests a rather sharp feature, which is impossible to re- produce with this simplistic model we used. In general sharp features require presence of external seed photons such as used e.g. by Böttcher et al. (2013) to model the SED of the source in a lower state. However, there is no observational evidence for such an external seed photon field from optical spectroscopy nor from the scenario we presented for the flaring behaviour within this epoch therefore no such a component was added to the mod- elling. There are also other possible two-zone model setups, such as a spine-sheath model (Ghisellini et al., 2005), where a slower sheath of the jet surrounds a faster spine. For the previous VHE γ-ray flaring epoch a spine-sheath mode (Ghisellini et al., 2005), was shown to provide a reasonable fit to the SED data (Anderhub et al., 2009; Tavecchio & Ghisellini, 2009). We tested the spine- sheath model for the spectral energy distributions shown here and found acceptable agreement with the SED data. This em- phasizes that SED data alone are not enough to separate different two-zone models, but must be combined with constraints from VLBA and light curve variability. According to the scenario described above, which is sup- ported by the dedicated VLBI study we performed in Sec. 4, we explain the extremely fast rotation of ∼ 360◦ as produced by turbulence in the interaction between a superluminal knot and a stationary feature near the core. Being dependent on the orienta- tion of the shock and the magnetic field threading it, EVPA pro- vides a unique tool to understand the acceleration mechanisms and behavior of the shocked plasma. Recent studies on EVPA swings larger than 180◦ simultaneous with HE γ-ray emission (Marscher et al., 2008, 2010; Abdo et al., 2010) have been inter- preted as additional evidence for a helical magnetic field struc- ture. The existing models focusing on the description of the syn- chrotron polarization features (e.g. Lyutikov et al., 2005) ap- ply a simple and time-independent power-law electron spectrum not taking into account possible predictions for the resulting HE emission. Our model, on the other hand, does not include a de- tailed geometry of the magnetic field and the angle-dependent synchrotron emissivity and polarization. At the moment only two models may represent the SED of blazars together with their synchrotron polarization features, including rotations of the EVPA: the HMFM (Helical Magnetic Field Model) (Zhang et al., 2014) and the TEMZ (Turbulent, Extreme Multi-Zone) (Marscher, 2014). In the HMFM large polarization angle rota- 13 M. L. Ahnen et al.: MWL characterization of S5 0716+714 during an unprecedented outburst tions by & 180◦ are explained with the passage of a moving shock through a region with a highly disordered field: the com- pression of the shock orders the field partially, but this order- ing is seen at different depths as time advances owing to light- travel delays. This leads to an apparent rotation of the polariza- tion of 180◦ per shock. In the TEMZ model randomness in the magnetic field direction in different turbulent cells can cause ob- served rotations in the linear polarization vector, even fast as the one observed in our case. Turbulence in general gives at different times ”clusters” with small EVPA variation, relatively smooth EVPA rotations, step-wise EVPA changes, and random fluctu- ations. The behavior of S5 0716+714 in Phase A is consistent with this scenario. 8. Summary and Conclusions The BL Lac object S5 0716+714 has been studied in a multi- wavelength frame from radio to the VHE γ-ray band. In January 2015 an unprecedented outburst of S5 0716+714 was registered in all energy bands, from low frequency radio to VHE, and after almost a month another high state was detected by the MAGIC and Swift-XRT instruments only. We divide the data into two phases (Phase A from 18 to 27 January 2015 – MJD 57040 to MJD 57050, and Phase B from 12 to 17 February 2015 – MJD 57065 to MJD 57070) that represent very different charac- teristics, allowing a deep study of the broadband SEDs. The broadband flaring activity period of Phase A coincides with the passage of a moving feature through a stationary feature (at ∼0.1 mas). We have found a very fast change in the electric vector position angle during Phase A. The > 400 degree swing in the optical EVPA is explained here as the passage of a su- perluminal knot through a stationary feature near the radio core. The VHE emission is then found originating in the entrance and exit of a superluminal knot in and out a recollimation shock in the inner jet. This suggests that shock-shock interaction in the jet seems to be responsible for the observed flares and EVPA swing. The jet behaviour, studied with VLBA-BU-BLAZAR data, is in agreement with the scenario described in Rani et al. (2015), suggesting a connection between jet kinematics and the observed broadband flaring activity. More precisely, the γ-ray emission in the HE and VHE bands is attributed to a shock in the helical jet downstream of the core, closely followed by an optical and X-ray outburst in the core. The presence of low radio activity, observed during Phase A, was not reported in April 2008, when MAGIC observed the source in the VHE range for the first time (Anderhub et al., 2009), but it could be a delayed response of a previous less intense flare, as observed in the past in the same source, when between optical/γ flares lagged the radio counter- parts almost two months (Rani et al., 2013, 2014). The first peak in the VHE γ-ray emission takes place ∼ 2 days after the very fast EVPA rotation and the second ∼ 18 days after the new knot has been emerged from the VLBA core. This is a strong indication that the VHE γ-ray emission is associated to a component entering and exiting the core region. The broadband SEDs, for the first time including MAGIC and Fermi-LAT simultaneous data and the quasi-simultaneous NuStar data, could not be described by a simple one-zone model. Instead we used a two-zone model, where two spherical blobs are co-spatial and provide seed photons to each other. This mod- elling setup provides an acceptable description of the spectral energy distributions in Phase A and B, even if it is certainly an over-simplified presentation of the true physical processes taking in place when superluminal knots enter and exit the recollima- tion shock region. Finally we also investigated the redshift of S5 0716+714. Using simultaneous data from MAGIC and Fermi-LAT the red- shift was calculated to be z = 0.31±0.02stat±0.05sys, confirming the value present in literature based on the photometric detection of the host galaxy (Nilsson et al., 2008) and the more recent up- per limit from a direct detection (Danforth et al., 2013). S50716+714 is an intermediate BL Lac object, and only a handful of these sources have been detected in VHE γ-rays. In almost all detections of VHE γ-rays, activity in other bands (op- tical and/or HE γ-rays) has been seen, but our very comprehen- sive dataset provided a unique insight on how these VHE γ-ray flares are connected to the activity in the jet. In the end of December 2017 S5 0716+714 was flaring again in VHE γ-rays (Mirzoyan et al., 2017). It will be interesting to see if the recognized patterns will repeat also during this ongoing flaring period. This will be studied in a future paper. Acknowledgements. 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 ac- knowledged. 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 num- ber 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 sci- ence 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 research was supported by an appointment to the NASA Postdoctoral Program at the Goddard Space Flight Center, administered by Universities Space Research Association through a contract with NASA.We thank the Swift team duty scientists and science planners. The Metsähovi team acknowledges the support from the Academy of Finland to our observing projects (numbers 212656, 210338, 121148, and others). The VLBA is an in- strument of the Long Baseline Observatory. The Long Baseline Observatory is a facility of the National Science Foundation operated by Associated Universities, Inc. The Submillimeter Array is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics and is funded by the Smithsonian Institution and the Academia Sinica. 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. The St. Petersburg University team acknowl- edges support from Russian Science Foundation grant 17-12-01029. The BU group acknowledges support by NASA under Fermi Guest Investigator grant NNX14AQ58G and by NSF under grant AST-1615796. Part of this work was done with funding by the UK Space Agency. The VLBA is an instrument of the National Radio Astronomy Observatory. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under co- operative agreement by Associated Universities, Inc. The PRISM (Perkins Re- Imaging SysteM) camera at Lowell Observatory was developed by K. Janes et al. at BU and Lowell Observatory, with funding from the NSF, BU, and Lowell Observatory. This paper makes use of data obtained with the 100 m Effelsberg radio-telescope, which is operated by the Max-Planck-Institut für Radioastronomy (MPIfR) in Bonn (Germany). Part of this work is based on archival data, software or online services provided by the ASI Science Data Center (ASDC). PYRAF is a product of the Space Telescope Science Institute, 14 M. L. Ahnen et al.: MWL characterization of S5 0716+714 during an unprecedented outburst which is operated by AURA for NASA. This paper is partly based on obser- vations carried out with the IRAM 30m. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). IA acknowledges support by a Ramon y Cajal grant of the Ministerio de Economia, Industria y Competitividad (MINECO) of Spain. The research at the IAA-CSIC was supported in part by the MINECO through grants AYA2016-80889-P, AYA2013-40825-P, and AYA2010- 14844, and by the regional government of Andalucia through grant P09-FQM- 4784. The Liverpool Telescope is operated by JMU with financial support from the UK-STFC. References Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJ, 710, 1271 Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, Nature, 463, 919 Abramowski, A., Acero, F., Aharonian, F., et al 2013, A&A, 550, A4 Acciari, V. A., Aliu, E., Aune, T., et al 2009, ApJ, 707, (1), 612-620 Ackermann, M., Ajello, M., Allafort, A., et al. 2011, ApJ, 743, 171 Ackermann, M., Ajello, M., Albert, A., et al. 2012, Science, 338, 1190 Acero, F., Ackermann, M., Ajello, M. et al. 2015, ApJS, 218, 23 Agudo, I., Thum, C., Molina, S. N., et al. 2018a, accepted in MNRAS Agudo, I., Thum, C., Ramakrishnan, V., et al. 2018b, MNRAS, 473, 1850 Akitaya, H., Moritani, Y., Ui, T. et al. 2014, Proc of the SPIE, 9147, 91474O Aleksić, J., Ansoldi, S., Antonelli, L. A., et al. 2014, A&A, 567, A135 Aleksić, J., Ansoldi, S., Antonelli, L. A., et al. 2015, A&A, 573, A50 Aleksić, J., Ansoldi, S., Antonelli, L. A., et al. 2016, Astroparticle Phys. 72,76 Anderhub, H., Antonelli, L. A., Antoranz, P., et al. 2009, ApJ, 704, L129 Angelakis, E., Fuhrmann, L., Marchili, N., et al. 2015, A&A, 575, A55 Antonucci, R. R. J., Hickson, P., Olszewski, E. W., & Miller, J. S. 1986, AJ, 92,1 Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009 ApJ, 697 1071 Bach, U., Krichbaum, T. P., Ros, E., et al. 2005, A&A, 433, 815 Bachev, R., et al. 2015, ATel 6957 Baars, J. W. M., Genzel, R., Paulini-Toth, I. I. K., & Witzel, A., 1977, A&A, 61, 99 Böttcher, M., Reimer, A., Sweeney, K., & Prakash, A., 2013, ApJ, 768, 54 Britzen, S., Kam, V. A., Witzel, A., et al. 2009, A&A, 508, 1205 Breeveld, A. A., Landsman, W., Holland S. T., et al. 2011, AIP Conf. Proc., 1358, 373 Burrows, D. N., Hill, J. E., Nousek, J. A.,et al. 2004, Proc. SPIE, 5165, 201 Carrasco, L., et al. 2005, ATel 6902 Chandra, S., Zhang, H., Kushwaha, P., et al. 2015a, ApJ, 809, 130 Chandra, S., et al. 2015b, ATel 6962 Chen, A. W., D’Ammando, F., Villata, M., et al. 2008, A&A, 489, L37 Chatterjee, R., Jorstad, S. G., Marscher, A. P., 2008, ApJ, 689, 79-94 Chatterjee, R., Bailyn, C. D., Bonning, E. W., et al. 2012, ApJ, 749, 191 Cohen, M. H., Meier, D. L., Arshakian, T. G., et al. 2014, ApJ, 787, 151 Danforth, C. W., Nalewajko, K., France, K., Keeney, B. A., 2013, ApJ, 764, 57 Domı́nguez, A., Primack, J. R., Rosario], D. J., et al. 2011, MNRAS, 410, 2556- 2578. Franceschini, A., Rodighiero, G., Vaccari, M. et al. 2008, A&A, 487, 837-852 Foschini, L., Tagliaferri, G., Pian, E., et al. 2006, A&A, 455, 871 Fuhrmann, L., Krichbaum, T. P., Witzel, A., et al. 2008, A&A, 490, 1019 Fuhrmann, L., Larsson, S., Chiang, J., et al. 2014, MNRAS, 441, 1899 Ghisellini, G., Tavecchio F., & Chiaberge, M., et al. 2005, A&A, 432, 401 Giommi, P., Massaro E., Chiappetti, L., et al. 1999, A&A, 351, 59 Gould R. J. & Schreder G. P., 1967, Phys. Rev., 155, 1404 Gurwell, M. A., Peck, A. B., Hostler S. R., et al. 2007, ASP Conference Series 375, 234 Hartman, R. C., Bertsch, D. L., Bloom, S. D., et al. 1999, ApJS, 123, 79 Harrison, F. A., Craig W. W., Christensen F. E., et al. 2013, ApJ, 770, 103 Jermak, H., Steele, I. A., Lindfors, E., et al. 2016, MNRAS, 462, 4267 Jorstad, S. G., Marscher, A. P., Mattox, J. R., et al. 2001, ApJ, 556, 738 Jorstad, S. G., Marscher, A. P., Lister, M. L., et al. 2005, AJ, 130, 1418 Jorstad, S. G., Marscher, A. P., Larionov, V. M., et al. 2010, ApJ, 715, 362 Jorstad, S. G., Marscher, A. P., Morozova, D. A, et al. 2017, ApJ, 846, 98 Kalberla, P. M. W., Burton, W. B., Hartmann, D., et al. 2005, A&A, 440, 775 Larionov, V. M., et al. 2008a, ATel 1502 Larionov, V. M., Jorstad, S. G., Marscher, A. P., et al. 2008b, A&A, 492, 389 Larionov, V. M., Jorstad, S. G., Marscher, A. P., et al. 2013, ApJ, 768, 40 Li, T. -P. & Ma,Y. -Q., 1983 ApJ, 272, 317L Lin, Y. C., Bertsch, D. L., Dingus, B. L., et al. 1995, ApJ, 442, 96 Lister, M. L., Aller, M. F., Aller, H. D., et al. 2013, ArXiv e-prints Liu, X., Song, H.-G., Marchili, N., et al. 2012, A&A, 543, A78 Lott, B., Escande, L., Larsson, S., and Ballet, J., 2012, A&A, 544, A6 Lyutikov, M., Pariev, V. I., & Gabuzda, D. C., 2005, MNRAS, 360, 869 Marscher, A. P., Jorstad S. G., D’Arcangelo F. D., et al. 2008, Nature, 452, 966 Marscher, A. P., Jorstad S. G., Larionov, V. M., et al. 2010, ApJ, 710, L126 Marscher, A. P. 2014, ApJ, 780,87 Mattox, J. R., Bertsch, D. L., Chiang, J., et al 1996 ApJ, 461, 396 Mirzoyan, R., et al. 2015, ATel 6999 Mirzoyan, R., et al. 2017, ATel 11100 Moderski, R., Sikora, M., and Błażejowski, M. 2003, A&A, 451, 435 Montagni, F., Maselli, A., Massaro, E., et al. 2006, A&A, 451, 435 Nilsson, K., Pursimo, T., Sillanpää, A., Takalo, L. O., & Lindfors, E. 2008, A&A, 487, L29 Nilsson, K., et al. 2017, submitted to A&A Paiano, S., Landoni, M., Falomo, R., et al. 2017, ApJ, 837, (2), 144 Poole, T. S., Breeveld, A. A., Page, M. J., et al. 2008, MNRAS383, 627 Poutanen, J., Zdziarski, A. A. and Ibragimov, A., 2008, MNRAS389, 1427 Prandini, E., Bonnoli G., Maraschi L., Mariotti M. and Tavecchio F., 2010, MNRAS405, L76 Prandini, E., Bonnoli G., Maraschi L., Mariotti M. and Tavecchio F., 2011, PoS(CRF 2010)012 Quirrenbach, A., Witzel, A., Wagner, S., et al. 1991, ApJ, 372, L71 Raiteri, C. M., Villata, M., Tosti, G., et al. 2003, A&A, 402, 151 Fallah Ramazani, V, Lindfors, E., & Nilsson, K., 2017, A&A, 608, A68 Rani, B., Gupta, A. C., Joshi, U. C., Ganesh, S., & Wiita, P. J. 2010, ApJ, 719, L153 Rani, B., Gupta, A. C., Joshi, U. C., Ganesh, S., & Wiita, P. J. 2011, MNRAS, 413, 2157 Rani, B., Krichbaum, T. P., Fuhrmann, L., et al. 2013, A&A, 552, A11 Rani, B., Krichbaum, T. P., Marscher, A. P., et al. 2014, A&A, 571, L2 Rani, B., Krichbaum, T. P., Marscher, A. P., et al. 2015, A&A, 578, A123 Rastorgueva, E. A., Wiik, K. J., Bajkova, A. T., et al. 2011, A&A, 529, A2+ Richards, J. L., Max-Moerbeck, W., Pavlidou, V., et al. 2011, ApJS, 194, 29 Sault R. J., Teuben P. J., & Wright M. C. H., 1995, ASP Conference Series, 77, 433-436. Schlafly, Edward F. & Finkbeiner, Douglas P. (2011), ApJ, 737, 103 Schulz, A. & Lenzen, R. (1983), A&A, 121, 158 Senturk, G. D., Errando, M., Boettcher, M., et al. 2013, ApJ, 764, 119 Spiridonova, O. I., et al. 2015, ATel 7004 Stecker F. W., DeJager O. C., Salamon M. H., 1992, ApJ, 390, L49 Steele, I. A., Kopač, D., Arnold, D. M., et al. 2017, ApJ, 843, 143 Tagliaferri, G., Ravasio, M., Ghisellini, G., et al. 2003, A&A, 400, 477 Tavecchio, F., & Ghisellini, G., 2009, MNRAS, 394, 13 Tavecchio, F., Becerra-Gonzalez, J., Ghisellini, G., et al. 2011, A&A, 534, A86 Teräsranta, H., Tornikoski, M., Mujunen, A. et al. 1998, A&AS132, 305. Thum, C., Agudo, I., Molina, S. N., et al. 2018, MNRAS, 473, 2506 Uttley, P., Edelson, R., McHardy, I. M., et al. (2003) ApJ, 584, L53 Vaughan, S., Edelson, R., Warwick, R. S. and Uttley, P. 2003, MNRAS, 345, 1271 Wagner, S. J., Witzel, A., Heidt, J., et al. 1996, AJ, 111, 2187 Wierzcholska, A. & Siejkowski, H. 2016, MNRAS, 458, 2350 Witzel, A., Schalinski, C. J., Johnston, K. J., et al. 1988, A&A, 206, 245 Zanin, R., Carmona, E., Sitarek, J., et al. 2013, Proc. of the 33th International Cosmic Ray Conference (ICRC) 2-9 July 2013, Rio de Janeiro, Brazil, 773 Zhang, H., Chen, X., and Böttcher, M., 2014, ApJ, 789, 66 15 M. L. Ahnen et al.: MWL characterization of S5 0716+714 during an unprecedented outburst Appendix A: Knot parameters for S5 0716+714 (July 2014-August 2015) In Table A.1 the parameters of components from the Section 4 are listed. They indicate the flux density, S , the relative RA and DEC with respect to the core, the distance from the core r, the position angle with respect to the core, Θ, and the size of the component, a (FWHM of the Gaussian) for every knot. Some knots present in Table A.1 do not have any assigned name be- cause they have been seen at 1-2 epochs only. 1 ETH Zurich, CH-8093 Zurich, Switzerland 2 Università di Udine, and INFN Trieste, I-33100 Udine, Italy 3 National Institute for Astrophysics (INAF), I-00136 Rome, Italy 4 Università di Padova 5 Technische Universität Dortmund, D-44221 Dortmund, Germany 6 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. 7 Saha Institute of Nuclear Physics, HBNI, 1/AF Bidhannagar, Salt Lake, Sector-1, Kolkata 700064, India 8 Max-Planck-Institut für Physik, D-80805 München, Germany 9 now at Centro Brasileiro de Pesquisas Fı́sicas (CBPF), 22290-180 URCA, Rio de Janeiro (RJ), Brasil 10 Unidad de Partı́culas y Cosmologı́a (UPARCOS), Universidad Complutense, E-28040 Madrid, Spain 11 Inst. de Astrofı́sica de Canarias, E-38200 La Laguna, and Universidad de La Laguna, Dpto. Astrofı́sica, E-38206 La Laguna, Tenerife, Spain 12 University of Łódź, Department of Astrophysics, PL-90236 Łódź, Poland 13 Deutsches Elektronen-Synchrotron (DESY), D-15738 Zeuthen, Germany 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ät Würzburg, D-97074 Würzburg, Germany 17 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 18 Departament de Fı́sica, and CERES-IEEC, Universitat Autónoma de Barcelona, E-08193 Bellaterra, Spain 19 Universitat de Barcelona, ICC, IEEC-UB, E-08028 Barcelona, 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; The University of Tokushima, 770-8502 Tokushima, Japan 21 Inst. for Nucl. Research and Nucl. Energy, Bulgarian Academy of Sciences, BG-1784 Sofia, Bulgaria 22 Università di Pisa, and INFN Pisa, I-56126 Pisa, Italy 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 INFN, I-35131 Padova, Italy 28 ASI Science Data Center and INFN, 06123 Perugia, Italy 29 CEN Bordeaux-Gradignan, 33170 Gradignan, France 30 NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA 31 Instituto de Astrofsica de Andaluca (CSIC), Apartado 3004, E- 18080 Granada, Spain 32 Max–Planck–Institut für Radioastronomie, Auf dem Hügel, 69, D– 53121, Bonn, Germany 33 Crimean Astrophysics Observatory, P/O Nauchny, Crimea, 298409, Russia 34 Astron. Inst., St. Petersburg State University, Russia 35 Harvard-Smithsonian Center for Astrophysics, MA 02138 Cambridge, USA 36 Metsähovi Radio Observatory, Aalto University, FI-02540 Kylmälä 37 Department of Electronics and Nanoengineering, Aalto University, FI-00076 Aalto, Finland 38 Tuorla Observatory, University of Turku, Väisäläntie 20, 21500 Piikkiö, Finland 39 Department of Physical Science, Hiroshima University, Higashi- hiroshima, 739-8526, Japan 40 Institute for Astrophysical Research, Boston University, Boston MA 02215 41 Mullard Space Science Lab., UCL, Dorking, RH5 6NT, UK 42 Pulkovo Observatory, St.-Petersburg, Russia 43 Kavli Institute for Particle Astrophysics and Cosmology, Stanford University and SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025 44 Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, Bartycka 18, 00-716 Warsaw, Poland 45 Owens Valley Radio Observatory, California Institute of Technology, Pasadena, CA 91125, USA 46 CePIA, Astronomy Department, Universidad de Concepcion, Casilla 160-C, Concepcion, Chile 47 Astrophysics Research Institute, Liverpool John Moores University, Brownlow Hill, Liverpool, L3 5RF, UK 16 M. L. Ahnen et al.: MWL characterization of S5 0716+714 during an unprecedented outburst Table A.1: Knot parameters for S5 0716+714 (July 2014-August 2015) Epoch MJD S (Jy) x(mas) y(mas) r(mas) Θ(deg) a(mas) Knot - 2014-Jul-28 - 2014.5726 56867 1.619 0.000 0.000 0.000 0.0 0.020 A0 2014.5726 56867 0.548 0.092 0.005 0.093 86.8 0.031 A1 2014.5726 56867 0.017 0.186 0.551 0.582 18.7 0.186 - 2014-Sep-23 - 2014.7288 56924 1.322 0.000 0.000 0.000 0.0 0.020 A0 2014.7288 56924 0.205 0.095 0.021 0.097 77.6 0.031 A1 2014.7288 56924 0.033 0.158 0.457 0.484 19.1 0.189 - 2014-Nov-15 - 2014.8740 56977 2.676 0.000 0.000 0.000 0.0 0.022 A0 2014.8748 56977 0.105 0.095 0.021 0.097 77.6 0.031 A1 2014.8740 56977 0.306 0.178 0.095 0.202 62.0 0.112 K14a - 2914-Dec-5 - 2014.9288 56997 1.717 0.000 0.000 0.000 0.0 0.019 A0 2014.9288 56997 0.571 0.078 0.055 0.095 54.8 0.038 A1 2014.9288 56997 0.183 0.206 0.139 0.248 56.1 0.094 K14a 2014.9288 56997 0.019 0.163 0.440 0.469 20.3 0.274 - 2014-Dec-29 - 2014.9945 57021 1.171 0.000 0.000 0.000 0.0 0.017 A0 2014.9945 57021 0.168 0.093 0.077 0.121 50.5 0.023 A1 2014.9945 57021 0.245 0.250 0.179 0.307 54.4 0.130 K14a 2014.9945 57021 0.021 0.246 0.602 0.651 22.2 0.416 - 2015-Feb14 - 2015.1233 57067 1.901 0.000 0.000 0.000 0.0 0.013 A0 2015.1233 57067 0.314 0.093 0.080 0.123 49.2 0.065 A1 2015.1233 57067 0.122 0.232 0.222 0.321 46.3 0.155 K14a 2015.1233 57067 0.036 0.236 0.527 0.577 24.1 0.412 - 2015-Apr-11 - 2015.2767 57123 2.094 0.000 0.000 0.000 0.0 0.017 A0 2015.2767 57123 0.373 0.069 0.056 0.089 51.3 0.057 A1 2015.2767 57123 0.045 0.163 0.169 0.235 44.1 0.100 K14b 2015.2767 57123 0.089 0.287 0.390 0.484 36.3 0.179 K14a - 2015-May-11 - 2015.3589 57153 1.581 0.000 0.000 0.000 0.0 0.023 A0 2015.3589 57153 0.309 0.074 0.065 0.098 48.7 0.063 A1 2015.3589 57153 0.032 0.204 0.165 0.262 50.9 0.095 K14b 2015.3589 57153 0.069 0.286 0.407 0.498 35.1 0.201 K14a - 2015-Jun-9 - 2015.4385 57182 0.947 0.000 0.000 0.000 0.0 0.018 A0 2015.4385 57182 0.308 0.106 0.098 0.144 47.5 0.067 A1 2015.4385 57182 0.039 0.229 0.191 0.298 50.3 0.069 K14b 2015.4385 57182 0.054 0.320 0.437 0.542 36.2 0.189 K14a - 2015-Jul-2 - 2015.5014 57205 1.612 0.000 0.000 0.000 0.0 0.021 A0 2015.5014 57205 0.124 0.020 0.060 0.063 18.7 0.000 K15 2015.5014 57205 0.144 0.075 0.096 0.122 37.9 0.084 A1 2015.5014 57205 0.041 0.253 0.204 0.325 51.1 0.084 K14b 2015.5014 57205 0.057 0.322 0.384 0.531 40.0 0.141 K14a 2015.5014 57205 0.016 0.376 0.665 0.764 29.5 0.401 - 2015-Aug-1 - 2015.5836 57235 1.865 0.000 0.000 0.000 0.0 0.019 A0 2015.5836 57235 0.238 0.045 0.090 0.101 26.3 0.033 K15 2015.5836 57235 0.245 0.087 0.126 0.153 34.6 0.083 A1 2015.5836 57235 0.032 0.241 0.278 0.368 41.0 0.107 K14b 2015.5836 57235 0.018 0.388 0.455 0.598 40.4 0.119 K14a 2015.5836 57235 0.020 0.315 0.667 0.738 25.3 0.398 17 1 Introduction 2 Observations and analysis 2.1 VHE -ray observations 2.2 HE -ray observations 2.3 X-ray and optical/UV observations 2.3.1 NuSTAR 2.3.2 Swift-XRT 2.3.3 Swift-UVOT 2.4 Optical observations and polarimetry data 2.5 Radio observations 3 Multi-wavelength light curves 3.1 Analysis of radio activity in a larger time-frame 3.2 Electric Vector Position Angle swing 4 Jets evolution study with very-long-baseline interferometry 5 Spectral fitting of NuSTAR and Swift-XRT data 6 VHE Differential Energy Spectrum and EBL deabsorbtion 6.1 Redshift estimation 7 Broadband Spectral Energy Distribution 8 Summary and Conclusions A Knot parameters for S5 0716+714 (July 2014-August 2015)