© ESO 2022

1. Introduction

Just after inflation, the Universe’s dark and baryonic matter density functions resembled realisations of nearly constant, isotropic, and stationary Gaussian random fields (GRFs) (Linde et al. 2008). Due to gravity alone, these fields evolved into the highly inhomogeneous, network-like large-scale structure (LSS) present today (Springel et al. 2005). The late-time Universe consists of two components: voids and the Cosmic Web, which can be further partitioned into sheets, filaments, and (galaxy) clusters¹. The initial density conditions, through the morphology of the Cosmic Web, determine the spatial distribution of galaxies and some of their internal properties, such as the magnetic field and spin. Hahn et al. (2007), inspired by the seminal work of Zel’dovich (1970), provided a² rigorous definition (known as the ‘T-web’) of the four canonical structure types (voids, sheets, filaments, and clusters) based on the number of positive eigenvalues (0, 1, 2, or 3, respectively) of the tidal field tensor. Their N-body simulations reveal that the prevalence of the structure types (quantified by e.g., volume-filling factors) is evolving, with filaments and sheets disappearing in favour of clusters and voids at the present day³.

Of the four structure types, galaxy clusters are most easily studied: in the X-ray, optical and radio bands, for example through thermal bremsstrahlung, gravitational lensing and synchrotron radiation. As the most massive gravitationally bound structures in the Universe thus far, they weigh up to ∼10¹⁵ M_⊙, contain up to ∼10³ galaxies, and are pervaded by a dilute (∼10³ m⁻³), hot (∼10⁸ K) and magnetised (∼1 μG) hydrogen- and helium-dominated plasma: the intra-cluster medium (ICM) (Cavaliere & Rephaeli 2011).

At the outskirts of clusters, the ICM transitions into the warm–hot intergalactic medium (WHIM) – the dominant baryonic constituent of filaments. The WHIM is a plasma of nearly primordial chemical composition, but is less dense (∼1 − 10 m⁻³), cooler (∼10⁵ − 10⁷ K) and less magnetised (∼10⁻³ − 10⁻¹ μG) than the ICM. Compared to clusters, filaments are therefore harder to detect in all three wavelength bands. Despite this, filaments are cosmologically relevant, as simulations predict that the WHIM harbours up to ∼90% of the Universe’s baryons (Cen & Ostriker 1999; Eckert et al. 2015; de Graaff et al. 2019; Tanimura et al. 2019).

The formation of filaments (and, indirectly, the galaxy clusters they fuel) has occurred primarily through the influx of dark matter (DM) and cold gas from sheets and voids. Once these free-falling pockets of gas reach a filament’s surface, they generate supersonic accretion shocks, with upstream Mach numbers ℳ_u ∼ 10⁰ − 10² (Ryu et al. 2003). Following Ensslin et al. (1998), Miniati et al. (2001), Hoeft & Brüggen (2007) have proposed that merger shocks in clusters boost electrons in the high-energy tail of the ICM’s Maxwell–Boltzmann distribution to ultrarelativistic velocities via diffusive shock acceleration (DSA). By extension, the Hoeft & Brüggen (2007) model could also describe (again via DSA) how accretion shocks in filaments boost electrons in the high-energy tail of the WHIM’s Maxwell–Boltzmann distribution to ultrarelativistic velocities. The DSA mechanism details how charges diffuse back and forth across the shock front, trapped in a magnetic mirror, and gain speed accordingly (Krymskii 1977; Axford et al. 1977; Bell 1978a,b; Blandford & Ostriker 1978; Drury 1983; Blandford & Eichler 1987; Jones & Ellison 1991; Baring 1997; Malkov & Drury 2001). Once released, these high-energy electrons subsequently spiral along the magnetic field lines of the intergalactic medium (IGM), glowing in synchrotron light. It is thought that this accretion-shock-based radiation mechanism provides the dominant contribution to the filaments’ synchrotron Cosmic Web (SCW) signal.

Although the Mach numbers of accretion shocks in filaments are higher than those of merger shocks in clusters (where they are ℳ_u ∼ 1 − 5) (Ryu et al. 2003), shocks in filaments remain fainter due to the aforementioned adverse density, temperature and magnetic field strength conditions⁴. This is why observing filaments through synchrotron radiation constitutes an ambitious, futuristic frontier.

Modern radio telescopes, such as the upgraded Giant Metrewave Radio Telescope (uGMRT), the Expanded Very Large Array (EVLA) and the Low-frequency Array High-band Antennae (LOFAR HBA), have enabled detailed studies of particle acceleration in the cluster IGM (e.g., Di Gennaro et al. 2018; Kale 2020; Locatelli et al. 2020; Mandal et al. 2020) by the detection of radio halos, phoenices and relics. In contrast, no single shock in filaments has hitherto been observed. Doing so would open up density, temperature, Mach number and magnetic field strength regimes different by orders of magnitude via which astrophysical shock models could be held to the test. For example, the top panel of Fig. 1 shows that DSA predicts an (angularly unresolved) synchrotron spectral index α = −1 for virtually all filament shocks (with a slight dependency on the adiabatic index γ). Significant deviations from α = −1 would falsify standard DSA⁵. A better understanding of astrophysical shocks has ramifications beyond the study of large-scale structure, as possibly similar shocks are found in, amongst other places, accreting X-ray binaries, stellar and pulsar winds, and supernova remnants.

Fig. 1. Dependence of three central (merger or accretion) shock quantities on upstream Mach number ℳu and adiabatic index γ, as derived from the Rankine–Hugoniot jump conditions in ideal gases. Radiation from ℳu < 10 shocks is predominantly generated by electrons that stem from the ICM, whilst radiation from ℳu > 10 shocks is predominantly generated by electrons that stem from the WHIM. Top: the spectral index α of angularly unresolved shocks, assuming standard DSA. Middle: the downstream-over-upstream plasma density ratio. Bottom: the downstream-over-upstream plasma temperature ratio.

Routine detection of the filament IGM in synchrotron would provide a novel way to address the missing baryon problem – the possible discrepancy between today’s mean baryon density as inferred from galaxy surveys versus that predicted by CMB-constrained ΛCDM models.

Direct imaging of the filament IGM in the low-frequency radio window adds spectral diversity to a growing list of methodologies that trace the WHIM, complementing X-ray observations of ionised oxygen (O VII) absorption along the line of sight to quasars (Nicastro et al. 2018), microwave measurements of the thermal Sunyaev–Zel’dovich effect due to hot gas between adjacent galaxies (de Graaff et al. 2019), millimetre searches for the hyperfine spin-flip transition of single-electron nitrogen ions (N VII) (Bregman & Irwin 2007) and dispersion measurements of localised fast radio bursts (FRBs) (Macquart et al. 2020). Just as direct imaging of the WHIM in the X-ray band (Eckert et al. 2015), detecting baryons through synchrotron emission does not necessitate a special (line-of-sight) geometry. Telescopes such as the LOFAR could thus corroborate the current baryon census, in which still ∼34% is not identified conclusively (de Graaff et al. 2019). Thirdly, the low galaxy number density in filaments means that its IGM, if largely untouched by galactic feedback, can retain information on its initial conditions for billions of years. In particular, MHD simulations (Vazza et al. 2015, 2017) demonstrate that different assumptions for the dominant physical process that drove the growth of magnetism in filaments, lead to different strengths and morphologies of the IGM’s magnetic fields today. The authors evaluate three different scenarios for cosmic magnetogenesis: primordial, dynamo, and galactic models, and calculate their evolution from a cosmological redshift of 38 to the present day. The models are calibrated by magnetic field strength measurements of the ICM. In primordial models, seed field fluctuations grow as LSS formation compresses and rarefies the magnetic field lines. In dynamo models, the seed field is much weaker (10⁻²² T vs. 10⁻¹³ T), but grows in strength as energy in solenoidal, turbulent gas motion is converted into magnetic energy through small-scale dynamos (Ryu 2008). In galactic (or ‘astrophysical’) models, finally, no seed field is assumed, with star formation and outflows from the jets of supermassive black holes (SMBHs) at galactic centres being the dominant contributors to the emergence of the magnetised cosmic web. Figure 2 illustrates that different magnetisation scenarios give rise to morphologically different low-frequency specific intensity functions. Thus, SCW observations could rule out at least some models of magnetogenesis – a puzzle widely considered to be amongst the most significant open problems in cosmology.

Fig. 2. Two simulated specific intensity functions at νobs = 150 MHz, assuming synchrotron radiation from merger and accretion shocks (as in Hoeft & Brüggen 2007) in LSS at z = 0.045. Top: primordial scenario for magnetogenesis, starting from BIGM = 1 nG in the Early Universe. (In contrast, the simulation underlying this article’s predictions starts at BIGM = 0.1 nG.) Bottom: galactic scenario for magnetogenesis, in which seeding occurs through AGN outflows and supernova winds. See Gheller & Vazza (2020) for details on these scenarios.

Hitherto, several groups have constrained the magnetic field strength of the filament IGM by means of Stokes I low-frequency radio observations. Using Murchison Widefield Array (MWA) Epoch of Reionisation Field 0 data and the equipartition energy condition, Vernstrom et al. (2017) derived a parameter-dependent upper limit of B_IGM < 0.03 − 1.98 μG, while Brown et al. (2017) used Parkes 64m Telescope S-PASS data to find a density-weighted upper limit of B_IGM < 0.13 μG at the present epoch. Finally, using MWA GLEAM and ROSAT RASS data, Vernstrom et al. (2021) suggest an average magnetic field strength of 30 nG < B_IGM < 60 nG, using both equipartition and inverse Compton arguments. In an upcoming publication (Oei et al., in prep.), we describe and present a LOFAR search, that includes new methodology.

Alternatively, O’Sullivan et al. (2019, 2020), Stuardi et al. (2020) constrain the properties of the IGM’s magnetic fields via rotation measure (RM) synthesis applied to LOFAR Two-metre Sky Survey (LoTSS) observations of the lobes of (giant) radio galaxies, finding B_IGM < 4 nG in filaments. Vernstrom et al. (2019) have performed a similar analysis with NVSS data, finding B_IGM < 40 nG in filaments.

The goal of this article is, first and foremost, to explain a new method for SCW prediction, that yields a probability distribution over the set Map(𝕊²,ℝ_≥0) of specific intensity functions on the 2-sphere 𝕊². The secondary goal is to demonstrate the method’s potential, by showing results for the modern, nearby Universe (z < 0.2) over half of the Northern Hemisphere – 25% of the full sky.

In Sect. 2, we first provide a general outline of our SCW signal prediction method. We then detail the methodology step-by-step, providing a background to the input data as we proceed. In Sect. 3, we develop a simple geometric model of the cosmic web that yields SCW redshift predictions. In Sect. 4, we analyse predictions of the main method, and compare these with predictions of the geometric model. Finally, in Sect. 5, we discuss caveats of the current work, and give recommendations for future extensions, before we present conclusions in Sect. 6.

We adopt a⁶ concordance inflationary ΛCDM cosmological model 𝔐 = (Ω_Λ, 0 = 0.728, Ω_DM, 0 = 0.227, Ω_BM, 0 = 0.045, h = 0.702, σ₈ = 0.807, n_s = 0.961), so that Ω_M, 0 := Ω_DM, 0 + Ω_BM, 0 = 0.272 and the Hubble constant H₀ := h ⋅ 100 km s⁻¹ Mpc⁻¹ (Jasche et al. 2015).

2. Methods

2.1. Overview

For structural clarity, we begin with an overview of our synchrotron Cosmic Web prediction approach.

1. Our starting point is a probability distribution over the total (i.e. baryonic and dark) matter density fields of the modern, nearby Universe.

2. Then, using a cosmological MHD simulation, we established the connection between the total matter density ρ and the merger and accretion shock synchrotron monochromatic emission coefficient (MEC) j_ν, in the form of the conditional probability distribution P(j_ν(r) | ρ(r)).

3. Next, to create a MEC field, one could independently realise the random variables (RVs) j_ν(r) | ρ(r) for all r in the volume of the inferred density fields. However, such an approach would disregard the fact that in reality, adjacent locations have similar physical conditions and thus similar MECs (i.e. the MEC field exhibits spatial correlation). To account for spatial correlation, we generated realisations of an isotropic, stationary, three-dimensional Gaussian random field (GRF) with zero mean and unit variance. By applying the cumulative density function (CDF) of the standard normal RV in point-wise fashion, the GRF transforms into a ‘percentile random field’ – with values strictly between 0 and 1. This new field inherits the spatial correlation present in the GRF. We then let the percentile random field determine the MEC field, by plugging the percentile scores into the inverse CDF of j_ν(r) | ρ = ρ(r) for all r in the volume.

4. After having obtained a merger and accretion shock synchrotron MEC field with spatial correlation, we ray traced through this three-dimensional structure, taking into account effects due to the Universe’s expansion. This, finally, yields a merger and accretion shock synchrotron Cosmic Web (MASSCW) specific intensity function on the sky.

5. By repeating the above steps over and over, picking different density field and GRF realisations every time, we generated a probability distribution on Map(𝕊²,ℝ_≥0) for the MASSCW signal.

