11institutetext: Leiden Observatory, Leiden University, Niels Bohrweg 2,
NL-2300 RA Leiden, the Netherlands 22institutetext: Department of Physics and
Astronomy, University of Bologna, Via Gobetti 93/2, I-40129 Bologna, Italy
33institutetext: INAF, Istituto di Radioastronomia, Via Gobetti 101, I-40129
Bologna, Italy 44institutetext: Hamburg Observatory, Hamburg University,
Gojensbergweg 112, 21029 Hamburg, Germany 55institutetext: Imperial Centre for
Inference and Cosmology, Imperial College London, 1012 Blackett Laboratory,
Prince Consort Road, London SW7 2AZ, United Kingdom 66institutetext: Anton
Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904,
NL-1098 XH Amsterdam, the Netherlands

# Filamentary Baryons and Where to Find Them

A forecast of synchrotron radiation from merger and accretion shocks in the
local Cosmic Web

Martijn S.S.L. Oei E-mail<EMAIL_ADDRESS>_In dearest memory of
Grandma Maria, who brightened my universe with so much love and laughter._ 11
Reinout J. van Weeren 11 Franco Vazza 223344 Florent Leclercq 55 Akshatha
Gopinath 66 Huub J.A. Röttgering 11

(Received February 1, 12021 H.E.; accepted March 7, 12022 H.E.)

###### Abstract

Context. The detection of synchrotron radiation from the intergalactic medium
(IGM) that pervades the filaments of the Cosmic Web, constitutes an upcoming
frontier to test physical models of astrophysical shocks and their radiation
mechanisms, trace the missing baryons, and constrain magnetogenesis: the
origin and evolution of extragalactic magnetic fields.

Aims. At the moment of writing, the first synchrotron detections of the IGM
within filaments have been claimed. Now is the time to develop a rigorous
statistical framework to predict sky regions with the strongest signal, and to
move from mere detection to inference: the identification of the most
plausible physical models and parameter values from the observations.

Methods. Current theory posits that the filament IGM lights up through shocks
that originate from large-scale structure formation. With Bayesian inference,
we generate a probability distribution on the set of specific intensity
functions that represent our view of the merger- and accretion-shocked
synchrotron Cosmic Web (MASSCW). We combine the BORG SDSS total matter density
posterior, which is based on spectroscopic observations of galaxies within
SDSS DR7, snapshots of Enzo MHD cosmological simulations, Gaussian random
fields (GRFs) and a ray tracing approach to arrive at the result.

Results. We present a physics-based prediction of the MASSCW signal, including
principled uncertainty quantification, for a quarter of the sky and up to
cosmological redshift $z_{\mathrm{max}}=0.2$. The super-Mpc 3D resolution of
the current implementation limits the resolution of the predicted 2D imagery,
so that individual merger and accretion shocks are _not_ resolved. The MASSCW
prior can be used to identify the most promising fields to target with low-
frequency radio telescopes, and to conduct actual detection experiments. We
furthermore calculate a probability distribution for the flux-density-weighted
mean (i.e. sky-averaged) redshift $\bar{\bar{z}}$ of the MASSCW signal up to
$z_{\mathrm{max}}$, and find a median of $\bar{\bar{z}}=0.077$. We construct a
low-parametric analytic model that produces a similar distribution for
$\bar{\bar{z}}$, with a median of $\bar{\bar{z}}=0.072$. Extrapolating the
model, we can calculate $\bar{\bar{z}}$ for _all_ large-scale structure in the
Universe (including what lies beyond $z_{\mathrm{max}}$) and show that, if one
only considers filaments, $\bar{\bar{z}}$ depends on virtually one parameter.
As case studies, we finally explore the predictions of our MASSCW specific
intensity function prior in the vicinity of three galaxy clusters: the
Hercules Cluster, the Coma Cluster, and Abell 2199; and in three deep LOFAR
HBA fields: the Lockman Hole, Abell 2255, and the Ursa Major Supercluster.

Conclusions. We describe and implement a novel, flexible and principled
framework for predicting the low-frequency, low-resolution specific intensity
function of the Cosmic Web due to merger and accretion shocks that arise
during large-scale structure formation. The predictions guide Local Universe
searches for filamentary baryons through half of the Northern Sky. Once
cosmological simulations of alternative emission mechanisms have matured, our
approach can be extended to predict additional physical pathways that
contribute to the elusive synchrotron Cosmic Web signal.

###### Key Words.:

Cosmology: miscellaneous – dark matter – large-scale structure of Universe –
Galaxies: clusters: intracluster medium – intergalactic medium – Magnetic
fields – Magnetohydrodynamics (MHD) – Methods: numerical – statistical –
Radiation mechanisms: non-thermal – Radio continuum: general – Shock waves

##  1 Introduction

###  1.1 The Cosmic Web

Just after inflation, the Universe’s dark and baryonic matter density
functions resembled realisations of nearly constant, isotropic and stationary
Gaussian random fields (Linde, 2008). By the influence of 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.111Some authors use
‘Cosmic Web’ to refer exclusively to filaments, but in this work, the term is
used to refer to _all_ of the late-time Universe excluding the voids. We
furthermore differentiate between ‘cosmic web’ (for the concept in an
arbitrary universe), and ‘Cosmic Web’ (for the concept in ours). 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 magnetic field and spin. Hahn et al. (2007), inspired by the seminal work
of Zel’dovich (1970), have provided a222Many other cosmic web classifiers
exist, such as DIVA (Lavaux & Wandelt, 2010), V-web (Hoffman et al., 2012),
and LICH (Leclercq et al., 2017). 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.333Forero-Romero et al. (2009) subsequently refined
Hahn’s definition, by counting eigenvalues above a tuneable threshold (rather
than zero) related to the gravitational collapse timescale. The partitioning
of the cosmic web into the four structure types depends sensitively on the
choice of this threshold, and it should therefore always be mentioned when
quantitative structure-type properties are stated, such as volume- and mass-
filling fractions.

###  1.2 Baryons in the ICM and WHIM

