Multimessenger Constraints on Intergalactic Magnetic Fields from the Flare of TXS 0506+056

The origin of magnetic fields in the universe is an open problem. Seed magnetic fields possibly produced in early times may have survived up to the present day close to their original form, providing an untapped window to the primeval universe. The recent observations of high-energy neutrinos from the blazar TXS 0506+056 in association with an electromagnetic counterpart in a broad range of wavelengths can be used to probe intergalactic magnetic fields via the time delay between the neutrinos and gamma rays as well as the time dependence of the gamma-ray fluxes. Using extensive three-dimensional Monte Carlo simulations, we present a novel method to constrain these fields. We apply it to TXS 0506+056 and, for the first time, derive constraints on both the magnetic-field strength and its coherence length, considering six orders of magnitude for each.


INTRODUCTION
A long-standing problem in cosmology concerns the origin of magnetic fields in the universe. Two broad classes of mechanisms to explain magnetogenesis exist. Primordial (or cosmological) mechanisms posit that global processes taking place in the early universe could give rise to seed magnetic fields. Examples of such processes are inflation and phase transitions such as the electroweak and the quantum chromodynamics phase transitions (see Durrer & Neronov 2013 for a review). Astrophysical mechanisms, on the other hand, suggest that small-scale processes during the formation of structures gave rise to magnetic fields.
A way to distinguish between primordial and astrophysical mechanisms is to look for signs of magnetisation in cosmic voids. A negative signal would favour an astrophysical origin, but a positive one would not necessarily provide unambiguous evidence of a cosmological origin, since winds/outflows could carry magnetised material to the intergalactic medium (IGM). Nevertheless, r.batista@astro.ru.nl, andrey.saveliev@desy.de the contamination of the IGM by winds and outflows is limited, and the seed magnetic fields far from structures, near the centre of voids, should remain in their pristine form (Furlanetto & Loeb 2001;Bertone et al. 2006). Moreover, primordial and astrophysical mechanisms lead to distinct magnetic power spectra and hence different coherence lengths.
Faraday rotation measures of distant objects (Vallee 1991) and the cosmic microwave background (CMB) (Jedamzik & Saveliev 2019) provide the respective upper limits of B 10 −9 G and B 10 −11 G on the strength of intergalactic magnetic fields (IGMFs) at large scales. Gamma-ray-induced electromagnetic cascades yield, in general, the lower limit B 10 −16.5 G (Neronov & Vovk 2010;Tavecchio et al. 2010;Taylor et al. 2011;Vovk et al. 2012;Ackermann et al. 2018), with some additional exclusion regions between 10 −15 and 10 −13 G from the absence of extended emission around objects (H. E. S. S. Collaboration 2014; VERITAS Collaboration 2017). The coherence length of IGMFs (L c ) is poorly constrained, lying in the range 10 −12 L c /Mpc 10 3.5 , the lower bound corresponding to the resistivity time scale of the universe, and the upper limit referring to the Hubble horizon.
The detection of the high-energy neutrino event IC-170922A by the IceCube Observatory (IceCube Collaboration 2018; IceCube et al. 2018) together with electromagnetic counterparts at various wavelengths (Ice-Cube et al. 2018) marked the dawn of neutrino astronomy. These observations have been associated with a flare of the blazar TXS 0506+056, located at redshift z ≈ 0.3365 ± 0.0010 (Paiano et al. 2018). The Major Atmospheric Gamma-ray Imaging Cherenkov (MAGIC) telescopes measured the flux at E > 9 × 10 10 eV indicating two periods of enhanced activity: MJD 58029.22 and MJD 58030.24 (Ansoldi et al. 2018). The hypothesis that the correlation between the electromagnetic and the neutrino signals happened by chance is rejected at a 3σ-level (IceCube et al. 2018).

