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

V. A. Acciari Instituto de Astrofísica de Canarias and Dpto. de Astrofísica,
Universidad de La Laguna, E-38200, La Laguna, Tenerife, Spain  T. Aniello
National Institute for Astrophysics (INAF), I-00136 Rome, Italy  S. Ansoldi
Università di Udine and INFN Trieste, I-33100 Udine, Italy  also at
International Center for Relativistic Astrophysics (ICRA), Rome, Italy  L. A.
Antonelli National Institute for Astrophysics (INAF), I-00136 Rome, Italy  A.
Arbet Engels Max-Planck-Institut für Physik, D-80805 München, Germany  M.
Artero Institut de Física d’Altes Energies (IFAE), The Barcelona Institute of
Science and Technology (BIST), E-08193 Bellaterra (Barcelona), Spain  K. Asano
Japanese MAGIC Group: Institute for Cosmic Ray Research (ICRR), The University
of Tokyo, Kashiwa, 277-8582 Chiba, Japan  D. Baack Technische Universität
Dortmund, D-44221 Dortmund, Germany  A. Babić Croatian MAGIC Group: University
of Zagreb, Faculty of Electrical Engineering and Computing (FER), 10000
Zagreb, Croatia  A. Baquero IPARCOS Institute and EMFTEL Department,
Universidad Complutense de Madrid, E-28040 Madrid, Spain  U. Barres de Almeida
Centro Brasileiro de Pesquisas Físicas (CBPF), 22290-180 URCA, Rio de Janeiro
(RJ), Brazil  J. A. Barrio IPARCOS Institute and EMFTEL Department,
Universidad Complutense de Madrid, E-28040 Madrid, Spain  I. Batković
Università di Padova and INFN, I-35131 Padova, Italy  J. Becerra González
Instituto de Astrofísica de Canarias and Dpto. de Astrofísica, Universidad de
La Laguna, E-38200, La Laguna, Tenerife, Spain  W. Bednarek University of
Lodz, Faculty of Physics and Applied Informatics, Department of Astrophysics,
90-236 Lodz, Poland  E. Bernardini Università di Padova and INFN, I-35131
Padova, Italy  M. Bernardos  Università di Padova and INFN, I-35131 Padova,
Italy  A. Berti Max-Planck-Institut für Physik, D-80805 München, Germany  J.
Besenrieder  Max-Planck-Institut für Physik, D-80805 München, Germany  W.
Bhattacharyya Deutsches Elektronen-Synchrotron (DESY), D-15738 Zeuthen,
Germany  C. Bigongiari National Institute for Astrophysics (INAF), I-00136
Rome, Italy  A. Biland ETH Zürich, CH-8093 Zürich, Switzerland  O. Blanch
Institut de Física d’Altes Energies (IFAE), The Barcelona Institute of Science
and Technology (BIST), E-08193 Bellaterra (Barcelona), Spain  H. Bökenkamp
Technische Universität Dortmund, D-44221 Dortmund, Germany  G. Bonnoli
Instituto de Astrofísica de Andalucía-CSIC, Glorieta de la Astronomía s/n,
18008, Granada, Spain  Ž. Bošnjak Croatian MAGIC Group: University of Zagreb,
Faculty of Electrical Engineering and Computing (FER), 10000 Zagreb, Croatia
G. Busetto Università di Padova and INFN, I-35131 Padova, Italy  R. Carosi
Università di Pisa and INFN Pisa, I-56126 Pisa, Italy  G. Ceribella Japanese
MAGIC Group: Institute for Cosmic Ray Research (ICRR), The University of
Tokyo, Kashiwa, 277-8582 Chiba, Japan  M. Cerruti⋆ Universitat de Barcelona,
ICCUB, IEEC-UB, E-08028 Barcelona, Spain  now at Université de Paris, CNRS,
Astroparticule et Cosmologie, F-75013 Paris, France  Y. Chai Max-Planck-
Institut für Physik, D-80805 München, Germany  A. Chilingarian Armenian MAGIC
Group: A. Alikhanyan National Science Laboratory, 0036 Yerevan, Armenia  S.
Cikota  Croatian MAGIC Group: University of Zagreb, Faculty of Electrical
Engineering and Computing (FER), 10000 Zagreb, Croatia  E. Colombo Instituto
de Astrofísica de Canarias and Dpto. de Astrofísica, Universidad de La Laguna,
E-38200, La Laguna, Tenerife, Spain  J. L. Contreras IPARCOS Institute and
EMFTEL Department, Universidad Complutense de Madrid, E-28040 Madrid, Spain
J. Cortina Centro de Investigaciones Energéticas, Medioambientales y
Tecnológicas, E-28040 Madrid, Spain  S. Covino National Institute for
Astrophysics (INAF), I-00136 Rome, Italy  G. D’Amico Max-Planck-Institut für
Physik, D-80805 München, Germany  now at Department for Physics and
Technology, University of Bergen, NO-5020, Norway  V. D’Elia National
Institute for Astrophysics (INAF), I-00136 Rome, Italy  P. Da Vela Università
di Pisa and INFN Pisa, I-56126 Pisa, Italy  now at University of Innsbruck  F.
Dazzi National Institute for Astrophysics (INAF), I-00136 Rome, Italy  A. De
Angelis Università di Padova and INFN, I-35131 Padova, Italy  B. De Lotto
Università di Udine and INFN Trieste, I-33100 Udine, Italy  A. Del Popolo INFN
MAGIC Group: INFN Sezione di Catania and Dipartimento di Fisica e Astronomia,
University of Catania, I-95123 Catania, Italy  M. Delfino Institut de Física
d’Altes Energies (IFAE), The Barcelona Institute of Science and Technology
(BIST), E-08193 Bellaterra (Barcelona), Spain  also at Port d’Informaciò
Científica (PIC), E-08193 Bellaterra (Barcelona), Spain  J. Delgado Institut
de Física d’Altes Energies (IFAE), The Barcelona Institute of Science and
Technology (BIST), E-08193 Bellaterra (Barcelona), Spain  also at Port
d’Informaciò Científica (PIC), E-08193 Bellaterra (Barcelona), Spain  C.
Delgado Mendez Centro de Investigaciones Energéticas, Medioambientales y
Tecnológicas, E-28040 Madrid, Spain  D. Depaoli INFN MAGIC Group: INFN Sezione
di Torino and Università degli Studi di Torino, I-10125 Torino, Italy  F. Di
Pierro INFN MAGIC Group: INFN Sezione di Torino and Università degli Studi di
Torino, I-10125 Torino, Italy  L. Di Venere INFN MAGIC Group: INFN Sezione di
Bari and Dipartimento Interateneo di Fisica dell’Università e del Politecnico
di Bari, I-70125 Bari, Italy  E. Do Souto Espiñeira Institut de Física d’Altes
Energies (IFAE), The Barcelona Institute of Science and Technology (BIST),
E-08193 Bellaterra (Barcelona), Spain  D. Dominis Prester Croatian MAGIC
Group: University of Rijeka, Department of Physics, 51000 Rijeka, Croatia  A.
Donini Università di Udine and INFN Trieste, I-33100 Udine, Italy  D. Dorner
Universität Würzburg, D-97074 Würzburg, Germany  M. Doro Università di Padova
and INFN, I-35131 Padova, Italy  D. Elsaesser Technische Universität Dortmund,
D-44221 Dortmund, Germany  V. Fallah Ramazani Finnish MAGIC Group: Finnish
Centre for Astronomy with ESO, University of Turku, FI-20014 Turku, Finland
now at Ruhr-Universität Bochum, Fakultät für Physik und Astronomie,
Astronomisches Institut (AIRUB), 44801 Bochum, Germany  L. Fariña Institut de
Física d’Altes Energies (IFAE), The Barcelona Institute of Science and
Technology (BIST), E-08193 Bellaterra (Barcelona), Spain  A. Fattorini
Technische Universität Dortmund, D-44221 Dortmund, Germany  L. Font
Departament de Física, and CERES-IEEC, Universitat Autònoma de Barcelona,
E-08193 Bellaterra, Spain  C. Fruck Max-Planck-Institut für Physik, D-80805
München, Germany  S. Fukami ETH Zürich, CH-8093 Zürich, Switzerland  Y.
Fukazawa Japanese MAGIC Group: Physics Program, Graduate School of Advanced
Science and Engineering, Hiroshima University, 739-8526 Hiroshima, Japan  R.
J. García López Instituto de Astrofísica de Canarias and Dpto. de Astrofísica,
Universidad de La Laguna, E-38200, La Laguna, Tenerife, Spain  M. Garczarczyk
Deutsches Elektronen-Synchrotron (DESY), D-15738 Zeuthen, Germany  S.
Gasparyan  Armenian MAGIC Group: ICRANet-Armenia at NAS RA, 0019 Yerevan,
Armenia  M. Gaug Departament de Física, and CERES-IEEC, Universitat Autònoma
de Barcelona, E-08193 Bellaterra, Spain  N. Giglietto INFN MAGIC Group: INFN
Sezione di Bari and Dipartimento Interateneo di Fisica dell’Università e del
Politecnico di Bari, I-70125 Bari, Italy  F. Giordano INFN MAGIC Group: INFN
Sezione di Bari and Dipartimento Interateneo di Fisica dell’Università e del
Politecnico di Bari, I-70125 Bari, Italy  P. Gliwny University of Lodz,
Faculty of Physics and Applied Informatics, Department of Astrophysics, 90-236
Lodz, Poland  N. Godinović Croatian MAGIC Group: University of Split, Faculty
of Electrical Engineering, Mechanical Engineering and Naval Architecture
(FESB), 21000 Split, Croatia  J. G. Green Max-Planck-Institut für Physik,
D-80805 München, Germany  D. Green Max-Planck-Institut für Physik, D-80805
München, Germany  D. Hadasch Japanese MAGIC Group: Institute for Cosmic Ray
Research (ICRR), The University of Tokyo, Kashiwa, 277-8582 Chiba, Japan  A.
Hahn Max-Planck-Institut für Physik, D-80805 München, Germany  T. Hassan
Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas,
E-28040 Madrid, Spain  L. Heckmann Max-Planck-Institut für Physik, D-80805
München, Germany  J. Herrera Instituto de Astrofísica de Canarias and Dpto. de
Astrofísica, Universidad de La Laguna, E-38200, La Laguna, Tenerife, Spain  J.
Hoang IPARCOS Institute and EMFTEL Department, Universidad Complutense de
Madrid, E-28040 Madrid, Spain  now at Department of Astronomy, University of
California Berkeley, Berkeley CA 94720  D. Hrupec Croatian MAGIC Group: Josip
Juraj Strossmayer University of Osijek, Department of Physics, 31000 Osijek,
Croatia  M. Hütten Japanese MAGIC Group: Institute for Cosmic Ray Research
(ICRR), The University of Tokyo, Kashiwa, 277-8582 Chiba, Japan  T. Inada
Japanese MAGIC Group: Institute for Cosmic Ray Research (ICRR), The University
of Tokyo, Kashiwa, 277-8582 Chiba, Japan  R. Iotov  Universität Würzburg,
D-97074 Würzburg, Germany  K. Ishio  University of Lodz, Faculty of Physics
and Applied Informatics, Department of Astrophysics, 90-236 Lodz, Poland  Y.
Iwamura  Japanese MAGIC Group: Institute for Cosmic Ray Research (ICRR), The
University of Tokyo, Kashiwa, 277-8582 Chiba, Japan  I. Jiménez Martínez
Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas,
E-28040 Madrid, Spain  J. Jormanainen  Finnish MAGIC Group: Finnish Centre for
Astronomy with ESO, University of Turku, FI-20014 Turku, Finland  L. Jouvin
Institut de Física d’Altes Energies (IFAE), The Barcelona Institute of Science
and Technology (BIST), E-08193 Bellaterra (Barcelona), Spain  D. Kerszberg
Institut de Física d’Altes Energies (IFAE), The Barcelona Institute of Science
and Technology (BIST), E-08193 Bellaterra (Barcelona), Spain  Y. Kobayashi
Japanese MAGIC Group: Institute for Cosmic Ray Research (ICRR), The University
of Tokyo, Kashiwa, 277-8582 Chiba, Japan  H. Kubo Japanese MAGIC Group:
Department of Physics, Kyoto University, 606-8502 Kyoto, Japan  J. Kushida
Japanese MAGIC Group: Department of Physics, Tokai University, Hiratsuka,
259-1292 Kanagawa, Japan  A. Lamastra National Institute for Astrophysics
(INAF), I-00136 Rome, Italy  D. Lelas Croatian MAGIC Group: University of
Split, Faculty of Electrical Engineering, Mechanical Engineering and Naval
Architecture (FESB), 21000 Split, Croatia  F. Leone National Institute for
Astrophysics (INAF), I-00136 Rome, Italy  E. Lindfors Finnish MAGIC Group:
Finnish Centre for Astronomy with ESO, University of Turku, FI-20014 Turku,
Finland  L. Linhoff Technische Universität Dortmund, D-44221 Dortmund, Germany
S. Lombardi National Institute for Astrophysics (INAF), I-00136 Rome, Italy
F. Longo Università di Udine and INFN Trieste, I-33100 Udine, Italy  also at
Dipartimento di Fisica, Università di Trieste, I-34127 Trieste, Italy  R.
López-Coto Università di Padova and INFN, I-35131 Padova, Italy  M. López-Moya
IPARCOS Institute and EMFTEL Department, Universidad Complutense de Madrid,
E-28040 Madrid, Spain  A. López-Oramas Instituto de Astrofísica de Canarias
and Dpto. de Astrofísica, Universidad de La Laguna, E-38200, La Laguna,
Tenerife, Spain  S. Loporchio INFN MAGIC Group: INFN Sezione di Bari and
Dipartimento Interateneo di Fisica dell’Università e del Politecnico di Bari,
I-70125 Bari, Italy  B. Machado de Oliveira Fraga Centro Brasileiro de
Pesquisas Físicas (CBPF), 22290-180 URCA, Rio de Janeiro (RJ), Brazil  C.
Maggio Departament de Física, and CERES-IEEC, Universitat Autònoma de
Barcelona, E-08193 Bellaterra, Spain  P. Majumdar Saha Institute of Nuclear
Physics, HBNI, 1/AF Bidhannagar, Salt Lake, Sector-1, Kolkata 700064, India
M. Makariev Inst. for Nucl. Research and Nucl. Energy, Bulgarian Academy of
Sciences, BG-1784 Sofia, Bulgaria  M. Mallamaci Università di Padova and INFN,
I-35131 Padova, Italy  G. Maneva Inst. for Nucl. Research and Nucl. Energy,
Bulgarian Academy of Sciences, BG-1784 Sofia, Bulgaria  M. Manganaro Croatian
MAGIC Group: University of Rijeka, Department of Physics, 51000 Rijeka,
Croatia  K. Mannheim Universität Würzburg, D-97074 Würzburg, Germany  M.
Mariotti Università di Padova and INFN, I-35131 Padova, Italy  M. Martínez
Institut de Física d’Altes Energies (IFAE), The Barcelona Institute of Science
and Technology (BIST), E-08193 Bellaterra (Barcelona), Spain  A. Mas Aguilar
IPARCOS Institute and EMFTEL Department, Universidad Complutense de Madrid,
E-28040 Madrid, Spain  D. Mazin Japanese MAGIC Group: Institute for Cosmic Ray
Research (ICRR), The University of Tokyo, Kashiwa, 277-8582 Chiba, Japan  Max-
Planck-Institut für Physik, D-80805 München, Germany  S. Menchiari  Università
di Siena and INFN Pisa, I-53100 Siena, Italy  S. Mender Technische Universität
Dortmund, D-44221 Dortmund, Germany  S. Mićanović Croatian MAGIC Group:
University of Rijeka, Department of Physics, 51000 Rijeka, Croatia  D. Miceli
Università di Udine and INFN Trieste, I-33100 Udine, Italy  now at Laboratoire
d’Annecy de Physique des Particules (LAPP), CNRS-IN2P3, 74941 Annecy Cedex,
France  T. Miener IPARCOS Institute and EMFTEL Department, Universidad
Complutense de Madrid, E-28040 Madrid, Spain  J. M. Miranda Università di
Siena and INFN Pisa, I-53100 Siena, Italy  R. Mirzoyan Max-Planck-Institut für
Physik, D-80805 München, Germany  E. Molina Universitat de Barcelona, ICCUB,
IEEC-UB, E-08028 Barcelona, Spain  A. Moralejo Institut de Física d’Altes
Energies (IFAE), The Barcelona Institute of Science and Technology (BIST),
E-08193 Bellaterra (Barcelona), Spain  D. Morcuende IPARCOS Institute and
EMFTEL Department, Universidad Complutense de Madrid, E-28040 Madrid, Spain
V. Moreno Departament de Física, and CERES-IEEC, Universitat Autònoma de
Barcelona, E-08193 Bellaterra, Spain  E. Moretti Institut de Física d’Altes
Energies (IFAE), The Barcelona Institute of Science and Technology (BIST),
E-08193 Bellaterra (Barcelona), Spain  T. Nakamori Japanese MAGIC Group:
Department of Physics, Yamagata University, Yamagata 990-8560, Japan  L. Nava
National Institute for Astrophysics (INAF), I-00136 Rome, Italy  V. Neustroev
Finnish MAGIC Group: Astronomy Research Unit, University of Oulu, FI-90014
Oulu, Finland  M. Nievas Rosillo Instituto de Astrofísica de Canarias and
Dpto. de Astrofísica, Universidad de La Laguna, E-38200, La Laguna, Tenerife,
Spain  C. Nigro Institut de Física d’Altes Energies (IFAE), The Barcelona
Institute of Science and Technology (BIST), E-08193 Bellaterra (Barcelona),
Spain  K. Nilsson Finnish MAGIC Group: Finnish Centre for Astronomy with ESO,
University of Turku, FI-20014 Turku, Finland  K. Nishijima Japanese MAGIC
Group: Department of Physics, Tokai University, Hiratsuka, 259-1292 Kanagawa,
Japan  K. Noda Japanese MAGIC Group: Institute for Cosmic Ray Research (ICRR),
The University of Tokyo, Kashiwa, 277-8582 Chiba, Japan  S. Nozaki Japanese
MAGIC Group: Department of Physics, Kyoto University, 606-8502 Kyoto, Japan
Y. Ohtani Japanese MAGIC Group: Institute for Cosmic Ray Research (ICRR), The
University of Tokyo, Kashiwa, 277-8582 Chiba, Japan  T. Oka Japanese MAGIC
Group: Department of Physics, Kyoto University, 606-8502 Kyoto, Japan  J.
Otero-Santos Instituto de Astrofísica de Canarias and Dpto. de Astrofísica,
Universidad de La Laguna, E-38200, La Laguna, Tenerife, Spain  S. Paiano
National Institute for Astrophysics (INAF), I-00136 Rome, Italy  M. Palatiello
Università di Udine and INFN Trieste, I-33100 Udine, Italy  D. Paneque Max-
Planck-Institut für Physik, D-80805 München, Germany  R. Paoletti Università
di Siena and INFN Pisa, I-53100 Siena, Italy  J. M. Paredes Universitat de
Barcelona, ICCUB, IEEC-UB, E-08028 Barcelona, Spain  L. Pavletić Croatian
MAGIC Group: University of Rijeka, Department of Physics, 51000 Rijeka,
Croatia  P. Peñil IPARCOS Institute and EMFTEL Department, Universidad
Complutense de Madrid, E-28040 Madrid, Spain  M. Persic Università di Udine
and INFN Trieste, I-33100 Udine, Italy  also at INAF Trieste and Dept. of
Physics and Astronomy, University of Bologna, Bologna, Italy  M. Pihet  Max-
Planck-Institut für Physik, D-80805 München, Germany  P. G. Prada Moroni
Università di Pisa and INFN Pisa, I-56126 Pisa, Italy  E. Prandini⋆ Università
di Padova and INFN, I-35131 Padova, Italy  C. Priyadarshi Institut de Física
d’Altes Energies (IFAE), The Barcelona Institute of Science and Technology
(BIST), E-08193 Bellaterra (Barcelona), Spain  I. Puljak Croatian MAGIC Group:
University of Split, Faculty of Electrical Engineering, Mechanical Engineering
and Naval Architecture (FESB), 21000 Split, Croatia  W. Rhode Technische
Universität Dortmund, D-44221 Dortmund, Germany  M. Ribó Universitat de
Barcelona, ICCUB, IEEC-UB, E-08028 Barcelona, Spain  J. Rico Institut de
Física d’Altes Energies (IFAE), The Barcelona Institute of Science and
Technology (BIST), E-08193 Bellaterra (Barcelona), Spain  C. Righi⋆ National
Institute for Astrophysics (INAF), I-00136 Rome, Italy  A. Rugliancich
Università di Pisa and INFN Pisa, I-56126 Pisa, Italy  N. Sahakyan⋆ Armenian
MAGIC Group: ICRANet-Armenia at NAS RA, 0019 Yerevan, Armenia  T. Saito
Japanese MAGIC Group: Institute for Cosmic Ray Research (ICRR), The University
of Tokyo, Kashiwa, 277-8582 Chiba, Japan  S. Sakurai Japanese MAGIC Group:
Institute for Cosmic Ray Research (ICRR), The University of Tokyo, Kashiwa,
277-8582 Chiba, Japan  K. Satalecka⋆ Deutsches Elektronen-Synchrotron (DESY),
D-15738 Zeuthen, Germany  F. G. Saturni National Institute for Astrophysics
(INAF), I-00136 Rome, Italy  B. Schleicher Universität Würzburg, D-97074
Würzburg, Germany  K. Schmidt Technische Universität Dortmund, D-44221
Dortmund, Germany  F. Schmuckermaier  Max-Planck-Institut für Physik, D-80805
München, Germany  T. Schweizer  Max-Planck-Institut für Physik, D-80805
München, Germany  J. Sitarek Japanese MAGIC Group: Institute for Cosmic Ray
Research (ICRR), The University of Tokyo, Kashiwa, 277-8582 Chiba, Japan  I.
Šnidarić  Croatian MAGIC Group: Ruder Bošković Institute, 10000 Zagreb,
Croatia  D. Sobczynska University of Lodz, Faculty of Physics and Applied
Informatics, Department of Astrophysics, 90-236 Lodz, Poland  A. Spolon
Università di Padova and INFN, I-35131 Padova, Italy  A. Stamerra National
Institute for Astrophysics (INAF), I-00136 Rome, Italy  J. Strišković Croatian
MAGIC Group: Josip Juraj Strossmayer University of Osijek, Department of
Physics, 31000 Osijek, Croatia  D. Strom Max-Planck-Institut für Physik,
D-80805 München, Germany  M. Strzys Japanese MAGIC Group: Institute for Cosmic
Ray Research (ICRR), The University of Tokyo, Kashiwa, 277-8582 Chiba, Japan
Y. Suda Japanese MAGIC Group: Physics Program, Graduate School of Advanced
Science and Engineering, Hiroshima University, 739-8526 Hiroshima, Japan  T.
Surić  Croatian MAGIC Group: Ruder Bošković Institute, 10000 Zagreb, Croatia
M. Takahashi Japanese MAGIC Group: Institute for Cosmic Ray Research (ICRR),
The University of Tokyo, Kashiwa, 277-8582 Chiba, Japan  R. Takeishi Japanese
MAGIC Group: Institute for Cosmic Ray Research (ICRR), The University of
Tokyo, Kashiwa, 277-8582 Chiba, Japan  F. Tavecchio National Institute for
Astrophysics (INAF), I-00136 Rome, Italy  P. Temnikov Inst. for Nucl. Research
and Nucl. Energy, Bulgarian Academy of Sciences, BG-1784 Sofia, Bulgaria  T.
Terzić Croatian MAGIC Group: University of Rijeka, Department of Physics,
51000 Rijeka, Croatia  M. Teshima  Max-Planck-Institut für Physik, D-80805
München, Germany  Japanese MAGIC Group: Institute for Cosmic Ray Research
(ICRR), The University of Tokyo, Kashiwa, 277-8582 Chiba, Japan  L. Tosti
INFN MAGIC Group: INFN Sezione di Perugia, I-06123 Perugia, Italy  S. Truzzi
Università di Siena and INFN Pisa, I-53100 Siena, Italy  A. Tutone National
Institute for Astrophysics (INAF), I-00136 Rome, Italy  S. Ubach  Departament
de Física, and CERES-IEEC, Universitat Autònoma de Barcelona, E-08193
Bellaterra, Spain  J. van Scherpenberg Max-Planck-Institut für Physik, D-80805
München, Germany  G. Vanzo Instituto de Astrofísica de Canarias and Dpto. de
Astrofísica, Universidad de La Laguna, E-38200, La Laguna, Tenerife, Spain  M.
Vazquez Acosta Instituto de Astrofísica de Canarias and Dpto. de Astrofísica,
Universidad de La Laguna, E-38200, La Laguna, Tenerife, Spain  S. Ventura
Università di Siena and INFN Pisa, I-53100 Siena, Italy  V. Verguilov Inst.
for Nucl. Research and Nucl. Energy, Bulgarian Academy of Sciences, BG-1784
Sofia, Bulgaria  I. Viale Università di Padova and INFN, I-35131 Padova, Italy
C. F. Vigorito INFN MAGIC Group: INFN Sezione di Torino and Università degli
Studi di Torino, I-10125 Torino, Italy  V. Vitale INFN MAGIC Group: INFN Roma
Tor Vergata, I-00133 Roma, Italy  I. Vovk Japanese MAGIC Group: Institute for
Cosmic Ray Research (ICRR), The University of Tokyo, Kashiwa, 277-8582 Chiba,
Japan  M. Will Max-Planck-Institut für Physik, D-80805 München, Germany  C.
Wunderlich Università di Siena and INFN Pisa, I-53100 Siena, Italy  T.
Yamamoto Japanese MAGIC Group: Department of Physics, Konan University, Kobe,
Hyogo 658-8501, Japan  D. Zarić Croatian MAGIC Group: University of Split,
Faculty of Electrical Engineering, Mechanical Engineering and Naval
Architecture (FESB), 21000 Split, Croatia  M. Hodges  Owens Valley Radio
Observatory, California Institute of Technology, Pasadena, CA 91125, USA  T.
Hovatta Finnish Center for Astronomy with ESO (FINCA), University of Turku,
FI-20014, Turku, Finland  Aalto University Metsähovi Radio Observatory,
Metsähovintie 114, 02540 Kylmälä, Finland  S. Kiehlmann Institute of
Astrophysics, Foundation for Research and Technology-Hellas, GR-71110
Heraklion, Greece  Department of Physics, Univ. of Crete, GR-70013 Heraklion,
Greece  I. Liodakis Finnish Center for Astronomy with ESO (FINCA), University
of Turku, FI-20014, Turku, Finland  W. Max-Moerbeck Departamento de
Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes,
Santiago, Chile  T. J. Pearson Owens Valley Radio Observatory, California
Institute of Technology, Pasadena, CA 91125, USA  A. C. S. Readhead  Owens
Valley Radio Observatory, California Institute of Technology, Pasadena, CA
91125, USA  R. A. Reeves Departamento de Astronomía, Universidad de
Conceptión, Concepción, Chile  A. Lähteenmäki Aalto University Metsähovi Radio
Observatory, Metsähovintie 114, 02540 Kylmälä, Finland  Aalto University
Department of Electronics and Nanoengineering, P.O. BOX 15500, FI-00076 AALTO,
Finland  M. Tornikoski Aalto University Metsähovi Radio Observatory,
Metsähovintie 114, 02540 Kylmälä, Finland  J. Tammi Aalto University Metsähovi
Radio Observatory, Metsähovintie 114, 02540 Kylmälä, Finland  F. D’Ammando⋆
INAF-Istituto di Radioastronomia, Via Gobetti 101, I-40129 Bologna, Italy  A.
Marchini Università di Siena, I-53100 Siena, Italy M. Cerruti, E. Prandini, C.
Righi,
K. Satalecka, N. Sahakyan, F. D’Ammando<EMAIL_ADDRESS>
###### Abstract

