Constraints on the intergalactic magnetic field strength from $\gamma$-ray observations of GRB 221009A

Characteristics of the cascade gamma-ray signal resulting from very-high-energy gamma-ray sources, such as gamma-ray bursts, can be used to constrain the strength and structure of intergalactic magnetic fields (IGMF). There has been a debate on whether GRB 190114C, the first gamma-ray burst with observed TeV photons, can constrain the IGMF. Recently, LHAASO detected the brightest-of-all-time GRB 221009A, which has much larger energy in TeV band and the spectrum extends to energy above 10 TeV, providing an unprecedented opportunity to studying IGMF. We perform a Monte-Carlo simulation of the cascade process with the public ELMAG code, considering the TeV data of GRB 221009A observed by LHAASO. By comparing the resulting cascade emission with the flux limit obtained from Fermi-LAT observations, we infer a limit of $B\ge 10^{-18.5}\rm G$ for IGMF. Though this limit may not be as strong as the limit from blazars, it serves as an independent constraint on IGMF from a new class of TeV sources.


INTRODUCTION
The magnetic fields in galaxies and galaxy clusters are thought to result from the amplification of seed magnetic fields, which might exist in their initial form in the intergalactic medium (Neronov & Vovk 2010).There are two broad classes of models for their origin of the seed magnetic fields (Durrer & Neronov 2013): (1) cosmological models, in which the seed fields are generated in the early universe before the structure formation; (2) astrophysical models, in which the seed fields are produced during the epochs accompanying the gravitational collapse leading to structure formation.Measurements of the strength of intergalactic magnetic fields (IGMF) can provide an important clue on the origin of the initial seed fields.
Very high energy (VHE, ≥100 GeV) transient sources such as flaring active galactic nuclei (AGNs) and gamma-ray bursts (GRBs) are viable tools to constrain IGMF (Plaga 1995;Dai & Lu 2002;Wang et al. 2004;Razzaque et al. 2004;Ichiki et al. 2008).During propagating in the universe, TeV photons emitting from these sources will interact with Extragalactic Background Light (EBL) and produce e + e − pairs via pair production process.The created pairs are deflected by the IGMF and radiate secondary GeV-TeV emission through inverse Compton scattering (ICS) off Cosmic Microwave Background (CMB) photons.When these cascade photons arrive at Earth, their properties, such as the spectrum and the time delay with respect to the primary emission, carry critical information of the IGMF, and can be thus used as a diagnosis of the IGMF (Dai & Lu 2002;Wang et al. 2004;Dermer et al. 2011;Taylor et al. 2011).
This method was first applied to blazars (Neronov & Vovk 2010), which give a lower bound on the IGMF at the lever of B≥ 10 −16 G from the non-detection of GeV gamma-ray emission from electromagnetic cascade initiated by the primary TeV gamma-rays in intergalactic medium.However, the constraints were subsequently found to be subject to significant systematic effects, such as the unknown duty cycle (Dermer et al. 2011;Taylor et al. 2011), the poorly constrained spectral properties of the source, and uncertainties in the EBL spectrum (Arlen et al. 2014;Finke et al. 2015).For example, restricting TeV activity of 1ES 0229+200 to a timescale of 3-4 years during which the source has been observed leads to a more robust lower limit of B≥ 10 −18 G (Dermer et al. 2011).Under such circumstances, GRBs, as an independent class of transient TeV sources, becomes crucial for constraining IGMF.
As a short-lived powerful TeV source, GRB is suitable for studying IGMF because the time-delayed cascade photons can be easily distinguished from the primary photons.GRB 190114C is the first GRB observed with TeV emission and a limit of B ≥ 10 −19.5 G on IGMF (for the coherence length of λ ≤ 1 Mpc) was inferred (Wang et al. 2020).On the other hand, Dzhatdoev et al. (2020) used Monte-Carlo code ELMAG to calculate the cascade emission for various EBL models, and found that the sensitivity of Fermi/LAT is not sufficient to constrain the IGMF.The discrepancy could be due to that Dzhatdoev et al. (2020) only take into account the primary TeV photons during the period from 62 s to 2454 s, while Wang et al. (2020) consider the power-law decay of the afterglow flux starting from 6 s (MAGIC Collaboration et al. 2019;Ravasio et al. 2019;Wang et al. 2019), which leads to a difference in the fluence of the primary TeV photons by a factor of 5. Recently, LHAASO detected the brightest-of-all-time GRB 221009A(LHAASO Collaboration 2023), which has much higher fluence in TeV band and the spectrum extends to energy above 10 TeV (Huang et al. 2022), thus it offers a precious opportunity for us to constrain IGMF.
LHAASO detected more than 60,000 photons with energies greater than ∼ 200 GeV from this GRB(LHAASO Collaboration 2023).The observed spectra at various time intervals show sharp steepening at high energies due to the EBL absorption.The intrinsic spectra after correcting for EBL absorption can be described by a single power-law extending to the highest observed energy.By integrating the time-resolved intrinsic spectra, the isotropic equivalent energy in TeV is E (0.3−15TeV) ≈ 1 × 10 53 erg.In addition, the clear steepening of the flux observed by LHAASO in the decay phase, which is consistent with a jet break, provides us the information regarding the half-opening angle (θ j = 0.8 • ).
The rest of the paper is organized as follows.In §2, we describe our analysis of Fermi-LAT data in order to set experimental upper limits on the pair echo intensity from GRB 221009A.In §3, we use the open code EL-MAG to study the electromagnetic cascades of primary TeV photons and obtain constraints on IGMF.Finally we give conclusions in §4.
The Fermi-LAT extended type data for the GRB 221009A are taken from the Fermi Science Support Center1 .Only the data within 30 We perform a binned maximum likelihood analysis for this GRB, and considering the LAT SOU RCE events between 100 MeV and 1 TeV.The corresponding instrument response function (IRF) (P 8R3 SOU RCE V 3)2 is used.A maximum zenith angle of 90 • is adopted to reduce the contamination from the γ-ray Earth limb.The LAT 12-year Source Catalog (4FGL-DR3) sources are included in our analysis.For the main background component, we consider the isotropic emission template ("iso P 8R3 SOU RCE V 3 v1.txt")and the diffuse Galactic interstellar emission template (IEM; gll iem v07.f its) in our analysis.The parameters of isotropic emission and IEM are left free.
Assuming a power-law spectrum (dN/dE = AE Γ ) for this burst, we obtained its spectral energy distribution (SEDs) of different time intervals with software f ermipy (version v1.1) (Wood et al. 2017).For the energy bins with a test statistic (TS3 ) value < 9 (TS=9 corresponding to 3σ significance), we extract the upper limits at 95% confidence level.The SEDs of different time intervals are shown as the black points in Figures 1, 2and 3.