ELECTROMAGNETIC CASCADES AND IGMFS
High-energy gamma rays may initiate electromagnetic cascades in the intergalactic medium. Aharonian et al. (1994) and Plaga (1995) suggested that they can be used to measure IGMFs. The underlying physics of this process is well known. A blazar emits high-energy gamma rays, which can interact with pervasive radiation fields such as the extragalactic background light (EBL) and the CMB, producing electron-positron pairs: γ + γ bg → e + + e − , with γ bg denoting the background (CMB/EBL) photon. Electrons and positrons upscatter CMB/EBL photons to high energies via inverse Compton (e ± + γ bg → e ± + γ). The high-energy photons produced will then restart the whole process, creating a cascade of particles. The cascade will stop when the energy of the photons drop below the kinematic threshold for pair production. The electrons and positrons are sensitive to the local magnetic field where they were produced. They are deflected in opposite directions, proportionally to the magnetic-field strength.
Most of the observed high-energy gamma rays from cosmologically distant objects are attenuated and the spectrum measured between GeV and TeV energies is a combination of both prompt and cascade gamma rays. The former are produced at the source and do not interact during propagation, whereas the latter are secondaries from high-energy primaries that underwent a cascade process in the intergalactic medium. The spectrum of TXS 0506+056 indicates that at least a fraction of the flux at ∼ GeV-TeV is comprised of secondary photons produced during propagation, as we will show later. We do not expect a significant emission in TeV and above because the high-energy part of the spectrum is absorbed by the EBL, given the distance of the object. Note that the spectrum should extend up to energies of at least ∼ 400 GeV, as observed by MAGIC.
The High-Energy Stereoscopic System (H.E.S.S.) and the Very Energetic Radiation Imaging Telescope Array System (VERITAS) have also observed this object, but only upper limits were provided (IceCube et al. 2018).
The flaring activity of TXS 0506+056 in a broad wavelength band started in 2017 June, lasted for about six months, and was preceded and succeeded by a quiescent period. The peak luminosity was reached around the time of detection of IC-170922A, decreasing slowly thereafter. It is tempting to directly correlate the veryhigh-energy gamma-ray signals observed by MAGIC with IC-170922A and to attempt to directly measure the strength of IGMFs as a function of the coherence length. However, multi-TeV gamma rays need not be produced simultaneously with the enhanced neutrino emission. In fact, they can be produced anytime during acceleration; nevertheless, the neutrino and TeV gamma-ray peaks lie within the same time window (Gao et al. 2019). Therefore, we fix the duration of the flare: ∆t flare ≈ 6 months.
In the absence of IGMFs, both the neutrino and the electromagnetic signals would be detected roughly within the same time interval, ∆t flare . Any time delay incurred by IGMFs would be shorter than the period of flaring activity, otherwise the light curves of gamma rays would be considerably different than those at other wavelengths. Thus, the observation of very-high-energy gamma rays in coincidence with the flare sets an upper bound on the maximum time delay due to IGMFs: ∆t IGMF < ∆t flare . This provides limits on the strength and coherence length of IGMFs.

MODEL AND SIMULATIONS
Here we employ three-dimensional Monte Carlo simulations. The development of electromagnetic cascades in the intergalactic medium is modelled using the CR-Propa 3 code (Alves Batista et al. 2016b) considering all relevant interactions and energy loss process, namely pair production, inverse Compton scattering, synchrotron emission, and adiabatic losses due to the expansion of the universe. The spectrum of gamma rays emitted by TXS 0506+056 is assumed to be a power law with spectral index −α and an exponential cut-off at E max , i.e., E −α exp(−E/E max ). Several scenarios were studied, fixing E max (10 10 ≤ E max /eV ≤ 10 14 ) and assuming 0 ≤ α ≤ 4. The following EBL models are used: Gilmore et al. (2012), Domínguez et al. (2011), and the upper and lower limits by Stecker et al. (2016). For all cases studied, secondary photons produced in the cascade contribute to the observed flux at ∼ GeV to at least a percent level. A number of scenarios for the magnetic field were considered, including all combinations of 10 −19 ≤ B/G ≤ 10 −14 and 10 −2 ≤ L c /Mpc ≤ 10 3 , both in logarithmic steps of 1, in addition to the case B = 0.
The neutrino emission takes place during a high state of the object. The spectral parameters of this period are not necessarily the same as the ones during the normal (low) state. Cascade photons stemming from gamma rays emitted during the low state may contribute to the observed flux at E 1 GeV, together with the flux from the high state. Thus, the spectrum of gamma rays effectively emitted by the source is: (1) wherein η denotes the flux enhancement in the high state (subscript 'h') with respect to the low state ('l'). This is computed within time interval ∆t flare , fixed by the duration of the neutrino flare. Another relevant parameter is the time scale over which TXS 0506+056 is a gammaray emitter at the low state, denoted by ∆t AGN . Typically, AGN activity times range between 10 6 and 10 8 years (Parma et al. 2002). We use ∆t AGN = 10, 10 4 , and 10 7 years.
We estimated how much of the total flux comprises secondary photons produced in the cascade. To this end, we assumed B = 0, and took only one of the components of equation 1, which is equivalent to setting η = 0. For the conservative case wherein E max = 10 11.5 eV and α < 3, our simulation results suggest that at 10 GeV, at least 10% of the flux correspond to cascade photons. This increases for stronger EBL models like the Stecker et al. (2016), upper limit. Therefore, with the observations by MAGIC extending up to E ∼ 400 GeV, we expect a sizeable contribution of secondary photons to the total flux, as shown by Saveliev & Alves Batista (2020).
Photohadronic and hadronuclear interactions in blazar jets can produce neutrinos, as well as gamma rays of similar energies. According to most models the maximum cosmic-ray energy that can explain the neutrino observations is E CR ∼ 10 16 eV (Ansoldi et al. 2018;Keivani et al. 2018;Murase et al. 2018;Liu et al. 2019;Gao et al. 2019). Consequently, neutrinos and gamma rays with E ∼ 10 15 eV should be produced. This energy can be significantly degraded if the source environment is opaque to high-energy gamma rays. We do not concern ourselves with this absorption; instead, our phenomenological model describes a gamma-ray flux that is injected into the intergalactic medium after escaping the object.