The blazar TXS 0506+056 (catalog ) 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ähovi, 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.

High energy astrophysics (739); Blazars (164); Jets (870); Gamma rays (637).

††facilities: MAGIC, $Fermi$/LAT, NuSTAR, Swift, KVA, REM, ASAS-SN, Metsähovi,
OVRO

##  1 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
$\gamma$ rays with both Fermi/LAT and MAGIC. Fermi/LAT detected the blazar in
a long-lasting (about six months long) flaring state, while MAGIC detected a
fast flare about 10 days after the IceCube neutrino. The chance probability
that the IceCube event and the flaring in the Fermi/LAT (MAGIC) energy band
happened simultaneously in space and time was estimated at the 3 $\sigma$ (3.5
$\sigma$) 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 $\sigma$ level (IceCube 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.

Figure 1: Multi-wavelength light curve of TXS 0506+056 during the 2017-2019
campaign. From top to bottom: MAGIC; $Fermi$/LAT; X-rays (Swift/XRT and
NuSTAR); UVOT; optical (ASAS-SN, KVA, REM, Siena); radio (OVRO and Metsähovi).
The light-blue dotted line shows the arrival time of the IceCube-170922A
neutrino, the grey dashed line marks the end of the data-set previously
published in (IceCube Collaboration et al., 2018a; 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.

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 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 $\gamma$-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-ray-neutral-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 $\gamma$-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 $\sigma$ in the neutrino sky (Aartsen et al.,
2020). The hottest spot in the all-sky scan (5 $\sigma$ pre-trials, reduced to
$1.3\sigma$ 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 $\sigma$ pre-trial; 2.9 $\sigma$ post-trials),
while TXS 0506+056 is significant at the 3.6 $\sigma$ 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
$Fermi$/LAT sample (Albert et al., 2020), significant at the 2.3 $\sigma$
level (combined space and time probability) and at the 2.6 $\sigma$ level when
considering information from IceCube, that also detected a high-energy
neutrino consistent with this source (IceCube-100608A).

Figure 2: MAGIC spectral energy distribution as measured during this campaign
(red). Low state is shown on the $left$ and the flare on the $right$. Overlaid
are the points from the low state ($left$, blue) and flares ($right$, blue and
green) observed in 2017 (Ansoldi et al., 2018).

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 spine-layer 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 $\gamma$-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) $\gamma$ 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,
$Fermi$/LAT, NuSTAR, Swift/XRT, Swift/UVOT, optical telescopes, Metsähovi 14
m, OVRO 40 m); in Section 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.