6. With a minor modification, this procedure naturally leads to another probability distribution on Map(𝕊²,ℝ_≥0), which captures the direction-dependent mean cosmological redshift of the MASSCW.

2.2. Modern-day total matter density field posterior

2.2.1. A brief history of Bayesian large-scale structure reconstruction

A decade of research in large-scale structure reconstruction from galaxy surveys has culminated in a suite of highly principled, physics-based Bayesian inference techniques that unveil the content, dynamics and history of the local Cosmic Web. These techniques represent the state-of-the-art of LSS reconstruction, and their applications are manifold (see e.g., Jasche & Lavaux 2019). The application relevant to our purposes is the reconstruction of the modern-day (i.e. z = 0) total (i.e. baryonic and dark) matter density field of the local Universe, in the form of a probability distribution over all possible density fields.

Several variations of the same general framework exist. Inspired by the map-level CMB inference methods of Wandelt et al. (2004), some theoretical groundwork common to all these techniques has been developed by Kitaura & Enßlin (2008), Jasche et al. (2010) and Jasche & Kitaura (2010), who also demonstrated the practical feasibility of Bayesian LSS reconstruction with the ARGO, ARES and HADES codes, respectively. In Jasche & Wandelt (2013), the authors described BORG (Bayesian Origin Reconstruction from Galaxies), the first algorithm to include physics of structure formation, and provided a proof-of-concept on simulated Sloan Digital Sky Survey Data Release 7 (SDSS DR7) data. In Jasche et al. (2015), the actual SDSS DR7 main galaxy sample was analysed with a refined version of BORG, yielding the publically-available BORG SDSS data release. These data form the basis for the SCW prediction method presented in this paper.

The BORG SDSS has hitherto been used for Cosmic Web classification (Leclercq et al. 2015), for the study of galaxy properties as a function of environment (Leclercq et al. 2016) and to unveil the dynamics of DM streams (Leclercq et al. 2017), all in a probabilistic and time-dependent fashion. Further refinements to BORG have been presented alongside applications to new datasets, namely 2M++ (Lavaux & Jasche 2016; Jasche & Lavaux 2019) and SDSS3-BOSS (Lavaux et al. 2019). Algorithms related to BORG are ELUCID (Wang et al. 2013, 2014), COSMIC BIRTH (Kitaura et al. 2021) and BARCODE (Bos et al. 2019).

2.2.2. Main principles

Only a brief description of the main ideas underlying these approaches falls within the scope of this paper. The methods all use the crucial insight that the statistical properties of the total density field in the Early Universe are well understood (i.e. Gaussian, with a theoretically predicted and observationally verified covariance function), and that the modern-day total density field relates to the initial field deterministically by means of (approximately Newtonian) gravity⁷. Thus, realisations from today’s highly non-Gaussian total density field can be generated by forward modelling the effect of gravity on a collisionless⁸ fluid, which is initialised as a Gaussian random field at the time of the CMB⁹. After a CMB-epoch total density GRF has been generated, and evolved into a modern total density field, a Poisson point process – with an intensity function that attempts to capture the relation between galaxy locations and the surrounding matter distribution – is used to calculate the likelihood of finding galaxies at their measured locations assuming the true underlying total density field is the one currently considered. As suggested before, spectroscopic galaxy surveys, with hundreds of thousands of galaxies pinpointed in three-dimensional comoving space, form the typical input data. Most notably, the likelihood for the modern total density field also fixes the likelihood for the initial total density field from which the modern one was evolved. Because the prior on the initial total density fields is Gaussian, the initial total density field can not only be assigned a likelihood, but also a posterior probability. Hamiltonian Monte Carlo (HMC) Markov chains provide the means to explore the high-dimensional posterior distribution of initial total density fields consistent with the galaxy data. A posterior distribution over modern-day total density fields, which is used in this article, is generated as a by-product of this process.

2.2.3. BORG SDSS

The BORG SDSS modern-day total density field posterior used in this article is based on optical galaxy data from the Sloan Digital Sky Survey (SDSS) (York et al. 2000; Strauss et al. 2002) Legacy Survey (Abazajian et al. 2009). This catalogue provides right ascensions ϕ, declinations θ, spectroscopic redshifts and Petrosian r-band apparent magnitudes m_r of about a million galaxies in the Northern Galactic Cap (NGC), which covers a region of 7646 deg² of Northern Sky, and three stripes of the Southern Galactic Cap (SGC), which together cover 386 deg². In these regions, the SDSS Legacy Survey is spectroscopically complete for galaxies with m_r < 17.77.

The BORG SDSS, which is based on 372 198 NGC galaxies (see Jasche et al. 2015 for selection details), infers structures in the density field of ∼3 Mpc h⁻¹ and larger. Each realisation of the posterior is a cube containing 256³ voxels, that represents a region ℛ ⊂ ℝ³ with comoving volume (750 Mpc h⁻¹)³. Naturally, ℛ is chosen to correspond to the half of the Northern Sky that contains the NGC (up to a cosmological redshift z ∼ 0.2).

2.3. Relating total matter density to synchrotron monochromatic emission coefficient

Next, to convert total matter density fields ρ to proper synchrotron monochromatic emission coefficient (MEC) fields j_ν, we established the relation between ρ and j_ν (including its variability). In radiative transfer theory, the MEC quantifies the amount of radiative energy released per unit of time, volume, frequency and solid angle (Rybicki & Lightman 1986).

2.3.1. The accretion shock Ansatz

Cosmological simulations demonstrate that accretion shocks during large-scale structure formation are ubiquitous, and dominate the thermalisation of kinetic energy of baryons falling onto filaments. Given the suite of known particle acceleration processes, and the fact that such shocks are almost certain to exist, shock acceleration is expected to be the biggest contributor to the synchrotron emission from filaments. This paper follows this hypothesis, by only considering the merger and accretion shock contribution to the SCW.

Even under the Ansatz that merger and accretion shocks drive the SCW signal, considerable uncertainty surrounding the correct physical description remains. For example, it is an open question whether the electrons that eventually radiate in synchrotron light originate from the high-energy tail of the thermal Maxwell–Boltzmann distribution, or are rather accelerated to cosmic ray (CR) energies by SMBHs in the centres of galaxies, and then flung out via jets into the IGM. Apart from uncertainty in the source of energetic electrons, the complex nature of plasma physics also makes it hard to establish how already-energetic electrons attain ultra-relativistic energies in magnetised shock fronts, with the theory of diffusive shock acceleration (DSA) being just one of multiple scenarios. In DSA, CRs are accelerated by repeated crossings of the shock front, which acts as a magnetic mirror (Malkov & Drury 2001; Xu et al. 2020). Thus, we stress that the functional form of the SCW MEC, which we require to establish a connection with the total matter density field, depends on ill-constrained assumptions surrounding the exact radiation mechanism initiated by the accretion shocks.

2.3.2. The Hoeft & Brüggen (2007) model

Hoeft & Brüggen (2007) have derived an analytic expression for the synchrotron power density of cosmological shock waves, assuming that the radiating electrons exclusively originate from the high-energy tail of the thermal Maxwell–Boltzmann distribution, and that the electron energy spectrum at the shock front is well-described by DSA. As the DSA-based formulae of Hoeft & Brüggen (2007) have been partially successful in explaining observations of shocks in the ICM (e.g., van Weeren et al. 2019; Locatelli et al. 2020), we postulate that the same formulae can describe synchrotron emission due to accretion shocks onto filaments¹⁰.

We assume that the power density P_ν of a single shock obeys Eq. (32) of Hoeft & Brüggen (2007). Let the shock surface area be A, and let the effective width of the downstream region be ⟨y⟩ (intuitively, this is the thickness of the shock in the direction perpendicular to the surface). The effective shock volume is V = A⟨y⟩. If shocks radiate isotropically, the MEC j_ν is direction-independent. However, like ρ, j_ν is a volume-averaged quantity, so that it depends on the scale on which the averaging occurs. Index all shocks in the universe, so that P_ν, i(z,ν) is the power density of shock i at cosmological redshift z and emission frequency ν, and A_i and ⟨y⟩_i are its surface area and effective width¹¹. Moreover, let ℛ_i(z) ⊂ ℝ³ be the region of space occupied by this shock at redshift z; the Lebesgue measure of ℛ_i is V_i = ∫_ℛidr = A_i⟨y⟩_i. Under the approximation that the shock emission is homogeneous within ℛ_i, the average total MEC within an arbitrary region with Lebesgue measure is

(1)

The MEC of a single shock (thus dropping indices) located at r at redshift z, and averaged over its own effective volume, is

(2)

Here, n_e, d is the downstream electron number density, T_d is the downstream plasma temperature, α is the integrated spectral index of the associated synchrotron emission, B_d is the downstream plasma magnetic field strength, B_CMB is the CMB magnetic field strength, and Ψ is a dimensionless quantity with a strong dependence on the upstream Mach number ℳ_u and a weak dependence on T_d (as can be seen in Fig. 4 of Hoeft & Brüggen 2007). Ψ approaches unity for high Mach numbers, and so for the WHIM, where the upstream Mach numbers are expected to be high, Ψ(r,z) ≈ 1. Like Ψ, α does not depend on the spacetime coordinate (r,z) directly, but rather via ℳ_u and the adiabatic index γ. Concretely, for ℳ_u > 1,

(3)

Note that α → −∞ as ℳ_u → 1+, and that α → 1/4(1−3γ) as ℳ_u → ∞. The top panel of Fig. 1 illustrates this behaviour. The electron acceleration efficiency ξ_e quantifies the fraction of the shock’s thermal energy that is used to accelerate suprathermal electrons. Unfortunately, ξ_e is unknown for WHIM shocks, but a comparison to supernova remnant (SNR) shocks suggests ξ_e = 0.05 (Keshet et al. 2004).

2.3.3. Enzo simulations

To capture, in a statistical sense, the relationship between the total matter density and the Hoeft & Brüggen (2007) monochromatic emission coefficient, we turned to MHD simulations. In particular, we used snapshots (i.e. 3D spatial fields at constant time) of the largest uniform-grid cosmological MHD simulation to date (Vazza et al. 2019), which is based on the Enzo code (Bryan et al. 2014). These cubic snapshots cover a comoving volume of (100 Mpc)³ with 2400³ voxels, yielding a (comoving) resolution of per voxel edge. The simulations recreate the evolution of the baryonic and dark matter density functions ρ_BM and ρ_DM – as well as (thermo)dynamic quantities¹² such as the gas temperature T, magnetic field strength B and gas velocity v – under Newtonian gravity. (However, the effects of the expansion of the Universe as predicted by general relativity are still incorporated.) No galactic physics is included. Shocks in the snapshots can be identified by searching for temperature and velocity jumps, as described in Vazza et al. (2009)¹³. For shocked voxels, we simply took ⟨y⟩ to be the voxel edge length, and calculated α and Ψ after establishing ℳ_u. With the exception of ξ_e, which we assumed to be a constant throughout, the simulation thus enables us to compute all factors on the RHS of Eq. (2).

Taking their product for the z = 0.025 snapshot, we obtained j_ν, and compared it to ρ := ρ_BM + ρ_DM to study their relationship¹⁴.

2.3.4. Conditional probability distribution

In order to convert the BORG SDSS total matter density fields into synchrotron MEC fields, we needed to predict j_ν from ρ – ideally including uncertainty. Because the BORG SDSS has a comoving resolution of 4.17 Mpc per voxel edge compared to Enzo’s 4.17 × 10⁻² Mpc per voxel edge, we blurred both j_ν and ρ with a Gaussian kernel whose standard deviation is half of the comoving resolution ratio. We then treated each (j_ν, ρ) pair of the blurred fields as a draw from the joint distribution P(j_ν,ρ). By binning the data, we performed a simplistic form of kernel density estimation (KDE). Next, we calculated the conditional P(j_ν | ρ) from the joint by dividing it by the marginal .

The result is the set of probability distributions shown in Fig. 3. We overplot the single shock scaling relation expected for the WHIM regime: j_ν ∝ ρ^23/6. See Appendix B for a derivation.

Fig. 3. Conditional probability distributions P(jν | ρ) of synchrotron monochromatic emission coefficient jν (at rest-frame (emission) frequency ν = 150 MHz) given total matter density ρ for the Hoeft & Brüggen (2007) formalism, as derived from z = 0.025 snapshots of cosmological MHD simulations by Vazza et al. (2019). The dashed lines indicate a jν ∝ ρ23/6 (single shock) scaling relation for various proportionality constants. ρ is shown relative to the current-day mean total matter density , while the weakly constrained electron acceleration efficiency ξe is divided out from jν. Each conditional is numerically approximated by a probability mass function with 1000 bins, whose edges vary by a constant factor and span 40 orders of magnitude in jν.

2.4. Generating monochromatic emission coefficient fields

In principle, we could – for each BORG SDSS sample ρ(r) – convert to j_ν(r) by drawing from P(j_ν | ρ=ρ(r)) independently for every r ∈ ℛ. However, this is clearly suboptimal, as accurate SCW MEC fields exhibit spatial correlation up to megaparsec scales, both under the merger and accretion shock Ansatz as well as in a turbulence scenario (Govoni et al. 2019; Brunetti & Vazza 2020).

