Investigating the blazar TXS 0506+056 through sharp multi-wavelength eyes during 2017-2019

The blazar TXS 0506+056 got into the spotlight of the astrophysical community in September 2017, when a high-energy neutrino detected by IceCube (IceCube-170922A) was associated at the 3 $\sigma$ level to a $\gamma$-ray flare from this source. This multi-messenger photon-neutrino association remains, as per today, the most significant one ever observed. TXS 0506+056 was a poorly studied object before the IceCube-170922A event. To better characterize its broad-band emission, we organized a multi-wavelength campaign lasting 16 months (November 2017 to February 2019), covering the radio-band (Mets\"ahovi, OVRO), the optical/UV (ASAS-SN, KVA, REM, Swift/UVOT), the X-rays (Swift/XRT, NuSTAR), the high-energy $\gamma$ rays (Fermi/LAT) and the very-high-energy (VHE) $\gamma$ rays (MAGIC). In $\gamma$ rays, the behaviour of the source was significantly different from the 2017 one: MAGIC observations show the presence of flaring activity during December 2018, while the source only shows an excess at the 4$\sigma$ level during the rest of the campaign (74 hours of accumulated exposure); Fermi/LAT observations show several short (days-to-week timescale) flares, different from the long-term brightening of 2017. No significant flares are detected at lower energies. The radio light curve shows an increasing flux trend, not seen in other wavelengths. We model the multi-wavelength spectral energy distributions in a lepto-hadronic scenario, in which the hadronic emission emerges as Bethe-Heitler and pion-decay cascade in the X-rays and VHE $\gamma$ rays. According to the model presented here, the December 2018 $\gamma$-ray flare was connected to a neutrino emission that was too brief and not bright enough to be detected by current neutrino instruments.


INTRODUCTION
On September 22, 2017, the IceCube neutrino observatory (IceCube Collaboration et al. 2006) detected a 290 TeV neutrino (IceCube-170922A) from a direction consistent with the blazar TXS 0506+056, which was found to be flaring in γ rays with both Fermi /LAT and MAGIC. Fermi /LAT detected the blazar in a longlasting (about six months long) flaring state, while MAGIC detected a fast flare about 10 days after the IceCube neutrino. The chance probability that the Ice-Cube event and the flaring in the Fermi /LAT (MAGIC) energy band happened simultaneously in space and time was estimated at the 3 σ (3.5 σ) level (IceCube Collaboration et al. 2018a), and is, as per today, the strongest evidence for joint photon and neutrino emission from an active galactic nucleus (AGN). This multi-messenger photon-neutrino emission is the smoking gun attesting for the presence of highly-relativistic hadrons in AGN jets. Its firm detection would unequivocally identify AGNs as (ultra-) high-energy cosmic-ray accelerators. Motivated by the detection of IceCube-170922A, the IceCube collaboration searched for additional neutrinos coming from TXS 0506+056 and identified a cluster of events in 2014-2015, significant at the 3.5 σ level (Ice-Cube Collaboration et al. 2018b). Interestingly, this six-months-long neutrino flare was not accompanied by an enhanced electromagnetic activity, although the multi-wavelength coverage was limited due to the absence of any specific alert at that time, and relies only on information from survey instruments.
The two events, offering the first chance to test hadronic radiative models on multi-messenger data, have been extensively studied and modeled by several research groups. Standard blazar lepto-hadronic models, in which the neutrino emission is associated with the decay of pions produced in proton-photon interactions  Ansoldi et al. 2018;Padovani et al. 2018). All optical and X-ray fluxes are corrected from absorption by neutral material. MAGIC data are shown as observed, and are not corrected for pair-production absorption on the extragalactic background light. in the jet, are able to fit the electromagnetic spectral energy distribution (SED) and at the same time reproduce a neutrino rate consistent with the single event seen by IceCube in 2017. While single-zone models face the disadvantage of requiring a high proton power (often super-Eddington) to fit the data (Cerruti et al. 2019;Gao et al. 2019), solutions that make use of external target photon fields seem more promising (Ansoldi et al. 2018;Keivani et al. 2018;Petropoulou et al. 2020). Another important conclusion from the IceCube-170922A modeling efforts is that pure hadronic models, in which the γ-ray emission is dominated by proton-synchrotron radiation, are excluded.
On the other hand, the 2014-2015 neutrino-only flare challenges current blazar emission models. Pion decay injects neutrinos together with photons and electrons/positrons in the emitting region. The latter radiate synchrotron and inverse-Compton photons (via a e ± -pair cascade that transfers the power to lower frequencies). This photon emission associated with the neutrino one has to be significantly suppressed in order to reproduce the 2014/2015 event. Alternatively, that photon emission can be channeled into the MeV band, where the upper limits are less constraining (see Rodrigues et al. 2019). This might suggest the existence of a dark neutrino emitting region, opaque to photons, and physically separated from the blazar emitting region that dominates the electromagnetic SED (Reimer et al. 2019;Xue et al. 2019).
A possible solution could also be the "cosmic-rayneutral-beam" scenario, in which neutrons produced in proton-photon interactions escape the region and travel further down in the jet before interacting again and producing neutrinos, while the associated cascade emission is isotropized in the larger-scale jet (see Murase et al. 2018;Zhang et al. 2020). Another alternative is represented by interactions of the relativistic jet with obstacles, such as clouds or stars: in this case neutrinos can be efficiently produced in proton-proton interactions (see Sahakyan 2018;Banik et al. 2020) without strong ambient photon fields triggering electron-positron cascades.
Although the first candidate neutrino source is a blazar, it is important to underline that it is unlikely that γ-ray blazars represent the bulk of the diffuse neutrino background (Aartsen et al. 2017;Hooper et al. 2019). With 10 years of IceCube data, there are still no sources individually significant at more than 5 σ in the neutrino sky (Aartsen et al. 2020). The hottest spot in the all-sky scan (5 σ pre-trials, reduced to 1.3σ post-trials) is consistent with the AGN NGC 1068 (Seyfert galaxy with starburst activity); in the source catalog search, the most significant one is again NGC 1068 (4.1 σ pre-trial; 2.9 σ post-trials), while TXS 0506+056 is significant at the 3.6 σ level (pre-trial). ANTARES, a neutrino telescope located in the northern hemisphere (Ageron et al. 2011), reports 3FGL J2255.1+2411 as the most significant blazar among the F ermi/LAT sample (Albert et al. 2020), significant at the 2.3 σ level (combined space and time probability) and at the 2.6 σ level when considering information from IceCube, that also detected a high-energy neutrino consistent with this source (IceCube-100608A).
TXS 0506+056 was a very poorly studied object before 2017. The interest triggered by its association with a high-energy neutrino prompted a measurement of its redshift (z=0.337, Paiano et al. 2018), as well as several multi-wavelength follow-ups. It is remarkable that the source stands out among the blazar population as an atypical high-luminosity/high-synchrotron-peak blazar (Padovani et al. 2019), with a synchrotron peak frequency located at about 10 15 Hz (right at the divide between intermediate and high-frequency-peaked blazars, see Abdo et al. 2010) and a peak luminosity of about 10 47 erg s −1 (a value more typical of flat-spectrum radio quasars). High-resolution VLBA radio images show possible signs of deceleration of the jet and/or a spinelayer structure within 1 mas from the mm-VLBI core of the source, as well as an apparent limb brightening (Ros et al. 2020). After the 2017 multi-messenger event, follow-up observations have been conducted in the TeV band (with the MAGIC and VERITAS telescopes, see Ansoldi et al. 2018;Abeysekara et al. 2018), and in the optical-infrared band (Hwang et al. 2020).
The absence of archival observations of TXS 0506+056 (besides survey instruments such as Fermi /LAT and OVRO) means that, in 2017, it was not possible to estimate duty cycles of the source at various wavelengths, or estimate, for example, how exceptional was the γ-ray flare detected by MAGIC, or if the X-ray flux during the neutrino event was also unusually high. To significantly improve our knowledge of this unique source, we started a multi-wavelength campaign during the 2017-2019 observing season, including data from radio to very-high-energy (VHE, E > 100 GeV) γ rays. Results from this campaign are presented in this paper and organized as follows. In Section 2 we present the data analysis, from high to low frequencies ( . MAGIC spectral energy distribution as measured during this campaign (red). Low state is shown on the lef t and the flare on the right. Overlaid are the points from the low state (lef t, blue) and flares (right, blue and green) observed in 2017 (Ansoldi et al. 2018).
3 we discuss the observed multi-wavelength variability of the source; in Section 4 we present the multi-messenger modeling of the observations; discussion of the results and conclusions are in Section 5.

OBSERVATIONS AND DATA ANALYSIS
This paper is based on an extended multi-wavelength data set collected from November 2017 till February 2019, i.e. giving continuation to the observations published in Ansoldi et al. (2018). Dedicated monitoring observations were performed with Swift/XRT (Burrows et al. 2005) and Swift/UVOT (Roming et al. 2005), as well as NuSTAR (Harrison et al. 2013). They were coordinated to maximise the simultaneity of the observations with MAGIC. MAGIC observations were usually accompanied by the KVA (Nilsson et al. 2018) optical telescope, and additional optical measurements were performed with REM (Zerbi et al. 2001;Covino et al. 2004), and the University of Siena telescope. TXS 0506+056 is also systematically monitored by the ASAS-SN project (Kochanek et al. 2017) and the radio telescopes OVRO (Richards et al. 2011) and Metsähovi. Fig.1 shows the observed light curves from the VHE γ rays to the radio band. To ease the comparison with the previously published results, we include in the same plot data from September and October 2017, i.e. the multi-wavelength campaign that followed the IceCube neutrino alert (IceCube Collaboration et al. 2018a;Ansoldi et al. 2018;Padovani et al. 2018). The following sub-sections describe the data analysis for each of the instruments involved in the observing campaign.

MAGIC
The observations at VHE were performed with the Major Atmospheric Gamma-ray Imaging Cherenkov (MAGIC) stereoscopic system of telescopes, located at the Observatorio del Roque de los Muchachos, on the Canary island of La Palma, Spain. Each of the two telescopes has a 17m-diameter mirror and the system can achieve a sensitivity of 0.66% of the Crab Nebula flux above 0.2 TeV in 50 h of observations (Aleksić et al. 2016).
Between November 2017 and February 2019 a total of ∼79 hours of good quality data, at zenith angles between 5 • and 50 • , were collected. The observations were performed only in "dark" conditions (i.e. not affected by moonlight 1 ) and in wobble mode (Fomin et al. 1994). The data have been analysed using the MAGIC Analysis and Reconstruction Software (MARS, Moralejo et al. 2009;Zanin et al. 2013). When necessary, the flux values are corrected for the atmospheric extinction due to clouds and aerosols using the LIDAR system at the MAGIC site (Fruck et al. 2014). The resulting daily binned light curve of integral fluxes or flux upper limits above 90 GeV is presented in the top panel of Fig.1. The upper limits are calculated for all observations with significance below 2 σ, at the 95% confidence level, using the method of Rolke et al. (2005), and assuming a 30% systematic error on the signal efficiency, as nuisance parameter. Table 1 shows the summary of the TXS 0506+056 flux levels as measured by MAGIC: two flaring episodes are highlighted: on December 1, and 3, 2018 (MJD 58453 and 58455); the remaining period is considered as low state. The flux levels in table 1 are calculated assuming a power-law spectrum, with photon indices of 3.7 ± 0.2 for the flare and 3.5 ± 0.4 for the low state period. The typical systematic uncertainties of MAGIC spectral measurements can be divided into: < 15% in energy scale, 11-18% in flux normalization and ±0.15 for the energy spectrum power-law slope Aleksić et al. (2016). The blazar remained in a low state during most of the monitored period (∼74 h). The average flux F(>90 GeV) = (2.7±2.1)×10 −11 cm −2 s −1 , calculated from the excess at a significance of 4 σ, is the lowest VHE γ-ray emission level observed from this source so far. Most notably, the difference between low state and high state integral photon flux spans an order of magnitude. The measured high state fluxes are comparable with the flare detected by MAGIC in October 2017, shortly after the neutrino alert IceCube-170922A.
The similarities between the 2018 and 2017 high states are also evident in the MAGIC spectral energy distributions shown in Fig. 2 (right). Both the photon index and flux levels are compatible within errors with the previous measurements.

Fermi/LAT
In this work, the publicly available F ermi/LAT (Atwood et al. 2009) data collected in the period from November 1, 2017, to March 1, 2019 (531187205 to 573091205 MET), are used. The Pass 8 SOURCE data were analysed with Fermi Science-Tools (1.2.1) with P8R3 SOURCE V2 instrument response functions. The 100 MeV -300 GeV events from a circular region of interest of 12 • around the γ-ray position of TXS 0506+056 (RA, Dec)= (77.359, 5.701) were downloaded and analyzed using the standard procedure described by the F ermi/LAT Collaboration. The data are binned into pixels of 0.1 • × 0.1 • and standard cuts were applied to select the good time intervals ("DATA QUAL>0" and "LAT CONFIG==1", z max < 90 • ). The model file was created based on the F ermi/LAT fourth source catalog (4FGL, Abdollahi et al. 2020) and it includes all point sources within 17 • from TXS 0506+056 and standard templates for Galactic diffuse emission model (gll iem v07) and isotropic diffuse emission (iso P8R3 SOURCE V2 v1).
The γ-ray light curves were calculated by performing an unbinned maximum likelihood analysis with the appropriate quality cuts as described above. The photon indexes of all sources, except TXS 0506+056, as well as the normalization of both background components were fixed to the best-fit values obtained for the whole time period. Only the normalizations of the sources within 12 • from TXS 0506+056 were left free during the analysis. Initially the light curve was computed using 7-day intervals and then using the adaptive binning method (Lott et al. 2012). In the latter case, the time bin widths are flexibly adapted to produce bins with constant flux uncertainty above the optimal energies (E 0 ). The adaptively binned light curve computed for 20% uncertainty and E 0 = 216.03 MeV is shown in Fig. 1. The γ-ray flux and photon index variations are investigated using the 7-day binned light curve ( Fig. 5).
In the low state, the γ-ray spectrum of the TXS 0506+056 has been computed by limiting the analysis only to the periods contemporaneous with the MAGIC low state, and the data observed during MJD 58451-58456 (November 29, to December 4, 2018) have been used for the active state. During the analysis, the spectrum of TXS 0506+056 has been modeled by a power-law function considering the normalization and index as free parameters. The γ-ray emission spectra in the low and active states are shown in Fig. 6.

NuSTAR
NuSTAR (Harrison et al. 2013) observed TXS 0506+056 with its two co-aligned X-ray telescopes with corresponding focal planes, focal plane module A (FPMA) and B (FPMB), five times between April 3, 2018, and January 7, 2019, for an exposure time varying between 21.3 ks and 31.7 ks (see Table 4). The level 1 data products were processed with the NuSTAR Data Analysis Software (nustardas) package (v2.0.0). Cleaned event files (level 2 data products) were pro-duced and calibrated using standard filtering criteria with the NUPIPELINE task and version 20201101 of the calibration files available in the NuSTAR CALDB and the OPTIMIZED parameter for the exclusion of the South Atlantic Anomaly passages. Spectra of the sources were extracted from the cleaned event files using a circle of 30 -radius, while the background was extracted from a nearby circular regions of 70 -radius on the same chip of the source. The ancillary response files were generated with the numkarf task, applying corrections for the point spread function losses, exposure maps and vignetting. The spectra were rebinned with a minimum of 20 counts per energy bin to allow for χ 2 spectrum fitting. All errors are given at the 90% confidence level.
The NuSTAR spectra have been fitted in the 3.0-78 keV energy range with an absorbed power-law with the Galactic absorption corresponding to a hydrogen column density of n H = 1.11 × 10 21 cm −2 (Kalberla et al. 2005). The results of the fits are presented in Table 4. The 3-78 keV flux varied during the NuSTAR monitoring up to a factor of 2, in particular from (4.54 ± 0.34) × 10 −12 to (7.15 ± 0.39) × 10 −12 erg cm −2 s −1 between December 8, 2018, and January 7, 2019. No significant variation of the photon index has been observed in the monitored period.

Swift/XRT
The Neil Gehrels Swift Observatory (Swift) has observed the source 41 times between November 8, 2017, and December 11, 2018. All XRT (Burrows et al. 2005) observations were performed in photon counting mode. The XRT spectra were generated with the Swift XRT data products generator tool at the UK Swift Science Data Centre 2 (for details see Evans et al. 2009). The X-ray spectra in the 0.3-10 keV energy range is fitted by an absorbed power-law model using the photoelectric absorption model tbabs (Wilms et al. 2000) with a HI column density consistent with the Galactic value in the direction of the source, as reported in Kalberla et al. (2005), i.e. 1.11 × 10 21 cm −2 . A large number of spectra show low number of counts (i.e. < 200), therefore not allowing us to use the χ 2 statistics. To maintain the homogeneity in the analysis, the spectra are grouped using the task grppha to have at least one count per bin and the fit has been performed with the Cash statistics (Cash 1979). We used the spec-2 http://www.swift.ac.uk/user objects tral redistribution matrices in the Calibration database maintained by HEASARC. The X-ray spectral analysis was performed using the XSPEC 12.9.1 software package (Arnaud 1996). Single observations with a number of counts < 15, for which it is not possible to constrain the spectral parameters, are not considered. The X-ray spectrum is variable with a photon index, Γ X , varying between 1.27 ± 0.37 and 2.71 ± 0.70.

Swift/UVOT
The UVOT telescope (Roming et al. 2005) on board the Swift satellite observed TXS 0506+056 simultaneously with XRT, guaranteeing a coverage in the optical and UV band. Fluxes for each of the six UVOT filters (V, B, U, UW1, UM2, UW2) are calculated with the Heasoft tool uvotmaghist, using a 5 -radius circular region for the source, and a 20 -radius circular region for the background. The data have been analyzed with the most recent (September 2020) UVOT calibration file. 3 The flux measurements are corrected for absorption assuming E B−V = 0.092 as value for the Galactic extinction (Schlafly & Finkbeiner 2011). To study the spectral shape in the optical/UV band, and its evolution with time, a power-law fit to the UVOT flux points has been performed for all observations including at least four different filters. All UVOT spectra are characterized by an index softer than 2.0, indicating that the peak of the synchrotron emission is located below the V band. All power-law indices from UVOT observations are shown in Fig. 4. 2.6. Optical telescopes TXS 0506+056 is automatically monitored by the ASAS-SN project (Kochanek et al. 2017), and the optical light curve has been extracted using the automatic online tool. 4 Figure 1 shows the g and V band light curve extracted with an aperture radius of 16 .
The KVA measurements were performed contemporaneously with MAGIC for most of the nights, and the data were analysed using the standard procedure, with a signal extraction aperture of 5 in the R-band (Nilsson et al. 2018).
The Rapid Eye Mount telescope (REM) observations were triggered after the MAGIC flare. REM is a 60-cm robotic telescope located at the ESO La Silla Observatory (Zerbi et al. 2001;Covino et al. 2004). It includes an optical camera with the Sloan filters g, r, i, z and a near-infrared camera equipped with J-H-K filters. Data reduction was carried out following the standard procedures, with the subtraction of an averaged bias frame dividing by the normalised flat frame. The photometric calibration was achieved by using the APASS catalogue. In order to minimise any systematic effect, we performed differential photometry with respect to a selection of non-saturated reference stars selecting the signal around the centroid of a source within a circle of 10 -radius. The monitoring of REM went on for 7 days in r band and the magnitudes observed, converted from the SDSS to Johnson-Cousins photometric system, are comparable with the other measurements in the R-band.
The Astronomical Observatory of the University of Siena observed TXS 0506+056 in the context of a program devoted to optical photometry of blazars as a follow-up of MAGIC observing campaigns. The instrumentation consists of a remotely operated 30 cm f/5.6 Maksutov-Cassegrain telescope installed on a Comec 10 micron GM2000-QCI German equatorial mount. The detector is a Sbig STL-6303 camera equipped with a 3072 x 2048 pixels KAF-6303E sensor; the filter wheel hosts a set of Johnson-Cousins BVRI filters. Multiple 300s images of TXS 0506+056 were acquired in the R band at each visit. After standard dark current subtraction and flat-fielding, images for each visit were averaged and aperture photometry was performed on the average frame by means of the MaximDL software package 5 , extracting the signal within 7 around the source centroids. The choice of reference and control stars was consistent with the one for the KVA data. The obtained magnitudes have been treated as indicated for the KVA data, as far as correction of galactic extinction, conversion to flux and subtraction of the host galaxy contribution to the measured flux.
In order to reproduce the light curve (Fig. 1), all the optical data are corrected for absorption assuming E B−V = 0.092.

Metsähovi
The 37 GHz observations were performed with the 13.7 m diameter Metsähovi radio telescope. 6 A typical integration time to obtain one flux density data point 5 https://diffractionlimited.com/product/maxim-dl/ 6 https://www.aalto.fi/en/services/metsahovis-main-instruments is between 1200 and 1800 s. The detection limit of our telescope at 37 GHz is on the order of 0.2 Jy under optimal conditions. 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. Sources NGC 7027, 3C 274 and 3C 84 are used as secondary calibrators. 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 contribution from the measurement rms and the uncertainty of the absolute calibration.

OVRO
The 15 GHz data presented here were taken in the frame of the large-scale, fast-cadence monitoring program with the 40 m telescope at the Owens Valley Radio Observatory (OVRO). The data analysis was performed according to a standard analysis described in Richards et al. (2011). The flux density of the source increases from about 0.6 Jy in November 2017 to more than twice this value in February 2019. We also note that an increase in the flux was also seen in the 2016 data, well before the neutrino event in September 2017 (Hovatta et al. 2021).

VARIABILITY STUDY
AGNs are known for their variability and TXS 0506+056 is no exception. Fig. 1 shows that a clear variability is observed in all wavelengths. We use the fractional variability parameter F var as defined by Vaughan et al. (2003) to quantify the observed trends. The different sensitivity and flux sampling of the differ-ent instruments introduce some biases in this method (see e.g. Aleksić et al. 2014). On the other hand, it is a relatively simple way to quantify and compare the flux variability among the observed energy bands. Table 5 and Fig 3 present the F var values for all the instruments contributing to this study, apart from MAGIC. The very low VHE γ-ray activity during these two years yielded low significance flux measurements for most of the single night observations. Therefore the fractional variability in this energy band cannot be estimated.
The most reliable variability estimates can be done using the data from monitoring instruments such as F ermi/LAT and ASAS-SN. A relatively large data set was also collected with the X-ray and UV instruments on board Swif t and the KVA optical telescope. High sensitivity NuSTAR observations complete the multi-wavelength picture. The most pronounced variability F var ∼0.40 is observed in the F ermi/LAT γ-ray band, while the X-ray variability is at a lower level of ∼0.25-0.30. The radio, optical and UV bands display a moderate variability of F var ∼0.15-0.22. A slight spread between the KVA R-band data and the ASAS-SN and UVOT observations at shorter wavelengths is visible. This is probably due to sparser KVA observations, mostly covering high emission states. The radio variability is at a marginally higher level than the optical and UV ones, which is a quite common blazar characteristic.
The VHE γ-ray flare observed in December 2018 is very similar to the one observed in October 2017, both in term of the flux level, as well as day-scale variability. In a matter of a few days, we see an order of magnitude flux increment and then a decay.

Spectral variability
An important aspect to investigate when evaluating blazar variability is the dependence of the synchrotron peak energy or the high-energy peak energy on the flux state. Very often in γ-ray blazars, a harder-whenbrighter behaviour, meaning a harder spectrum during enhanced flux states, is observed in both X-and γ-ray bands. This trend is related to the fact that the enhanced flux state can be explained with an injection of high-energy particles, shifting the peaks of the SED towards higher energies and higher flux levels (MAGIC Collaboration et al. 2020). However, some exceptions to this trend are reported in literature (e.g., Liodakis et al. 2018). Table 2. Results of the Spearman correlation study between the photon index and the flux level.

Band
R p-value UV (V band) -0.17 0.19 X-ray (2-10 keV) -0.90 2 × 10 −15 X-ray (0.3-10 keV) -0.23 0.14 X-ray (3-78 keV) -0.25 0.59 γ-ray (0.2-300 GeV) 0.36 0.006 To investigate if TXS 0506+056 data in different bands follow this trend, a correlation study between the photon flux and the photon index was performed, using the Spearman coefficient R. The results of our study are reported in Table 2. The last column of the table reports the p-value of the hypothesis of no correlation: smaller p-values indicate a higher probability of (anti)correlation.
In our analysis, no correlation is found in UVOT data, Figure 4, left. The corresponding Spearman coefficient is R = −0.17, and the p-value is 0.19. In soft X-ray, Figure 4 right panel, a strong anti-correlation between the photon flux in the 2-10 keV band and the photon index emerges, as supported by a value of R = −0.90 (p-value=2 × 10 −15 ). This correlation is lost when considering data in the full XRT energy range, 0.3-10 keV, gray markers in Figure 4, right panel. No significant spectral variability has been observed in hard X-rays during the NuSTAR monitoring of the source.
At higher energies, only Fermi/LAT data were considered for the study, as the faintness of the signal in the MAGIC energy range prevented a detailed study in this band. The flux-vs-index correlation plot obtained with the Fermi/LAT data, averaged over 7-days bins, is displayed in Fig. 5. A linear fit to the data suggests a correlation between the two values, indicating a softerwhen-brighter trend in γ rays. The Spearman coefficient value, however, is only R = 0.36 (p-value = 0.006), so with the available data set we can only conclude that a hint of correlation is deducible in the GeV band. gray, the whole range of optical-to-X-ray observations). For the VHE high-state, the only instrument performing contemporaneous observations with MAGIC were F ermi/LAT and ASAS-SN. We include nonetheless the nearest multi-wavelength observations in radio, optical, UV and X-rays (on December 8, 2018). Blazar SEDs can be fit by both leptonic and hadronic models. Neutrino detection from a blazar can indeed be the key to discriminate between the two emission scenarios. Fully aware that in absence of neutrinos the SED of TXS 0506+056 can be solely described by leptons, we model the SEDs in a lepto-hadronic scenario inspired by the September 2017 photon-neutrino event: the dominant radiative process in the γ-ray band is external-inverse-Compton, with hadronic radiative components (Bethe-Heitler and pion-decay pair-cascade) emerging in the X-ray and VHE bands. The dominant low-energy photon field for proton-photon interactions is synchrotron radiation from the jet layer that encompasses the jet spine in which the plasmoid moves (Ghisellini et al. 2005). This geometry is also supported by recent spatially resolved radio observations described in Ros et al. (2020).
The numerical simulation is performed with the code described in Cerruti et al. (2015), expanded to allow for hadronic interactions over arbitrary external photon fields. Photo-meson interactions in the code are implemented via the Monte-Carlo code SOPHIA (Mücke et al. 2000), while Bethe-Heitler pair-production is implemented following Kelner & Aharonian (2008). Absorption of VHE γ rays on the extragalactic background light is implemented using the model by Franceschini et al. (2008). The pair-cascades associated with hadronic interactions are computed iteratively generation-bygeneration. This assumes that the cascades are not self-supported, and that the dominant photon fields are either the external target photon field, or the synchrotron emission by primary electrons, but never the photon emission from the cascades themselves. The radiative code computes the steady-state emission from a spherical plasmoid in the jet, filled by a tangled, homogeneous magnetic field, and interacting with external photons from the jet layer. The plasmoid is fully characterized by three free parameters: the bulk Lorentz factor Γ (translated into the Doppler factor δ assuming an angle to the line of sight θ view ), the radius R , and the magnetic flux density B . 7 The emitting region is filled with a population of primary electrons and protons, without any assumption on the underlying acceleration mechanism. The primary particle distributions are parametrized by broken-power-law functions with exponential cut-offs at their maximum Lorentz factor γ Max , and are thus characterized by six free parameters each (the three Lorentz factors γ min , γ break , γ Max , the two indexes α 1 , and α 2 , and the normalization K ), for a total of nine free parameters. In addition, the external photon field from the jet layer carries other eight free parameters: the layer region is a hollow cylinder (parametrized by its inner and outer radii R in,layer and R out,layer , and its height H layer ), moving with bulk Lorentz factor Γ layer and filled with a homogeneous magnetic field B layer , and an electron population parametrized by a broken power-law function with its six free parameters. Such a model is clearly degenerate, and several assumptions have to be introduced to reduce the number of free parameters: the primary electrons and protons are considered to be co-accelerated, sharing the same index α 1,e = α 1,p ; the primary protons are not cooled (i.e. all cooling processes become relevant at Lorentz factors higher than γ p,Max ) and thus the proton distribution is parametrized by a simple power-law function; the minimum proton Lorentz factor is fixed to 1, to avoid biasing the total luminosity estimate; the Doppler factor of the plasmoid is fixed to a rather typical value of 40 (Tavecchio et al. 2010;Zhang et al. 2012), with an angle to the line of sight θ view = 0.8 • ; the geometry of the layer is R in,layer = R , R out,layer = 1.5R , while the height of the layer in the frame of the spine is H layer = H layer /(ΓΓ layer (1 − ββ layer )) = R (see Ghisellini et al. 2005); the photon field from the layer is computed using δ layer = 4, B layer = 1.8 G, α e,1,layer = 2.0, α e,2,layer = 2.9, γ e,min,layer = 290, γ e,break,layer = 2000, and γ e,Max,layer = 3.9 × 10 5 (identical to the values used in Ansoldi et al. 2018). The normalization of the synchrotron emission from the layer is left free to vary and is expressed as photon energy density in the reference frame of the plasmoid, u ph,layer . The external inverse-Compton component is corrected by a factor δ/δ layer following Dermer (1995). These physical constraints reduce the number of free parameters of the model to eleven: R , B , γ e,min , γ e,break , γ e,Max , γ p,Max , α 1 , α 2,e , K e , K p , and u ph,layer . This number remains larger than the number of independent observables, so no numerical fitting of the model to the data is performed. We limit 7 Here and in the following, primed quantities are given in the reference frame of the plasmoid. Double primed quantities are in the reference frame of the layer.
ourselves to identify a lepto-hadronic solution (among many) that can reproduce the observations.
The multi-wavelength observations and the variability studies presented in the previous sections provide further constraints to the modeling. The position of the synchrotron peak is located in the optical band or lower, as determined by UVOT. The X-ray band marks the transition between two different emission components, as clearly shown by the different spectral variability properties in Swift/XRT and NuSTAR data. In the lepto-hadronic scenario, this is naturally explained by the transition from the synchrotron component to the inverse-Compton one, with additional contribution from the Bethe-Heitler pair cascade. We start by studying the low-state emission of TXS 0506+056. The free parameters of the model are first optimized to get a leptonic modeling of the optical-UV and γ-ray part of the SED. The proton population is then increased until the emission from Bethe-Heitler and photo-meson cascades starts violating X-ray and VHE observations. An additional constraint comes from IceCube observations, that identify a hotspot at the position of TXS 0506+056 corresponding to 12.3 events in ten years. We compute the expected neutrino rate by convolving the neutrino spectrum from the model with the IceCube effective area for point-source searches. The lepto-hadronic solution, shown in Figure 6 (with model parameters provided in Table 3), has a hadronic contribution emerging in the X-rays and at VHE as Bethe-Heitler and pion-decay cascade. The model for the low state predicts an Ice-Cube neutrino rate of 1.4 events in ten years (estimated using the IceCube point-source effective area), much less than the recent results by Aartsen et al. (2020, 12.3 neutrinos in ten years). This supports the idea that the hotspot at the position of TXS 0506+056 in IceCube data is dominated by flaring events (i.e. the flare which occurred in 2014-15 and the high-energy neutrino of September 2017).
It is important to underline that the hadronic model remains degenerate unless explicit assumptions on the acceleration/cooling processes (and thus the maximum proton Lorentz factor γ p,Max ) are made. The maximum proton energy is related to the peak energy of the neutrino emission: lower values of γ p,Max would imply a higher neutrino rate (because the spectrum would then peak closer to the maximum of IceCube sensitivity) at the expense of a higher total jet power. In an opposite way, the total jet power (in our model equal to 7.6 × 10 47 erg s −1 , or 7.6 − 76 L Edd for a super-massive black-hole with mass M • = 10 9−8 M ) can be reduced if the maximum proton energy is increased, resulting in a shift towards higher energies of the neutrino peak. The total jet power is very dependent on the index of the proton distribution and on γ p,min . The model presented here, with γ p,min = 1 and α p = 2.0 is explicitly conservative in this sense: a harder particle index, or γ p,min 1 will significantly lower the total jet power.
Once we obtained a satisfactory modeling of the lowstate emission of TXS 0506+056, we investigated the December 2018 VHE flare. In this case the only truly simultaneous observations are from MAGIC, Fermi /LAT, and ASAS-SN: the behaviour of the synchrotron component of the SED during the flare is unknown. On December 3, 2018 the MAGIC collaboration issued an alert (Mirzoyan 2018) to encourage further multiwavelength observations of TXS 0506+056. Together with the X-ray and optical instruments, both IceCube and ANTARES set upper limits on the neutrino emission in the week/10 days around the VHE γ-ray flare (Coleiro & Dornic 2018;Vandenbroucke 2018). They are also shown on the SED plot. The model shown in Figure 6 (with model parameters in Table 3) is obtained through a minimal variation of the better constrained low-state model: the only parameters that change are the ones from the primary particle populations: the electron break Lorentz factor is increased to 10 4 , while the maximum Lorentz factor is increased to 2 × 10 4 (while maintaining the same electron energy density); the primary proton distribution is unchanged in spectral shape, but increased in overall normalization by about a factor of two. The model for the γ-ray flaring state predicts the presence of an optical/UV/soft-X-ray flare that was missed due to absence of simultaneous observations in that band. We highlight that this statement is very model-dependent, and there exist other solutions (obtained assuming that other model parameters are driving the variability) with a very different behaviour in the synchrotron and neutrino components. In this scenario the expected neutrino rate in IceCube data (computed using the GOLD event effective area, Blaufuss et al. 2019) is 2.1 × 10 −4 events in a day, or 1.5 × 10 −3 in a week, completely consistent with the non-detection of neutrinos associated with this γ-ray flare.