Table 1: MAGIC measurements of TXS 0506+056 Period | Exposure | Significance | Flux${>90\ \textrm{GeV}}$
---|---|---|---
| [h] |  | [10-11 cm-2 s-1]
MJD 58453 | 2.5 | 3.8 $\sigma$ | 9.8$\pm$2.5
MJD 58455 | 1.8 | 5.4 $\sigma$ | 18.0$\pm$3.4
MJD 58453-55 | 4.3 | 6.7 $\sigma$ | 12.5 $\pm$ 1.8
Rest | 74.4 | 4.0 $\sigma$ | 2.7$\pm$2.1

##  2 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 $\gamma$ 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.

###  2.1 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 $\sim$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
moonlight111Moonlight is defined as camera illumination resulting in the
measured current higher than 2000 $\mu$A) 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 $\sigma$, 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 $\pm$ 0.2 for the flare and 3.5 $\pm$ 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
$\pm 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 ($\sim$74
h). The average flux F($>$90 GeV) = (2.7$\pm$2.1)$\times$10-11 cm-2 s-1,
calculated from the excess at a significance of 4 $\sigma$, is the lowest VHE
$\gamma$-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.

###  2.2 Fermi/LAT

In this work, the publicly available $Fermi$/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^{\circ}$
around the $\gamma$-ray position of TXS 0506+056 (RA, Dec)= (77.359, 5.701)
were downloaded and analyzed using the standard procedure described by the
$Fermi$/LAT Collaboration. The data are binned into pixels of
$0.1^{\circ}\times 0.1^{\circ}$ and standard cuts were applied to select the
good time intervals (”DATA_QUAL$>$0” and ”LAT_CONFIG==1”, $z_{\rm
max}<90^{\circ}$). The model file was created based on the $Fermi$/LAT fourth
source catalog (4FGL, Abdollahi et al., 2020) and it includes all point
sources within $17^{\circ}$ 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 $\gamma$-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^{\circ}$ 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_{\rm 0}$). The adaptively binned light curve
computed for 20% uncertainty and $E_{\rm 0}$ = 216.03 MeV is shown in Fig. 1.
The $\gamma$-ray flux and photon index variations are investigated using the
7-day binned light curve (Fig. 5).