To generate spatially correlated draws from our conditional probability distribution, we used realisations of a Gaussian random field 𝒵 over ℛ. In general, the statistical properties of a GRF are determined by its mean, and its covariance function or kernel. In our case, we set the mean to 0, and chose an isotropic and stationary kernel, as suggested by the Cosmological Principle. From the variety of remaining choices commonly used in machine learning (ML), we picked¹⁵ the squared-exponential (SE) kernel K_SE : ℛ × ℛ → ℝ_≥0 (Rasmussen & Williams 2006):

(4)

The variance allows one to endow the GRF with dimensionality and scale it at will. As our GRF only served to draw spatially correlated MEC samples, we set σ_SE = 1.

The lengthscale l_SE is the characteristic scale of spatial correlation, and forms an important model choice. Some of the longest coherent shocks observed thus far are the Toothbrush Relic in galaxy cluster 1RXS J0603.3+4214 (van Weeren et al. 2012) and the Sausage Relic in galaxy cluster CIZA J2242.8+5301 (Di Gennaro et al. 2018), both with a largest linear size (LLS) of ∼2 Mpc. Although the majority of observed shocks seem smaller (suggesting l_SE < 2 Mpc), the presence of noise in radio imagery obfuscates the full extent of shock fronts, biasing measurements towards small LLSs. As a compromise, we chose l_SE = 2 Mpc.

As is customary in cosmological simulations – where GRFs with isotropic and stationary kernels are used to initialise the Early Universe matter density fields – we generated our zero mean, unit covariance GRF realisations on the BORG SDSS voxel grid using Fourier analysis. This necessitated calculating the power spectrum of the kernel, which in the case of the 3D SE kernel is

(5)

Finally, Fourier-generated GRFs on a finite (numerical) grid require an appropriate normalisation; we used Parseval’s theorem to find the correct factor.

This procedure generated a ‘standard normal’ GRF 𝒵, in the sense that each location’s RV is standard normal: 𝒵(r) ∼ 𝒩(0,1) for all r ∈ ℛ. Next, we applied the CDF Φ of the standard normal to the GRF in point-wise fashion (i.e. voxel-wise, in practice), thereby creating a ‘percentile random field’ 𝒫:

(6)

This field has values strictly between 0 and 1, and inherits the spatial correlations present in 𝒵.

Let Ξ_ρ(x) : ℝ_≥0 → [0, 1] be the CDF of the RV j_ν | ρ, and let be its inverse. Then the final MEC field is given by

(7)

Because ρ and 𝒫 are random fields on ℛ, so is j_ν. A graphical summary of the procedure described in this subsection, and defined in Eq. (7), is shown in Fig. 4.

Fig. 4. Overview of the procedure for converting a probability distribution over total (i.e. baryonic and dark) matter density fields into a probability distribution over monochromatic emission coefficient (MEC) fields. We show a fixed slice through ℛ, of 750 Mpc h−1 × 750 Mpc h−1 for the density and MEC fields, and a subregion of 150 Mpc h−1 × 150 Mpc h−1 for the Gaussian and percentile random fields (to more clearly illustrate the bijective mapping between the two). Left: BORG SDSS density ρ sample, in multiples of the current-day mean total matter density . Middle left: Gaussian random field 𝒵 sample using a squared exponential (SE) kernel with lengthscale lSE = 2 Mpc. Middle right: the corresponding percentile random field 𝒫 sample. Right: the corresponding MEC jν sample at an emission frequency of ν = 150 MHz, with the weakly-constrained electron acceleration efficiency ξe divided out.

2.5. Generating specific intensity functions

To find the specific intensity of the SCW on the sky, I_ν, we simulated the passage of light rays through our 3D MEC fields via ray tracing. Conveniently, projecting a ray’s 4D null geodesic in a pure Friedmann–Lemaître–Robertson–Walker (FLRW) metric onto 3D comoving space results in a straight line. We provide a brief derivation in Appendix C. By assuming an exact FLRW metric, we neglected all spacetime deformations due to local (large-scale structure) energy density fluctuations. Our results thus do not feature gravitational lens effects around massive clusters.

We estimated I_ν in a sky patch by sampling many directions within it, and simulating the passage of light through our reconstructed SCW in comoving space for each such direction (or ‘ray’) .

2.5.1. From proper monochromatic emission coefficient to observer’s specific intensity

Next, we establish the relation between the observed specific intensity I_ν of a ray that has travelled through the SCW to Earth from direction , and the MEC j_ν along its path.

Be j_ν : ℝ³ × ℝ_≥0 × ℝ_{> 0} → ℝ_≥0 the function that assigns to each 3-tuple (r,z,ν) containing a comoving locus r, a cosmological redshift z (which represents time) and an EM wave frequency ν, the proper (rather than comoving) MEC j_ν(r,z,ν) (with SI units W m⁻³ Hz⁻¹ sr⁻¹). Conceptually, we must differentiate between the specific intensity of large-scale structure at the time of emission (‘there and then’), and at today’s observing epoch (‘here and now’). Today’s quantity I_ν for direction and observing frequency ν_obs is

(8)

We provide a brief derivation, alongside explicit expressions for r(z) and E(z), in Appendix D.

For each specific intensity function that we wished to generate, we then drew M = 10⁶ ray directions uniformly from the π sr lune that covers the SDSS-constrained half of the Northern Sky. For each ray, we calculated the corresponding specific intensity by combining Eqs. (2) and (8).

Note that this requires evaluating the MEC field at a range of emission frequencies. Numerically, we realised this by building both the Enzo-derived conditional probability distribution P(j_ν | ρ) and the MEC field sample (using the same density field sample and GRF sample) at two emission frequencies ν: ν = ν_obs and ν = ν_obs(1+z_max). Next, for each voxel, we calculated a spectral index from the two MEC fields by assuming a power-law spectrum, and used it to find the MEC at the emission frequency ν = ν_obs(1+z) needed given the voxel’s cosmological redshift z.

We repeated this process many times by selecting 1000 density samples from the BORG SDSS posterior (we discarded the first 2110 samples due to burn-in, and thinned the Monte Carlo Markov chain by a factor 10) and generating one independent and identically distributed (IID) GRF sample for each. (There is no compelling reason to thin the chain and we could have used multiple GRF samples per density sample; the current approach is merely to limit data storage, compute time and energy usage.) Thus, we numerically realised 1000 ⋅ 1 = 1000 specific intensity functions. Together, they form a probability distribution over MASSCW specific intensity functions over a quarter of the sky.

2.6. Generating redshift functions

Our methodology can also be used to characterise the redshift properties of the MASSCW.

2.6.1. Specific-intensity–weighted mean redshift

A ray’s total specific intensity is found by summing up the specific intensity contributions from all voxels along its path, each of which has a known redshift. Therefore, a ray’s mean redshift is found by weighing each voxel’s redshift by the corresponding specific intensity contribution, before dividing by the sum of such contributions. Concretely, we define the specific-intensity–weighted mean redshift at observing frequency ν_obs of a ray passing through MEC field j_ν to be

(9)

By iterating over total matter density and GRF samples, we generated a probability distribution over the function , in exact analogy to the generation of the distribution over specific intensity functions .

One might wonder what the specific-intensity–weighted mean redshift of the MASSCW signal is for a randomly picked direction on the sky. The set , with rays drawn from a uniform distribution over the sky, can be viewed as a random sample from a specific-intensity–weighted mean redshift random variable . An RV is fully characterised by its CDF. Each MASSCW prior realisation generates another empirical CDF (ECDF) for .¹⁶

2.6.2. Flux-density–weighted mean redshift

One might also wonder what the mean redshift of the MASSCW signal is – not for a single direction, but overall. Due to strong attenuation with redshift, the MASSCW specific intensity is usually ‘high’ for directions with a ‘low’ specific-intensity–weighted mean redshift, and vice versa. This means that, although the signal could originate from far away for most of the sky (given sufficiently small cluster and filament volume-filling fractions), the sky-averaged redshift can still be low. We define the flux-density–weighted mean redshift of the MASSCW signal at observing frequency ν_obs to be

(10)

where is shorthand for .

3. Redshift predictions from geometric cosmic web model

3.1. Overcoming the redshift limitation

Although we can calculate the specific intensity function , the specific-intensity–weighted mean redshift function , the single-direction redshift RV ECDF and the flux-density–weighted mean redshift for each realisation of the MASSCW prior, all four quantities suffer from the fact that the BORG SDSS reconstructions stop at a redshift z_max = 0.2. It is of prime interest to know to what extent such redshift limitations, which we anticipate will become less stringent in the future, affect our inferences. For example, if most of the real MASSCW signal originates from z > 0.2, the predictive power of our specific intensity function distribution would be limited.

Here we introduce a simple analytical model that allows us to calculate both the distribution of , and , for an arbitrary value of z_max. The model reproduces the results calculated from the MASSCW prior, if it is equally limited to z_max = 0.2. This provides tentative evidence that the model captures the essential elements of the MASSCW signal, so that it might be used for extrapolation. We thus use the model to predict the distribution of and , for the case that describes the actual Universe: z_max → ∞.

3.2. Model formulation

To retrieve MASSCW redshift properties, we propose a simple geometric model. The key idea is that the dominant MASSCW redshift for an arbitrary sightline can be calculated via a weighted sum of the redshifts of the clusters and filaments it passes through, where the weights are set by the geometry and MECs of the LSS pierced. The redshifts of the structures a sightline pierces through depend on the particular LSS realisation the observer is immersed in. However, the cosmological principle dictates that LSS is statistically similar everywhere, and thus typical geometric cosmic web parameters must exist (and can be retrieved from numerical simulations; see e.g., Gheller et al. 2015).

Here we assume that sightlines pass through an arrangement of identical cubic unit cells – ignoring morphological variation among clusters, filaments, sheets and voids – with comoving edge length l. The edges are surrounded by square cuboids, which represent filaments, of typical comoving width w_f and length l_f ≔ l − w_f (and thus volume )¹⁷. Galaxy clusters reside where filaments meet, and – in this simplistic model – are cubes with comoving edge length w_f (and thus volume ). By far the largest component of a unit cell is its central void, a cubical region of comoving edge length l_f (and thus volume ). Finally, unit cells contain sheets, which fill the regions bounded by filaments and voids. They are of typical comoving length l_f and thickness w_f (and thus volume ). Figure 5 depicts a face of a unit cell. If we were to assume a perfect crystal structure for the cosmic web, some sightlines would never encounter clusters or filaments whilst others would consistently do so at regular (comoving) intervals. To avoid this unphysical scenario, we assume that our sightline-of-interest enters every unit cell it encounters at a random position on the face of incidence. For additional simplicity, we assume that the sightline always hits such faces perpendicularly.

Fig. 5. Geometry of the simple model, in which we treat the cosmic web as a collection of randomly displaced cubic unit cells. We show one face (grey square) of a unit cell, with edges of comoving length l. Filaments of typical comoving width wf and length lf ≔ l − wf cover the light-blue–shaded region of the face, and extend to the faces of neighbouring unit cells (grey dashed lines). Where filaments meet, galaxy clusters reside (orange-shaded regions), each of which is also connected to a filament oriented away from the observer. The central part of the face represents a sheet (dark blue), behind which a large void looms (not depicted). The figure is not to scale: filament widths are exaggerated compared to realistic filament lengths.

If one would choose a point-of-incidence on the face shown in Fig. 5 in uniform fashion, the sightline would hit a sheet and then a void (dark-blue-shaded area) with probability

(11)

Note that p_s − v → 1 when , and p_s − v → 0 when , as required. Likewise, the probability that the sightline hits a filament only (light-blue-shaded area) is

(12)

Finally, the probability that the sightline hits a cluster first, and then the filament behind it (orange-shaded area) is

(13)

Figure 6 depicts these pierce probabilities as a function of the ratio between filament width and unit cell edge length . For a comparison between the LSS volume-filling factors (VFFs) of this simple model and those from cosmological simulations, see Appendix E.

Fig. 6. When a sightline enters a cubic unit cell of the cosmic web perpendicularly to a face (as in Fig. 5) and at a random position on it, one of three different LSS piercing events occurs. We show the probabilities for each event as a function of the ratio of the typical filament width and cubic unit cell edge length (see Eqs. (11)–(13)).

Pick a sightline, and successively label each unit cell boundary this sightline crosses with a natural number n ∈ ℕ_≥1. If the first boundary crossing occurs at a comoving distance d₁, then the nth crossing happens at comoving distance d_n = d₁ + (n−1)l. Let 𝔐 be the cosmological model of preference, and let z_𝔐(d) denote the function that converts comoving distance to cosmological redshift under this model. The unit cell boundary crossings occur at redshifts z_n ≔ z_𝔐(d_n) for all n ∈ ℕ_≥1.

Because we want to be able to calculate the dominant redshift of the MASSCW signal produced by filaments and clusters up to a given cosmological redshift z_max only, we introduce a parameter N. Intuitively, N ∈ ℕ_≥1 is the label of the last boundary crossing within the LSS considered. More formally, N ≔ max {n ∈ ℕ_≥1 | z_n < z_max}.