Monte Carlo Simulation
Very high energy photons are absorbed by EBL when they are propagating in the IGM, creating electrons and positrons through pair production.These pairs then scatter CMB photons to the GeV domain via inverse Compton (IC) radiation.Concurrently, the IGMF de- flects the pairs, causing these secondary GeV photons to reach the observer with different directions and arrive later compared to the TeV photons.The characteristics of these secondary GeV photons, in terms of their duration and strength, can be thus used to constrain the IGMF.It is worth noting that the angular spread of pair production and IC emission induces an intrinsic time delay even in the absence of the intervening magnetic fields.As estimated by Vovk (2023) for GRB 190114C, this time spread can span a range from 10 3 to 10 5 seconds depending on the energy.However, in this study, the effect of the intrinsic time spread is negligible due to the much longer time intervals used in this study.
The open code ELMAG (Kachelrieß et al. 2012) is a Monte Carlo simulation program designed to study electromagnetic cascades on the EBL, including the deflections of charged particles in IGMF.Due to the introduction of some additional features, this program provides an accurate description of the particle trajectories.These features include the turbulence of extragalactic magnetic fields, the opening angle of jet and the calculation of three-dimensional trajectories of the secondary electrons and positrons by solving the Lorentz force equation.In this work, we use the newest edition ELMAG 3.03 (Blytt et al. 2020) to perform a full threedimensional simulation.
We choose three time intervals: T 0 +4000 seconds to 3 days, T 0 +3 days to 30 days and T 0 +30 days to 8 months to probe an early stage and two longer observation windows.Here, the 4000 s is the moment when the Fermi-LAT detector entered the field of view for the second cycle.Furthermore, the 3rd day represents the transition time from Fermi-LAT detection to non-detection of GeV photons (Stern & Tkachev 2023).Moreover, to obtain a more stringent constraint on the IGMF, we set the observation limit at 30 days, considering that the theoretical cascade flux diminishes with time (i.e., decreases as t −1 ) while the Fermi-LAT upper limit flux follows t −1/2 for long-term observations.The final time bin extending to 8 months corresponds to the data accumulated as of writing of this paper.
We use a turbulent magnetic field with a Kolmogorov spectrum4 , where the minimum and maximum spatial scales are set as the default values L min = 5 × 10 −4 Mpc and L max = 5 Mpc, respectively.A coherence length of 1 Mpc is assumed, consistent with the analytic calculation.
Due to the potential systematic effect existing in different EBL models, we follow the approach described in LHAASO Collaboration (2023) to take three EBL models into consideration, namely, the models of weak attenuation, standard situation and strong attenuation in the Saldana-Lopez et al. ( 2021) model.The corresponding intrinsic spectra, after the correction for these EBL models, have been obtained in LHAASO Collaboration (2023).To simplify the calculation, the spectral index is adopted from the value during the time interval when the TeV flux is the highest 5 .Then we perform the simulations for these three EBL models using the input parameters shown in Table .1.In our calculation, we treat the primary emission as an instantaneous injection to simplify the calculation.But for the first time interval, considering that the duration of primary TeV emission is comparable to this time interval, such simplification is less accurate.Therefore we calculate the start time for each time-resolved spectrum of the primary emission (LHAASO Collaboration 2023) and integrate their cascade emission over the time interval.The results are displayed in Figure 1.Notably, the lower and upper ranges of the color bands in Figure 1 represent the results obtained from the low and high EBL intensity models, respectively.For IGMF ranging from B = 10 −20 G to B = 10 −18.5 G, the theoretical flux significantly exceeds the Fermi-LAT upper limits in some energy ranges, while for the case of B = 10 −18 G (as shown in the bottom panel of Figure 1), the theoretical flux is below the upper limits in the low EBL model case.
We thus conclude that a conservative lower limit of the IGMF is B ≥ 10 −18.5 G.This limit on IGMF is much more stringent than that derived from GRB 190114C (Wang et al. 2020;Dzhatdoev et al. 2020).

