The Gamma-Ray Window to Intergalactic Magnetism

One of the most promising ways to probe intergalactic magnetic fields (IGMFs) is through gamma rays produced in electromagnetic cascades initiated by high-energy gamma rays or cosmic rays in the intergalactic space. Because the charged component of the cascade is sensitive to magnetic fields, gamma-ray observations of distant objects such as blazars can be used to constrain IGMF properties. Ground-based and space-borne gamma-ray telescopes deliver spectral, temporal, and angular information of high-energy gamma-ray sources, which carries imprints of the intervening magnetic fields. This provides insights into the nature of the processes that led to the creation of the first magnetic fields and into the phenomena that impacted their evolution. Here we provide a detailed description of how gamma-ray observations can be used to probe cosmic magnetism. We review the current status of this topic and discuss the prospects for measuring IGMFs with the next generation of gamma-ray observatories.


Introduction
The advent of imaging air Cherenkov telescopes (IACTs) enabled the study of very-high-energy (VHE; E 1 TeV) processes involving gamma rays with unprecedented precision. With small angular resolutions (θ psf ∼ 0.1 • ), IACTs such as the High-Energy Stereoscopic System (H.E.S.S.) [1], the Major Atmospheric Gamma Imaging Cherenkov (MAGIC) [2,3], and the Very Energetic Radiation Imaging Telescope Array System (VERITAS) [4,5] provide a unique view of the gamma-ray universe above TeV energies. These observations are supplemented, at higher energies, by measurements with water-Cherenkov detectors such as the High Altitude Water Cherenkov Experiment (HAWC) [6], the Astrophysical Radiation with Ground-based Observatory at YangBaJing (ARGO-YBJ) [7], the Large High Altitude Air Shower Observatory (LHAASO) [8], and the Tibet Air Shower Experiment (Tibet-ASγ) [9]. The launch of the Fermi Large Area Telescope [10] (Fermi-LAT) in 2008 was undoubtedly one of the most important landmarks in gamma-ray astrophysics. The complementarity between Fermi-LAT and IACTs has been crucial to glimpse into extreme cosmic accelerators, and to shed light on large-scale properties of the Universe, including the topic of this review: intergalactic magnetic fields (IGMFs).
In the last decades, active galactic nuclei (AGNs) were observed across the whole electromagnetic spectrum, from radio to gamma rays (for reviews see, e.g., Refs. [11,12]). Ever since the Whipple Telescope observed the first BL Lac-type AGN at very-high energies, Mrk 421, in 1992 [13], blazars -a sub-class of AGNs -have been extensively studied. Because their relativistic jets point approximately towards Earth, their emission can probe the Universe over vast distances. At gamma-ray energies, in particular, they can be used to probe the extragalactic background light (EBL) [14][15][16][17][18][19][20] and IGMFs [21][22][23], as well as fundamental physics [24][25][26]. It is thanks to combined observations of blazars by IACTs and Fermi-LAT that the first studies aiming to constrain IGMFs were possible.
The process whereby high-energy gamma rays emitted by blazars initiate electromagnetic cascades in the intergalactic space has been known for over half a century [27][28][29][30][31]. It follows immediately from the idea of cascades that magnetic fields can interfere with their development. Despite the numerous works on the topic (e.g., [32-37]), it was not until later that the potential of electromagnetic cascades as a method to probe IGMFs was fully realised by Plaga [38], though there have been considerations of the underlying concept even before that [39].
The seminal work of Neronov & Vovk [40] sparked an avalanche of subsequent investigations along the same lines in the following years (e.g., [41][42][43][44][45][46][47]), most of which derived lower bounds on the strength of IGMFs ranging from B 10 −18 G to B 10 −16 G, depending on the specific details of the analysis, which is in line with more recent works [47][48][49]. So far, only constraints on IGMFs have been derived, as opposed to actual measurements. Apart from very general considerations, the coherence length (L B ) of IGMFs has filaments, whose fields are B ∼ 0. 1-10. nG [52,57,65]. They compose the cosmic web, whose magnetic properties are poorly known. This is, to a large extent, due to the scarcity of observational data, owing to the intrinsic difficulties in measuring magnetic fields at scales larger than clusters of galaxies. For this reason, numerical simulations play an important role providing the full picture of how magnetic fields are distributed in the cosmic web. For further details, the reader is referred to, e.g., Ref. [65].
A natural question that arises is how magnetic fields in galaxies reached the µG level we observe today. One possible explanation is that astrophysical dynamos can amplify seed magnetic fields by many orders of magnitude. In this context, these seed fields are required to have strengths of at least 10 −22 to 10 −15 G [69,70], the actual value depending on the particular model. However, if the seeds are strong enough (B 10 −11 G), one does not need to invoke dynamos. In this case, adiabatic compression [55] (potentially together with some stretching and shearing of flows [71]) is sufficient.
Given the distance scales involved in typical IGMF studies using particle probes, only the large-scale distribution of magnetic fields is relevant. In this case, clusters of galaxies fill 10 −3 of the Universe's volume (the so-called volume filling factor), such that their magnetic field is virtually negligible if we are studying how particles propagate over large distances and the effects of magnetic fields upon them. Filaments are believed to have filling factors ∼ 10 −3 -10 −1 , whereas cosmic voids are the most important contribution, with a filling factor 10 −1 . Therefore, magnetic fields in the voids are, to first order, the dominant component that determines how particles propagate over cosmological baselines.
The origin of the seed magnetic fields is one of the main open questions in astrophysics today. In Sec. 2.2, we briefly mention some of the main mechanisms for magnetogenesis, focusing on providing some estimates of the relevant observables -field strength, coherence length, and helicity -that can be probed with high-energy gamma-ray observations. Before that, in Sec. 2.1, we provide some important definitions and the mathematical framework required for understanding cosmic magnetism. Conclusively, in Sec. 2.3 we present techniques other than gamma-ray astrophysics used to constrain IGMFs and the results originating from them.

Statistical Observables
Stochastic magnetic fields can be characterised by a number of observables which correspond to different statistical averages. The first one is the average magnetic-field strength (B). When considering this quantity, it is a common misconception to talk about the mean of B because at cosmological scales it is expected that B i = 0 for each individual component i, in particular for a Gaussian distribution, which is the typical first-order assumption. The relevant quantity, in this context, is the root mean square of the magnetic field, defined from where V is the considered volume. Another related quantity, the magnetic helicity H B , is given by where A(r) is the vector potential, i.e. B = ∇ × A. Originally, H B had been defined for a vanishing normal magnetic field component everywhere at the boundary of V, even though it is possible to drop this condition in a more general case [72]. As the name suggests, H B is directly related to the topology of the magnetic field, more precisely to whether the magnetic field on average is left-or right-handed. It should be noted that magnetic helicity is a well-defined quantity as it is gauge-invariant if the aforementioned requirement of a vanishing normal magnetic field at the border is fulfilled. Furthermore, it is conserved in ideal MHD and plays an important role for the time evolution of the (energy content of the) magnetic field in general (cf. Sec. 2.2.1).
To define the other quantities, we need to introduce the Fourier transform, which for a given (magnetic) field B(r) is given by and represents the mode for the wave vector k. We can then determine the statistical connection between any two modes, represented by wave vectors k and k , by calculating the corresponding ensemble average (denoted as ... ), given by Assuming that the magnetic field is homogeneous and isotropic, the general form of P ab is where δ ab and abc are the Kronecker delta and the Levi-Civita symbol, respectively, M k is the spectral magnetic energy, H k is the spectral magnetic helicity density, and c H is a numerical constant which depends on the convention used. It is important to mention here that there is a fundamental relation between M k and H k : one can show (see, for example, [73]) that for a given k the value of H k is limited by the corresponding value of M k , such that the right hand side of Eq. 6 is also called maximal (spectral) helicity, and the actual spectral helicity density may be expressed as a fraction f H of it, with −1 ≤ f H ≤ 1.
In general, one assumes M k ∝ k α B −1 , which means that the spectrum is given by a power law for which the spectral index α B defines its type. It is assumed that at small scales (i.e., for large values of k), IGMFs have a Kolmogorov (α B = −2/3, see [74,75]) or an Iroshnikov/Kraichnan (α B = −1/2, see [76,77]) spectrum at present time. Both values for the spectral index were derived from dimensional considerations, with the latter one assumed to be valid under the assumption of a strong mean magnetic field. Still, due to the fact that the numerical values of these two spectral indices are very close to each other, up to the present day it has not been possible to distinguish between them experimentally [78,79]. For large scales (i.e. small k) one expects a Batchelor spectrum (α B = 5), as predicted using general causality arguments in [80] and confirmed by semi-analytical simulations in [81]. Other works also considered a white noise spectrum (α B = 3) [82,83]. On the other hand, IGMFs produced during Inflation (cf. Sec. 2.2.1) are expected to be scale-invariant, which corresponds to α B = 0 (see, for example, [84]). Note that there are different ways to define the spectral index, such that the numerical values in other publications might differ from the one used here while describing the same kind of spectrum. The last essential characteristic statistical observable considered here is the correlation length (L B ) which is given by 6 of 57 In a simplified way, L B can be understood as the average size of the eddies of the magnetic field. Again it should be noted that several different ways to define the correlation length are found in the literature, such that small differences (for example by a factor of a few) are possible. This is discussed, e.g., in [85], where also the power-law case relevant here is addressed in more detail.

Origin
While the origin of IGMFs is still unknown, there are two classes of scenarios for their magnetogenesis present in the literature [55,58,64]. In cosmological scenarios strong seed magnetic fields were created during the early Universe and later decayed to their present state. In astrophysical scenarios, on the other hand, weak seed magnetic fields emerged due to local effects (for example, a battery process) in astrophysical objects, being subsequently amplified by dynamo mechanisms. In the following, we will present possible mechanisms for both of these classes of scenarios. Note that this separation is done here purely for reasons of clarity and comprehensibility. In reality, the situation may be more complex, as the particular mechanism may be the result of (yet) unknown physics, or the actual origin of IGMFs may turn out to be a combination of multiple processes, astrophysical and/or cosmological.

Cosmological Scenarios
The seed magnetic fields of cosmological scenarios, i.e., the PMFs, are thought to be created by some major cosmological effect, such that they permeate the whole Universe. Without any claim to completeness (more details may be found in [58,59,86]), we list some of these possibilities below.
Inflation. Magnetic fields may have been produced during inflation (for a review on Inflation in general, see, e.g., [87]). However, if Maxwellian conformal invariance is preserved, these fields are predicted to be exceedingly weak (B 10 −50 G at the epoch of galaxy formation [86]), being negligible for all practical purposes [88]. Models for inflationary magnetogenesis that are of astrophysical relevance must generate much stronger fields. Because conformally invariant fields are not produced in an expanding conformally-flat spacetime, one has to introduce a coupling of the electromagnetic field with the inflaton, and/or an additional coupling which breaks the conformal or gauge invariance, mainly of the form R µναβ F µν F αβ or R µν A µ A ν , respectively [58] (where F αβ is the electromagnetic field tensor, R µν is the Ricci tensor and R µναβ is the Riemann curvature tensor), even though other terms are also possible [88,89]. After the seminal publications in the field [88,90] the follow-up works (see, for example, [91][92][93][94][95][96][97][98][99]) then further explored the idea or investigated more exotic scenarios. Due to a large parameter space the resulting magnetic field strength estimations, even in the simplest models, range from 10 −65 to 10 −9 G [92].