Let X_n(ν_obs) be the RV denoting the contribution to the MASSCW specific intensity at observing frequency ν_obs picked up by the sightline during the nth full unit cell crossing (that is whilst travelling through the nth newly entered unit cell; the unit cell the observer resides in does not contribute). To retain low complexity, we assign the complete specific intensity contribution from the nth crossing to the redshift of incidence z_n. It is useful to regard X_n as an instance of a more general random variable: X_n(ν_obs) ≔ X(ν_obs,z_n). X(ν_obs,z) is the RV with support (0,ΔI_ν, f(ν_obs,z),ΔI_{ν, c − f}(ν_obs,z)), and corresponding probabilities (p_s − v,p_f,p_c − f). Note that we assume that sheets and voids have a vanishing contribution to the MASSCW signal. Here,

(14)

(15)

The factor (1+z)⁻⁴ follows from Eq. (D.7): the specific intensity as an integral over comoving radial distance. As we assume LSS to occur at regularly-spaced radial comoving distances d_n, this integral can be approximated by a Riemann sum, where each term equals the integrand multiplied by Δr = l.

Although redshift-dependent, we assume the MEC to be constant (at a given emission frequency) within – but different between – clusters and filaments of the same unit cell. Respectively, j_ν, f(ν,z) and j_ν, c(ν,z) represent the typical proper filament and cluster SCW MECs at emission frequency ν and cosmological redshift z. We use the dimensionless parameter 𝒞 to denote the typical cluster-to-filament-SCW MEC ratio at some reference emission frequency ν = ν_ref and z = 0: j_ν, c(ν_ref,0) ≔ 𝒞j_ν, f(ν_ref,0). To retain minimal complexity, we propose to describe the spectral and temporal dependencies of both j_ν, f(z) and j_ν, c(z) as power laws in and (1+z):

(16)

(17)

with α_f, β_f, α_c and β_c constants. Under this choice of MEC function parametrisation, the expectation value of X(ν_obs,z) becomes

(18)

(19)

3.3. Specific-intensity–weighted mean redshift

Consider the random vector [X₁,X₂,...,X_N]^T. Whenever max{X₁, X₂, ..., X_N}> 0, we can define the MASSCW specific-intensity–weighted mean redshift RV

(20)

We stress that , being a deterministic function of the RVs {X₁, X₂, ..., X_N}, is itself an RV. It quantifies the variety of specific-intensity–weighted mean redshifts that could occur in different LSS realisations (with the same model parameters) for a fixed direction (sightline) of the sky. Conversely, it could be interpreted as representing the redshift variety that occurs in a fixed LSS realisation for different directions of the sky, provided that the sky contains many (almost) independent patches of LSS. This second interpretation invites for a comparison between the specific-intensity–weighted mean redshift distribution of our MASSCW prior skies, and the distribution of .

Because X_n(ν_obs) features in both numerator and denominator of the fraction in Eq. (20), redshift-independent multiplicative factors cancel out. Therefore, j_ν, f(ν_ref,0) need not be specified; also, becomes ν_obs-independent if α_f = α_c. In general, the distribution of is determined by the cosmological model 𝔐 and 9 additional parameters (w_f, l, d₁, z_max, α_f, β_f, α_c, β_c and 𝒞). Using d₁, we can force the filaments and clusters nearest to the observer to occur at a fixed comoving distance. Alternatively, one could randomise the observer’s position with respect to the unit cells by choosing d₁ ∼ Uniform(0,l), making it an RV. In such case, 8 parameters remain that determine the distribution of .

3.4. Flux-density–weighted mean redshift

The geometric model also allows to calculate the flux-density–weighted mean (that is sky-averaged) MASSCW redshift . We did so by considering M ∈ ℕ_≥1 sightlines (instead of one) and their associated specific-intensity–weighted mean redshift RVs in the sense of Eq. (20). In total, these RVs depend on a set of N ⋅ M discrete RVs X_nm, which can be partitioned into N subsets of M IID RVs: X_nm ∼ X_n. We find by summing the elements of , each weighted by the corresponding specific intensity, and dividing the result by the sum of these weights; all whilst we let M → ∞:

(21)

(22)

(23)

Here we use the fact that the limit of a ratio of two sequences equals the ratio of the sequences’ limits (provided that the denominator sequence does not converge to 0), and that the sample mean as M → ∞, for all n ∈ {1, 2, ..., N}. Like before, as 𝔼(X_n(ν_obs)) features in both numerator and denominator of the fraction in Eq. (23), redshift-independent multiplicative factors cancel out. Thus, is ν_obs-independent when α_f = α_c.

Equation (19) makes clear that clusters dominate over filaments at cosmological redshift z when

(24)

At ν_obs = ν_ref and z = 0, this inequality reduces to

(25)

which is amply satisfied (see Sect. 3.5 for parameter estimates). Near ν_obs = ν_ref and for z of order unity, the inequality continues to hold as long as α_f − α_c + β_f − β_c is a number of order unity at most. In addition, most redshifts high enough that the inequality is violated, do not contribute meaningfully to . Therefore, around ν_obs = ν_ref is typically dominated by clusters, and . As long as Eq. (24) holds, becomes insensitive to variations in 𝒞, w_f, α_f and β_f; setting z_max = ∞, just 3 significant parameters remain: l, d₁ and α_c + β_c.

Finally, it is of significant interest to calculate for the filaments’ SCW signal only, discarding the dominant SCW signal from merger and accretion shocks in galaxy clusters. To find a filament-only with the geometric model, we just set 𝒞 = 0. This simplifies the expression for 𝔼(X(ν_obs,z)), and reduces the number of relevant parameters considerably. When 𝒞 = 0, becomes not only independent of α_c and β_c, but also of w_f and ν_obs:

(26)

As before, if we set z_max = ∞, is a function of just three parameters: l, d₁ and α_f + β_f, in this case.

3.5. Parameter estimates

This subsection provides concrete parameter estimates, which are necessary to evaluate the model in practice.

Cosmological simulations indicate w_f ∼ 10⁰ − 10¹ Mpc, while l ∼ 10¹ − 10² Mpc. The ratio of these parameters is constrained too: Table E.1 shows that a VFF comparison between the geometric model and cosmological simulations favours .

From Fig. 3, we estimated 𝒞 ∼ 10⁵ − 10⁷ for ν_ref = 150 MHz, assuming filament environments are characterised by , and cluster environments are characterised by (at the ∼3 Mpc h⁻¹ resolution, at least). Technically, Fig. 3 does not show P(j_ν | ρ) for z = 0 (the redshift at which 𝒞 is defined to be the cluster-to-filament SCW MEC ratio), but for z = 0.025. However, for our purposes, the evolution of this relationship is likely of negligible importance.

To find appropriate values for α_f and α_c, we calculated the MEC-weighted mean integrated spectral index from our Enzo simulation snapshot as a function of shock environment. Figure 7 shows the result when weighing by MECs at emission frequency ν = 150 MHz; in general, . Nevertheless, it seems that α_f = α_c = −3/2 are reasonable choices.

Fig. 7. Monochromatic-emission-coefficient–weighted mean integrated spectral index at an emission frequency of ν = 150 MHz, including all shocks that occur at locations where the dark matter (DM) density ρDM is below some threshold. Larger thresholds incorporate shocks from a wider range in ρDM, eventually including shocks from both filament and cluster environments. To provide an idea of the uncertainty in , we use jackknife realisations that each contain 80% of the shocked voxels present in our Enzo snapshot at z = 0.025.

Next, to find β_f and β_c, we revisited the single-shock MEC expression of Eq. (2). Inverse Compton (IC) scattering to the CMB contributes a factor (1+z)⁻⁴ to the single-shock MEC as and B_CMB ∝ (1+z)² (see Appendix B). Further factors of (1+z) follow by considering the typical redshift evolution of the proper BM density ρ_BM, proper magnetic field strength B, proper shock velocity v relative to the upstream IGM and comoving shock number density; all for both filament and cluster environments. In the linear regime of structure formation within an Einstein–de Sitter universe, the density contrast is proportional to (1+z)⁻¹. As the proper mean matter density is proportional to (1+z)³, the proper BM density in filaments is expected to be proportional to (1+z)². In clusters, the density contrast grows more strongly; numerical simulations suggest that the density contrast is roughly proportional to (1+z)⁻². Likewise, this would imply that the proper BM density in clusters is proportional to (1+z). Assuming that the magnetic field strength only evolves due to field line compression, we have (see Appendix B)¹⁸. In filaments then, B_f ∝ (1+z)^4/3, whilst in clusters B_c ∝ (1+z)^2/3. As the single-shock MEC is proportional to B^1 − α, the proper magnetic field strength approximately contributes a factor (1+z)^10/3 for filaments, and a factor (1+z)^5/3 for clusters. Lacking further knowledge, we applied Occam’s razor and assumed that the proper shock velocity and comoving shock number density do not evolve with redshift.

Collecting factors of (1+z), we found β_f = 4/3 and β_c = −4/3. We stress that these values are highly uncertain and should be used as indications only.

4. Results

4.1. MASSCW specific intensity function distribution: general results

All the results given are for ν_obs = 150 MHz. At the BORG SDSS LSS reconstruction resolution, realisations from our MASSCW specific intensity function distribution exhibit a factor ∼10³ of specific intensity variation over the sky, stemming from the highly variable power and localised nature of merger and accretion shocks in the Cosmic Web. These variations clearly appear in Fig. 8a, where we show a single realisation of over the full lune. Some sharp-edged structures are visible: these are due to the voxelised nature of the BORG SDSS density field samples. The effect is most pronounced for high-MEC voxels close to the observer, which result in bright patches that span hundreds of square degrees.

Fig. 8. Realisation of the merger- and accretion-shocked synchrotron Cosmic Web (MASSCW) priors at observing frequency νobs = 150 MHz. We show exactly 1/4 of the total sky. Top: specific intensity Iν. Bottom: specific-intensity–weighted mean (cosmological) redshift .

In Fig. 8b, we show the corresponding MASSCW specific-intensity–weighted mean redshift function realisation . All redshifts are within [0, 0.2], because the BORG SDSS reconstructions are limited to z_max = 0.2. A comparison with Fig. 8a reveals that generally, sky patches of high MASSCW specific intensity are due to structures at low redshift, and vice versa.

In Figs. 9a,b, we show the mean specific intensity function and mean specific-intensity–weighted mean redshift function at ν_obs = 150 MHz for each of the 10⁶ ray-traced directions, pooling all 10³ realisations. Both of these functions represent a summary statistic calculated from each ray’s marginal I_ν- and -distribution (that is the distribution obtained by marginalising out – from the joint (prior) distribution – the RVs of all directions but one). We calculated the mean specific intensity after removing, for each marginal distribution separately, the lowest 1% and highest 1% of values.

Fig. 9. Per-direction mean of the MASSCW priors at observing frequency νobs = 150 MHz, exhibiting less variability than the single realisation of Figs. 8a,b (note the smaller colour ranges). Three famous galaxy clusters stand out; these are (in increasing order of declination) the Hercules Cluster, the Coma Cluster, and Abell 2199. Top: specific intensity mean. Bottom: specific-intensity–weighted mean (cosmological) redshift mean.

Three bright spots stand out; these are (in increasing order of declination) the Hercules Cluster, the Coma Cluster, and Abell 2199. Note that we have not included the specific intensity contribution of radio halos around galaxy clusters, which are of different origin: turbulent reacceleration (Brunetti et al. 2001; Petrosian 2001). As observations suggest that the radio halo contribution usually dominates over the merger and accretion shock contribution, our results cannot be directly compared to actual galaxy cluster images.

The median specific intensity function (not shown) strongly resembles the mean specific intensity function, but is smaller over the whole lune. This is indicative of skewed (single-direction) marginal specific intensity distributions. The 5th to 95th percentile mean-to-median ratios span the interval (1.3, 2.1); the median mean-to-median ratio is 1.5.

Probabilistic approaches also allow for quantification of prediction uncertainty. In Figs. A.1a,b, we present both an absolute and relative measure of spread, again calculated from the MASSCW specific intensity prior marginals. Figure A.1a shows the standard deviation after performing the same filtering as was done for the mean. The resulting function is highly similar to the mean: directions that are brighter on average also tend to have larger absolute prediction uncertainties. Of natural interest is also the relative prediction uncertainty , which we show in Fig. A.1b. In this sense, the specific intensity function can be most accurately predicted around galaxy clusters and superclusters, reaching σ_Iν ∼ 60% μ_Iν. As expected, the relative prediction uncertainty is highest outside the SDSS DR7 coverage (i.e. near the ϕ = 90° and ϕ = 270° edges of the lune). It is also high in regions within the SDSS DR7 coverage that have low average brightness, inverting the trend that characterises the absolute prediction uncertainty.

4.2. MASSCW specific intensity function distribution: outskirts of massive galaxy clusters

