Gamma-ray Cosmology and Tests of Fundamental Physics

The propagation of gamma-rays over cosmological distances is the subject of extensive theoretical and observational research at GeV and TeV energies. The mean free path of gamma-rays in the cosmic web is limited above 100 GeV due to the production of electrons and positrons on the cosmic optical and infrared backgrounds. Electrons and positrons cool in the intergalactic medium while gyrating in its magnetic fields, which could cause either its global heating or the production of lower-energy secondary gamma-rays. The energy distribution of gamma-rays surviving the cosmological journey carries observed absorption features that gauge the emissivity of baryonic matter over cosmic time, constrain the distance scale of $\Lambda$CDM cosmology, and limit the alterations of the interaction cross section. Competitive constraints are in particular placed on the cosmic star-formation history as well as on phenomena expected from quantum gravity and string theory, such as the coupling to hypothetical axion-like particles or the violation of Lorentz invariance. Recent theoretical and observational advances offer a glimpse of the multi-wavelength and multi-messenger path that the new generation of gamma-ray observatories is about to open.


The extragalactic gamma-ray sky
The total electromagnetic emission from astrophysical objects in the Universe is dominated by thermal emission from stars, gas and dust. Their optical and infrared emission, originating from nuclear processes and thermal radiation, dominates the cosmic background radiation at wavelengths ranging from 0.1 µm to 100 µm [1]. These optical and infrared backgrounds are collectively referred to as the extragalactic background light (EBL). While millions of galaxies have been observed over the full-sky in the optical band, only a fraction of them also emits non-thermal radiation at radio and X-ray wavelengths [2]. Such non-thermal emission originates from charged particles, electrons and positrons or protons and nuclei, which are accelerated by violent processes in astrophysical environments. These charged particles radiate away their energy by gyrating in magnetic fields or by interacting with particles and lower-energy photons populating the surrounding medium [3]. The frontier of the electromagnetic spectrum up to which such astrophysical accelerators are observed is the gamma-ray domain, which is defined as the electromagnetic band covering energies larger than a few MeV.
Gamma-rays from astrophysical sources are observed with dedicated particle-physics instruments either in space or on the ground. The most successful technique at energies from 100 MeV to 1 TeV is based on detectors orbiting Earth, such as Fermi-LAT [4], that combine a tracker and a calorimeter in which gamma-rays generate showers of electrons and positrons. The hits and charge deposits from the secondary leptons are extracted to reconstruct the arrival direction and energy of the primary gamma-ray. The effective area of such space detectors, on the order of a square meter, is a limiting factor beyond few hundreds of GeV. The non-thermal emission from astrophysical sources follows to first order a power-law spectrum arXiv:2202.00523v1 [astro-ph.CO] 1 Feb 2022 of energy, with a differential number of photons per infinitesimal energy band that strongly decreases with energy, E, typically as E −2 . In the very-high-energy range, from 100 GeV to 100 TeV, the atmosphere of Earth itself is used as a calorimeter, so that ground-based telescopes separated by one hundred meters, such as H.E.S.S., MAGIC and VERITAS, build up an effective collection area ten thousand times larger than that of gamma-ray satellites [5].
The field of high-energy astrophysics, dedicated to the study of non-thermal emission, is relatively young compared to thermal astrophysics (see, e.g., Refs. [6,7] and references therein for a historical perspective). The development of the radar technique during World War II enabled in the 1950s the first detections of radio galaxies such as Centaurus A (at a luminosity distance d ∼ 4 Mpc), M 87 (d ∼ 16 Mpc) and Cygnus A (d ∼ 230 Mpc). Such radio observations, combined with optical spectroscopy, triggered the development of multiwavelength extragalactic astronomy in the early 1960s, with the first measurements of redshifts of quasars at cosmological distances, namely 3C 273 and 3C 48 at redshifts z = 0.16 and z = 0.37, i.e. d ∼ 770 Mpc and d ∼ 2000 Mpc, respectively [8]. The late 1960s and early 1970s revealed that some of these radio galaxies and quasars, such as M 87, Centaurus A and 3C 273, also emit X rays. The quasar 3C 273 was first detected in the gamma-ray band in 1978 by the CoS-B satellite [9] and remained the only known extragalactic gamma-ray source for over a decade. The non-thermal luminosity of such quasars, nowadays encompassed together with radio galaxies in the more general class of active galactic nuclei (AGNs [2]), goes up to ∼ 10 41 W which is equivalent to hundreds of trillions of solar luminosities, that is thousands of times larger than the stellar emission from the host galaxy. The tremendous luminosities of AGNs like 3C 373 are enabled by the relativistic beaming and energy shift of photon fields from jetted outflows. If the jet is closely aligned with the line of sight, the AGN is called a blazar. The observed photon energy is then enhanced with respect to the emitted photon energy in the jet frame, by a Doppler factor δ ∼ 10. The observed bolometric luminosity is further enhanced by a factor δ 4 [3]. The gamma-ray band is of particular interest in the study of jetted AGN for two reasons: not only does it probe the highest energies per photon but it also often encompasses more than half of the electromagnetic luminosity of AGNs, thus providing an intense and energetic beam to study cosmological processes and fundamental physics.
Extragalactic gamma-ray astronomy was established as a disciplinary field in its own right during the 1990s, with the successful operation of the Whipple telescope on Mount Hopkins (Arizona, USA) and the launch of EGRET on board the CGRO satellite. Figure 1 illustrates the growth of the number of gamma-ray sources identified outside of the Milky Way since that epoch. At the turn of the millennium, space-borne EGRET observations had enabled the firm identification over the entire sky of nearly a hundred extragalactic sources above 100 MeV (see 3EG Catalog [10]). Cross-matches with lower-energy multiwavelength observations revealed that most of these accelerators are AGN, although without specific sub-classification (see AGN, AGU, BCU in 3EG pie chart of Figure 1). The only non-AGN extragalactic source was the Large Magellanic Cloud, which is the Milky Way's largest satellite galaxy at d ∼ 0.05 Mpc. Meanwhile, the first ground-based gamma-ray telescopes, with their limited field of view of a few degrees, were pointed to the most promising targets selected among radio, optical and X-ray extragalactic sources. The first extragalactic gamma-ray source observed above few hundreds of GeV is the blazar Mrk 421, at z = 0.031 or d = 140 Mpc, which was detected in 1992 [11]. A single other gamma-ray source had been discovered by the Whipple telescope three years earlier at such energies, the Crab Nebula located within the Milky Way at a distance of about 2 kpc from Earth [12]. With a flux in March-June 1992 estimated to 30% that of the Crab, even a relatively low emission state 1 of Mrk 421 proved blazars to be extremely bright gamma-ray beacons distributed on extragalactic scales. Seven other blazars and a single radio galaxy, M 87, had been detected by 2005 with the Whipple telescope and its successors such as HEGRA (La Palma, Canary Islands) and CAT (Pyrenees, France) [14]. 2 The late 2000s saw the emergence of a new generation of gamma-ray telescopes, H.E.S.S. (Khomas Highlands, Namibia), MAGIC (La Palma, Canary Islands) and VERITAS (Arizona, USA), with lower energy thresholds down to 100 GeV, as well as improved background rejection and imaging capabilities (68% containment angle better than 5 arcmin above 1 TeV). They revealed a population of TeV extragalactic sources dominated by blazars of the BL Lac type (BLL). Unbeamed and steady TeV emission was also discovered from a few galaxies: the radio galaxy Centaurus A and the starburst galaxies NGC 253 and M 82, all of which are located in the Council of Giants that surrounds the Milky Way at d = 3 − 6 Mpc [15]. The launch of the Fermi-LAT satellite in 2008, with a gain of about a factor of five in effective area and in field of view with respect to EGRET, triggered a tremendous growth in the number of extragalactic gamma-ray sources detected above 100 MeV. With an angular resolution (68% containment angle) better than a degree above 1 GeV [16], cross-identification with radio, optical and X-ray catalogs established the prominence of blazars in the gamma-ray sky, although with marked differences with the emerging population uncovered by ground-based telescopes at higher energies. About half of the blazars categorized in the first Fermi-LAT catalog ( [17], 1FGL in Figure 1) were found to be flat-spectrum radio quasars (FSRQs), at variance with the prevalence of BLLs at higher energies. The relatively high accretion rate and bright photon fields in the environment of FSRQs with respect to BLLs is nowadays understood as the probable cause for the fainter but higher-energy peak luminosity of BLLs [18].
As shown in Figure 1, the number of extragalactic sources identified above 100 MeV and above 100 GeV increased by nearly order-of-magnitude in the 2000s. After an initial phase with a high discovery rate, a saturation in number of discovered sources was prevented in the 2010s thanks to improved observing strategies and analysis techniques, continued operation of Fermi-LAT and upgrades of H.E.S.S., MAGIC and VERITAS to lower the energy thresholds of the observatories below 100 GeV. Among the many discoveries enabled by relentless efforts from the gamma-ray community, one could highlight the discovery of the gravitationally-lensed FSRQ B2 0218+35 located at z = 0.954 (d ∼ 6000 Mpc), a distance record for AGNs observed above 100 GeV, only recently superseded by the detection of an FSRQ announced at z = 0.991 [25]. For comparison, blazars observed by Fermi-LAT have been firmly detected above 100 MeV from beyond z = 3 and tentative detections have even been reported out to z > 4 [26]. Other interesting discoveries lie in the few dozen objects above 100 MeV that do not belong to the classical types of extragalactic gamma-ray sources, namely blazars and radio galaxies. For example, narrow-line Seyfert 1 (NLSY1) galaxies can display jets with an inferred power similar to that of BLLs but with presumably brighter photon fields, as expected in FSRQs [27,28]. The inferred black-hole masses within NLSY1 are estimated around ten million solar masses, rather than billions for numerous blazars and radio galaxies, which would suggest a high accretion rate for NLSY1. Together with compact steep spectrum sources and steep-spectrum radio quasars, NLSY1 could provide one of the missing links in the understanding of jet formation around super massive black holes.
Last but not least, the 2010s and 2020s have seen the emergence of gamma-ray bursts (GRBs) as an entirely new extragalactic population at very-high energies. While the emission of other extragalactic sources is often thought to stem from galactic-scale outflows, be they jets for active galaxies or collective winds for starforming ones, the power engine of GRBs is of Satellites at E > 100 MeV Telescopes at E > 100 GeV Figure 1. Evolution of the number of extragalactic gamma-ray sources at high energies as a function of time. The extragalactic sources are usually identified as such by matching them with multiwavelength counterparts. The dashed gray line shows significant detections by the EGRET and Fermi-LAT gammaray satellites at energies larger than 100 MeV. The pie charts display the distribution of sub-classes of extragalactic sources, as extracted from the 1EG [19], 3EG [10], 1FGL [17], 3FGL [20], 1FLGC [21], 4FGL-DR2 [22], 1FLT [23] and 2FLGC [24] catalogs. The solid black line shows significant detections by gamma-ray telescopes above 100 GeV as extracted from TeVCat [14] after the 2021 international cosmic-ray conference (ICRC 2021 stellar size and their outbursts are understood as being induced by the collapse of short-lived massive stars or by a binary merger. Short GRBs, with a duration below a couple of seconds, are subdominant in the gamma-ray band (less than 10% above 100 MeV, none above 100 GeV) with respect to long GRBs [24]. The latter are now detected up to z = 4.35 (GRB 080916C [29]) by gamma-ray satellites above 100 MeV and up to z = 1.1 (tentative redshift of GRB 201216C [30]) by gamma-ray telescopes above 100 GeV.
Multi-wavelength observations of extragalactic gamma-ray sources provide remarkable insights into the non-thermal processes at play in astrophysical environments, as discussed in other reviews of this Special Issue on Extragalactic TeV Astronomy. Of particular interest to the present discussion are source redshifts, whose spectroscopic measurements can be challenging when the emission lines from the host galaxy are overwhelmed by non-thermal emission from the outflows. Despite these difficulties, dedicated spectroscopic campaigns, preceding, coordinated with, or following up the gamma-ray discoveries (see Ref. [31] and references therein), have enabled to constrain the distances of over 80% of the extragalactic gamma-ray sources detected above 100 GeV. These redshift measurements and gamma-ray detections act as the starting blocks of the journey of gamma-rays over cosmological distances, whose multiple paths are discussed in the following sections.

A gamma-ray journey through cosmic ages
As described in the introduction, gamma-rays with energies of at least 100 GeV have now been observed from sources (mostly blazars) over a large range of distances, out to z ≈ 1 from the ground and z ≈ 4 from space-borne instruments. The observation of gamma-rays from such distant emitters provides us with a unique opportunity to study the propagation of the highest-energy electromagnetic radiation over cosmological distances [32]. In this section, we review the interactions that gamma-rays can undergo along their journey to Earth. Just as the absorption of optical emission can be used to infer properties of, e.g., stellar or planetary atmospheres [33], the absorption of gamma-rays in the intergalactic medium is a complementary tool to study the cumulative energy release by stars and dust grains over the history of the Universe. Furthermore, the high energies and long distance scales involved might enable us to observe effects of physics beyond the standard model, such as oscillations into hypothetical scalar particles or the departure from Lorentz invariance. The processes discussed in this section are shown in Figure 2, which illustrates the path of a gamma-ray from its creation in the jet of a blazar (we take here the example of an FSRQ) until it enters the Milky Way. For illustration, we follow the path of three gamma-rays produced in the blazar jet (left-most panel). The top gamma-ray travels unperturbed until it is detected on Earth. The second gamma-ray in the center is absorbed in intergalactic space, whereas the bottom gamma-ray gets transformed into a hypothetical axion-like particle (ALP) and converts back into a gamma-ray. We start our discussion with the absorbed gamma-ray.
It was realized in the 1960s that gamma-rays produced in distant sources should undergo absorption in collisions with photons from background radiation fields [39] such as intergalactic starlight [40], radio-frequency fields [41], the cosmic microwave background (CMB) [42,43], or radiation fields within the source [44]. The interaction of a high-energy gamma-ray, γ, with a photon from a background radiation field, γ bkg , produces an electron-positron pair, γ + γ bkg → e + + e − if the following threshold condition is met [45] E ≥ 2 m e c 2 2 1 − cos θ .
In Eq. (1), m e is the mass of the electron (and positron), E and are the energies of the gamma-ray and background photon, respectively, which are denoted with a prime in the comoving cosmological frame, and θ is the angle between the momenta of the two photons. and a radius equal to the Schwarzschild radius 2 r g with r g = GM • /c 2 ≈ 1.5 × 10 14 cm ≈ 4.8 × 10 −5 pc. The inner and outer radii of the accretion disk, broad line region (BLR), and dusty torus are based on typical values from the literature with R disk, in = 6 r g , R disk, out = 200 r g , R BLR, in = 0.01 pc, R BLR, out = 0.1 pc, R torus, in = 5 pc, R torus, out = 23 pc (e.g., Ref. [34]). The jet is assumed to have a parabolic base with a transition into a conical shape at 2 × 10 5 r g (see Eq. (13) and Figure 1 in Ref. [35]). Not shown here is the possibility of internal absorption of gamma-rays on the photon fields within the source (e.g., on photons originating from the BLR or the dusty torus). Several observations of FSRQs indicate that gamma-rays are produced outside the BLR (see, e.g., Refs. [36,37] and references therein). Center: Once produced, the gamma-rays propagate through the host galaxy (r host = 20 kpc [38]) and potentially through a Galaxy cluster (with r cluster = 1 Mpc here) before entering cosmic voids. The upper gamma-ray propagates all the way towards the observer whereas the central gamma-ray interacts with an EBL photon to form an electron-positron pair at a distance equal to the mean free path, Γ γγ ≈ 250 Mpc, of a 1 TeV gamma-ray. The electron and positron gyrate in the IGMF with r gyro ≈ 540 kpc (γ/10 6 )(B/10 −15 G) −1 and inverse-Compton scatter off CMB photons after a mean free path of Γ eγ ≈ 730 kpc (γ/10 6 ) −1 . The lower gamma-ray converts into an ALP in the magnetic field of the galaxy cluster. Right: The surviving gamma-ray enters the Milky Way and the ALP converts back to a gamma-ray in the magnetic field of our Galaxy.
The Breit-Wheeler cross section for pair production, σ γγ , is a function of the velocity of the electron (positron) in the center of mass frame with β 2 = 1 − 2(m e c 2 ) 2 /[E (1 − cos θ )] and peaks for β ≈ 0.7. Averaging over the (1 − cos θ ) term, one finds that the peak of the cross section is reached when ≈ 1.0 eV(E/TeV) −1 or for a wavelength λ ≈ 1.2 µm(E/TeV). As a consequence, gamma-rays at energies E > 100 GeV most likely interact with photons at optical to infrared wavelengths. The absorption process is shown in the central panel of Figure 2 for the central gamma-ray line. For a specific photon density dn/d (in units of eV −1 cm −3 ), the mean-free path, Γ γγ , for pair production occurring at a redshift z is given by where Θ is the Heavyside step function which ensures that the threshold condition in Eq. (1) is met. In Section 3, we discuss the EBL, which is the radiation field dominating the absorption of gamma-rays. For reference, the mean free path of a 1 TeV gamma-ray interacting with EBL photons is Γ γγ ≈ 250 Mpc for a source at d = 1 Gpc (z ≈ 0.2) assuming the EBL model in Ref [46]. This is the distance at which the pair production occurs in the central panel in Figure 2. The gamma-ray absorption, quantified by the optical depth τ γγ , is then found as a line-of-sight integral up to the redshift z 0 of the source over the mean-free path, where d /dz is the Jacobian for the transformation of the line-of-sight integral over distance to redshift, which depends on the adopted cosmology. Throughout this review, we assume a flat ΛCDM cosmology with reference values for the matter density Ω M = 0.3, dark-energy density Ω Λ = 0.7 and Hubble constant at current epoch H 0 = 70 km s −1 Mpc −1 , so that The initial flux emitted by the blazar, F emitted , is attenuated and we observe the flux F obs = F emitted exp(−τ γγ ). As discussed in Section 3, we can use the observation of sources at different energies and redshifts to reconstruct τ γγ and thus reconstruct the EBL photon density, dn/d .
In the pair-creation process, the electrons and positrons acquire about half of the initial gamma-ray energy, so that their Lorentz factors reach γ = E/(2m e c 2 ) ≈ 10 6 (E/TeV). If we follow the path of the electrons and positrons in Figure 2, we see that they can undergo an inverse-Compton scattering process with background radiation fields. The most abundant photon field is the CMB with an integrated energy density of U CMB ≈ 0.26 eV cm −3 and an average photon energy of CMB ≈ 630 µeV. For such CMB photons and for electrons at γ ≈ 10 6 , the inverse-Compton scattering occurs in the Thomson regime, i.e. the electrons only loose a small fraction of their energy in each scattering event (see, e.g., [47]). In the Thomson regime, the mean free path Γ eγ and the average photon energy after scattering, ˜ , are where σ T ≈ 6.65 × 10 −25 cm 2 is the Thomson cross section and t cool = γ/|γ| is the cooling time. Two observations can be made from these equations. On the one hand, for a sufficiently large primary gamma-ray energy, E, the upscattered CMB photon can reach an energy ˜ that is in the gamma-ray regime, so that the secondary gamma-ray could again pair produce on a background photon and initiate a full electromagnetic cascade. On the other hand, Eq. (4) indicates that the electron-positron pairs can travel over hundreds of kpc before they scatter. As they are charged, they are deflected in magnetic fields along the line of sight. A sufficiently large magnetic field would then broaden the pair beam. The resulting gamma-ray cascade would arrive with a time delay with respect to the primary gamma-ray and could even form an extended halo around the otherwise point-like blazars [48,49]. These effects have been suggested as potential probes of the strength of the intergalactic magnetic field (IGMF), which is extremely difficult to observe directly (see, e.g., Ref. [50]). The current constraints are discussed in Section 3.
We now turn to processes beyond the standard model that can affect gamma-ray propagation. One possibility is that gamma-rays oscillate into axions in the presence of external electromagnetic fields. The axion was originally proposed to explain the apparent absence of the violation of parity symmetry (P) and time reversal (T) in QCD [51][52][53]. It was soon realized that this axion is also a dark-matter candidate [54][55][56][57] and that particles with similar properties, so-called axion-like particles (ALPs), arise in extensions of the standard model, most notably string theory (see, e.g., Ref. [58] for a review). The interaction between axions (or ALPs) and photons arises through the Lagrangian density L aγ = g aγ E · Ba, 3 where g aγ is the coupling between axions and photons (in units of inverse energy), E is the electric field associated to the photon, B is an external magnetic field, and a is the axion field strength [55,60]. This Lagrangian term describes the photon-axion Primakoff effect which, similarly to neutrino flavor mixing, results in an oscillation between photons and axions in the presence of external fields.
If produced through the so-called misalignment mechanism in order to make up a fraction Ω a, DM of the total dark matter density (as determined, e.g., with the Planck satellite to be Ω Planck DM = 0.26 [61]), the photon-ALP coupling must fulfill the relation (adopted from Ref. [57]), where θ 1 is the initial misalignment angle of the ALP field in the early Universe (which takes values between −π and π) and N is a model dependent parameter; in the simplest case, we take θ 1 N ∼ 1. The Feynman diagram for the photon-axion interaction is shown in the bottom central panel of Figure 2. We used the example of a gamma-ray converting into an ALP in the magnetic field of a galaxy cluster that might harbor the blazar. The ALP can then convert back into a gamma-ray in, e.g., the magnetic field of the Milky Way, as discussed in Section 4. Photon-axion oscillations could be detected in gamma-ray observations of blazars mainly in two ways. On the one hand, energy-dependent oscillations could be visible around two specific energies E low crit and E high crit (see, e.g., Ref. [62,63] with the addition from Ref. [64]), where m a is the mass of the axion; ω plasma is the plasma frequency of the traversed medium with electron density n el ; 13 G is the critical magnetic field and α is the fine-structure constant; χ photon describes the strength of the dispersion of the photon from background radiation fields. If the only relevant field is the CMB with [64]. Below E low crit , the oscillations are suppressed due to a momentum mismatch between the (massive) axion and the photon with an effective mass equal to the plasma frequency. Above E high crit , the mixing is suppressed as well: photon-photon dispersion (i.e. higher order corrections to the photon propagator) with the magnetic field and / or background radiation fields becomes more likely than photon-ALP oscillations. On the other hand, for gamma-ray energies E low crit < E < E high crit the mixing becomes maximal and independent of energy. For a path length L in a homogeneous magnetic field, the oscillation probability is simply given by P aγ ≈ (2πL/λ osc ) 2 with the oscillation length λ osc ∼ 4π/(g aγ B) ∼ 0.4 Mpc (B/µG) −1 (g aγ /10 −11 GeV −1 ) −1 . Interestingly, in this strong mixing regime, a fraction of the gamma-rays could avoid pair production by converting into axions if λ osc Γ γγ . If these axions re-convert into photons close to Earth, the observed photon flux would be enhanced in comparison to the no-axion case and the Universe would appear more transparent to high-energy gamma-rays than expected from conventional physics only. The status for axion searches using both effects is reviewed in Section 4.
Lastly, gamma-ray propagation could also be altered if Lorentz invariance is violated. Such a Lorentz invariance violation (LIV) is expected from theories that try to unify quantum field theory with gravity (see Ref. [65] and references therein). A recent review of the astrophysical effects of LIV is given in Refs. [66,67], which we loosely follow hereafter. One of the consequences of LIV would be that the photon velocity is not always the speed of light, c, but depends on the energy of the photon. Some effective-field theories predict that LIV can lead to a modified dispersion relation for particles at leading order n [68]. For photons with energy E and momentum p, this can be written as where ξ n is a coefficient of the underlying theory and M is the energy scale at which new physics is expected, e.g., the Planck energy scale E Pl = √h c 5 /G ≈ 1.22 × 10 19 GeV. Tight bounds have been placed on the superluminal case (plus sign in Eq. (9)), in particular by recent observations by HAWC and LHAASO of gamma-rays beyond 100 TeV from sources in the Milky Way [69,70]. We focus the discussion here on subluminal searches, i.e. a negative sign in Eq. (9), which correspond to velocities below the speed of light for photons. Within specific models predicting LIV, such as the Standard-Model Extension [71], corrective terms with an even order n emerge from operators that are invariant under charge/parity/time-reversal (CPT) symmetries, contrarily to terms with an odd order n [72,73]. Quadratic modifications of the dispersion relation may thus be favored over linear modifications within CPT-invariant theories.
In terms of gamma-ray propagation, the modified dispersion relation has two important consequences. On the one hand, the threshold for pair production is modified by adding the additional term |ξ n |E n+2 /(4M n ) to the right hand side of Eq. (1), where ξ n = 1 if LIV only affects photons and ξ n = (1 − 2 n ) −1 if it affects both photons and electrons [74,75]. As a result, pair production will cease once the ratio E n+2 /M n becomes large compared to (2m e c 2 ) 2 . The latter two terms are comparable at E ∼ 20 TeV for M = E Pl and n = 1, so that the optical depth τ γγ would decrease towards zero at gamma-ray energies larger than E . Thus, if one can measure blazars beyond these energies, it should be possible to probe LIV.
On the other hand, an energy-dependent speed of light would lead to a time delay ∆t between two photons emitted simultaneously with energy difference ∆E > 0 propagating from a source at redshift z (e.g., Ref. [76] but see also Ref. [77] for an alternative formulation), (10) We have approximated the integral with polynomials F n (z) = ∑ n+1 k=1 p k z k of the (n + 1)-th order for 0.01 z 4 which show an absolute error of less then 4 % compared to the numerical integration. In the case n = 1, F 1 (z) is almost linear with F 1 (z) ≈ 1.2 z whereas, for n = 2, F 2 (z) increases from ≈ 10 −4 for z = 0.01 to ≈ 9 for z = 4. 4 Sources as distant as possible and electromagnetic radiation over a wide (linear) energy range are preferred to maximize the potentially observable effect. As discussed in Section 4, gamma-ray bursts and blazar flares are best suited for these kinds of searches.

Probing the content of the intergalactic medium
The baryonic and electromagnetic contents of the Universe can be probed through observations of extragalactic gamma-ray sources at various distances, as illustrated in Figure  3. The first entity of interest consists of all ultra-violet to infrared photons emitted since reionization, from a time when the Universe was about 1 Gyr old to the current epoch, about 13 Gyr later. Such photons accumulate in the intergalactic medium to build up the main constituent of astrophysical origin in the "spectrum" of the Universe [1]: the EBL. The EBL energy density, dn/d , is composed of two humps each spread over nearly two decades in energy: the cosmic infrared background (CIB) at wavelengths 8 µm < λ < 1000 mm or energies 0.16 eV > > 1.2 meV, which mainly arises from emission by dust grains, and the cosmic optical background (COB) at wavelengths 0.09 µm < λ < 8 µm or energies 13.6 eV > > 0.16 eV, which mainly stems from stellar radiation that escaped its environment without absorption on dust grains. As such, the evolution of the EBL with redshift provides an integrated measure of the evolving star-formation rate density, also called the cosmic star-formation history (CSFH [90]). The second entity of interest is the magnetic field which pervades the space between filaments and clusters of galaxies: the IGMF [50]. The IGMF is a crucial ingredient to understand cosmic magnetism on both large and small spatial scales, particularly as it could have provided the seed later amplified through compression and 10 6 10 0  Figure 3. Evolution within a Hubble radius of the two classical entities relevant to gamma-ray propagation, the IGMF and the EBL, as a function of cosmic age. The redshift range over which constraints on gamma-ray propagation have been established by ground-based telescopes (E > 100 GeV) is illustrated by the upper white box in the right-hand-side panel. The latter observations are exemplified by the extreme blazar of BLL type 1ES 0229+200 (z = 0.14), which has been seminal for IGMF studies [78,79], and by the most distant TeV blazar used to constrain the EBL, the FSRQ PKS 1441+25 [80,81]. The redshift range over which constraints have been established by satellite-based observations (E > 100 MeV) is illustrated by the lower white box, exemplified by the FSRQ-type PKS 1502+106 (z = 1.84) which is the most-distant blazar detected by Fermi-LAT above 30 GeV [82] with a firm spectroscopic redshift [83]. The emission of the EBL from a cosmic age of about 1 Gyr is illustrated by the colored area, which figures the star-formation rate density inferred from gamma-ray observations [84], as parameterized by Ref. [85]. Earlier phases of the expansion of the Universe are illustrated in the middle panel, which features the decoupling that resulted in the CMB emission, and in the left-hand-side panel, which includes the quantum chromodynamics (QCD) phase transition (or cross-over) [86] while the electroweak (EW) phase transition (or cross-over) [87] occurred at earlier times. As illustrated by the magnetic cells shown in the first two panels, both the EW and QCD eras could have seen the IGMF emergence, but it should be noted that the IGMF could also have been seeded at later epochs during the formation of large-scale structures. The reader is referred to Figure 3 in Ref. [88] and Figure 9 in Ref. [89] for more specific constraints on the EBL spectrum and IGMF parameter space, respectively. dynamo effects into the µG fields observed in galaxies and clusters [91]. The IGMF could originate both from astrophysical outflows (AGN jets, starburst winds) or primordial phase transitions, e.g., at the break of the electroweak symmetry [87]. The reader is referred to Ref. [88] for a didactic review on these matters. As discussed in Section 2, gamma-rays can interact with photons from the EBL en route to the observer, which results in an absorption of the emitted spectrum. The magnitude of the absorption depends both on the cross section of interaction and on the density of target EBL photons in the wavelength range of interest. For example, spectral observations of nearby blazars such as Mrk 421 and Mrk 501 (z ∼ 0.03) at 20 TeV probe the EBL density at typical wavelengths ranging in 10-100 µm, that is a fraction of the CIB, while observations at 100 GeV of sources at z = 1 typically probe the COB density at wavelengths ranging in 0.1-1 µm. Gamma-ray observations of sufficiently bright sources at different redshifts thus enable the study of the EBL intensity over nearly three decades in energy (0.1-100 µm). The main limitation of the indirect gamma-ray measurement technique lies in the unknowns on the emitted spectrum, which were historically solved as follows.
The first studies of the EBL based on observed gamma-ray spectra were developed after the discovery of Mrk 421 by the Whipple observatory. The TeV spectrum was jointly modeled with the unabsorbed spectral measurement from EGRET at MeV-GeV energies, assuming that the emitted spectrum follows a straight power-law dependence with energy from MeV to TeV energies [92]. As soon noted by Ref. [93], such an approach neglects any possible intrinsic curvature in the emitted spectrum, which could be attributed either to absorption inside the source itself rather than on the line of sight or to the interplay between acceleration and radiative processes that would reduce the flux at very-high energies. Although only weak upper limits on the EBL density could be placed with the first observations, constraining limits were derived a decade later with the observations of the extreme blazars H 2356+304 (z = 0.165) and 1ES 1101-232 (z = 0.186) by H.E.S.S. [94]. At the lack of MeV-GeV observations to constrain the maximum hardness of the emitted spectrum, the authors of Ref. [94] used a theoretical bound on the index of the power-law spectrum and managed to solve a longstanding discrepancy between direct observations of dark patches of the night sky and the integrated galaxy light (IGL), which is determined from counting galaxies in deep-field surveys. Direct observations encompass not only the EBL but also one-to-two order-ofmagnitude brighter foregrounds (zodiacal light from dust grains in the Solar system, stellar light from our Galaxy), so that a percent bias in foreground subtraction could result in a 100% overestimation of the EBL [95]. Galaxy counts, on the other hand, rely on well-calibrated but limited samples of galaxies. These surveys could fail to account for low-surface-brightness galaxies below the sensitivity threshold of the instrument or miss any diffuse component of the EBL [96]. Up to the late 2000s, the most optimistic direct measurements exceeded the IGL by over an order of magnitude. This "optical controversy" [96] appeared to be solved by indirect constraints from H.E.S.S., with gamma-ray upper limits suggesting that most of the COB was resolved in known galaxies.
The "optical controversy", or more specifically the question of the actual intensity of the COB between 0.5 µm and 1 µm, is not simply a matter of discrepant measurement techniques (see Refs. [96,97] for a recent status of the scientific debate). This region of the EBL spectrum could contain contributions from the sources of reionization (UV photons at ∼ 0.1 eV emitted at z ∼ 10), from gas or stars stripped from their host galaxies, e.g., during mergers, or radiative signatures from exotic particles such as decaying or interacting dark-matter candidates. Moreover, the limitations of the gamma-ray indirect technique, namely the unknowns in the emitted spectrum, left trying questions open. Is the theoretical bound on intrinsic hardness sufficiently motivated? With only gamma-ray upper limits at hand, would an inferred EBL level below the IGL suggest a too large gamma-ray transparency, as expected from exotic processes?
The conceptual ditch was crossed after the launch of Fermi-LAT in 2008. Broad-band measurements above 100 MeV revealed that the unabsorbed gamma-ray spectra of blazars do show intrinsic curvature and that it can be disentangled from line-of-sight processes [98]. The techniques employed to account for intrinsic curvature in the first gamma-ray measurements of the EBL [98,99] slightly differed, as discussed in Ref. [100]. In Ref. [98], the Fermi-LAT Collaboration measured intrinsic curvature through a successful modeling of the unabsorbed part of each spectrum in their sample with a log-parabola function (parabola in log-log space, when a power-law model is a linear function in log-log space). The authors then modeled the absorbed part as F obs = F emitted exp[−aτ ref (E, z 0 )], by fixing the parameters of the emitted spectrum based on the lower-energy observations and letting free the normalization, a, of reference EBL models. In Ref. [99], the H.E.S.S. Collaboration, having a more limited access to the emitted spectrum, left both the EBL normalization and intrinsic parameters free to vary. The optimal level of curvature was determined by testing several curved variations of powerlaw models and by adopting the preferred one in a frequentist approach (maximization of a goodness-of-fit estimator). The absorption feature, which shows a characteristic dependence on energy and redshift, could thus be reconstructed by combining observations from multiple sources at different redshifts, with a 6σ signal driven by blazars at z ≈ 0.5 − 1.6 observed by Fermi-LAT up to ∼ 500 GeV and a 9σ signal driven by blazars at z ≈ 0.1 − 0.2 observed by H.E.S.S. up to ∼ 5 TeV. The specific intensity of the COB, νI ν = c 4π 2 dn/d (units of nW m −2 sr −1 ), was determined to be at a level compatible with the IGL within statistical and systematic uncertainties, on the order of 20 − 30% each. H.E.S.S. and Fermi-LAT confirmed, for the first time through a measurement, previous exclusions of the most optimistic direct measurements [94,101]. The analysis techniques have been refined over the past decade and larger datasets have been scrutinized. Ref. [102] compiled a database of 86 archival spectra from 30 sources observed by all ground-based instruments, together with satellitebased constraints on the emitted spectrum, and detected the absorption signature at the 11σ level based on an EBL scaling approach similar to those employed by H.E.S.S. and Fermi-LAT. A model-independent approach was also developed in Ref. [102] to measure the EBL intensity in four distinct wavelength bands spanning 0.3 − 100 µm. Imposing that the inferred EBL spectrum does not drop below the IGL level provided a refined measurements with 8 bins over the same wavelength range. Both EBL-model-dependent and -independent techniques were subsequently employed by the H.E.S.S. [103], VERITAS [104] and MAGIC [105] Collaborations, the latter also including Fermi-LAT constraints on the intrinsic spectrum. Although the overall EBL intensity can now be estimated from ground-based gamma-ray observations with a statistical accuracy of 10-20% and a systematic one of 20-30%, wavelengthresolved measurements come with a poorer resolution, particularly in the region of the "optical controversy". The status of the controversy at the time of writing is the following. The New Horizons mission recently performed a direct measurement of the COB beyond Pluto's orbit [97], estimating the EBL intensity at 0.6 µm to 15.9 ± 4.2 nW m −2 sr −1 (stat. + sys.). This value can be compared to the most accurate IGL estimate to date at a comparable wavelength (r-band), 8.11 ± 0.33 nW m −2 sr −1 (stat. + sys.), with an impressive accuracy on the order of 5% obtained through the combination of the GAMA and DEVILS deep-field surveys with observations from the Hubble space telescope [85]. The tension between galaxy counts and direct measurements has thus been reduced from a factor of ten in the early 2010s down to a factor of two, with a low-level significance if we loosely treat the uncertainties as random ones. For comparison, the latest wavelength-resolved measurements from H.E.S.S. and VERITAS limit the maximum EBL intensity at comparable wavelengths to 15 − 25 nW m −2 sr −1 (sys. limited), with estimates that, although favoring galaxy counts, do not settle the debate and leave ample room for unresolved contributions.
An alternative approach to the forward-folding gamma-ray spectral analysis described above has been developed over the past years, particularly to explore the evolution of the EBL with redshift and to constrain cosmological parameters. The EBL density indeed results from the accumulated luminosity density or specific emissivity, j (units of cm −3 s −1 ), of each redshift layer since the ignition of the first stars, so that dn/d = (1 + z) 3 /c ∞ z dz d /dz j( , z )/ , where = (1 + z ) [84]. For a given luminosity density, the mean dust extinction, A , and the amount of light emitted per star-formation-rate unit, K (units of M ), provides an estimate of the CSFH, ρ(z) = K 10 0.4A j( , z) (units of M yr −1 Mpc −3 ). With an estimate of the optical depth τ γγ (see Eq. (3)), gamma-ray spectral observations could constrain both the evolution of light emission through cosmic ages, ρ(z), and parameters of the ΛCDM cosmological model, in particular the Hubble constant via the distance element d /dz (see Ref. [88] for a review). Given the complexity of the emissivity-based calculation of the optical depth, which involves a double integral over redshift and an integral over EBL photon energies, the authors of Refs. [84,106] summarized ensembles of gamma-ray spectra (either from satellite-and/or ground-based observations) into estimates of the optical depth in successive energy bins, by fitting a scaling factor with respect to reference EBL models. The Fermi-LAT Collaboration in particular obtained the first gamma-ray indirect measurement of the CSFH through this technique [84]. Their CSFH and EBL estimates show good agreement with galaxy surveys (see also Ref. [85] for an independent discussion). While the gamma-ray constraints mostly come from blazars at z < 3, the integral nature of the gamma-ray measurements provides sensitivity to EBL emission at higher redshifts. The Fermi-LAT Collaboration estimated a UV luminosity density above z = 4 in line with the lowest values from Lyman-break galaxy surveys and some room remains at z ∼ 6 for a sufficiently bright UV field to drive reionization. In Ref. [107], the authors also employed estimates of the optical depth to constrain the Hubble constant and reported an accuracy of 3-6 km s −1 Mpc −1 on H 0 , where the bounds correspond either to a fixed or free matter density, Ω M . Such accuracy is relevant to test current tensions between late and early-type measurements of H 0 , typically distant by 5 km s −1 Mpc −1 with a tension at a significance level beyond 4σ [108]. The latest gamma-ray measurements are enticing as they open new scientific avenues to gamma-ray cosmology. Nonetheless, it can be noted that direct estimates of the optical depth thus far rely on scaled versions of EBL models that in turn depend on the CSFH and cosmological parameters, so that estimates of the latter are in principle not free from potential biases in the EBL models themselves (see, e.g., Refs. [109][110][111] for recent models). One of the main question in this line of research during the upcoming decade, independently from the adopted method, lies in the control on systematic uncertainties, which will require improvement upon the typical 20% resolution on the overall EBL intensity achieved as of today.
In a nutshell, the study of the EBL absorption signature in extragalactic gamma-ray spectra has shifted, over the past decade, from a detection quest to the pursuit of precision measurements. The observational proof of gamma-ray interactions with the EBL naturally leads to the question of how the produced electron-positron pairs dissipate their energy and of the potential observational signatures of such an energy release. As discussed in Section 2, a prime candidate process for dissipation is inverse-Compton scattering with background photon fields. With a bolometric intensity about twenty times larger than that of the EBL, the CMB provides copious amounts of target photons that can be upscattered to GeV energies by the secondary charged particles (see Eq. (5)). As shown in Eq. (4), such charged secondaries are produced at distances tens to hundreds of times larger than the typical Mpc scale that corresponds to the radius of the cluster or thickness of the filament hosting the gamma-ray source, so that it was early noted [49] that the detection of a GeV secondary signal could be used as a probe of the IGMF. The latter denotes the largely unknown magnetic fields in cosmic voids, whose volume filling factor is dominant in the cosmic web (see Ref. [112] for a didactic review on voids). An alternative candidate process to inverse-Compton scattering lies in the development of plasma instabilities in the beam of electrons and positrons (see dedicated overview in Ref. [89] and references therein). The nature and growth of these instabilities depend on the properties of the intergalactic medium (temperature, density) and on the intensity, energy distribution and angular distribution of the pairs in the beam (see Ref. [113] for state-of-the-art simulations and analytical treatment). The present understanding of the competition between the two cooling mechanisms remains limited. If plasma instabilities were to develop faster than inverse-Compton cooling and over a sufficiently long time (hundreds of years based on Eq. (52) in Ref. [113]), the relaxation of the beam could contribute to the heating of the intergalactic medium [114][115][116]. Observational constraints on the temperature history of the intergalactic medium nonetheless tightly limit the contribution of the secondary pairs, e.g., to the reionization of helium at z > 2 − 3 [117]. Clear-cut observational signatures of energy losses through plasma instabilities thus currently remain out of reach.
On the other hand, the dissipation of the energy of the pairs through inverse-Compton cooling would suggest the presence of a bump at GeV energies in the spectrum of sufficiently hard gamma-ray emitters. The best candidates for searches of GeV bumps are extreme blazars as their emitted flux is expected to show the largest TeV-to-GeV ratio among known extragalactic sources [118]. The detection of the extreme blazar 1ES 0229+200 (z = 0.14 or d ≈ 650 Mpc) up to 10 TeV [119] triggered a series of joint GeV-TeV studies, with up-to-date reviews presented in Refs. [88,89,120]. After a decade of searches, no secondary GeV signal has been observed [121], be it as a spectral bump expected for a focused pair beam or as an angularly extended signal should the electrons and positrons be deflected by the IGMF. Current GeV-TeV constraints have thus been interpreted as placing a lower limit on the IGMF strength, at the level of 10 −15 G for a coherence length larger than 10 kpc. 5 It should be noted that the limits imposed on the magnetic-field strength display order-of-magnitude variations when changing the source intrinsic parameters, in particular the duration of the active emission from the blazar, as well as the jet opening angle and orientation with respect to the line of sight. Despite being model-dependent, such gamma-ray constraints provide inputs to theoretical studies of the origin of cosmic magnetism [50,87], which could arise either from astrophysical outflows or from a primordial phase transition. The IGMF strength expected from the latter scenarios is furthermore bound on the upper side by the CMB spectrum and anisotropies, at the level of 10 −11 G for the most recent constraints [124]. More general limits, which are independent from the origin of cosmic magnetism, are placed on the maximum strength of the IGMF with rotation measures, at the level of 10 −9 G [125].
To conclude this section, while the pair-production process shown in the upper insert of Figure 2 can nowadays be claimed as one of the observational discoveries of gamma-ray astronomy, the fate of the pairs in the lower insert remains elusive. No definite conclusion on the main dissipation process of electron-positron beams has been reached on theoretical grounds. The detection of a secondary gamma-ray signal at GeV energies would provide crucial inputs to IGMF models and could guide the identification of the seeds of cosmic magnetism in large-scale structures.

Deviations from the standard model at ultra-high and ultra-low masses
We have seen in Section 2 how processes beyond the standard model can alter cosmological gamma-ray propagation. These processes involve either the oscillation of photons in axions and ALPs (with masses below m a 10 −6 eV) in astrophysical magnetic fields or LIV. The latter is expected at energy scales around the Planck energy E Pl . In this section, we review how gamma-ray observations can put constraints on axions, ALPs, and LIV. The constrained mass (or energy) scales are shown in Figure 4. For completeness, although not related to gamma-ray propagation, the masses of hypothetical weakly interacting massive particles (WIMPs) are also shown as they could self-annihilate (or decay) leading to the production of gamma-rays. We refer the reader to Ref. [129] for a recent overview of how ground-based gamma-ray telescopes can probe WIMPs, ALPs, and other beyond-the-standard-model phenomena.
We first turn to the searches for axions and ALPs. The interest in these particles was spurred by a tentative observation with the PVLAS experiment of an optical rotation generated . Energy and mass scales of physics beyond the standard model that can be probed with gamma-ray observations with satellites (blue boxes), ground-based telescopes (red boxes), or both (purple boxes). The range over which axions and ALPs are constrained lies below 10 eV. The mass range for self-annihilating WIMPs, between 10 6 and 10 14 eV, could be probed down to the thermal cross-section with future CTA observations of the Galactic center (Gal. Cen. [126]) or Fermi-LAT observations of dwarf spheroidal galaxies (dSphs [127]). LIV is probed up to 10 28 eV (n = 1) with both temporal and spectral observations as labelled in the Figure. The reader is referred to the repository in Ref. [128] and Figure 1 in Ref. [66] for more specific constraints on the parameter spaces of ALPs and LIV, respectively.
in vacuum by a magnetic field, which could be interpreted as evidence for ALPs [130]. 6 It was soon realized that ALPs over a broad range of masses, m a , and photon couplings, g aγ , would also affect the propagation of gamma-rays from extragalactic sources [138,139]. Since then, mainly two observational signatures have been identified: (a) spectral features-an energy dependent dimming of the flux in the form of (chaotic) oscillations-which occur around the critical energies given in Eqs. (7) and (8); (b) an enhancement of the gamma-ray flux in comparison to the case of pure absorption on background radiation fields. This enhancement should occur in the strong mixing regime, i.e. for gamma-ray energies E low crit < E < E high crit . Various astrophysical magnetic fields along the line of sight have been studied in this context. Starting from the production site of gamma-rays within AGNs, several authors considered photon-ALP conversion in the magnetic field of the jet, the termination lobes of the jet and the host galaxy [62,[139][140][141][142][143][144][145]. In particular, a coherent toroidal component of the jet magnetic field, which can be of the order of O(G) at pc distances from the central black hole (e.g., Ref. [146]), could lead to spectral features in the gamma-ray range, as can be inferred from Eq. (8). The region of ALP masses that can be probed through the conversion in these magnetic fields is roughly 1 neV m a 100 neV (region labeled "AGN jet" in Figure 4). However, the exact morphology of the jet magnetic fields is still unknown, as is the location of the gamma-ray production site along the jet. Ideally these parameters are constrained through multi-wavelength observations and / or are left as additional nuisance parameters in an ALP hypothesis test [145]. Radiation fields could also have a significant contribution to E high crit in particular in FSRQs [147].
Following the line of sight of the photon-ALP beam, the turbulent magnetic fields of galaxy clusters can act as efficient photon-ALP converters [143,148]. The fields could reach tens of µG in cluster centers [149] and induce chaotic oscillations in gamma-ray spectra around E low crit . Observations with ground-based telescopes and Fermi-LAT of bright gammaray emitting AGN in galaxy clusters provide world-leading constraints on the photon-ALP coupling for ALP masses around 1 neV m a 1 µeV (the region labeled "Gal. Cl.", for Galaxy Clusters, in Figure 4 [150][151][152][153][154][155]). One of the most studied sources in this regard is NGC 1275 at the center of the Perseus cluster. It should be noted though that the exact strength and configuration of the regular components of the large-scale magnetic field could change the ALP bounds considerably [156].
After leaving the source and a potential galaxy cluster, photon-ALP conversion could also occur in the IGMF (see Section 3 for a discussion of gamma-ray constraints on this magnetic field). This scenario has the appeal that the IGMF should be ubiquitously present in cosmic voids and, thus, one could expect photon-ALP conversion signatures in all AGN spectra. Using a cell-like turbulent IGMF with a field strength of the order of 1 nG, which is close to the maximum value allowed by rotation measures [125], it was realized that photon-ALP mixing could lead to an increased transparency of the Universe [157][158][159][160]. A more realistic IGMF modeling, based on large-scale cosmological simulations of the magnetic field evolution, led to the same result [161][162][163]. 7 Interestingly, evidence for such a reduced opacity has been suggested by several authors and interpreted as evidence for ALPs [148,157,[165][166][167][168][169][170][171]. This conclusion was reached mainly by observing an increasing gamma-ray flux in blazars at high energies and redshifts (and correspondingly large values of τ γγ ) when correcting the spectra for EBL absorption, which is not expected in standard blazar emission scenarios. Lower limits on the photon-ALP coupling to explain these evidences were derived in Ref. [172]. In contrast, other studies have confirmed that standard EBL absorption describes the observed spectra well when large samples of spectra from ground-based telescopes and satellites are used [102,173] and when the response of the instrument, in particular the energy resolution, is taken into account [174]. Furthermore, a recent analysis of the highest energy photons recorded by the Fermi-LAT showed no preference for the scenario of photon-ALP mixing in the IGMF [175].
Lastly, the photon-ALP beam enters the Milky Way. Its magnetic field can be modeled as a coherent and turbulent component (e.g., Ref. [176]). The coherence length for the turbulent component is too small compared to λ osc to cause any significant oscillations. In contrast, depending on the line of sight, the coherent field can be of the order of O(µG) over kpc length scales leading to a potentially large oscillation probability [133,134,140,148]. The lineof-sight dependence through the Galactic magnetic field (GMF) could lead to a correlated detection of distant blazars (as ALPs reconvert back in the GMF to gamma-rays) [177] and an anisotropy in the diffuse gamma-ray background [178,179]. Furthermore, several authors have searched for ALP-induced spectral features around E low crit in Galactic sources [180][181][182]. A preference for the existence of ALPs was for example claimed in the analysis of pulsar spectra measured with the Fermi-LAT [180]. Such a claim is, however, at odds with results from the CAST experiment, which searches for ALPs produced in the Sun [183]. 8 Interestingly, the photon-ALP conversion in the GMF could also enable the detection of ALP-induced bursts of gamma-rays from core-collapse supernovae, which should arrive simultaneously with the neutrinos produced in the collapse. The non-detection of such a burst from SN 1987A has lead to stringent bounds on ALPs with masses below m a 1 neV [133,134,184]. ALP production in the supernova population has similarly been used to place constraints on the photon-ALP coupling [185]. These constraints would disfavor the ALP-induced reduction of gamma-ray opacity discussed above. 9 Heavier axions and ALPs could also be probed indirectly by gamma-ray observations. The decay time of axions and ALPs into two photons is (e.g., Ref. [187]) τ aγγ = 64πg −2 aγ m −3 a ≈ 1.3 × 10 27 s (g aγ /10 −11 GeV −1 ) −2 (m a /eV) −3 , which is ten orders of magnitude larger than the present age of the Universe. Decay of light ALPs discussed so far is thus irrelevant to observational probes. However, dark-matter axions and ALPs with masses around 1 eV m a 10 eV could decay into optical and UV photons and contribute to the EBL photon density (region labelled "EBL" in Figure 4) [188]. This would in turn increase the opacity of the Universe. First analyses indicate that ground-based gamma-ray observations are promising to constrain these comparatively heavy ALPs [189][190][191].
Besides constraining the ALP parameter space, the propagation of gamma-rays provides constraints on LIV (see, e.g., Refs. [66,67,88] for recent reviews). Similarly to hints for a reduced opacity discussed above for ALPs, the study of modifications of the pair-production threshold induced by LIV gained traction when the first blazars, namely Mrk 421 and Mrk 501, were observed in flaring states at energies up to ∼ 20 TeV (e.g., Refs. [192,193] for historical references). Contrarily to ALP spectral modifications, the LIV effect is expected above a fixed energy, around 20 TeV for n = 1 and M = E Pl . Such an effect emerges from the modified dispersion relation, to order n, of the particles involved in pair production, which has been treated either by deriving the modified threshold energy [194] or by considering an effective mass for photons [195]. Both approaches are purely based on kinematics and assume no modification of the pair-production cross section, which can only be predicted in specific extensions of the standard model that break Lorentz invariance (see discussion in Ref. [196] and references therein).
Using state-of-the-art EBL models, there is currently no indication for a reduction of the Universe transparency to gamma-rays at very-high energies. Ground-based gammaray observations led to exclusions using searches in single observations [197] or combined analysis of multiple sources [102,198], where the latter ones provide the strongest constraints up to 1 − 10 E Pl (the probed range of M is labeled "Blazar spectra" in Figure 4). Interestingly, constraints on the LIV energy scale could also be derived from the non-observation of ultrahigh-energy photons at E > 10 18 eV with the Pierre Auger Observatory [199]. For an unaltered pair-production threshold, a photon flux at these energies would be strongly absorbed in the interaction with CMB and radio photons. Not unlike gamma-ray constraints on the IGMF, ultra-high-energy constraints on LIV depend on the source properties which remain a limiting unknown.
The other channel through which gamma-ray propagation can be used as a probe of LIV lies in the study of time delays as a function of energy and distance [68]. Ground-based observations of blazar flares have led to constraints on the LIV scale (box labeled "∆t Blazars" in Figure 4) to be above 0.2 E Pl [197,200,201]. Ground-based gamma-ray telescopes have also recently discovered TeV emission from GRBs [202][203][204]. In particular, the detection by MAGIC of GRB 190114C (z = 0.42) at energies above 0.2 TeV ruled out M < 0.5 E Pl [205]. The most stringent time-delay constraints come from Fermi-LAT observation of GRB 090510 at a redshift of z ∼ 0.90 which rules out M < 2 E Pl [206]. The values of M probed with GRBs are shown in Figure 4 in the region labelled "∆t GRBs". It should be noted that even though these time-delay analyses achieve strong constraints on the LIV scale M, assumptions have to made about the intrinsic time spread of the emission. This problem is starting to be addressed both for AGN flares [207] and GRBs [208,209] in a systematic way, but further work appears to be needed to disentangle intrinsic from possible line-of-sight effects [210].

Knowns, expectations and hopes from gamma-ray cosmology
The exploration of the extragalactic gamma-ray sky around TeV energies has unveiled some of the most extreme accelerators in the Universe: from blazars detected up to 20 TeV at light travel times of ∼ 500 Myrs (z ≈ 0.03) out to distant ones detected above 200 GeV at ∼ 8 Gyrs (z ≈ 1). Gamma-ray observations of jetted AGNs at very-high energies provide powerful probes of the content of intergalactic medium, as discussed in Section 3, and of physics beyond the standard model, as discussed in Section 4. The recent discovery of emission beyond 100 GeV from long GRBs also provides enticing prospects to complement the explorations performed with AGNs.
The long-sought signal of gamma-ray absorption on the EBL, predicted in the 1960s, has been uncovered by observations from both gamma-ray satellites and ground-based telescopes. Measurements of the EBL imprint thus far primarily constrain the present-day COB at optical and near-infrared wavelengths, with a precision better than 30%. Gamma-ray astronomy at intermediate and high redshifts (z > 1) from space has further provided its first competitive constraints on the EBL evolution with redshift, which is dictated by the cosmic history of star formation, dust consumption and AGN formation. The study of the baryonic and light contents of the ionized Universe, and of their co-evolution, is expected to strongly benefit from continued operation of Fermi-LAT in space during the early-science phase of the new-generation gamma-ray observatory, the Cherenkov Telescope Array (CTA [211]). Through its key-science programs, CTA will be able to probe the EBL evolution out to at least z = 2 with a precision two-to-three times better than current-generation telescopes [120]. CTA will in particular probe with unprecedented accuracy the diffuse or unresolved emissions in the region of the "optical controversy", in conjunction with deepand wide-field observations from JWST, the Rubin Observatory, Euclid and the Roman Space Telescope [211]. The possible identification of distant (re-)ionizing sources with dedicated optical and near-infrared observations together with CTA constraints on emission from the overall population(s) can let us hope for major steps forward in constraining the low-end of the galaxy luminosity function and in the understanding of the intergalactic medium. At higher EBL wavelengths, around 100 µm, CTA observations of high emission states from nearby extragalactic sources are expected to place tight constraints on the CIB, not only assessing the cosmic-dust emissivity but also resolving an important source of uncertainties in the propagation of ultra-high-energy nuclei (E > 10 18 eV) as they photodissociate on this infrared background [212]. As ultra-high-energy nuclei are charged, their propagation is also affected by magnetic fields on galactic, intrahalo, intracluster and intergalactic scales, which could impact the interpretation of anisotropies detected on large-angular scales at E > 8 × 10 18 eV [213] and suggested on intermediate angular scales at E > 40 × 10 18 eV [214]. New insights on cosmic magnetism are expected from radio observations of clusters and filaments with SKA and LOFAR (see, e.g., the MAGCOW project [215]). Such observations could guide numerical simulations of large-scale structures and inform gamma-ray searches for IGMF signatures. CTA in particular offers the potential to probe the IGMF in voids up to strengths of at least 0.3 pG [120]. A strong interplay of theoretical/numerical advances and astroparticle/astronomical observations offers exciting prospects in our understanding of the intergalactic medium in the upcoming decade.
As we have seen in Section 4, searches for physics beyond the model affecting gammaray propagation have also become a vibrant field of study. CTA will in particular have unprecedented sensitivity searching for spectral irregularities [120] and a reduced gamma-ray opacity (e.g., Refs. [162,216]) as induced by ALPs. With the foreseen energy range of CTA, it will be possible to probe higher ALP masses potentially up to hundreds of neV for couplings 10 −12 GeV −1 for cluster magnetic fields of the order of 10 µG [120]. At these masses and couplings, CTA observations will begin to probe the parameter space where ALPs could constitute the entire cold dark matter content of the Universe. A similar region of the ALP parameter space could be probed across the energy ranges of the TAIGA detector [217] and the LHAASO water Cherenkov detector [218]. For high values of g aγ (or in optimistic magnetic field scenarios), LHAASO could further detect reconverted ALPs at the high-energy end of the diffuse gamma-ray background [219]. Space-based gamma-ray observatories could also detect ALP-induced bursts of gamma-rays from supernovae, either in the Milky Way [161] or in close-by galaxies [220]. Additionally, theoretical advancements open up the parameter space of the QCD axion towards lower masses which could be probed with gamma-ray instruments (e.g., Refs. [221][222][223]). The development of open-source software [224,225] facilitates the computation of photon-ALP oscillation probabilities in various magnetic field environments for young scientists entering the domain. Possible evidences or constraints derived from future gamma-ray observations should be easily testable with the next-generation of axion experiments, also coming online in the next decade (e.g., Ref. [226]). In that sense, the search for axions and ALPs with gamma-ray observations nicely complements the experimental efforts to search for this kind of particle.
Searches for LIV will greatly benefit from deep-field observations with CTA as well as from wide-field observations beyond tens of TeV with LHAASO and SWGO. In light of the recent transient events detected with already operating ground-based telescopes, detections of GRBs with CTA appear to be guaranteed, presumably at a rate larger than 0.1 − 1 per year that was initially suggested [227][228][229]. While LHAASO and SWGO will continuously monitor a large fraction of the sky and possibly detect some serendipitous events [230], CTA will have to follow-up on external triggers provided by multi-wavelength facilities such as Swift, SVOM or Fermi-GBM [211]. For the linear subluminal LIV scenario, the sensitivity of CTA could reach 10 E Pl depending on the time delay between the GRB trigger and the start of the observation [229]. LHAASO could even reach a sensitivity to LIV of the order of ∼ 20 E Pl [218]. These limits will of course depend on a better understanding of the intrinsic time lags, which need to be incorporated in the data analysis [209,210]. Furthermore, optical follow-up observations will be of paramount importance in order to determine the redshift and hence the distance of GRBs. CTA observations of blazars should also improve considerably the sensitivity to alterations of the pair-production threshold. The authors of Ref. [120] used simulated CTA observations of a flare from Mkn 501 as well as long-term observations of the hard-spectrum blazar 1ES 0229+200 to estimate the possible constraints on M in the linear and quadratic LIV scenario. With single observations of these sources, energy scales M ∼ E Pl could be achieved in the linear case and a source stacking should push the sensitivity a factor two-to-three beyond current constraints.
To conclude, the study of gamma-ray propagation has placed competitive constraints on cosmological entities and on fundamental physics over the past decade. Such constraints and discoveries have been enabled by the capablities of instruments such as Fermi-LAT, H.E.S.S, MAGIC and VERITAS, which have left us with challenging questions: Have we resolved only half of the optical light in the Universe? Where is cosmic magnetism stemming from? Could dark matter consist in light particles such as axions? Is Lorentz symmetry broken beyond the Planck scale? We look forward to CTA and other upcoming gamma-ray observatories to help the multi-wavelength and multi-messenger communities take up the challenge.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

1EG
First