FIT RESULTS
We are now able to constrain the strength of IGMFs using information from both messengers, gamma rays and neutrinos. We first fit the spectrum for the low state. We find that, in this case, the spectral parameters remain virtually unaltered regardless of the magneticfield properties: α l = 2.2 and E max,l = 250 GeV if we only consider combinations of α l and E max,l for which the fit produces P-values p > 10 −3 . The second step is to use these values to scan the all the combinations of the parameters E max,h , α h , B, and L c . One example of the fitted spectrum is shown in figure 1.
In the scope of this work we are interested in the case of a non-vanishing magnetic field, such that as a first step we carried out a hypothesis test for each of the four considered EBL models, with the null hypothesis being B = 0. To do so, we marginalized over all other quantities, obtaining a probability distribution for the seven values of B simulated. We find that only for the models by Domínguez et al. (2011) and the lower limit model by Stecker et al. (2016) the null hypothesis can be rejected. For completeness, we also included the EBL models by Gilmore et al. (2012) and the upper limit by Stecker et al. (2016) in our considerations, which, based on the likelihood ratio analysis we carried out, do not disfavour the B = 0 case, but still allow for nonzero magnetic fields. This caveat should be borne in mind hereafter.
To constrain the magnetic field (B) and coherence length (L c ), we first marginalise our results over the spectral parameters. Then we derive two-dimensional marginalised confidence regions for B and L c , as shown in figure 2. The results for other values of ∆t AGN are very similar. In fact, we found this parameter to have little influence on the constraints. A summary with the best-fit intervals is shown in table 1. Note, again, that not all models allow us to reject the null hypothesis (B = 0). Our results shown in figure 2 provide seemingly weak constraints on the parameter space. For instance, with the EBL by Domínguez et al. (2011), the 90% contour disfavours very large coherence lengths (L c 100 Mpc) for fields weaker than B ∼ 10 −18 G. If IGMFs have galactic scales (L c 10 − 100 kpc), then B 10 −18 G (at 90% C.L.). Similarly, for the Stecker et al. (2016) lower-limit model, L c 10 kpc is not contained within the 95% confidence region.