SUMMARY AND CONCLUSIONS
We have performed an extensive multi-wavelength campaign from radio to VHE γ rays on the "neutrino blazar" candidate TXS 0506+056. The observations cover the time span from November 2017 (right after the neutrino event of September 2017) to February 2019. The source has been very poorly studied before 2017, and this work represents the most complete characterization of the behaviour of TXS 0506+056 on a Note-The quantities flagged with a star ( ) are derived quantities, and not model parameters.
The luminosity of the emitting region has been calculated as L = 2πR 2 cΓ 2 (u B + u e + u p ), where u B , u e , and u p are the energy densities of the magnetic field, the electrons, and the protons, respectively. For the high-state model, only the parameters that change with respect to the low-state one are indicated.

time-scale of a few years.
In the VHE band, the source showed another daylong flare on December 1-3, 2018, observationally very similar (both in flux level and photon index) to the one of September 2017. In the rest of the multi-wavelength campaign, TXS 0506+056 was not detected neither on single days, nor in short periods of time. When aggregating all data collected with MAGIC, which amounts to 75 hours, one gets a VHE γ-ray excess at the level of 4 σ. No neutrino flares associated with the VHE flare in December 2018 have been reported. In the HE γ-ray band, TXS 0506+056 was in a very different state compared to that from September 2017: while at that time the source was in a long-term (month-scale) brightening, during this campaign the source showed a much lower average state with much faster (days-to-week) flaring activity. No significant flares are detected in X-rays or optical/UV. In the radio band, OVRO observations at 15 GHz revealed the presence of a long-term brightening of the source which increased its flux by around a factor of two throughout the campaign. The same brightening is seen by Metsähovi at 37 GHz, with the same approximate two-fold increase in flux during the discussed time period. While the radio behaviour seems uncorrelated with the other bands on the time-scale of the observing campaign, it is important to underline that the physical zones sampled by radio observations are different from the ones at higher energies, and increases in radio flux densities have been associated to γ-ray flaring activity (Jorstad et al. 2001;Lähteenmäki & Valtaoja 2003;León-Tavares et al. 2011), and, tentatively, with neutrinos (Hovatta et al. 2021). Spectral studies reveal only marginal spectral variability in UVOT and F ermi/LAT, indicating that the position of the SED peaks was stable during the campaign. On the other hand, a very interesting spectral variability is seen in Xrays: while the hard-X-ray band (3-70 keV, sampled by NuSTAR) shows no spectral variability, the soft-X-rays are much more variable, with a clear harder-whenbrighter behaviour detected in the 2-10 keV band with XRT observations. When studying the whole XRT band (0.3-10 keV), the anti-correlation between the photon index and the flux is lost, suggesting the presence of a transition between different radiation mechanisms in the keV band, that in a lepto-hadronic model corresponds to the transition from synchrotron radiation by primary electrons, through the one by Bethe-Heitler pairs, to inverse-Compton.
We perform a lepto-hadronic modeling of TXS 0506+056, starting from the low-state of the source which is now, for the first time, very well characterized. A scenario in which the emission is dominated by the leptonic component (inverse-Compton in the γ-ray band), with hadronic interactions happening between protons in the blazar region and photons from the jet-layer, can provide a good description of the SED. The model predicts a detection rate by IceCube of 0.14 neutrinos per year, or 1.4 neutrinos in ten years of observations, consistent with the most recent estimates by IceCube (12.3 neutrinos in ten years). This modeling strengthens the idea that the neutrino emission from AGNs is dominated by their flaring events. The VHE flare of December 2018 does not have simultaneous observations at lower energies, which prevents us to efficiently constrain the theoretical model. We present a tentative modeling of these observations assuming that the only parameters that changed with respect to the low-state of the source were the electron and proton primary distributions. The conclusion is that the γ-ray flare had to be accompanied by a brightening in the optical/UV/X-ray band. Within this model, the associated neutrino emission would have been too rapid and not bright enough to be detected by IceCube (predicted rate of 2.1 × 10 −4 neutrinos per day).
The association of the neutrino IceCube-170922A with the flaring blazar TXS 0506+056, has promoted this once vanilla blazar to one of the most important γ-ray sources in the sky. To fully understand the multi-messenger gamma-neutrino association, a thorough characterization of the source is needed. In this work we presented the largest multi-wavelength campaign on TXS 0506+056 to date, with 16 months of data-taking from radio to VHE, but future observations are definitely needed. It is of particular importance the characterization of the variability properties of the source and its duty cycle in the various energy bands, to ultimately quantify how exceptional the 2017 γ-ray flare was. Several questions remain open, and new one have arisen from this campaign. While the source's behaviour in the GeV band changed significantly from September 2017, we detected a new VHE flare very similar to the previous one. Another highlight of this campaign is the radio behaviour of TXS 0506+056, which is undergoing a long-term brightening completely uncorrelated to any other wavelength. Future multi-wavelength and multimessenger observations of TXS 0506+056 will help cast light on hadronic acceleration in AGNs jets, and on the role of AGNs as cosmic rays and neutrino sources.