Of the four structure types, galaxy clusters are most easily studied: in the
X-ray, optical and radio bands, e.g. through thermal Bremsstrahlung,
gravitational lensing and synchrotron radiation. As the most massive
gravitationally bound structures in the Universe thus far, they weigh up to
${\sim}10^{15}\ M_{\odot}$, contain up to ${\sim}10^{3}$ galaxies, and are
pervaded by a dilute (${\sim}10^{3}\ \mathrm{m}^{-3}$), hot (${\sim}10^{8}\
\mathrm{K}$) and magnetised (${\sim}1\ \mathrm{\mu 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 (${\sim}1-10\ \mathrm{m}^{-3}$), cooler (${\sim}10^{5}-10^{7}\
\mathrm{K}$) and less magnetised (${\sim}10^{-3}-10^{-1}\ \mathrm{\mu 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 $\sim 90\%$ of
the Universe’s baryons (Cen & Ostriker, 1999; Eckert et al., 2015; de Graaff
et al., 2019; Tanimura et al., 2019).

###  1.3 Synchrotron radiation from merger and accretion shocks

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 $\mathcal{M}_{\mathrm{u}}\sim 10^{0}-10^{2}$ (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
$\mathcal{M}_{\mathrm{u}}\sim 1-5$) (Ryu et al., 2003), shocks in filaments
remain fainter due to the aforementioned adverse density, temperature and
magnetic field strength conditions.444Precisely how much fainter shocks in
filaments are compared to those in clusters, is unknown, because the typical
filament IGM magnetic field strength $B_{\mathrm{IGM}}$ and the electron
acceleration efficiency $\xi_{e}$ of shocks — and especially weak ones — is
highly uncertain. This is why observing filaments through synchrotron
radiation constitutes an ambitious, futuristic frontier.

###  1.4 Scientific prospects of detecting synchrotron radiation from the
filament IGM

####  1.4.1 Physics of astrophysical shocks

Figure 1: The dependence of three central (merger or accretion) shock
quantities on upstream Mach number $\mathcal{M}_{\mathrm{u}}$ and adiabatic
index $\gamma$, as derived from the Rankine–Hugoniot jump conditions in ideal
gases. Radiation from $\mathcal{M}_{\mathrm{u}}<10$ shocks is predominantly
generated by electrons that stem from the ICM, whilst radiation from
$\mathcal{M}_{\mathrm{u}}>10$ shocks is predominantly generated by electrons
that stem from the WHIM. Top: the spectral index $\alpha$ of angularly
unresolved shocks, assuming standard DSA. Middle: the downstream-over-upstream
plasma density ratio. Bottom: the downstream-over-upstream plasma temperature
ratio.

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 Figure 1 shows that DSA predicts an (angularly unresolved)
synchrotron spectral index $\alpha=-1$ for virtually all filament shocks (with
a slight dependency on the adiabatic index $\gamma$). Significant deviations
from $\alpha=-1$ would falsify standard DSA.555Caprioli & Haggerty (2019)
describe recent advances in DSA theory, which suggest steeper electron energy
and synchrotron spectra for the filaments’ strong shocks (i.e. $\alpha<-1$ for
$\gamma=\nicefrac{{5}}{{3}}$ and $\mathcal{M}_{\mathrm{u}}>10$). 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.

####  1.4.2 Missing baryons

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 $\Lambda$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 ${\sim}34\%$
is not identified conclusively (de Graaff et al., 2019).

Figure 2: Two simulated specific intensity functions at
$\nu_{\mathrm{obs}}=150\ \mathrm{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
$B_{\mathrm{IGM}}=1\ \mathrm{nG}$ in the Early Universe. (In contrast, the
simulation underlying this article’s predictions starts at
$B_{\mathrm{IGM}}=0.1\ \mathrm{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.

####  1.4.3 Magnetogenesis

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^{-22}\
\mathrm{T}$ vs. $10^{-13}\ \mathrm{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.
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_{\mathrm{IGM}}<0.03-1.98\ \mathrm{\mu G}$, while
Brown et al. (2017) used Parkes 64m Telescope S-PASS data to find a density-
weighted upper limit of $B_{\mathrm{IGM}}<0.13\ \mathrm{\mu 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\
\mathrm{nG}<B_{\mathrm{IGM}}<60\ \mathrm{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_{\mathrm{IGM}}<4\ \mathrm{nG}$ in
filaments. Vernstrom et al. (2019) have performed a similar analysis with NVSS
data, finding $B_{\mathrm{IGM}}<40\ \mathrm{nG}$ in filaments.

###  1.5 Article structure

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
$\mathrm{Map}\left(\mathbb{S}^{2},\mathbb{R}_{\geq 0}\right)$ of specific
intensity functions on the 2-sphere $\mathbb{S}^{2}$. 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 Section 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 Section 3, we develop a simple geometric
model of the cosmic web that yields SCW redshift predictions. In Section 4, we
analyse predictions of the main method, and compare these with predictions of
the geometric model. Finally, in Section 5, we discuss caveats of the current
work, and give recommendations for future extensions, before we present
conclusions in Section 6.
We adopt a666Rather fortunately, the BORG SDSS and Enzo data products combined
in this work have been made assuming almost identical cosmological parameters;
the set reported is of the BORG SDSS. Authors of future work who strive to
achieve _full_ self-consistency should use the same cosmological model for the
LSS reconstructions as for the MHD simulations. concordance inflationary
$\Lambda$CDM cosmological model
$\mathfrak{M}=(\Omega_{\Lambda,0}=0.728,\Omega_{\mathrm{DM},0}=0.227,$
$\Omega_{\mathrm{BM},0}=0.045,h=0.702,\sigma_{8}=0.807,n_{s}=0.961)$, so that
$\Omega_{\mathrm{M},0}\coloneqq\Omega_{\mathrm{DM},0}+\Omega_{\mathrm{BM},0}=0.272$
and the Hubble constant $H_{0}\coloneqq h\cdot 100\ \mathrm{km\ s^{-1}\
Mpc^{-1}}$ (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. 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. 2.

Then, using a cosmological MHD simulation, we establish the connection between
the total matter density $\rho$ and the merger and accretion shock synchrotron
monochromatic emission coefficient (MEC) $j_{\nu}$, in the form of the
conditional probability distribution $P\left(j_{\nu}\left(\mathbf{r}\right)\
|\ \rho\left(\mathbf{r}\right)\right)$.

  3. 3.

Next, to create a MEC field, one could independently realise the random
variables (RVs) $j_{\nu}\left(\mathbf{r}\right)\ |\
\rho\left(\mathbf{r}\right)$ for all $\mathbf{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 generate 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. Note that 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_{\nu}\left(\mathbf{r}\right)\ |\
\rho=\rho\left(\mathbf{r}\right)$ for all $\mathbf{r}$ in the volume.

  4. 4.

After having obtained a merger and accretion shock synchrotron MEC field with
spatial correlation, we ray trace 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. 5.

By repeating the above steps over and over, picking different density field
and GRF realisations every time, we generate a probability distribution on
$\mathrm{Map}\left(\mathbb{S}^{2},\mathbb{R}_{\geq 0}\right)$ for the MASSCW
signal.

  6. 6.

With a minor modification, this procedure naturally leads to another
probability distribution on $\mathrm{Map}\left(\mathbb{S}^{2},\mathbb{R}_{\geq
0}\right)$, 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., 2019) 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.777On super-Mpc scales, baryonic effects
such as gas and radiation pressure play a subdominant role in structure
formation, and can therefore be safely ignored. Thus, realisations from
today’s highly non-Gaussian total density field can be generated by forward
modelling the effect of gravity on a collisionless888The assumption of a
collisionless fluid is apt for DM, but only approximate for BM. fluid, which
is initialised as a Gaussian random field at the time of the CMB.999Because
simulating gravity with an $N$-body simulation from the CMB to the present day
is computationally expensive, and the method requires this process to be
repeated thousands of times, all authors resort to methods that approximate
Newtonian gravitational evolution. In order of increasing accuracy (and
numerical cost), they invoke either first-order Lagrangian perturbation theory
(i.e. the Zel’dovich approximation (ZA) (Zel’dovich, 1970)), second-order
Lagrangian perturbation theory (2LPT) or particle mesh (PM) models, which
approach the $N$-body solution. 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 $\phi$, declinations $\theta$,
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\ \mathrm{deg}^{2}$ of Northern Sky, and three stripes of the
Southern Galactic Cap (SGC), which together cover $386\ \mathrm{deg}^{2}$. 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
${\sim}3\ \mathrm{Mpc}\ \mathrm{h}^{-1}$ and larger. Each realisation of the
posterior is a cube containing $256^{3}$ voxels, that represents a region
$\mathcal{R}\subset\mathbb{R}^{3}$ with comoving volume $(750\ \mathrm{Mpc}\
\mathrm{h}^{-1})^{3}$. Naturally, $\mathcal{R}$ is chosen to correspond to the
half of the Northern Sky that contains the NGC (up to a cosmological redshift
$z\sim 0.2$).

###  2.3 Relating total matter density to synchrotron monochromatic emission
coefficient

Next, to convert total matter density fields $\rho$ to proper synchrotron
monochromatic emission coefficient (MEC) fields $j_{\nu}$, we establish the
relation between $\rho$ and $j_{\nu}$ (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.101010Araya-Melo et al. (2012) also took this
approach and used the MareNostrum simulation to establish that, under these
assumptions, filaments at a redshift of $0.15$ should produce a flux density
of $10^{-1}\ \mathrm{\mu Jy}$ at $150\ \mathrm{MHz}$.
We assume that the power density $P_{\nu}$ of a single shock obeys Equation 32
of Hoeft & Brüggen (2007). Let the shock surface area be $A$, and let the
effective width of the downstream region be $\langle y\rangle$ (intuitively,
this is the thickness of the shock in the direction perpendicular to the
surface). The effective shock volume is $V=A\langle y\rangle$. If shocks
radiate isotropically, the MEC $j_{\nu}$ is direction-independent. However,
like $\rho$, $j_{\nu}$ 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_{\nu,i}\left(z,\nu\right)$ is the power density of shock $i$ at
cosmological redshift $z$ and emission frequency $\nu$, and $A_{i}$ and
$\langle y\rangle_{i}$ are its surface area and effective width.111111Because
the (cosmology-dependent) function $z\left(t\right)$ is strictly decreasing
and thus invertible, we can use $z$ as a time coordinate. Moreover, let
$\mathcal{R}_{i}\left(z\right)\subset\mathbb{R}^{3}$ be the region of space
occupied by this shock at redshift $z$; the Lebesgue measure of
$\mathcal{R}_{i}$ is
$V_{i}=\int_{\mathcal{R}_{i}}\mathrm{d}\mathbf{r}=A_{i}\langle y\rangle_{i}$.
Under the approximation that the shock emission is homogeneous within
$\mathcal{R}_{i}$, the average _total_ MEC within an arbitrary region
$\bar{\mathcal{R}}\subset\mathbb{R}^{3}$ with Lebesgue measure
$\bar{V}=\int_{\bar{\mathcal{R}}}\mathrm{d}\mathbf{r}$ is

$\displaystyle j_{\nu}\left(z,\nu\right)=\frac{1}{4\pi\
\mathrm{sr}}\frac{1}{\bar{V}}\sum_{i}P_{\nu,i}\left(z,\nu\right)\frac{\int_{\mathcal{R}_{i}\left(z\right)\
\cap\ \bar{\mathcal{R}}}\mathrm{d}\mathbf{r}}{V_{i}}.$ (1)

The MEC of a _single_ shock (thus dropping indices) located at $\mathbf{r}$ at
redshift $z$, and averaged over its own effective volume, is

$\displaystyle j_{\nu}\left(\mathbf{r},z,\nu\right)=\frac{1}{4\pi\
\mathrm{sr}}\frac{P_{\nu}}{V}$ $\displaystyle=\frac{445\
\mathrm{Jy}}{\mathrm{deg^{2}\ Mpc}}\ \left(\frac{\langle
y\rangle}{\mathrm{Mpc}}\right)^{-1}\frac{n_{e,\mathrm{d}}\left(\mathbf{r},z\right)}{10\
\mathrm{m}^{-3}}\frac{\xi_{e}}{10^{-2}}\left(\frac{\nu}{150\
\mathrm{MHz}}\right)^{\alpha\left(\mathbf{r},z\right)}$
$\displaystyle\left(\frac{T_{\mathrm{d}}\left(\mathbf{r},z\right)}{10^{8}\
\mathrm{K}}\right)^{\frac{3}{2}}\frac{\left(\frac{B_{\mathrm{d}}}{\mathrm{\mu
G}}\right)^{1-\alpha\left(\mathbf{r},z\right)}\left(\mathbf{r},z\right)}{\left(\frac{B_{\mathrm{d}}}{\mathrm{\mu
G}}\right)^{2}\left(\mathbf{r},z\right)+\left(\frac{B_{\mathrm{CMB}}}{\mathrm{\mu
G}}\right)^{2}\left(z\right)}\left(\frac{150}{1400}\right)^{\alpha\left(\mathbf{r},z\right)}\Psi\left(\mathbf{r},z\right).$
(2)

Here, $n_{e,\mathrm{d}}$ is the downstream electron number density,
$T_{\mathrm{d}}$ is the downstream plasma temperature, $\alpha$ is the
integrated spectral index of the associated synchrotron emission,
$B_{\mathrm{d}}$ is the downstream plasma magnetic field strength,
$B_{\mathrm{CMB}}$ is the CMB magnetic field strength, and $\Psi$ is a
dimensionless quantity with a strong dependence on the upstream Mach number
$\mathcal{M}_{\mathrm{u}}$ and a weak dependence on $T_{\mathrm{d}}$ (as can
be glanced from Figure 4 of Hoeft & Brüggen (2007)). $\Psi$ approaches unity
for high Mach numbers, and so for the WHIM, where the upstream Mach numbers
are expected to be high, $\Psi\left(\mathbf{r},z\right)\approx 1$. Like
$\Psi$, $\alpha$ does not depend on the spacetime coordinate
$\left(\mathbf{r},z\right)$ _directly,_ but rather via
$\mathcal{M}_{\mathrm{u}}$ and the adiabatic index $\gamma$. Concretely, for
$\mathcal{M}_{\mathrm{u}}>1$,

$\displaystyle\alpha=\alpha\left(\mathcal{M}_{\mathrm{u}}\left(\mathbf{r},z\right),\gamma\left(\mathbf{r},z\right)\right)=\frac{\tfrac{1}{4}\left(1-3\gamma\right)\mathcal{M}_{\mathrm{u}}^{2}-1}{\mathcal{M}_{\mathrm{u}}^{2}-1}.$
(3)

Note that $\alpha\to-\infty$ as $\mathcal{M}_{\mathrm{u}}\to 1+$, and that
$\alpha\to\nicefrac{{1}}{{4}}\left(1-3\gamma\right)$ as
$\mathcal{M}_{\mathrm{u}}\to\infty$. The top panel of Figure 1 illustrates
this behaviour. The electron acceleration efficiency $\xi_{e}$ quantifies the
fraction of the shock’s thermal energy that is used to accelerate suprathermal
electrons. $\xi_{e}$ is unknown for WHIM shocks, but a comparison to supernova
remnant (SNR) shocks suggests $\xi_{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
turn to MHD simulations. In particular, we use 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
$\left(100\ \mathrm{Mpc}\right)^{3}$ with $2400^{3}$ voxels, yielding a
(comoving) resolution of $41\tfrac{2}{3}\ \mathrm{kpc}$ per voxel edge. The
simulations recreate the evolution of the baryonic and dark matter density
functions $\rho_{\mathrm{BM}}$ and $\rho_{\mathrm{DM}}$ — as well as
(thermo)dynamic quantities121212As long as one considers a single shock, space
(at a fixed time) can be classified into an upstream region and a downstream
region, with associated temperatures ($T_{\mathrm{u}}$ and $T_{\mathrm{d}}$,
respectively) and magnetic field strengths ($B_{\mathrm{u}}$ and
$B_{\mathrm{d}}$, respectively). This distinction is less meaningful once one
considers multiple shocks at the same time: a given location could then be
upstream for some shocks, and downstream for others. Therefore, cosmological
simulations do not evolve an upstream or downstream temperature field, but
just a _general_ temperature field $T$. Analogously, cosmological simulations
maintain only _general_ magnetic field component fields $B_{x}$, $B_{y}$ and
$B_{z}$ (and thus a magnetic field strength field
$B\coloneqq\sqrt{B_{x}^{2}+B_{y}^{2}+B_{z}^{2}}$), without upstream and
downstream distinction. 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).131313From the middle panel of Figure 1,
it is clear that shocks also induce a jump in gas density. However, this jump
is modest, saturating at e.g. a factor $4$ for $\gamma=\nicefrac{{5}}{{3}}$.
By contrast, the jump in temperature, seen in the bottom panel of the same
figure, can be several orders of magnitude — in the WHIM, at least. For
shocked voxels, we simply take $\langle y\rangle$ to be the voxel edge length,
and calculate $\alpha$ and $\Psi$ after establishing
$\mathcal{M}_{\mathrm{u}}$. With the exception of $\xi_{e}$, which we will
assume to be a constant throughout, the simulation thus enables us to compute
all factors on the RHS of Equation 2.
Taking their product for the $z=0.025$ snapshot, we obtain $j_{\nu}$, and
compare it to $\rho\coloneqq\rho_{\mathrm{BM}}+\rho_{\mathrm{DM}}$ to study
their relationship.141414By repeating our analysis for snapshots of other
redshifts, one could study the time evolution of this relationship. However,
as the Hoeft & Brüggen (2007) model is likely only a rough description of the
actual synchrotron emission mechanism in filaments, we currently consider such
level of detail superfluous.

####  2.3.4 Conditional probability distribution

In order to convert the BORG SDSS total matter density fields into synchrotron
MEC fields, we need to predict $j_{\nu}$ from $\rho$ — ideally including
uncertainty. Because the BORG SDSS has a comoving resolution of $4.17\
\mathrm{Mpc}$ per voxel edge compared to Enzo’s $4.17\cdot 10^{-2}\
\mathrm{Mpc}$ per voxel edge, we blur both $j_{\nu}$ and $\rho$ with a
Gaussian kernel whose standard deviation is half of the comoving resolution
ratio. We then treat each $(j_{\nu},\rho)$ pair of the blurred fields as a
draw from the joint distribution $P\left(j_{\nu},\rho\right)$. By binning the
data, we perform a simplistic form of kernel density estimation (KDE). Next,
we calculate the conditional $P\left(j_{\nu}\ |\ \rho\right)$ from the joint
by dividing it by the marginal
$P\left(\rho\right)\coloneqq\int_{0}^{\infty}P\left(j_{\nu},\rho\right)\mathrm{d}j_{\nu}$.
The result is the probability distribution shown in Figure 3.

Figure 3: Conditional probability distributions $P\left(j_{\nu}\ |\
\rho\right)$ of synchrotron monochromatic emission coefficient $j_{\nu}$ (at
rest-frame (emission) frequency $\nu=150\ \mathrm{MHz}$) given total matter
density $\rho$ 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_{\nu}\propto\rho^{\nicefrac{{23}}{{6}}}$
(single shock) scaling relation for various proportionality constants. $\rho$
is shown relative to the current-day mean total matter density
$\bar{\rho}_{0}=\Omega_{\mathrm{M,0}}\ \rho_{c,0}$, while the weakly
constrained electron acceleration efficiency $\xi_{e}$ is divided out from
$j_{\nu}$. 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_{\nu}$.

We overplot the single shock scaling relation expected for the WHIM regime:
$j_{\nu}\propto\rho^{\nicefrac{{23}}{{6}}}$. See Appendix B for a derivation.

Figure 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 $\mathcal{R}$, of $750\ \mathrm{Mpc}\ h^{-1}\times 750\
\mathrm{Mpc}\ h^{-1}$ for the density and MEC fields, and a subregion of $150\
\mathrm{Mpc}\ h^{-1}\times 150\ \mathrm{Mpc}\ h^{-1}$ for the Gaussian and
percentile random fields (to more clearly illustrate the bijective mapping
between the two). Left: BORG SDSS density $\rho$ sample, in multiples of the
current-day mean total matter density
$\bar{\rho}_{\mathrm{0}}=\Omega_{\mathrm{M,0}}\ \rho_{\mathrm{c,0}}$. Middle
left: Gaussian random field $\mathcal{Z}$ sample using a squared exponential
(SE) kernel with lengthscale $l_{\mathrm{SE}}=2\ \mathrm{Mpc}$. Middle right:
the corresponding percentile random field $\mathcal{P}$ sample. Right: the
corresponding MEC $j_{\nu}$ sample at an emission frequency of $\nu=150\
\mathrm{MHz}$, with the weakly-constrained electron acceleration efficiency
$\xi_{e}$ divided out.

###  2.4 Generating monochromatic emission coefficient fields

In principle, we could — for each BORG SDSS sample
$\rho\left(\mathbf{r}\right)$ — convert to $j_{\nu}\left(\mathbf{r}\right)$ by
drawing from $P\left(j_{\nu}\ |\ \rho=\rho\left(\mathbf{r}\right)\right)$
independently for every $\mathbf{r}\in\mathcal{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 use realisations of a Gaussian random field (GRF)
$\mathcal{Z}$ over $\mathcal{R}$. 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 choose an isotropic and stationary kernel,
as suggested by the Cosmological Principle. From the variety of remaining
choices commonly used in machine learning (ML), we pick151515The choice of
kernel should ideally reflect the morphology of merger and accretion shocks.
However, given the current spatial resolution of the BORG SDSS density fields
and our requirement that the kernel be isotropic (disregarding the relation
between shock morphology and local LSS orientation), the approach is already
so approximate that no particular kernel is clearly preferred over the others.
the squared-exponential (SE) kernel
$K_{\mathrm{SE}}:\mathcal{R}\times\mathcal{R}\to\mathbb{R}_{\geq 0}$
(Rasmussen & Williams, 2006):

$\displaystyle
K_{\mathrm{SE}}\left(\mathbf{r},\mathbf{r}^{\prime}\right)=\sigma_{\mathrm{SE}}^{2}\exp{\left(-\frac{||\mathbf{r}-\mathbf{r}^{\prime}||_{2}^{2}}{2l_{\mathrm{SE}}^{2}}\right)},\
\textrm{so that}$
$\displaystyle\mathrm{Cov}\left(\mathcal{Z}\left(\mathbf{r}\right),\mathcal{Z}\left(\mathbf{r}^{\prime}\right)\right)=K_{\mathrm{SE}}\left(\mathbf{r},\mathbf{r}^{\prime}\right)\
\textrm{with}\ \mathbf{r},\mathbf{r}^{\prime}\in\mathcal{R}.$ (4)

The variance $\sigma_{\mathrm{SE}}^{2}$ allows us to endow the GRF with
dimensionality and scale it at will. As our GRF will only serve to draw
spatially correlated MEC samples, we set $\sigma_{\mathrm{SE}}=1$.
The lengthscale $l_{\mathrm{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 ${\sim}2\ \mathrm{Mpc}$. Although the majority of observed shocks
seem smaller (suggesting $l_{\mathrm{SE}}<2\ \mathrm{Mpc}$), the presence of
noise in radio imagery obfuscates the full extent of shock fronts, biasing
measurements towards small LLSs. As a compromise, we choose
$l_{\mathrm{SE}}=2\ \mathrm{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 generate our zero mean, unit covariance GRF realisations on the
BORG SDSS voxel grid using Fourier analysis. This necessitates calculating the
power spectrum of the kernel, which in the case of the 3D SE kernel is

$\displaystyle P_{\mathrm{SE}}\left(\mathbf{k}\right)=\left(2\pi
l_{\mathrm{SE}}^{2}\right)^{\frac{3}{2}}\exp{\left(-2\pi||\mathbf{k}||_{2}^{2}l_{\mathrm{SE}}^{2}\right)}.$
(5)

Finally, Fourier-generated GRFs on a finite (numerical) grid require an
appropriate normalisation; we use Parseval’s theorem to find the correct
factor.
This procedure generates a ‘standard normal’ GRF $\mathcal{Z}$, in the sense
that each location’s RV is standard normal:
$\mathcal{Z}\left(\mathbf{r}\right)\sim\mathcal{N}\left(0,1\right)$ for all
$\mathbf{r}\in\mathcal{R}$. Next, we apply the CDF $\Phi$ of the standard
normal to the GRF in point-wise fashion (i.e. voxel-wise, in practice),
thereby creating a ‘percentile random field’ $\mathcal{P}$:

$\displaystyle\mathcal{P}\left(\mathbf{r}\right)\coloneqq\Phi\left(\mathcal{Z}\left(\mathbf{r}\right)\right),\
\textrm{where}\
\Phi\left(x\right)\coloneqq\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{x}\exp{\left(\frac{-y^{2}}{2}\right)}\
\mathrm{d}y.$ (6)

This field has values strictly between $0$ and $1$, and inherits the spatial
correlations present in $\mathcal{Z}$.
Let $\Xi_{\rho}\left(x\right):\mathbb{R}_{\geq 0}\to[0,1]$ be the CDF of the
RV $j_{\nu}\ |\ \rho$, and let
$\Xi_{\rho}^{-1}\left(y\right):[0,1]\to\mathbb{R}_{\geq 0}$ be its inverse.
Then the final MEC field is given by

$\displaystyle
j_{\nu}\left(\mathbf{r}\right)\coloneqq\Xi_{\rho\left(\mathbf{r}\right)}^{-1}\left(\mathcal{P}\left(\mathbf{r}\right)\right),\
\textrm{for all}\ \mathbf{r}\in\mathcal{R}.$ (7)

Because $\rho$ and $\mathcal{P}$ are random fields on $\mathcal{R}$, so is
$j_{\nu}$. A graphical summary of the procedure described in this subsection,
and defined in Equation 7, is shown in Figure 4.

###  2.5 Generating specific intensity functions

To find the specific intensity of the SCW _on the sky_ ($I_{\nu}$), we
simulate 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 neglect 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 estimate $I_{\nu}$ in a sky patch by sampling many directions $\hat{r}_{i}$
within it, and simulating the passage of light through our reconstructed SCW
in comoving space for each such direction (or ‘ray’) $\hat{r}_{i}$.

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

Next, we establish the relation between the observed specific intensity
$I_{\nu}$ of a ray that has travelled through the SCW to Earth from direction
$\hat{r}$, and the MEC $j_{\nu}$ along its path.
Be $j_{\nu}:\mathbb{R}^{3}\times\mathbb{R}_{\geq
0}\times\mathbb{R}_{>0}\to\mathbb{R}_{\geq 0}$ the function that assigns to
each 3-tuple $\left(\mathbf{r},z,\nu\right)$ containing a comoving locus
$\mathbf{r}$, a cosmological redshift $z$ (which represents time) and an EM
wave frequency $\nu$, the _proper_ (rather than comoving) MEC
$j_{\nu}\left(\mathbf{r},z,\nu\right)$ (with SI units $\mathrm{W\ m^{-3}\
Hz^{-1}\ sr^{-1}}$). 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_{\nu}$
for direction $\hat{r}$ and observing frequency $\nu_{\mathrm{obs}}$ is

$\displaystyle
I_{\nu}\left(\hat{r},\nu_{\mathrm{obs}}\right)=\frac{c}{H_{0}}\int_{0}^{\infty}\frac{j_{\nu}\left(r\left(z\right)\hat{r},z,\nu_{\mathrm{obs}}\left(1+z\right)\right)}{\left(1+z\right)^{4}E\left(z\right)}\mathrm{d}z.$
(8)

We provide a brief derivation, alongside explicit expressions for
$r\left(z\right)$ and $E\left(z\right)$, in Appendix D.
For each specific intensity function that we wish to generate, we now draw
$M=10^{6}$ ray directions uniformly from the $\pi\ \mathrm{sr}$ lune that
covers the SDSS-constrained half of the Northern Sky. For each ray, we
calculate the corresponding specific intensity by combining Equations 2 and 8.
Note that this requires evaluating the MEC field at a range of emission
frequencies. Numerically, we realise this by building both the Enzo-derived
conditional probability distribution $P\left(j_{\nu}\ |\ \rho\right)$ and the
MEC field sample (using the same density field sample and GRF sample) at two
emission frequencies $\nu$: $\nu=\nu_{\mathrm{obs}}$ and
$\nu=\nu_{\mathrm{obs}}\left(1+z_{\mathrm{max}}\right)$. Next, for each voxel,
we calculate a spectral index from the two MEC fields by assuming a power-law
spectrum, and use it to find the MEC at the emission frequency
$\nu=\nu_{\mathrm{obs}}\left(1+z\right)$ needed given the voxel’s cosmological
redshift $z$.

We repeat this process many times by selecting $1000$ density samples from the
BORG SDSS posterior (we discard the first $2110$ samples due to burn-in, and
thin the Monte Carlo Markov chain by a factor $10$) and generating $1$
independent and identically distributed (IID) GRF sample for each. (There is
no compelling reason to thin the chain and we could use multiple GRF samples
per density sample; the current approach is merely to limit data storage,
compute time and energy usage.) Thus, we numerically realise $1000\cdot
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 weighting 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 $\bar{z}$ at observing frequency
$\nu_{\mathrm{obs}}$ of a ray $\hat{r}$ passing through MEC field $j_{\nu}$ to
be

$\displaystyle\bar{z}\left(\hat{r},\nu_{\mathrm{obs}}\right)\coloneqq
I_{\nu}^{-1}\left(\hat{r},\nu_{\mathrm{obs}}\right)\ \int_{0}^{\infty}z\
\frac{\mathrm{d}I_{\nu}}{\mathrm{d}z}\left(\hat{r},\nu_{\mathrm{obs}},z\right)\
\mathrm{d}z.$ (9)

By iterating over total matter density and GRF samples, one generates a
probability distribution over the function
$\bar{z}\left(\hat{r},\nu_{\mathrm{obs}}\right)$, in exact analogy to the
generation of the distribution over specific intensity functions
$I_{\nu}\left(\hat{r},\nu_{\mathrm{obs}}\right)$.
One might wonder, _If I were to pick a random direction on the sky, what would
be the specific-intensity-weighted mean redshift of its MASSCW signal?_ The
set $\\{\bar{z}\left(\hat{r}_{i},\nu_{\mathrm{obs}}\right)\ |\
i\in\\{1,2,...,M\\}\\}$, with rays $\hat{r}_{i}$ drawn from a uniform
distribution over the sky, can be viewed as a random sample from a specific-
intensity-weighted mean redshift _random variable_ $\bar{Z}$. An RV is fully
characterised by its CDF. Each MASSCW prior realisation generates another
empirical CDF (ECDF) for $\bar{Z}$.161616The probability distribution over
ECDFs of $\bar{Z}$ thus obtained is less informative than the probability
distribution over functions $\bar{z}$, because it disposes of the directional
correlations in specific-intensity-weighted mean redshift
$\bar{z}\left(\hat{r},\nu_{\mathrm{obs}}\right)$.

####  2.6.2 Flux-density-weighted mean redshift

One might also wonder, _What is the mean redshift of the MASSCW signal — 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 $\nu_{\mathrm{obs}}$
to be

$\displaystyle\bar{\bar{z}}\coloneqq\left(\int_{\mathbb{S}^{2}}I_{\nu}\left(\hat{r},\nu_{\mathrm{obs}}\right)\mathrm{d}\Omega\right)^{-1}\int_{\mathbb{S}^{2}}\int_{0}^{\infty}z\
\frac{\mathrm{d}I_{\nu}}{\mathrm{d}z}\left(\hat{r},\nu_{\mathrm{obs}},z\right)\
\mathrm{d}z\ \mathrm{d}\Omega,$ (10)

where $\bar{\bar{z}}$ is shorthand for
$\bar{\bar{z}}\left(\nu_{\mathrm{obs}}\right)$.

##  3 Redshift predictions from geometric cosmic web model

###  3.1 Overcoming the redshift limitation

$l$$l$$w_{\mathrm{f}}$$w_{\mathrm{f}}$$w_{\mathrm{f}}$$w_{\mathrm{f}}$ Figure
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 $w_{\mathrm{f}}$ and length $l_{\mathrm{f}}\coloneqq
l-w_{\mathrm{f}}$ 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.

Although we can calculate the specific intensity function
$I_{\nu}\left(\hat{r},\nu_{\mathrm{obs}}\right)$, the specific-intensity-
weighted mean redshift function
$\bar{z}\left(\hat{r},\nu_{\mathrm{obs}}\right)$, the single-direction
redshift RV $\bar{Z}$ ECDF and the flux-density-weighted mean redshift
$\bar{\bar{z}}\left(\nu_{\mathrm{obs}}\right)$ for each realisation of the
MASSCW prior, all four quantities suffer from the fact that the BORG SDSS
reconstructions stop at a redshift $z_{\mathrm{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 $\bar{Z}$, and $\bar{\bar{z}}$, for an arbitrary value of
$z_{\mathrm{max}}$. The model reproduces the results calculated from the
MASSCW prior, if it is equally limited to $z_{\mathrm{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 $\bar{Z}$ and $\bar{\bar{z}}$, for the
case that describes the actual Universe: $z_{\mathrm{max}}\to\infty$.

###  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_{\mathrm{f}}$ and length $l_{\mathrm{f}}\coloneqq l-w_{\mathrm{f}}$ (and
thus volume $w_{\mathrm{f}}^{2}l_{\mathrm{f}}$).171717Often, filaments are
modelled as cylinders. In such cases, the length of the path of a sightline
through a filament depends on the exact point of incidence on the unit cell
boundary. However, in our simplistic geometric model, this effect does not
arise. Galaxy clusters reside where filaments meet, and — in this simplistic
model — are cubes with comoving edge length $w_{\mathrm{f}}$ (and thus volume
$w_{\mathrm{f}}^{3}$). By far the largest component of a unit cell is its
central void, a cubical region of comoving edge length $l_{\mathrm{f}}$ (and
thus volume $l_{\mathrm{f}}^{3}$). Finally, unit cells contain sheets, which
fill the regions bounded by filaments and voids. They are of typical comoving
length $l_{\mathrm{f}}$ and thickness $w_{\mathrm{f}}$ (and thus volume
$l_{\mathrm{f}}^{2}w_{\mathrm{f}}$). 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.
If one would choose a point-of-incidence on the face shown in Figure 5 in
uniform fashion, the sightline would hit a sheet and then a void (dark-blue-
shaded area) with probability

$\displaystyle
p_{\mathrm{s-v}}\left(\frac{w_{\mathrm{f}}}{l}\right)=\frac{\left(l-w_{\mathrm{f}}\right)^{2}}{l^{2}}=\left(1-\frac{w_{\mathrm{f}}}{l}\right)^{2}.$
(11)

Note that $p_{\mathrm{s-v}}\to 1$ when $\frac{w_{\mathrm{f}}}{l}\to 0$, and
$p_{\mathrm{s-v}}\to 0$ when $\frac{w_{\mathrm{f}}}{l}\to 1$, as required.
Likewise, the probability that the sightline hits a filament only (light-blue-
shaded area) is

$\displaystyle
p_{\mathrm{f}}\left(\frac{w_{\mathrm{f}}}{l}\right)=\frac{2w_{\mathrm{f}}\left(l-w_{\mathrm{f}}\right)}{l^{2}}=2\frac{w_{\mathrm{f}}}{l}\left(1-\frac{w_{\mathrm{f}}}{l}\right).$
(12)

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

$\displaystyle
p_{\mathrm{c-f}}\left(\frac{w_{\mathrm{f}}}{l}\right)=\left(\frac{w_{\mathrm{f}}}{l}\right)^{2}.$
(13)

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

Figure 6:  When a sightline enters a cubic unit cell of the cosmic web
perpendicularly to a face (as in Figure 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 Equations 11, 12 and 13).

Pick a sightline, and successively label each unit cell boundary this
sightline crosses with a natural number $n\in\mathbb{N}_{\geq 1}$. If the
first boundary crossing occurs at a comoving distance $d_{1}$, then the $n$-th
crossing happens at comoving distance $d_{n}=d_{1}+\left(n-1\right)l$. Let
$\mathfrak{M}$ be the cosmological model of preference, and let
$z_{\mathfrak{M}}\left(d\right)$ denote the function that converts comoving
distance to cosmological redshift under this model. The unit cell boundary
crossings occur at redshifts $z_{n}\coloneqq
z_{\mathfrak{M}}\left(d_{n}\right)$ for all $n\in\mathbb{N}_{\geq 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_{\mathrm{max}}$ only, we introduce a parameter $N$. Intuitively,
$N\in\mathbb{N}_{\geq 1}$ is the label of the last boundary crossing within
the LSS considered. More formally, $N\coloneqq\max\ \\{n\in\mathbb{N}_{\geq
1}\ |\ z_{n}<z_{\mathrm{max}}\\}$.
Let $X_{n}\left(\nu_{\mathrm{obs}}\right)$ be the RV denoting the contribution
to the MASSCW specific intensity at observing frequency $\nu_{\mathrm{obs}}$
picked up by the sightline during the $n$-th full unit cell crossing (i.e.
whilst travelling through the $n$-th _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 $n$-th 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}\left(\nu_{\mathrm{obs}}\right)\coloneqq
X\left(\nu_{\mathrm{obs}},z_{n}\right)$. $X\left(\nu_{\mathrm{obs}},z\right)$
is the RV with support $\left(0,\Delta
I_{\nu,\mathrm{f}}\left(\nu_{\mathrm{obs}},z\right),\Delta
I_{\nu,\mathrm{c-f}}\left(\nu_{\mathrm{obs}},z\right)\right)$, and
corresponding probabilities
$\left(p_{\mathrm{s-v}},p_{\mathrm{f}},p_{\mathrm{c-f}}\right)$. Note that we
assume that sheets and voids have a vanishing contribution to the MASSCW
signal. Here,

$\displaystyle\Delta I_{\nu,\mathrm{f}}\left(\nu_{\mathrm{obs}},z\right)$
$\displaystyle\coloneqq
w_{\mathrm{f}}j_{\nu,\mathrm{f}}\left(\nu_{\mathrm{obs}}\left(1+z\right),z\right)\left(1+z\right)^{-4}$
(14) $\displaystyle\Delta
I_{\nu,\mathrm{c-f}}\left(\nu_{\mathrm{obs}},z\right)$
$\displaystyle\coloneqq(w_{\mathrm{f}}j_{\nu,\mathrm{c}}\left(\nu_{\mathrm{obs}}\left(1+z\right),z\right)$
$\displaystyle+l_{\mathrm{f}}j_{\nu,\mathrm{f}}\left(\nu_{\mathrm{obs}}\left(1+z\right),z\right))\left(1+z\right)^{-4}.$
(15)

The factor $\left(1+z\right)^{-4}$ follows from Equation 52: 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 $\Delta 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_{\nu,\mathrm{f}}\left(\nu,z\right)$ and
$j_{\nu,\mathrm{c}}\left(\nu,z\right)$ represent the typical proper filament
and cluster SCW MECs at emission frequency $\nu$ and cosmological redshift
$z$. We use the dimensionless parameter $\mathcal{C}$ to denote the typical
cluster-to-filament-SCW MEC ratio at some reference emission frequency
$\nu=\nu_{\mathrm{ref}}$ and $z=0$:
$j_{\nu,\mathrm{c}}\left(\nu_{\mathrm{ref}},0\right)\coloneqq\mathcal{C}j_{\nu,\mathrm{f}}\left(\nu_{\mathrm{ref}},0\right)$.
To retain minimal complexity, we propose to describe the spectral and temporal
dependencies of both $j_{\nu,\mathrm{f}}\left(z\right)$ and
$j_{\nu,\mathrm{c}}\left(z\right)$ as power laws in
$\frac{\nu}{\nu_{\mathrm{ref}}}$ and $\left(1+z\right)$:

$\displaystyle j_{\nu,\mathrm{f}}\left(\nu,z\right)$
$\displaystyle=j_{\nu,\mathrm{f}}\left(\nu_{\mathrm{ref}},0\right)\left(\frac{\nu}{\nu_{\mathrm{ref}}}\right)^{\alpha_{\mathrm{f}}}\left(1+z\right)^{\beta_{\mathrm{f}}}$
(16) $\displaystyle j_{\nu,\mathrm{c}}\left(\nu,z\right)$
$\displaystyle=j_{\nu,\mathrm{c}}\left(\nu_{\mathrm{ref}},0\right)\left(\frac{\nu}{\nu_{\mathrm{ref}}}\right)^{\alpha_{\mathrm{c}}}\left(1+z\right)^{\beta_{\mathrm{c}}},$
(17)

with $\alpha_{\mathrm{f}}$, $\beta_{\mathrm{f}}$, $\alpha_{\mathrm{c}}$ and
$\beta_{\mathrm{c}}$ constants. Under this choice of MEC function
parametrisation, the expectation value of $X\left(\nu_{\mathrm{obs}},z\right)$
becomes

$\displaystyle\mathbb{E}\left(X\left(\nu_{\mathrm{obs}},z\right)\right)=p_{\mathrm{f}}\Delta
I_{\nu,\mathrm{f}}\left(\nu_{\mathrm{obs}},z\right)+p_{\mathrm{c-f}}\Delta
I_{\nu,\mathrm{c-f}}\left(\nu_{\mathrm{obs}},z\right)$ (18)
$\displaystyle=\left(1+z\right)^{-4}\left(\frac{w_{\mathrm{f}}}{l}\right)^{2}j_{\nu,\mathrm{f}}\left(\nu_{\mathrm{ref}},0\right)\cdot$
$\displaystyle\left(3l_{\mathrm{f}}\left(\frac{\nu_{\mathrm{obs}}}{\nu_{\mathrm{ref}}}\right)^{\alpha_{\mathrm{f}}}\left(1+z\right)^{\alpha_{\mathrm{f}}+\beta_{\mathrm{f}}}+\mathcal{C}w_{\mathrm{f}}\left(\frac{\nu_{\mathrm{obs}}}{\nu_{\mathrm{ref}}}\right)^{\alpha_{\mathrm{c}}}\left(1+z\right)^{\alpha_{\mathrm{c}}+\beta_{\mathrm{c}}}\right).$
(19)

###  3.3 Specific-intensity-weighted mean redshift

Consider the random vector $\left[X_{1},X_{2},...,X_{N}\right]^{T}$. Whenever
$\max\\{X_{1},X_{2},...,X_{N}\\}>0$, we can define the MASSCW specific-
intensity-weighted mean redshift RV

$\displaystyle\bar{Z}\left(\nu_{\mathrm{obs}}\right)\coloneqq\frac{\sum_{n=1}^{N}z_{n}X_{n}\left(\nu_{\mathrm{obs}}\right)}{\sum_{n=1}^{N}X_{n}\left(\nu_{\mathrm{obs}}\right)}.$
(20)

We stress that $\bar{Z}$, being a deterministic function of the RVs
$\\{X_{1},X_{2},...,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 $\bar{Z}$.
Because $X_{n}\left(\nu_{\mathrm{obs}}\right)$ features in both numerator and
denominator of the fraction in Equation 20, redshift-independent
multiplicative factors cancel out. Therefore,
$j_{\nu,\mathrm{f}}\left(\nu_{\mathrm{ref}},0\right)$ need not be specified;
also, $\bar{Z}$ becomes $\nu_{\mathrm{obs}}$-independent if
$\alpha_{\mathrm{f}}=\alpha_{\mathrm{c}}$. In general, the distribution of
$\bar{Z}$ is determined by the cosmological model $\mathfrak{M}$ and 9
additional parameters ($w_{\mathrm{f}}$, $l$, $d_{1}$, $z_{\mathrm{max}}$,
$\alpha_{\mathrm{f}}$, $\beta_{\mathrm{f}}$, $\alpha_{\mathrm{c}}$,
$\beta_{\mathrm{c}}$ and $\mathcal{C}$). Using $d_{1}$, 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_{1}\sim\mathrm{Uniform}\left(0,l\right)$, making it an RV. In such case, 8
parameters remain that determine the distribution of $\bar{Z}$.

###  3.4 Flux-density-weighted mean redshift

The geometric model also allows to calculate the flux-density-weighted mean
(i.e. sky-averaged) MASSCW redshift
$\bar{\bar{z}}\left(\nu_{\mathrm{obs}}\right)$. We do so by considering
$M\in\mathbb{N}_{\geq 1}$ sightlines (instead of one) and their associated
specific-intensity-weighted mean redshift RVs
$\\{\bar{Z}_{1},\bar{Z}_{2},...,\bar{Z}_{M}\\}$ in the sense of Equation 20.
In total, these RVs depend on a set of $N\cdot M$ discrete RVs $X_{nm}$, which
can be partitioned into $N$ subsets of $M$ IID RVs: $X_{nm}\sim X_{n}$. We
find $\bar{\bar{z}}$ by summing the elements of
$\\{\bar{Z}_{1},\bar{Z}_{2},...,\bar{Z}_{M}\\}$, each weighted by the
corresponding specific intensity, and dividing the result by the sum of these
weights; all whilst we let $M\to\infty$:

$\displaystyle\bar{\bar{z}}\left(\nu_{\mathrm{obs}}\right)$
$\displaystyle\coloneqq\lim_{M\to\infty}\frac{\sum_{m=1}^{M}\sum_{n=1}^{N}z_{n}X_{nm}\left(\nu_{\mathrm{obs}}\right)}{\sum_{m=1}^{M}\sum_{n=1}^{N}X_{nm}\left(\nu_{\mathrm{obs}}\right)}$
(21)
$\displaystyle=\lim_{M\to\infty}\frac{\sum_{n=1}^{N}\left(z_{n}M^{-1}\sum_{m=1}^{M}X_{nm}\left(\nu_{\mathrm{obs}}\right)\right)}{\sum_{n=1}^{N}\left(M^{-1}\sum_{m=1}^{M}X_{nm}\left(\nu_{\mathrm{obs}}\right)\right)}$
(22)
$\displaystyle=\frac{\sum_{n=1}^{N}z_{n}\mathbb{E}\left(X_{n}\left(\nu_{\mathrm{obs}}\right)\right)}{\sum_{n=1}^{N}\mathbb{E}\left(X_{n}\left(\nu_{\mathrm{obs}}\right)\right)}.$
(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
$M^{-1}\sum_{m=1}^{M}X_{nm}\to\mathbb{E}\left(X_{n}\right)$ as $M\to\infty$,
for all $n\in\\{1,2,...,N\\}$. Like before, as
$\mathbb{E}\left(X_{n}\left(\nu_{\mathrm{obs}}\right)\right)$ features in both
numerator and denominator of the fraction in Equation 23, redshift-independent
multiplicative factors cancel out. Thus, $\bar{\bar{z}}$ is
$\nu_{\mathrm{obs}}$-independent when
$\alpha_{\mathrm{f}}=\alpha_{\mathrm{c}}$.
Equation 19 makes clear that clusters dominate over filaments at cosmological
redshift $z$ when

$\displaystyle\mathcal{C}\gg
3\frac{l_{\mathrm{f}}}{w_{\mathrm{f}}}\left(\frac{\nu_{\mathrm{obs}}}{\nu_{\mathrm{ref}}}\right)^{\alpha_{\mathrm{f}}-\alpha_{\mathrm{c}}}\left(1+z\right)^{\alpha_{\mathrm{f}}-\alpha_{\mathrm{c}}+\beta_{\mathrm{f}}-\beta_{\mathrm{c}}}.$
(24)

At $\nu_{\mathrm{obs}}=\nu_{\mathrm{ref}}$ and $z=0$, this inequality reduces
to

$\displaystyle\mathcal{C}\gg 3\left(\frac{l}{w_{\mathrm{f}}}-1\right),$ (25)

which is amply satisfied (see Section 3.5 for parameter estimates). Near
$\nu_{\mathrm{obs}}=\nu_{\mathrm{ref}}$ and for $z$ of order unity, the
inequality continues to hold as long as
$\alpha_{\mathrm{f}}-\alpha_{\mathrm{c}}+\beta_{\mathrm{f}}-\beta_{\mathrm{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
$\bar{\bar{z}}$. Therefore, $\bar{\bar{z}}$ around
$\nu_{\mathrm{obs}}=\nu_{\mathrm{ref}}$ is typically dominated by clusters,
and $\bar{\bar{z}}\left(\nu_{\mathrm{obs}}\right)\approx\bar{\bar{z}}$. As
long as Equation 24 holds, $\bar{\bar{z}}$ becomes insensitive to variations
in $\mathcal{C}$, $w_{\mathrm{f}}$, $\alpha_{\mathrm{f}}$ and
$\beta_{\mathrm{f}}$; setting $z_{\mathrm{max}}=\infty$, just _3_ significant
parameters remain: $l$, $d_{1}$ and $\alpha_{\mathrm{c}}+\beta_{\mathrm{c}}$.
Finally, it is of significant interest to calculate $\bar{\bar{z}}$ for the
filaments’ SCW signal _only_ , discarding the dominant SCW signal from merger
and accretion shocks in galaxy clusters. To find a filament-only
$\bar{\bar{z}}$ with the geometric model, we just set $\mathcal{C}=0$. This
simplifies the expression for
$\mathbb{E}\left(X\left(\nu_{\mathrm{obs}},z\right)\right)$, and reduces the
number of relevant parameters considerably. When $\mathcal{C}=0$,
$\bar{\bar{z}}$ becomes not only independent of $\alpha_{\mathrm{c}}$ and
$\beta_{\mathrm{c}}$, but also of $w_{\mathrm{f}}$ and $\nu_{\mathrm{obs}}$:

$\displaystyle\bar{\bar{z}}\left(\nu_{\mathrm{obs}}\right)=\bar{\bar{z}}=\frac{\sum_{n=1}^{N}z_{n}\left(1+z_{n}\right)^{-4+\alpha_{\mathrm{f}}+\beta_{\mathrm{f}}}}{\sum_{n=1}^{N}\left(1+z_{n}\right)^{-4+\alpha_{\mathrm{f}}+\beta_{\mathrm{f}}}}.$
(26)

As before, if we set $z_{\mathrm{max}}=\infty$, $\bar{\bar{z}}$ is a function
of just _3_ parameters: $l$, $d_{1}$ and
$\alpha_{\mathrm{f}}+\beta_{\mathrm{f}}$, in this case.

###  3.5 Parameter estimates

Figure 7: The monochromatic-emission-coefficient-weighted mean integrated
spectral index $\bar{\alpha}$ at an emission frequency of $\nu=150\
\mathrm{MHz}$, including all shocks that occur at locations where the dark
matter (DM) density $\rho_{\mathrm{DM}}$ does not exceed some threshold.
Larger thresholds incorporate shocks from a wider range in
$\rho_{\mathrm{DM}}$, eventually including shocks from both filament and
cluster environments. To provide an idea of the uncertainty in $\bar{\alpha}$,
we use jackknife realisations that each contain $80\%$ of the shocked voxels
present in our Enzo snapshot at $z=0.025$.

This subsection provides concrete parameter estimates, which are necessary to
run the model in practice.
Cosmological simulations indicate $w_{\mathrm{f}}\sim 10^{0}-10^{1}\
\mathrm{Mpc}$, while $l\sim 10^{1}-10^{2}\ \mathrm{Mpc}$. The ratio of these
parameters is constrained too: Table 1 shows that a VFF comparison between the
geometric model and cosmological simulations favours
$\frac{w_{\mathrm{f}}}{l}\sim 10^{-1}$.
From Figure 3, we estimate $\mathcal{C}\sim 10^{5}-10^{7}$ for
$\nu_{\mathrm{ref}}=150\ \mathrm{MHz}$, assuming filament environments are
characterised by $\rho\sim\bar{\rho}_{0}$, and cluster environments are
characterised by $\rho\sim 10^{2}\ \bar{\rho}_{0}$ (at the ${\sim}3\
\mathrm{Mpc}\ h^{-1}$ resolution, at least). Technically, Figure 3 does not
show $P\left(j_{\nu}\ |\ \rho\right)$ for $z=0$ (the redshift at which
$\mathcal{C}$ 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 $\alpha_{\mathrm{f}}$ and
$\alpha_{\mathrm{c}}$, we calculate the MEC-weighted mean integrated spectral
index $\bar{\alpha}$ from our Enzo simulation snapshot as a function of shock
environment. Figure 7 shows the result when weighing by MECs at emission
frequency $\nu=150\ \mathrm{MHz}$; in general,
$\bar{\alpha}=\bar{\alpha}\left(\nu\right)$. Nevertheless, it seems that
$\alpha_{\mathrm{f}}=\alpha_{\mathrm{c}}=-\nicefrac{{3}}{{2}}$ are reasonable
choices.
Next, to find $\beta_{\mathrm{f}}$ and $\beta_{\mathrm{c}}$, we revisit the
single-shock MEC expression of Equation 2. Inverse Compton (IC) scattering to
the CMB contributes a factor $\left(1+z\right)^{-4}$ to the single-shock MEC
as $j_{\nu}\propto B_{\mathrm{CMB}}^{-2}$ and
$B_{\mathrm{CMB}}\propto\left(1+z\right)^{2}$ (see Appendix B). Further
factors of $\left(1+z\right)$ follow by considering the typical redshift
evolution of the proper BM density $\rho_{\mathrm{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
$\left(1+z\right)^{-1}$. As the proper _mean_ matter density is proportional
to $\left(1+z\right)^{3}$, the proper BM density in filaments is expected to
be proportional to $\left(1+z\right)^{2}$. In clusters, the density contrast
grows more strongly; numerical simulations suggest that the density contrast
is roughly proportional to $\left(1+z\right)^{-2}$. Likewise, this would imply
that the proper BM density in clusters is proportional to $\left(1+z\right)$.
Assuming that the magnetic field strength only evolves due to field line
compression, we have $B\propto\rho_{\mathrm{BM}}^{\nicefrac{{2}}{{3}}}$ (see
Appendix B).181818Ignoring magnetogenesis by outflows from AGN and supernovae,
this Ansatz likely underestimates the actual magnetic field strength growth
with time. In filaments then,
$B_{\mathrm{f}}\propto\left(1+z\right)^{\nicefrac{{4}}{{3}}}$, whilst in
clusters $B_{\mathrm{c}}\propto\left(1+z\right)^{\nicefrac{{2}}{{3}}}$. As the
single-shock MEC is proportional to $B^{1-\alpha}$, the proper magnetic field
strength approximately contributes a factor
$\left(1+z\right)^{\nicefrac{{10}}{{3}}}$ for filaments, and a factor
$\left(1+z\right)^{\nicefrac{{5}}{{3}}}$ for clusters. Lacking further
knowledge, we apply Occam’s razor and assume that the proper shock velocity
and comoving shock number density do not evolve with redshift.
Collecting factors of $\left(1+z\right)$, we find
$\beta_{\mathrm{f}}=\nicefrac{{4}}{{3}}$ and
$\beta_{\mathrm{c}}=-\nicefrac{{4}}{{3}}$. We stress that these values are
highly uncertain and should be used as indications only.

##  4 Results

(a)

(b)

Figure 8: Realisation of the merger- and accretion-shocked synchrotron Cosmic
Web (MASSCW) priors at observing frequency $\nu_{\mathrm{obs}}=150\
\mathrm{MHz}$. We show exactly $\nicefrac{{1}}{{4}}$ of the total sky. Top:
specific intensity $I_{\nu}$. Bottom: specific-intensity-weighted mean
(cosmological) redshift $\bar{z}$.

###  4.1 MASSCW specific intensity function distribution: general results

All the results given are for $\nu_{\mathrm{obs}}=150\ \mathrm{MHz}$. At the
BORG SDSS LSS reconstruction resolution, realisations from our MASSCW specific
intensity function distribution exhibit a factor ${\sim}10^{3}$ 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 Figure 8(a), where we show a single realisation
of $I_{\nu}\left(\hat{r},\nu_{\mathrm{obs}}\right)$ 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.
In Figure 8(b), we show the corresponding MASSCW specific-intensity-weighted
mean redshift function realisation
$\bar{z}\left(\hat{r},\nu_{\mathrm{obs}}\right)$. All redshifts are within
$[0,0.2]$, because the BORG SDSS reconstructions are limited to
$z_{\mathrm{max}}=0.2$. A comparison with Figure 8(a) reveals that generally,
sky patches of high MASSCW specific intensity are due to structures at low
redshift, and vice versa.
In Figure 9(a) and 9(b), we show the mean specific intensity function
$\mu_{I_{\nu}}\left(\hat{r},\nu_{\mathrm{obs}}\right)$ and mean specific-
intensity-weighted mean redshift function
$\mu_{\bar{z}}\left(\hat{r},\nu_{\mathrm{obs}}\right)$ at
$\nu_{\mathrm{obs}}=150\ \mathrm{MHz}$ for each of the $10^{6}$ ray-traced
directions, pooling all $10^{3}$ realisations.

(a)

(b)

Figure 9: Per-direction mean of the MASSCW priors at observing frequency
$\nu_{\mathrm{obs}}=150\ \mathrm{MHz}$, exhibiting less variability than the
single realisation of Figures 8(a) and 8(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.

Both of these functions represent a summary statistic calculated from each
ray’s marginal $I_{\nu}$\- and $\bar{z}$-distribution (i.e. 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.
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
$m_{I_{\nu}}\left(\hat{r},\nu_{\mathrm{obs}}\right)$ (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 $5^{\mathrm{th}}$ to $95^{\mathrm{th}}$
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 Figure 14(a) and Figure 14(b), we present both an absolute and
relative measure of spread, again calculated from the MASSCW specific
intensity prior marginals. Figure 14(a) shows the standard deviation
$\sigma_{I_{\nu}}\left(\hat{r},\nu_{\mathrm{obs}}\right)$ 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 $\sigma_{I_{\nu}}\mu_{I_{\nu}}^{-1}$, which we
show in Figure 14(b). In this sense, the specific intensity function can be
most accurately predicted around galaxy clusters and superclusters, reaching
$\sigma_{I_{\nu}}\sim 60\%\ \mu_{I_{\nu}}$. As expected, the relative
prediction uncertainty is highest outside the SDSS DR7 coverage (i.e. near the
$\phi=90\degree$ and $\phi=270\degree$ 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.

Figure 10: Three realisations (rows) of the merger- and accretion-shocked
synchrotron Cosmic Web (MASSCW) specific intensity prior at observing
frequency $\nu_{\mathrm{obs}}=150\ \mathrm{MHz}$, showing zoom-ins around
three massive Northern Sky galaxy clusters (columns). The zoom-ins show caps
of the celestial sphere of $10\degree$ radius. The dashed circles are at
$3\degree$, $6\degree$ and $9\degree$ 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 ${\sim}3\ \mathrm{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 ${\sim}1\ \mathrm{Mpc}$ MEC - total matter density
(anti)correlations are not resolved.

###  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 Figure 10, we show $I_{\nu}$ over spherical domes with an angular
radius of $12\degree$, whilst Figure 15 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.
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 ($\geq$ 50-hr) observations have been
conducted that complement the wide-field, 8-hr 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 Figure 11, we show the MASSCW specific
intensity prior single-direction (i.e. marginal) medians for three such deep
fields.191919For these close-ups, we present the median instead of the mean to
emphasise that our probabilistic approach enables the calculation of a variety
of useful summary statistics. 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.
The Lockman Hole (for which a 100-hr dataset is available (Tasse et al., in
prep.)) 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-hr
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-hr 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 Figure 10 and 15. For this reason, it
has been selected for deep LOFAR HBA observations.

Figure 11: Merger- and accretion-shocked synchrotron Cosmic Web (MASSCW)
specific intensity prior marginal medians at observing frequency
$\nu_{\mathrm{obs}}=150\ \mathrm{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.

###  4.4 Specific-intensity-weighted mean redshift

In the top panel of Figure 12, we put our geometric model to the test by
comparing a predicted specific-intensity-weighted mean redshift CDF
$F_{\bar{Z}}\left(\bar{z}\right)$ (green dotted line) to $10^{2}$ randomly
drawn ECDFs from our $10^{3}$ MASSCW prior realisations (translucent black
lines), at $\nu_{\mathrm{obs}}=150\ \mathrm{MHz}$. Both prior realisations and
the model reach a cumulative probability (‘fraction of the sky’) of $1$ at
$\bar{z}=0.2$: the BORG SDSS samples do not feature LSS beyond
$z_{\mathrm{max}}=0.2$ and the model is restricted accordingly. The model CDF
is constructed from $10^{5}$ draws from RV $\bar{Z}$ (see Equation
20).202020In a strict sense, given the finite number of realisations in our
numerical approximation, this function is also an ECDF. However, for all
practical purposes, it can be regarded as a CDF. We adopt parameters suggested
by Section 3.5: $w_{\mathrm{f}}=5\ \mathrm{Mpc}$, $l=50\ \mathrm{Mpc}$,
$d_{1}=25\ \mathrm{Mpc}$, $\mathcal{C}=10^{6}$ at $\nu_{\mathrm{ref}}=150\
\mathrm{MHz}$, $\alpha_{\mathrm{f}}=\alpha_{\mathrm{c}}=-\nicefrac{{3}}{{2}}$,
$\beta_{\mathrm{f}}=\nicefrac{{4}}{{3}}$ and
$\beta_{\mathrm{c}}=-\nicefrac{{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 $\bar{Z}$ for the
_true_ sky, which features LSS beyond $z_{\mathrm{max}}=0.2$. For the
observational study of filaments, it is of greatest interest to calculate
$F_{\bar{Z}}\left(\bar{z}\right)$ when the cluster contribution to the MASSCW
is ignored, as it would otherwise dominate. In the bottom panel of Figure 12,
we therefore show $F_{\bar{Z}}\left(\bar{z}\right)$ according to the geometric
model for $z_{\mathrm{max}}=\infty$ and $\mathcal{C}=0$. When $\mathcal{C}=0$,
$\alpha_{\mathrm{c}}$ and $\beta_{\mathrm{c}}$ become irrelevant. Four
parameters remain: $\alpha_{\mathrm{f}}+\beta_{\mathrm{f}}$, $l$,
$w_{\mathrm{f}}$ and $d_{1}$. The first of these has by far the most effect on
$F_{\bar{Z}}\left(\bar{z}\right)$, as evidenced by the differently coloured
curves (we keep $\alpha_{\mathrm{f}}=-\nicefrac{{3}}{{2}}$ constant and vary
$\beta_{\mathrm{f}}$, but varying $\alpha_{\mathrm{f}}$ and keeping
$\beta_{\mathrm{f}}$ constant would lead to identical results). We also vary
the purely geometric parameters $l$, $w_{\mathrm{f}}$ and $d_{1}$, but in such
a way that $w_{\mathrm{f}}=\tfrac{1}{10}l$ (as suggested by comparing the
model’s VFFs to those from cosmological simulations) and $d_{1}=\tfrac{1}{2}l$
(i.e. the observer is always put at the centre of a cubic unit cell). When
$l\in\\{50\ \mathrm{Mpc},75\ \mathrm{Mpc},100\ \mathrm{Mpc}\\}$ is increased,
and $w_{\mathrm{f}}$ and $d_{1}$ accordingly, modest changes in
$F_{\bar{Z}}\left(\bar{z}\right)$ occur (dash-dotted and dotted lines): the
distribution of $\bar{Z}$ attains larger spread.
$F_{\bar{Z}}\left(\bar{z}\right)$ provides detailed information. For example,
the model with $\beta_{\mathrm{f}}=\nicefrac{{4}}{{3}}$ and $l=50\
\mathrm{Mpc}$ predicts that for $\sim 50\%$ of the sky (${\sim}21,000\
\mathrm{sq.\ deg.}$), the MASSCW signal has a mean redshift of $0.35$ or
lower, whilst for $\sim 80\%$ of the sky (${\sim}33,000\ \mathrm{sq.\ deg.}$),
the mean redshift is $0.42$ or lower.

Figure 12: MASSCW specific-intensity-weighted mean redshift RV $\bar{Z}$
results. Top: distributions of $\bar{Z}$, 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 $z_{\mathrm{max}}=0.2$ only. The
other geometric model parameters are: $w_{\mathrm{f}}=5\ \mathrm{Mpc}$, $l=50\
\mathrm{Mpc}$, $d_{1}=25\ \mathrm{Mpc}$, $\mathcal{C}=10^{6}$,
$\alpha_{\mathrm{f}}=\alpha_{\mathrm{c}}=-\nicefrac{{3}}{{2}}$,
$\beta_{\mathrm{f}}=\nicefrac{{4}}{{3}}$ and
$\beta_{\mathrm{c}}=-\nicefrac{{4}}{{3}}$. Bottom: geometric model CDFs of
$\bar{Z}$ for filaments only ($\mathcal{C}=0$), though including all of them
($z_{\mathrm{max}}=\infty$). By far the most influential parameter is also the
most uncertain one; we therefore vary $\beta_{\mathrm{f}}$ over a plausible
range. We also vary $l\in\\{50\ \mathrm{Mpc},75\ \mathrm{Mpc},100\
\mathrm{Mpc}\\}$, and by extension also $w_{\mathrm{f}}$ and $d_{1}$, by
forcing $w_{\mathrm{f}}=\tfrac{1}{10}l$ and $d_{1}=\tfrac{1}{2}l$. As $l$
increases, the unit cell number density decreases; we symbolise this with more
sparsely drawn curves. We keep $\alpha_{\mathrm{f}}=-\nicefrac{{3}}{{2}}$
constant.

###  4.5 Flux-density-weighted mean redshift

The flux-density-weighted mean redshift $\bar{\bar{z}}$ constitutes the most
concise measure of the typical MASSCW signal epoch. In the top panel of Figure
13, we show a distribution (shaded grey) over $\bar{\bar{z}}$ generated via
KDE (Gaussian kernel, $\sigma_{\mathrm{KDE}}=4\cdot 10^{-4}$) from our
$10^{3}$ MASSCW prior samples. The median is $\bar{\bar{z}}=0.077$ (solid
grey); however, we stress that this is because only LSS up to
$z_{\mathrm{max}}=0.2$ is included. We present an example model (solid green)
that reproduces the median $\bar{\bar{z}}$, with main parameters
$z_{\mathrm{max}}=0.2$, $l=75\ \mathrm{Mpc}$, $d_{1}=60\ \mathrm{Mpc}$,
$\alpha_{\mathrm{c}}=-\nicefrac{{3}}{{2}}$ and
$\beta_{\mathrm{c}}=-\nicefrac{{4}}{{3}}$. (Adopting $\mathcal{C}=10^{6}$ at
$\nu_{\mathrm{ref}}=150\ \mathrm{MHz}$ as suggested by Section 3.5, clusters
dominate over filaments, and the other parameters play a very minor role.) To
explore the sensitivity of $\bar{\bar{z}}$ on the parameters, we generate
$5\cdot 10^{5}$ parameter sets by drawing from wide uniform distributions. We
draw $l\sim\mathrm{Uniform}\left(50\ \mathrm{Mpc},100\ \mathrm{Mpc}\right)$,
$d_{1}\sim\mathrm{Uniform}\left(0,l\right)$,
$w_{\mathrm{f}}\sim\mathrm{Uniform}\left(5\ \mathrm{Mpc},10\
\mathrm{Mpc}\right)$,
$\log_{10}\left(\mathcal{C}\right)\sim\mathrm{Uniform}\left(5,7\right)$,
$\alpha_{\mathrm{f}},\alpha_{\mathrm{c}}\sim\mathrm{Uniform}\left(-\nicefrac{{3}}{{2}}-\nicefrac{{1}}{{4}},-\nicefrac{{3}}{{2}}+\nicefrac{{1}}{{4}}\right)$,
$\beta_{\mathrm{f}}\sim\mathrm{Uniform}\left(\nicefrac{{4}}{{3}}-2,\nicefrac{{4}}{{3}}+2\right)$
and
$\beta_{\mathrm{c}}\sim\mathrm{Uniform}\left(-\nicefrac{{4}}{{3}}-2,-\nicefrac{{4}}{{3}}+2\right)$.
However, the resultant distribution for $\bar{\bar{z}}$ is mostly restricted
to the relatively small range of $\left(0.06,0.08\right)$. We stress that the
two distributions are conceptually distinct and are not meant to be compared
directly.
In the bottom panel of Figure 13, we show $\bar{\bar{z}}$ for filaments only
($\mathcal{C}=0$) and $z_{\mathrm{max}}=\infty$. The prediction then becomes a
function of three parameters: $\alpha_{\mathrm{f}}+\beta_{\mathrm{f}}$, $l$
and $d_{1}$. The dependencies on both $l$ and $d_{1}$ are weak, and
$\bar{\bar{z}}$ 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_{1}$ on $\bar{\bar{z}}$
increases when $l$ is larger. As the unit cells grow (so that their number
density drops $\propto l^{-3}$), $\bar{\bar{z}}$ becomes more sensitive to
contributions of individual boundary crossings (filaments), including the
first. Filament specific intensity contributions decrease quickly (especially
for low values of $\alpha_{\mathrm{f}}+\beta_{\mathrm{f}}$) and monotonically
with distance, so that the distance of the first filament to the observer
(e.g. nearby, $d_{1}=10\%\ l$, or far away, $d_{1}=90\%\ l$) is of some
importance.

Figure 13: MASSCW flux-density-weighted mean redshift $\bar{\bar{z}}$ results.
Top: distribution of $\bar{\bar{z}}$ calculated via KDE from our prior’s
$10^{3}$ realisations of $I_{\nu}$ and $\bar{z}$ (shaded grey), which are
limited to $z_{\mathrm{max}}=0.2$. The median is $\bar{\bar{z}}=0.077$ (solid
grey). Furthermore, we show the variety of geometric model predictions for
$\bar{\bar{z}}$ when similarly limited to $z_{\mathrm{max}}=0.2$, assuming
flat priors on all model parameters (shaded green; see text for details). The
median is $\bar{\bar{z}}=0.072$. 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\ \mathrm{Mpc}$, $d_{1}=60\
\mathrm{Mpc}$ and $\beta_{\mathrm{c}}=-\nicefrac{{4}}{{3}}$. Bottom:
predictions for $\bar{\bar{z}}$ considering filaments _only_ , varying all
three relevant geometric model parameters:
$\alpha_{\mathrm{f}}+\beta_{\mathrm{f}}$, $l$ and $d_{1}$. The two
$d_{1}=50\%\ l$ curves are virtually indistinguishable.

##  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 structures, or ‘bridges’, 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 length ${\sim}2\ \mathrm{Mpc}$)
(Botteon et al., 2018, 2020b), whilst the first unambiguous detection was in
Abell 399–401 (of length ${\sim}3\ \mathrm{Mpc}$) (Govoni et al., 2019). The
latter ridge features higher densities (${\sim}3\cdot 10^{2}\
\mathrm{m}^{-3}$), temperatures (${\sim}7\cdot 10^{7}\ \mathrm{K}$), and
magnetic field strengths (${\sim}1\ \mathrm{\mu 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\ \mathrm{Mpc}$
extent are smaller than the current BORG SDSS resolution of ${\sim}3\
\mathrm{Mpc}\ h^{-1}$, 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
$\rho_{\mathrm{BM}}\left(\rho\right)=\frac{\Omega_{\mathrm{BM,0}}}{\Omega_{\mathrm{M,0}}}\rho$
and
$T\left(\rho_{\mathrm{BM}}\right)=T_{\mathrm{ref}}\frac{\rho_{\mathrm{BM}}}{\rho_{\mathrm{BM,ref}}}$
(with $T_{\mathrm{ref}}=10^{6}\ \mathrm{K}$ and $\rho_{\mathrm{BM,ref}}=m_{p}\
\mathrm{m^{-3}}$, 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
Figure 16, for exactly the same sky regions as in Figure 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 (${\sim}10\degree$) 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 Section 2.4. Currently, the
BORG SDSS voxel length is larger than the shock correlation length
$l_{\mathrm{SE}}=2\ \mathrm{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_{\mathrm{SE}}$, instead of by the density reconstruction resolution.
A straightforward next step would be to use density reconstructions of the
remaining $\nicefrac{{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_{\mathrm{max}}=0.1$ instead of
$z_{\mathrm{max}}=0.2$.
Reconstructions that push beyond $z_{\mathrm{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
$F_{\bar{Z}}\left(\bar{z}\right)$ for $z_{\mathrm{max}}=\infty$. However, the
single most influential parameter governing $F_{\bar{Z}}\left(\bar{z}\right)$,
$\alpha_{\mathrm{f}}+\beta_{\mathrm{f}}$, is currently ill-
constrained.212121Interestingly, its determination could be done with existing
cosmological simulations. 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\ \mathrm{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. (2020) 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 $\left(100\
\mathrm{Mpc}\right)^{3}$ 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 Figure 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

Based on our predictions, what are the observational prospects of detecting
the filament SCW — and the filament _MASSCW_ in particular? As evinced by
radio bridge detections (Botteon et al., 2018; Govoni et al., 2019; Botteon et
al., 2020b) 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_{\mathrm{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.222222As discussed in Section 5.1, the MASSCW signal is
unlikely to be the dominant contributor to the SCW signal _in radio bridges_ ;
in the case of the all-sky search of Vernstrom et al. (2021), the dominant
emission mechanism remains unknown. For associated MASSCW predictions, see
Figures 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 Section 1.4.3); secondly, the filling factor of shocks of
appropriate strength and obliquity to trigger DSA; and, thirdly, the magnitude
of the electron acceleration efficiency $\xi_{e}$ in filaments (see Section
2.3.2). From Figure 10, we see that the specific intensity function in the
direction of filaments around massive low-redshift galaxy clusters, at
$\nu_{\mathrm{obs}}=150\ \mathrm{MHz}$ and degree-scale resolution, reaches
$I_{\nu}\sim 10^{-1}\ \mathrm{Jy\ deg^{-2}}$ for $\xi_{e}=1$. If $\xi_{e}\sim
10^{-2}$ in filaments, as in Keshet et al. (2004)’s SNR shocks, the actual
specific intensity is $I_{\nu}\sim 10^{-3}\ \mathrm{Jy\ deg^{-2}}$. A
conservative estimate of the uncertainty in $\xi_{e}$ is an order of
magnitude, yielding a range of specific intensity estimates $I_{\nu}\sim
10^{-4}\text{--}10^{-2}\ \mathrm{Jy\ deg^{-2}}$.
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 $\nu_{\mathrm{obs}}\ {\sim}10^{2}\
\mathrm{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^{\prime\prime}$ and $90^{\prime\prime}$) images with the LoTSS DR2
(Shimwell et al., 2022) that reveal the Milky Way’s specific intensity
function at $\nu_{\mathrm{obs}}=144\ \mathrm{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 $\sigma_{I_{\nu}}\sim 10^{0}\text{--}10^{1}\
\mathrm{Jy\ deg^{-2}}$ 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-hr 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,232323This estimate follows from considering the spherical harmonics
angular power spectrum $C_{\ell}$ of the interferometrically observed (and
thus large-angular-scale–deprived) synchrotron Milky Way. We model $C_{\ell}$
as a power law (La Porta et al., 2008) from degree $\ell_{0}$ onwards:
$\displaystyle C_{\ell}\left(\ell\right)=\cases{0}$ ℓ<ℓ0; (27) (28) where
$\beta<-1$. Generally, given an angular power spectrum, the total power $P$ is
$\displaystyle P=\sum_{\ell=0}^{\infty}C_{\ell}\left(\ell\right).$ (29) We
model the imposition of an additional inner $(u,v)$-cut as the removal of
power on all angular scales up to (but excluding) $\ell_{0}^{\prime}$, with
$\ell_{0}^{\prime}\geq\ell_{0}$. The (negative) relative change in total power
caused by the $(u,v)$-cut is
$\displaystyle\frac{P^{\prime}-P}{P}=\frac{\sum_{\ell_{0}^{\prime}}^{\infty}\ell^{\beta}}{\sum_{\ell_{0}}^{\infty}\ell^{\beta}}-1=\frac{\zeta\left(-\beta,\ell_{0}^{\prime}\right)}{\zeta\left(-\beta,\ell_{0}\right)}-1,$
(30) where $\zeta$ is the Hurwitz zeta function. For $\beta\in[-3,-2]$ (e.g.
La Porta et al., 2008; Ghosh et al., 2012; Sims et al., 2016), $\ell_{0}\sim
10^{2}$ and $\frac{\ell_{0}^{\prime}}{\ell_{0}}\sim 10^{0}$, we find relative
total power changes of $-90\%$ to $-99\%$.  leaving a Milky Way contamination
of $\sigma_{I_{\nu}}\sim 10^{-2}\text{--}10^{0}\ \mathrm{Jy\ deg^{-2}}$. (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^{-4}$ in a very
pessimistic case, $10^{-1}$ in a fairly optimistic case, and $10^{0}$ 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^{1}$ degree diameter to achieve a signal-to-
noise ratio of order unity (e.g. three). For example, the Ursa Major
Supercluster (see the bottom panel of Figure 11) is a region of roughly the
required extent. In case the MASSCW specific intensity is an order of
magnitude weaker, i.e. of order $I_{\nu}\sim 10^{-3}\ \mathrm{Jy\ deg^{-2}}$,
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
$I_{\nu}\left(\hat{r},\nu_{\mathrm{obs}}\right)$ can be generated from a
percentile random field $\mathcal{P}\left(\mathbf{r}\right)$ (or,
equivalently, a Gaussian random field $\mathcal{Z}\left(\mathbf{r}\right)$)
and a total matter density random field $\rho\left(\mathbf{r}\right)$;
$\hat{r}\in\mathbb{S}^{2}$, $\nu_{\mathrm{obs}}\in\mathbb{R}_{>0}$, and
$\mathbf{r}\in\mathcal{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_{\nu}=f\left(\mathcal{P},\rho\right)$. Implicitly, our approach has been to
sample $\rho$ from the BORG SDSS posterior, and $\mathcal{P}$ from another,
_independent_ distribution (which is fixed by the distribution of
$\mathcal{Z}$, and thus specified by the covariance function
$K_{\mathrm{SE}}$). Under what physical scenario is this justified? Consider
the joint distribution for the input random fields
$P\left(\mathcal{P},\rho\right)$, which can be written as the product of a
conditional and marginal: $P\left(\mathcal{P},\rho\right)=P\left(\mathcal{P}\
|\ \rho\right)P\left(\rho\right)$. Thus, the central assumption that underlies
our sampling approach is that $P\left(\mathcal{P}\ |\ \rho\right)\approx
P\left(\mathcal{P}\right)$, so that $P\left(\mathcal{P},\rho\right)\approx
P\left(\mathcal{P}\right)P\left(\rho\right)$. 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 Section 2.5.1, we calculate $I_{\nu}$ by generating $j_{\nu}$ 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
$[\nu_{\mathrm{obs}},\nu_{\mathrm{obs}}\left(1+z_{\mathrm{max}}\right)]$. Does
this assumption correspond to a plausible physical scenario? Let the MECs of
two shocks at some reference emission frequency $\nu=\nu_{\mathrm{ref}}$ be
$j_{1}$ and $j_{2}$, and let $j_{1}>j_{2}$ without loss of generality. Assume
the existence of some function $\alpha=\alpha\left(j_{\nu}\right)$, that
assigns integrated spectral indices to shocks based on their MEC at
$\nu=\nu_{\mathrm{ref}}$. At emission frequency $\nu$, the MEC ranks of the
shocks are the same as at $\nu_{\mathrm{ref}}$ if and only if

$\displaystyle
j_{1}\left(\frac{\nu}{\nu_{\mathrm{ref}}}\right)^{\alpha\left(j_{1}\right)}>j_{2}\left(\frac{\nu}{\nu_{\mathrm{ref}}}\right)^{\alpha\left(j_{2}\right)}\
\mathrm{for\ all}\ j_{1},j_{2}\in\mathbb{R}_{\geq 0},\ j_{1}>j_{2}.$ (31)

The rewritten inequality

$\displaystyle\frac{j_{1}}{j_{2}}>\left(\frac{\nu}{\nu_{\mathrm{ref}}}\right)^{\alpha\left(j_{2}\right)-\alpha\left(j_{1}\right)}\
\mathrm{for\ all}\ j_{1},j_{2}\in\mathbb{R}_{\geq 0},\ j_{1}>j_{2}$ (32)

suggests $\alpha\left(j_{2}\right)-\alpha\left(j_{1}\right)\leq 0$, or
$\alpha\left(j_{1}\right)\geq\alpha\left(j_{2}\right)$ for
$\nu>\nu_{\mathrm{ref}}$: $\alpha\left(j_{\nu}\right)$ must be a monotonically
_increasing_ function. Analogously, for $\nu<\nu_{\mathrm{ref}}$, we find that
$\alpha\left(j_{\nu}\right)$ must be a monotonically _decreasing_ function.
These scenarios are visualised in Figure 17. Using the Enzo simulation data
used in this work, we explore the spectral index - MEC relation at
$\nu_{\mathrm{ref}}=180\ \mathrm{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 $\alpha\left(j_{\nu}\right)$; 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.242424A conceptually correct way to address this problem would be
to realise both a spectral index and a MEC at some fixed emission frequency
($\nu=\nu_{\mathrm{obs}}$, say), using $P\left(\alpha,j_{\nu}\ |\
\rho\right)=P\left(\alpha\ |\ j_{\nu},\rho\right)P\left(j_{\nu}\ |\
\rho\right)$. The conditional probability distribution $P\left(\alpha\ |\
j_{\nu},\rho\right)$ could be learnt from the Enzo simulation data used in
this work, too. Our formalism allows for the generation of the function
$I_{\nu}\left(\hat{r},\nu_{\mathrm{obs}}\right)$ 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

  1. 1.

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.

  2. 2.

Using BORG SDSS total matter density reconstructions and Enzo cosmological
simulations, we build a prior distribution that is informative over half of
the Northern Sky, and that has a ${\sim}0.6\degree$ 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 $I_{\nu}\
\xi_{e}^{-1}\sim 10^{-1}\ \mathrm{Jy\ deg^{-2}}$ at $\nu_{\mathrm{obs}}=150\
\mathrm{MHz}$; $\xi_{e}$ is the highly uncertain electron acceleration
efficiency. Even at the ${\sim}3\ \mathrm{Mpc}\ h^{-1}$ reconstruction
resolution, our merger and accretion shock SCW (MASSCW) prior indicates that
the specific intensity for a given direction is highly uncertain (with a
typical standard deviation being ${\sim}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.

  3. 3.

With a simple geometric model of cubic unit cells, we calculate both the
distribution of the specific-intensity-weighted mean redshift RV $\bar{Z}$, as
well as the flux-density-weighted mean redshift $\bar{\bar{z}}$ of the MASSCW
signal for large-scale structure (LSS) reconstructions up to
$z_{\mathrm{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, we present filament-only geometric model predictions for
$\bar{\bar{z}}$ that include _all_ LSS (i.e. $z_{\mathrm{max}}=\infty$). These
predictions are highly insensitive to plausible variations in model parameters
$l$ and $d_{1}$, demonstrating that $\bar{\bar{z}}$ is effectively determined
by a single parameter: the sum of the typical MEC-weighted filament spectral
index $\alpha_{\mathrm{f}}$ and the MEC - cosmological redshift power law
exponent $\beta_{\mathrm{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.

###### Acknowledgements.

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.

## References

  * Abazajian et al. (2009) Abazajian, K. N., Adelman-McCarthy, J. K., Agüeros, M. A., et al. 2009, ApJS, 182, 543
  * Araya-Melo et al. (2012) Araya-Melo, P. A., Aragón-Calvo, M. A., Brüggen, M., & Hoeft, M. 2012, MNRAS, 423, 2325
  * Axford et al. (1977) Axford, W. I., Leer, E., & Skadron, G. 1977, in International Cosmic Ray Conference, Vol. 11, International Cosmic Ray Conference, 132
  * Baring (1997) Baring, M. G. 1997, in Very High Energy Phenomena in the Universe; Moriond Workshop, ed. Y. Giraud-Heraud & J. Tran Thanh van, 97
  * Bell (1978a) Bell, A. R. 1978a, MNRAS, 182, 147
  * Bell (1978b) Bell, A. R. 1978b, MNRAS, 182, 443
  * Blandford & Eichler (1987) Blandford, R. & Eichler, D. 1987, Phys. Rep, 154, 1
  * Blandford & Ostriker (1978) Blandford, R. D. & Ostriker, J. P. 1978, ApJ, 221, L29
  * Bos et al. (2019) Bos, E. G. P., Kitaura, F.-S., & van de Weygaert, R. 2019, MNRAS, 488, 2573
  * Botteon et al. (2020a) Botteon, A., Brunetti, G., van Weeren, R. J., et al. 2020a, ApJ, 897, 93
  * Botteon et al. (2018) Botteon, A., Shimwell, T. W., Bonafede, A., et al. 2018, MNRAS, 478, 885
  * Botteon et al. (2020b) Botteon, A., van Weeren, R. J., Brunetti, G., et al. 2020b, MNRAS, 499, L11
  * Bregman & Irwin (2007) Bregman, J. N. & Irwin, J. A. 2007, ApJ, 666, 139
  * Brown et al. (2017) Brown, S., Vernstrom, T., Carretti, E., et al. 2017, MNRAS, 468, 4246
  * Brunetti et al. (2001) Brunetti, G., Setti, G., Feretti, L., & Giovannini, G. 2001, MNRAS, 320, 365
  * Brunetti & Vazza (2020) Brunetti, G. & Vazza, F. 2020, Phys. Rev. Lett., 124, 051101
  * Bryan et al. (2014) Bryan, G. L., Norman, M. L., O’Shea, B. W., et al. 2014, ApJS, 211, 19
  * Caprioli & Haggerty (2019) Caprioli, D. & Haggerty, C. 2019, in International Cosmic Ray Conference, Vol. 36, 36th International Cosmic Ray Conference (ICRC2019), 209
  * Cavaliere & Rephaeli (2011) Cavaliere, A. & Rephaeli, Y. 2011, Astrophysics of Galaxy Clusters, Vol. 172
  * Cen & Ostriker (1999) Cen, R. & Ostriker, J. P. 1999, ApJ, 514, 1
  * de Graaff et al. (2019) de Graaff, A., Cai, Y.-C., Heymans, C., & Peacock, J. A. 2019, A&A, 624, A48
  * Di Gennaro et al. (2018) Di Gennaro, G., van Weeren, R. J., Hoeft, M., et al. 2018, ApJ, 865, 24
  * Drury (1983) Drury, L. O. 1983, Reports on Progress in Physics, 46, 973
  * Eckert et al. (2015) Eckert, D., Jauzac, M., Shan, H., et al. 2015, Nature, 528, 105
  * Ensslin et al. (1998) Ensslin, T. A., Biermann, P. L., Klein, U., & Kohle, S. 1998, A&A, 332, 395
  * Feretti et al. (1997) Feretti, L., Boehringer, H., Giovannini, G., & Neumann, D. 1997, A&A, 317, 432
  * Forero-Romero et al. (2009) Forero-Romero, J. E., Hoffman, Y., Gottlöber, S., Klypin, A., & Yepes, G. 2009, MNRAS, 396, 1815
  * Gheller & Vazza (2020) Gheller, C. & Vazza, F. 2020, MNRAS, 494, 5603
  * Gheller et al. (2016) Gheller, C., Vazza, F., Brüggen, M., et al. 2016, MNRAS, 462, 448
  * Gheller et al. (2015) Gheller, C., Vazza, F., Favre, J., & Brüggen, M. 2015, MNRAS, 453, 1164
  * Ghosh et al. (2012) Ghosh, A., Prasad, J., Bharadwaj, S., Ali, S. S., & Chengalur, J. N. 2012, MNRAS, 426, 3295
  * Govoni et al. (2019) Govoni, F., Orrù, E., Bonafede, A., et al. 2019, Science, 364, 981
  * Hahn et al. (2007) Hahn, O., Porciani, C., Carollo, C. M., & Dekel, A. 2007, MNRAS, 375, 489
  * Hodgson et al. (2021) Hodgson, T., Vazza, F., Johnston-Hollitt, M., & McKinley, B. 2021, Publications of the Astronomical Society of Australia, 38, e047
  * Hoeft & Brüggen (2007) Hoeft, M. & Brüggen, M. 2007, MNRAS, 375, 77
  * Hoffman et al. (2012) Hoffman, Y., Metuki, O., Yepes, G., et al. 2012, MNRAS, 425, 2049
  * Jasche & Kitaura (2010) Jasche, J. & Kitaura, F. S. 2010, MNRAS, 407, 29
  * Jasche et al. (2010) Jasche, J., Kitaura, F. S., Wandelt, B. D., & Enßlin, T. A. 2010, MNRAS, 406, 60
  * Jasche & Lavaux (2019) Jasche, J. & Lavaux, G. 2019, A&A, 625, A64
  * Jasche et al. (2015) Jasche, J., Leclercq, F., & Wandelt, B. D. 2015, J. Cosmology Astropart. Phys., 2015, 036
  * Jasche & Wandelt (2013) Jasche, J. & Wandelt, B. D. 2013, MNRAS, 432, 894