In the low state, the $\gamma$-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 $\gamma$-ray
emission spectra in the low and active states are shown in Fig. 6.

###  2.3 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
produced 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 $\chi^{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 nH = 1.11 $\times$ 1021 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 $\pm$
0.34) $\times$ 10-12 to (7.15 $\pm$ 0.39) $\times$ 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.

###  2.4 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 Centre222http://www.swift.ac.uk/user_objects (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 $\times$ 1021 cm-2. A large number of spectra show low
number of counts (i.e. $<$ 200), therefore not allowing us to use the
$\chi^{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 spectral 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,
$\Gamma_{\rm\,X}$, varying between 1.27 $\pm$ 0.37 and 2.71 $\pm$ 0.70.

###  2.5 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.333see https://www.swift.ac.uk/analysis/uvot/index.php 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.444https://asas-sn.osu.edu/ 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 package555https://diffractionlimited.com/product/maxim-dl/,
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$.

###  2.7 Metsähovi

The 37 GHz observations were performed with the 13.7 m diameter Metsähovi
radio telescope.666https://www.aalto.fi/en/services/metsahovis-main-
instruments A typical integration time to obtain one flux density data point
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.

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

Figure 3: Fractional variability parameter for each instrument. All points are
calculated using daily averaged fluxes, apart from the point at 1 GeV which
uses a 7-days average of $Fermi$/LAT flux.

##  3 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 Fvar as defined by Vaughan et al. (2003) to
quantify the observed trends. The different sensitivity and flux sampling of
the different 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 Fvar values for all the instruments contributing
to this study, apart from MAGIC. The very low VHE $\gamma$-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 $Fermi$/LAT and ASAS-SN. A relatively large
data set was also collected with the X-ray and UV instruments on board $Swift$
and the KVA optical telescope. High sensitivity NuSTAR observations complete
the multi-wavelength picture. The most pronounced variability
F${}_{var}\sim$0.40 is observed in the $Fermi$/LAT $\gamma$-ray band, while
the X-ray variability is at a lower level of $\sim$0.25-0.30. The radio,
optical and UV bands display a moderate variability of
F${}_{var}\sim$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 $\gamma$-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.

Figure 4: Left: Correlation between the optical flux in the V band and the
photon index in the optical-to-UV band measured by UVOT. Right: Correlation
between the X-ray flux in the 0.3-10, and 2-10 keV energy bands and the photon
indexes measured with Swift-XRT and that in the 3-78 keV band as measured with
NuSTAR. A clear break emerges between the two energetic regimes. In both
panels, the dashed lines show the average trend of the data and were obtained
with a linear fit.

###  3.1 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 $\gamma$-ray blazars, a harder-when-brighter
behaviour, meaning a harder spectrum during enhanced flux states, is observed
in both X- and $\gamma$-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).

