Absorption of Gamma-Ray Burst X-Ray Afterglows by the Missing Baryons: Confronting Observations with Cosmological Simulations

Anomalous absorption of the X-ray afterglow of gamma-ray bursts (GRBs) is prevalent and has been known since the days of BeppoSAX (Stratta et al. 2004). Beyond the expected Galactic absorption, there is an excess column density of ∼10²¹ cm⁻²–10²³ cm⁻² that remains to be explained. This excess absorption is detected through attenuation of the soft X-rays that increase toward low energies and is measured down to ∼0.3 keV. Such attenuation is typical of photoelectric (ionization) absorption. Identifying the origin of this excess absorption is challenging, because the spectra are void of any discrete features, such as absorption lines or edges. The common practice is to assume an absorber that is (a) neutral, (b) with solar abundances, and (c) at the redshift z of the GRB, and to quote an equivalent column density N_H. Lacking any supporting evidence in the data for these assumptions, one needs to keep in mind that the actual column density of the absorber depends strongly on all three of them and may be significantly different. For example, partial ionization or subsolar abundances will increase the derived N_H.

The X-ray telescope (XRT) on board the Swift observatory (Gehrels 2004) considerably increased the sample size of GRB afterglow spectra. The equivalent N_H increases dramatically with z, from ∼10²¹ cm⁻² at low z up to ∼10²³ cm⁻² at z > 5 (Campana et al. 2010). This result holds regardless of the Galactic absorption model (Willingale et al. 2013; Arcodia et al. 2016), the sample selection (Campana et al. 2012), or the variability of the X-ray afterglow (Valan et al. 2023). Rahin & Behar (2019) later found the same trend with a larger sample. This trend with z is curious, as no other evolution of GRBs with redshift is known. Additionally, the equivalent X-ray N_H can be orders of magnitude higher than the Lyα column density readily measured at the GRB host (Watson et al. 2007; Campana et al. 2010), which has no redshift evolution (Rahin & Behar 2019).

The directly measured quantity in GRB X-ray afterglows is not N_H, but the optical depth τ = σ N_H. The fact that for z ≳ 2, τ(z) tends to a constant τ(0.5 keV) ≈ 0.4, while N_H increases strongly with z, points to a varying absorption cross section (σ) with z, which is exactly what is expected from a photoionization cross section that decreases strongly with energy (z). This physical effect led Behar et al. (2011) to argue that the observed X-ray absorption may be due to the missing baryons in the intergalactic medium (IGM), while the host contribution to the 0.5 keV opacity decreases with z, becoming negligible at z > 2. Most of the IGM opacity accumulates at low z, eventually saturating at z > 3. Assuming a cosmologically uniform neutral absorber, and an approximate cosmic cross section of σ ∝ E^−2.5, they showed that the asymptotic X-ray optical depth at high z depends only on metallicity:

$$\begin{eqnarray}&&{\tau }_{\mathrm{IGM}}(z\to \infty )=\displaystyle \frac{2{Z}_{0}}{1+k},\end{eqnarray} \tag{ 1 }$$

where Z₀ is the metallicity at z = 0 in solar units and k represents the metallicity evolution $Z{(z)={Z}_{0}(1+z)}^{-k}$.

Starling et al. (2013) further accounted for ionization in the IGM at 10⁵–10^6.5 K and treated H and He separately from metals. They showed that the IGM could explain the observed X-ray opacity at z > 3, provided it is metal-enriched to ≈0.2 solar metallicity. However, they confirmed that intervening cold gas could also account for the observed opacity. Wang (2013) found a correlation between the N_H of GRB afterglows and that of IGM intervening systems. Dalton et al. (2021) confirmed the significant contribution of the IGM to N_H. The IGM interpretation is further supported by the lack of absorption variability despite drastic changes in the X-ray afterglow flux (Valan et al. 2023). Alternative explanations for the discrepancy between N_H and host opacity have been suggested (Watson 2011; Krongold & Prochaska 2013; Watson et al. 2013; Tanga et al. 2016; Heintz et al. 2018); however, none of these explain the z dependence of N_H.

