key: cord-0332888-11de1qeo
authors: Locatelli, Nicola; Vazza, Franco; Bonafede, Annalisa; Banfi, Serena; Bernardi, Gianni; Gheller, Claudio; Botteon, Andrea; Shimwell, Timothy
title: New constraints on the magnetic field in filaments of the cosmic web
date: 2021-01-15
journal: nan
DOI: 10.1051/0004-6361/202140526
sha: 5ce50d75bd589a2cc29220162ad06c93611c78c6
doc_id: 332888
cord_uid: 11de1qeo

Strong accretion shocks are expected to illuminate the warm-hot inter-galactic medium encompassed by the filaments of the cosmic web, through synchrotron radio emission. Given their high sensitivity, low-frequency large radio facilities may already be able to detect signatures of this extended radio emission from the region in between two close and massive galaxy clusters. In this work we exploit the non-detection of such diffuse emission by deep observations of two pairs of relatively close ($simeq 10$ Mpc) and massive ($M_{500}geq 10^{14}M_odot$) galaxy clusters using the LOw-Frequency ARray (LOFAR). By combining the results from the two putative inter-cluster filaments, we derive new independent constraints on the median strength of inter-galactic magnetic fields: $B_{rm 10 Mpc}<2.5times 10^2,rm nG,(95%, rm CL)$. Based on cosmological simulations and assuming a primordial origin of the B-fields, these estimates can be used to limit the amplitude of primordial seed magnetic fields: $B_0leq10,rm nG$. We advise the observation of similar cluster pairs as a powerful tool to set tight constraints on the amplitude of extragalactic magnetic fields.

At the largest scales of the Universe (≥ 10 Mpc), galaxy groups and clusters are connected by elongated distributions of galaxies called filaments and sheets which are believed to be also permeated by diffuse gas, and possibly by magnetic fields. Until now, a straightforward and direct detection of inter-galactic medium (IGM) and magnetic field (IGMF) has been prevented by the very low density of the plasma (n IGM ≤ 10 −4 cm −3 ) and its relatively low temperature (T IGM ≤ 10 7 K). However, increasing evidence (Nicastro et al. 2018; Macquart et al. 2020; Tanimura et al. 2020 ) is recently confirming the long-lived expectations for the warm-hot gas phase of the IGM (WHIM, with T WHIM ∼ 10 5 −10 7 , n WHIM ∼ 10 −5 −10 −4 ) to contain up to half of the baryon content at low redshift (Cen & Ostriker 1999; Davé et al. 2001; Reiprich et al. 2021) .