Figure 5: Correlation between the photon flux in the $\gamma$-ray band and the
photon index measured by Fermi/LAT.

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.

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\times 10^{-15}$
X-ray (0.3–10 keV) | -0.23 | 0.14
X-ray (3–78 keV) | -0.25 | 0.59
$\gamma$-ray (0.2–300 GeV) | 0.36 | 0.006

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\times 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 softer-when-
brighter trend in $\gamma$ 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.

##  4 Multi-messenger modeling

The SEDs of TXS 0506+056 during the 2017-2019 multi-wavelength campaign are
shown in Fig. 6, separately for the VHE flare of December, 2018, and the
average low-state. For the latter, we use the low-state spectra from Fermi/LAT
and MAGIC, and NuSTAR and Swift data from October 16, 2018 as an exemplary
spectrum for the low-state (together with, in gray, the whole range of
optical-to-X-ray observations). For the VHE high-state, the only instrument
performing contemporaneous observations with MAGIC were $Fermi$/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
$\gamma$-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 $\gamma$ 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-by-generation.
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
$\Gamma$ (translated into the Doppler factor $\delta$ assuming an angle to the
line of sight $\theta_{\rm view}$), the radius $R^{\prime}$, and the magnetic
flux density $B^{\prime}$.777Here 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. 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 $\gamma^{\prime}_{\rm{Max}}$, and are thus
characterized by six free parameters each (the three Lorentz factors
$\gamma^{\prime}_{\rm{min}}$, $\gamma^{\prime}_{\rm{break}}$,
$\gamma^{\prime}_{\rm{Max}}$, the two indexes $\alpha_{1}$, and $\alpha_{2}$,
and the normalization $K^{\prime}$), 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^{\prime\prime}_{\rm in,layer}$ and
$R^{\prime\prime}_{\rm out,layer}$, and its height $H^{\prime\prime}_{\rm
layer}$), moving with bulk Lorentz factor $\Gamma_{\rm layer}$ and filled with
a homogeneous magnetic field $B^{\prime\prime}_{\rm 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 $\alpha_{1,\rm{e}}=\alpha_{1,\rm{p}}$; the primary protons are not
cooled (i.e. all cooling processes become relevant at Lorentz factors higher
than $\gamma^{\prime}_{\rm{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
$\theta_{\rm view}=0.8^{\circ}$; the geometry of the layer is
$R^{\prime\prime}_{\rm in,layer}=R^{\prime}$, $R^{\prime\prime}_{\rm
out,layer}=1.5R^{\prime}$, while the height of the layer in the frame of the
spine is $H^{\prime}_{\rm layer}=H^{\prime\prime}_{\rm
layer}/(\Gamma\Gamma_{\rm layer}(1-\beta\beta_{\rm layer}))=R^{\prime}$ (see
Ghisellini et al., 2005); the photon field from the layer is computed using
$\delta_{\rm layer}=4$, $B^{\prime}_{\rm layer}=1.8\ \textrm{G}$, $\alpha_{\rm
e,1,layer}=2.0$, $\alpha_{\rm e,2,layer}=2.9$, $\gamma^{\prime}_{\rm
e,min,layer}=290$, $\gamma^{\prime}_{\rm e,break,layer}=2000$, and
$\gamma^{\prime}_{\rm e,Max,layer}=3.9\times 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^{\prime}_{\rm ph,layer}$. The
external inverse-Compton component is corrected by a factor
$\delta/\delta_{\rm layer}$ following Dermer (1995). These physical
constraints reduce the number of free parameters of the model to eleven:
$R^{\prime}$, $B^{\prime}$, $\gamma^{\prime}_{\rm{e,min}}$,
$\gamma^{\prime}_{\rm{e,break}}$, $\gamma^{\prime}_{\rm{e,Max}}$,
$\gamma^{\prime}_{\rm{p,Max}}$, $\alpha_{1}$, $\alpha_{2\rm{,e}}$,
$K_{\rm{e}}$, $K_{\rm{p}}$, and $u^{\prime}_{\rm 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 ourselves to identify
a lepto-hadronic solution (among many) that can reproduce the observations.

Figure 6: Spectral energy distribution of TXS 0506+056 during the 2017-2019
campaign, modeled in a lepto-hadronic scenario. Left: low-state (average)
during the multi-wavelength campaign. Black data points show the $Fermi$/LAT
average spectrum and MAGIC upper limits, and Swift and NuSTAR spectra on
October 16, 2018. Gray data points show the whole range of optical-to-X-ray
spectra during the campaign. Right: high-state in VHE $\gamma$ rays during
December 2018. Black data points show the $Fermi$/LAT, MAGIC, and ASAS-SN
contemporaneous spectra. Gray data points show the nearest observations in
radio, optical, UV and X-rays on December 8, 2018. In both plots, the solid
lines represent synchrotron emission by primary electrons (red); inverse-
Compton emission over photons from the jet layer (green); Bethe-Heitler
cascade (orange); cascade from $\pi^{0}$ (violet) and $\pi^{\pm}$ (pink)
decay; neutrino emission (blue). In the left plot we show the IceCube
$5\sigma$ sensitivity curve for point-source searches (for eight years
exposure and $0^{\circ}$ declination, see Aartsen et al., 2019). In the right
plot we show the IceCube and ANTARES flux upper limits (at lower and higher
energies, respectively) following the detection of a VHE flare in December
2018.
Table 3: Parameters of the lepto-hadronic models

| Low-state | High-state
---|---|---
$z$ | $0.337$ | ′′
$\theta_{\rm view}$ [deg] | $0.8$ | ′′
$\delta$ | $40$ | ′′
${}^{\star}\Gamma$ | $22$ | ′′
$R^{\prime}$ [1016 cm] | $1.1$ | ′′
${}^{\star}\tau_{\rm var}$ [hours] | $3.4$ | ′′
$B^{\prime}$ [G] | $1$ | ′′
${}^{\star}u^{\prime}_{B}$ [erg cm-3] | $0.04$ | ′′
$\gamma_{\rm e,min}^{\prime}$ | $800$ | ′′
$\gamma_{\rm e,break}^{\prime}$ | $2.2\times 10^{3}$ | $1\times 10^{4}$
$\gamma_{\rm e,max}^{\prime}$ | $10^{4}$ | $2\times 10^{4}$
$\alpha_{\rm e,1}=\alpha_{p,1}$ | $2.0$ | ′′
$\alpha_{\rm e,2}$ | $4.0$ | ′′
$K^{\prime}_{\rm e}$ [cm-3] | $10^{3}$ | $520$
${}^{\star}u^{\prime}_{\rm e}$ [erg cm-3] | $1.0\times 10^{-3}$ | $1.0\times 10^{-3}$
$\gamma_{\rm p,min}^{\prime}$ | $1$ | ′′
$\gamma_{\rm p,max}^{\prime}[10^{7}]$ | $1$ | ′′
$K_{\rm p}^{\prime}$ [cm-3] | $3.0\times 10^{3}$ | $5.0\times 10^{3}$
${}^{\star}u^{\prime}_{\rm p}$ [erg cm-3] | $70$ | $120$
$\delta_{{\rm layer}}$ | $4$ | ′′
${}^{\star}\Gamma_{{\rm layer}}$ | $2.2$ | ′′
$R^{\prime}_{{\rm in,layer}}$ [1016 cm] | $1.1$ | ′′
$B_{\rm layer}^{\prime\prime}$ [G] | $1.8$ | ′′
$\gamma_{\rm e,min,layer}^{\prime\prime}$ | $290$ | ′′
$\gamma_{\rm e,break,layer}^{\prime\prime}$ | $2000$ | ′′
$\gamma_{\rm e,max,layer}^{\prime\prime}$ | $3.9\times 10^{5}$ | ′′
$\alpha_{\rm e,1,layer}$ | $2.0$ | ′′
$\alpha_{\rm e,2,layer}$ | $2.9$ | ′′
$u^{\prime}_{\rm ph,layer}$ [erg cm-3] | $1.5\times 10^{-2}$ | ′′
${}^{\star}u^{\prime}_{\rm e}/u^{\prime}_{\rm B}[10^{3}]$ | $0.025$ | $0.024$
${}^{\star}u^{\prime}_{\rm p}/u^{\prime}_{\rm B}[10^{3}]$ | $1.7\times 10^{3}$ | $2.9\times 10^{3}$
${}^{\star}L$ [1047 erg s-1] | $7.6$ | $12.7$
${}^{\star}\nu$ [yr-1] | $0.14$ | -
${}^{\star}\nu_{\rm{GOLD}}$ [day-1] | - | $2.1\times 10^{-4}$

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\pi R^{\prime 2}c\Gamma^{2}(u^{\prime}_{\rm B}+u^{\prime}_{\rm
e}+u^{\prime}_{\rm p})$, where $u^{\prime}_{\rm B}$, $u^{\prime}_{\rm e}$, and
$u^{\prime}_{\rm 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.

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 $\gamma$-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 IceCube 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 $\gamma_{\rm p,Max}^{\prime}$) are made. The
maximum proton energy is related to the peak energy of the neutrino emission:
lower values of $\gamma_{\rm p,Max}^{\prime}$ 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\times 10^{47}$ erg s-1, or
$7.6-76$ $L_{Edd}$ for a super-massive black-hole with mass
$M_{\bullet}=10^{9-8}M_{\odot}$ ) 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 $\gamma_{\rm p,min}^{\prime}$. The model presented here,
with $\gamma_{\rm p,min}^{\prime}=1$ and $\alpha_{\rm{p}}=2.0$ is explicitly
conservative in this sense: a harder particle index, or $\gamma_{\rm
p,min}^{\prime}\gg 1$ will significantly lower the total jet power.

Once we obtained a satisfactory modeling of the low-state 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 multi-wavelength 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
$\gamma$-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\times 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 $\gamma$-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\times 10^{-4}$
events in a day, or $1.5\times 10^{-3}$ in a week, completely consistent with
the non-detection of neutrinos associated with this $\gamma$-ray flare.

##  5 Summary and conclusions

We have performed an extensive multi-wavelength campaign from radio to VHE
$\gamma$ 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 time-scale of a few
years.

In the VHE band, the source showed another day-long 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 $\gamma$-ray excess at the level of 4 $\sigma$. No
neutrino flares associated with the VHE flare in December 2018 have been
reported. In the HE $\gamma$-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 $\gamma$-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 $Fermi$/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 X-rays: 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-when-brighter 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 $\gamma$-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 $\gamma$-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\times 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
$\gamma$-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 $\gamma$-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
multi-messenger 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.

We would like to thank the Instituto de Astrofísica de Canarias for the
excellent working conditions at the Observatorio del Roque de los Muchachos in
La Palma. The financial support of the German BMBF, MPG and HGF; the Italian
INFN and INAF; the Swiss National Fund SNF; the ERDF under the Spanish
Ministerio de Ciencia e Innovación (MICINN) (PID2019-104114RB-C31,
PID2019-104114RB-C32, PID2019-104114RB-C33,
PID2019-105510GB-C31,PID2019-107847RB-C41, PID2019-107847RB-C42,
PID2019-107847RB-C44, PID2019-107988GB-C22); the Indian Department of Atomic
Energy; the Japanese ICRR, the University of Tokyo, JSPS, and MEXT; the
Bulgarian Ministry of Education and Science, National RI Roadmap Project
DO1-400/18.12.2020 and the Academy of Finland grant nr. 320045 is gratefully
acknowledged. This work was also supported by the Spanish Centro de Excelencia
“Severo Ochoa” (SEV-2016-0588, SEV-2017-0709, CEX2019-000920-S), the Unidad de
Excelencia “María de Maeztu” (CEX2019-000918-M, MDM-2015-0509-18-2) and by the
CERCA program of the Generalitat de Catalunya; 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. E.P. acknowledges funding from
Italian Ministry of Education, University and Research (MIUR) through the
”Dipartimenti di eccellenza” project Science of the Universe. M.C. has
received financial support through the Postdoctoral Junior Leader Fellowship
Programme from la Caixa Banking Foundation, grant No. LCF/BQ/LI18/11630012.
This research has made use of data from the OVRO 40-m monitoring program which
was supported in part by NASA grants NNX08AW31G, NNX11A043G and NNX14AQ89G,
and NSF grants AST-0808050 and AST-1109911, and private funding from Caltech
and the MPIfR. This publication makes use of data obtained at Metsähovi Radio
Observatory, operated by Aalto University in Finland.

##  Appendix A Multiwavelength data tables

Table 4: Log and fitting results of NuSTAR observations of TXS 0506+056 using a power-law model model with $N_{\rm H}$ fixed to Galactic absorption, i.e. 1.11$\times$1021 cm-2. Date | MJD | Net exposure time | Flux${}_{\rm\,3.0-78\,keV}$ | Photon index
---|---|---|---|---
[UT] |  | [s] | [10-12 erg cm-2 s-1] | [$\Gamma_{\rm\,X}$]
2018-04-03 | 58211 | 26030 | 5.82 ${}^{+0.33}_{-0.39}$ | 1.57 $\pm$ 0.09
2018-10-16 | 58407 | 31868 | 6.21 ${}^{+0.39}_{-0.35}$ | 1.52 $\pm$ 0.08
2018-11-15 | 58437 | 26880 | 4.24 ${}^{+0.26}_{-0.28}$ | 1.64 $\pm$ 0.11
2018-12-08 | 58460 | 25899 | 4.54 ${}^{+0.34}_{-0.34}$ | 1.68 $\pm$ 0.10
2019-01-07 | 58490 | 25889 | 7.15 ${}^{+0.38}_{-0.40}$ | 1.73 $\pm$ 0.08
Table 5: Fractional variability of TXS 0506+056 Instrument | Fermi/LAT | NuSTAR | Swift/XRT | Swift/UVOT | KVA | ASAS-SN | OVRO
---|---|---|---|---|---|---|---
Fvar |  0.410 $\pm$ 0.045 |  0.25$\pm$0.02 |  0.30$\pm$0.06 |  0.128$\pm$0.007 (V) |  0.194$\pm$0.002 |  0.143$\pm$0.001 (g) |  0.216 $\pm$ 0.002
|  |  |  |  0.120$\pm$0.005 (B) |  |  0.136$\pm$0.004 (V) |
|  |  |  |  0.138$\pm$0.005 (U) |  |  |
|  |  |  |  0.157$\pm$0.005 (W1) |  |  |
|  |  |  |  0.166$\pm$0.005 (M2) |  |  |
|  |  |  |  0.156$\pm$0.004 (W2) |  |  |

## References

  * Aartsen et al. (2017) Aartsen, M. G., Abraham, K., Ackermann, M., et al. 2017, ApJ, 835, 45, doi: 10.3847/1538-4357/835/1/45
  * Aartsen et al. (2019) Aartsen, M. G., Ackermann, M., Adams, J., et al. 2019, European Physical Journal C, 79, 234, doi: 10.1140/epjc/s10052-019-6680-0
  * Aartsen et al. (2020) —. 2020, Phys. Rev. Lett., 124, 051103, doi: 10.1103/PhysRevLett.124.051103
  * Abdo et al. (2010) Abdo, A. A., Ackermann, M., Agudo, I., et al. 2010, ApJ, 716, 30, doi: 10.1088/0004-637X/716/1/30
  * Abdollahi et al. (2020) Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33, doi: 10.3847/1538-4365/ab6bcb
  * Abeysekara et al. (2018) Abeysekara, A. U., Archer, A., Benbow, W., et al. 2018, ApJ, 861, L20, doi: 10.3847/2041-8213/aad053
  * Ageron et al. (2011) Ageron, M., Aguilar, J. A., Al Samarai, I., et al. 2011, Nuclear Instruments and Methods in Physics Research A, 656, 11, doi: 10.1016/j.nima.2011.06.103
  * Albert et al. (2020) Albert, A., André, M., Anghinolfi, M., et al. 2020, arXiv e-prints, arXiv:2012.15082.  https://arxiv.org/abs/2012.15082
  * Aleksić et al. (2014) Aleksić, J., Ansoldi, S., Antonelli, L. A., et al. 2014, A&A, 572, A121, doi: 10.1051/0004-6361/201424254
  * Aleksić et al. (2016) Aleksić, J., et al. 2016, Astropart. Phys., 72, 76, doi: 10.1016/j.astropartphys.2015.02.005
  * Ansoldi et al. (2018) Ansoldi, S., Antonelli, L. A., Arcaro, C., et al. 2018, ApJ, 863, L10, doi: 10.3847/2041-8213/aad083
  * Arnaud (1996) Arnaud, K. A. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes, 17
  * Atwood et al. (2009) Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071, doi: 10.1088/0004-637X/697/2/1071
  * Banik et al. (2020) Banik, P., Bhadra, A., Pandey, M., & Majumdar, D. 2020, Phys. Rev. D, 101, 063024, doi: 10.1103/PhysRevD.101.063024
  * Blaufuss et al. (2019) Blaufuss, E., Kintscher, T., Lu, L., & Tung, C. F. 2019, in International Cosmic Ray Conference, Vol. 36, 36th International Cosmic Ray Conference (ICRC2019), 1021.  https://arxiv.org/abs/1908.04884
  * Burrows et al. (2005) Burrows, D. N., Hill, J. E., Nousek, J. A., et al. 2005, Space Sci. Rev., 120, 165, doi: 10.1007/s11214-005-5097-2
  * Cash (1979) Cash, W. 1979, ApJ, 228, 939, doi: 10.1086/156922
  * Cerruti et al. (2019) Cerruti, M., Zech, A., Boisson, C., et al. 2019, MNRAS, 483, L12, doi: 10.1093/mnrasl/sly210
  * Cerruti et al. (2015) Cerruti, M., Zech, A., Boisson, C., & Inoue, S. 2015, MNRAS, 448, 910, doi: 10.1093/mnras/stu2691
  * Coleiro & Dornic (2018) Coleiro, A., & Dornic, D. 2018, The Astronomer’s Telegram, 12274, 1
  * Covino et al. (2004) Covino, S., Stefanon, M., Sciuto, G., et al. 2004, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 5492, Ground-based Instrumentation for Astronomy, ed. A. F. M. Moorwood & M. Iye, 1613–1622, doi: 10.1117/12.551532
  * Dermer (1995) Dermer, C. D. 1995, ApJ, 446, L63, doi: 10.1086/187931
  * Evans et al. (2009) Evans, P. A., Beardmore, A. P., Page, K. L., et al. 2009, MNRAS, 397, 1177, doi: 10.1111/j.1365-2966.2009.14913.x
  * Fomin et al. (1994) Fomin, V. P., Stepanian, A. A., Lamb, R. C., et al. 1994, Astroparticle Physics, 2, 137, doi: 10.1016/0927-6505(94)90036-1
  * Franceschini et al. (2008) Franceschini, A., Rodighiero, G., & Vaccari, M. 2008, A&A, 487, 837, doi: 10.1051/0004-6361:200809691
  * Fruck et al. (2014) Fruck, C., Gaug, M., Zanin, R., et al. 2014, arXiv e-prints, arXiv:1403.3591.  https://arxiv.org/abs/1403.3591
  * Gao et al. (2019) Gao, S., Fedynitch, A., Winter, W., & Pohl, M. 2019, Nature Astronomy, 3, 88, doi: 10.1038/s41550-018-0610-1
  * Ghisellini et al. (2005) Ghisellini, G., Tavecchio, F., & Chiaberge, M. 2005, A&A, 432, 401, doi: 10.1051/0004-6361:20041404
  * Harrison et al. (2013) Harrison, F. A., Craig, W. W., Christensen, F. E., et al. 2013, ApJ, 770, 103, doi: 10.1088/0004-637X/770/2/103
  * Hooper et al. (2019) Hooper, D., Linden, T., & Vieregg, A. 2019, J. Cosmology Astropart. Phys, 2019, 012, doi: 10.1088/1475-7516/2019/02/012
  * Hovatta et al. (2021) Hovatta, T., Lindfors, E., Kiehlmann, S., et al. 2021, A&A, 650, A83, doi: 10.1051/0004-6361/202039481
  * Hwang et al. (2020) Hwang, S., Im, M., Taak, Y. C., et al. 2020, arXiv e-prints, arXiv:2011.14049.  https://arxiv.org/abs/2011.14049
  * IceCube Collaboration et al. (2006) IceCube Collaboration, Achterberg, A., Ackermann, M., et al. 2006, Astroparticle Physics, 26, 155, doi: 10.1016/j.astropartphys.2006.06.007
  * IceCube Collaboration et al. (2018a) IceCube Collaboration, Aartsen, M. G., Ackermann, M., et al. 2018a, Science, 361, eaat1378, doi: 10.1126/science.aat1378
  * IceCube Collaboration et al. (2018b) —. 2018b, Science, 361, 147, doi: 10.1126/science.aat2890
  * Jorstad et al. (2001) Jorstad, S. G., Marscher, A. P., Mattox, J. R., et al. 2001, ApJ, 556, 738, doi: 10.1086/321605
  * Kalberla et al. (2005) Kalberla, P. M. W., et al. 2005, A&A, 440, 775, doi: 10.1051/0004-6361:20041864
  * Keivani et al. (2018) Keivani, A., Murase, K., Petropoulou, M., et al. 2018, ApJ, 864, 84, doi: 10.3847/1538-4357/aad59a
  * Kelner & Aharonian (2008) Kelner, S. R., & Aharonian, F. A. 2008, Phys. Rev. D, 78, 034013, doi: 10.1103/PhysRevD.78.034013
  * Kochanek et al. (2017) Kochanek, C. S., Shappee, B. J., Stanek, K. Z., et al. 2017, PASP, 129, 104502, doi: 10.1088/1538-3873/aa80d9
  * Lähteenmäki & Valtaoja (2003) Lähteenmäki, A., & Valtaoja, E. 2003, ApJ, 590, 95, doi: 10.1086/374883
  * León-Tavares et al. (2011) León-Tavares, J., Valtaoja, E., Tornikoski, M., Lähteenmäki, A., & Nieppola, E. 2011, A&A, 532, A146, doi: 10.1051/0004-6361/201116664
  * Liodakis et al. (2018) Liodakis, I., Romani, R. W., Filippenko, A. V., et al. 2018, Monthly Notices of the Royal Astronomical Society, 480, 5517, doi: 10.1093/mnras/sty2264
  * Lott et al. (2012) Lott, B., et al. 2012, A&A, 544, A6, doi: 10.1051/0004-6361/201218873
  * MAGIC Collaboration et al. (2020) MAGIC Collaboration, Acciari, V. A., Ansoldi, S., et al. 2020, A&A, 638, A14, doi: 10.1051/0004-6361/201935450
  * Mirzoyan (2018) Mirzoyan, R. 2018, The Astronomer’s Telegram, 12260, 1
  * Moralejo et al. (2009) Moralejo, A., Gaug, M., Carmona, E., et al. 2009, arXiv e-prints, arXiv:0907.0943.  https://arxiv.org/abs/0907.0943
  * Mücke et al. (2000) Mücke, A., Engel, R., Rachen, J. P., Protheroe, R. J., & Stanev, T. 2000, Computer Physics Communications, 124, 290, doi: 10.1016/S0010-4655(99)00446-4
  * Murase et al. (2018) Murase, K., Oikonomou, F., & Petropoulou, M. 2018, ApJ, 865, 124, doi: 10.3847/1538-4357/aada00
  * Nilsson et al. (2018) Nilsson, K., , et al. 2018, A&A, 620, A185, doi: 10.1051/0004-6361/201833621
  * Padovani et al. (2018) Padovani, P., Giommi, P., Resconi, E., et al. 2018, MNRAS, 480, 192, doi: 10.1093/mnras/sty1852
  * Padovani et al. (2019) Padovani, P., Oikonomou, F., Petropoulou, M., Giommi, P., & Resconi, E. 2019, MNRAS, 484, L104, doi: 10.1093/mnrasl/slz011
  * Paiano et al. (2018) Paiano, S., Falomo, R., Treves, A., & Scarpa, R. 2018, ApJ, 854, L32, doi: 10.3847/2041-8213/aaad5e
  * Petropoulou et al. (2020) Petropoulou, M., Murase, K., Santander, M., et al. 2020, ApJ, 891, 115, doi: 10.3847/1538-4357/ab76d0
  * Reimer et al. (2019) Reimer, A., Böttcher, M., & Buson, S. 2019, ApJ, 881, 46, doi: 10.3847/1538-4357/ab2bff
  * Richards et al. (2011) Richards, J. L., Max-Moerbeck, W., Pavlidou, V., et al. 2011, ApJS, 194, 29, doi: 10.1088/0067-0049/194/2/29
  * Rodrigues et al. (2019) Rodrigues, X., Gao, S., Fedynitch, A., Palladino, A., & Winter, W. 2019, ApJ, 874, L29, doi: 10.3847/2041-8213/ab1267
  * Rolke et al. (2005) Rolke, W. A., López, A. M., & Conrad, J. 2005, Nuclear Instruments and Methods in Physics Research A, 551, 493, doi: 10.1016/j.nima.2005.05.068
  * Roming et al. (2005) Roming, P. W. A., Kennedy, T. E., Mason, K. O., et al. 2005, Space Sci. Rev., 120, 95, doi: 10.1007/s11214-005-5095-4
  * Ros et al. (2020) Ros, E., Kadler, M., Perucho, M., et al. 2020, Astronomy & Astrophysics, 633, L1, doi: 10.1051/0004-6361/201937206
  * Sahakyan (2018) Sahakyan, N. 2018, ApJ, 866, 109, doi: 10.3847/1538-4357/aadade
  * Schlafly & Finkbeiner (2011) Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103, doi: 10.1088/0004-637X/737/2/103
  * Tavecchio et al. (2010) Tavecchio, F., Ghisellini, G., Ghirlanda, G., Foschini, L., & Maraschi, L. 2010, MNRAS, 401, 1570, doi: 10.1111/j.1365-2966.2009.15784.x
  * Teräsranta et al. (1998) Teräsranta, H., Tornikoski, M., Mujunen, A., et al. 1998, Astron. Astrophys. Suppl. Ser., 132, 305, doi: 10.1051/aas:1998297
  * Vandenbroucke (2018) Vandenbroucke, J. 2018, The Astronomer’s Telegram, 12267, 1
  * Vaughan et al. (2003) Vaughan, S., Edelson, R., Warwick, R. S., & Uttley, P. 2003, MNRAS, 345, 1271, doi: 10.1046/j.1365-2966.2003.07042.x
  * Wilms et al. (2000) Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914, doi: 10.1086/317016
  * Xue et al. (2019) Xue, R., Liu, R.-Y., Petropoulou, M., et al. 2019, ApJ, 886, 23, doi: 10.3847/1538-4357/ab4b44
  * Zanin et al. (2013) Zanin, R., Carmona, E., Sitarek, J., et al. 2013, in International Cosmic Ray Conference, Vol. 33, International Cosmic Ray Conference, 2937
  * Zerbi et al. (2001) Zerbi, R. M., Chincarini, G., Ghisellini, G., et al. 2001, Astronomische Nachrichten, 322, 275, doi: 10.1002/1521-3994(200112)322:5/6<275::AID-ASNA275>3.0.CO;2-N
  * Zhang et al. (2020) Zhang, B. T., Petropoulou, M., Murase, K., & Oikonomou, F. 2020, ApJ, 889, 118, doi: 10.3847/1538-4357/ab659a
  * Zhang et al. (2012) Zhang, J., Liang, E.-W., Zhang, S.-N., & Bai, J. M. 2012, ApJ, 752, 157, doi: 10.1088/0004-637X/752/2/157