Forming intracluster gas in a galaxy protocluster at a redshift of 2.16
Cosmology
In this work, we consider a spatially flat Λ Cold Dark Matter cosmological model, with Ω_{M} = 0.30, Ω_{Λ} = 0.70 and H_{0} = 70.0 km s^{−1} Mpc^{−1}. At the redshift of the Spiderweb complex (z = 2.156), 1″ corresponds to 8.29 kpc.
ALMA observations and reduction
An extensive observational campaign was performed during Cycle 6 to obtain a detailed view of the Spiderweb complex in Band 3 (project code: 2018.1.01526.S, principal investigator A. Saro). The data comprise measurements from the main 12m ALMA^{31} array in three different configurations (C431, C433 and C436), aimed at providing a highdynamicrange view of the structure, as well as from the 7m ACA^{32} (also known as the Morita Array), complementing the ALMA observations over SZrelevant scales (that is, over a uv range of 2.2–17.5 kλ, corresponding to scales 767–97 kpc at the redshift of the Spiderweb galaxy). The spectral setup for all the configurations was tuned to cover the frequency range 94.5–110.5 GHz, split over four 2GHzwide spectral bands centred at 95.5 GHz, 97.5 GHz, 107.5 GHz and 109.5 GHz, respectively. In particular, the last window targets the line emission resulting from the J = 3–2 transition from the carbon monoxide (CO; rest frequency 345.796 GHz). As we are interested in modelling only the continuum component of the observed signal, we conservatively exclude all the visibilities from the spectral window expected to contain the redshifted CO J = 3–2 line. In fact, excluding only the channels corresponding to the specific emission line would provide a slight improvement in the overall statistics of the analysed data. Previous studies^{33,34,35} of the molecular content of the nuclear region around the Spiderweb galaxy and of protocluster galaxies have shown, however, that the cold molecular gas within the Spiderweb complex is characterized by large velocity dispersion as well as broad differences in the systemic velocities of the member galaxies. Faint tails and offset components may hence contaminate the continuum signal in any channels in the proximity of the emission line, potentially affecting the model reconstruction.
Data calibration was performed in the Common Astronomy Software Applications^{36} (CASA; https://casa.nrao.edu/) package version 5.4.0 using the standard reduction pipeline provided as part of the data delivery. The direct inspection of output visibility tables highlighted no clear issues with the outcome of the pipeline calibration and we therefore did not perform any extra flagging or postprocessing tuning. The resulting rootmeansquare (RMS) noise levels of the observations are estimated to be 4.10 μJy beam^{−1}, 29.1 μJy beam^{−1} and 14.2 μJy beam^{−1} for the C431, C433 and C436 ALMA measurements, respectively, and 32.0 μJy beam^{−1} for the ACA data.
To obtain better knowledge of the spectral properties of the measured signals, we further include in our analyses archival Band 4 ALMA (project code: 2015.1.00851.S, principal investigator B. Emonts) and ACA (project code 2016.2.00048.S, principal investigator B. Emonts) observations. In both cases, we use the calibrated measurement sets provided by the European ALMA Regional Centre network^{37} through the calMS service^{38}. The achieved RMS noise levels amount to 99.0 μJy beam^{−1} and 6.29 μJy beam^{−1} for the ACA and ALMA data, respectively. As with the Band 3 measurements, we exclude from our analyses the spectral windows in Band 4 covering the CO J = 4–3 and [CI] ^{3}P_{1}–^{3}P_{0} emission lines^{33}.
All interferometric images presented in this work are generated using CASA package version 6.3.0.
Nested posterior sampling
To obtain a statistically robust detection of the potential SZ signal in the direction of the Spiderweb complex, we use the approach already used in refs. ^{16,39,40,41} in the context of ALMA+ACA studies of the SZ signal. In brief, we perform a visibilityspace analysis, which allows for exactly accounting for the nonuniform radiointerferometric transfer function, as well as taking advantage of the Gaussian properties of noise in the native Fourier space. Any extended model component is created in image space, taking into account the proper frequency scaling and primarybeam attenuation for any fields and spectral windows used in the analysis, and is then projected onto the visibility points by means of a nonuniform fast Fourier transform algorithm based on convolutional gridding (as implemented in the finufft library^{42}). Instead, in Fourier space, pointlike sources are trivially represented by a constant function with amplitude equal to the source flux corrected for the primarybeam attenuation at the source position and with a phase term defined by the offset between the point source and the phase centre of the observations. To allow for Bayesian model selection and averaging, the sampling of the posterior distribution is performed by means of the nested sampling algorithm^{43,44}. We specifically exploit the implementation provided in the dynesty (ref. ^{45}) package, which allows for robustly extending the sampling problem to moderatedimensional and highdimensional models. For any further details on the model reconstruction, we refer to the discussion provided by Di Mascolo et al.^{16,39}.
Obtaining a thorough description of the smallscale, complex morphology of the extended radio signal from the Spiderweb galaxy would require performing a pixellevel model inference (see ‘Sparse imaging’ section below), resulting in a posterior probability function with extreme dimensionality. Nested sampling techniques generally show better performances than Monte Carlo Markov chain algorithms in the case of moderatedimensional problems (mitigating the impacts of the socalled ‘curse of dimensionality^{44,46}). Still, sampling from highdimensional posterior distributions may easily become computationally intractable owing to the complexity of estimating highdimensional marginal likelihoods. Therefore, we decide to perform a first modelling run only on a largescale subset of the available data. The giant Lyα nebula observed to surround the Spiderweb galaxy is in fact expected to be confined within a diffuse halo of hot intracluster gas^{8,9,47,48,49}. In turn, the SZ footprint of potential intracluster gas within the Spiderweb protocluster should be expected to extend over characteristic scales ≳10″. We therefore introduce an upper cut in the visibility space at a uv distance of 65 kλ, whose corresponding angular scale is roughly twice the transverse width of the jet structure (that is, the size measured along the direction perpendicular to the jet direction) observed when imaging only the highresolution ALMA data from the C433 and C436 observations. Such a choice makes the jet signal spatially resolved only along the jet axis and allows for describing this as a limited collection of pointlike sources. As a result, the number of parameters required to model the observations remains limited and the analysis computationally manageable, while allowing for using nested sampling to obtain statistically meaningful information on our model inference.
The selection of the specific number of compact components required to model the radio source was performed by means of Bayesian model selection. In particular, we consider as the most favoured model set the one for which the introduction of an extra pointlike term would have caused a degradation or only a marginal improvement in the logevidence^{50} (that is, ((log ,{Z}_{n+1}log ,{Z}_{n})lesssim 0.50), in which n denotes the number of model components). For all of the components, as already described above, we assume the spatial morphology to be described by a Dirac δ function (that is, a constant function in Fourier space with nonnull phase term) and the source fluxes by a powerlaw spectral scaling. The priors for all the parameters—right ascension, declination, flux and spectral index—are described by uniform probability distributions. In particular, the right ascension and declination are allowed to vary within the area of r ≃ 1′ enclosed within the first null of the antenna pattern for the Band 3 ALMA data at the highest available frequency (that is, 107.5 GHz). To avoid label switching and force mode identifiability, we impose a further ordering prescription^{51} to the right ascension parameters. The flux and spectral index parameters are instead allowed to vary within uniform prior probability distributions. In particular, we constrain the source fluxes to be nonnegative and assume an upper limit of 10 mJy on each amplitude, around an order of magnitude larger than the emission peak in the ALMA map (1.54 mJy; see top panel of Extended Data Fig. 1). For the spectral indices, we consider a range [−5, 10], arbitrarily wide and set to cover both the cases of power spectra with negative and positive slopes, consistent with synchrotronlike and dustlike spectral properties, respectively. The prior limits are intentionally set to extend well beyond the values expected for such cases, to avoid overconstraining of the source spectral properties while allowing for a quick diagnostics of the actual constraining power of the available data with respect to spectral information.
To describe the pressure distribution of the ICM, we instead use a generalized Navarro–Frenk–White (gNFW) profile^{52}, widely shown to provide an accurate description of the average pressure distribution of the intracluster gas. In particular, we use different gNFW formulations from the literature. A summary is provided below.

–The universal profile (hereafter, A10 UP) by Arnaud et al.^{28} derived from the reconstruction of the pressure distribution in galaxy clusters from the REXCESS^{53} sample and the equivalent model obtained from the subset of systems with clear evidence of a disturbed morphology (A10 MD). Although this profile is calibrated on massive local systems (10^{14} M_{⊙} < M_{500} < 10^{15} M_{⊙}, z < 0.2), it is the base model used for the mass reconstruction in largescale SZ surveys and allows for a straightforward comparison with the literature.

–The pressure model^{54} (M14 UP) obtained from the Xray analysis of a highredshift (0.6 < z < 1.2) subsample of galaxy clusters detected by the South Pole Telescope^{55,56,57} (SPT). As for the previous case, we further consider the pressure profile (M14 NCC) reconstructed by excluding all the systems with clear evidence for the presence of a cool core, generally indicative of a more relaxed dynamical state. To our knowledge, this model represents the highestredshift, observationally motivated pressure profile available at present, in turn potentially providing a better description of the pressure distribution in the Spiderweb system than the A10 models.

–The median pressure profiles reconstructed from the OverWhelmingly Large Simulations (OWLS) suite of cosmological hydrodynamical simulations, cosmoOWLS^{29}, and considering different prescriptions for the physical model of AGNdriven heating of the intracluster gas. In particular, we consider the OWLS^{58} reference model (L15 REF), as well as the two AGN models (L15 8.0 and L15 8.5; we refer to Le Brun et al.^{29} for details). These models present two main advantages. First, they are built on simulated halos whose mass and redshift ranges (2 × 10^{12} M_{⊙} ≲ M_{500} ≲ 3 × 10^{15} M_{⊙}, 0 < z < 3) broadly overlap with the properties of the Spiderweb protocluster. Second, at z ≃ 2, the Spiderweb complex sits in a transitional phase for AGN feedback^{59,60} and the different flavours of the cosmoOWLS pressure model allow for directly using the SZ effect to explore different AGN scenarios.

–The massdependent and redshiftdependent extended pressure model^{61} (G17 EXT) calibrated on simulated galaxy clusters from the set of Magneticum Pathfinder hydrodynamical simulations (http://www.magneticum.org/). Although this profile is computed on massive galaxy clusters (M_{500} > 1.4 × 10^{14} M_{⊙}), it provides an explicit model for taking into account the departure from universality and selfsimilarity as a function of mass and redshift.
In all cases, the free parameters defining the gNFW profiles are the planeofsky coordinates of the model centroid and the mass parameter M_{500}. As for the radio model, the right ascension and declination parameters are prescribed to vary within the region encompassed by the first null of the Band 3 ALMA primary beam. For the mass M_{500}, we consider a loguniform distribution over the range [10^{12} M_{⊙}, 10^{15} M_{⊙}], to facilitate the posterior exploration over such a wide range of order of magnitudes.
We also tried fitting the coolcore versions of the A10 and M14 profiles above but found these to be systematically disfavoured ((Delta log ,Zlesssim 15.5)) in comparison with the listed models. This is, however, not surprising, as the presence of a wellformed cool core would be hardly consistent with the inherently disturbed nature of a protocluster complex.
Finally, we account for any potential systematics with data calibration by considering a scaling parameter for each of the measurement sets used in our analysis. For these, we assume normal prior distributions, with unitary central value and standard deviation of 5%, as reported in the ALMA Technical Handbook for the considered observing cycles.
Results
With regards to the extended signal from the Spiderweb galaxy, the criterion introduced above for the selection of the number of pointlike components supports the case for a total of eight distinct components over the entire search area (we summarize the key information on the results of the pointlike modelling in Extended Data Table 1). In particular, two components (ID1 and ID2) are found to be spatially consistent with the position of known protocluster members^{25,34,35,49,62,63,64}—that is, ERO 284 (ref. ^{12}) and HAE 229 (refs. ^{65,66}), respectively—located around 250 kpc west of the Spiderweb galaxy. The remaining components instead uniquely describe the radio emission associated with the Spiderweb galaxy, with one pointlike component (ID8) specifically corresponding to the bright lobe of the eastern radio jet and another (ID6) being nearly coincident with the Spiderweb galaxy itself (Extended Data Fig. 1). All the components describing the extended signal exhibit a negative spectral index, consistent with the synchrotron origin of the emission. The bestfit estimates highlight a spatial variation consistent with what is observed at lower frequencies in the Karl G. Jansky Very Large Array (VLA) data^{6,7,8,9}, which show that the spectrum of the eastern lobe is, on average, steeper than that of the Spiderweb galaxy and the western jet. The specific values are also in rough agreement with the results from the VLA analyses. A onetoone comparison is, however, not practicable, owing to the inherent highfrequency spectral steepening induced by radiative losses and the different modelling approaches. The two offset sources, on the other hand, are both characterized by positive spectral indices. Such a trend is possibly indicating a dominant contribution from thermal dust emission already at about 100 GHz and is in agreement with the potential presence of massive dust reservoirs in the galaxies, as already verified for HAE 229 (refs. ^{65,66}).
A summary of the inferred parameters for the different SZ models is instead provided in Extended Data Table 2. Despite the more or less substantial differences in the reconstructed parameters, all the adopted pressure profiles resulted in statistically consistent SZ models, not allowing a statistically motivated selection of a specific description. All the assumed models, with exception of the L15 8.0 and L15 8.5 cases, provide mass estimates M_{500} ≃ 3 × 10^{13} M_{⊙}. In fact, we note that the L15 8.5 profile from Le Brun et al.^{29} results in a mass of about 8 × 10^{13} M_{⊙}, consistent with many of the dynamical estimates reported in the literature^{11,13,15} for the Spiderweb protocluster. Considering that the profile is based on simulations with the prescription for an intense heating of the ICM owing to AGN feedback, the result might hint to a crucial role of the active core of the Spiderweb galaxy in shaping the intracluster/circumgalactic medium. Nevertheless, this model is statistically equivalent to many others in our sample, limiting the statistical relevance of the above considerations. For the same reason, any attempts of performing a comparison with past ICM studies of similarly highredshift systems—for example, XLSSC 122 (ref. ^{67}) or Cl J1449+0856 (ref. ^{68}), the clusters at the highest redshift known so far with a direct SZ measurement—would not be statistically meaningful.
Finally, we note that, in our analyses, all the scaling parameters are found to be broadly consistent with unity. In particular, the scaling factors are measured to be equal to (1.0{4}_{0.04}^{+0.04}) and (1.0{1}_{0.04}^{+0.03}) for the Band 3 ACA and ALMA observations, respectively, whereas—in the case of the Band 4 data—we find parameters of (1.0{0}_{0.07}^{+0.06}) for ACA and (1.0{3}_{0.04}^{+0.04}) for ALMA (the reported values are given by the Bayesian model averages for all models considered in this work; see the ‘Dependence of the SZ significance on the number of pointlike components’ section below for a discussion).
Comparison with masses from previous studies
Performing a proper comparison of the SZderived mass estimates (Extended Data Table 2) with independent measurements from the literature is nontrivial. The dynamical masses for the Spiderweb protocluster are in fact based on velocitydispersion estimates that might trace specific, yet not wellidentified substructures and that span almost an order of magnitude: from 204 km s^{−1} for one peak in the velocity distribution of Lyα emitters in the Spiderweb field^{11} and up to 1,360 km s^{−1} as measured for all satellites within 60 kpc from the Spiderweb galaxy^{14}. Nevertheless, the SZbased mass estimates presented in this work are much lower (a factor of about 2–4, depending on the model) than the dynamical values^{11,14,15} found in the literature for the mass of the whole protocluster structure. In fact, the same studies reported evidence for a complex velocity structure within the central region of the Spiderweb complex, hinting at the possibility for the system to be experiencing a major merger and still accreting large amounts of material from surrounding filaments. The fact that the Spiderweb system is embedded in a largescale filamentary structure was confirmed by widefield CO J = 1–0 mapping^{35}. As such, even the very core of the Spiderweb protocluster might not be fully virialized. Accordingly, the integrated SZ signal Y_{SZ} = Y_{SZ} (<5r_{500}) we measure from the ALMA+ACA observations—({Y}_{{rm{SZ}}}^{{rm{ALMA}}}=(1.6{8}_{0.32}^{+0.35})times 1{0}^{6},{{rm{Mpc}}}^{2}) (see Extended Data Table 2)—is, for instance, a factor of 3.5 times lower than that expected, taking the mass inferred from the velocitydispersion measurement σ_{v} ≃ 683 km s^{−1} reported by Shimakawa et al.^{15} for the galaxies within the 0.53Mpc region surrounding the Spiderweb galaxy.
This large difference between the expected SZ signal from velocitydispersion measurements and the observed ALMA+ACA SZ integrated flux could be the result of a scenario in which the measured SZ signal is dominated by the contribution from the most prominent (sub)halo. In fact, the integrated SZ flux Y_{SZ} scales steeply as a function of mass M, meaning that the SZ signal from a single halo would be larger than the sum of the SZ flux from a complex system of subhalos whose masses amount overall to the same value M. Under the assumptions that the Spiderweb protocluster is composed of several interacting substructures and the measured lineofsight velocity dispersion σ_{v} is providing an unbiased estimate of the total mass of the system, we can exploit the Y_{SZ}–M relation to obtain an estimate of the number n_{halos} of subhalos populating the Spiderweb complex.
First, the dispersion σ_{v} is converted into a dynamical mass estimate using the scaling relation calibrated by Saro et al.^{69} on a mock galaxy catalogue from the Millennium Simulation^{70}. This is rescaled to M_{500} assuming the conversion relation in Ragagnin et al.^{71} between masses at different overdensities. We then iterate over n_{halos} and compute the integrated SZ flux expected for each set of subhalos. For the sake of simplicity, we consider all the n_{halos} subcomponents to have equal mass ({M}_{500}^{{sigma }_{{rm{v}}}},/{n}_{{rm{halos}}}). Both the masses from our SZ analysis ({M}_{500}^{{rm{ALMA}}}) and from the velocitydispersion measurement ({M}_{500}^{{sigma }_{{rm{v}}}},/{n}_{{rm{halos}}}) are then used to obtain a measurement of the spherically integrated SZ flux Y_{SZ} by means of numerical integration of the pressure profiles over a volume of radius 5r_{500}. In our calculations, we consider the universal formulation by Arnaud et al.^{28} to describe the pressure distribution of the intracluster electrons. Finally, to derive the number of equalmass subhalos within the Spiderweb protocluster whose individual SZ fluxes would match the measured value ({Y}_{{rm{SZ}}}^{{rm{ALMA}}}={Y}_{{rm{SZ}}}left({M}_{500}^{{rm{ALMA}}}right)), we simply estimate the value n_{halos} for which the equality ({Y}_{{rm{SZ}}}left({M}_{500}^{{rm{ALMA}}}right)={Y}_{{rm{SZ}}}left({M}_{500}^{{sigma }_{{rm{v}}}},/{n}_{{rm{halos}}}right)) is satisfied (see Extended Data Fig. 2). We note that adopting an underlying pressure model different from the A10 profile used for producing Extended Data Fig. 2 is not causing a relevant variation in the total number of subcomponents n_{halos}, which are overall constrained to range between two to a maximum of four (except for the L15 8.5 case, resulting in n_{halos} ≃ 1). We further note that we made several attempts in using a physically motivated subhalo mass function (for example, the theoretical model provided by Giocoli et al.^{72}) instead of our simple equalmass distribution. However, we found that removing the strong constraint of having subhalos with the same masses makes the derivation of n_{halos} or of the subhalo mass parameters unconstrained, resulting in unstable and heavily degenerate results. We note, in any case, that all these considerations are derived a posteriori of the SZ modelling and, thus, do not affect the significance of the reported SZ detection.
This emerging multihalo picture is also consistent with the identification reported in the literature^{14,20} of double peaks in the velocity distribution. Nevertheless, we do not identify in the posterior distribution for the SZ model any notable secondary peaks that would be indicative of several pressure components in the ICM^{40,41} in the direct surroundings of the Spiderweb galaxy (see however the discussion below on the results of a multicomponent analysis). This might be caused by the chance lineofsight alignment of separate halos, but neither the spectroscopic information on the protocluster members nor the ICM constraints allow us to disentangle any distinct contributions from superimposed substructures. Similarly, any subhalos with similar masses would also be characterized by comparable pressure distributions, in turn causing their SZ signal to be barely distinguishable. Overall, the above result suggests that the SZ effect is tracing a minor portion of the larger Spiderweb structure in which the ICM has started building up and pressurizing, whereas the rest of system, extending over scales of tens of Mpc (refs. ^{11,12,20,25,35,49,65}) and tracing the region encompassed by the turnaround radius of the overdensity, has yet to undergo virialization. Any SZ signals associated with further subhalos or more extended structure are not constrained by the observations, which have limited sensitivity and may suffer from largescale interferometric filtering.
Multiple SZ components
The synchronous multiellipsoidal sampling^{73} typical of the main nested sampling algorithms—and, in particular, of dynesty (refs. ^{45,74}), the library used for our analysis—would naturally break into separate posterior modes in the presence of several peaks in the posterior density function. This would be the case, for instance, in the presence of several halos, resulting in distinct SZ components (as reported in, for example, refs. ^{40,41} in the case of merging cluster systems). However, as mentioned above, we do not find such evidence in the posterior probability distribution for any of the main modelling runs presented in this work. Nevertheless, we tested for the potential presence of any extra SZ features by performing a multicomponent analysis. In particular, we consider the same model description as in the singlehalo case above but introduce further SZ terms. As in the case of the radiosource modelling, to avoid label switching, we introduce an ordering condition on the centroid coordinates of the SZ components. This is applied first to the right ascension parameters and then to the declination direction, to test against any bias potentially introduced by the specific prior choice. Nevertheless, the results are found to be consistent between the two modelling sets.
Independently of the model used to describe the underlying pressure distribution, the sampler identifies a secondary SZ feature (27,.,{3}_{3.8}^{+2.{4}^{{primeprime} }}) southeast of the Spiderweb galaxy, corresponding to approximately 226 kpc at the protocluster redshift. This falls right at the r_{500} boundary of the main SZ component, implying that this secondary feature, if real, would be associated with a halo distinct from the one in which the Spiderweb galaxy is embedded. The actual existence of such a structure however cannot be firmly assessed. The images produced with the highresolution algorithm (see the ‘Sparse imaging’ section below) or after subtracting the radiosource model from the lowresolution set do not provide any clear evidence for any offcentre SZ structure. A lack of spatial correspondence is also noted with respect to the protocluster members, as the secondary SZ component cannot be clearly associated with any specific concentration of member galaxies, indicative of a distinct collapsed halo. The absence of a correspondence with protocluster galaxies further limits (if not excludes) the chances for the SZ component to be associated with the secondary velocity peak mentioned above. Above all this, the Bayesian evidence of the model comprising two SZ components improves on the onecomponent case by a factor of only ((log ,{Z}_{{{rm{n}}}_{{rm{s}}{rm{z}}}=2}log ,{Z}_{{{rm{n}}}_{{rm{s}}{rm{z}}}=1})lesssim 1.84), corresponding to an effective significance of σ_{eff} ≲ 1.90. As such, the ALMA+ACA data available at present are not able to support the unequivocal, statistically significant identification of a secondary pressure component.
In fact, increasing the flexibility of the SZ model beyond the twohalo scenario does not provide any effective improvements in the overall reconstruction. In particular, the inclusion of a third component induces the Bayesian evidence to degrade, with a reduction with respect to the twocomponent case of ((log ,{Z}_{{{rm{n}}}_{{rm{s}}{rm{z}}}=3}log ,{Z}_{{{rm{n}}}_{{rm{s}}{rm{z}}}=2})lesssim 1.67) and a limited improvement with respect to the reference singlehalo model ((log ,{Z}_{{{rm{n}}}_{{rm{s}}{rm{z}}}=3}log ,{Z}_{{{rm{n}}}_{{rm{s}}{rm{z}}}=1})lesssim 0.17) (which converts to σ_{eff} ≲ 0.60). At the same time, we observe that allowing for an ellipsoidal pressure distribution provides a marked improvement in the significance of the model (σ_{eff} ≲ 4.41). However, this concurrently makes the sampling converge to a hardly physical solution, with a planeofsky eccentricity (varepsilon =0.9{6}_{0.04}^{+0.02})—that is, corresponding to a planeofsky minor axis being only 4% of the respective major axis. This is mainly a consequence of the strong degeneracy between the mass (that is, the parameter controlling the overall amplitude and scale radius of the SZ signal) and the lineofsight extent of the ICM distribution, governed by the eccentricity parameter ε. In fact, SZ data alone cannot provide information on the lineofsight distribution of the optically thin ICM. We thus have to force the lineofsight scale radius to be equal to the geometric mean of the major and minor axes of the threedimensional ellipsoid, assumed for simplicity to lie on the plane of the sky.
Overall, the main consequence for the main SZ detection with data available at present is that a single spherically symmetric halo is sufficient to provide a statistically exhaustive description of the SZ footprint of the Spiderweb protocluster. The result of the elliptical modelling can only be interpreted as a marginal indication of an underlying morphological complexity, without however providing a conclusive and meaningful answer on the spatial properties of the SZ signal. At the same time, as demonstrated for the radiosource model, any resolvable asymmetry should be naturally traced by an ordered multicomponent model. Gaining a better understanding of the morphological properties of the forming ICM would require achieving improved quality from the observational side, in terms of both sensitivity and frequency coverage. Most importantly, though, the results reported above suggest that the robust identification of the SZ signal already with the simple spherical model provides only a lower limit to the actual significance of the detection and that this could only be enhanced when including in our models the description for any irregular and asymmetric features.
Dependence of the SZ significance on the number of pointlike components
First, we note that we do not observe a substantial variation in the fluxes of the compact components between the modelling runs with and without a SZ component. This provides a straightforward test of the robustness of our model reconstruction, in particular assuring against being driven in the SZ identification by the oversubtraction of the radio source. To properly assess whether the significance of the SZ signal is, however, dependent on the specific assumption on the number of pointlike components, we rerun the SZ modelling for the entire sample of model setups considered for finding the optimal set of compact sources. A summary of such a test for our reference model (A10 UP) is provided in Extended Data Fig. 3. We note that we consider here only the case n > 4, as this is found to represent the minimal condition for observing a sensible improvement in the imagespace residuals and for the sampler not to suffer from slow convergence.
The first outcome is that increasing the number of pointlike components beyond our reference model (n = 8) induces a drop in the Bayesian evidence, thus not justifying any further extension of the radio source to n ≥ 9 (second panel). This implies that, despite all the parameters remaining practically unvaried beyond the optimal set with n = 8, increasing the number of pointlike sources to n ≥ 9 makes the modelling incur data overfitting. On the other hand, before n = 6, the models exhibit a rapid increase in the overall significance in comparison with the null case, that is, the dataonly run with no model components except for the crossdata calibration parameters (first panel). Concurrently, for n ≤ 5, the right ascension and declination of the SZ centroid (fourth and fifth panels) are observed to roughly collapse on the values inferred for the secondary SZ feature identified at low significance in the multicomponent run discussed above. Corresponding to a region with low primarybeam amplitude (≲0.50, depending on the specific array and band), this is compensated by an abrupt increase in M_{500} when moving to n ≤ 5 (sixth panel), in turn resulting in a more extended (and thus more severely filtered) SZ signal. Despite being naively favoured by statistical reasoning on ({sigma }_{{rm{eff}}}^{{rm{sz}}}) (third panel), we note that the overall effective significance of the SZ models for n ≤ 5 is, however, substantially lower than the stable n ≥ 6 cases (first panel), thus limiting the validity of this reconstruction. Further, the cases n ≤ 5 correspond to radiosource models that substantially underfit the data and fail in describing the complex morphology of the extended emission from the Spiderweb galaxy (for this, we refer to Extended Data Fig. 4). As such, the identification of such an offset SZ feature cannot be reliably associated with the actual presence of any physical component and might be induced by spurious systematic features in the visibility data.
Despite overall supporting the effectiveness of our modelling and choice for the reference model, the above results clearly highlight a marginal level of variance introduced by the different assumptions on the number of compact sources. To take this into account and to limit any bias consequent to the choice of a specific model as reference, we thus decide to use Bayesian model averaging^{75,76} to generate modelmarginalized profiles when comparing our reconstruction with observations (see Fig. 1). In practice, we computed crossmodel posterior probability distributions by applying a weighted average to the original model posteriors, with the weight of each sampling point set equal to the respective Bayesian evidence. To fully account for any degeneracies between different model parameters, we generate radio and SZ uv models for each posterior sample and each pressure model, and then apply the Bayesian model averaging reduction to the resulting collection of uv profiles.
Systematics in the Bayesian analysis
We note that, because the dynamical mass estimates are computed assuming virial equilibrium for the entire structure, their values represent upper limits to the real mass of the Spiderweb protocluster. On the other hand, our SZderived mass estimates might be affected by nontrivial systematic biases associated with the assumptions used in our modelling.
First, the conversion between SZ signal and total mass is derived under hydrostatic equilibrium considerations. Given the disturbed nature of a galaxy protocluster, we can instead expect notable nonthermal contributions to the overall pressure support to the ICM—ranging from turbulent motion to dynamical effects caused by recent or continuing merger events. At the same time, hydrodynamical simulations^{77,78,79} show that the growth of the mass of a system resulting from a merger event does not correspond to a concomitant increase of the thermal SZ signal, owing to the temporal offset between mass evolution and gas thermalization. If the Spiderweb protocluster is observed while experiencing a merging process, the intrinsic SZ signal would thus be lower than expected from standard massobservable scaling relations.
Second, the reconstruction of the SZ signal relies on the assumption of selfsimilarity across cosmic time of the halo properties. Although this is observed^{80} to hold on large scales for galaxy clusters even up to z ≃ 1.9, we still lack an exhaustive description of the average thermodynamic properties for the ICM within protocluster complexes, as well as any observational information for the potential selfsimilar appearance of protocluster halos with the ones in their z ≲ 2 descendants. Also, for the sake of computational feasibility, we rely on the relatively strong assumption of spherically symmetric distribution for the electron pressure within the ICM. Nevertheless, the overall depth of the available data limits the possibility of obtaining better constraints on the SZ effect from the Spiderweb complex or achieving an improved separation of the signals from the radio source and the underlying SZ effect able to highlight any diffuse ICM halo with low surface brightness. Similarly, asymmetries in the morphology of the pressure distribution and, therefore, in the resulting SZ distribution, as well as to any secondary ICM components populating the Spiderweb complex, cannot be firmly assessed (see, however, ‘Multiple SZ components’ above for a discussion).
Third, given the limited information across the millimetre/submillimetre window with sufficient sensitivity to constrain the SZ spectrum, we are not able to disentangle any contribution to the overall SZ signal in the direction of the Spiderweb protocluster from the kinetic term. Similarly, we cannot constrain the relativistic correction to the SZ spectrum. Such effect is however generally subdominant at the virial temperature expected for such a lowmass system (for example, the relativistic correction to the thermal SZ effect at 2 keV is on the order of 1.2%; we refer to Mroczkowski et al.^{3} for a review of the different contributions to the SZ effect). Therefore, we assume the measured SZ signal to be entirely caused by its nonrelativistic thermal component.
Sparse imaging
The standard tools available at present for performing radiointerferometric imaging and deconvolution are not optimally suited to the joint reconstruction of the signal from an extended radio source superimposed over a diffuse SZ decrement. Crosscontamination may in fact cause both an underestimation of the flux of the former when using CLEANlike algorithms^{81}, and—at the same time—introduce a notable suppression of the underlying SZ footprint. Also, common approaches exploiting scale separation between the signature from the radio sources and the ICM in galaxy clusters^{82,83} would not be able to provide a robust characterization of the two signals, as their characteristic scales exhibit a broad overlap in the case of the Spiderweb protocluster.
Building on the extensive literature on compressed sensing^{84,85,86} and sparsitybased component separation^{87,88} and imaging^{89,90,91,92}, we thus developed an algorithm for taking full advantage of both the different and, in the case of the SZ effect, wellconstrained spectral behaviour of the measured signals, as well as the information on the different spatialcorrelation properties. In particular, we assume the total surface brightness I_{ν} in a given direction on the plane of sky (x, y) and at a given frequency ν to be described as
$${I}_{nu }(x,y)={I}_{{rm{RS}}}(x,y)cdot {(nu /{nu }_{0})}^{alpha (x,y)}+{I}_{{rm{SZ}}}(x,y)cdot {g}_{{rm{SZ}}}(nu )$$
(1)
Here I_{RS}(x, y) is the surface brightness of the radio source computed at the reference frequency ν_{0}, whereas α(x, y) is the corresponding spatially varying spectral index. The terms I_{SZ}(x, y) and g_{SZ}(ν) denote, instead, the amplitude in units of Compton y and spectral scaling for the thermal SZ effect^{2,3}. Because, as already discussed above, we are not including in our SZ model any corrections owing to relativistic terms, the spectral SZ scaling is determined solely by the frequency and therefore does not contribute to the overall model through any specific free parameter. To achieve a highfidelity reconstruction of the complex morphology of the structure of the radio jet from the Spiderweb galaxy, we model I_{RS}(x, y) and α(x, y) as a collection of pixels with angular scale matched to the longest baseline according to the Nyquist sampling theorem (given the maximum uv distance of u_{max} ≃ 1.3 Mλ, we set the pixel scale equal to 0.5/u_{max} ≃ 0.07″). On the other hand, the signaltonoise ratio of the SZ effect from the Spiderweb complex is not high enough to allow for adopting an analogous approach to constrain the SZ decrement. To exploit the expected largescale correlation, we hence describe I_{SZ}(x, y) through an isothermal, circular β model^{93}, with free centroid coordinates, amplitude normalization and core radius, whereas we arbitrarily fix β to a value of 1 (commonly used in ICM studies and adopted, for example, by the SPT collaboration as the base for their matchedfilter cluster template^{56,57,94}). We tested against potential biases introduced by the specific choice of β but did not observe a dependence of the solution on this parameter. Clearly, this model is less refined than the gNFW profile broadly used in the literature and in the lowresolution analysis presented above. Nevertheless, here, we are mostly interested in capturing the average properties of the faint SZ distribution, while providing a good estimate of the base level on top of which to measure the surface brightness of the radio structure and avoid its undersubtraction when imaging the SZ effect.
Using a nonparametric approach for reconstructing the radio jet has the clear advantage of allowing for constraining the signal on the finest scales accessible in the ALMA+ACA measurements. However, this causes the model to comprise (O(1{0}^{4})) free parameters, making the inference problem intractable through posterior sampling methods (for example, Markov chain Monte Carlo, nested sampling). At the same time, solving simultaneously for all the free parameters in the model of equation (1) through an optimization scheme would imply facing the shortcomings of a nonconvex likelihood function. We hence use an iterative, blockcoordinate optimization process^{95} with proximal point update^{96,97,98}. In this scheme, we alternatively solve at a time only for one set of parameters—namely, the radiosource amplitude, its spatially varying spectral index and the SZ centroid coordinates, normalization and core radius—through a boxconstrained limitedmemory Broyden–Fletcher–Goldfarb–Shanno algorithm^{99,100}. We use the implementation provided in the SCIPY^{101} package and assume all the standard values for controlling the convergence of the algorithm. The selection of the specific coordinate set to solve for is performed by means of the Gauss–Southwell rule^{102,103}, which introduces an ordering in the optimization blocks based on the norm of the Jacobian of the corresponding likelihood function.
When modelling the radiosource amplitude, we consider a LASSOlike regression^{104} step to promote sparsity,
$${L}_{{rm{R}}{rm{S}}}=,frac{1}{2}{chi }_{{rm{R}}{rm{S}}}^{2}+{lambda }_{{ell }_{1}}parallel {I}_{{rm{R}}{rm{S}}}{parallel }_{{ell }_{1}}$$
(2)
in which ({chi }_{{rm{RS}}}^{2}) is computed when fixing all the model parameters but I_{RS}, whereas the constant ({lambda }_{{{ell }}_{1}}) controls the strength of the ℓ_{1}norm regularization. To mitigate the bias and amplitude truncation^{105} introduced by the ℓ_{1} norm, we further iterate over the full blockcoordinate loop after reweighting the regularization constant according to the previous solution for the radio source. At every main step k of the optimization algorithm, we scale the regularization hyperparameter as
$${lambda }_{{{ell }}_{1}}^{(k)}=frac{{{epsilon }}^{(k1)}}{{{epsilon }}^{(k1)}+parallel {I}_{{rm{RS}}}^{(k1)}{parallel }_{{{ell }}_{1}}}{lambda }_{{{ell }}_{1}}^{(k1)}.$$
(3)
in which ϵ^{(k)} is an auxiliary coefficient aimed at ensuring numerical stability. We assume this to decrease with the number of reweighting iterations, ({{epsilon }}^{(k)}=max left({{epsilon }}^{(0)},1{0}^{k},{sigma }_{{rm{pix}}}right)), with ϵ^{(0)} given by the standard deviation per pixel of the signal as measured directly from the dirty image and σ is the RMS of the imagespace noise per pixel^{105}. Similarly, we can express the initial value ({lambda }_{{{ell }}_{1}}^{(0)}) in units of σ_{pix} as^{106}({lambda }_{{{ell }}_{1}}^{(0)}=n{sigma }_{{rm{pix}}}{Sigma }_{i}{w}_{i}), in which w_{i} are the weights of the individual visibilities and n is a userdefined quantity controlling the overall depth of the modelling process (that is, the optimization algorithm neglects any imagespace features with significance lower than about nσ_{pix}). Each blockcoordinate optimization loop is then run until the relative variation in the value of the likelihood function between two full sets of updates falls below an arbitrary tolerance of 10^{−8}. The same tolerance on the relative change of the likelihood function between two consecutive solutions is further used for interrupting the iteration over the reweighting steps.
For the starting values for the model parameters, we consider ({I}_{{rm{RS}}}^{(0)}) to be null everywhere and set ({alpha }_{{rm{RS}}}^{(0)}) to the spectral index of −0.70 typical for synchrotron emission. The coordinates of the SZ distribution are initially set on the phase centre of the Band 3 ALMA observations, roughly coinciding with the Spiderweb galaxy, whereas the initial amplitude and scale radius are initially set equal to zero and to an arbitrary scale of 30″, respectively. We tested for such choices by injecting random values for each parameter block at the beginning of the optimization run and found no substantial deviations in the final results.
In our analysis, we attempted to impose an extra smoothing to the radiosource amplitude through a total squared variation regularization^{107}, found to guarantee, in astronomical imaging, better results than the sole ℓ_{1}norm or standard isotropic total variation. However, we observed only a minor impact on the final solution and hence decided not to include any total squared variation term in the likelihood function. Similarly, the limited significance of the available data did not allow for successfully exploiting any waveletbased algorithms, generally demonstrated^{89,106,108,109,110} to provide an improved performance compared with simple imagebased techniques.
In our analysis, we adopt a conservative threshold n = 4. In our tests, a lower parameter showed improved subtraction capabilities but also deteriorated purity in the solution. As an internal crosscheck, we further tested against the recovered flux of the radio source, finding a maximumaposteriori estimate from the sparse modelling of 2.606 mJy against the value (2.60{2}_{0.007}^{+0.008},{rm{mJy}}) inferred using the nested sampling approach. Finally, the reconstructed source manifests a morphology that matches closely the signal measured by the VLA in Xband (Extended Data Fig. 5).
Insights from multiwavelength information
The combination of the ALMA+ACA data with the wealth of the multiwavelength information available on the Spiderweb galaxy allows us to propose several potential scenarios that could explain the intensity and morphology of the SZ signal and the associated offset with respect to the Spiderweb galaxy itself. Although mild spatial deviations between ICM halos and central galaxies are not a surprise even in moderately disturbed systems, the identification of such a displacement can inform our observational understanding of the dynamical state of a structure (see also ‘Comparison with simulated protoclusters’ below). Nevertheless, the available ALMA+ACA observations are not allowing to gain any insights on the physical and thermodynamic state of the forming ICM beyond the SZ detection—for example, the fraction of the total ICM reservoir that has already thermalized, the overall virialization state, the relative contribution of largescale accretion, mergers and AGN feedback in heating the assembling ICM, the multiscale morphology of the protoICM and its connection (or absence thereof) with the physical properties of protocluster galaxies. Obtaining a robust characterization of the intertwined processes determining the complex physical state of the forming ICM will thus necessarily require to further improve the observational coverage of the protocluster.
Given the inherently disturbed nature of the Spiderweb protocluster, the system could be characterized by marked contributions from the kinetic SZ effect in the direction of intracluster substructures with nonnegligible relative velocity with respect to the average bulk motion of the whole system. In particular, several studies^{15,19} report that the protocluster galaxies in the region southwest of the Spiderweb galaxy exhibit a relative blueshift compared with that of the central galaxy. This was interpreted as evidence for the presence of an approaching background feature, potentially tracing a filament infalling onto the protocluster core. Assuming that such structure hosts a nascent intracluster gas, this would thus be associated with a reduction of the total SZ flux owing to the infilling of the thermal SZ signal by a kinetic SZ contribution. An approaching ICM component with the same mass as that corresponding to the observed SZ feature would require, to completely cancel out its own thermal SZ effect, a lineofsight velocity of ({v}_{{rm{LoS}}}=,1,,92{8}_{699}^{+641},{rm{km}},{{rm{s}}}^{1}) when assuming the temperature estimate ({T}_{{rm{e}}}=2.0{0}_{0.40}^{+0.70},{rm{keV}}) from the Chandra Xray data^{23}. Alternatively, assuming the empirical mass–temperature scaling by Lovisari et al.^{111} and the resulting temperature of ({T}_{{rm{e}}}=0.8{5}_{0.25}^{+0.42},{rm{keV}}), the required velocity would be ({v}_{{rm{LoS}}}=79{6}_{266}^{+356},{rm{km}},{{rm{s}}}^{1}). Notably, the former velocity value roughly corresponds to the largest velocity reported by Miley et al.^{19}, whereas both estimates are consistent with the average velocity for both the Lyαemitting satellite members^{19} or the blueshifted components^{14} measured by the Spectrograph for INtegral Field Observations in the Near Infrared^{112} (SINFONI) on the Very Large Telescope (VLT). The same observations further show the presence of a small overdensity of redshifted members in the region immediately northeast of the Spiderweb galaxy. Similarly, the asymmetry in the SZ effect with respect to the central galaxy could be alternatively caused by an enhancement of the total SZ decrement owing to a negative kinetic contribution.
Although a notable contribution of the kinetic SZ effect can thus explain some of the multifrequency observations, we argue that this is probably not the only process responsible for the observed asymmetry. In fact, if the kinetic SZ effect were solely responsible for the skewed SZ decrement in the direction of the Spiderweb core, the Xray emission should not be affected. We do instead observe a marginal asymmetry between the extended halo and the Spiderweb galaxy also in the Chandra Xray data of the Spiderweb field, from which the bulk of the extended thermal emission is reported^{23,25} to have a mild shift towards the northeast with respect to the Xray signal from the central AGN (central panel of Extended Data Fig. 6; we refer to a forthcoming, dedicated paper (M. Lepore et al., manuscript in preparation) for an extended description of the Xray data processing). However, a further suppression of the overall Xray signal is observed in the region southwest of the main core—analogous to the SZ case—hence suggesting a local depletion in the forming intracluster halo. Recently, Anderson et al.^{8} and Carilli et al.^{9} reported direct evidence for the interaction of the extended radio jets with the Spiderweb environment, which could potentially result in the inflation of cavities in the diffuse protocluster gas. As commonly observed in the cores of local coolcore clusters, the presence of buoyant bubbles would induce a reduction of both the Xray and SZ signals over the limited regions covered by the cavities—provided that the bulk of their pressure is because of nonthermal components^{113}—and thus contribute to the observed asymmetry in the ICM emission.
Although potentially important, none of the processes discussed above seem to effectively explain why the protoICM centroid is offset with respect to the Spiderweb galaxy (Fig. 3) and the centres of the Lyα halo (right panel of Extended Data Fig. 6) and of the Hα nebula^{49} in which the Spiderweb galaxy is embedded. In a relatively relaxed system, the central galaxy should in fact sit deep in the centre of the gravitational well of the protocluster complex, along with most of the gas undergoing virialization, whereas the warm gas nebula is expected^{8,9,47,48} to be embedded in an extended halo of hot plasma. However, the lack of strong evidence for dominant SZ or Xray signatures from such a largescale halo, in combination with the numerous measurements^{33,63,114} showing that the central approximately 70 kpc in the Spiderweb complex hosts a massive reservoir of atomic and molecular gas, hints at a substantial, possibly predominant presence in the core region of warm and cold gas phases rather than ionized gas at the virial temperature. Emonts et al.^{115} further report the detection of a tail in the CO J = 1–0 line emission from the surroundings of the Spiderweb galaxy, interpreted as being associated with preferential accretion of cold gas along a filament. Notably, such emission is oriented along the same direction as the asymmetry in the SZ and Xray morphologies. At the same time, past observations provide evidence for relative motion of the Spiderweb galaxy within its own molecular halo^{33,114}, as well as a largescale systematic asymmetry in the velocity distribution of the protocluster members, denoting ordered flows of galaxies towards the central galaxy^{14,19}. Such motion could also naturally explain the difference observed in the radio morphology of the two opposite jets from the central radio galaxy^{8,9,48} (see Fig. 3). The western lobe is in fact found to be experiencing a stronger interaction with the environment than its eastern counterpart. Although recent studies^{8,9} based on VLA observations of the Spiderweb galaxy satisfactorily explained this behaviour as resulting from the inhomogeneity of the density field in the surroundings of the Spiderweb galaxy, it is not possible to exclude that the relative motion of the radio galaxy within the diffuse intracluster atmosphere is in fact contributing to the overall ram pressure enhancement. All these hints, taken together, suggest that the shift between the bulk distribution of the protoICM and the Spiderweb galaxy may be the net result of the merging activity experienced by the system, with a multiphase subhalo or filament infalling into the cold, compact core from the east–northeast and causing a local enhancement of the SZ and Xray signals through shock or adiabatic compression of the intracluster gas.
It is worth mentioning that the clear displacement between the protoICM halo and the Spiderweb galaxy might be indicative of the subdominant role of AGN feedback in providing the energy support required to heat the assembling intracluster gas. As reported by Tozzi et al.^{23}, the radio galaxy might be observed in the midst of the establishment of the radiomode feedback regime and potentially providing a subefficient heating mechanism. We can use the estimates of the nonthermal pressure and the volume of the radiobright regions provided by Carilli et al.^{9} to roughly evaluate the nonthermal energy provided by the radio jet, E_{nth} ≃ 6 × 10^{60} erg. Models of jet propagation through the ICM^{116} suggest that a comparable amount should already be deposited into the shocked intracluster gas. In fact, the jet energy E_{nth} is comparable with the ICM thermal energy within the approximately 60kpc sphere encapsulating the radio jets. This means that the mechanical energy might have a nonnegligible impact on the thermodynamics of the gas (and SZ flux) from the same region. However, to effectively affect the thermal energy of the ICM in a few 10^{13} M_{⊙} cluster up to r_{500}, the jet has to operate for an orderofmagnitude longer time. All in all, despite the broad evidence^{7,8,9,20,25,117} for a major crossinteraction between the central AGN and the assembling environment, the available data are still not allowing for inferring a conclusive answer on the relative contribution of AGN heating (versus, for example, gravitational process as infall and major mergers) to the overall thermal energy budget of the system. An investigation of the energetics and thermodynamics of the AGN–ICM system will thus be deferred to a future study.
Finally, it is important to note that the correspondence of the dominant features and, in particular, of the lack of thermal signal to the west of the Spiderweb galaxy in both the Xray and SZ data (see central panel of Extended Data Fig. 6) exclude that any potential undersubtraction of the signal from the extended radio source could have a dominant effect on our results (see, for example, ‘Nested posterior sampling’ above for a collection of tests on the model reconstruction). Nevertheless, improvements to our reconstruction would require a more extensive observational coverage, including deeper Xray and multifrequency SZ measurements, and higher angular resolution in particular for the latter. It is worth noting that the measured integrated Compton Y parameter—(Y( < 5{r}_{500})=(1.6{8}_{0.32}^{+0.35})times 1{0}^{6},{{rm{Mpc}}}^{2})—corresponds to a SZ flux density of about 0.6 mJy at 100 GHz. Such measurements are however barely in reach of currentgeneration subarcminuteresolution millimetrewave facilities.
Comparison with simulated protoclusters
The parametric modelling presented in the previous sections has the fundamental advantage of providing a straightforward and highly informative characterization of the average properties of the Spiderweb protocluster (for example, the halo mass). Nevertheless, owing to the many simplifying assumptions required to perform the model reconstruction using sensitivitylimited measurements, the analytic models used in our analysis are not able to describe the finer details of the intracluster gas. Furthermore, all the models for the ICM pressure profile available in the literature at present and explored in this work are calibrated on systems on a different mass or redshift range than the Spiderweb protocluster. In fact, we still lack an analytic description derived from observational considerations of the average thermodynamic distributions and, in particular, of the characteristic pressure profile for galaxy clusters and protocluster structures at z ≳ 2. Finally, because these models provide mean pressure profiles, they cannot capture the diversity and the possibly complex morphology of the pressure structure of the ICM in protoclusters.
The shortcomings of the parametric modelling motivates the use of hydrodynamical simulations to perform a more thorough exploration of the properties of the SZ signal measured in the ALMA+ACA observations of the Spiderweb system. We stress again that this analysis has only the scope of presenting a qualitative comparison with numerical predictions of objects evolved in a cosmological environment, without any pretence to identifying and analysing an exact numerical equivalent of the Spiderweb complex. The simulated protoclusters analysed in this work are taken from the ‘Dianoga’ simulation set^{18}. In particular, these comprise five regions that at z = 0 have sizes in the range 25–40 comoving Mpc h^{−1}, selected from a parent 1 h^{−1} Gpc^{3} darkmatteronly cosmological volume around five massive objects (with M_{500} at z = 0 in the range (5.636–8.811) × 10^{14} M_{⊙} h^{−1}). Such regions are resimulated with the Nbody hydrodynamical code OpenGadget3 at high resolution (with darkmatter and gasparticle masses of 3.4 × 10^{7} h^{−1} M_{⊙} and 6.1 × 10^{6} h^{−1} M_{⊙}, respectively). Besides hydrodynamics, which is described by an improved version^{118} of the smoothedparticle hydrodynamics, these simulations include radiative cooling, star formation, chemical enrichment and energy feedback from both supernovae and AGN (we refer to Planelles et al.^{17} and Bassini et al.^{18} for details on previous versions of the simulations used here and to a forthcoming paper by S. Borgani et al. (manuscript in preparation) for this specific sample).
The highresolution regions are large enough to resolve several halos, some of which—but not all—are progenitors of the main clusters selected at z = 0. For the purpose of this comparison, we select all systems with a mass M_{500} ≳ 1.3 × 10^{13} M_{⊙} at z = 2.16, for a total sample of 27 distinct objects. SZ maps are generated at the corresponding snapshot using the mapmaking software SMAC (ref. ^{119}) and projected onto the visibility plane. We note that, to simplify the comparison, the projected SZ maps are generated considering only the thermal component and they do not encompass any kinetic SZ term. Despite the obvious discrepancies associated with the specific differences in the morphological properties and dynamical states of the simulated protoclusters with respect to the Spiderweb complex, Fig. 2 shows that the mock uv profiles for simulated halos with mass in the same range of the parametric estimates (approximately 3 × 10^{13} M_{⊙}) are broadly consistent with the ALMA+ACA measurements.
To further gain visual insight into the simulated protoclusters, we produce synthetic maps by projecting the simulated ALMA+ACA visibilities back into the image space, consistent with the imagingreconstruction process of Extended Data Fig. 6. Three out of the full sample of simulated clusters are shown in Extended Data Fig. 7, with one specifically selected as exhibiting a SZ footprint resembling that observed for the Spiderweb protocluster, one as the most massive halo in the simulated set and another as exhibiting a displacement between the SZ peak and the halo centre analogous to what is observed between the SZ signal and the central galaxy in the Spiderweb complex. We checked that, for the entire sample, the physical displacement between the darkmatter halos and the corresponding dominant galaxies were never larger than about 6 kpc, more than an order of magnitude smaller than the observed SZgalaxy displacement in the Spiderweb protocluster. In the context of the discussion on such an offset, the halo centroid and galaxy for each of the simulated halos can thus be considered to be consistent. To introduce realistic noise and optimally reflect the noise properties of the available observations, we jackknife the original measurements by changing the sign of a randomly selected half of the visibilities. This noiselike realization is finally added to the idealized SZ data. Consistent with the uv result above, we can observe an overall agreement in terms of the average amplitude of the SZ surface brightness in the reconstructed images between the SZ signal measured in the ALMA+ACA data (Extended Data Fig. 6) and the estimate for the simulated lowmass system. We note that the selection of a specific projection axis has only a minor effect on the measured SZ amplitude. Over the whole sample, the change in both the peak and the average SZ signal associated with a different projection exhibits a RMS variation of about 10%, reaching a maximum of around 30% in only three cases, corresponding to objects that manifest morphological disturbances.
In any case, what is most important in such a comparison is the possibility of gaining a sensible understanding of the effect of radiointerferometric filtering on the observed distribution of the SZ signal. In the specific case shown in the bottom row of Extended Data Fig. 7, we can observe that the displacement between the SZ distribution and the centre of the corresponding darkmatter halo (in turn roughly coincident with the position of the most massive galaxy in the structure; see above) is exaggerated by the radiointerferometric suppression of the bulk, relatively symmetric SZ signal on large scales. The net result is an enhancement of local SZ asymmetries that finally drive the definition of the observed ICM morphology. In fact, the observed offset is not indicative of an absolute displacement between the ICM and the Spiderweb galaxy but rather of a positional mismatch between the peak of the pressure distribution and the radio galaxy, typical of disturbed systems.
It is clear that, to move this comparison beyond the sole qualitative analysis and to perform a robust crosscharacterization of the observed and simulated SZ signals, a broader exploration of the model parameter space in the available hydrodynamical simulations would be required. Such a detailed investigation will however be deferred to a future work.