Upon approaching the virial radius of a galaxy cluster from a connected filament, the WHIM transitions into the ICM, and the IGM’s magnetic field grows stronger. Especially in models such as Hoeft & Brüggen (2007), in which the synchrotron-emitting electrons originate from the thermal pool, these cluster outskirts constitute the most promising targets to find synchrotron emission from the filament IGM. Therefore, we also show zoom-ins of three realisations from our MASSCW specific intensity function prior near three massive galaxy clusters: the Coma Cluster, the Hercules Cluster, and Abell 2199. In Fig. 10, we show I_ν over spherical domes with an angular radius of 12°, whilst Fig. A.2 shows the specific-intensity–weighted mean redshift for the same realisations and sky regions. For the Coma Cluster, two filaments are discernible: one in northeastern, and one running in western direction. For the Hercules cluster, the realisations suggest one northern and two southbound filaments. Finally, for Abell 2199 a southbound filament is evident.

Fig. 10. Three realisations (rows) of the merger- and accretion-shocked synchrotron Cosmic Web (MASSCW) specific intensity prior at observing frequency νobs = 150 MHz, showing zoom-ins around three massive Northern Sky galaxy clusters (columns). The zoom-ins show caps of the celestial sphere of 10° radius. The dashed circles are at 3°, 6° and 9° from the cluster centre. Note that the usual radio halos that permeate galaxy clusters are not shown; these are caused by turbulent reacceleration, and we only show the merger and accretion shock contribution. The ∼3 Mpc h−1 resolution of the 3D total matter density and monochromatic emission coefficient (MEC) fields limits the resolution of the 2D specific intensity fields, so that individual shocks and ∼1 Mpc MEC – total matter density (anti)correlations are not resolved.

By inspecting the mean redshift functions, one can verify that the identified structures indeed lie at the cluster redshift, rather than being structures closer by or further away that appear connected to the clusters in chance alignments.

4.3. MASSCW specific intensity function distribution: deep LOFAR HBA fields

With the LOFAR HBA, several deep (≥50-h) observations have been conducted that complement the wide-field, 8-h approach of the LOFAR Two-metre Sky Survey (LoTSS) (Shimwell et al. 2019). Under thermal-noise–limited conditions, deep fields are of prime interest to search for a signature of the filament SCW. For this reason, in Fig. 11, we show the MASSCW specific intensity prior single-direction (i.e. marginal) medians for three such deep fields¹⁹. We generate imagery on arbitrary Flexible Image Transport System (FITS) grids by first calculating the marginal medians for each traced ray, and then applying a Voronoi tessellation to achieve a prediction for every pixel.

Fig. 11. Merger- and accretion-shocked synchrotron Cosmic Web (MASSCW) specific intensity prior marginal medians at observing frequency νobs = 150 MHz for three deep LOFAR HBA fields. Individual shocks should not be discernible in these statistical aggregates. The colour scales share the same lower bound. Top: Lockman Hole. Middle: Abell 2255. Bottom: Ursa Major Supercluster.

The Lockman Hole (for which a 100-h dataset is available Tasse et al. 2021) is a field known for its relatively low Milky Way column densities of neutral hydrogen and dust (Lockman et al. 1986), making it favourable for study in the extreme UV and soft X-ray bands, amongst others. The potential for multi-wavelength synergy has made deep observations of the Lockman Hole a LOFAR HBA priority (Mahony et al. 2016). Abell 2255 (for which a 75-h dataset is available Botteon et al., in prep.) is a merging galaxy cluster (Feretti et al. 1997; Botteon et al. 2020a) at z = 0.08 that is part of the North Ecliptic Pole (NEP) Supercluster (Mullis et al. 2001; Shim et al. 2011). In this dynamic environment, merger and accretion shocks could light up the filaments connected to Abell 2255 (Pizzo et al. 2008). Lastly, the Ursa Major Supercluster (for which a 50-h dataset is being assembled) at z = 0.06 stands out as the most prominent structure in our MASSCW specific intensity prior after the galaxy clusters shown in Figs. 10 and A.2. For this reason, it has been selected for deep LOFAR HBA observations.

4.4. Specific-intensity–weighted mean redshift

In the top panel of Fig. 12, we put our geometric model to the test by comparing a predicted specific-intensity–weighted mean redshift CDF (green dotted line) to 10² randomly drawn ECDFs from our 10³ MASSCW prior realisations (translucent black lines), at ν_obs = 150 MHz. Both prior realisations and the model reach a cumulative probability (‘fraction of the sky’) of 1 at : the BORG SDSS samples do not feature LSS beyond z_max = 0.2 and the model is restricted accordingly. The model CDF is constructed from 10⁵ draws from RV (see Eq. (20))²⁰. We adopt parameters suggested by Sect. 3.5: w_f = 5 Mpc, l = 50 Mpc, d₁ = 25 Mpc, 𝒞 = 10⁶ at ν_ref = 150 MHz, α_f = α_c = −3/2, β_f = 4/3 and β_c = −4/3. Without further parameter tuning, the model CDF reproduces the trend revealed by the majority of prior ECDFs. This correspondence is evidence that the geometric model in general (and not just for some highly specific choice of parameters) captures the main features of the MASSCW, and motivates calculating the distribution of for the true sky, which features LSS beyond z_max = 0.2. For the observational study of filaments, it is of greatest interest to calculate when the cluster contribution to the MASSCW is ignored, as it would otherwise dominate. In the bottom panel of Fig. 12, we therefore show according to the geometric model for z_max = ∞ and 𝒞 = 0. When 𝒞 = 0, α_c and β_c become irrelevant. Four parameters remain: α_f + β_f, l, w_f and d₁. The first of these has by far the most effect on , as evidenced by the differently coloured curves (we keep α_f = −3/2 constant and vary β_f, but varying α_f and keeping β_f constant would lead to identical results). We also vary the purely geometric parameters l, w_f and d₁, but in such a way that (as suggested by comparing the model’s VFFs to those from cosmological simulations) and (i.e. the observer is always put at the centre of a cubic unit cell). When l ∈ {50 Mpc, 75 Mpc, 100 Mpc} is increased, and w_f and d₁ accordingly, modest changes in occur (dash-dotted and dotted lines): the distribution of attains larger spread.

Fig. 12. MASSCW specific-intensity–weighted mean redshift RV results. Top: distributions of , showing 100 randomly selected ECDFs from our prior (grey curves), and a geometric model CDF (green curve). Both data and model consider LSS up to redshift zmax = 0.2 only. The other geometric model parameters are: wf = 5 Mpc, l = 50 Mpc, d1 = 25 Mpc, 𝒞 = 106, αf = αc = −3/2, βf = 4/3 and βc = −4/3. Bottom: geometric model CDFs of for filaments only (𝒞 = 0), though including all of them (zmax = ∞). By far the most influential parameter is also the most uncertain one; we therefore vary βf over a plausible range. We also vary l ∈ {50 Mpc, 75 Mpc, 100 Mpc}, and by extension also wf and d1, by forcing and . As l increases, the unit cell number density decreases; we symbolise this with more sparsely drawn curves. We keep αf = −3/2 constant.

provides detailed information. For example, the model with β_f = 4/3 and l = 50 Mpc predicts that for ∼50% of the sky (∼21 000 sq. deg), the MASSCW signal has a mean redshift of 0.35 or lower, whilst for ∼80% of the sky (∼33 000 sq. deg), the mean redshift is 0.42 or lower.

4.5. Flux-density–weighted mean redshift

The flux-density–weighted mean redshift constitutes the most concise measure of the typical MASSCW signal epoch. In the top panel of Fig. 13, we show a distribution (shaded grey) over generated via KDE (Gaussian kernel, σ_KDE = 4 × 10⁻⁴) from our 10³ MASSCW prior samples. The median is (solid grey); however, we stress that this is because only LSS up to z_max = 0.2 is included. We present an example model (solid green) that reproduces the median , with main parameters z_max = 0.2, l = 75 Mpc, d₁ = 60 Mpc, α_c = −3/2 and β_c = −4/3. (Adopting 𝒞 = 10⁶ at ν_ref = 150 MHz as suggested by Sect. 3.5, clusters dominate over filaments, and the other parameters play a very minor role.) To explore the sensitivity of on the parameters, we generate 5 × 10⁵ parameter sets by drawing from wide uniform distributions. We draw l ∼ Uniform(50 Mpc,100 Mpc), d₁ ∼ Uniform(0,l), w_f ∼ Uniform(5 Mpc,10 Mpc), log₁₀(𝒞) ∼ Uniform(5,7), α_f, α_c ∼ Uniform(−3/2−1/4,−3/2+1/4), β_f ∼ Uniform(4/3 − 2, 4/3 + 2) and β_c ∼ Uniform(−4/3−2,−4/3+2). However, the resultant distribution for is mostly restricted to the relatively small range of (0.06,0.08). We stress that the two distributions are conceptually distinct and are not meant to be compared directly.

Fig. 13. MASSCW flux-density–weighted mean redshift results. Top: distribution of calculated via KDE from our prior’s 103 realisations of Iν and (shaded grey), which are limited to zmax = 0.2. The median is (solid grey). Furthermore, we show the variety of geometric model predictions for when similarly limited to zmax = 0.2, assuming flat priors on all model parameters (shaded green; see text for details). The median is . As a concrete example, we show the flux-density–weighted mean redshift of a single model that reproduces the data median (solid green). This model is given by l = 75 Mpc, d1 = 60 Mpc and βc = −4/3. Bottom: predictions for considering filaments only, varying all three relevant geometric model parameters: αf + βf, l and d1. The two d1 = 50% l curves are virtually indistinguishable.

In the bottom panel of Fig. 13, we show for filaments only (𝒞 = 0) and z_max = ∞. The prediction then becomes a function of three parameters: α_f + β_f, l and d₁. The dependencies on both l and d₁ are weak, and is thus, to good approximation, determined by just a single parameter. Unfortunately, this parameter is currently highly uncertain, and should be constrained with upcoming cosmological simulations. The various coloured curves suggest that the impact of d₁ on increases when l is larger. As the unit cells grow (so that their number density drops ∝ l⁻³), becomes more sensitive to contributions of individual boundary crossings (filaments), including the first. Filament specific intensity contributions decrease quickly (especially for low values of α_f + β_f) and monotonically with distance, so that the distance of the first filament to the observer (e.g., nearby, d₁ = 10% l, or far away, d₁ = 90% l) is of some importance.

5. Discussion

5.1. Inter-cluster radio bridges: evidence for turbulence rather than shocks

In the last few years, LOFAR HBA observations around 140 MHz have revealed massive structures connecting clusters of galaxies at the onset of merging. These ‘bridges’, or ‘ridges’, likely are compressed filaments of the Cosmic Web. The first tentative discovery (which has since been confirmed) of an inter-cluster bridge was in Abell 1758 (of projected length ∼2 Mpc) (Botteon et al. 2018, 2020b), whilst the first unambiguous detection was in Abell 399–401 (of projected length ∼3 Mpc) (Govoni et al. 2019). The latter bridge features higher densities (∼3 × 10² m⁻³), temperatures (∼7 × 10⁷ K), and magnetic field strengths (∼1 μG) than are expected to exist in typical filaments.

These two inter-cluster bridges cannot be, and are not expected to be, faithfully reproduced by our MASSCW prior. Firstly, both bridges do not lie in the reconstructed volume: Abell 1758 lies at a redshift of z = 0.279, beyond the current redshift range (z < 0.2) of BORG SDSS reconstructions; Abell 399–401 has a favourable redshift of z = 0.07, but falls outside of the SDSS DR7 footprint. Another complication is that bridges of 2 − 3 Mpc extent (assuming a minor length component along the line-of-sight) are smaller than the current BORG SDSS resolution of ∼3 Mpc h⁻¹, and would therefore only barely be identifiable if they fell within the reconstructed volume.

Notwithstanding practical difficulties, Govoni et al. (2019) show that the Hoeft & Brüggen (2007) model alone cannot explain the LOFAR HBA data. Significantly, this discovery has cast doubt on the widespread hypothesis – also adopted in this article – that accretion shocks generate the dominant contribution to the filaments’ SCW signal. Other mechanisms, like turbulence (Brunetti & Vazza 2020), could play an important role. Cosmological simulations that model turbulence in the sense of (Brunetti & Vazza 2020) that would allow a recalculation of the SCW prior under this alternative emission mechanism, do not yet exist. We stress that our methods, although presently employed to generate MASSCW predictions, are of general nature, and can be used to explore alternative (combinations of) SCW emission mechanisms in the future.

5.2. Generalisation to other wavelength bands

We have presented our methodology in the context of predicting the Cosmic Web’s contribution to the radio sky. However, the methodology does not rely on any special property of radio emission. Apart from comparing different emission mechanisms for the same wavelength band (as discussed in the previous subsection), the methodology could also be extended to predicting the Cosmic Web’s contribution to the sky in other wavelength bands. With missions such as XRISM and Athena on the horizon, an extension into the X-ray window is of prime interest. As a demonstration, we have tentatively calculated the thermal bremsstrahlung component of the Cosmic Web’s contribution to the X-ray sky. For plain simplicity, we do not invoke snapshots from cosmological simulations, but rather assume the proportionalities and (with T_ref = 10⁶ K and ρ_BM, ref = m_p m⁻³, where m_p is the proton rest mass). Furthermore, we assume a hydrogen–helium plasma of primordial chemical composition, and use the Gaunt factor regimes from Novikov & Thorne (1973). We show results in Fig. A.3, for exactly the same sky regions as in Fig. 11. Note that although the existence of a thermal bremsstrahlung component from the WHIM to the X-ray sky is uncontested and features relatively low uncertainties, it is likely that the oxygen line emission component dominates. We have not considered this component due to its relative complexity and distance to our work’s main focus, but propose its calculation as a promising direction for future research.

