1. Introduction
The term ‘cosmic web’ has been used to evoke the distribution of matter in our Universe on the very largest of scales. In this model, perturbations shortly after inflation have been amplified by gravitational instability to drive a process of hierarchical structure formation: low-density regions have evolved into giant voids, high-density regions have seeded galaxies, groups, and clusters, and these are connected by a network of low-density filaments and sheets (e.g. Baugh et al. Reference Baugh2004). Outside of galaxies and cluster cores, some 40–50% of the baryonic mass of the Universe is believed to trace this cosmic web structure, existing as a diffuse, highly-ionised, warm–hot intergalactic medium (e.g. Cen and Ostriker Reference Cen and Ostriker1999; Davé et al. Reference Davé2001).
Until recently, the existence of the WHIM and its distribution in a cosmic web was primarily inferred from simulations of Universe evolution. Tentative empirical results have begun to support this model (e.g. Eckert et al. Reference Eckert2015; Nicastro et al. Reference Nicastro, Kaastra, Krongold, Borgani, Branchini, Cen, Dadina, Danforth, Elvis, Fiore, Gupta, Mathur, Mayya, Paerels, Piro, Rosa-Gonzalez, Schaye, Shull, Torres-Zafra, Wijers and Zappacosta2018; de Graaff et al. Reference de Graaff, Cai, Heymans and Peacock2019; Tanimura et al. Reference Tanimura, Aghanim, Douspis, Beelen and Bonjean2019). Most recently, Macquart et al. (Reference Macquart2020) provided compelling evidence in support of this hypothesis by tracing the dispersion measure of a small collection of fast radio bursts, the origins of each having been traced to a known host galaxy. Despite the small sample size, this result has provided the strongest evidence yet for these missing baryons residing along the line of sight in the intracluster medium.
The synchrotron cosmic web is the expected radio component emitted by this large-scale structure (Brown Reference Brown2011). As part of the ongoing process of large-scale structure formation, simulations such as Vazza et al. (Reference Vazza, Ferrari, Brüggen, Bonafede, Gheller and Wang2015, 2019) have modelled large-scale accretion processes and shown them capable of producing shocks surrounding filaments and the outermost regions of galaxy clusters many times the local speed of sound, with Mach numbers as high as
${\mathcal{M}} \sim$
10–100. Such shocks are capable of producing high-energy electrons by way of diffusive shock acceleration (DSA) and, in the presence of large-scale intracluster magnetic fields, this suprathermal population will in turn produce synchrotron emission (e.g. Keshet et al. Reference Keshet, Katz, Spitkovsky and Waxman2009). The strength of this emission is predicted to be extremely weak, however, and in previous modelling, it was predicted to be at or below the level of detectability of the current generation of low-frequency radio instruments (Vazza et al. Reference Vazza, Ferrari, Brüggen, Bonafede, Gheller and Wang2015). Using more recent cosmological numerical simulations as guidance, Gheller and Vazza (Reference Gheller and Vazza2020) tested the detectability of a number observables that trace the gas in the cosmic web using a cross-correlation technique. They reported that the observables with the strongest cross-correlation with the underlying distribution of galaxies should be those involving the Sunyaev–Zeldovich effect and, in spite of the difficulty of its detection, the diffuse synchrotron radio emission from the shocked cosmic web.
Two separate papers by Brown et al. (Reference Brown2017) and Vernstrom et al. (Reference Vernstrom, Gaensler, Brown, Lenc and Norris2017) both attempted a statistical detection using a cross-correlation analysis. Vernstrom et al. (Reference Vernstrom, Gaensler, Brown, Lenc and Norris2017) used a deep, 180 MHz observation of a 21.8° × 21.8° field using the Murchison Widefield Array (MWA; Tingay et al. Reference Tingay2013) and cross-correlated this (residual) map with galaxy density maps at various redshifts, smoothed to scales ranging from 1 to 4 Mpc. Their expectation was that there would be a peak in the cross-correlation at 0° offset, and this expectation was rooted in an assumption that cosmic web emission broadly traced the large-scale mass distribution of the Universe. However, other radio sources such as active galactic nuclei (AGNs) and star-forming galaxies (SFGs) also correlate with galaxy density maps, making it necessary to accurately model these related emission populations to distinguish their signals. Indeed, despite the authors detecting a peak in the correlation at 0° offset, it was this confounding factor that prohibited any claim to a positive detection.
Brown et al. (Reference Brown2017) similarly attempted to use a cross-correlation analysis using the S-Band Polarization All-Sky Survey observed with Parkes at 2.3 GHz (Carretti et al. Reference Carretti, Staveley-Smith, Govoni, Feretti, Murgia, Giovannini, Vacca and Brown2013). However, rather than using galaxy density as a proxy for the cosmic web, they used cosmological simulations that were constrained to reproduce the local large-scale structure and which tracked the evolution of thermal gas and magnetic field strengths. The resulting synchrotron cosmic web emission S was modelled as a function of thermal electron density n e and magnetic field strength B in the form S ∝ n e B 2. From this, they produced a large-scale, low-resolution map of the local synchrotron cosmic web showing that it broadly and smoothly traced out the underlying mass density of the local Universe. The cross-correlation showed no statistically significant detection.
Both papers point to the great difficulty of making a detection of the synchrotron cosmic web. In particular, they point to the need for future detection attempts to be able to accurately model how the radio emission traces the underlying matter density, as well as to understand the much brighter population of AGN and SFG radio sources, how these cluster with respect to the underlying synchrotron cosmic web, and how their emission may produce confounding signals.
Most recently, Vernstrom et al. (Reference Vernstrom, Heald, Vazza, Galvin, West, Locatelli, Fornengo and Pinetti2021) have claimed the first definitive detection of the synchrotron cosmic web. By using luminous red galaxies as tracers of cluster cores, Vernstrom et al. stacked nearby cluster pairs found in low-frequency radio data produced both by the MWA as well as Owens Valley Radio Observatory Long Wavelength Array (OVRO-LWA Eastwood et al. Reference Eastwood2018) and identified a residual signal produced in the spanning intracluster medium. The authors argue against alternative explanations for the signal such as intervening cluster emission or overdense AGN and SFG emission spanning the filaments; whilst the results of this experiment are promising, the work also points to the need to accurately understand and model the combined extragalactic population alongside cosmic web emission so as to provide robust constraints on possible contamination within the stacking signal.
In response to the need for such simulations, we present the first sky model providing both the synchrotron cosmic web alongside a realistically clustered AGN and SFG radio population, the ‘FIlaments & GAlactic RadiO’ (FIGARO) simulation. We provide this model in the form of 10 4° × 4° light cones out to a redshift of z = 0.8, and valid for observing frequencies ranging from 100 to 1 400 MHz. To do this, we have combined cosmic web emission extracted from the magnetohydrodynamic (MHD) simulations by Vazza et al. (Reference Vazza, Ettori, Roncarelli, Angelinelli, Brüggen and Gheller2019) and have populated the light cone with AGN and SFG radio sources using the Tiered Radio Extragalactic Continuum Simulation (T-RECS) codebase produced by Bonaldi et al. (Reference Bonaldi, Bonato, Galluzzi, Harrison, Massardi, Kay, De Zotti and Brown2019). Uniquely, this latter population are positioned and clustered realistically with respect to the underlying mass density of the cosmological simulation. We expect this simulation to be important in developing observing and detection strategies for detecting the cosmic web with both current as well as upcoming low to mid-frequency radio telescopes.
This paper is structured as follows. We begin by discussing the construction of this simulation as well as the verification steps taken during that process. In Section 2, we discuss the underlying MHD simulation, the extraction of synchrotron cosmic web emission from the snapshots and the calibration of this emission against the small known population of radio relics. In Section 3, we discuss the extraction and validation of dark matter halos from the simulation, taking care to verify their number density and clustering properties as this population are a key input for the T-RECS simulation. Section 4 details the light cone construction process itself, stacking the cosmic web and dark matter halos out to a redshift depth of z = 0.8. Finally, in Section 5, we discuss how we use the halo light cones to position the AGN and SFG radio population using the T-RECS simulation codebase, before presenting the completed simulation catalogue in Section 6. We finish in Section 7 with a discussion oriented around the possibility of detection of the cosmic web with the latest generation of radio telescopes.
Throughout this paper, and including the original cosmological simulation, a ΛCDM cosmological model is assumed, with density parameters ΩBM = 0.0478 (baryonic matter), ΩDM = 0.2602 (dark matter), and ΩΛ = 0.692, and the Hubble constant H 0 = 67.8 km s−1 Mpc−1. All stated observing resolutions refer to the full width at half maximum (FWHM) of a circular Gaussian.
2. Cosmological simulation and snapshots
In beginning to build FIGARO, we must start with an underlying and evolving cosmological simulation. We have used the MHD simulation detailed in Vazza et al. (Reference Vazza, Ettori, Roncarelli, Angelinelli, Brüggen and Gheller2019), which was produced using enzo.Footnote a , Footnote b This simulation encompassed a comoving volume of 1003 Mpc3 with a uniform grid of 2 4003 cells and 2 4003 dark matter particles, each with a fixed mass set at 8.62 × 106 M⊙. This gives a spatial resolution of 41.63 kpc3 (comoving) per cell. The simulation was initialised at z = 45 with a simple uniform magnetic field of B 0 = 0.1 nG, and these fields were evolved in time using the MHD method of Dedner et al. (Reference Dedner, Kemm, Kröner, Munz, Schnitzer and Wesenberg2002).
Whilst larger, purely dark matter simulations do exist, a full MHD treatment is significantly more computationally expensive. The simulation used here is the largest of its kind to date, and presents a careful computational trade-off between cellular resolution and total volume. In the latter case, the volume is large enough to allow the simulation of a small population of ∼1014 M⊙ clusters and galaxy groups.
Figure 1. A comparison of the radio relic luminosity function (RRLF) from Nuza et al. (Reference Nuza, Hoeft, van Weeren, Gottlöber and Yepes2012) (dashed grey) with the the measured relics in each snapshot volume for HB07 (blue) and HB07 with additional fossil electrons (red). Relic statistics were calculated by summing emission in the annulus around dark matter halos with radii
$0.5 \cdot r_{200} < r < 1.5 \cdot r_{200} $
.
An important caveat to note is that the simulation does not include either radiative gas cooling or feedback processes due to star formation and AGN, both of which are important to the evolution of cluster interiors, nor was the spatial resolution sufficiently high to capture turbulent dynamo effects of the dense cluster regions that can significantly magnify magnetic fields. Our regions of interest for this catalogue, however, are the filaments and cluster peripheries and these regions are sufficiently distant from core cluster environs that these effects, at least to first order, are not especially relevant.
2.1. Radio
Synchrotron radio emission is produced by electrons at relativistic energies interacting with background magnetic fields. In this simulation we trace synchrotron emission solely as a result of diffusive shock acceleration (DSA) typically associated with accretion shocks, cluster mergers, and other large-scale structure formation processes (e.g. Ryu et al. Reference Ryu, Kang, Hallman and Jones2003).
In DSA, a small fraction of ambient electrons are accelerated to relativistic energies, which then radiate due to their interaction with non-neglible intracluster magnetic fields. We model the resulting radio power based on Hoeft and Brüggen (Reference Hoeft and Brüggen2007), henceforth HB07:
\begin{equation} P_\nu \propto S \cdot n_d \cdot \xi\left(M,T\right) \cdot \nu^{\alpha} \cdot T_d^{{3}/{2}} \cdot \frac{B^{1 - \alpha}}{B^2_{{\textrm{CMB}}} + B^2},\end{equation}
where S is the shock surface area; n d is the downstream electron density; ξ(M, T) is the electron acceleration efficiency which is a function of Mach number and temperature; ν is the frequency and α is the spectral index of the radio emission; T d is the downstream electron temperature; B is the magnetic field in each cell; and B CMB is the equivalent magnetic field strength of the CMB where B CMB ≈ 3.25(1 + z)2 μG.
Both B and ξ(M, T) are poorly understood, especially in the highly diffuse filaments and cluster outskirts, yet are key parameters in calculating the radio power. The magnetic field strength along filaments is primarily the result of adiabatic gas compression, and in these low-density regions, the field strength is closely related to the magnetic field seeding scenario (e.g. Vazza et al. Reference Vazza, Ferrari, Brüggen, Bonafede, Gheller and Wang2015). The seed strength, however, is unknown and thus predictions of these magnetic fields vary between cosmological simulations in the range of 10−4−0.1 nG. At best, we have upper limits provided by Planck from CMB observations which limit the magnetic field strengths to a few nG on scales of 1 Mpc (Planck Collaboration 2016). In this simulation, the magnetic field was seeded uniformly at 0.1 nG—an order of magnitude lower than these limits—but there is latitude in the choice of this parameter. Moreover, for values of B ≪ B CMB, the radio power scales P ∝ B 2, making the simulation particularly sensitive to the seed strength.
The electron acceleration efficiency ξ(M, T) is a second, crucial parameter that is difficult to model. This parameter estimates the fraction of thermal electrons that are accelerated to relativistic energies as a result of a shock. The model used here is based on HB07 which depends upon the strength of the shock and the thermal temperature of the downstream electrons. However, this model does not exceed ∼10−3 for reasonable Mach values, and this value is insufficient to account for many observed radio relics (e.g. Botteon et al. Reference Botteon, Brunetti, Ryu and Roh2020). For these events, we now believe ‘fossil’ electron populations—those which have been previously accelerated by, for example, AGN activity or previous DSA shocks—allow for a much higher effective electron acceleration efficiency (e.g. Pinzke, Oh, & Pfrommer Reference Pinzke, Oh and Pfrommer2013). As shocks were here calculated in post-processing, this simulation does not have a ‘memory’ of previous shocks nor does it model AGN activity; at each snapshot, the electron population is therefore always assumed to be at thermal equilibrium resulting in underestimated radio emission, especially in dense cluster environments.
As a mitigation, we introduce a modified HB07 model (herein: ‘HB07 + fossil’) for the acceleration efficiency. This mitigation builds upon HB07 with a special case for weak shocks in dense environments: for shocks with
${\mathcal{M}} < 5$
and thermal temperature T
$ > $
107K, we arbitrarily set ξ(M, T) = 10−2. The final HB07+fossil model was then a weighted sum of the original HB07 model with this special case, with weightings 0.95 and 0.05, respectively.
These weightings were chosen to best reproduce the radio relic luminosity function (RRLF) derived by Nuza et al. (Reference Nuza, Hoeft, van Weeren, Gottlöber and Yepes2012). Radio relics are a known class of radio source produced by DSA at the periphery of clusters and is driven by shocks from large-scale structure formation processes. Radio relics should be well modelled by our simulation and provide a means by which to calibrate the radio emission of our models. In Figure 1, we compare the relic counts for both the HB07 (blue) and HB07+fossil (red) models for each of the snapshot volumes of our simulationFootnote c with the Nuza et al. (Reference Nuza, Hoeft, van Weeren, Gottlöber and Yepes2012) RRLF in dashed grey.Footnote d We can observe that the standard HB07 model results in significantly under-powered relics by about one order of magnitude, whilst the modified HB07 + fossil weighted model brings the relic counts into good agreement.
2.2. Snapshots
We have extracted four ‘snapshots’ of the simulation at redshifts ranging from z = 0.025 to z = 0.6, as detailed in Table 1, and these snapshots have formed the basis for constructing our light cones. This selection of snapshots was driven purely by the data still available. Amongst the observables produced by the simulation, the three of note to this work are: dark matter density cubes; the calculated luminosity of the synchrotron emission at 1 400 MHz; and the shock values (i.e. cell Mach numbers).
Table 1. The simulation snapshots used to construct the light cone and associated halo catalogues.
3. Halo finding
Alongside the radio emission, the mass distribution, and in particular the dark matter halo positions, are a key ingredient in building the FIGARO simulation. This arises because T-RECS relies on dark matter halos to position AGN and SFG and thus accurately mimic the clustering properties of these sources with respect to the cosmic web emission. Thus for each snapshot, we needed to identify dark matter halos in post-processing using each snapshot’s associated dark matter cube. The processed outputs of the simulation have been smoothed using a Cloud-In-Cell (CIC) kernel, prohibiting traditional ‘friend of friend’ algorithms. Instead, to identify dark matter halos we used the following simple algorithm: a sphere was progressively grown around the most massive voxel until the mean density of the enclosed volume reduced beneath a threshold density
$\bar{\rho}$
, in this case
$\bar{\rho} = 200 \rho_c$
, where ρ
c
is the critical density of the Universe at that respective redshift. From this, we could then interpolate the virial radius r
200 and mass M
200 of the halo. The mass of enclosed voxels was set to zero and the process was repeated on the next most massive voxel in the snapshot, and so on, until exhausting the volume of resolved halos.
To provide enough potential source positions, T-RECS requires that we find halos down to ∼109.5 M⊙. Given the relatively coarse spatial and mass resolution of our simulation, detecting these low-mass halos is challenging. To achieve this we set the minimum allowable halo size as having a radius of 1 voxel, that is, incorporating just 7 voxels in total. At this limit, however, we are especially wary of the introduction of errors or bias due to discretisation effects, and we have therefore performed a number of sanity checks upon our halo catalogue.
In Figure 2 we show the halo counts, binned by mass, for each respective snapshot. We compare these counts with the halo mass function by Angulo et al. (Reference Angulo, Springel, White, Jenkins, Baugh and Frenk2012) (produced using HMFCalc by Murray et al. Reference Murray, Power and Robotham2013) and we observe good agreement down to about 1010.5 M⊙, below which our catalogue becomes increasingly incomplete. The lowest mass halos range from 109.87 M⊙ at z = 0.025 to 109.57 M⊙ at z = 0.6, where this lower limit steadily decreases due to the decreasing cosmological scale factor and corresponding increased resolution of the simulation at earlier epochs.
Figure 2. Halo counts for each snapshot, compared with the halo mass function by Angulo et al. (Reference Angulo, Springel, White, Jenkins, Baugh and Frenk2012) (produced using HMFCalc by Murray et al. Reference Murray, Power and Robotham2013) at the respective redshifts. The halos are binned by mass into intervals of 0.1 dex.
As a second sanity check, in Figure 3, we compare the clustering properties of the FIGARO halos with those from the Millennium II (MII) and Planck Millennium (PM) simulations by calculating the two-point correlation functions of each respective halo population (for details, see Appendix A). Both the MII and PM simulations are of sufficiently high resolution so as to be able to properly resolve low-mass halos on the order of 109.5 M⊙ and therefore provide a good benchmark for comparison. What we observe is good agreement between the simulations across all mass ranges for scales greater than 0.3 Mpc. For scales less than this, both MII and PM turnover slightly, whereas FIGARO continues its power law behaviour. This discrepancy exists for all mass bins and so is not directly related to the discretisation of the simulation, although its cause is unknown.
Figure 3. The two-point correlation of halos for the FIGARO (blue), Millennium II (MII; red), and Planck Millennium (PM; green) simulations. The results are shown at redshifts z ≈ 0.025 and z ≈ 0.309 and are binned by mass in 0.5 log10(M⊙) increments. Additionally, an all mass result is also shown. In all, the halo clustering properties of FIGARO are consistent with both MII and PM simulations with the exception of scales under ∼0.3 Mpc.
These sanity checks satisfy to us that our halo catalogue is sufficiently robust, despite the coarse resolution of our simulation.
4. Light cone construction
With both the synchrotron cosmic web and the dark matter halos catalogued for each of our underlying snapshot volumes, we can now proceed to tile these volumes so as to construct the FIGARO light cones. This cone was constructed spanning a 4° × 4° angular extent and with a depth extending to a redshift of z = 0.8. The simulation volume, with sides of length 100 Mpc (comoving) only extends out to z = 0.023, and the field of view exceeds a comoving transverse distance of 100 Mpc at z = 0.35. This necessitates replicating the simulation volume along the length of the light cone a total of 34 times, as well as laterally tiling the volume beyond z = 0.35. We discuss here the construction of this light cone.
The web and halo light cones were constructed identically. In both cases, we began with a catalogue of voxel positions and properties. For the web, these properties were the voxel luminosity and shock value; for the halos these were the halo mass and radius. Each cone was constructed by appending the snapshot catalogue in 100 Mpc increments along the length of the light cone. To mitigate repeating structures along the line of sight, we introduced a random transverse offset on each iteration as well as sequentially rotating the volume through each of its three axes. Then, for each entry in the snapshot catalogue, we calculated its redshift, latitude and longitude and, if these latter values fit within the 4° × 4° field of view, they would be appended to the final catalogue. Note that throughout, we assume the field of view is small enough that we can safely assume a flat screen projection.
As noted, at redshifts of about z = 0.35 it becomes necessary to tile the simulation volume so as to fully cover the field of view. Since the simulation volume wraps at its edges, this tiling does not cause discontinuities or edges in the final light cone. The volume was tiled in a 2 × 2 arrangement from z = 0.35.
The random transverse offset allows us to additionally produce multiple light cones or ‘realisations’ which, especially at low redshifts, can be significantly different to each other. In one realisation, for example, the low redshift cone may be largely empty whilst in another, just by chance, there may be a massive galaxy cluster within the field of view. This reflects the kind of cosmic variance we should expect. We have produced 10 realisations in total.
In the case of the halos, their properties were simply appended to a catalogue. In the case of the web emission, however, some additional processing was required. For each web voxel, we calculated its flux S ν from the simulation luminosity values L ν using the standard radio luminosity function with k-correction:
\begin{equation} S_{\nu}\!\left(z\right) = L_{\nu} \frac{{(1 + z)}^{\alpha}}{4 \pi D_L^2},\end{equation}
where D
L
is the luminosity distance at the redshift z of the voxel and α is the voxel’s spectral index. In turn, the spectral index α was calculated from the voxel’s associated Mach number
${\mathcal{M}}$
using the relation:
\begin{equation} \alpha = \frac{{\mathcal{M}}^2}{1 - {\mathcal{M}}^2},\end{equation}
where the power law relation uses the positive sign convention S ∝ ν α . This latter relation is also used to scale the luminosity values from the 1 400 MHz output of the simulation to the observing frequency. Finally, this flux value was appended to a map spanning the field of view with 3′′ × 3′′ resolution.
The relatively small volume of the simulation raises issues as we go deeper in redshift space. Whilst the simulation volume is large enough to give a good representative sample of the Universe for small and medium mass clusters, it is, however, just large enough to simulate a single massive galaxy cluster, on the order of 1015 M⊙. As we have shown in Table 2, at low redshifts, we sample only a very small fraction of the total simulation volume; since this sample will primarily be low or medium mass clusters, it is a reasonable representation of the Universe as whole. Moreover, at these low redshifts, the differences between realisations should give an accurate picture of the effect of cosmic variance. However, as we go deeper in redshift and as the fractional volume increases towards unity, the different realisations will increasingly sample from the same parts of the simulation, thus becoming self-similar. Moreover, the particular statistics associated with the small population of high mass clusters in our volume, which also host the most luminous cosmic web emission, will begin to dominate statistics derived from these redshift slices. Finally, at redshifts of around z = 0.35, where tiling becomes necessary, the most massive galaxy cluster in the volume will be replicated in every single 100 Mpc increment. This single cluster and its unique evolution will dominate most statistics derived from these deeper redshift slices. This is a limitation we will need to consider when we later discuss the statistics of cosmic web emission.
Table 2. The comoving volume enclosed by regular Δz = 0.05 redshift slices and the 4° × 4° field of view. For higher redshifts
$(z \gtrapprox 0.15)$
, the volume as a fraction of the total simulation volume is sufficiently large enough that different realisations will be increasingly similar. For redshift slices greater than approximately z = 0.25, it becomes necessary to duplicate the simulation volume more than once.
5. T-RECS
T-RECS (Bonaldi et al. Reference Bonaldi, Bonato, Galluzzi, Harrison, Massardi, Kay, De Zotti and Brown2019) is a simulation of the continuum radio sky from 150 MHz to 20 GHz that models radio emission from AGN and SFG using a number of empirically derived relations, and forms the final ingredient required for FIGARO. The simulation provides a thorough range of properties for each radio source, ranging from general properties of intrinsic luminosity, redshift, and physical size to more specific properties such as AGN classification, radio jet angle, and SFG ellipticity. Each of the two general populations were further broken down into subpopulations. The AGNs consisted of: steep-spectrum sources (SS-AGNS), flat-spectrum radio quasars (FSRQ), and BL Lac. The SFGs were also further subdivided into three subpopulations: late-type, spheroidal, and lensed spheroidal galaxies.
The T-RECS simulation was constructed with careful attention paid to the spatial clustering of radio sources. This clustering was implemented by associating each AGN and SFG radio source with a dark matter halo extracted from a 5° × 5° light cone that Bonaldi et al. (Reference Bonaldi, Bonato, Galluzzi, Harrison, Massardi, Kay, De Zotti and Brown2019) had constructed, originally derived from the PM simulation. Bonaldi et al. (Reference Bonaldi, Bonato, Galluzzi, Harrison, Massardi, Kay, De Zotti and Brown2019) describe in detail the way in which each radio population was associated with dark matter halos of a particular mass range, which we briefly summarise here. In the first case, the AGN subpopulations were associated with dark matter halos by first relating a stellar mass M * to each halo mass M h (i.e. M * = F(M h )), and then in turn calculating the fraction of galaxies hosting an AGN (radio loud or radio quiet) as a function of that host galaxy mass. The SFGs, on the other hand, used an abundance matching process to relate the known distribution of halo masses and the known distribution of SFG luminosities. This allowed for associating the most luminous SFGs with dark matter halos, whilst for SFGs where the luminosity implied a dark matter halo mass smaller than allowed for by the resolution of the PM simulation, a random distribution on the sky was instead assumed.
Crucially for our purposes, Bonaldi et al. (Reference Bonaldi, Bonato, Galluzzi, Harrison, Massardi, Kay, De Zotti and Brown2019) also made the simulation code publicly available. The has allowed us to run the simulation ourselves, and in the process input our own dark matter halo catalogue for each of the respective realisation light cones. Thus the output extragalactic catalogue positions AGN and SFG sources with respect to the underlying matter density of our lightcones, and in this way we accurately model how the cosmic web and the embedded radio population cluster with respect to each other.
We made use of the T-RECS codebase largely without modification, with the exception of a bug fix that corrected a random number generation routine that incorrectly produced strongly biased results. We provided the halo catalogues of each light cone realisation in the format T-RECS expected and ran the simulation for a 4° × 4° field of view, with a lower flux threshold of S 1 400MHz = 0.1μJy, and for a redshift range that encompassed the full cosmic web light cone. The choice of lower flux threshold was selected so as to fully simulate classical confusion noise in the frequency range 150–1 400 MHz assuming a highest resolution of 8′′ (see Appendix B).
Finally, we also note the ease with which future improvements to T-RECS modelling of AGN and SFG populations can readily be regenerated into this simulation, especially as deeper, large-field surveys provide more accurate statistics of these faint sources (e.g. Prandoni 2021).
5.1. Clustering of the radio population
To validate the spatial distribution of the resulting radio catalogue, we turn to the the angular two-point correlation function. This function ω(θ) describes the apparent two-dimensional clustering of radio sources based on their angular separation, without reference to their redshift. Here, we compare the results of two recent empirical measurements of this value, by Magliocchetti et al. (Reference Magliocchetti, Popesso, Brusa, Salvato, Laigle, McCracken and Ilbert2017) and Hale et al. (Reference Hale, Jarvis, Delvecchio, Hatfield, Novak, Smolčić and Zamorani2018), against that of our simulation.
Magliocchetti et al. (Reference Magliocchetti, Popesso, Brusa, Salvato, Laigle, McCracken and Ilbert2017) reported a measurement of the angular two-point correlation function based on a catalogue of sources derived from observations at 1.4 GHz of a 2 deg2 region of the COSMOS field. By setting a lower flux threshold of 0.15 mJy, above which value the catalogue was considered complete, and assuming ω(θ) to be a power law of the form ω(θ) = Aθ 1−γ where γ was set to 2, they were thereby able to derive a value for the constant of proportionality as A = 2.2 ± 0.4 × 10−3. Hale et al. (Reference Hale, Jarvis, Delvecchio, Hatfield, Novak, Smolčić and Zamorani2018) similarly used observations of a 2 deg2 region of the COSMOS field, but at the higher frequency range of 2–4GHz and using a much lower flux threshold 12.65μJy (5.5 times the median image noise). They assumed the same power law form but instead fixed γ = 1.8 and thereby derived a value for the constant of proportionality as log 10 A = −2.8 ± 0.1.
Figure 4. The angular two-point correlation of all radio sources (AGN, SFG, and all subtypes) for FIGARO (realisation 1; blue), the original TRECS (red), and other FIGARO realisations (grey crosses), calculated across the redshift range
$ 0 < z < 0.8$
. Left: Comparison to Magliocchetti et al. (Reference Magliocchetti, Popesso, Brusa, Salvato, Laigle, McCracken and Ilbert2017) using their minimum flux threshold of 0.15 mJy at 1 400 MHz. Right: Comparison to Hale et al. (Reference Hale, Jarvis, Delvecchio, Hatfield, Novak, Smolčić and Zamorani2018) using their minimum flux threshold of 12.65μJy at 3 GHz.
In Figure 4, we compare these measurements against those of our simulation. We estimate the angular version of the two-point correlation function using the equation described in Appendix A with the exception that the Euclidean metric r is replaced with the apparent angular separation θ between sources. All error estimates are calculated using bootstrap sampling with 100 iterations. The function was calculated for angular separations 10−3.25
$ < $
z
$ < $
10−0.25 in equal logarithmically spaced bins of width log 10Δθ = 0.5.
In the left panel, we compare the measurement of Magliocchetti et al. (Reference Magliocchetti, Popesso, Brusa, Salvato, Laigle, McCracken and Ilbert2017) (shaded blue region) with those of realisation 1 from FIGARO (blue, error bars indicate three standard deviations). We also calculate ω(θ) for the other realisations (grey crosses), which indicate the variation in this function purely as a result of cosmological variance; however, for clarity we have not included their associated error bars. For reference, we also plot ω(θ) of T-RECS using its PM cosmology for the matching redshift range 0
$ < $
z
$ < $
0.8.Footnote
e
The simulated catalogues have had a flux threshold of S
1.4GHz
$ > $
0.15 mJy applied. On the right, we compare the measurement of Hale et al. (Reference Hale, Jarvis, Delvecchio, Hatfield, Novak, Smolčić and Zamorani2018) with the same set of simulations except subject to a flux threshold of S
3GHz
$ > $
12.65μJy.
In both cases, FIGARO is reasonably consistent with these empirical measurements for
$\theta \gtrapprox 10^{-1.75}$
, especially when we take into account the spread of values measured from different realisations, as well as with the original T-RECS. FIGARO (and the original T-RECS), however, appears to be under-clustered on angular scales smaller than this. One partial explanation for this is a result of our coarse simulation volume. The resolution of our simulation (and therefore the halo catalogue) was 41.7 kpc, which at a redshift of z = 0.1 corresponds to a minimum angular resolution of 19.9′′ and even at a redshift of z = 0.3 the minimum angular resolution of 7.0′′ is still greater than the smallest bin for which we have calculated ω(θ). As a result, for a given redshift slice, we should expect ω(θ) will decline for values of θ smaller than the respective minimum angular scale, and this will introduce an under-clustering bias to our results on these small angular scales. The second contributing factor is the absence of higher redshift (z
$ > $
0.8) sources, although a comparison of the two-point statistics for the full T-RECS catalogue against the redshift limited catalogue shows this effect to be small. Nonetheless, we consider our results to be at least as good as the original T-RECS catalogue and sufficient for our present purposes.
6. Results
Combining the final cosmic web light cones and T-RECS catalogues, we produce 10 4° × 4° realisations of the radio sky encompassing the redshift range 0
$ < $
z
$ < $
0.8. We additionally provide a background catalogue consisting only of T-RECS sources in the redshift range 0.8
$ < $
z
$ < $
8; these sources were not clustered but were instead positioned randomly across the field of view using a uniform distribution, which was uniquely generated for each realisation. Figure 5 is illustrative of the final cosmic web catalogue, where in the top image, we present realisation 5 encompassing the full simulated redshift range of cosmic web emission (0
$ < $
z
$ < $
0.8) and T-RECS sources (0
$ < $
z
$ < $
8) at 900 MHz, and convolved to a beam resolution of 20′′. The density fluctuations of the extragalactic radio populations are subtly visible, and using this colour scale we can also identify a small handful of peaks in the cosmic web emission. In the bottom panel, we present the cosmic web emission only, presented using a logarithmic colour scale.
Figure 5. A full-field (4° × 4°) image of realisation 5 at 900 MHz encompassing the redshift range
$0 < z < 0.8 $
for the cosmic web and
$ 0 < z < 8$
for T-RECS sources, and convolved to a beam resolution of 20′′. Top: The combined simulation with T-RECS extragalactic sources and faint, background cosmic web emission. The colour scale ranges from 0 to saturate at 200μJy. Bottom: Cosmic web emission only, using a logarithmic colour scale.
Accompanying this paper, we also make available a full data release of the simulation for each realisation, which includes catalogues describing their respective halos, cosmic web maps, and AGN and SFG populations along the length of the redshift cone. For each realisation, we provide the full data in Hierarchical Data Format 5 (HDF5), with five named datasets. These dataset names and their corresponding schemas are: halos (Table 3), web (Table 4), sfg (Table 5), and agn (Table 6).
Table 3. Halo dataset schema which describes the properties of dark matter halos along the length of a realisation’s redshift cone.
Table 4. Cosmic web array with four dimensions. The values for x_coord and y_coord denote the centre of the pixel, with each pixel having a value of Jy and occupying an area ∼3′′ × 3′′. The formulas above assume zero indexing (i.e. i ∈ [0, 1, . . .]).
Table 5. SFG dataset schema which describes the properties of SFG radio sources along the length of a realisation’s redshift cone.
Table 6. AGN dataset schema which describes the properties of AGN radio sources along the length of a realisation’s redshift cone.
Figure 6. The cumulative radio luminosity contained within the spherical volume surrounding each of the 100 most massive dark matter halos as a function of radius, for the snapshot volume at z = 0.025. The luminosity has been normalised to the value at r = 2.5 · r 200 for each curve. The median value across all halos is indicated in blue.
The cosmic web radio emission itself is provided as a four-dimensional array, with the first axis spanning integrated redshift slices in 0.05 increments, the next two axes spanning the full 4° × 4° field of view, and the final two a tuple containing the flux value and spectral index. The field of view is approximated as a flat screen, with each pixel coordinate offset by 3′′; with this approximation, each pixel occupies approximately 3′′ × 3′′ with a maximum error at the edge of 0.06%.
6.1. Morphology and the ubiquity of relics
In Figure 5, we showed the full 4° × 4° image of a typical realisation, with the lower figure showing just the cosmic web emission with logarithmic scaling. The phrase ‘cosmic web’ has been evoked to describe the large-scale structure of the Universe as a kind of sponge, and the synchrotron component has often been expected to follow the same hierarchy of structures with sheets, filaments and dense clusters. However, this is not what we observe. Instead, we observe ubiquitous relic-like, shocked shells of emission that surround dark matter halos. Indeed, these are the primary source of emission in the volume. These shocked shells envelop dark matter halos and in the lowest redshift snapshot, for example, ∼96% of the total power in the volume is located in the spherical shell (r
$ < $
1.5 · r
200) of the 100 most massive dark matter halos. The connecting filaments or inflows are not easily discernible, even in a logarithmically scaled image.
In Figure 6, we show the distribution of power for the z = 0.025 snapshot by plotting the cumulative radio luminosity contained within the spherical volume surrounding each of the 100 most massive dark matter halos, as a function of radius. The cumulative luminosity is plotted in grey for each of the halos, and in blue we plot the mean cumulative luminosity. The universal absence of emission at the core reminds us that in this simulation we do not model Fermi II shock processes such as radio halos. Instead, in general, we observe DSA processes occurring and most luminious in the range 0.75 · r
200
$ < $
r
$ < $
1.5 · r
200, although with significant spread in this range from as low as 0.25 · r
200 out to greater than 2 · r
200. The shells themselves are far from isotropic, having complex, filamentary structures with knots and shock fronts orders of magnitude more luminous than elsewhere in the shell.
The absence of significant DSA processes in the innermost cores of dark matter halos is primarily due to the increased matter density and corresponding increase to the speed of sound, as well as the proportionally small area of the shock front; shocks in these regions are therefore not effective electron accelerators despite these regions containing the highest density electron population and magnetic field strengths (Vazza et al. Reference Vazza, Brüggen, van Weeren, Bonafede, Dolag and Brunetti2012). As the shock proceeds significantly far away from the cluster core it loses energy, encounters an increasingly sparse electron environment, and magnetic field strengths decline; thus similarly these outermost environments are ineffective at producing synchrotron emission. In between, however, there exists a sweet spot which maximises the synchrotron output, precisely as we observe in the simulation.
Moreover, the majority of radio power within the volume surrounds just a handful of interacting clusters. In Figure 7, we plot the power contained within the spheres centred on dark matter halos with radius r = 1.5 · r 200 as a function of halo mass. We observe a power law correlation between mass and the cumulative radio power (L ∝ M 3.6; dashed grey line), although we observe significant scattering around this trend that shows the dependence on the specific interaction and merger histories of individual clusters. The plot also makes it clear that just handful of dark matter halo environments account for the majority of the power output in the volume: 90% of total power is located within the spheres of radius r = 1.5 · r 200 surrounding just 12 dark matter halos, which cumulatively account for 0.45% of the total volume.
Figure 7. Cosmic web power across the snapshot volume z = 0.025 within the dark matter halo spheres
$ r < 1.5 \cdot r_{200}$
as a function of dark matter halo mass. There is a power law trend indicated by dashed line (L ∝ M
3.6), but significant scattering occurs primarily as a result of the interaction and merger histories of specific clusters.
6.2. A brief survey of the brightest, most detectable features
To further illustrate the typical emission morphologies, in Figure 8, we provide a small sample of emission regions drawn from redshift slices 0.10
$ < $
z
$ < $
0.15 and 0.15
$ < $
z
$ < $
0.20 of realisations 1 and 2, observed at 900 MHz and with a 20′′ resolution. These sources were chosen as they are amongst the brightest emission and provided a good range of morphologies that are present throughout the catalogue. In contrast to previous figures, these are presented using a linearly-scaled colour map, ranging from 0 to the 99.5th percentile pixel, to make clear the morphology and angular size of the emission structures as would actually be observed. This is to correct for the impression of broad, diffuse emission that may be taken from logarithmic colour scales, when in fact the emission landscape we observe is one with bright ‘knotty’ peaks of emission and otherwise extremely faint surrounding islands of emission. The contours, however, are logarithmically scaled, indicating 0.1, 0.01, and 0.001 of the 99.5th percentile value. The r
200 radius of nearby dark matter halos are indicated by red dashed circles.
Figure 8. A sample of emission features at 900 MHz from from various realisations showing both familiar relic formations as well as more unusual shock morphologies. Each image has been convolved to a resolution of 20′′ (beam size indicated by white circle in top right), and the colour map (Jy beam−1) is varyingly scaled from 0 Jy beam−1 to the value of the 99.5th percentile pixel. Red circles indicate the r
200 dark matter halo extent for halos with M
200
$ > $
1012 M⊙. Note that all emission occurs outside the core region of dark matter halos, and only appears to be centrally located due to projection effects.
The most recognisable features in this collection are the very traditional double relic morphologies, which can be seen in Figure 8(a), (b), (c), and (d). These examples clearly show matching pairs of arced emission structures as shocked waves travel outwards from the dark matter halo centre, with very classic bow wave morphology. Figure 8(e) is more complex but still recognisable as a relic system, whilst 8f is asymmetric and appears as single relic system.
The remaining examples demonstrate the effect of projection angle on what are otherwise similarly structured shells of emission about dark matter halos. Figure 8(g), (h), and (i), for example, appear in projection to have centrally located emission. Figure 8(j) has two apparently point-like peaks of emission that are amongst the brightest features in this collection, with the brightest peaking at 8 μJy beam−1; these ‘points’ are aided in their apparent brightness by summing long, narrow emission structures in the radial direction.
Figure 8(k) and (l) are the most complex in this set of examples, and both involve numerous interacting dark matter halo systems. They are also the most luminous, presumably aided by these cluster interactions. The two Southern peaks in Figure 8(k) peak at approximately 12 μJy beam−1, whilst the central peak of emission in Figure 8(l) measures 52 μJy beam−1.
Crucially, note the absence of broad emission features, connecting filamentary bridges between dark matter halos, and in general, emission that isn’t associated with shocked shells surrounding dark matter halos.
6.3. The embedded radio population
We now consider the addition of the T-RECS sources. In Figure 9, we show a 50′ × 50′ field of view showing the redshift range 0.15
$ < $
= z
$ < $
0.2 extracted from realisation 3 and centred on a massive cluster system. This redshift range has a comoving depth of 203 Mpc, and so this extraction incorporates the full simulation volume stacked twice in the radial direction. The panels allow us to compare the distribution of the synchrotron cosmic web in comparison to the underlying mass distribution as well as the embedded radio population.
Figure 9. A 50 × 50′ field of view showing redshift range
$0.15 \leq z < 0.2 $
extracted from realisation 3. This redshift range has a comoving depth of 203 Mpc, and so this extraction incorporates the full simulation volume stacked twice in the radial direction. The 5 Mpc scale has been calculated at the mean redshift z = 0.175. Left: The synchrotron cosmic web emission at 900 MHz with resolution 20′′. Middle: The associated mass distribution with halos of mass
$M > 10^{12.5} $
M⊙ indicated by dashed red circles of radii r
200. Right: The combined cosmic web emission and TRECS radio population for this redshift range at 900 MHz with resolution 20′′. The TRECS radio population are modelled as simple point sources.
The left panel shows the cosmic web emission component of this redshift range observed at 900 MHz and at a resolution of 20′′. The colour scale in the image spans seven orders of magnitiude from the plausibly detectable 10 μJy beam−1 at the heart of the most massive dark matter halo in the field of view, to 10 nJy beam−1 on the virial periphery, to under 1 nJy beam−1 along the filaments. Around the most massive cluster, we observe a number of shocked shells of radio emission; these surround the cluster core, however, projection effects mean the brightest emission appears to align with the core region of the cluster.
The central panel shows the mass distribution within the field of view, where we have also overlaid the dark matter halos detected by our halo finding algorithm with mass M
$ > $
1012.5 M⊙. We can observe that this is in fact a merging system of two massive systems at the centre of the field, and it is clear that a bridge of increased mass density extends between the two systems. As we observed prior, the northern, most massive system is enshrouded by shocked shells of cosmic web emission, whilst the southern is comparatively quiet in this regard, and there is no significant emission associated with the bridge.
In the rightmost panel of Figure 9 we show the associated T-RECS radio populations, with both AGN and SFG modelled as simple point sources. The T-RECS sources are not uniformly distributed and display clustering coincident with the most massive portions of the field of view, as well as less dense regions coincident with the dark matter voids. A number of bright T-RECS sources overlap with the peak of the cosmic web emission. Whilst we have modelled the T-RECS sources as simple point sources, at this redshift and at a resolution of 20′′ the lobes of radio loud AGN can be resolved. In this sense, the rightmost panel represents a best-case view of the cosmic web, whereas a more sophisticated representation of the T-RECS catalogue would likely occlude the underlying cosmic web emission much more significantly.
7. Discussion
We orient the discussion of FIGARO primarily around practical questions of the detectability of the cosmic web with current and future radio instruments. This discussion will include its apparent flux distribution, characteristic angular scales, its correlation with the much brighter AGN and SFG populations and, ultimately, the possibility of its detection with the cross-correlation method discussed in the introduction.
To make these comparisons, we will make reference to three idealised instruments that map approximately to: the MWA in its phase 2 configuration (Wayth et al. Reference Wayth2018); the Australian Square Kilometre Array Pathfinder (ASKAP; Hotan et al. Reference Hotan2021); and the proposed Square Kilometre Array Low (SKA Low).Footnote f For our purposes, we characterise these instruments simply by an observing frequency, resolution, and noise limit, as shown in Table 7.
7.1. Flux comparison
In this section, we quantify flux differences between the embedded radio population and the synchrotron cosmic web. We begin by examining the flux sum across various redshift slices. In Table 8, we show the flux sum at 150 MHz for both of the embedded radio population and the cosmic web emission, computed across all 10 realisations. Across the depth of the redshift cone out to z = 0.8, the total flux attributed to AGN and SFG sources is 3.95 Jy deg−2, whilst the cosmic web emission over this depth is almost a factor of 150 lower, at 0.027 Jy deg−2. The full flux of the AGN and SFG sources out to z = 8 is 8.95Jy, for which we have no commensurate cosmic web flux values. If we bin the length of the light cones in Δz = 0.05 slices, we see that the cosmic web emission peaks in the two redshift slices 0.05
$ < $
z
$ < $
0.1 and 0.1
$ < $
z
$ < $
0.15, each providing a fractional signal of about 0.09% of the total flux of the simulation.
Table 7. Idealised observing configurations that approximately map to the MWA, SKA Low, and ASKAP radio interferometers. The resolution refers to the FWHM of a circular Gaussian beam.
Table 8. Flux statistics across all 10 realisations at 150 MHz. For average flux density values, we also provide the flux-weighted mean spectral index (as the exponent) allowing for extrapolation up to 1 400 MHz. Whilst the flux sums only depend on frequency, the final two columns, the 100th and 99.9th percentile values, are calculated with respect to the idealised MWA configuration.
Table 8 also shows the flux-weighted spectral index, which appears as the exponent value, and which allows for extrapolation of the flux sum up to 1 400 MHz. These spectral index values do not appear to be redshift dependent across the nearby redshift range we have considered, and show a consistent α ≈ −0.8 for extragalactic sources and α ≈ −1.25 for cosmic web emission.Footnote g Thus, whilst at 150 MHz, the total cosmic web emission is a factor of ∼150 times fainter than the embedded extragalactic emission, at 900 MHz this ratio increases to ∼300. All things being equal, the cosmic web signal is most detectable at these lower frequencies. It’s also important to note that the presence of fossil electron populations, which are not part of our model, may bias the cosmic web emission steeper still.
To further draw out the flux differences between these two populations, we next consider the distribution of flux values on the sky; this is dependent on the observing configuration, and in particular the beam resolution. In Figure 10, we bin pixels by flux density (log 10ΔS = 0.25) and show the proportion of the sky covered by these values, for each of the cosmic web and extragalactic sources when mapped according to the observing configurations in Table 7. Vertical dotted lines indicate the noise threshold of each observing configuration. As before, extragalactic sources are modelled simply as point sources. The peak flux distribution of the two populations differs by about two–three orders of magnitude for all configurations, with the vast majority of the cosmic web area of emission more than four–five orders of magnitude fainter than the bulk of the extragalactic emission. Only a small proportion of the cosmic web is either directly or statistically (e.g. via cross-correlation) detectable. If we refer back to Table 8, the final two columns show the peak (i.e. 100th percentile) and 99.9th percentile flux density values when observed in the MWA configuration. The peaks are bright enough that they are directly detectable with the MWA; their morphology, however, tends to resemble the aforementioned apparently point-like knots, making them hard to correctly identify. The extended emission surrounding these knots rapidly declines in brightness, and this can be seen in the 99.9th percentile flux values which are already a factor of ∼10 lower than the current best noise limits of the MWA.
Figure 10. The sky coverage of the cosmic web (blue) and embedded extragalactic (red) emission as a function of flux density binned in log 10 S = 0.25 increments. We show the sky coverage both as a function of redshift slice, as well as idealised observing configurations for MWA (150 MHz, 60′′), SKA Low (150 MHz, 10′′), and ASKAP (900 MHz, 20′′). Vertical dotted lines indicate the noise threshold for each configuration.
What is not apparent here is the variability of the cosmic web flux distribution between different realisations. In the lowest redshift slice, the relatively small volume contained within the 4° × 4°, 0
$ < $
z
$ < $
0.05 cone slice results in a highly variable flux distribution which reflects the real degree of cosmic variance of the very nearby Universe; thus for this nearest redshift slice, the realisations differ in total flux by five orders of magnitude depending on whether they sample from a void, a filament or, by chance, the core region of a massive cluster. By z
$ > $
0.1, this variability has significantly reduced as the enclosed volume of each redshift slice encompasses greater proportions of the simulation volume; the redshift slice 0.1
$ < $
z
$ < $
0.15, for example, differs by just one order of magnitude across all realisations.
7.2. Cosmic web angular scale
We next consider the characteristic angular extent of cosmic web sources, quantified through the radial autocorrelation function. This is described in Appendix C with the modification that we cross-correlate the map against itself (i.e. we set map A = B). Figure 11 shows the radial autocorrelation as a function of angular offset, R(θ), for a number of redshift slices out to redshift z = 0.3. This redshift range is typical of previous detection attempts, for example, Brown (Reference Brown2011) which only extended as deep as z = 0.05 and Vernstrom et al. (Reference Vernstrom, Gaensler, Brown, Lenc and Norris2017) which had a mean redshift depth of z = 0.2. The autocorrelation of each specific realisation is shown in grey, and the mean across all 10 realisations is shown in blue; the large spread of results, especially at low redshifts, is primarily a result of cosmic variance.
Figure 11. The radial autocorrelation as a function of angular offset, R(θ), for redshift slices (Δz = 0.05) out to redshift z = 0.3. The autocorrelation of each specific realisation is shown in grey, and the mean across all 10 realisations in shown in blue.
To compare results, Table 9 compiles the minimum, mean, and maximum FWHM values for each redshift slice. In the lowest redshift bin (0
$ < $
z
$ < $
0.05), θ
FWHM varies from as small as 55′′ to as much as 149′′; for redshift bin (0.05
$ < $
z
$ < $
0.1) the range is 28–58′′; and for the redshift bin (0.1
$ < $
z
$ < $
0.15) the range is much more narrow at 28–43′′. Given that these volumes are much smaller than the simulation volume (refer to Table 2), we believe that these ranges are likely representative of the real spread. The deepest three redshift bins in Figure 11 have mean values for θ
FWHM of 28′′, 22′′, and 22′′, respectively, however, these values are increasingly likely to be affected by systematics arising from the limited simulation volume.
Table 9. The characteristic angular scale of cosmic web emission, measured here by the FWHM of the autocorrelation of the cosmic web maps, for redshift slices of Δz = 0.05 out to z = 0.3. The minimum, mean, and maximum are calculated across each of the 10 realisations and are indicative of the expected cosmic variance between 4° × 4° fields.
These angular sizes are significantly less than the typical angular length of intracluster filaments at their respective redshifts, and this latter scale has often been used as a proxy for the expected angular extent of cosmic web emission (e.g. Vernstrom et al. Reference Vernstrom, Gaensler, Brown, Lenc and Norris2017). This smaller than expected angular extent arises as a result of the ‘knotty’ morphology of most emission structures, combined with our previous observation that the emission does not, in general, bridge cluster pairs but rather is restricted to halo shells. Whilst previous detection attempts, such as Vacca et al. (Reference Vacca2018), have deliberately chosen observing strategies to ensure sensitivity to extended and diffuse emission structures presumed to have angular scales that are multiple arcminutes in extent, these results suggest that a much finer resolution on the order of 20–40′′ is best for nearby cosmic web emission.
7.3. Cross-correlating the cosmic web
The cosmic web signal is too faint to detect directly, with the exception of outlier emission knots, and so one promising method is the radial cross-correlation method used by both Brown et al. (Reference Brown2017) and Vernstrom et al. (Reference Vernstrom, Gaensler, Brown, Lenc and Norris2017). In this method, a kernel image is first constructed and weighted based on where cosmic web emission is believed to concentrate across a map. This ‘best guess’ map for the cosmic web is then cross-correlated and radially averaged in an effort to bring out the cosmic web signal that is otherwise hidden beneath the image noise; a peak at or very near 0° offset suggests a possible detection of the cosmic web. In Vernstrom et al. (Reference Vernstrom, Gaensler, Brown, Lenc and Norris2017), the correlation kernel was produced by using galaxy density maps, which were believed to be a proxy for large-scale mass density and in turn for the synchrotron component of the cosmic web. A challenge with this approach is that the extragalactic radio population also correlates with this density map, producing an excess cross-correlation signal that is difficult to separate from any potential cosmic web signal.
We here attempt to reproduce this ‘false’ correlation by cross-correlating FIGARO with the mass density maps pertaining to redshift slices of depth Δz = 0.05, for the redshift range 0.05
$ < $
z
$ < $
0.3. Prior to cross-correlation, the mass density maps were smoothed with a Gaussian kernel having a FWHM of 1 Mpc pertaining to the mean redshift of the respective slice. In Figure 12, we show the cross-correlation of FIGARO for each of the idealised observing configurations with these mass density maps (solid line), as well as a ‘null’ result where the underlying cosmic web emission has been spatially flipped (dashed line). All AGN and SFG sources (spanning the full redshift range 0
$ < $
z
$ < $
8) have been represented as simple point sources and have been cleaned down to the noise threshold of the respective observing configuration. We observe that all results peak at zero, indicating a positive correlation with the mass density maps. In the case of the null result (dashed line), this correlation is solely the result of the AGN and SFG populations clustering similarly to the underlying mass distribution; this ‘false’ signal is precisely the issue Vernstrom et al. (Reference Vernstrom, Gaensler, Brown, Lenc and Norris2017) encountered. However, we also observe an excess correlation in the solid line, and this is due to the additional presence of the cosmic web. The degree of this excess is quite significant, accounting for more than half of the signal in the closest redshifts, and reducing somewhat at higher redshifts. We should note that the presence of other cluster emission that is not modelled in our simulation, such as radio halos, faint and diffuse AGN remnants, and radio phoenix, would further strengthen the ‘null’ result signal and reduce the relative excess; these populations are not well understood and it is very difficult at this time to model their additional contribution to the cross-correlation signal.
Figure 12. Cross-correlation of mass density with the FIGARO simulations, for a variety of redshift slices (solid lines), and compared with the ‘null’ case where the underlying cosmic web emission has been flipped (dashed). The excess correlation versus the null result therefore shows the cosmic web component of the cross-correlation result.
This spatial alignment of unrelated cluster emission and cosmic web emission makes detection extremely challenging. Nonetheless, we can still ask the question: if we had perfect knowledge of the distribution of the cosmic web, can we avoid the cross-correlation from being polluted by extragalactic population? To answer this, we have computed the cross-correlation of the full FIGARO maps (cosmic web, AGN, and SFG sources combined) with a kernel that is the cosmic web emission located in webshift slices of depth Δz = 0.05. We do this for all three observing configurations and have cleaned each map down to its respective noise threshold to produce a residual image. In Figure 13, we present the results of this cross-correlation and also additionally a ‘null’ result where we have spatially flipped the cosmic web signal in the map. The cosmic web kernel correlates strongly with itself, as expected, even in amongst the additional AGN and SFG residual sources; and moreover, this correlation seems to peak in the redshift slice 0.1
$ < $
z
$ < $
0.15 suggesting this is a sweet spot that maximises cosmological volume against the fading brightness that arises with increasing luminosity distance. The strongest correlation is observed in the MWA configuration, which by observing at 150 MHz takes advantage of the steep spectral index of cosmic web emission, but which is also advantaged by its lower resolution beam that is approximately at the characteristic scale of the cosmic web emission. Despite the disadvantage of observing at a higher frequency, the ASKAP configuration also produces a strong correlation signal suggesting this might well be a valid observing configuration.
Figure 13. Cross-correlation of the all FIGARO realisations with kernels of the known cosmic web signal for a variety of redshift slices (solid line), in comparison with a ‘null’ result where the cosmic web signal has been spatially flipped (dashed line). The cosmic web emission is easily detectable even in amongst the extragalactic sources. Note the lack of correlation peak for the null result, showing the extragalactic sources do not cluster in the same way as the cosmic web emission. Left: MWA, 150 MHz, 60′′ resolution. Middle: SKA Low, 150 MHz, 10′′ resolution. Right: ASKAP, 900 MHz, 20′′ resolution.
Most significantly, however, we observe no significant signal in the null results in Figure 13, and in fact in many cases this signal is not only weak, it is also weakly negative indicating anti-correlation. This is an important result: on average, the cosmic web signal is spatially separate and distinct from the extragalactic radio population, and they do not cluster in the same way. Whilst this provides for the possibility of devising a kernel that does not correlate with other cluster emission, of course any real detection attempt will not be blessed with perfect knowledge of the cosmic web distribution. The much more difficult task will be in devising kernels based only on approximate knowledge and proxies such as the mass distribution, but which do not strongly correlate with central cluster emission. We leave the development of kernels based only on limited knowledge of the Universe to future work.
8. Conclusion
We have released the FIGARO simulation: the first combined simulation of the cosmic web and its embedded extragalactic radio population spanning 10 unique 4° × 4° fields, out to a redshift of z = 0.8, and valid over a frequency range 100–1 400 MHz. In doing this, we have brought together the largest MHD simulation to date encompassing 1003 Mpc3—from which we have derived a model of synchrotron cosmic web emission, calibrated to match observed radio relic population statistics—and combined this simulation with the previously released T-RECS codebase to model the embedded extragalactic AGN and SFG embedded populations. Uniquely, these latter populations have been positioned such that they follow the underlying mass distribution of the MHD simulation, and thus cluster realistically alongside the associated cosmic web emission. We make these simulations publicly available as we believe they can be helpful in developing cosmic web detection techniques, as well as in modelling and constraining confounding signals that may arise in the course of these detection attempts.
In addition, we have provided an early analysis of the FIGARO simulation and its properties. Foremost amongst these results is to emphasise the spatial distribution and morphology that the cosmic web takes: as opposed to a sponge-like, filamentary structure that traces the underlying mass distribution, we find that the radio cosmic web is composed of ubiquitous, relic-like shells, principally in the spherical annuli 0.75 · r
200
$ < $
r
$ < $
1.5 · r
200 enshrouding dark matter halos. The filaments proper are largely empty of any significant emission. Moreover, the actual distribution of power in these shells is irregular and ‘knotty’. The brightest of these knots rise to the level of detection with the current generation of instruments whilst the surrounding emission rapidly declines in brightness, with a characteristic angular extent of
$\lessapprox 58^{\prime\prime}$
already by
$z \gtrapprox 0.05$
, much more compact than has been previously assumed. Indeed, some of these bright, knotty peaks of cosmic web emission are of an angular size that they could be easily confused with more compact emission sources, and indeed may already be present in the current generation of sky surveys.
Our simulation gives hope to the future detection of the cosmic web, primarily since the cosmic web emission and the extragalactic radio sources do not cluster in the same way. Our cross-correlation analysis shows the cosmic web emission is clearly detectable both at 150 and 900 MHz by way of cross-correlation, even in the presence of the much more luminous AGN and SFG sources; and that these embedded sources show negligible correlation, and in many cases, weak anti-correlation, with the cosmic web emission. In this case, the kernel used in the cross-correlation was the known cosmic web map itself; in a real detection attempt, a kernel will need to be devised with only approximate knowledge of the distribution of cosmic web flux. We hope the present simulation can aid in the development of such a kernel.
Besides the AGN and SFG populations, there are additional confounding factors that may impact detection techniques. One such factor in particular is the radio halo population which consists of large, extended low-surface brightness features that are centrally located in cluster and group cores, and which without careful consideration could contribute to false or exaggerated signals. The introduction of this population into FIGARO could prove a useful direction in future work.
Acknowledgements
F.V. acknowledges financial support from the ERC Starting Grant ‘MAGCOW’, no. 714196. The cosmological simulations on which this work is based have been produced using the ENZO code (http://enzo-project.org), running on Piz Daint supercomputer at CSCS-ETHZ (Lugano, Switzerland) under project s805 (with F.V. as PI, and the collaboration of C. Gheller and M. Brüggen). We also acknowledge the usage of online storage tools kindly provided by the INAF Astronomical Archive (IA2) initiative (http://www.ia2.inaf.it).
A. Two-point correlation function
To compare clustering properties, we make use of the two-point correlation function ξ(r). For a given set of points in a volume, this function estimates the probability of finding two points separated by some distance r with respect to a set of points that were randomly (uniformly) distributed. We have calculated this function using the estimator provided by Landy and Szalay (Reference Landy and Szalay1993):
\begin{equation} \xi\left(r\right) = \frac{ \langle D,D \rangle -2f \langle D,R \rangle + f^2 \langle R,R \rangle }{f^2 \langle R,R \rangle} \end{equation}
In this equation the notation
$\langle A, B\rangle$
indicates the number of unique pairs of points (a, b), where a ∈ A and b ∈ B, for which the condition r
$ < $
|a − b|
$ < $
r + dr is satisfied. The set D is our data, and contains the set of points for which we wish to calculate the two-point correlation function; the set R is a set of points randomly assigned to the same volume according to a uniform distribution; f = length(D)/length(R) is a normalisation constant to account for any differences in the number of points in each set. To calculate the error in this estimate, we use the ‘bootstrap resampling’ method (see, e.g. Ling et al. Reference Ling, Frenk and Barrow1986) whereby we repeat the calculation for a number of independently generated R, over which we can calculate the mean and standard deviation of the resulting values.
In the case of the angular two-point correlation function, the Euclidean distance is replaced by the apparent angular distance in projection on the celestial sphere.
B. Confusion limit calculation
The choice of lower flux threshold when generating the T-RECS catalogue was selected so as to fully simulate classical confusion noise in the frequency range 150–1 400 MHz assuming a maximum resolution of 8′′. We make use of the ‘probability of deflection’ or P(D) technique that is described in Vernstrom et al. (Reference Vernstrom2014), and references therein, to calculate the classical confusion noise. This technique allows for us to calculate, for any given point on the sky, the probability distribution of the flux density per beam based upon two inputs: a differential source count (dN/dS) and a synthesised beam shape. The classical confusion noise can then be estimated from the width of the peak in this distribution. In this case, we calculated the differential source count from the T-RECS catalogue itself, calculated independently at 150 and 1 400 MHz, and for the synthetic beam shape we assumed a circular Gaussian with FWHM of 8′′. In Figure B.1, we compare the P(D) distribution that results from a range of lower flux thresholds, that is to say, by setting
$\frac{dN}{dS_{1\,400}}(S < S_0) = 0$
for a range of S
0. At 150 MHz, for all the lower flux thresholds considered, we can see the central peak is well formed indicating that the confusion noise (at ∼18 μJy beam−1) is well simulated. At 1 400 MHz, on the other hand, we can see that for the lower flux thresholds of 5 × 10−7Jy and 10−6Jy the peaks are malformed, and only at 10−7Jy and below is the confusion noise (at ∼3 μJy beam−1) well simulated.
Figure B.1. The probability of deflection for the T-RECS catalogue at 150 and 1 400 MHz with an 8′′ circular Gaussian beam. Each curve shows the distribution for a different lower 1 400 MHz threshold cutoffs, showing that simulating sources down to 0.1μJy at 1 400 MHz is sufficient to simulate the classical confusion noise across this frequency range.
C. Radial cross-correlation
The radial cross-correlation of discrete maps A and B is constructed by first calculating its 2D cross-correlation function, which is defined as
\begin{equation} R(\Delta x, \Delta y) = \sum_{i, j} \frac{ \left(A(i, j) - \bar{A} \right)\left(B(i + \Delta x, j + \Delta y) - \bar{B} \right) }{ N(\Delta x, \Delta y) \cdot \sigma_A \sigma_B },\end{equation}
where A(i, j) is the (i, j)th component of map A,
$\bar{A}$
is the map mean, σ
A
is the standard deviation across the map, and N(Δx, Δy) is a normalisation function. In essence, the maps are offset from each other by Δx, Δy and we sum the product of all overlapping values, where the normalisation function simply counts the number of overlapping cells. The radial autocorrelation function is simply the radial average of this function, with
$r = \sqrt{\Delta x^2 + \Delta y^2}$
, and where in practice, we discretise radial values into bins and average over the 2D values of the function that fall within the bin.