The above IGM hypothesis begs for confrontation with cosmological simulations. Campana et al. (2015) used the RAMSES code (Teyssier 2002) to show that 20% of the metals in the universe can reside in the IGM and around galaxies and can explain the X-ray opacity toward distant X-ray sources. In this paper, we use the IllustrisTNG simulations (Nelson et al. 2018a) to compute X-ray opacities at cosmic distances. Illustris simulations include star formation, metal enrichment, and radiative feedback that were not available for Campana et al. (2015). In the following sections, we describe the simulations (Section 2), the computed opacities (Section 3), and the computational results (Section 4). We then compare these opacities and their z dependence with Swift GRB N_H measurements in Section 5, and draw our conclusions in Section 6.

We employ the publicly available IllustrisTNG simulation (Nelson et al. 2018a, 2019; Pillepich et al. 2019) to compute optical depths toward sources at different redshifts. We worked with TNG504, which evolves a cube of 350 comoving Mpc (cMpc h^–1) on the side. It has a dark mass particle resolution of 2.8 × 10⁸ M_⊙.

Higher-resolution and larger simulations are available, but require computing resources that are prohibitive for us. TNG504 has 100 snapshots; each snapshot captures a specific redshift. We use 96 of them between z = 0 and 10. Metallicity and neutral H fraction are available for only 18 of the 96 snapshots. These 18 snapshots are distributed all along the redshift range 0 ≤ z ≤ 10, in steps of approximately Δz = 0.1 for 0 ≤ z ≤ 1 and Δz = 1 for 1 < z ≤ 10. ¹

In the other snapshots, we linearly interpolated. For column density measurements, we collapse each snapshot to two dimensions. Six such collapsed snapshots are shown in Figure 1, demonstrating structure formation from z = 0–10.

Zoom In Zoom Out Reset image size

In order to compute the column densities and optical depths, we need to extract the number density from the simulations and integrate over cosmic sightlines. To obtain the average mass density ρ in each bin of the 2D-collapsed (64 × 64) grid (Figure 1), we take its total mass and divide by its comoving volume 350*(350/64)² ∼10⁵(cMpc h^–1)³. The physical number density of each species s (e.g., H, He, ...) is then

$$\begin{eqnarray}&&{n}_{s}=\displaystyle \frac{\rho {A}_{s}{\left(1+z\right)}^{3}}{{M}_{s}},\end{eqnarray} \tag{ 2 }$$

where A_s is the average abundance of the s species in the bin and M_s is the mass of that species. The factor of (1 + z)³ gives the physical density at the snapshot of redshift z. Fully ionized atoms do not contribute to the photoelectric optical depth. The ionization of H is readily given by the simulations and features a sharp ionization phase at z = 5–6, leaving only ∼1% of H neutral. The ionization of He and metals is not provided by the simulation. We assume that He is ionized once at z ≈ 4.5 (its first ionization potential is higher than H) and is fully ionized at z ≈ 2.5. This is roughly where observations of IGM absorption systems indicate He is fully ionized (Furlanetto & Oh 2008; Dixon & Furlanetto 2009; Upton Sanderbeck & Bird 2020). At z < 2.5, we assume that the low neutral He fraction follows that of H from the simulation. The neutral fractions of H and He used in the present computations are summarized in Figure 2. Metal ionization is neglected, since partly ionized metals still absorb X-rays by K-shell (1s electron) photoionization. Indeed, most of the IGM in IllustrisTNG is in small-scale filaments dominated, e.g., by Li-like O and lower charge states (Nelson et al. 2018b).

Zoom In Zoom Out Reset image size

In order to compute the optical depth in the simulation, we sum separately over the neutral atomic species of H, He, and all metals (Z) treated as a whole:

$$\begin{eqnarray}&&\tau ({E}_{0},z)=\displaystyle \sum _{s\in \left\{{\rm{H}},\mathrm{He},{\rm{Z}}\right\}}{\int }_{0}^{z}{n}_{s}(z^{\prime} ){\sigma }_{s}[E(z^{\prime} )]c\displaystyle \frac{{dt}}{{dz}^{\prime} }{dz}^{\prime} ,\end{eqnarray} \tag{ 3 }$$