Accretion shocks are believed to reside along and within the filaments of the cosmic web as well as at the outskirts of galaxy clusters. These shocks are expected to amplify the magnetic fields and to accelerate particles up to relativistic energies (Ryu et al. 2008) . Their presence might then enable the detection of the WHIM through its synchrotron emission signature at radio wavelengths, and indeed the direct observation of the tip of the iceberg of this diffuse emission has already been obtained at radio frequencies (Govoni et al. 2019; Botteon et al. 2020b) . In these few cases the plasma conditions are still hotter and denser than what expected for the WHIM and the detected emission lays within the clusters virial radii. Further works in Send offprint requests to: E-mail: nlocat@mpe.mpg.de this direction will be helped in the near future by the promising new and upcoming radio facilities (the next generation Very Large Array ngVLA, the Karoo Array Telescope MeerKAT, the square kilometer array SKA-mid) and especially at very low frequencies (LOFAR, the Murchison Widefield Array MWA, SKAlow) where the emission should be brighter up to further out the clusters virial radii thanks to the expected spectral behaviour as S ν ∝ ν −1 with respect to frequency ν (Vazza et al. 2015) . A way to overcome sensitivity limitations is to quantify the Faraday effect induced by magneto-ionized plasma along the line of sight to a polarized background radio source and build a tomography of the WHIM by means of a grid of background sources (see Akahori et al. 2014; Vacca et al. 2016) . A thorough exploitation of this method currently suffers from the lack of large and dense grids of polarized sources, however it is expected to provide important results thanks to the upcoming radio facilities (Locatelli et al. 2018) . Complementary, recent upper limits on the IGMF intensity and scale have been derived from the cross-correlation of diffuse radio synchrotron emission with the underlying galaxy distribution Brown et al. 2017) , by crosscorrelating the difference in rotation measures of physically related pairs of extended radio galaxies, compared with the one derived from randomly paired and close lobes (Vernstrom et al. 2019; O'Sullivan et al. 2020; Stuardi et al. 2020 ) and by stacking full sky low-frequency radio images of close red luminous galaxy pairs (Vernstrom et al. 2021) .

Why is it important to assess the IGMF properties in the cosmic web at late times? The magnetic fields in galaxies and galaxy clusters, commonly observed today, arise from strong amplification from efficient MHD small-scale mechanisms (Ryu et al. 2008) which are responsible of a fast saturation of the fields, thus erasing information on their initial conditions, and in turn of their origin (Beresnyak & Miniati 2016) . Instead, in the WHIM Fig. 1 . Diagram of the method outlined in Sect. 2. Thick boxes highlight the most computationally expensive steps. Links labelled with the letter "i" are computed iteratively over the simulated pairs. environment, the amplification of primordial magnetic fields is found in simulations to be mainly driven by the field compression as its lines freeze into the plasma plus the contribution of small scale shocks. These mechanisms do not bring the field to saturation and provide a tool of assessing the history and original conditions of the field by means of the level of magnetisation observed today (Vazza et al. 2014 (Vazza et al. , 2015 Donnert et al. 2018 ). For the above reasons it is crucial to constrain the magnetic field in the WHIM in order to determine the original scenario for the large scale magnetic field origin and evolution in the Universe.

Cosmological MHD simulations predict the intensity of the IGMF at low redshift to range within 1 and 100 nG (Dolag et al. 1999; Brüggen et al. 2005; Vazza et al. 2017) . In this paper we introduce a novel method for a robust inference of an upper limit on the initial B 0 and current B values of the IGMF within the large-scale filaments of the cosmic web. The method explores the amount of diffuse emission detected at 144 MHz with the LOw Frequency ARray (LOFAR) telescope along the direction connecting two different pairs of galaxy clusters. We outline the method used to explore the upper limits on the IGMF into cosmological filaments in the following Sect. 2; we show its results in Sect. 3 and discuss their assumptions and implications in Sect. 4; we draw our conclusions in Sect. 5. Throughout this work we assumed a ΛCDM cosmological model, with baryonic and dark matter and dark energy density parameters Ω BM = 0.0455, Ω DM = 0.2265, Ω Λ = 0.728 respectively and a Hubble constant H 0 = 70.2km s −1 Mpc −1 .

In order to look for large scale emission from the cosmic web, we observed pairs of galaxy clusters and the putative inter-cluster filaments connecting them. The cluster pairs have been selected from the Meta-Catalogue of X-Ray Detected Clusters of Galaxies 1 (MCXC Piffaretti et al. 2011 ) by applying cuts in declination (δ ≥ 10 deg), redshift (z ≤ 0.3) and maximum angular separation (θ ≤ 5 deg). These values have been tailored to the proposed LOFAR observations. The two most promising pairs that, according to cosmological simulations, maximise the probability of a physical connection between the clusters in terms of total mass and separation (real and projected), have been proposed and observed at the LOFAR during Cycle 9 (Proposal 1 http://heasarc.gsfc.nasa.gov/W3Browse/all/mcxc.html Id:LC9_020). The most important properties of the two observed pairs of clusters are given in Table 1 .

In a nutshell, after calibrating, imaging and removing contaminating sources from the LOFAR data, we quantify the confidence of having observed (or not) diffuse emission from the inter-cluster filaments by injecting simulated diffuse emission produced by large scale (≥Mpc) accretion shocks, into the original radio visibility data (uvw) for a large (O(100)) subset of simulated filaments/cluster pairs and by imaging them as done for the real observations.

In this section we provide further details on the analysis performed on the actual observations, the simulated data set prepared for the injection and the injection procedure (also sketched by the diagram in Fig. 1 ).

We have observed the two cluster pairs RXCJ1659.7+3236-RXCJ1702.7+3403 (hereafter RXC_J1659-J1702) and RXCJ1155.3+2324-RXCJ1156.9+2415 (hereafter RXC_J1155-J1156) using LOFAR. The fields containing these two targets have been co-observed together with two pointings of the LOFAR Two-metre Sky Survey (LoTSS; Shimwell et al. 2017) taking advantage of the multi-beam capabilities of LOFAR. The observing setup of our observations thus follows that of LoTSS, namely 8 hr on-source time book-ended by two 10 min scans on the flux density calibrator in the frequency range 120-168 MHz using LOFAR in HBA_DUAL_INNER mode (see Shimwell et al. 2017, for details) . A first calibration and imaging run has been performed adopting the pipelines developed to analyse LoTSS pointings (Shimwell et al. 2017 (Shimwell et al. , 2019 , that aim to correct both for direction-independent and direction-dependent effects exploiting PREFACTOR van Weeren et al. 2016; de Gasperin et al. 2019) , KILLMS (Tasse 2014a,b; Smirnov & Tasse 2015) , and DDFACET (Tasse et al. 2018) softwares. In particular, we have used the improved version of the direction-dependent data reduction pipeline (v2.2 2 , the same used for the forthcoming second LoTSS data release DR2, Shimwell et al. in preparation) to produce images of the full LOFAR field-of-view at the central frequency of 144 MHz at high (6 ) and low resolution (20 , shown in Fig. 2 ) using a Briggs weighting scheme (robust=-0.5). We refer the reader to Tasse et al. (2021) for a thorough description of the steps performed by the pipeline. Using the sky models derived from the pipeline, we have subtracted the sources from the uv-data in two different fashions: either we have subtracted all sources (by means of their clean components) found in the high and low resolution maps or we have subtracted all sources detected in the high resolution image. The model components have been determined during the high and low resolution images deconvolution making use of the PYthon Blob Detector and Source Finder (PYBDSF; Mohan & Rafferty 2015) . The subtraction of the model components has been performed in the visibility domain, corrupting the model components by the directionindependent antenna gains obtained from the calibration. We have then produced dirty images from the subtracted data. The images include the residual contribute to the surface brightness resulting from model approximation plus artefacts associated with imperfect model and solutions (Im empty ) plus patches of faint extended emission in the case in which only sources from the high resolution model have been subtracted (Im diffuse ). The subtracted (dirty-)images Im empty at low resolution (20 ) have a rms noise floor of ∼ 160 and 240 µJybeam −1 for the two fields, respectively. The noise difference is consistent with the amount of flagged (i.e. discarded) data in the two observations.

We have extracted the simulated inter-cluster filaments from the suite of simulations of the cosmic web properties described in Vazza et al. (2019) performed with the cosmological MHD code ENZO 3 (Bryan et al. 2014) . They consist in a comoving 100 3 Mpc 3 box with a uniform grid of 2400 3 cells (and 2400 3 dark matter particles) with linear (comoving) resolution of 41.6 kpc per cell and dark matter mass m dm = 8.62 × 10 6 M per dark matter particle. Magnetic fields have been initialized at z = 45 as a uniform background of of B 0 = 0.1 nG and evolved at run-time using the MHD method of Dedner (Dedner et al. 2002) . We note that a uniform initial magnetic field here corresponds to a scaleinvariant spectrum in the models used for Cosmic Microwave Background (CMB) analysis (Aghanim et al. 2019; Paoletti & Finelli 2019) . We also note that the run is non-radiative and does not include any treatment for star formation or feedback from active galactic nuclei (AGN). While to a first approximation, these processes are not very relevant for the radio and X-ray properties of the peripheral regions of galaxy clusters and filaments,(e.g. Gheller & Vazza 2019) , the effect of outflows from AGN and galaxies can mix metals and magnetic fields at least out to the virial radius of clusters (e.g. Biffi et al. 2018) . Simulations by our group have showed that such effects is negligible on the radio emission on the scale of filaments (Vazza et al. 2017) , however only future and more refined simulations including galaxy 3 www.enzo-project.org formation-related effects in filaments, at a much higher resolution, will be able to assess with more certainty whether the limits obtained from observations (such as the one obtained in this very work) can be straightforwardly related to a limit on primordial seed fields.

We have produced synthetic maps of synchrotron radio emission assuming that diffusive shock acceleration (DSA, e.g. Kang et al. 2012 and references therein) accelerates a small fraction of thermal electrons swept by structure formation shocks up to relativistic energies, as in Vazza et al. (2019) . We have computed the radio emission from electrons in the downstream cooling region of shocks using the model of Hoeft & Brüggen (2007) and based on the shocks identified in post-processing in the simulation. The total acceleration efficiency at shocks, ξ e (M) (with M the Mach number) is assumed to be the combination of two variables: the kinetic energy flux dissipated onto the acceleration of cosmic rays, ψ(M), and the fraction going into electron acceleration, ξ e , giving ξ e (M) = ξ e · ψ(M). Following Hoeft & Brüggen (2007) , the radio emission in the downstream of each shock is directly linked to the power-law energy distributions N γ ∝ γ −p of electrons accelerated by the shock front during a cooling time, through the integrated radio spectrum of I(ν) ∝ ν −s , where s = (p − 1)/2 + 1/2, with p = 2(M 2 + 1)/(M 2 − 1). With this approach and for the range of M 5 shocks usually found within and around simulated filaments, as well as for the ≤ µG magnetic fields in filaments (e.g. Vazza et al. 2017) , the radio emission thus scales as I(ν) ∝ ξ e B 2 ν −2 . The baseline model used in this work assumes ξ e = 10 −2 , which is in line with DSA ex- Table 1 . Main parameters of the two pairs of galaxy clusters observed in this work, based on the Meta-Catalogue of X-Ray Detected Clusters of Galaxies (MCXC Piffaretti et al. 2011 ). In the last three columns, we give the 3-dimensional distance of between the two cluster centres, normalized by the radius of the two clusters), the angular separation of the two cluster centres, and the 3-dimensional length of the filament (considering the cluster to cluster distance).

pectations for the maximal acceleration efficiency of relativistic electrons by strong shocks (e.g. Hoeft & Brüggen 2007; Kang et al. 2012; Bykov et al. 2019b ) and is also compatible with the modelling of supernova remnants (e.g. Uchiyama et al. 2007; Bykov et al. 2019a ). However, we notice that in typically weak (M ≤ 4) shocks leading to radio relics in galaxy clusters (e.g. van Weeren et al. 2019), the acceleration efficiencies implied by the observed relic radio fluxes can be much larger (ξ e ∼ 0.1 − 1, see e.g. Stuardi et al. 2019 and Botteon et al. 2020a ), thus making our maximal value of ξ e = 10 −2 a conservative one. We notice however that very recent particle-in-cell (PIC) simulations (albeit in 1D and with some limiting assumptions) have derived a ∼ 5% electron acceleration efficiency by strong shocks (Xu et al. 2020) . For the remainder of the paper, the reader must thus bare in mind that our limits on B 10Mpc must be accordingly rescaled if a different value of ξ e is adopted. Fig. 3 gives an example of filaments connecting a massive cluster to other groups in its surrounding, in an ENZO cosmological simulation. A small but significant fraction of radio emission from shocks running on filaments connecting some of the pairs (e.g. M1-M2, M1-M4 and M1-M6 in Fig. 3) is above the detection threshold in LOFAR-HBA for our baseline model. In all cases, the detectable emission comes from relatively small and localised patches, extended a few ∼ 10 at most, with irregular shapes. Detecting the radio signal from cosmic filaments is indeed made challenging by the fact that the detectable fraction is just the tip of the iceberg of the wider 'radio cosmic web', which makes a morphological classification of the emission often ambiguous. Indeed advanced Deep Learning techniques have been proposed for the detection of the cosmic web in next radio surveys (Gheller et al. 2018 ). In the following section we use this model to constrain the amplitude of the ξ e B 2 combination based on our real LOFAR observations.

For each observed pair of clusters we have selected simulated pairs holding individual cluster masses M 500 > 10 13 M and linear (comoving) and projected angular distance of clusters in the pair within 20% deviation from the values of the observed pair (see the last and second-last columns in Table 1 ). We have obtained 171 simulated cluster pairs selected for RXC_J1659-J1702 and 139 pairs for RXC_J1155-J1156. The cluster pairs selected above mirror the separation selection criteria of our observations but include less massive clusters that may not involve a physical connection within one pair. We have thus analysed in addition a sub-sample of high-mass clusters (M 500 > 3×10 13 M ) for which a physical connection and the presence of a intercluster filament was verified either manually and through a high temperature cut T WHIM > 10 5 K of the WHIM within the intercluster filament. We refer to this sub-sample as best.

Taken a simulated cluster pair, a box has been drawn and extracted along the direction connecting the pair. The box has been rescaled to match the angular scale and comoving transverse distance indicated by the pixel size and redshifts of the observed clusters, and the intensity of synchrotron emission has been scaled to match its luminosity distance by conserving the total power (see the lower left panel in Fig. 4 for an example).

The flux density has been also multiplied by a constant factor

Under the assumption that the amplification of the magnetic field into the simulated filaments is affected by negligible small scale dynamo (Ryu et al. 2008 ) and that it is thus mainly driven by the adiabatic compression of the magnetic field lines following from flux freezing 4 into the plasma condensing during structure formation (Vazza et al. 2017) , the magnetic field B in the simulation is scalable with respect to B 0 and in turn the synchrotron emissivity is scalable with respect to f B ∝ B 2 0 (assuming a constant ξ e = 10 −2 ).

The rescaled simulated image of each mock pair of galaxy clusters has been injected into the source-subtracted Measurement Set (MS), following the procedure sketched in Fig.1 . In detail, each rescaled image has been first Fourier-transformed, then written into the MS and finally added to the visibilities of the source-subtracted sky using WSCLEAN (Offringa et al. 2014) .

We note that such injection does not take into account direction dependent effects that may act on the MS radio data. With the same software, the resulting data-set has been imaged and deconvolved with a 20 uv-taper and synthesized beam, and Briggs weighting scheme (Briggs 1995) with robust=-0.25 and 2000 minor cycles (see Fig. 4 for output examples). For realistic values of the normalisation parameter f B , the detectable emission is fragmented into small and sparse patches, associated with shocks internal to filaments. Therefore, we resort to statistical methods to assess the likelihood of each mock image to be compatible with our observed LOFAR fields. We compute the integral of the image power spectrum P S ≡ log 10 k max k min P(k ) dk where P(k ) is the power spectrum and k min and k max are determined by the image and beam size respectively (see Fig. 9 for an example). The sky model subtraction can leave bright residual artefacts depending on the goodness of the model used. These residuals can be as bright as ∼ 0.1 Jy beam −1 around point-like sources and they may dominate the integral of the image power spectrum P S . A zero-padding mask has been manually generated for each N. Locatelli et al.: Constraints on the cosmic web magnetic field pair of clusters in order to exclude those artefact from the computation of P S in all images.

From the cumulative probability distributions of the statistic P S resulting from all the source injections, we can access how likely is for a model to provide an expected value smaller than the one recovered from the observations. From the image Im empty , a 2.5 deg ×2.5 deg square centered on the axis of the cluster pair, in which all the sources (point-like plus extended) have been subtracted, we computeP S ≡ P S (Im empty ).P S corresponds to the total power in Im empty distributed over all scales from twice the beam size (k min ) up to half the image size (k max ). All images resulting from injection thus have P S equal or larger than the one computed for Im empty (red dotted lines in Fig. 5 ). The statistic P S resulting from the image in which diffuse emission has not been subtracted Im diffuse are indicated by the blue dashed vertical lines in Fig. 5 . Totals of 6 and 15 mJy of diffuse emission have been found in Im diffuse in excess of Im empty for RXC_J1659-J1702 and RXC_J1155-J1156 respectively. We outline the probabilities for the different models in Table 2. The Table reports also results from injection performed in the image plane (instead that in the uvw-plane) found in general to produce different probabilities of non-detection with respect to injection through the uvw-plane. We discuss this alternative method in Sect. 4.

The overall probability of a magnetic field model is simply the product of the probabilities (of the model to produce lower statistics) of the two cluster pairs, since the experiments have been run independently on each pair. We compute these probabilities in the "all" bottom lines in Table 2 . A model is more likely to be discarded, when its probability of having a smaller P S than in our LOFAR observations is very small (or very high alternatively).

Our main results can be so summarised: the primordial scenario with a seed magnetic field of B 0 10 nG has a small probability P(< P S ) 0.05 of explaining the small power excess in our observation of the RXC_J1659-J1702 and RXC_J1155-J1156 pairs, we then reject it with a confidence level (CL) of > 95%.

The models with a lower seed magnetic fields B 0 < 10 nG yield non-negligible (≥ 0.1) probabilities to produce a statistic equal to (o smaller than) the one observed.

By tightening the constraint on individual cluster masses and on the presence of a inter-cluster filament connecting the clusters, the model with B 0 = 10 nG can be rejected with even higher confidence CL> 97%.

If any of the patches of diffuse emission observed is produced by shocks in the WHIM, then the B 0 = 0.1 nG model is highly disfavoured, as it is basically unable to produce any detectable emission (i.e the probability in Table 2 of this model to produce less diffuse emission than what found in Im empty are always ≈ 1; we note that they have not been plotted in Fig. 5 ). Although we do not reject this scenario, we consider it implausible (see Sect. 4).

From the original simulation holding B 0 = 0.1 nG, we extract the probability distribution functions (PDF) of the magnetic field values B 10Mpc across the mock filaments selected according to the properties of the observed cluster pairs. We plot the resulting PDF(log B 10Mpc ) in Fig. 6 . Given the expected lack of dynamo amplification in the WHIM, the magnetic field distributions PDF(log B 10Mpc ) corresponding to the other B 0 models can easily be rescaled linearly with the input seed field. We find a skewed distribution encompassing B 10Mpc = 1.0 − 7.4 nG values (90% confidence range) with median B 10Mpc = 2.5 nG (equivalent to log B 10Mpc = −8.6) for the full sample. We note that the value of the magnetic field that produces the simulated synchrotron emission lays in the high part of the B 10Mpc distribution, as can be seen from the emission-weighted B 10Mpc dis- tribution in Fig. 6 . We also give in Fig.7 the average profiles of mass-weighted magnetic field strength for all simulated filaments extracted with the procedure above, for the two cluster pairs. On average, the profile of magnetic field is very uniform across 10 − 20 Mpc, with an average magnetic field along the line of sight for these objects of ∼ 2 − 3 nG and a tail of rare and massive filaments that can reach ∼ 10 nG (Fig. 7) . We find that the magnetic field along the direction connecting the clusters is on average enhanced with respect to volume-filling values. The enhancement is significant regardless whether the gas reaches temperatures typical of filaments (T > 10 5 K) or not, but the enhancement is larger for higher temperatures, as expected from compression of the magnetic field lines during the structures growth (e.g. Gheller et al. 2016) .

To interpret the results provided in Fig. 5 and Table 2 , we postulate three different assumptions that exploit the different type of source-subtraction performed in the analysis, and that can be used to derive different priors from our data (vertical lines in Fig. 5 ):

I : none of the residual diffuse emission after the point-like source subtraction (Im empty ) is produced by the shocked cosmic web;

II : all of the residual diffuse emission in excess of Im empty (i.e. Im diffuse ) is produced by the shocked cosmic web; III : at least some of the excess diffuse emission present in Im diffuse comes from the cosmic web.

We analyse their implications separately in the following paragraphs.

Provided that we can fix the ξ e acceleration efficiency at strong shocks (ξ e ≈ 10 −2 ), the assumption that none of the observed emission comes from cosmological shocks (I), produces in principle tighter constraints on B 10Mpc (and B 0 ), since P(P S (Im empty )) ≤ P(P S (Im diffuse )) always. In practice, the constraints are just slightly tighter due to the small amount of diffuse emission found into Im diffuse with respect to Im empty . Thus, under the first hypothesis that we did not observe the cosmic web emission, by scaling the B 10Mpc distribution obtained from B 0 = 0.1 nG to match the B 0 = 10 nG model (i.e. a factor ×100) we infer an upper limit to the current median IGMF into filaments of B < 0.25 µG with 95% confidence (the same confidence level that applies to the rejection of B 0 ≥ 10 nG models of Table 2 . Probabilities of obtaining a statistic P S lower than the one observed in Im empty and Im diffuse , computed from source injection performed in the uvw-plane and in the image-plane.

the primordial magnetic field scenario). By considering the best sub-sample of cluster pairs with higher masses and connected by an inter-cluster filament we can further improve the constraints on B by rejecting with a higher CL> 97% the B 0 = 10 nG model. Not all pairs of clusters are physically connected by a cosmic filaments, and in the absence of a detection in other wavelengths (e.g. via Sunyaev-Zeldovich or X-ray analysis, e.g. Govoni et al. 2019 ) one needs to resort to linking probabilities as a function of distance, which can be derived from cosmological simulations. For example, early Dark-Matter only simulations estimated that ∼ 80% of pairs in the same mass range of our clusters, and separated by ∼ 15 Mpc/h, are physically connected by filaments (Colberg et al. 2005) . With more recent and resolved simulations, also including gas physics, we can revise this number and tailor it to the exact mass difference and separation of our LO-FAR pairs. Using the algorithm outlined in Banfi et al. (2021) , we reconstructed the network of filaments connecting halos in our simulation, by checking for the actual presence of a matter bridge between pairs of clusters with a 25 or a 15 Mpc separation. This allows us to associate to each LOFAR observation a probability for the presence of a filament, through the ratio between the number of physical filaments in the simulation (best sample) over the total pairs of galaxy clusters found at a given distance. This results in 35% and 65% probability of having a filament between clusters at 25 and 15 Mpc separation, respectively. It shall be remarked, that even in the case without an actual gas connection between clusters, the region in between is far from being empty, because massive clusters at a relatively short separation also are indicators of a large cosmic overdensity, which is often associated with the presence of other filaments or threads of the cosmic web along the line of sight, which may account for a non-negligible radio emission. This still allows to exclude B 0 > 10 nG with high enough confidence. So, the absence of an actual filament does not dramatically affect the validity of the limits inferred from the full sample in this work.

The assumption that the excess diffuse emission present in Im diffuse with respect to Im empty is entirely due to shocked plasma of the WHIM can be readily tested by looking in detail at the diffuse emission patches which have been detected. In Fig. 8 we present close-up clippings of the diffuse patches found close to the pair RXC_J1659-J1702, taken from the low resolution LOFAR images (before source-subtraction). They are meant to help in assessing the nature of some of the diffuse emission, indicated by the dashed red circles in the panels of Fig. 8 . We also mark with green X symbols the position of sources already known from the VLA Faint Images of the Radio Sky at Twentycentimeters (FIRST) survey (Becker et al. 1995) . Most of the diffuse emission is plausibly linked to the lobes of radio-galaxies already detected at higher frequencies. Panel 7 (panels are numbered from 1 to 9 from left to right, top to bottom) shows what looks like either a radio lobe or an artefact linked to a lowfrequency point source. Panels 2, 5 and 6 instead show diffuse emission which is neither obviously linked to radio galaxies, nor to deconvolution artefacts. However, looking at their coordinates, sources 5 and 6 are found to be distant from the axis connecting the clusters, albeit within the imaged portion of the sky around the pair (see also Fig. 2 ). This makes their physical connection to the putative inter-cluster filament unlikely, even if the shocked cosmic web is expected to fill the space in between clusters in a non-trivial way, as shown in Fig.3 . Furthermore, the point-like source embedded into the diffuse emission in panel 6 is also found to be at a different redshift with respect to the cluster pair. For the above reasons, these patches can hardly be used in the comparison with the simulated inter-cluster filaments. The diffuse patches in panel 2 instead embed optical galaxies with redshift z = 0.087, 0.093, consistent with the cluster pair z = 0.095 − 0.101, however they are likely dying faint radio lobes, with no FIRST counterpart. Interestingly, the positions of Sloan Digital Sky Survey (SDSS) sources (orange crosses) in the panels of Fig. 8 seem to cluster close to the diffuse radio emission. This trend is actually expected if the radio signal is caused by radio galaxies and lobes which show the tendency to be found in clusters or groups of galaxies (Prestage & Peacock 1988; Allington-Smith et al. 1993; Zirbel 1997; Wing & Blanton 2011) . This behaviour additionally disfavours the link between the excess diffuse emission and the cosmic web. All in all, since there is no easy way to cross check all the different patches, we still compute the statistics for the most conservative scenario by assuming that the level of observed diffuse emission in excess of Im empty is entirely due to the cosmic web. In this case, the level of confidence associated to the rejection of the same models loosens. However, the models rejected by casting hypothesis II are the same ones resulting from hypothesis I, though with slightly lower or even equal CL. (e.g. for the B 0 = 10 nG best model the CL for its rejection decreases from 98% to 97% while for the full sample remains unchanged. Furthermore, since hypothesis II has been falsified already by the examples described above and shown in Fig. 8 , hypothesis I is strengthen in favour of hypothesis II and we thus refer to the former in order to draw our conclusions.

For completeness, a third additional and interesting way of interpreting our data with respect to the simulation results is the complementary hypothesis to the first one: we assume that at least some of the diffuse emission in excess of Im empty comes from the cosmic web. The associated probabilities is then trivially P(> P S ) = 1 − P(< P S ). In this case we are not interested in the level of diffuse emission in Im diffuse , since we want to produce at least the one in Im empty . This scenario, though disfavored, cannot be discarded a priori since this would imply checking (e.g. through cross-correlations) all the different patches of diffuse emission in Im diffuse and prove that all of them are not connected to the emission from the Cosmic Web, it is then instructive to inspect its implications. Under the assumption that we did see the cosmic web emission at least in part, then the B 0 = 0.1 nG model is ruled out with high confidence ≥ 99% since it is not able to produce any observable emission brighter than the noise level of our LOFAR observation. In this scenario, B 0 > 0.1 nG can be set as a lower limit to the primordial mag- netic field intensity and in turn B > 2 nG as the median value for the magnetic field into filaments today.

While checking that the source injection procedure (presented in Sect. 2 and sketched in Fig. 1 ) is actually needed in order to derive robust limits on B 10Mpc and B 0 we also demonstrate that the method is essential to interpret observations in details by means of the outcome of simulations when dealing with radio data. With this respect, we have produced the same statistic P S for the simulations directly added to the residual image Im empty in terms of simple image sum, rather than following the central FFT + write + sum procedure involving visibilities. This procedure is much simpler and faster (shortening the computing time of a factor ∼ 600). In Fig. 9 we plot the power spectra resulting from the injection in the RXC_J1659-J1702 field of one source as example (images of the same source are shown in Fig. 4 ) in order to inspect differences between the injection through the uvw-(black lines) and image-plane (grey lines), for the different B 0 models. As can be seen by comparing the black and grey lines, when the injection is performed within the image-plane the level of simulated emission at large scales is generally underestimated. As a consequence the models are consistent with the data with different probability up to ±30%P(< P S ) (see the values in Table 2 ). We interpret the difference in the results as due to the lack of model convolution with the instrument's Point Spread Function (PSF). In addition, the lack of convolution of the emission with a visibility weighting scheme able to maximise the evidence for extended diffuse emission into the data may also play a similar role. As far as an upper cut on the scales of the emission (corresponding to a lower bound on the baseline length in radio interferometers) is taken into account, and detailed power spectrum information does not constitute the largest budget of uncertainty in one analysis (in our case is the scatter in the properties of the -unknown-inter-cluster WHIM), the image sum is a much faster approach than the source injection through the uvwplane, however it shall be used with caution as results are biased by a different sampling of the scales. The strength of the bias depends either on the sampling (window) function and the source power spectrum.

As a final caveat, our analysis assumes that for strong shocks in and around filaments, the acceleration efficiency of electrons is the one suggested by DSA, i.e. ξ e ∼ 10 −2 . This assumes, in turn, that despite the rather low particle density and magnetisation, shocks can form and undergo particle acceleration similar to what is already observed for the outer regions of galaxy clusters in form of giant radio relics (see van Weeren et al. 2019, for a review). Moreover, our analysis assumes that the acceleration of electrons at shocks can proceed independently on the obliquity between the upstream magnetic field and the shock normal. However, recent numerical works by Banfi et al. (2020) have shown that shocks surrounding the cosmic web are more often quasi-perpendicular than random chance, as an effect to the peculiar gas velocity flow following the formation of filaments. In this case, the vast majority of shocks in filaments are quasiperpendicular and thus likely to be suitable for efficient electron acceleration (Xu et al. 2020) . Furthermore, Masters et al. (2017) have recently reported a significant electron acceleration by the strong quasi-parallel shock while crossing the Saturn bow shock by the Cassini space mission, i.e. in plasma conditions similar to the intra-cluster medium. The acceleration seems to occur in the portion of the shock where upstream cosmic-ray streaming instabilities generate perpendicular small-scale magnetic field components, leading to particle acceleration.

In this work, we attempt for the first time to combine dedicated LOFAR-HBA observations of inter-cluster filaments and numer-ical simulations of the magnetic cosmic web, in order to derive upper limits on the magnetisation of the WHIM.

While our LOFAR observations do detect patches of diffuse emission of unclear origin their morphology does not allow us to firmly associate the origin of the most prominent ones to the cosmic web. However, the presence of a faint diffuse large scale excess in comparison with numerical models allows us to derive inferences on the average magnetisation of such filaments, and possibly on the allowed initial amplitude of primordial seed magnetic fields. As a main outcome of our work, by fixing ξ e = 0.01 for strong shocks, we derive an upper limit for the median magnetic field strength in filaments connecting massive galaxy clusters: B 10Mpc < 0.25 µG. Based on the dynamical evolution of magnetic fields given by present simulations (which is mostly dominated by simple compression of magnetic field lines), this also implies an upper limit of B 0 < 10 nG on the amplitude of primordial seed fields. The estimates above rely on the assumption that our observations constitute non-detections of diffuse emission from the cosmic web (hypothesis I).

As a mutually exclusive interpretation of our data, if some of the detected emission partially came from the shocked WHIM (hypothesis III), this would imply a median magnetic field of order of B 10Mpc ≥ nG (see e.g. Fig. 6 ). This would be an important outcome as it would also possibly indicate primordial magnetic fields with intensity B 0 ≥ 0.1 nG.

The hypothesis that all of the excess diffuse emission detected in our maps came from the cosmic web (hypothesis II) can be easily rejected even by a visual check of some of the sources. Given the uncertainties connected to our method and the limited statistics of 'detections' in our sample, we favour the first interpretation of our results (hypothesis I, setting B 10Mpc < 0.25 µG and B 0 < 10 nG).

To put our new limits in comparison with other recent works (Hackstein et al. 2016; Pshirkov et al. 2016; Vernstrom et al. 2017; Brown et al. 2017; Vernstrom et al. 2019; O'Sullivan et al. 2020 ; Natwariya 2020, and Paoletti & Finelli 2019 for joint BI-CEP2/Keck -Planck 2018 updated results), we show them in Fig. 10 separating the limits inferred for the IGMF or the magnetic field intensity of the WHIM (red arrows) and for the primordial magnetic field intensity B 0 (blue arrows). We note that our limits are still in agreement with the recent limits 0.134 < B 0 /nG < 0.316 set by the level of excess diffuse emission observed by ARCADE2 and EDGES 21cm line experiments (Natwariya 2020). Furthermore, an apparent tension seems to arise between our lower limit to B 10Mpc into filaments and the one derived from other probes such as the level of anisotropy in the arrival direction of charged ultra-high-energy cosmic rays used to limit the average amplitude of magnetic fields in voids to ≤ 1nG (Hackstein et al. 2016) , or from the non-detection of a trend of rotation measures from distant radio sources with respect to redshifts (Pshirkov et al. 2016) . Although computed over similar linear scales ≥ Mpc and globally refer to the IGM, they can still hardly be directly compared since referred to different IGM environments (e.g. voids, filaments, averaged). Interestingly a recent work suggests that primordial magnetic fields with amplitude ∼ 0.1 nG would possibly alleviate the existing tension between cosmological and standard candle-based estimates of H 0 (Jedamzik & Pogosian 2020) .

While it is hard to derive conclusive limits from these data, as no robust detection (although tentative) of the diffuse emission from the cosmic web can be claimed, this first attempt stresses the potential of low-frequency radio observations in constraining extragalactic magnetic fields, and its relevance to the study of cosmic magnetogenesis. With the analysis and the values ob- Fig. 10 . Summary of the current upper and lower limits to the volumeaveraged IGMF (red arrows), B in the filaments' WHIM (green arrows) and B 0 (blue arrows). tained in this work, we can forecast to produce tighter constraint than the ones posed by CMB experiments by covering a ∼ ×10 larger sample of cluster pairs similar to the ones analysed here even in the case of other non-detections.

The code used to produce the simulations discussed in this paper is public (enzo-project.org). Significant subsets of the simulations are publicly available at https://cosmosimfrazza.myfreesites.net/the_magnetic_cosmic_web, while larger data set can be shared upon request

American Astronomical Society Meeting Abstracts

Proceedings of the International Symposium on Grids and Clouds (ISGC)

PyBDSF: Python Blob Detection and Source Finder Natwariya, P. K. 2020, arXiv e-prints

A&A proofs: manuscript no. aa

We thank Klaus Dolag for providing corrections and helpful scientific feedback. NL, and FV acknowledge financial support from the ERC Starting Grant "MAGCOW", no. 714196. ABon acknowledges financial support from the ERC Starting Grant "DRANOEL", no. 714245. ABot acknowledges support from the VIDI research programme with project number 639.042.729, which is financed by the Netherlands Organisation for Scientific Research (NWO). Radio imaging made use of WSClean v2.6 (Offringa et al. 2014 ) and DDF (Tasse et al. 2018) . NL wishes to thank S. Carozzi and A. Crotti for psychological aid during and after the lock-down that followed the covid-19 pandemic situation, when this work took shape. This paper is based (in part) on data obtained with the International LOFAR Telescope (obs. ID LC9_020, PI F.V.) and analysed using LOFAR-IT infrastructure. LOFAR (van Haarlem et al. 2013) is the Low Frequency Array designed, constructed by ASTRON and collectively operated by the ILT foundation. These data were (partly) processed by the LOFAR Two-Metre Sky Survey (LoTSS) team. This team made use of the LOFAR direction independent calibration pipeline (https://github.com/lofar-astron/prefactor) which was deployed by the LOFAR e-infragroup on the Dutch National Grid infrastructure with support of the SURF Co-operative through grants e-infra 160022 and e-infra 160152 (Mechev et al. 2017 , PoS(ISGC2017)002). The cosmological simulations were performed with the ENZO code (http://enzo-project.org), which is the product of a collaborative effort of scientists at many universities and national laboratories.FV acknowledges the usage of computational resources on the Piz Daint supercomputer at CSCS-ETHZ (Lugano, Switzerland) under project s805, and he also acknowledges the usage of online storage tools kindly provided by the Inaf Astronomica Archive (IA2) initiave (http://www.ia2.inaf.it).