5.3. Total matter density reconstructions: resolution, coverage and depth

The voxelised nature of density reconstructions like the BORG SDSS causes large-scale (∼10°) blocky shapes in specific intensity and specific-intensity–weighted mean redshift function realisations that do not represent plausible real-life morphologies. The angular scales at which these discontinuities occur, depend on the distances between the responsible (high-MEC) voxels and the observer. Future density reconstructions that are run at higher resolution will contain fewer of these problematic patterns. Furthermore, as hinted at in the previous paragraph, a modest resolution improvement (i.e. by a factor of order unity) would allow MASSCW predictions to contain inter-cluster bridges.

An improved resolution would also give relevance to the method of generating spatially correlated MEC fields as described in Sect. 2.4. Currently, the BORG SDSS voxel length is larger than the shock correlation length l_SE = 2 Mpc, so that the MEC fields obtained by our method are only marginally different from MEC fields where each voxel’s MEC is drawn independently from the appropriate conditional probability distribution.

A promising idea that is not pursued in this work, is to interpolate the density fields (by e.g., doubling or tripling the number of voxels along each dimension), so that it becomes possible to generate merger and accretion shocks of the appropriate size. Note that this approach would not add new small-scale structure to the density reconstructions, but merely ensures that the size of shocks generated stochastically are determined by the length scale l_SE, instead of by the density reconstruction resolution.

A straightforward next step would be to use density reconstructions of the remaining 3/4 of the sky to complete the MASSCW prior: generating predictions for the Southern Sky is relevant for SKA searches of the filament SCW. This could be done with the already available BORG 2M++ (Jasche & Lavaux 2019), although these reconstructions remain shallower than those of the BORG SDSS, with z_max = 0.1 instead of z_max = 0.2.

Reconstructions that push beyond z_max = 0.2 would also improve our predictions. The yield of such an extension depends on the filament-only specific-intensity–weighted mean redshift CDF for z_max = ∞. However, the single most influential parameter governing , α_f + β_f, is currently ill-constrained²¹. We therefore do not yet know to what extent deeper reconstructions would improve our MASSCW priors.

5.4. Cosmological simulations

The MASSCW predictions presented in this work are based on the statistical relationship between total matter density and the MASSCW MEC as inferred from Enzo cosmological simulations. However, in absence of tight observational constraints, these cosmological simulations must assume one of many magnetogenesis scenarios. The snapshots used in this work assume a primordial magnetogenesis scenario, starting from z = 45 with a seed magnetic field with a uniform comoving strength of 0.1 nG. More complex spectral energy distributions of primordial magnetic fields are however not excluded by present constraints from the CMB (e.g., Vazza et al. 2021 and references therein). Future work should explore the effect of different choices in the magnetogenesis scenario landscape on the MASSCW predictions, which could see a systematic change in specific intensity by an order of magnitude for filament-dominated directions.

A minor additional uncertainty comes from the fact that the (100 Mpc)³ cube used in this work is not yet large enough to fully capture the density–MASSCW MEC relation in a statistically exhaustive manner. This modest problem could be alleviated by appending the joint probability distribution from which the conditional probability distributions shown in Fig. 3 are derived with data from more simulation runs.

Finally, in this work, RAM limitations have necessitated discarding shocks with upstream Mach numbers below 2 in our shock identification procedure. This Mach number cut means that the MECs assigned to voxels are lower bounds. Future inclusion of the shocks now omitted could increase the MASSCW specific intensity functions by a (direction-dependent) factor of order unity.

5.5. Observational considerations

Using our predictions, we can evaluate the observational prospects of detecting the filament SCW – and the filament MASSCW in particular. As evinced by radio bridge detections (Botteon et al. 2018, 2020b; Govoni et al. 2019) and a statistical all-sky (or close to all-sky) detection (Vernstrom et al. 2021), both special-geometry and global observations of the filament SCW are already possible with modern low-frequency radio telescopes. The more interesting question therefore is whether detections are possible on an intermediate level – that of individual regions of large-scale structure – so that spatially resolved measurements of the intergalactic magnetic field strength B_IGM in filaments come within reach. In this section, we therefore explore the observational prospects and challenges of detecting the filament SCW around individual massive galaxy clusters, such as the Hercules Cluster, the Coma Cluster, Abell 2199 and Abell 2255, and around larger LSS complexes, such as the cluster triple in the Lockman Hole and the Ursa Major Supercluster. We assume that the MASSCW signal is the dominant contributor (in fact, the only contributor) to the SCW signal of these LSS regions²². For associated MASSCW predictions, see Figs. 10 and 11.

Making firm observability forecasts is in the first place hampered by the fact that there are three major unknowns in MASSCW predictions: firstly, the strength of the seed magnetic field in the Early Universe combined with the dominant process by which this field has evolved into magnetic fields in filaments today (see Sect. 1); secondly, the filling factor of shocks of appropriate strength and obliquity to trigger DSA; and, thirdly, the magnitude of the electron acceleration efficiency ξ_e in filaments (see Sect. 2.3.2). From Fig. 10, we see that the specific intensity function in the direction of filaments around massive low-redshift galaxy clusters, at ν_obs = 150 MHz and degree-scale resolution, reaches I_ν ∼ 10⁻¹ Jy deg⁻² for ξ_e = 1. If ξ_e ∼ 10⁻² in filaments, as in Keshet et al. (2004)’s SNR shocks, the actual specific intensity is I_ν ∼ 10⁻³ Jy deg⁻². A conservative estimate of the uncertainty in ξ_e is an order of magnitude, yielding a range of specific intensity estimates I_ν ∼ 10⁻⁴ − 10⁻² Jy deg⁻².

In this work, we have considered the MASSCW signal without contaminants. However, in most directions, the Milky Way is the dominant contributor to the specific intensity function at ν_obs  ∼ 10² MHz, and actual observational attempts to detect the SCW should therefore avoid the Galactic Plane, the North Polar Spur (Salter 1983), and other bright synchrotron Milky Way features. The extent to which the Milky Way hampers a SCW detection, depends on the typical angular scales of the synchrotron Milky Way, the SCW, and those measurable by the interferometer. Oei et al. (in prep.) have made compact-source–subtracted, low-resolution (60″ and 90″) images with the LoTSS DR2 (Shimwell et al. 2022) that reveal the Milky Way’s specific intensity function at ν_obs = 144 MHz up to degree scales. (By the lack of baselines shorter than 68 metres, larger scales are resolved out.) These images have most power on the degree scale and show specific intensity variations of σ_Iν ∼ 10⁰ − 10¹ Jy deg⁻² in the off-Galactic plane region. Meanwhile, the scales at which a sky region’s SCW specific intensity function has most power depend on the distances to the region’s most massive large-scale structures and vary per dominant emission mechanism assumed, with power on larger scales for turbulence compared to merger and accretion shocks.

A discrepancy in dominant scales between the specific intensity functions of the synchrotron Milky Way and the SCW can be leveraged to bolster the prospects of a SCW detection. Doing so appears important, because already at the 8-h depth of the LoTSS DR2, the Milky Way’s specific intensity dominates over the thermal noise. By using an inner (u, v)-cut that removes most Milky Way emission (but that retains most SCW emission), Milky Way–induced specific intensity variations can be reduced by one to two orders of magnitude²³, leaving a Milky Way contamination of σ_Iν ∼ 10⁻² − 10⁰ Jy deg⁻². (The optimal (u, v)-cut choice for a particular sky region can be derived by comparing the angular power spectrum of the synchrotron Milky Way with that of the region’s SCW predictions.) For depths such that the Milky Way remains the dominant noise source, the signal-to-noise ratio for a solid angle of a square degree centred around a massive filament thus is 10⁻⁴ in a very pessimistic case, 10⁻¹ in a fairly optimistic case, and 10⁰ in a very optimistic one. The signal-to-noise ratio grows with the square root of the number of such solid angles considered, and thus linearly with the angular diameter of the observed region. Notably, in the fairly optimistic case, we need to observe a region of 10¹ degree diameter to achieve a signal-to-noise ratio of order unity (for example three). For example, the Ursa Major Supercluster (see the bottom panel of Fig. 11) is a region of roughly the required extent. In case the MASSCW specific intensity is an order of magnitude weaker, that is of order I_ν ∼ 10⁻³ Jy deg⁻², the sky region required to detect the MASSCW is roughly the entire sky.

Galaxies populate filaments and generate synchrotron radiation. If they are not masked or removed from the observed imagery, their presence could mimic a SCW detection signal in cross-correlation experiments with low-spatial-resolution MASSCW predictions. (However, at high spatial resolution, Hodgson et al. 2021 show that the SCW and synchrotron emission from galaxies can be separated, as these signals trace the LSS matter distribution in different ways.) This issue underlines the importance of deep observations that feature low thermal noise levels, although – in contrast to the usual situation – in order to minimise systematic rather than random errors. As noted before, for LoTSS DR2 observations, Milky Way contamination dominates over thermal noise, rendering thermal noise largely irrelevant – at least prior to the suppression of degree scales. However, low thermal noise levels do allow for a more thorough subtraction of the galaxy population in filaments and thus help control an important systematic effect.

5.6. Independence of random fields

In this work, we have outlined how a specific intensity random field can be generated from a percentile random field 𝒫(r) (or, equivalently, a Gaussian random field 𝒵(r)) and a total matter density random field ρ(r); , ν_obs ∈ ℝ_{> 0}, and r ∈ ℛ. The function f that maps the two input random fields to the output random field is deterministic, informed by conditional probability distributions extracted from cosmological simulations, and non-linear in both arguments. Symbolically, I_ν = f(𝒫,ρ). Implicitly, our approach has been to sample ρ from the BORG SDSS posterior, and 𝒫 from another, independent distribution (which is fixed by the distribution of 𝒵, and thus specified by the covariance function K_SE). There is a physical scenario in which this is justified. The joint distribution for the input random fields P(𝒫,ρ) can be written as the product of a conditional and marginal: P(𝒫,ρ) = P(𝒫 | ρ)P(ρ). Thus, the central assumption that underlies our sampling approach is that P(𝒫 | ρ) ≈ P(𝒫), so that P(𝒫,ρ) ≈ P(𝒫)P(ρ). Our approximation is thus that the density field does not inform where, given a set of locations with the same density, high (or, equivalently, low) shock emission is more likely to occur.

The real world will violate this assumption to some extent. For example, a point in the outskirts of a galaxy cluster could be as dense as a point along the central axis of a prominent filament; however, the cluster point still likely has a different MASSCW MEC probability distribution than the filament point. One reason could be the presence of passing merger shocks in cluster outskirts; another could be the higher typical speed by which accretion shocks crash onto clusters, compared to the typical speed by which they hit filaments.

5.7. Spectral indices in 3D and 2D

In Sect. 2.5.1, we calculate I_ν by generating j_ν at different emission frequencies with the same percentile random field. This approach implicitly assumes that if shocks were to be ordered by their MEC, the ordering remains invariant over emission frequency range [ν_obs, ν_obs(1+z_max)]. Does this assumption correspond to a plausible physical scenario? Let the MECs of two shocks at some reference emission frequency ν = ν_ref be j₁ and j₂, and let j₁ > j₂ without loss of generality. Assume the existence of some function α = α(j_ν), that assigns integrated spectral indices to shocks based on their MEC at ν = ν_ref. At emission frequency ν, the MEC ranks of the shocks are the same as at ν_ref if and only if

(30)

The rewritten inequality

(31)

suggests α(j₂) − α(j₁) ≤ 0, or α(j₁) ≥ α(j₂) for ν > ν_ref: α(j_ν) must be a monotonically increasing function. Analogously, for ν < ν_ref, we find that α(j_ν) must be a monotonically decreasing function. These scenarios are visualised in Fig. A.4. Using the Enzo simulation data used in this work, we explore the spectral index - MEC relation at ν_ref = 180 MHz, and find a general downward trend for spectral index as a function of MEC at this frequency. However, the relation is scattery and therefore not fully described by a monotonically decreasing function α(j_ν); this implies that shock MEC ranks do change when varying the emission frequency. Although our approach to determining spectral indices remains approximate, compared to other uncertainties in our methodology, the error thus introduced is likely of minor importance²⁴.

Our formalism allows for the generation of the function at two (or more) different observing frequencies. In turn, this enables the calculation of spectral indices for specific intensity rather than MEC (i.e. ‘in 2D’ instead of ‘in 3D’), emulating the type of spectral analysis routinely performed by observational astronomers. We caution that the procedure for MEC spectral index assignment used in this work does not respect the full diversity of spectral behaviour present in the Enzo simulation, and instead forces MEC spectral indices to approach the MEC-weighted mean. In turn, this also causes specific intensity spectral index variations to be biased low. Future work should adapt this procedure so that a plausible specific intensity spectral index prior can be added to the potent suite of predictions that follow simultaneously from our methodology.