where σ_s is the photoelectric cross section. For the metals, we use a cross section weighted by relative solar abundances, for absorption at $z^{\prime} ,{\sigma }_{s}(E)$ needs to be evaluated at the de-redshifted energy of $E(z^{\prime} )={E}_{0}(1+z^{\prime} )$. Since σ_s strongly decreases with energy, as can be seen in Figure 3, higher-z absorbers require higher (column) density to match the observed τ. c is the speed of light.

Zoom In Zoom Out Reset image size

In terms of the cosmological model, τ is written as

$$\begin{eqnarray}&&\begin{array}{l}\tau ({E}_{0},z)\\ =\displaystyle \frac{c}{{H}_{0}}\displaystyle \sum _{s\in \left\{{\rm{H}},\mathrm{He},{\rm{Z}}\right\}}{\int }_{0}^{z}\displaystyle \frac{{n}_{s}(z^{\prime} ){\sigma }_{s}[{E}_{0}(1+z^{\prime} )]}{(1+z^{\prime} )\sqrt{{{\rm{\Omega }}}_{M}{\left(1+z^{\prime} \right)}^{3}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}}}{dz}^{\prime} ,\end{array}\end{eqnarray} \tag{ 4 }$$

where H₀ = 67.7 km s⁻¹ Mpc⁻¹ is the Hubble constant. Ω_M = 0.31 and Ω_Λ = 0.69 are the mass and dark energy fractions of the present-day energy density of the universe. The integral over $z^{\prime}$ in Equation (4) is carried out as a discrete sum over the simulation snapshots, in ${dz}^{\prime}$ steps determined by the ${\rm{\Delta }}z^{\prime}$ between snapshots. Random lines of sight are chosen in each snapshot, to avoid the line of sight going through the same structure over and over again (Figure 1). Structures larger than the simulation box therefore cannot be accounted for. We repeat the calculation of Equation (4) 300 times, in order to study the statistics over many random lines of sight, similar to real GRB samples. A consistency check of the simulation confirms an average comoving baryon density of 2 × 10⁻⁷ cm⁻³ and hydrogen column densities (computed by Equation (3), but without σ) of N_H = 4 × 10²¹ up to z = 1 and 10²³ cm⁻³ up to z = 10, which are consistent with mean analytical values (e.g., Behar et al. 2011).

Figure 4 shows the results of our calculations for τ(z) according to Equation (4) and averaged over 300 random lines of sight. The τ(z) values in the random lines of sight represent the opacity of the IGM expected for random GRB locations. We tested two cases. In the first case, we assumed all atoms were neutral. The results are shown in the upper part of Figure 4. It can be seen that τ at high z approaches 0.9, which is higher than what the GRB afterglow observations indicate. Also, it can be seen that the He opacity dominates τ over H and the small contribution of metals. When we include the effects of ionization (Figure 2), the contribution of H and He becomes negligible for z < 3, while the metals dominate. The values obtained at z = 10 are τ_H ≈ 0.008, τ_He ≈ 0.05, and τ_Z ≈ 0.09, for a total of τ_total ≈ 0.15. This can be seen in the lower part of Figure 4. The figure shows that τ is not sensitive to the exact epoch in which He fully ionized, unless it ionized at z < 2.5, which could increase τ.

Zoom In Zoom Out Reset image size

The result above includes ionization, which is more physical, and shows that τ is determined almost entirely by the metallicity, and mostly the metallicity at low z. Interestingly, τ ∝ Z₀ agrees well with the crude approximation of a homogeneous IGM presented in Equation (1). The average metallicity, however, in the simulation is surprisingly low, peaking at Z₀ ≈ 0.06 at z = 0 and decreasing below 0.01 at z = 3, as shown in Figure 5. In terms of Equation (1), the metallicity index is k = 1.7. These metallicities are low compared to the measurements of metallicity of ∼0.3 (median) in circumgalactic media (Prochaska et al. 2017) and a metallicity approaching solar in the IGM along lines of sight to GRBs and quasars (see the compilation in Figure 2 of Savaglio 2009).

Zoom In Zoom Out Reset image size

The factor that goes into the calculation of τ in Equation (4) is not the metallicity, but the physical metal density ${n}_{M}(z^{\prime} )$. Figure 6 shows n_M as a function of z averaged over 300 lines of sight. n_M increases with star formation, peaking between z = 2 and 3, and it decreases at low z due to the expansion of the universe. This behavior—along with the decreasing σ_Z (E)—explains why τ increases in Figure 4 up to z ∼ 2 and not beyond. The dominance of metals in the bottom panel of Figure 4 deems Equation (1) to be a good approximation, demonstrating that τ_IGM(z → ∞ ) is determined predominantly by the metallicity at low z.