Constraining IGMF using analytic approach
For a cross-check, we perform an analytic calculation of the cascade emission and then derive the constraint on IGMF.A comprehensive description of the analytic approach is provided in the Appendix.We find that the constraint on IGMF through the analytic process is ≥ 10 −18 Gauss.The analytic approach does not take into account of the uncertainties in EBL models.Ad-5 The effect of the the uncertainties in the TeV flux normalization and slope and their time variation on the cascade flux is found to be unimportant as it is smaller than that caused by the uncertainty in EBL models.ditionally, we explore the case where the maximum energy is 7 TeV, which is the maximum energy reported by WCDA observations.We find that the constraint on IGMF is almost the same, as shown in the appendix Figure 3.
3.3.On the origin of the 398 GeV photon Xia et al. (2022) proposed that a 398 GeV photon arriving at 0.4 day could come from the cascade emission.Under this hypothesis, it requires the strength of IGMF to be ∼ 10 −17 G.In order to verify this scenario, we adopt the precise parameters of the TeV emission measured by LHAASO to study the possibility.We use ELMAG to calculate the cascade flux within the energy range of 200 to 800 GeV over a period of from 0.2 to 0.8 days after the burst for three magnetic fields (10 −18 G,10 −17 G and 10 −16 G).For simplicity, we only consider the standard case of the EBL model in Saldana-Lopez et al. (2021).Figure .2presents a comparison between the theoretical flux of the cascade emission and the observed data by Fermi-LAT.As we can see, the theoretical flux is much lower than the observed one at the energy around 398 GeV.By considering the corresponding exposure time and effective area of LAT during this event, we estimate the expected photon number.The Poisson probabilities for detecting one such photon are 0.65%, 0.82%, 0.17% for the three magnetic field cases, respectively.Although the probability for such an event is low, we cannot conclusively rule out the scenario at a confidence level of ≳ 3σ.Nevertheless, should this scenario hold, it would require a magnetic field of approximately 10 −17 G, which does not conflict with our lower boundary on IGMF.GRB 221009A, as a once-in-ten-thousand-year event (Burns et al. 2023), has the largest amount of energy in the TeV band and its spectrum extends to at least 10 TeV (Huang et al. 2022).In this work, we obtain the constraints on the IGMF based on the LHAASO and Fermi-LAT observations of GRB 221009A.We used a Monte-Carlo code ELMAG 3.03 to calculate the cascade emission.The constraints on magnetic field from the analytic calculation are consistent with that from the Monte-Carlo code.The Monte-Carlo approach takes into account the uncertainty in the EBL intensity and we obtain a conservative lower limit, which is B ≥ 10 −18.5 G for the coherence length of λ ≤1 Mpc.This result is much more stringent than that derived from GRB 190114C (Wang et al. 2020;Dzhatdoev et al. 2020).While this limit may be not as strong as that from blazars (see the recent work Acciari et al. (2023), which obtained a limit of ≥ 10 −17 G ), it serves as an independent constraint on IGMF from a new class of TeV sources.
One significant systematic effect originates from uncertainties in the EBL model.Following the treatment in LHAASO Collaboration (2023), we adopted two additional EBL intensity models to testify the systematic error they may bring, corresponding to the lower (weak attenuation) and upper (strong attenuation) boundary of the error range in the Saldana-Lopez et al. ( 2021) model.Most of other EBL models are compatible with this error band.We would like to point out that our limits are derived from the lowest EBL intensity in the model of Saldana-Lopez et al. (2021), which represnts a conservative constraint on IGMF.
Another parameter affecting the constraints on the IGMF strength is the coherence length, λ coh .In this work, we set the coherence length as 1 Mpc (approximately equal to the cooling length of electrons λ IC which upscatter CMB photons to 1 GeV).As indicated by the equation of θ B (λ coh ≤ λ IC ) in the Appendix, the constraint will become more stringent as λ coh decreases, scaling as λ −1/2 coh (Neronov & Vovk 2010), while the IGMF constraint will be independent of coherence length if λ coh > 1Mpc.
We thank the anonymous referee for valuable suggestions, and we thank Ievgen Vovk for helpful discussions.This assumes that all photons that were emitted by electron pairs during the IC cooling time t IC have been received by observer for a typical delay time ∆t B .
However, a more sophisticated treatment should consider the effect of cooling on the the energy spectrum of pairs, which leads to a shift in population of electrons from high energy to low energy.In this study, we take into account the evolution of electron distribution (∂N (γ e , t)/∂t = −∂[ γe N (γ e , t)]/∂γ e ).Consequently the distribution of electrons and time delay become time-dependent (dN e,0 (ϵ)/dγ e → dN e,t (ϵ)/dγ e ; ∆t B (ϵ, γ e , B) → ∆t B (ϵ, γ e , B, t)).As a result, the electron distribution is modified to be dN dγ e = dt dϵ dN e,t (ϵ) dγ e 1 max (∆t B (ϵ, γ e , B, t) , ∆t obs ) (A5) and the spectrum of IC photons is given by where f (x) = 2x ln x + x + 1 − 2x 2 and x = ϵ IC /4γ 2 e ϵ 0 .Additionally, it is worthy noting that the resulting photon flux represents an average effect over the entire observation period ∆t obs , which spans from t start to t end .If the time delay of the secondary photon emission is less than the start time of the observation period (i.e., ∆t B (ϵ, γ e , B, t) < t start ), then all corresponding secondary photons have already arrived, so we remove the contribution from these particles in the calculation.