DISCUSSION
IGMFs are commonly constrained using TeV-emitting extreme blazars (see e.g. Neronov & Vovk 2010, Taylor et al. 2011, Arlen et al. 2014, Ackermann et al. 2018). Due to EBL absorption, TeV fluxes from very distant objects are not observed, and only objects at redshifts z 0.20 are used for this purpose. TXS 0506+056 is located at z ≈ 0.34, so that any flux at TeV energies would be suppressed. Consequently, the high-energy flux would be shifted to lower energies, retaining some information about the injection spectrum and the intervening magnetic fields. Alone, these constraints would likely be weak, but the picture changes with the temporal information provided by neutrinos.
Our results, shown in fig. 2, are compatible with bounds derived by other authors (Neronov & Vovk 2010;Tavecchio et al. 2010;Taylor et al. 2011;Arlen et al. 2014;Ackermann et al. 2018). This is not surprising given that we could not exclude most of the parameter space probed. Nevertheless, we successfully derived some constraints on the coherence length with our de-tailed three-dimensional treatment of gamma-ray propagation, whereas other works commonly rely in the smallangle approximation for calculating angular and temporal profiles (see e.g. Ackermann et al. 2018). Moreover, in these works the magnetic fields have a cell-like structure, whereas we consider a more realistic Kolmogorov turbulent field (details can be found in Alves Batista et al. 2016a).
Similar analyses could be performed in a more clear fashion using GRBs (Takahashi et al. 2008(Takahashi et al. , 2012. In this case, there would be no low state and the number of free parameters in the fit would decrease. This analysis was carried out recently for GRB 190114C (Wang et al. 2020).
The IceCube Collaboration (2018) reported activity near TXS 0506+056 in 2014/2015. An analysis of Fermi-LAT data does not provide any indications of enhanced activity around this period (Garrappa et al. 2019). Otherwise, we could have applied this same method to constrain IGMFs.
Plasma instabilities may quench electromagnetic cascades propagating in the intergalactic medium (Broderick et al. 2012;Sironi & Giannios 2014;Broderick et al. 2018). This is, however, under dispute (see e.g. Miniati & Elyiv 2013). Nevertheless, even if plasma instabilities are taken into account, this should not significantly affect our estimates because of the transient nature of the phenomenon, which would not allow the instability to grow fast enough over the brief period of enhanced activity (Alves Batista et al. 2019).
Coherence lengths of ∼ 10 kpc are disfavoured at a 90% confidence level for the two EBL models shown in fig. 2. This would disfavour models in which the IGM is magnetised by galactic winds (e.g. Bertone et al. (2006)), since they predict L c ∼ 1 − 10 kpc. Models in which IGMFs were generated by cosmic rays escaping galaxies prior to reionisation (e.g., Miniati & Bell (2011)) are only marginally compatible with our results. Fields seeded by AGNs are expected to have ∼Mpc scales (Durrer & Neronov 2013), and are well within the estimated bounds. We have not probed L c 10 kpc, a region of the parameter space that would allow us to constrain some cosmological magnetogenesis models (Durrer & Neronov 2013).
The detection of high-energy neutrinos correlating with the electromagnetic signal favours hadronic models (Ansoldi et al. 2018;Gao et al. 2019;Keivani et al. 2018) and strengthens the case for multi-TeV gammaray production in the blazar. The bounds presented here are rather robust and rely only on the assumption that the gamma rays at E 1 GeV observed by Fermi-LAT, and the highest-energy bin observed by MAGIC at E ≈ 400 GeV, are produced during the neutrino flaring period. Therefore, the reliability of our limits reflects the significance of this correlation.
The multi-messenger approach used to constrain IGMFs based on time delays between high-energy gamma rays and their counterparts (in neutrinos and other wavelengths) is powerful. It differs from the methods commonly found in the literature (Plaga 1995;Neronov & Semikoz 2009;Neronov et al. 2013) in that we do not attempt to correlate the GeV and TeV emissions with each other. Instead, we use the time delay between the neutrino and the gamma-ray signals. For this reason, this method enables us to probe the universe up to high redshifts, since no TeV signal is required to perform these estimates and we can evade the limitations placed by the EBL attenuation of the very-high-energy part of the spectrum.
The same idea used here could be applied whenever a high-energy gamma-ray signal correlates with the arrival of high-energy neutrinos, since this guarantees the production of very energetic gamma rays in the source, though there is no guarantee that they can escape the environment. In the absence of neutrinos, the time delays between the cascade photons and another messenger such as gravitational waves (e.g. from mergers of compact objects) would rely on the existence of a putative multi-TeV emission which, in the case of high-redshift sources, cannot be observed. Either way, gamma rays in the ∼GeV band, together with another messenger (neutrinos, gravitational waves, or ∼TeV gamma rays) can be used to place limits on both the strength and the coherence length of IGMFs.

SUMMARY AND OUTLOOK
We have here derived combined constraints on both the coherence length and the strength of IGMFs using three-dimensional simulations. On one hand, the constraints we derived are relatively weak compared to the existing ones, given that we used only one object. On the other hand, this is the first time that bounds on the coherence length are obtained, which is of utmost importance for understanding the origin of IGMFs.
We showed that the intrinsic spectral parameters of the object are degenerate with respect to the magneticfield ones. It follows that magnetic-field effects may be important when fitting the high-energy region of spectral energy distributions. We investigate this in more detail in another work (Saveliev & Alves Batista 2020).
Our results exclude the hypothesis of null IGMFs for only two EBL models, out of the four tested. For these models, considering the range of parameters studied, we could obtain bounds on the coherence length. The upper bound was inferred to be L c 300 Mpc. The lower bound depends on the EBL: L c 30 kpc for the more general case and, more interestingly, L c 300 kpc for a weak EBL, at a 90% C.L.