Post-Inflationary.
It is also possible that magnetic fields emerged between Inflation and the EWPT, for example during or before reheating. The general idea is that the coupling between the electromagnetic and a scalar field breaks the Maxwellian conformal invariance. In particular, the scalar field in question may be an oscillating inflaton, which decays into radiation and reheats the Universe [100], resulting in IGMFs with B 10 −15 G on ∼ Mpc scales. In another scenario [101], Majorana neutrino decays result in lepton asymmetries, and ultimately in baryon asymmetries via anomalous processes, subsequently leading to the violation of lepton/baryon numbers. This then may produce relic hypercharge magnetic fields which are converted to electromagnetic fields during the EWPT, giving ∼ 10 −18 G field strength with L B 10 pc today. More recently, the idea of a Weibel instability emerging and subsequently amplifying a possible inflationary magnetic field during this era has been considered [102].
Electroweak Phase Transition. Within the SM, the EWPT is assumed to be rather smooth [103], such that in order to realize a first order transition, mechanisms beyond the Standard Model (BSM) have to be considered [104]. The basic idea of magnetogenesis during the EWPT was first laid out by [105]. Due to the restrictions of possible values of the vacuum expectation value of the Higgs field, Φ, which breaks the electroweak symmetry, and the fact that it varies with the position in space, we have ∂ µ Φ = 0, such that the electromagnetic field strength does not necessarily compensate effects arising from the Higgs field. Hence, we expect a non-vanishing magnetic field after the phase transition. Note, however, that the magnetic field depends on gradients of Φ. Other possible scenarios may be found in [106], with the general conclusion that magnetic fields of up to ∼ 10 −11 G on scales of ∼ 10 kpc are possible [59].
Quantum Chromodynamics Phase Transition. In a similar fashion to the EWPT, it should be mentioned here that within the Standard Model of particle physics (SM) the quantum chromodynamics phase transition (QCDPT) is considered to be of the second order or crossover type [107][108][109][110], such that, for it to be of first order, a SM extension has to be invoked [111]. Several works [112][113][114] discuss magnetogenesis due to the growth of bubbles of the hadronic phase and, subsequently, charge separation, which ultimately leads to the creation of electric currents and consequently of magnetic fields with an estimated field strength of the order of ∼ 10 −16 G on ∼ kpc scales [113].
It is usually assumed that immediately after magnetogenesis most of the magnetic-field energy is concentrated on a characteristic length called the integral scale. The basic idea, as described in [115], is that throughout the evolution of the Universe up to the present day, the magnetic energy decays starting with small scales, such that the integral scale is increasing until it reaches the coherence length of IGMFs today. Throughout the years, there has been a large number of simulations, both numerical and (semi-)analytic, which modelled this time evolution for different magnetogenesis scenarios [81,83,[115][116][117][118][119][120].
As a final remark it should be pointed out that, due to the fact that magnetic helicity is (nearly) conserved, it plays an important role in the time-evolution of magnetic fields, in particular by causing the so-called inverse cascade of energy, i.e. the transfer of magnetic energy from small to large scale fluctuations [73,115,[121][122][123]. These inverse cascades, however, do not seem to be exclusive to helical fields, as shown in recent simulations [124].
It is well possible that PMFs actually were helical. One of the first works along these lines was [125], suggesting the creation of a left-handed PMF due to a change of the Chern-Simons number. Other possible mechanisms include extra dimensions [126], the coupling to the cosmic axion field [127] or an axion-like coupling [128], the Riemann tensor [129], a spectator field [84], or an inherently helical coupling [130] during Inflation in the first place. Recently, also the possibility of helicity generation via a chiral cosmological medium around the EWPT has been considered [131], however the authors found the effect to be suppressed due to the value of the baryon to entropy ratio.

Astrophysical Scenarios
A number of possible mechanisms also exists for the astrophysical scenario. They all have in common that magnetic fields are created locally due to some astrophysical process. Some of them are concisely described below.
Biermann Battery. It is manifestly difficult to create magnetic fields from scratch due to the fact that in classical MHD, if B(r) = 0 at some instant in time, then this is true for all later times. A way to evade this limitation is through the Biermann battery mechanism [132,133], for which the basic idea is that the misalignment of temperature and density gradients induces an electric field which ultimately results in the generation of a magnetic field. Prior to Reionisation, this process produces exceedingly weak fields in the intergalactic space (B 10 −24 G) [134]. In protogalaxies, these fields can reach B ∼ 10 −22 -10 −17 G [135][136][137]. For other astrophysical and cosmological settings see also [56,57,[138][139][140]. Galactic Outflows. One obvious candidate to produce IGMFs are the galaxies themselves, as they eject matter and energy into the intergalactic space. Most authors assume that this can be driven by stars, in particular magnetised winds, or cosmic rays [66,[141][142][143][144]. However, other possibilities exist as well, including the magnetisation of voids by giant radio lobes or bubbles from AGNs, even though energetics requirements generally do not allow for such a substantial effect over the age of the Universe [68,[145][146][147].
Cosmic-Ray Return Currents. In addition to the outflow scenario discussed above, cosmic rays escaping from a galaxy create a charge imbalance resulting in electric fields and, subsequently, return currents. Ultimately, these return currents can produce magnetic fields on scales which are sufficiently large to provide the seed for IGMFs [148,149].
Photoionisation during the Reionisation Era. During Reionisation high-energy photons are able to escape from objects like population III stars, protogalaxies, and quasars into the (then) neutral intergalactic medium (IGM). This causes photoionisation which ultimately causes the generation of radial currents (and electric fields), inducing magnetic fields with strengths B ∼ 10 −25 − 10 −20 G on scales between ∼ 1 kpc and 10 Mpc [150][151][152]. Remarkably, this mechanism can generate global magnetic fields through astrophysically-initiated mechanisms. This seeding scheme agrees with results of large-scale cosmological MHD simulations by Garaldi et al. [153], although they could, in principle, be subdominant with respect to seeds produced through the Biermann battery.
Primordial Vorticity. In a seminal paper by Harrison [154], it was suggested that due to relativistic effects electromagnetic fields are coupled to vorticity 2 , such that rotating protogalaxies could create primordial vorticity that could generate magnetic fields in the radiation-dominated era. However, vorticity is predicted to decay rather fast in the early Universe, such that more advanced theories based on the same idea, but with vorticity appearing at later stages or using higher-order effects, had to be introduced [155][156][157][158].
Several of the mechanisms listed above require a dynamo mechanism in order to amplify the magnetic field strength to the observed present-day values. Especially the small-scale dynamo has attracted major interest in this context (see [65,73,159,160] for some recent results), even though simulating it numerically poses a challenge due to the size of the scales which have to be resolved.

General Constraints
In this review we focus on constraints on IGMFs from gamma-ray observations. However, since IGMFs can interact through various electromagnetic phenomena throughout the Universe, there are other ways to derive bounds on them. In this section we present the general ideas to do so, based on Ref. [161] and including some more recent developments.
First, there is a generic lower and upper limit on the coherence length (L B ) [161]. The latter is given by the size of the observed Universe, i.e., the Hubble radius. On the other hand, the IGMF decays due to magnetic diffusion, such that the lower limit on L B is given by the length scale equivalent to the corresponding decay time, i.e., L B 2 × 10 11 m [55].
As for the magnetic-field strength (B), measurements of the Zeeman splitting of H I lines can be used to set upper bounds on this quantity. This can be done either for our own galaxy or for the radiation from distant quasars [164][165][166], both consistently giving a result of the order of ∼ µG. In the latter case, any 2 This mechanism can be viewed as cosmological, since it involves density perturbations. However, the necessary conditions for the vorticity generation involve protogalaxies, so we chose to classify it as an astrophysical magnetogenesis model.   [161]. The upper bounds are due to Zeeman splitting and Faraday rotation observations of extragalactic objects [161]. The 'early magnetic dissipation' bound indicates the region of the parameter space excluded by freely decaying MHD in the early Universe [58,115]. Other limits from cosmology come from CMB observations (spectrum [58,162] and anisotropies [163]); the currently strongest limit [60], labelled 'JS19', is shown for the case of a scale-invariant spectrum (α B = 0) which leads to the most conservative bounds.
stronger IGMFs along the line of sight to the object would have a measurable impact on the observations, thus giving a robust upper limit for the IGMF strength.
Another constraint on IGMFs is derived from Faraday rotation measurements of polarised radio emission from quasars and other extragalactic sources. Faraday rotation describes the (wavelength-dependent) rotation of the polarisation plane of polarised electromagnetic radiation when it traverses a magnetised medium. Therefore, the value of the relevant observable, the so-called rotation measure (RM), may be subdivided into contributions from the host galaxy, the IGM, and the Milky Way. With a rigorous statistical analysis of RM data, one can then identify the impact of the IGMF, and hence derive limits on the IGMF strength which in general depend on L B . There are many studies on the topic [167][168][169][170][171][172][173][174], all of which give upper limits ranging from nG to a few µG. This is also confirmed by other methods, like the interpretation of radio observations as the result of shock acceleration in galaxy clusters [175,176]. In this context, fast radio bursts (FRBs) can play an important role [177], delivering both rotation and dispersion measures. As a consequence, magnetic fields along the line of sight can be better inferred because the use of these two observables reduce the reliance on models on models of the electron density distribution [178]. RMs can be related to the magnetogenesis model, as shown in Ref. [179].
In addition, an important set of limits can be derived from cosmology considering the cosmological scenarios for magnetogenesis (cf. Sec. 2.2.1). An indirect, theoretical approach is to consider a given mechanism of magnetic-field generation and derive the corresponding limits on the initial magnetic-field strength and correlation length, and then calculate their time evolution via freely-decaying MHD up to the present day. A detailed description is given in [115], while in Fig. 1 we present the region which contains most of these constraints, following [58]. In general, one can state that these limits bound the field strength from above and the coherence length from below, the latter due to the fact that in cosmological scenario IGMFs are generated at small scales (see above).
From an observational point of view, most of the limits on IGMFs from cosmology are derived using the cosmic microwave background (CMB), as there is a large range of effects through which magnetic fields can impact the background radiation. The most basic idea, developed already in [180], is to assume a homogeneous field throughout the Universe and then to derive the temperature anisotropies expected from that. Comparing this dataset to data from COBE [163] or, more recently, from Planck [181], gives an upper limit of around 4 nG (marked in Fig. 1 as "CMB anisotropies"). Since then the upper limit has been dramatically improved by using CMB observations in combination with such effects as spectral distortions (in Fig. 1 we present a limit stemming from this phenomenon based on [162], denoted "CMB spectrum"), temperature anisotropies, polarisation, non-Gaussianity and Reionisation (for a concise review see [60]). The best upper limit so far, B ∼ 10-50 pG [60], was derived by considering the change of the Recombination process itself via density fluctuations due to the presence of PMFs. In addition, CMB observations are also interesting because they may also be used to derive constraints on magnetic helicity [182][183][184].
Finally, ultra-high-energy cosmic rays (UHECRs), i.e. nuclei with energies above 10 18 eV, may be used to constrain IGMFs [161,[185][186][187][188][189]. The general principle used here is that, since UHECRs are charged, they are deflected. Hence, once their sources are identified, the corresponding deflection angle can be measured, providing a direct measure of the magnetic-field strength orthogonal to the line of sight. Ref. [190] used the observed excess of UHECRs with energies ∼ 10 20 eV in the direction of Centaurus A to constrain the local extragalactic magnetic field, obtaining B 10 −8 G. This local constrain evidently serves as an upper limit for IGMFs. Ref. [191] used the anisotropy reported by the Pierre Auger Observatory to associate UHECR detections with extragalactic objects and to derive upper limits of B ∼ 10 −9 G for L B < 100 Mpc and ∼ 10 −10 G for L B > 100 Mpc. More recently, Ref. [192] found that for B > 6 × 10 −10 G the Auger anisotropy measurements are in good agreement with the local density of star-forming galaxies. On the other hand, if the local density is treated as a model parameter, the authors found a conservative upper limit of BL 1/2 B < 24 nG Mpc 1/2 . In principle UHECRs may also be used to constrain the helicity of IGMFs, as argued in [193] and demonstrated numerically in [194].
Note, however, that with the direct UHECR observations available today, it is rather difficult to derive IGMF constraints, as their sources would have to be known (see also Secs. 3.1 and 4.5 for an indirect gamma-ray-based approach on how to use UHECRs for deriving IGMF constraints). Moreover, the statistic of events at the highest energies (E 4 × 10 19 eV) is fairly limited, while the composition (and thus the charge) of UHECRs is only known statistically, and not on an event-by-event basis [195], posing severe challenges for any attempt to constrain IGMFs with UHECRs. Finally, the distribution of magnetic fields in the cosmic web is more complex than that in cosmic voids, and much more uncertain [65]. Numerical studies of UHECR propagation in magnetic fields lead to very discrepant results, such that the prospects for UHECR astronomy (and thus IGMF constraints using UHECRs) are far from clear (cf. [196][197][198][199]).

Electromagnetic Cascades
In this section we lay out the theoretical foundations for understanding how electromagnetic cascades develop in intergalactic space. We start off, in Sec. 3.1, by describing some classes of astrophysical objects that can emit particles that initiate the electromagnetic cascades. We describe two scenarios, depending on the type of particle that initiate the cascade which ultimately leads to the observed gamma-ray signal. In the first, the cascades are triggered by high-energy gamma rays (or electrons), whereas in the second, they are initiated by ultra-high-energy cosmic rays. After describing how electromagnetic cascades originate, we proceed to Sec. 3.2, where we give a detailed account of how they develop, how they interact with the photon fields that pervade the Universe, and how IGMFs can affect them. In Sec. 3.3 we present approximate analytical descriptions for the cascade process, which can also be treated in more detail with numerical codes, as described in Sec. 3.6 and illustrated in Sec. 3.7. In Sec. 3.4 we chime into the debate surrounding the role played by plasma instabilities in the development of cascades. Other potentially relevant propagation effects are concisely mentioned in Sec. 3.5.

Origin
The most common source of high-energy gamma rays used in IGMF studies are blazars. These objects are a sub-class of active galactic nuclei (AGNs) whose relativistic jets point approximately towards Earth [200]. Their spectral energy distribution is characterised by a low-energy hump corresponding to synchrotron emission by relativistic electrons [11,12]. There is also a second notable hump which, in the case of high and extreme synchrotron-peaked objects, is of interest for IGMF studies, peaks at ∼ TeV energies [12]. These objects are excellent cosmological probes as the very-high-energy emission assures the production of a substantial electronic component in the cascade, which can be used to probe IGMFs and the EBL [201].
In general, AGNs are active over time scales of T ∼ 10 6 -10 8 years [218], which makes it hard to use temporal information for constraining IGMFs since they are, for all practical purposes, quasi-steady sources. However, some objects, such PKS 2155-304 [214], Mrk 421 [13,219] and Mrk 501 [211,220,221], display short-time variability [222]. This information can, in principle, be used together with light curves in other wavelengths in the context of multimessenger campaigns to improve the constraints on IGMFs via time delays (see Sec. 4). Interestingly, for blazars that are slightly misaligned with respect to the line of sight, the GeV gamma rays stemming from the TeV emission could still be observed today over angular scales of ∼ 1 • even if the objects are no longer high-energy emitters [223].
Another class of objects that can potentially be used to probe IGMFs are gamma-ray bursts (GRBs). They emit highly collimated relativistic jets of high-energy radiation within a short time. GRBs are the most luminous events known, reaching isotropic-equivalent luminosities of ∼ 10 54 erg s −1 (see, e.g., [224][225][226] for reviews). Only recently were GRBs observed at very-high energies, with the detection of a bright flash from GRB 190114C [227], which was used for IGMF studies [228,229].
GRBs are interesting cosmological probes because they can be used exactly in the same manner as blazars, while in general providing more accurate temporal information. In this case, the highand very-high-energy components depend strongly on the properties of intervening IGMFs [230][231][232]. Moreover, if their HE light curve were known, in principle it would be possible to reconstruct a possible TeV emission even in the absence of VHE measurements, based only on the cascade signal at ∼ GeV energies, up to high redshifts [233][234][235]. Note that this argument only holds if the TeV light curve is known from theoretical models, which is not the case [224,226], or if there are well-defined relations between the GeV and TeV light curves.
The shape of the intrinsic spectrum of the sources of interest for this work, whether a blazar or a GRB, is not precisely known. In general, it is assumed to be a power law of the form where f cut (E) denotes a function that suppresses the spectrum above a given energy E max , which depends on the mechanism responsible for particle acceleration (and consequently for gamma-ray emission). This function is typically an exponential, log-parabola, or similar [236][237][238][239][240]. Interestingly, the value of E max that could be inferred with observations depends on the opacity of the Universe to gamma-ray propagation, i.e., the distribution of photon fields such as the EBL, as well as on the properties of the intervening IGMFs [241].
Cosmic-ray-induced electromagnetic cascades in the intergalactic medium may lead to observational signatures that resemble those initiated by gamma rays. These cascades are evidently affected by intervening IGMFs, as discussed in, e.g., Refs. [242][243][244][245][246][247]. Therefore, gamma rays from cosmic rays can, in principle, also be used to probe IGMFs. In the case of GRBs, this was suggested by the authors of Ref. [248]. Similarly, blazars are prominent contenders to emit UHECRs that can induce electromagnetic cascades in the IGM [249,250]. For highly collimated jets, it could be even possible to distinguish this hadronic scenario from the standard picture wherein gamma rays from the sources induce the cascades [251].
The cosmic rays relevant for this type of analysis are UHECRs, since they can produce electromagnetic cascades during intergalactic propagation, via photonuclear or hadronuclear interactions. In fact, this type of scenario has been suggested to explain observations of some blazars, as they lead to better agreement with the measurements [247,249,[252][253][254][255][256].
One process that creates electrons and positrons that trigger cascades is Bethe-Heitler pair production:

wherein
A Z X denotes an arbitrary cosmic-ray nucleus X of atomic mass A with Z protons interacting with a background photon (γ bg ).
Nuclear interactions also produce electrons and photons, starting with the photodisintegration of cosmic-ray nuclei (e.g., The most important hadronic channel for the generation of cascade-inducing particles (electrons and photons) is photopion production. For a cosmic-ray proton, p + γ bg → ∆ + → p + π 0 and p + γ bg → ∆ + → n + π + . The decay of the neutral pion produces photons (π 0 → γ + γ) and the decay of charged pions 3 lead to the generation of leptons (π + → µ + + ν µ ), including muons, whose decays produce electrons (µ + → ν e +ν µ + e + ). Note that the cascades stemming from the by-products of pion decays also occur for an arbitrary nucleus A Z X. In this case, the production rate depends on the number of each nucleonic species (see, e.g., Ref. [257] for further details).
While it is, in principle, possible to constrain IGMFs with UHECR-produced gamma rays, this is not straightforward. Firstly, the sources of UHECRs are not known. Secondly, they are deflected by intervening Galactic and extragalactic magnetic fields, potentially spoiling any correlation between the source direction and the gamma rays. For more details on the cosmic-ray-gamma-ray connection, the reader is referred to some reviews on the topics: [195,258].

Theory of Propagation
The particle physics aspects relevant for the propagation of high-energy gamma rays are well known. At energies E 400 GeV high-energy gamma rays interact with background photon fields predominantly at infrared frequencies, generating electron-positron pairs: γ + γ bg → e + + e − . The mean free path for this process is typically of the order of tens to hundreds of Mpc. These pairs up-scatter photons from (mostly) the CMB to high energies via inverse Compton scattering (e ± + γ bg → e ± + γ). These new photons, in turn, can either travel straight to Earth or, if their energy is above the threshold for pair production, restart this process, leading to an electromagnetic cascade in the intergalactic medium.
The picture outlined in the previous paragraph is theoretically simple, but there are uncertainties that complicate the modelling of electromagnetic cascades in the IGM. The most important one is the distribution of the EBL, which is not precisely known. At extremely high gamma-ray energies (E 10 17 eV), the contribution of the cosmic radio background (CRB) starts to become relevant. A comparison of EBL and CRB models, as well as the density of CMB photons, is illustrated in Figure 2.
In general, the inverse of the mean free path λ for a particle of energy E and mass m interacting with isotropically-distributed photons of differential number density 4 dn(ε,z) where z is the redshift (see below), ε refers to the energy of the background photon, and F depends on the process of interest, with kinematic limits s min and s max . For pair production, F = sσ PP (s), with s min = 4m 2 e c 4 and s max = 4Eε. For inverse Compton scattering, F = σ IC (s − m 2 e c 4 )/β, wherein β denotes the speed of the electrons, in units of the speed of light. The kinematic limits, in this case, are s min = m 2 e c 4 and s max = m 2 e c 4 + 2Eε(1 + β). Here, σ PP and σ IC denote, respectively, the cross sections for pair production and for inverse Compton scattering. Note that the minimum and maximum energies are, in principle, unbounded, i.e., ε min → 0 and ε max → ∞, but in practice they quickly vanish outside a given energy range. In the case of the EBL, for example, for purposes of calculations, ε min 10 −4 eV and ε max 10 eV (see Fig. 2). 3 There are other decay channels. For the purposes of this review, we present only the most relevant one. One example is the electronic mode (π + → e + + ν e ) that occurs much more rarely ( 10 −3 ) than the main one. 4 The differential photon number density is defined in a way such that  Following Ref. [161], we can approximate Equation 10 as for pair production, where κ is a model-specific parameter of the order of κ ∼ 1, and as for inverse Compton scattering. Typically, gamma rays with ∼ TeV energies produce pairs after travelling distances larger than ∼ 100 Mpc. For inverse Compton scattering, the typical distance TeV electrons travel before they undergo interactions is ∼ 30 kpc. In Figure 3 we show the inverse of the mean free path for these processes, as obtained from Eq. 10.
Higher-order processes to pair production and inverse Compton scattering are important for the propagation of gamma rays of E 10 15 eV. In particular, for scenarios wherein UHECRs induce electromagnetic cascades (see Sec. 3.1), they are an essential ingredient to understand gamma-ray production. The higher-order equivalent of the Breit-Wheeler pair production is the double pair production [267,268] (γ + γ bg → e + + e − + e + + e − ). This process has been extensively studied in various astrophysical contexts, including the propagation of high-energy photons [269,270]. Inverse Compton  scattering can also occur as a second-order process called triplet pair production (e ± + γ bg → e ± + e + + e − ). Its role in the propagation of high-energy photons has long been recognised [271][272][273][274]. This process starts to become important at E 10 17 eV, for cosmological distances. The (inverse) mean free paths for double and triplet pair production are also shown in Fig. 3.
The cosmological propagation of any particle is subject to adiabatic energy losses due to the expansion of the Universe. The change in redshift (dz) for a small propagated distance (dx) can be written as is the Hubble parameter, which in a flat Lambda cold dark matter (ΛCDM) universe is given by with H 0 denoting the Hubble "constant", i.e., the value of the Hubble parameter at present time. Here, Ω m and Ω Λ parameters represent the fraction of the total energy density of the Universe corresponding to matter and dark energy, respectively. According to recent measurements, H 0 ≈ 67.37 km s −1 Mpc −1 , Ω m ≈ 0.3147, and Ω Λ ≈ 0.6853 [275]. Generally, in the presence of a magnetic field the charged component of the electromagnetic cascade (electrons) loses energy through synchrotron emission. However, synchrotron losses are small for intergalactic gamma-ray propagation, since B 10 −9 G (see Figure 1).
The last theoretical ingredient missing for understanding how electromagnetic cascades propagate is the interaction of its charged component with magnetic fields. The equation of motion for a particle of charge q with velocity v in a magnetic field B can be written as where p is the particle momentum. As a consequence of this equation, the electrons and positrons will deflect away from each other. This deflection consists of a circular/helical movement (around the magnetic-field lines), characterized by the Larmor radius r L , given by where p is the absolute value of the particle momentum.
We can now draw a general picture of how gamma rays can be used to constrain IGMFs. Consider an object located at a distance D from Earth, corresponding to a redshift z src , emitting a jet of high-energy gamma rays with an opening angle Θ jet , as sketched in Figure 4. Let Θ los denote the angle between the jet axis and the line of sight, i.e., the angle of misalignment. The primary gamma rays are generated at a redshift z src and produce pairs at z PP , travelling for a distance dictated by the mean free path for pair production (λ PP ) for the energy and redshift of interest (see Eqs. 10 and 11). The pairs produced are deflected by intervening IGMFs, forming an angle δ with the direction of the parent gamma ray. The distance the electrons travel is typically of the order of the mean free path for inverse Compton scattering (λ IC ; see Eqs. 10 and 12). The up-scattered gamma rays can restart the cascade depending on their energy, such that the cascade would have multiple generations of particles. In Figure 4 only one generation is shown. Finally, the secondary gamma rays are detected at Earth forming an angle θ obs with respect to the line of sight, i.e., the line connecting the observer and the object. With this scheme in mind, we can estimate the relevant gamma-ray observables, namely the spectrum, arrival directions, and light curves, either analytically (see Sec. 3.3) or numerically (see Sec. 3.6).
The secondary photons resulting from the electrons deflected in the presence of IGMFs will be delayed compared to primary gamma rays emitted at the same time. Depending on the distance to the source, the duration of the emission and the properties of the IGMF, some of the secondary gamma rays from the cascade, produced from electrons via inverse Compton scattering, will not be able to arrive at Earth within one Hubble time, leading to an energy-dependent decrease in the flux. Therefore, the three main gamma-ray observables relevant for IGMF studies are • spectral effects; • angular distribution; • time delays.
Naturally, quite often combinations of these strategies are employed, as will be presented in more detail in Sec. 4.

Analytical Description of Propagation and Observables
Neronov & Semikoz [161] presented a pedagogical model describing how gamma-ray telescopes can be used to probe IGMFs. This model is a suitable approximation for the energy range of interest, between GeV and tens of TeV. Making use of some simplifying assumptions, they derived analytical expressions for the expected signatures of specific combinations of magnetic-field strength (B) and coherence length (L B ). It is beyond the scope of this review to derive the formulae, but it is certainly worth transcribing the main results and some of the steps required to obtain them.
One can distinguish two regimes of propagation for the charged component of the electromagnetic cascades. They are determined by an interplay between the characteristic scales of inverse Compton scattering (λ IC ) and the coherence length of the magnetic field (L B ). For λ IC L B , the propagation is quasi-rectilinear (ballistic), whereas for λ IC L B , the electrons diffuse before they produce the secondary photons via IC scattering. In the former case, the electrons can be seen as effectively moving in a homogeneous magnetic field, such that in the small-angle approximation δ L B /r L , wherein r L is given by Equation 15. In the latter case, we have δ √ λ IC L B /r L . Together with equations 12 and 15 we then obtain the estimate where E e is the electron energy at redshift z PP . Note that a more detailed investigation [276] shows that the deflection angle also weakly depends on the spectral index α B of the magnetic field (see Sec. 2.1) for L B λ IC . For distant sources, the pairs are produced closer to the source than to Earth (λ PP D). If δ 1, then we can adopt the approximation z z src z PP , which allows us to derive an analytic expression for θ obs : where τ θ is the ratio between the angular diameter distance from the observer to the source and the mean free path for pair production, λ PP . Morphologically, this corresponds to a "halo" of secondary photons around the point-like source. Note that, while the morphology of the arrival directions does resemble a halo in the axi-symmetric case, this is not always the case. Depending on the geometry of the jet (Θ los > 0 • ; see Fig. 4) and properties of the intervening magnetic fields (e.g., helical fields), more complex shapes arise. We continue to use the term 'halo' nonetheless. An interesting and somewhat more accurate approach to estimate the size of such haloes was presented in Ref. [277], in which the moments of the halo distribution are calculated from diffusion-cascade equations. This method is applicable whenever the distribution of gamma rays emitted by the source is isotropic or the jet opening angle (Θ jet ) is sufficiently large.
Another important quantity when determining IGMFs from electromagnetic cascades is the time delay ∆t B , defined as the difference between the following two quantities: the cumulative propagation time of the "reprocessed" gamma rays resulting from the cascades (see Fig. 4), consisting of the lifetime t PP of the primary gamma ray until it results in pair production and of the duration t sec of the cascade from the secondary gamma rays; and the light-travel time (t prim ) of primary gamma rays. Therefore, one can write the equation For the standard consideration of IGMFs we have z PP z src = z 1 and δ 1, such that Eq. 18 becomes The last observable we describe here concerns the probing of magnetic helicity of IGMFs using gamma rays and was first suggested in [278]. Since then, it has been further extended and investigated in a significant number of publications [279][280][281][282][283][284][285][286][287][288]. There it was shown that the helical part of the magnetic field spectrum (see Eq. 5) has a direct impact on the morphology of the halo around the gamma-ray source. In particular, when the magnetic field is helical, the halo becomes "twisted", i.e. instead of an (elongated) circular or oval halo, as one would expect from considering the simple analytic formulas derived above, the result is a spiral-like pattern (see Fig. 8).
This twisted pattern or, more specifically, its handedness, can be measured by the quantity Q introduced in [278] (and summarised in [59]) as Q(n 1 ,n 2 ,x los ) = (n 1 ×n 2 ) ·x los , wheren 1 andn 1 are the unit vectors of the arrival directions of two particles with the respective energies E 1 and E 2 (with E 1 < E 2 ), andx los is the unit vector along the line of sight from the observer to the source. Using this one can calculate the so-called Q-statistics, given by Q(θ max obs ) = Q(n 1 ,n 1 ,x los ) θ obs ≤θ max obs , i.e. the average over all photons with angles θ obs up to a value of θ max obs . If the direction of the line of sight is not known, the arrival directionn 3 of a third particle with an energy E 3 (with E 3 > E 2 > E 1 ) may be considered instead ofx los . In fact, by generally considering such triplets of particles from any direction in the sky, one can calculate the generalised Q-quantity (and, subsequently, the corresponding statistics) as [280] Q(E 1 , E 2 , E 3 , θ max ) = 1 where for every particle with the arrival directionn 3 (and given energy E 3 ) the second summation is carried out over all particles with the given energies E 1 and E 2 (with E 3 > E 2 > E 1 ) with arrival directionŝ n 1 andn 2 , respectively, which lie inside "patches" of angular size θ max aroundn 3 . Finally, the values N 1 , N 2 , and N 3 in Eq. 22 are the corresponding total numbers of particles for each of the three energies. The final step in connecting the Q-statistics (and therefore the handedness of particle arrival directions) with the handedness of the magnetic field (and therefore its helicity) is to consider the case θ max → π/2. As shown in [280], this is proportional to the helical part of the spectrum H k , as defined in Eq. 5.
An alternative to the Q-statistics, introduced in [284], is the S-statistics which, for a single source, can be used to quantify the spiral shape of the halo.

Plasma Instabilities
The physics of electromagnetic cascades described above is well understood, but it neglects the back-reaction of the intergalactic medium on the cascades. This is a common assumption adopted in most IGMF studies, but if it turns out to be a poor approximation, plasma effects may become dominant. It was suggested [289] that the electrons in the cascade interact with the IGM and lead to the generation of plasma instabilities, losing their energy and consequently heating the IGM. Due to the extreme parameters of the interacting components (for example, a factor of up to 10 24 between the density of the electron beam and the background plasma [290]), it is practically impossible to exactly calculate the impact of the instabilities on the development of the cascade. Nevertheless, one can rely on approximations and/or extrapolations.
The IGM parameters relevant for plasma instabilities are its temperature, which is typically T IGM ∼ 10 4 K [289], and the density, which in the cosmic voids is n IGM ∼ 0.1 m −3 [291]. Another important parameter is density of the gamma-ray beam, which is related to its luminosity.
As mentioned above, there is no general agreement on whether plasma instabilities are important for the propagation of electromagnetic cascades. Even if one accepts this assumption, it is not clear which kind of instability could be dominant. In fact, the modulation [290,[292][293][294], oblique [289,295,296], kinetic [297], and longitudinal [298] instabilities, as well as non-linear Landau Damping [299] have been considered in the literature. On the other hand, Ref. [300] found that even if they are present, the effect of plasma instabilities is too small to cause a significant impact on observations. A comparison of the energy-loss length for different types of instabilities is shown in Fig. 5.
Several authors subsequently published results of actual simulations of gamma-ray propagation including possible plasma instabilities effects and compared them to actual observations [301][302][303]. The results show that, while the instabilities can, indeed, lead to appreciable deviations from the paradigmatic picture of cascade development, they may not be sufficient to render gamma-ray constraints of IGMFs completely ineffective. In this case, all that would be required is a more detailed modelling of the electromagnetic cascades -which, understandably, would be more susceptible to uncertainties due to the inclusion of an additional and poorly-understood effect.
There is also another window of opportunity to evade plasma instabilities, even if they majorly disrupt electromagnetic cascades. The growth rates of plasma instabilities are often estimated using simplifying assumptions like a continuous and constant stream of particles. However, if the object in question emits gamma rays in flares, the temporal structure of the resulting charged beam should be  Cooling rates due to plasma instabilities computed at z = 0, according to different models [289,292,294,295,299]. This example is for a typical scenario for an IGM density of n IGM = 10 −1 m −3 and temperature of 10 4 K, for a blazar beam with luminosity L = 10 38 J/s. We also present the inverse mean free path for inverse Compton scattering in the CMB for comparison.
considered. In particular, if the duration of the flare is short enough, the instability might not have enough time to fully develop, consequently having no significant impact on the electrons.

Other Propagation Phenomena
So far we have discussed how electromagnetic cascades propagate in the Universe in light of a standard picture entirely contained within the framework of quantum electrodynamics (Sec. 3.2). We also briefly discussed how plasma instabilities could quench electromagnetic cascades propagating in the IGM (Sec. 3.4). In this subsection, we briefly describe how other physical phenomena could affect the development of electromagnetic cascades and, consequently, observations of high-energy gamma-ray sources, with direct implications for IGMF studies.
One phenomenon that can interfere with the propagation of gamma rays and consequently compromise IGMF constraints is gravitational lensing. Massive objects can significantly deform the space-time surrounding them, altering the path along which particles travel. As a result, in the context of gamma rays, gravitational lenses can significantly deform the morphology of haloes and, in the case of flaring objects, increase the time delays due to this gravitationally-induced contribution. The first source for which gravitational lensing has been observed at gamma-ray energies (up to 30 GeV) was PKS 1830-211 [304]. Since then, the phenomenon has been detected for this and other gamma-ray-emitting objects [305][306][307][308]. Ref. [309] investigated how macrolenses could compromise estimates of the optical depth for pair production, concluding that this effect would not lead to any measurable changes in this observable. This results is corroborated by the more detailed study of [310].
Other potentially important phenomena arise in the context of BSM models. The most widely studied BSM processes that could interfere with the gamma-ray-IGMF framework we present here involve Lorentz invariance violation (LIV) and interactions with axion-like particles (ALPs). Because of the potentially important role played by these phenomena in determining the gamma-ray signatures of sources used for IGMF constraints, we briefly touch upon these issues. Lorentz invariance violation is a possible consequence of various BSM approaches, especially in the context of quantum gravity. The standard approach, from the field theory side, is to create a minimal SM extension, in particular by introducing additional terms to the SM Lagrangian, resulting in a effective field theory with LIV [311].
In terms dynamics, the main effect when it comes to the propagation of particles is the modification of the dispersion relation, given by where E LIV and E is the particle energy with and without LIV, respectively, η is a dimensionless parameter measuring the strength of the LIV, and M Pl is the Planck mass. This, on the one hand, changes the threshold of a given reaction, and, as a consequence, also changes the corresponding propagation length, as it modifies the limits of the integral in Eq. 10. In addition, new reactions, which are not possible without LIV, may then be kinematically allowed, such as, to name the ones most relevant in the context of this review, spontaneous photon decay into pairs/photons, the vacuum Cherenkov effect for electrons and charged UHECRs, as well as spontaneous photodisintegration of multi-nucleon nuclei. For an overview on the modifications of the processes in the electromagnetic cascades and in UHECR propagation, see [312][313][314] and [312,315,316], respectively. All this may result in a significant modification of particle propagation and, therefore, impact the corresponding observations. In particular, [317] showed that LIV might dramatically increase the interaction length of pair production for energies above 100 TeV and therefore suppress the cascade development. On the other hand, LIV might also imply that the photons' speed is energy-dependent, thus resulting in energy-dependent time delays [314].
Axion-like particles, or ALPs, appear in extensions of the SM. They are pseudo Nambu-Goldstone bosons associated with a broken symmetry U(1). They were originally introduced by Peccei and Quinn [318,319] as a solution to the strong CP problem. ALPs couple to standard-model particles via the Lagrangian [320] where E γ denotes the electric field of the photon itself, and B ext represents an external magnetic field (in the context of this review, IGMFs). The coupling constant g aγ determines how strongly photons, in our case gamma rays, will interact with the ALP field. For a distance x, the probability of a photon to convert into an ALP (and vice-versa) is Here the ∆ terms refer to the solution of the equation of motion derived from the Lagrangian (Eq. 24). They describe: the coupling between photons and ALPs (∆ aγ = g aγ B T ) between a photon of energy E and the ALP field for a an external magnetic field B T transverse to the direction of propagation of the photon; the kinetic term (∆ a = m 2 a /2E) for an ALP of mass m a ; the polarisation states of the photon (∆ and ∆ ⊥ ), which, in our case, encompass the contribution of the IGM plasma (∆ pl = −ω pl /2E) for a plasma frequency ω pl , and the QED vacuum polarisation (∆ QED ∝ B 2 T ), which depends on the direction (∆ QED,⊥ = 7∆ QED /2 and ∆ QED, = 2∆ QED ). For more details, the reader is referred to, for example, Ref. [320].
Two regimes of propagation can be identified [321], depending on whether the gamma-ray energy is larger or smaller than the critical energy (E c ), given by: The limit E E c corresponds to the so-called strong mixing. In this case, the probability of conversion (see Eq. 25) does not depend on the energy. If E E c , the energy dependence becomes salient leading to an effective low-energy cut-off.
The propagation of gamma rays will be affected by ALPs in multiple ways. Firstly, the magnetic fields in the sources will contribute to the total ALP-photon mixing. Secondly, once the gamma rays are injected into the intergalactic space, they may initiate electromagnetic cascades as described in Sec. 3.2. Upon entering the Galaxy, ALP-photon mixing may also occur due to the Galactic magnetic field. The oscillation probability from Eq. 25 will then be a combination of the probabilities in each of these environments, as discussed in, e.g., Refs. [322][323][324][325].
In the case of gamma rays propagating over cosmological distances, the oscillation probability from Equation 25 implies deviations from the expected transparency of the Universe, since gamma rays will be able to travel longer without undergoing pair production. A number of works investigated the possibility that this "pair production anomaly" could be related to ALPs (e.g. [324,[326][327][328][329]).
The effects of IGMFs on gamma-ray-ALP interconversion has been studied adopting several methods ranging from semi-analytical approaches to more detailed simulations. While the first studies on the topic assumed relatively simple magnetic-field configurations, later studies [330,331] improved the treatment including turbulent fields (see Sec. 2.1 for details). Investigations considering the actual distribution of magnetic fields in the magnetised cosmic web have also been performed [332] While ALPs are an important ingredient that could play a leading role on the intergalactic propagation of gamma rays in a magnetised Universe, they have not been observed and only constraints exist. Some limits were derived using gamma-ray observations [328,[333][334][335], but much of the parameter space is excluded due to observations of photons in wavelengths other than gamma rays. Interestingly, some works have been obtaining combined limits on IGMFs and ALPs together [336]. For reviews on the status of the field, see, e.g., Refs. [337][338][339].

Propagation Codes
The simulation of the propagation of electromagnetic cascades in the intergalactic medium is often done numerically, employing some approximations to enable (semi-)analytical solutions (e.g., [48,340,341] Elmag [342,345] is a Fortran code that tracks the development of electromagnetic cascades. In the first two versions of the code [342], the effects of magnetic fields on the charged cascade component, namely time delays and deflections, were taken into account using the small-angle approximation. Therefore, this version was limited to low magnetic-field strengths. The newest version, Elmag 3.01 [345], adds a Lorentz-force solver that enables three-dimensional simulations assuming turbulent magnetic fields generated based to Eqs. 5 and 7, following Refs. [346,347], as well as custom grids. CRPropa [343,348,349] is a well-known code for ultra-high-energy cosmic-ray propagation, written in C++ and with Python bindings (since version 3). The original CRPropa [348] and CRPropa 2 [349] made use of the numerical methods from Ref. [340], namely using transport equations to treat the development of electromagnetic cascades. The newest version include a full treatment of electromagnetic cascades [286,343,350]. A variety of magnetic-field configurations are available, and the code is flexible enough to handle customisations and arbitrary magnetic-field grids. Moreover, it can generate turbulent magnetic fields on the grid or using grid-less methods [351], with improvements from [352]. Earlier releases of CRPropa 3 supported the propagation of gamma rays with energies 10 17 eV through the Monte Carlo EleCa code [353]. Due to the computational limitations, EleCa is restricted to the highest energies, but in CRPropa a hybrid approach using the transport-equation treatment of [340] was available. Recent developments enable a full Monte Carlo treatment of photons from ultra-high down to GeV energies, which is useful for exploring the UHECR-induced cascade scenarios (see Sec. 3.1).
A Fortran code for cascade propagation was developed by Fitoussi et al. [344]. It does not rely on any approximations, performing the full three-dimensional propagation of the cascades. In this code, the magnetic field is composed of cells with randomly oriented strengths. A semi-analytical treatment of the cascade development in Mathematica is implemented in γ-Cascade [354].
Plasma instabilities are often neglected in simulations, or treated within a dedicated MHD computational framework. In [303], grplinst was presented. It is a module for the CRPropa code that implements plasma effects on the electrons as an additional energy-loss term of the form where x is the length of the trajectory described by an electron (or positron) of energy E e , z is the redshift, and τ is the electron energy-loss time due to the plasma instability. Within this simplified treatment, the time scale in which the instability grows (T ) is taken to be the electron cooling time (τ). Therefore, Eq. 27 is overestimated, since in reality T ≤ τ. More recently, these same parametrisations of grplinst [303] were implemented in Elmag 3.02.

Examples
To illustrate the effects of IGMFs, we present the results of Monte Carlo simulations of the development of electromagnetic cascades in the IGM. To this end, we use the CRPropa code [343], but similar results could have been derived with, e.g., Elmag [342,345] or the code presented in Ref. [344].
We first select the archetypical blazar 1ES 0229+200 [202], used in several IGMF studies (cf. Sec. 3.1). This object is located at a distance corresponding to z 0.14. We fix the coherence length to L B = 1 Mpc to illustrate the formation of the haloes around the source. This is shown in Fig. 6. Note that these plots are shown in the coordinate system of the simulation, as we would observe from Earth, but they can be immediately converted to another coordinate system, such as Galactic or equatorial.
It is evident from Fig. 6 that a significant fraction of the flux is not contained within a finite-sized containment radius centred at the source. This causes spectral changes with respect to the point-like source flux, as shown in Fig. 7.
If magnetic fields are maximally helical, then the halo shape shown in Fig. 6 changes considerably. In fact, the changes can be so drastic that the morphology of the arrival direction pattern is no longer a standard axi-symmetric halo. For a source pointing straight at Earth (Θ los ), we expect a spiral-like pattern, as shown in Fig. 8 for a hypothetical source at z = 0.08. In this case, the handedness of the halo reflects the sign of the helicity: left-handed for H B > 0, and right-handed for H B < 0.
The arrival directions of gamma rays can be quantified through the calculation of the Q-factors (see Eq. 22). The effects of the helicity of IGMFs are more pronounced for large coherence lengths, hence the choice of L B = 250 Mpc in the example of Fig. 8. The smaller the ratio between the source distance and the coherence length, the more diluted the signal is, which would reflect in the Q-factors (see [284,288]).

Results
Fluxes of distant objects, like the ones used in IGMF studies, are normally computed for a point-like source. The magnetically-induced broadening of the electromagnetic cascade will naturally affect the measured point-like flux, especially at lower (typically E 10 GeV) energies, since these are predominantly secondary gamma rays if the intrinsic spectrum extends beyond ∼ TeV energies. In this case, the larger the angular broadening caused by IGMFs, the more pronounced is the suppression of the gamma-ray flux from a point-like source, since a fraction of the events will leak outside the point spread function (PSF) of the detector.
Gamma-ray sources are observed during a given time. If IGMFs are such that the incurred time delays (T) exceed this window of observation, then the measured flux will be affected. The suppression will, in general, be stronger as energy decreases because this contribution is likely produced by lower-energy electrons whose Larmor radii increase with energy (see Eq. 15). The relevance of this effect depends on an interplay between the duration of the emission, which depends on the type of object (see Section 3.1), the time window of observation, and the magnetically-induced time delay.
A number of studies attempted to constrain IGMFs based on gamma rays using different methods. One possible way to classify these studies is the number of sources used, such that in Sec. 4.1 we describe the results of analyses of individual gamma-ray sources, and in Sec. 4.2 those of multiple stacked sources.  [202]. The source parameters are the same as in Figure 6.
In general, the results are for the magnetic-field strength. However, there have also been attempts to constrain the coherence length and helicity of IGMFs with gamma rays, which we discuss in Sec. 4.3 and 4.4, respectively, followed by results for IGMFs considering that the cascades are induced by UHECRs in Sec. 4.5. Finally, in Sec. 4.6 we discuss the prospects for IGMF measurements with gamma-ray observatories.

Analyses of Individual Sources
The first constraints on IGMFs using gamma rays were derived by Neronov Fig. 9. This dependence of the lower limit of IGMFs on the coherence length, L B , follows from the simplified approach commonly used (see Sec. 3.3), and is adopted in most of the works to which we refer below, with a few exceptions. Most importantly, this work was the first to firmly exclude the case B = 0. An earlier investigation [355] of the blazar 1ES 1101-232 concluded that an exceedingly hard intrinsic spectrum for this object would be required to account for observations, unless the EBL was more intense and the IGMF were stronger (B 10 −15 G). Ref. [356] argued along the same lines, when interpreting observations of the blazar H1426+428. . The source, assumed to be located at z = 0.08, emits gamma rays with a spectrum E −1.5 and an exponential cutoff at E max = 100 TeV (see Eq. 9). Its jet has an opening angle Θ jet = 5 • and it points directly at Earth (Θ los = 0 • ; see Fig. 4). The colour scale indicates the normalised spectrum-weighted number of detected events in the angular bin.
the authors perform additional checks about the energy range of the Fermi-LAT data used in the analysis, demonstrating that the results are the same regardless of whether a dataset containing gamma rays with energies starting from 100 MeV or 1 GeV is used. The significance of these results decreases slightly for other EBL models (Refs. [259,357]). Note that a more detailed treatment of the cascade interactions would increase the flux at lower energies, so that these estimates are actually conservative.  Fig. 9. In principle, because the coherence length was assumed to be L B = 1 Mpc in both works, one could be tempted to extrapolate the conclusions to L B > 1 Mpc. However, because the objects used to derive the constraints (1ES 1218+304 and PKS 2155-304) show some intrinsic variability [358][359][360], care should be taken extrapolating the bounds to larger values of the coherence length.
The Fermi-LAT Collaboration [47] compiled a catalogue of sources that was used to constrain IGMFs and performed a detailed analysis of this sample of blazars. They found no evidence of extended emission, neither around individual objects, nor in the stacked analysis. This way, they could constrain the allowed values for the strength of IGMFs: B 10 −16.5 G for L B 10 kpc. This result is conservative, assuming T 10 yr. If this condition is relaxed, the bounds are even stronger: B 10 −14 G and B 10 −12.5 G, for T 10 4 yr and T 10 7 yr, respectively, as shown in Fig. 9. While these limits were derived for a jet with half-opening angle Θ jet ≤ 10 • (see Fig. 4), no misalignment was considered (Θ los = 0 • ). These results were derived for a combination of sources which include, among others, 1ES 0229+200 and 1ES 1218+304. Because there are indications that they could be variable [359,361], if these sources are removed from the analysis, limits which are more conservative are obtained. for the EBL models labelled D11 [261] and the lower-limit S16l model [263]. Note that the regions plotted refer exclusively to the region of the parameter space reported in the corresponding references, without extrapolations to high/low values of the coherence length. The grey region are the combined excluded regions from Fig. 1, obtained via other methods.
The time scale over which a given source emits gamma rays (T) influences the bounds one can derive. For instance, Ref.
[43] analysed VERITAS and Fermi-LAT data. The lower bounds they obtained for L B 1 Mpc were B 3 × 10 −18 G and B 2 × 10 −16 G, for source activity periods T 3 yr and T → ∞, respectively. The former is evidently more conservative, as it encompasses exclusively the period for which there are observations of the object (RBG J0710+591). Similar considerations about T were discussed in [341], which obtained B 10 −18 G for L B 1 Mpc, for 1ES 0229+200 (see Fig. 9). These results are order-of-magnitude compatible with the lower limits by Ref. [362], which are B 10 −16 -10 −18 G.
The flaring object Mrk 501 drew much interest for IGMF constraints, due to its variability (see, e.g., [21,23,47,[363][364][365][366][367][368]). The prospect for detecting pair echoes from this same object was studied in [369]. With an analysis of observations of its 2009 flare by MAGIC, VERITAS, and Fermi-LAT, the authors of [366] argued that its spectrum and time profile could be explained by IGMFs with B 10 −17 -10 −16 G (for L B 1 Mpc). Ref. [367] studied the pair echoes from this same blazar, concluding that B 10 −20 G, assuming L B 1 kpc, at a 90% confidence level. Making use of similar methods, Ref. [370] analysed data from ARGO-YBJ and Fermi-LAT for Mrk 421, excluding B 10 −20.5 G for L B 1 kpc, at a 4σ level. The results of the latter analysis are particularly interesting because they do not make assumptions about the intrinsic spectrum of the source during periods when it is not observed.
A somewhat elaborate treatment of the cascade development was adopted in [44], in which Monte Carlo simulations were used to derive bounds on IGMFs for a sample of three blazars (1ES 0229+200, RBG J0710+591, and 1ES 1218+304). The limits obtained depend on the strategy adopted for the analysis: B 10 −15 G considering the absence of haloes, as observed by Fermi-LAT, and B 10 −17 G considering time delays.
An important factor that significantly affects IGMF estimates is the model of the EBL. Ref.
[45] studied this dependence for the archetypical extreme blazar 1ES 0229+200. They found B 10 −17 G, which can increase by nearly two orders of magnitude depending on the EBL model. In fact, the EBL is one of the main intrinsic uncertainties that hinders the exclusion of the scenario B = 0, as argued in [371]. In this analysis, among the seven blazars considered, only one led to B > 0 irrespective of the choice of EBL model. However, the uncertainties in the intrinsic source spectrum and EBL model might be unrealistic, as noted in refs. [58,59].
The analysis by Dolag et al. [46] is interesting because it employed magnetic fields obtained from cosmological magnetohydrodynamical simulations from [196]. These simulations are constrained, i.e., they roughly reproduce the distribution of large-scale structures up to hundreds of Mpc. At larger distances, this cosmological volume was replicated up to the distance of 1ES 0229+200. The authors showed that more than ∼ 60% of the Universe along the line of sight of this object have magnetic fields with strength B 10 −16 G. Interestingly, this analysis also showed that haloes can be used to probe the maximal energy of the gamma rays emitted by a source, E max (see Equation 9). In fact, there is a considerable correlation between the value of E max that could be inferred from fits in the presence and in the absence of IGMFs [241].
Ref. [372] employed a Monte Carlo code to model the development of electromagnetic cascades initiated by GRBs. However, because GRBs had not been detected at TeV energies until recently, when MAGIC observed GRB 190114C [227], the authors extrapolated the GeV flux of GRB 130427A measured by Fermi-LAT up to TeV energies. If the extrapolation is correct, the lower limit obtained is B ∼ 10 −17.5 G (for L B 1 Mpc).
With the first observation of VHE emission by a GRB, it was possible to effectively constrain IGMFs with gamma-ray observations. By combining Fermi-LAT and MAGIC data, Ref. [229] obtained a lower limit of B 10 −19.5 G for L B 100 kpc. Ref. [228] performed a similar analysis for GRB 190114C using Monte Carlo simulations. They concluded that Fermi-LAT is not sensitive enough to detect the cascade signal from this GRB on time scales of one month. The discrepancy between these two works is due to a combination of factors. Firstly, the former employed a simpler semi-analytical method, whereas the latter performed detailed Monte Carlo simulations using the three-dimensional version of the Elmag code. Moreover, the authors of [228] reconstructed the intrinsic spectrum following [255], while in Ref. [229] a fixed spectral index α = 2 (for a power-law distribution ∝ E −α ) was assumed. Yet another difference is the treatment of the time information of the photons. While [228] accounted only for photons detected more than 62 s after the burst, [229] adopted 6 s. This issue is far from simple, as it requires knowledge of the inner workings of gamma-ray bursts. For further details, the interested reader is referred to the works on GRB 190114C by the MAGIC Collaboration [227,373].

Stacked and Diffuse Analyses
It is rather difficult to observe magnetically-induced haloes from individual sources, as they are normally not bright enough to be detected [37, [374][375][376]. Hence, techniques that are more sensitive are needed. Analyses of stacked samples of blazars could be useful for this purpose, since the signal-to-background ratio increases, easing the identification of any excess to the detector's PSF.
The authors of [377] performed a stacked analysis of 170 AGNs using 11 months of Fermi-LAT data. They claim to have found an excess to the PSF of the detector 0.5 − 0.8 • , at a 3.5σ level. Haloes of these sizes are caused by B 10 −15 G (see Eq. 17). Nevertheless, it was later shown that these results could be attributed to instrumental effects associated to different treatments of photons measured in different parts of the detector [378] (see also Ref. [379]). This was not included in the PSF used in Ref. [377], thus leading to an incorrect estimate of the strength of IGMFs.
The stacked analysis from Ref. [380] considered 24 selected high-synchrotron-peaked BL Lacs, at z < 0.5. Using Fermi-LAT data, the authors found indications of extended emission, consistent with B ∼ 10 −17 -10 −15 G. However, an updated analysis using 12 objects of the same population resulted in no compelling evidence for an extended emission, with only a modest 2σ significance for B ∼ 10 −15 G [381].
Another stacked analysis of a sample of 394 AGNs, 158 of which present flaring activity, was performed in Ref. [382]. Interestingly, the method employed considers temporal information of the sources by comparing the fluxes during low quiescent states and during flaring periods. No evidence for pair haloes was found. The recent analysis by Fermi-LAT [47] corroborates this result, finding no indications for extended emission in the stacked source samples of high-synchrotron peaked BL Lacs.
Using the method introduced in [383], Ref. [376] identified misaligned blazars from a catalogue of radio-loud AGNs, and performed a search for pair haloes around the stacked sample of these objects. They showed that a magnetic field with BL 1/2 B 10 −15 G Mpc 1/2 would lead to specific halo anisotropy patterns that are not observed, thus providing an upper limit on the strength of IGMFs. Note, however, that the assumptions about the intrinsic properties of the considered sources are subject to uncertainties. Considering the available lower limits derived in the works discussed above, the parameter space that would remain available for IGMFs would be tiny or, considering the stronger constrains from the recent results by Fermi-LAT [47], inexistent. The authors of [376] then conclude that, if there is indeed no room for IGMFs that can explain the observations, then some other process might be at play that quenches the electromagnetic cascades. They claim that this could be due to, for instance, plasma instabilities (see Section 3.4). More recently, the same group claimed to have found convincing evidence of the non-existence of pair haloes. Using the same method, they exclude B ∼ 10 −16 -10 −15 G with L B > 100 Mpc at a 3.9σ level, and B ∼ 10 −17 -10 −14 G at 2σ [49] (see Fig. 9).
An interesting idea to constrain IGMFs is to study their possible imprints on the diffuse gamma-ray background (DGRB) [384,385], even though the validity of the assumptions used in these analyses is unclear. The presence of IGMFs may suppress the lower-energy diffuse gamma rays measured. Interestingly, the authors of Ref. [384] claim that, in fact, the observations by Fermi-LAT already disfavours the scenario with null IGMF. This agrees with [385], who found that the contribution of cascade gamma rays from blazars to the DGRB changes significantly in the presence of IGMFs, such that for B 10 −12 G the blazar contribution to the spectrum of the DGRB changes.
In the context of diffuse searches, it is important to keep in mind that, in addition to the uncertainties in the EBL models (see Fig. 2), there are fluctuations correlated to the processes that produce the EBL photons. This leads to inhomogeneities in the EBL distribution that can affect the propagation of electromagnetic cascades. However, as shown in [386], this effect is small ( 1%), so it should have little impact on IGMF measurements using diffuse gamma-ray observations.

Bounds on the Coherence Length
A method to measure the coherence length was suggested in Ref. [387]. In this case, the slope of the light curve of secondary gamma rays would provide an upper limit on L B . More specifically, the time dependence of the flux would be ∝ 1/ √ ∆t B for coherence lengths much larger than the mean free path for inverse Compton scattering (L B λ IC ), and approximately constant if L B λ IC . Similarly, the angular profile of the haloes can also retain information about the coherence length. For L B λ IC , the surface brightness profile is roughly uniform, whereas for L B λ IC it decays as the angular distance to the centre of the source increases.
With the first multimessenger observations of high-energy neutrinos from the flaring blazar TXS 0506+056 in coincide with electromagnetic radiation [388,389], Ref.
[50] used the cascade signal delayed with respect to the neutrino emission to constrain IGMFs. The derived limits depend on the EBL model, such that the hypothesis of null IGMFs could only be rejected for two out of the four models tested (Domínguez et al. [261] and the lower-limit model by Stecker et al. [263]). Interestingly, while the bounds are not robust, this work derived, within the investigated parameter space , limits on the coherence length of IGMFs for the first time: 30 kpc L B 300 Mpc, at a 90% C.L., shown in Fig. 9. Naturally, the significance of this result depends on the reliability of the neutrino-gamma-ray correlation and on the assumptions made, namely that the IGMF has a Kolmogorov power spectrum, and that the intrinsic spectrum of TXS 0506+056 both during the flaring and quiescent periods can be described by a power-law with an exponential cut-off.

Constraints on the Magnetic Helicity
There has been a growing interest in probing the helicity of IGMFs, given its importance for understanding magnetogenesis. All-sky analyses of Fermi-LAT data employing the parity-odd correlators described in Section 3 found indications of IGMFs with B ∼ 10 −14 G at L B ∼ 10 Mpc and an overall negative (left-handed) helicity [279,281]. More recently, a re-analysis of a larger data set showed this result to be a fluctuation stemming from a miscalculation of the statistical significance that neglected the look-elsewhere effect [288]. The same publication also claims that it is currently challenging to detect helicity, both in the fluxes of individual sources and in the diffuse gamma-ray background. In addition, the authors of [287] claim that they did not find any handedness using the Q-statistics for Fermi data, being, however, unable to state definitively whether there is actually no handedness present or whether the Q-statistics is not sensitive enough for measuring it.
The signatures of helical IGMFs on the shape of haloes are unique, with significant deviations from axial symmetry, as illustrated in Fig. 8. Moreover, the sign of the helicity directly correlates with the handedness of the morphology of the arrival directions of gamma rays. In Ref. [282], the authors employed a semi-analytical method to show that spiral-like patterns are the natural shape of the arriving gamma rays for helical fields. Nevertheless, within this simple framework, IGMFs were assumed to be homogeneous, which is not realistic unless the coherence lengths involved are exceedingly high, comparable to the distance of the source. In the more realistic case of turbulent magnetic fields with coherence lengths possibly shorter than the distance from Earth to the gamma-ray source in question, the spiral pattern could vanish, being diluted into something closer to a typical axisymmetric halo. This was, indeed, observed in a more detailed study using three-dimensional Monte Carlo simulations [284]. This work, however, does support the measurement of the helicity of IGMFs for L B 50 Mpc, for sources at redshifts z 0.10.

Constraints from UHECR-produced Gamma Rays
In Sec.3.1 the model proposed in Refs. [249,252] was presented. In this scenario, the flux of some extreme blazars could be attributed to cosmic-ray interactions along the line of sight. Essey et al. [253] constrained IGMFs considering that the gamma rays observed are a combination of those emitted by the blazars and those stemming from CR interactions. In this case, the combined limits from all three blazars analysed favour 10 −17 B/G 10 −14.5 , at 95% C.L. This result is robust with respect to the choice of EBL model. It is also shown in Fig. 9. Other authors also performed similar investigations (e.g., [256,[390][391][392][393]).
Interestingly, within the UHECR-cascade framework, photons with E 10 TeV could be detected even if the sources are very distant (z 0.1). Nevertheless, for IGMFs with B 10 −14 G, significant deviations from a point-like flux would be expected due to magnetically-induced deflections, compromising any constraints that one could derive in the context of this hadronic model [253].
For blazars, the investigations of the role of line-of-sight interactions in gamma-ray measurements [249,250] were also shown to lead to time variabilities that are characteristic for specific magnetic-field properties, of the order of years for B = 10 −15 G [390]. Nevertheless, the variabilities cannot be too short since even for weak IGMFs of the order of 10 −18 G cascade photons with E = 10 GeV would be magnetically-delayed by ∼ 10 yr [391]. In the purely leptonic scenario, this timescale is shorter by a ten fold [251].
A detailed account of the effect of magnetic fields on both the electromagnetic cascades as well as on their progenitor UHECRs was presented in [256]. Using three-dimensional simulations of the magnetised cosmic web from [196] and detailed numerical methods, the authors found that the cascade broadening could be detected with next-generation gamma-ray telescopes and possibly some of the ones in operation today.
Note that the propagation of cosmic rays is not trivial. There are many uncertainties involved (see, e.g., [394,395] for a discussion), which might compromise the production of gamma rays. Moreover, depending on the location of the blazar in the cosmic web, local magnetic fields (e.g. in filaments) might significantly deflect cosmic rays away from the line of sight [196][197][198][199] (see the discussion at the end of Section 2).

Prospects for Measurements of IGMFs
From the discussion so far a general picture of IGMFs emerges, wherein gamma rays play a fundamental role in excluding part of the parameter space shown in Fig. 1, as summarised in Fig. 9. It is important to bear in mind that there are many factors that could compromise the derived limits shown in the latter figure. This includes uncertainties regarding the intrinsic source spectrum and possible variability, the knowledge of the EBL, the distribution of magnetic fields in the Universe, the contribution of a putative hadronic component to the cascade, etc. Moreover, plasma instabilities may quench electromagnetic cascades, even if this effect is minor. The central question that arises is, therefore, if the next generation of gamma-ray telescopes, whether ground-or space-based (or a combination of both), will be able to unambiguously detect them. In this subsection we briefly revisit the theory that can be directly connected to the experiments. In particular, we highlight here the requirements for next-generation detectors to be sensitive enough to detect IGMFs.
In general, the detection of haloes depends on two factors. First, the size of the extended emission should be such that it is fully contained within the field of view (FoV) of the detector, of size θ fov . Second, this extension must exceed the angular resolution of the detector (θ psf ). In other words, the signal can be observed if the PSF and FoV of the instruments satisfy θ psf < θ obs < θ fov .
Current-generation IACTs (VERITAS, MAGIC, H.E.S.S.) can resolve scales of ∼ 0.08 • , with a typical FoV of 6 • . The upcoming Cherenkov Telescope Array [396] will reach angular resolutions as high as ∼ 0.02 • with a field of view of 20 • , improving the possible constraints on IGMFs [397][398][399][400]. For instance, 50 hours of observations of the blazar 1ES 0229+200 could be used to probe magnetic-field strengths of B ∼ 10 −13.5 G at a 5σ-level, for L B 1 Mpc [401]. In particular, with angular resolutions of 0.13 • at E 100 GeV, a combination of CTA and Fermi-LAT observations could also be used to probe magnetic helicity [284], although it is not clear if any measurable signal could be extracted [288]. Fig. 10 shows the region of the parameter space that can be probed with the halo strategy. Note that, while the PSF must necessarily be smaller than the size of the halo for any extended emission to be identified, the condition θ fov > θ obs is not a strict requirement. Nevertheless, if the halo is not entirely contained within the field of view of the instrument, it becomes difficult (but not impossible) to reconstruct the image, due to uncertainties stemming from the reconstruction procedure and the motion of the telescope to scan the region surrounding the source. Typically, IACTs have higher angular resolutions near the centre of the FoV, decreasing radially from that point. . This example was calculated using the approximation given by Equation 18 [161]. The source is assumed to be located at a distance corresponding to redshift z = 0.10. Fig. 10 was obtained using simplifying assumptions, in particular Equation 17, derived in [161]. If these estimates were improved using detailed Monte Carlo simulations and instrument response functions were accounted for, then the picture could change slightly. Nevertheless, a recent work by the CTA Consortium [401] using simulations obtained with the CRPropa code is in qualitative agreement with Fig. 10.
More generally, IGMF constraints based on gamma-ray observations employing the halo strategy depend on the point-source sensitivity of the instruments, shown in Fig. 11. A simple comparison with the simulations of Fig. 7 demonstrates that instruments like CTA will be able to probe IGMFs stronger than ∼ 10 −14 G, as shown in Ref. [401]. The sensitivities shown in Fig. 11 are a useful guide for a first assessment of the instrumental capabilities in IGMF studies through simple comparisons with theoretical expectations (Fig. 7). Nevertheless, there are multiple conceivable ways to probe IGMFs with halo searches. The simplest one is the direct search for an extended emission, as we discussed in the preceding paragraphs, but one could also employ methods involving the fit of the halo profile and comparison with the background, for example. This would lead to differences in the sensitivity curves, as discussed in detail in Ref. [397] for the case of CTA.
It is worth stressing that facilities operating at slightly higher or lower energies can play an important role in this type of study, despite being seldom considered for IGMF studies. They can be used to constrain putative PeV gamma rays as well as cascade photons in the GeV band. The current and upcoming facilities operating at higher energies, like LHAASO [8] and the planned Southern Wide-field Gamma-ray Observatory (SWGO) [402], formerly known as the Southern Gamma-ray Survey Observatory (SGSO), can help in the precise determination of the intrinsic spectrum of the sources and consequently lead to better models. Observatories such as the planned e-ASTROGAM [403] and the All-sky Medium Energy Gamma-ray Observatory (AMEGO) [404] can detect secondary (cascade) photons in the MeV-GeV band, thus providing additional insights into IGMFs. For the extreme blazars with hard spectra (see discussion in Sec. 3.1), in particular, this will ultimately reduce the uncertainties when constraining IGMFs. Interestingly, observations around GeV energies may also probe spectral features expected from some plasma instability models (e.g. [303]; see also Sec. 3.4).  Figure 11. Comparison of the point-source sensitivity for various gamma-ray observatories. The Fermi-LAT band encompasses sources at various positions in the sky, for the P8R3_SOURCE_V2 instrument response function [405]. The sensitivities for the IACTs, namely VERITAS [406], MAGIC [3], H.E.S.S. [407], and CTA [396] are given for 50 hours of observations. For SWGO [402] and LHAASO [8], the curves shown are for five and one year, respectively. One year is also the observation time used to derive the sensitivity for AMEGO [404]. For HAWC [408], the curve corresponds to 507 days (which is equivalent to approximately 3000 hours) of observations. The thick black lines correspond to the simulations from Fig. 7. Note that the instrument response functions of each detector are not folded into the simulations; the corresponding sensitivities are shown here just for the sake of comparison.
All currently operating instruments can resolve short-duration events from sources at distances closer than z ∼ 1, probing magnetic fields with strengths B 10 −17 G for L B 1 Mpc; note that the exact value of B that can be probed depends on the distance to the source. For stronger magnetic fields, however, it becomes difficult to detect time delays if they are larger than a few years or a decade. In fact, according to Eq. 18, the expected time delay for 10 GeV gamma rays assuming B 10 −17 G and L B 100 kpc would already be ∆t B 100 yr, posing obstacles for measurement within a reasonable time window of a few decades. Fig. 12 shows the region of the parameter space that can be probed with the time-delay strategy. It is clearly favourable for probing the region of the parameter space corresponding to weaker magnetic fields (compared to Fig. 10). This particular example is for a source at redshift z = 0.42, the same as GRB 190114C [227]. . This example was calculated using the approximation given by Equation 18 [161]. The source is assumed to be located at a distance corresponding to redshift z = 0.42.
In a recent work [368], the prospects for measuring strong IGMFs (B 10 −12 G) were analysed using the constrained cosmological simulations of the cosmic web from Ref. [409], based on gamma-ray observations from both Mrk 421 and Mrk 501. The authors argue that, at least for the latter object, IGMFs with B 10 −12 -10 −11 G and L B 10 kpc could be measured in the energy range between 1 and 10 TeV via halo searches. Such strong IGMFs could, in principle, be invoked to resolve the Hubble tension [61].

Outlook
Following on the footsteps of pioneer ground-based gamma-ray detectors, in particular IACTs like the Whipple Observatory and HEGRA, currently-operating facilities such as H.E.S.S., VERITAS, MAGIC have made outstanding progress in studying the VHE universe. Complemented by space-borne detectors like Fermi-LAT and AGILE at energies below ∼ 100 GeV, and by ground-based particle detectors such as HAWC, Tibet-ASγ, and ARGO-YBJ at higher energies (∼ 100 TeV), we have in recent years made significant progress towards understanding the Universe at high energies. At the dawn of the multimessenger era, the discovery potential of ground-based gamma-ray facilities can be maximised by working with other observatories across the whole electromagnetic spectrum, as well partners measuring cosmic rays, neutrinos, and gravitational waves. In this context, joint studies through multimessenger networks such as the Astrophysical Multimessenger Observatory Network (AMON) [410] are extremely useful to orchestrate campaigns of follow-up observations. It is through coordinated efforts of multiple of these facilities that we can pave new roads to fully exploit the potential of gamma rays as probes of cosmology and fundamental physics (e.g., Lorentz invariance violation, axion-like particles; see Sec. 3.5).
Within this landscape, we identify a unique opportunity for measuring IGMFs using gamma rays as well as other messengers.
Many challenges lie ahead in the coming decade. Firstly, it is possible that next-generation IACTs like CTA will still not be sensitive enough to enable measurements of magnetically-induced haloes. This limitation is certainly true for other ground-based instruments given the lower angular resolution of water-Cherenkov detectors compared to IACTs. Secondly, there are theoretical issues that need to be addressed, including the issue of plasma instabilities (see Sec. 3.4). Moreover, future studies should start relying on more detailed magnetic-field models, capturing also the magnetisation of the cosmic web wherein the gamma-ray sources used as "lighthouses" to probe the cosmos are embedded. We are entering an era of precision measurements and, therefore, also require more accurate tools to model the three-dimensional propagation of electromagnetic cascades if we wish to exploit the data as much as possible. Finally, there is room for novel methods to be devised to measure IGMFs, involving, among other messengers, gamma rays.
New insights into cosmic magnetism will be obtained with the Square Kilometre Array (SKA) [411,412]. Through measurements of Faraday of polarised extragalactic sources (e.g., FRBs, GRBs, quasars) SKA will deliver a tomographic map of extragalactic magnetic fields, disentangling part of the IGMF component, and offering clues on the structure and evolution of IGMFs [413]. Figures 1 and 9 neatly summarise the space of parameters for IGMFs allowed by measurements. However, the landscape of cosmic magnetism is more complex than this simple two-dimensional parameter space. Besides the magnetic-field strength (B) and the coherence length (L B ), IGMFs may be helical, such that a third dimension (H B ) should be added to this plot. Moreover, the magnetic power spectrum (α B ) can also play a role in the development of electromagnetic cascades, adding a fourth dimension. It is manifestly difficult to scan over all these parameters (B, L B , H B , and α B ) simultaneously. Still, these caveats should be borne in mind when constraining IGMFs, since there might be degeneracies. In this context, observation of gamma-ray sources can play an important role, given its ability to probe all these parameters. Nevertheless, besides more sensitive gamma-ray observatories, theoretical efforts in this direction are needed.
With the promising prospects for measuring IGMFs using next-generation gamma-ray observatories, we can invert the reasoning presented in Sec. 2.3: from the measurements, assuming IGMFs have a cosmological origin, we could constrain certain aspects of cosmology by inferring the specific parameters that characterise them. In fact, all the IGMF parameters mentioned in the previous paragraph may, in principle, be used for this purpose. With the measurement of B, one directly obtains the overall energy content of IGMFs. Like any other form of energy permeating the Universe, this would have an immediate impact in its global evolution, such that it could be necessary to consider this contribution as an addition to the standard ΛCDM model. Moreover, measurements of both the spectral index (α B ) and the coherence length (L B ) could be used to constrain the major processes from which they originate, like Inflation, QCDPT and EWPT. In the case of a phase transition, in particular, these measurements could allow us to infer its order. Finally, the measurement of a non-zero magnetic helicity (H B ) would strongly hint at a general CP violation in the Universe, with clear implications for various aspects of particle cosmology.
In summary, it is fair to say that gamma rays represent a unique observational window into cosmic magnetism. With the advances in gamma-ray astronomy, we could already capitalise on this window of opportunity to better understand IGMFs and to start constraining the B-L B parameter space. In the coming decades, the next generation of instruments might improve our understanding of cosmic magnetism more than ever, probing magnetism at cosmological scales and providing us a glimpse into the mechanisms whereby magnetic fields originated.