A.2. Limits on IGMF using GRB 221009A
Here, we consider four cases of the IGMF strength, ranging from 10 −20 to 10 −18 Guass.With the analytic description provided, we have better understanding on the evolution trend of the cascade emission spectrum for varying magnetic field strengths as discussed in Section 3.According to Eq.A3, the arrival time is proportional to magnetic field strength and inversely proportional to the energy of primary photons.As a result, the peak energy of the cascade spectra varies systematically with the magnetic field strength and observing time, as shown in Fig. 3.The favorable choice for studying the magnetic strengths is to ensure that the peak energy falls within the Fermi-LAT energy window.Specifically, choosing a relatively later time interval can provide a better constraint for the high IGMF cases.
In our calculation, the primary emission is treated as instantaneous injection.We assume that the shape of intrinsic spectrum is a power-law (i.e., dN/dE = A(E/TeV) −Γ ), which is consistent with the results reported in LHAASO Collaboration (2023).The normalization factor (A = 1.56 × 10 −7 TeV −1 cm −2 s −1 ) is the average value of the timeresolved spectra over the full TeV emission duration 1749 s, and the spectral index (Γ = 2.36) is adopted from the value during the time interval when the TeV flux is dominant (from T 0 + 248 s to T 0 + 326 s).As there was no clear cutoff or break in the intrinsic spectrum below 7 TeV (LHAASO Collaboration 2023), we set the maximum energy at 15 TeV, which we believe to be a reasonable choice.The coherence length of IGMF is set to be 1 Mpc. Figure 3 presents a comparison between the expected spectrum of the echo emission and the upper limit imposed by the Fermi/LAT observations.We can find that the predicted flux exceeds the upper limits imposed by Fermi/LAT for IGMF with B ≤ 10 −18 G, so we obtain a limit on IGMF of B ≥ 10 −18 G for our analytic approach.