6. Conclusions

In this work, we describe and implement the first methodology to produce a (prior) probability distribution over specific intensity functions representing the synchrotron cosmic web (SCW) of the Local Universe. We assume merger and accretion shocks to be the main generators of the SCW, and assume a primordial magnetogenesis scenario for the evolution of magnetic fields in the IGM. However, the methodology is general enough to explore alternative physical hypotheses in the future. Our prior can be used to guide and verify observational attempts to detect the SCW with low-frequency radio telescopes such as the LOFAR and the SKA.

1. Using BORG SDSS total matter density reconstructions and Enzo cosmological simulations, we have built a prior distribution that is informative over half of the Northern Sky, and that has a ∼0.6° resolution for LSS at z = 0.1. Although not a fundamental limitation of the methodology, the current resolution is not high enough to resolve individual merger and accretion shocks. Typically, filaments near massive structures give at ν_obs = 150 MHz; ξ_e is the highly uncertain electron acceleration efficiency. Even at the ∼3 Mpc h⁻¹ reconstruction resolution, our MASSCW prior indicates that the specific intensity for a given direction is highly uncertain (with a typical standard deviation being ∼100% of the mean) due to uncertainty regarding the presence and highly variable nature of shock emission along the line-of-sight. We present (marginal median) MASSCW specific intensity predictions for three deep LOFAR HBA fields: the Lockman Hole, Abell 2255, and the Ursa Major Supercluster.

2. With a simple geometric model of cubic unit cells, we have calculated both the distribution of the specific-intensity–weighted mean redshift RV , as well as the flux-density–weighted mean redshift of the MASSCW signal for LSS reconstructions up to z_max = 0.2. We obtain results that closely resemble those found numerically from our data-driven MASSCW prior, whose construction is much more involved. Encouraged by this, we present filament-only geometric model predictions for that include all LSS (i.e. z_max = ∞). These predictions are highly insensitive to plausible variations in model parameters l and d₁, demonstrating that is effectively determined by a single parameter: the sum of the typical MEC-weighted filament spectral index α_f and the MEC – cosmological redshift power law exponent β_f. Its future determination will characterise the completeness of the MASSCW predictions put forth in this work.

In an optimistic case, our prior already reveals a great deal about filamentary baryons – and where to find them.

Acknowledgments

M.S.S.L. Oei and R.J. van Weeren acknowledge support from the VIDI research programme with project number 639.042.729, which is financed by the Netherlands Organisation for Scientific Research (NWO). F. Vazza acknowledges support from the ERC STG MAGCOW (714196) from the H2020. F. Leclercq acknowledges funding from the Imperial College London Research Fellowship Scheme. The cosmological simulations used in this work were produced with the Enzo code (enzo-project.org) and run on the Piz-Daint supercluster at LSCS (Lugano) under project ‘s701’ with F. Vazza as P.I.. M.S.S.L. Oei warmly thanks Andrea Botteon, Matthias Hoeft, Vincent Icke, Josh Albert, Lara Anisman, Jacob Bakermans and Jesse van Oostrum for helpful discussions. In dearest memory of Grandma Maria, who brightened my universe with so much love and laughter.

References

Appendix A: Additional figures

In this appendix, we list four figures referenced in the main text, but relegated for structural clarity.

Fig. A.1. MASSCW priors allow for a quantification of prediction uncertainty. Here we show both an absolute and a relative measure of spread for the single-direction specific intensity distributions (i.e. marginals) at νobs = 150 MHz. For each direction, we discard data outside the 1–99% percentile range. Top: marginal standard deviation (absolute uncertainty), which closely resembles the marginal mean of Fig. 9a. Bottom: marginal standard deviation over mean (relative uncertainty), which reveals an inverted trend.

Fig. A.2. Same as Fig. 10, but now for specific-intensity–weighted mean redshift instead of specific intensity.

Fig. A.3. Thermal bremsstrahlung specific intensity prior marginal medians at Eobs = 1 keV for three deep LOFAR HBA fields. For merger and accretion shock synchrotron predictions of the same regions, see Fig. 11. Top: Lockman Hole. Middle: Abell 2255. Bottom: Ursa Major Supercluster.

Fig. A.4. Three scenarios for simple power-law synchrotron spectra of shocks that are consistent with our methodology. Due to cosmological redshifting, the calculation of the specific intensity Iν at observing frequency νobs necessitates knowing the monochromatic emission coefficient (MEC) jν at a range of emission frequencies ν (see Eq. 8). In this work, we generate jν at emission frequencies ν = νobs and ν = νobs(1+zmax), where zmax is the maximum cosmological redshift of the large-scale structure reconstructions. We do so using the same Gaussian random field realisation 𝒵 (and thus percentile random field realisation 𝒫), implicitly assuming that shocks retain their MEC percentile rank over the frequency range [νobs, νobs(1+zmax)]. This assumption holds in the three scenarios sketched above. Each solid line represents a shock with constant percentile rank (over the range shown, at least); the graphs are drawn with logarithmic scaling. Regarding the sign of the spectral index, we use the convention jν ∝ να. Left: the integrated spectral index α is a monotonically increasing function of jν at ν = νobs (‘brighter shocks have flatter spectra’). Middle: the integrated spectral index α is a monotonically decreasing function of jν at ν = νobs(1+zmax) (‘brighter shocks have steeper spectra’). Right: the integrated spectral index α is a jν-independent constant (at all frequencies); this is a limiting case of the previous two.

Appendix B: Single-shock synchrotron MEC - total matter density scaling relation

In this appendix, we derive a synchrotron power density - total matter density scaling relation for individual shocks with high upstream Mach numbers in cluster outskirts and filaments, assuming the Hoeft & Brüggen (2007) model and γ = 5/3. Because the power density and MEC of a single shock are proportional, this immediately also yields the desired single-shock synchrotron MEC - total matter density scaling relation.

B.1. Temperature and the speed of sound

Note that , where v_u is the shock velocity relative to the upstream plasma, and c_s, u is the speed of sound in the upstream plasma. The Newton–Laplace equation for an ideal gas predicts , where T_u is the upstream plasma temperature.²⁵ The upstream Mach number of a shock incident on the WHIM would therefore be higher than that of a shock incident on the ICM if the shocks arrive at the same velocity v_u relative to these media.

B.2. The filament regime: low magnetic field strengths

One of the prime reasons for pursuing SCW detections is to gauge the unknown strength of the Universe’s largest magnetic fields. Numerical simulations by Vazza et al. (2015, 2017) that reproduce the observed magnetic field strengths in galaxy clusters, predict magnetic field strengths in filaments that depend strongly on the magnetogenesis scenario considered, ranging between 10⁻¹ − 10² nG. For the purposes of finding a power density - matter density scaling relation, the relevant quantity to compare the filament IGM magnetic field strength B_IGM with at cosmological redshift z, is the CMB magnetic field strength B_CMB(z).

As the CMB is well modelled by a blackbody²⁶, the CMB magnetic field strength B_CMB(z) is derived by equating the electromagnetic energy density u_EM of a blackbody of temperature T to the electromagnetic energy density of a magnetic field of magnitude B:

(B.1)

where σ is the Stefan–Boltzmann constant and μ₀ is the vacuum permeability. Upon rearranging, and for T = T_CMB and B = B_CMB, we find

(B.2)

Let a be the scale factor and let a₀ be its present-day value. As due to the expansion of the Universe, it follows that T_CMB ∝ 1 + z, and thus

(B.3)

Using T_CMB(0) = 2.725 K yields B_CMB(0) = 3.238 μG.

Thus, under all plausible scenarios of magnetogenesis, B_IGM(z) ≪ B_CMB(0) ≤ B_CMB(z).

B.3. Magnetic field strength and baryon density

A scaling relation between B and ρ_BM follows from considering the conservation of magnetic flux as the Universe expands. The magnetic flux through a surface is the product of the surface area and the magnetic field strength (and the cosine of the angle between the surface normal and the magnetic field). Over time, the surface area increases ∝a², so that the magnetic field strength must follow ∝a⁻² if conservation of magnetic flux is to hold. Finally, as ρ_BM ∝ a⁻³, one obtains . Fig. 4 of Vazza et al. (2017) compares this scaling relation with simulated magnetic field strengths and baryon densities under various scenarios of magnetogenesis.

B.4. The power density expression simplifies

The Hoeft & Brüggen (2007) power density folded into Eq. 2 simplifies appreciably if high Mach numbers and low magnetic field strengths are assumed. In such a regime, α ≈ −1 and Ψ ≈ 1, while . Thus,

(B.4)

From the Rankine–Hugoniot jump conditions, one can derive that the compression factor (Fig. 1, central panel), as a function of ℳ_u and γ, is

(B.5)

The same equations dictate that the temperature increase (Fig. 1, bottom panel), as a function of ℳ_u and γ, is

(B.6)

Thus, for ℳ_u ≫ 1, and . Returning to the power density scaling relation, we find

(B.7)

To arrive at the second line, we use the definition of the upstream Mach number and the upstream sound speed - upstream plasma temperature scaling relation.

B.5. Upstream velocity and total density

We investigate the upstream velocity - total (i.e. dark and baryonic matter) density relation for three simple geometries. We invoke Gauss’ law for gravity to find expressions for the gravitational field, derive the gravitational potential using the fact that the gravitational force is conservative, and equate, for a test particle, the loss in gravitational potential energy to the gain in kinetic energy. We assume that the structures considered have hard edges and are equidense (with total matter density ρ ≔ ρ_BM, u + ρ_DM) within.

Upon impact, the velocity of a test particle starting from rest at a distance d, and falling onto an isolated spherical galaxy cluster with radius R, is

(B.8)

Upon impact, the velocity of a test particle starting from rest at a distance d, and falling onto an isolated cylindrical filament with radius R, is

(B.9)

Upon impact, the velocity of a test particle starting from rest at a distance d, and falling onto an isolated, thick planar sheet of half-thickness R, is

(B.10)

In all three cases, R determines the size of the structure types. Assuming no relation between total density ρ and R, we find , irrespective of the geometry.

B.6. Baryon density and dark matter density

Structure formation theory predicts that after decoupling, and if gas pressure is ignored, ρ_BM, u ∝ ρ_DM. Including gas pressure, this proportionality is expected to remain valid on large scales only. For filaments in particular, Fig. 6 of Gheller et al. (2016) shows that the baryon fraction in filaments remains close to f_cosmic ≔ Ω_BM, 0(Ω_BM, 0+Ω_DM, 0)⁻¹(=0.167) over four orders of magnitude of total baryonic mass, and for redshifts from 1 to 0; this is consistent with ρ_BM, u ∝ ρ_DM.

Using , and ρ_BM, u ∝ ρ_DM, the final scaling relation becomes

(B.11)

Note that, to obtain a scaling relation for the total (rather than single-shock) MEC, one should also consider how the shock number density relates to the total matter density. The total MEC - matter density and the single-shock MEC - matter density scaling relation exponents are only the same when no such relationship exists.

Appendix C: Ray tracing in the cosmological setting

Projecting a ray’s 4D null geodesic in a pure Friedmann–Lemaître–Robertson–Walker (FLRW) metric onto 3D comoving space results in a straight line. This follows readily from the FLRW metric in hyperspherical coordinates

(C.1)

where ds² is the spacetime line element, c is the speed of light in vacuo, t is physical time since the Big Bang, and a(t) is the scale factor (with a₀ ≔ a(t₀) being its present-day (t = t₀) value). Also, r is the radial comoving distance, k is the Universe’s Gaussian curvature (with SI units m⁻²), is the transverse comoving distance and dΩ² ≔ dθ² + cos²θdϕ². (The sinc function follows the mathematical (unnormalised) convention.) Finally, let the location of present-day Earth be the spatial origin. An initially radial (dΩ² = 0) null (ds² = 0) geodesic thus satisfies

(C.2)

regardless of k. As time progresses, such rays maintain their direction and only change in r; this justifies considering the path ℒ of a light ray with direction in comoving space as the set of points .

Appendix D: Observer’s specific intensity

Our aim is to derive an expression for the specific intensity in direction at observing frequency ν_obs. We follow Chapter 12 of Peacock (1999), but generalise to arbitrary ΛCDM models (by allowing Λ ≠ 0), and recast the results in terms of the MEC instead of the emissivity.

As in Appendix C, consider a FLRW metric with arbitrary Gaussian curvature k. A comoving volume element dV_c and the corresponding proper volume element dV_p at comoving radial distance r and cosmological redshift z = z(r) that cover a solid angle dΩ on the sky are given by

(D.1)

(D.2)

Recall that we have defined j_ν as the proper MEC, and assume that the filament IGM radiates isotropically. The luminosity density dL_ν of the volume, seen in direction and at emission frequency ν = (1+z)ν_obs, then equals

(D.3)

The corresponding observer’s flux density dF_ν of the volume in direction at observing frequency ν_obs is