Zoom In Zoom Out Reset image size

The results of the IllustrisTNG simulation are compared in Figure 7 with 353 GRB Swift/XRT measurements of τ, including 123 upper limits, taken from Rahin & Behar (2019). Overall, the measured GRB opacities are higher than the average values from the simulation. The two come close at high z; both data sets have significant scatter. On average at high z, the simulation results tend to τ = 0.15 ± 0.07 (95%), while the observed data tend to τ ≈ 0.4, starting at z ∼ 2–3. Only a few GRBs are available for z > 5.

Zoom In Zoom Out Reset image size

The main reason for the low τ value in the simulation is its low average metallicity at low z (Figures 5 and 6). Complicating the comparison is the absorption in the host galaxy. We do not attempt to compute this contribution to τ, but note that for a fixed host column density (e.g., N_H = 10²¹ cm⁻²; Behar et al. 2011), it will strongly diminish with z, due to the decreasing cross section (Figure 3).

Indeed, the difference between the measured τ values and those from the simulation decreases in Figure 7 with z. However, the high-z measured values remain about three times higher than in the simulation. In order for the IGM to explain the observed τ ≈ 0.4 at high z, one needs a z = 0 metallicity of Z₀ = 0.2, rather than Z₀ = 0.06 from the simulation. Indeed, we compute τ with Z₀ = 0.2 and obtain z = 10, τ_Z ≈ 0.3, and τ_total ≈ 0.36, within the measured range. Note that the median τ values in Figure 7 include the upper limits. Hence, the true median of τ(z) of GRB afterglows is below the plotted value and closer to the simulation of the IGM. This reduces the required host contribution.

Another difference between the measured τ values and those from the simulation is the distribution around the mean. The 95% range of the lines of sight in the simulation is represented in Figure 7 by the blue shaded area. The distribution about the mean is rather normal and bounded by a factor of a few. On the other hand, the distribution of observed values (green points and red triangles) is much broader and flatter, exceeding an order of magnitude at most redshifts. At low z, this can be ascribed to the variety of host galaxies and the GRB environment. At z > 3, where the IGM becomes important, the distribution narrows, yet suffers from a deficiency of measurements.

We employed the IllustrisTNG simulations to estimate the total IGM X-ray optical depth τ toward high-z sources and compared it with measured GRB afterglows. The conclusions can be summarized as follows:

• 1.  
  Ionization effects in the IGM are important. Neglecting them results in an overestimation of τ by a factor of 6 (Figure 4).
• 2.  
  The early ionization of H and He between 2.5 < z < 5 (Figure 2) leaves the metals as the dominant (>60%) X-ray absorbers in the IGM (Figure 4, bottom).
• 3.  
  The 0.5 keV IGM opacity in the simulation is τ = 0.15 ± 0.07, which is sufficient to explain the observed values of ≈0.4 (Figure 7). The opacity difference is explained by a residual host contribution, by the low metallicity in the simulation (Figure 5), and by the numerous GRBs with τ upper limits (Figure 7).
• 4.  
  The broad distribution of GRB afterglow τ values at low z must be a result of the variety of host environments. At high z, the distribution narrows and better agrees with simulations (Figure 7).
• 5.  
  The simulations show that the IGM can account for a major fraction of the missing baryons. The exact contribution depends strongly on the IGM metallicity and ionization. These need to be better constrained by observations in order to benchmark the simulations.

This work was supported by a Center of Excellence of Israel Science Foundation (grant No. 1937/19). The authors thank Adi Nusser for useful discussions. The IllustrisTNG simulations were undertaken with computer time awarded by the Gauss Centre for Supercomputing (GCS) under GCS Large-Scale Projects GCS-ILLU and GCS-DWAR on the GCS share of the supercomputer Hazel Hen at the High Performance Computing Center Stuttgart (HLRS), as well as on the machines of the Max Planck Computing and Data Facility (MPCDF) in Garching, Germany. We thank an anonymous referee for useful comments that improved the manuscript.