Figure 1 .
Figure 1.Comparison between the expected flux of the echo emission and the observed data by Fermi/LAT in Monte-Carlo simulation.The dotted lines represent the echo emission using the standard flux of the EBL model of Saldana-Lopez et al. (2021), while the color bands represent the echo emission considering the uncertainties in the EBL intensity.

Figure 2 .
Figure 2. Comparison between the expected flux of the echo emission and the observed data by Fermi/LAT in the time interval of T0 + 0.2-T0 + 0.8 day, during which the 397 GeV photon was detected.
The work is supported by the National Key R&D Program of China under grant No. 2022YFF0711404, the NSFC Grants No.12121003, No. 12203022 and No. U2031105, and China Manned Spaced Project (CMS-CSST-2021-B11).

Figure 3 .
Figure 3.Comparison between the expected flux of the echo emission and the observed data by Fermi/LAT in the analytic approach.The SEDs of the echo emission are averaged over the observation time 4000s-3 days (top panel), 3-30 days (middle panel) and 1-8 months (bottom panel).The dashed and solid lines represent the cases assuming a maximum energy of 7 TeV and 15 TeV, respectively.The black points denote the Fermi-LAT detection while the black arrows denote the upper limits.

Table 1 .
Input parameters used in the Monte-Carlo simulation for three EBL model.The intrinsic spectrum is assuming as a power-law function (i.e., dN/dE = A(E/TeV) −Γ ).