(D.4)

and so the observer’s specific intensity dI_ν in direction at observing frequency ν_obs is

(D.5)

(D.6)

We neglect absorption so that the specific intensity of the ray only accumulates as the ray travels through LSS to the observer: the Universe is mostly optically thin for ν < 1 GHz; we assume this holds perfectly. So, by collecting all contributions along the ray’s path, we obtain the equivalent of Eq. 12.12 of Peacock (1999):

(D.7)

Alternatively, one can express as an integral over z. Because

(D.8)

with the dimensionless Hubble parameter being and today’s curvature density parameter being Ω_K, 0 = 1 − Ω_R, 0 − Ω_M, 0 − Ω_Λ, 0, we have

(D.9)

Combining Eq. D.7, D.8 and D.9,

(D.10)

(Barring notational differences, the Ω_Λ, 0 = 0 (and Ω_R, 0 = 0) limit of this formula is Eq. 12.10 of Peacock (1999).) Finally, we can read off that

(D.11)

Appendix E: Volume-filling fractions

In this appendix, we compute the volume-filling fractions (VFFs) of the four canonical structure types (clusters, filaments, sheets, and voids) as predicted by the cubic unit cell geometric model developed in Sect. 3. A single cubic unit cell features two typical lengthscales: a large scale (l_f), and a small scale (w_f), which can be interpreted as the typical filament length and width, respectively. The cube obtained by raising the large scale to the third power represents a void. Similarly, the three rectangular cuboids obtained by taking the product of the square of the large scale, and the small scale, resemble three sheets. The three rectangular cuboids obtained by taking the product of the large scale, and the square of the small scale, resemble three filaments. Finally, the cube obtained by raising the small scale to the third power resembles a cluster. The natural volume-filling fractions (VFFs) suggested by this geometry are thus

(E.1)

(E.2)

(E.3)

(E.4)

See Fig. E.1. To see that the VFFs sum to 1, we rewrite 1 using the binomial theorem:

(E.5)

Fig. E.1. Cosmic web structure type volume-filling fractions (VFFs) as predicted by the cubic unit cell model of Sect. 3. These VFFs depend on just one parameter: the ratio between the (comoving) filament width and the (comoving) cubic unit cell edge length .

and recognise the VFFs as the four terms in this expansion.

The VFFs obtained from this one-parameter model for , are similar to those retrieved by Forero-Romero et al. (2009) from cosmological simulations for eigenvalue threshold λ_th = 1 and effective smoothing scales R_eff ∼ 10⁰ Mpc (see Table E.1).

Appendix F: Notation

This paper adopts the SI system of units (for formulae concerning electromagnetism), and the following symbols. We list dimensionalities in SI base units with the radian ‘rad’ appended. Current-day quantities are subscripted with a zero: for example a₀ is today’s scale factor, while a is the scale factor for arbitrary times. Upstream and downstream quantities are subscripted with a ‘u’ or ‘d’, respectively: for example T_u is the upstream plasma temperature, whilst T is the general plasma temperature.

All Tables

All Figures

Fig. 1. Dependence of three central (merger or accretion) shock quantities on upstream Mach number ℳu and adiabatic index γ, as derived from the Rankine–Hugoniot jump conditions in ideal gases. Radiation from ℳu < 10 shocks is predominantly generated by electrons that stem from the ICM, whilst radiation from ℳu > 10 shocks is predominantly generated by electrons that stem from the WHIM. Top: the spectral index α of angularly unresolved shocks, assuming standard DSA. Middle: the downstream-over-upstream plasma density ratio. Bottom: the downstream-over-upstream plasma temperature ratio.

In the text

Fig. 2. Two simulated specific intensity functions at νobs = 150 MHz, assuming synchrotron radiation from merger and accretion shocks (as in Hoeft & Brüggen 2007) in LSS at z = 0.045. Top: primordial scenario for magnetogenesis, starting from BIGM = 1 nG in the Early Universe. (In contrast, the simulation underlying this article’s predictions starts at BIGM = 0.1 nG.) Bottom: galactic scenario for magnetogenesis, in which seeding occurs through AGN outflows and supernova winds. See Gheller & Vazza (2020) for details on these scenarios.

In the text

Fig. 3. Conditional probability distributions P(jν | ρ) of synchrotron monochromatic emission coefficient jν (at rest-frame (emission) frequency ν = 150 MHz) given total matter density ρ for the Hoeft & Brüggen (2007) formalism, as derived from z = 0.025 snapshots of cosmological MHD simulations by Vazza et al. (2019). The dashed lines indicate a jν ∝ ρ23/6 (single shock) scaling relation for various proportionality constants. ρ is shown relative to the current-day mean total matter density , while the weakly constrained electron acceleration efficiency ξe is divided out from jν. Each conditional is numerically approximated by a probability mass function with 1000 bins, whose edges vary by a constant factor and span 40 orders of magnitude in jν.

In the text

Fig. 4. Overview of the procedure for converting a probability distribution over total (i.e. baryonic and dark) matter density fields into a probability distribution over monochromatic emission coefficient (MEC) fields. We show a fixed slice through ℛ, of 750 Mpc h−1 × 750 Mpc h−1 for the density and MEC fields, and a subregion of 150 Mpc h−1 × 150 Mpc h−1 for the Gaussian and percentile random fields (to more clearly illustrate the bijective mapping between the two). Left: BORG SDSS density ρ sample, in multiples of the current-day mean total matter density . Middle left: Gaussian random field 𝒵 sample using a squared exponential (SE) kernel with lengthscale lSE = 2 Mpc. Middle right: the corresponding percentile random field 𝒫 sample. Right: the corresponding MEC jν sample at an emission frequency of ν = 150 MHz, with the weakly-constrained electron acceleration efficiency ξe divided out.

In the text

Fig. 5. Geometry of the simple model, in which we treat the cosmic web as a collection of randomly displaced cubic unit cells. We show one face (grey square) of a unit cell, with edges of comoving length l. Filaments of typical comoving width wf and length lf ≔ l − wf cover the light-blue–shaded region of the face, and extend to the faces of neighbouring unit cells (grey dashed lines). Where filaments meet, galaxy clusters reside (orange-shaded regions), each of which is also connected to a filament oriented away from the observer. The central part of the face represents a sheet (dark blue), behind which a large void looms (not depicted). The figure is not to scale: filament widths are exaggerated compared to realistic filament lengths.

In the text

	Fig. 6. When a sightline enters a cubic unit cell of the cosmic web perpendicularly to a face (as in Fig. 5) and at a random position on it, one of three different LSS piercing events occurs. We show the probabilities for each event as a function of the ratio of the typical filament width and cubic unit cell edge length (see Eqs. (11)–(13)).
In the text

Fig. 7. Monochromatic-emission-coefficient–weighted mean integrated spectral index at an emission frequency of ν = 150 MHz, including all shocks that occur at locations where the dark matter (DM) density ρDM is below some threshold. Larger thresholds incorporate shocks from a wider range in ρDM, eventually including shocks from both filament and cluster environments. To provide an idea of the uncertainty in , we use jackknife realisations that each contain 80% of the shocked voxels present in our Enzo snapshot at z = 0.025.

In the text

	Fig. 8. Realisation of the merger- and accretion-shocked synchrotron Cosmic Web (MASSCW) priors at observing frequency νobs = 150 MHz. We show exactly 1/4 of the total sky. Top: specific intensity Iν. Bottom: specific-intensity–weighted mean (cosmological) redshift .
In the text

Fig. 9. Per-direction mean of the MASSCW priors at observing frequency νobs = 150 MHz, exhibiting less variability than the single realisation of Figs. 8a,b (note the smaller colour ranges). Three famous galaxy clusters stand out; these are (in increasing order of declination) the Hercules Cluster, the Coma Cluster, and Abell 2199. Top: specific intensity mean. Bottom: specific-intensity–weighted mean (cosmological) redshift mean.

In the text

Fig. 10. Three realisations (rows) of the merger- and accretion-shocked synchrotron Cosmic Web (MASSCW) specific intensity prior at observing frequency νobs = 150 MHz, showing zoom-ins around three massive Northern Sky galaxy clusters (columns). The zoom-ins show caps of the celestial sphere of 10° radius. The dashed circles are at 3°, 6° and 9° from the cluster centre. Note that the usual radio halos that permeate galaxy clusters are not shown; these are caused by turbulent reacceleration, and we only show the merger and accretion shock contribution. The ∼3 Mpc h−1 resolution of the 3D total matter density and monochromatic emission coefficient (MEC) fields limits the resolution of the 2D specific intensity fields, so that individual shocks and ∼1 Mpc MEC – total matter density (anti)correlations are not resolved.

In the text

	Fig. 11. Merger- and accretion-shocked synchrotron Cosmic Web (MASSCW) specific intensity prior marginal medians at observing frequency νobs = 150 MHz for three deep LOFAR HBA fields. Individual shocks should not be discernible in these statistical aggregates. The colour scales share the same lower bound. Top: Lockman Hole. Middle: Abell 2255. Bottom: Ursa Major Supercluster.
In the text

Fig. 12. MASSCW specific-intensity–weighted mean redshift RV results. Top: distributions of , showing 100 randomly selected ECDFs from our prior (grey curves), and a geometric model CDF (green curve). Both data and model consider LSS up to redshift zmax = 0.2 only. The other geometric model parameters are: wf = 5 Mpc, l = 50 Mpc, d1 = 25 Mpc, 𝒞 = 106, αf = αc = −3/2, βf = 4/3 and βc = −4/3. Bottom: geometric model CDFs of for filaments only (𝒞 = 0), though including all of them (zmax = ∞). By far the most influential parameter is also the most uncertain one; we therefore vary βf over a plausible range. We also vary l ∈ {50 Mpc, 75 Mpc, 100 Mpc}, and by extension also wf and d1, by forcing and . As l increases, the unit cell number density decreases; we symbolise this with more sparsely drawn curves. We keep αf = −3/2 constant.

In the text

Fig. 13. MASSCW flux-density–weighted mean redshift results. Top: distribution of calculated via KDE from our prior’s 103 realisations of Iν and (shaded grey), which are limited to zmax = 0.2. The median is (solid grey). Furthermore, we show the variety of geometric model predictions for when similarly limited to zmax = 0.2, assuming flat priors on all model parameters (shaded green; see text for details). The median is . As a concrete example, we show the flux-density–weighted mean redshift of a single model that reproduces the data median (solid green). This model is given by l = 75 Mpc, d1 = 60 Mpc and βc = −4/3. Bottom: predictions for considering filaments only, varying all three relevant geometric model parameters: αf + βf, l and d1. The two d1 = 50% l curves are virtually indistinguishable.

In the text

Fig. A.1. MASSCW priors allow for a quantification of prediction uncertainty. Here we show both an absolute and a relative measure of spread for the single-direction specific intensity distributions (i.e. marginals) at νobs = 150 MHz. For each direction, we discard data outside the 1–99% percentile range. Top: marginal standard deviation (absolute uncertainty), which closely resembles the marginal mean of Fig. 9a. Bottom: marginal standard deviation over mean (relative uncertainty), which reveals an inverted trend.

In the text

	Fig. A.2. Same as Fig. 10, but now for specific-intensity–weighted mean redshift instead of specific intensity.
In the text

	Fig. A.3. Thermal bremsstrahlung specific intensity prior marginal medians at Eobs = 1 keV for three deep LOFAR HBA fields. For merger and accretion shock synchrotron predictions of the same regions, see Fig. 11. Top: Lockman Hole. Middle: Abell 2255. Bottom: Ursa Major Supercluster.
In the text

Fig. A.4. Three scenarios for simple power-law synchrotron spectra of shocks that are consistent with our methodology. Due to cosmological redshifting, the calculation of the specific intensity Iν at observing frequency νobs necessitates knowing the monochromatic emission coefficient (MEC) jν at a range of emission frequencies ν (see Eq. 8). In this work, we generate jν at emission frequencies ν = νobs and ν = νobs(1+zmax), where zmax is the maximum cosmological redshift of the large-scale structure reconstructions. We do so using the same Gaussian random field realisation 𝒵 (and thus percentile random field realisation 𝒫), implicitly assuming that shocks retain their MEC percentile rank over the frequency range [νobs, νobs(1+zmax)]. This assumption holds in the three scenarios sketched above. Each solid line represents a shock with constant percentile rank (over the range shown, at least); the graphs are drawn with logarithmic scaling. Regarding the sign of the spectral index, we use the convention jν ∝ να. Left: the integrated spectral index α is a monotonically increasing function of jν at ν = νobs (‘brighter shocks have flatter spectra’). Middle: the integrated spectral index α is a monotonically decreasing function of jν at ν = νobs(1+zmax) (‘brighter shocks have steeper spectra’). Right: the integrated spectral index α is a jν-independent constant (at all frequencies); this is a limiting case of the previous two.

In the text

	Fig. E.1. Cosmic web structure type volume-filling fractions (VFFs) as predicted by the cubic unit cell model of Sect. 3. These VFFs depend on just one parameter: the ratio between the (comoving) filament width and the (comoving) cubic unit cell edge length .
In the text