© ESO 2022

1 Introduction

The formation and evolution of the large-scale structure (LSS) has long been a key target of astrophysical and cosmological investigations. Simulations predict the existence of a thin filamentary structure, the so-called cosmic web (Bond et al. 1996), connecting the knots in which galaxy clusters reside. Observationally, a detailed study of the cosmic web structure and evolution can be pursued by investigating the properties of the visible matter tracing the underlying dark matter (DM) distribution. This has promoted several investigations aiming for a census of cosmic baryons (Persic & Salucci 1992; Fukugita et al. 1998). Compared to observational estimates, theoretical studies predict that about half of the baryons in the Universe must indeed be undetected (Cen & Ostriker 1999; Davé et al. 2001). In particular, independently of the identification tools used to identify the cosmic web constituents (e.g. Libeskind et al. 2018), namely filaments, sheets, and voids, cosmological simulations consistently predict that a significant fraction (~40%) of the cosmic baryon budget is in the form of cool intergalactic medium and warm-hot intergalactic medium (WHIM; Cui et al. 2019; Martizzi et al. 2019), with most of it located in filaments (Tuominen et al. 2021). The properties of cosmic filament populations have also been studied from a statistical point of view in order to characterise their evolution and global properties in correlation with their tracers, namely gas and galaxies (Gheller et al. 2015, 2016; Cui et al. 2018; Martizzi et al. 2019). Recent investigations of various recent cosmological simulations, including the Magneticum Pathfinder suite we analyse here, show that statistically different populations exist: longer, thinner filaments that typically connect smaller structures, and shorter bridgelike filaments that are usually denser and connected to massive objects (Galárraga-Espinosa et al. 2020, 2021).

Observational reconstructions of the cosmic web are therefore tightly connected to the extensive search for the diffuse cold gas and WHIM in the filaments, which still largely eludes our detections (Bregman 2007). Reconstructions of the LSS through galaxy mapping have been successfully performed over the past decades (see e.g. Tempel et al. 2014; Malavasi et al. 2020). The gaseous content of the LSS has been observed mainly in the far-UV (see Bregman 2007; Péroux & Howk 2020, for a review), whereas in X-rays, direct observations of the gaseous cosmic web have proven to be very challenging (Nicastro et al. 2008). The main reason for this lies in the difficulty of observing the tenuous warm plasma residing in the filaments, whose low densities (<10⁻⁴cm⁻³) and temperatures (T ~ 10⁵-10⁷ K) make its X-ray emission very faint. Only recently has significant progress been made in this direction, starting with the first detection of the WHIM reported by Nicastro et al. (2018), which was obtained from the absorption lines of highly ionised oxygen (OVII) in spectra with a high signal-to-noise ratio of a quasar at z > 0.4, using the XMM-Newton Reflection Grating Spectrometer (RGS). Overall, direct observations of warm gas in filaments have mainly been limited to special cases, such as specific lines of sight in the direction of high-redshift quasars for the detection of WHIM in intervening absorbers (Nevalainen et al. 2015, 2019; Nicastro et al. 2018; Kovács et al. 2019) or regions in between close pairs of clusters or close to cluster outskirts, in which the WHIM emission could be enhanced (e.g. Briel & Henry 1995; Finoguenov et al. 2003; Durret et al. 2004). Longer gaseous filamentary structures up to tens of h⁻¹ cMpc have instead been investigated mostly through statistical approaches, such as stacking of observations of the thermal Sunyaev-Zeldovich (SZ) effect and of X-ray data (de Graaff et al. 2019; Tanimura et al. 2020a,b; Lim et al. 2020). Recently, Khabibullin & Churazov (2019) investigated the potential of eROSITA stacking approaches in characterising the WHIM. Specifically, accounting for the contribution of the resonantly scattered cosmic X-ray background and comparing the WHIM imprints in X-ray absorption and emission, they showed that the truly diffuse gas that fills cosmic filaments can be significantly characterised and better distinguished from clumps of denser matter.

A promising strategy for direct detections of the WHIM X-ray emission is to focus on specific regions encompassing multiple cluster systems (such as merging clusters and super-clusters) or on the prolongation of cluster outskirts (Fujita et al. 2008; Planck Collaboration Int. VIII 2013; Eckert et al. 2015; Bulbul et al. 2016; Hattori et al. 2017; Akamatsu et al. 2017; Parekh et al. 2017; Connor et al. 2018, 2019; Ghirardini et al. 2021). Several cluster binary systems in the local Universe have been thoroughly investigated to find evidence of cosmic diffuse gas in the bridges connecting the member clusters. Well-known case studies in the low-redshift Universe are A3556 and A3558 in the Shapley supercluster (z = 0.048, Mitsuishi et al. 2012; Ursino et al. 2015), RXC J1825.3+3026 (z ~ 0.065) and CIZA J1824.1+3029 (z ~ 0.071; Botteon et al. 2019), the A399-A401 system (at redshifts z ~ 0.0724 and z ~ 0.0737, respectively; Sakelliou & Ponman 2004; Fujita et al. 2008; see also thermal SZ observations by Bonjean et al. 2018), the merging pair A222-A223 (z ~ 0.21; Werner et al. 2008), or the multiple merging system Abell 1758 (z ~ 0.28; Durret et al. 2011).

The cluster system A3391-A3395 belongs to this group and has been widely studied in the literature, with special attention to the gaseous bridge between the pair members (Tittley & Henriksen 2001). This system comprises two main clusters, the northern cluster Abell 3391, and the southern cluster Abell 3395 (a double system itself; Reiprich & Böhringer 2002) at z ~ 0.05, separated by about 50′ (approximately 3Mpc) on the sky. Both systems have X-ray temperatures of T_X ~ 5 keV (Reiprich & Böhringer 2002; Vikhlinin et al. 2009). Early results by Reiprich & Böhringer (2002), combining ROSAT and ASCA data, indicated masses of M₅₀₀ ~ 5 × 10¹⁴M_⊙ for A3391 and M₅₀₀ ~ 8.8 × 10¹⁴M_⊙ for both A3395s and A3395n (assuming an Einstein-de Sitter cosmology and a Hubble constant H₀ = 50km s⁻¹Mpc⁻¹). Especially for A3395 and despite the uncertainties, these values were likely overestimated, probably because the system is in a merging phase. More recent estimates by Piffaretti et al. (2011), based on ROSAT X-ray luminosities, report lower masses of about M₅₀₀ ~ 2 × 10¹⁴M_⊙ for both A3391 and A3395.

From previous studies, debated conclusions have been drawn on the dynamical status of the pair, with different interpretations of the nature of the gas observed between the two main clusters (Tittley & Henriksen 2001; Planck Collaboration Int. VIII 2013; Sugawara et al. 2017; Alvarez et al. 2018). Sugawara et al. (2017) and Alvarez et al. (2018), for instance, reported the presence of hot gas in the bridge region, concluding that this could have been tidally stripped or heated by the interaction between A3391 and A3395. This is in contrast with earlier ASCA and ROSAT results from Tittley & Henriksen (2001), rather favouring a scenario in which the filament connecting the two systems is in fact longer and aligned lengthwise with the line of sight (l.o.s.). The geometry of the system with respect to the l.o.s. and the temperature of the interconnecting gas are key aspects to interpret the nature of the bridge, where the detection of cold-warm gas with T ≲ 1 keV can indicate the presence of true filament gas.

More recently, the A3391/95 system has also been targeted by the German-built X-ray telescope on board the Spectrum-Roentgen-Gamma (SRG) mission, the extended ROentgen Survey with an Imaging Telescope Array (eROSITA; Predehl et al. 2021), during the performance verification (PV) phase observational campaign (Reiprich et al. 2021). Thanks to the high sensitivity in the soft X-ray band, large field of view, and good spatial resolution (Merloni et al. 2012), these eROSITA observations allowed finding indications of emission from warm gas in the pair bridge, in addition to the known hotter gas emission. In the eROSITA observation, the bridge emission qualitatively spans ~3 Mpc (projected) (Reiprich et al. 2021), which cannot be entirely attributed to the known galaxy group (ESO 161-IG 006) that is located in between the two main clusters. Furthermore, the emission from hot and warm diffuse gas is detected well beyond the projected virial boundaries of the main clusters, which for the first time allows mapping a continuous warm-hot emission filament in the X-rays, spanning from north to south 4 degrees across the observed field. This WHIM filament extends north of A3391 and south of A3395 for a total projected length of ~15 Mpc at the median redshift of the system, and its presence is confirmed through the Planck SZ map and DECam (optical) galaxy density distribution as well. With this unprecedented X-ray detection of a WHIM filament in the observation of a single system reported by Reiprich et al. (2021), the new eROSITA PV observations also capture the formation of the LSS in the region in which the system is located. Several groups and clusters are discovered in its local environment (i.e. at the same redshift), and there are indications that the substructures move towards the A3391/95 system (see Veronica et al. 2022, Ramos-Ceja et al., in prep.).

With the primary goal of accompanying the observational results obtained by the eROSITA PV observations of the Abell 3391/95 system and its surroundings and to help interpreting them, we resort here to cosmological hydrodynamical simulations. Specifically, we focus on a simulated pair of galaxy clusters that was extracted from the Magneticum cosmological hydrodynamical simulation and shares several similarities with the observed A3391/95 system. The simulated analogue allows us to investigate the intrinsic properties of an A3391/95-like system and its formation and evolution. In particular, we aim at reconstructing the system assembly history by considering the main clusters and the additional structures in the surrounding environment. Furthermore, we can make use of the simulations to study the thermal and chemical properties of the diffuse gas component in the system at the selection redshift (z = 0.07) and throughout the cosmic evolution. This will allow us to explore the expected properties of the WHIM, mainly populating filaments, and of the bridge gas, in comparison to the ICM characteristics within the clusters. Especially, we can directly investigate the origins of the gas in the bridge region between the pair clusters to determine whether this can be filament gas or rather tidally stripped intracluster gas.

The paper is organised as follows. In Sect. 2 we present the simulation suite we used for the analysis, namely the Mag-neticum Pathfinder cosmological hydrodynamical simulations. Interested in multiple cluster systems in the local Universe, we focus on one simulation snapshot at low redshift (z = 0.07) and explore the cluster pair candidates in Sect. 3, depending on three-dimensional and projected separation, as well as on the mass and mass ratio of the members. From these, we select and introduce the cluster pair that best represents the theoretical analogue of the observed A3391/95 system (Sect. 3.1). In Sect. 4 we illustrate the main results of our study. The three-dimensional assembly of the pair system in the cosmic web, with its local environment, is presented in Sect. 4.1 together with the results on the origin of the diffuse gas residing in the bridge between the member clusters of the pair. In Sects. 4.2 and 4.3 we discuss the thermal and chemical properties and their time evolution of the diffuse gas in various thermal phases and spatial regions in and around the pair. In Sect. 4.4 we then investigate the groups that are aligned along large-scale filaments connected to the pair in order to find the signatures of their infall towards the overdense cosmic web knot in which the pair is located. Finally, we summarise and conclude in Sect. 6.

2 Magneticum Pathfinder simulations

We employed a large cosmological volume that was simulated to include a variety of physical processes to study the properties of multiple galaxy cluster systems and the evolution of the local cosmic web in which they are embedded. For this scope, we considered the Magneticum Pathfinder simulations¹, which are a set of advanced cosmological smoothed-particle hydrodynamics (SPH) simulations, comprising boxes of different volumes and resolution. The simulations were performed with the TreePM/SPH code P-Gadget3, an extended version of P-Gadget2 (Springel 2005). This code includes several improvements of the SPH formulation, such as the treatment of viscosity and artificial conduction, and the use of higher-order kernels (Dolag et al. 2005; Beck et al. 2016). The code also accounts for a large variety of physical processes describing the evolution of the baryonic components. These comprise radiative cooling and heating from a uniform time-dependent ultraviolet (UV) background (Haardt & Madau 2001), as well as a sub-resolution model for star formation(Springel & Hernquist 2003). Following detailed stellar evolution models, a description of chemical enrichment is also included as in Tornatore et al. (2004, 2007). Specifically, metals are produced from stellar sources depending on their mass and typical lifetimes, assuming an initial mass function (IMF) according to Chabrier (2003) and the mass-dependent lifetime function from Padovani & Matteucci (1993). The three main enrichment channels are supernovae type Ia (SNIa) and type II (SNII) and low- and intermediate-mass stars undergoing the asymptotic giant branch (AGB) phase. Eleven different chemical elements are explicitely traced (i.e. H, He, C, Ca, O, N, Ne, Mg, S, Si, and Fe) assuming stellar yields from van den Hoek & Groenewegen (1997), for AGB stars, Thielemann et al. (2003) for SNIa, and Woosley & Weaver (1995) for SNII. The radiative gas cooling depends on the local gas metallicity in a self-consistent way, as described in Wiersma et al. (2009). SNII also contribute to thermal and kinetic energy feedback by driving galactic winds, with a mass-loading rate proportional to the star formation rate (SFR) and a resulting wind velocity of v_w = 350 km s⁻¹ (Springel & Hernquist 2003). The simulations further account for black hole (BH) growth and gas accretion, powering energy feedback from active galactic nuclei (AGN), based on the implementation by Springel et al. (2005) and Di Matteo et al. (2005) and with modifications as in Fabjan et al. (2010). As shown in previous studies, the Magneticum simulations successfully reproduce many observed properties of cosmic structures, such as kinematical and morphological properties of galaxies (Teklu et al. 2017; Remus et al. 2017; Schulze et al. 2018; Remus & Forbes 2021), chemical properties of both galaxies and clusters (Dolag et al. 2017), and statistical properties of the AGN population at various redshifts (Hirschmann et al. 2014; Steinborn et al. 2016; Biffi et al. 2018a). At cluster scales, the Magneticum set has been employed to constrain the scatter in the observable X-ray luminosity-temperature relation through ICM velocity diagnostics (Biffi et al. 2013) as well as the pressure profile of the ICM (Gupta et al. 2017), and to explore the level of contamination of ICM emission due to the central AGN that is expected for synthetic eROSITA observations (Biffi et al. 2018a).

Because several well-known observed cluster pairs are at low redshift, and in particular, motivated by the eROSITA PV observation of the A3391-A3395 system, we focused on the simulation snapshot corresponding to redshift z = 0.07. The specific simulation we considered is one of the larger ones, namely the Magneticum "Box2" cosmological box at high resolution (hr). This box comprises a comoving volume of (352 h^-1 cMpc)³ and is resolved with 2 × 1584³ particles, corresponding to a mass resolution of m_DM = 6.9 × 10⁸ h⁻¹ M_⊙ and m_gas = 1.4 × 10⁸ h⁻¹ M_⊙ for dark matter and gas particles, respectively (see also Biffi et al. 2018a). Because the gas can spawn up to four stellar particles, which decreases its mass, or can accrete metals from chemical pollution, the mass of the gas particles can vary in time. At this resolution, the softening length for DM and gas particles is ε_DM.gas = 3.75 h⁻¹ kpc, and for stellar particles, it is ε_stars = 2 h⁻¹ kpc. The simulations assume a standard ACDM cosmological model with the Hubble parameter set to h = 0.704, the density parameters for matter, dark energy, and baryons equal to Ω_M = 0.272, Ω_Α = 0.728, and Ω_b = 0.0451, and σ₈ = 0.809 for the normalisation of the fluctuation amplitude at 8 Mpc (according to the seven-year results of the Wilkinson Microwave Anisotropy Probe, WMAP, Komatsu et al. 2011). The data of the simulation volume we used for the analysis, with several pre-computed halo properties and mock X-ray data (generated with the PHOX X-ray photon simulator, Biffi et al. 2012, Biffi et al. 2013), are publicly available at the Cosmological Web Portal² (Ragagnin et al. 2017).

In order to identify subhaloes, we employed the SUBFIND substructure finding algorithm (Springel et al. 2001; Dolag et al. 2009). The main halo detection is based on a standard friends-of-friends algorithm (FoF; Davis et al. 1985). Taking the presence of baryons into account as well (Dolag et al. 2009), SUB-FIND also identifies all self-bound substructures within the main haloes around local density peaks. In particular, the subhaloes of a cluster are defined as all the substructures identified within its virial radius (R_vir, computed with the spherical-overdensity approach using the top-hat model by Eke et al. 1996). SUBFIND also allows identifying the corresponding content of gravitationally bound DM, gas, stars, etc. for each substructure. The centre of each halo is identified with the position of the minimum of the potential well, considering the member particles.

For each main structure, characteristic quantities were computed, such as characteristic radii and enclosed masses. For data comparison purposes, we computed the R₅₀₀ radius, that is, the radius of the sphere encompassing an average density that is 500 times the critical density of the Universe (ρ_cr), and the mass within it, M₅₀₀. The estimate of R₂₀₀ used in the following corresponds to 1.5 × R₅₀₀ for simplicity, and we verified that this matches the value computed at the actual overdensity of Δ = 200 (with respect to ρ_cr) within a few percent.

For all metallicity values reported in this work, we refer to the solar abundance pattern by Anders & Grevesse (1989) as a reference for simplicity. Because the solar unit is purely used to normalise our estimates, this can be easily re-scaled for other solar reference values. For instance, in order to adopt the Asplund et al. (2009) solar pattern instead of Anders & Grevesse (1989), the Fe abundances reported in the following sections should be multiplied by a factor $Z_{\text{Fe,}\odot }^{\text{AG}89}/Z_{\text{Fe,}\odot }^{\text{Asp}109}=4.68\times {10}^{-5}/3.16\times {10}^{-5}~1.5$.

3 Galaxy cluster pair candidates

In the Magneticum Box2/hr box at z = 0.07, a total of 10429 haloes are identified with masses in the range 10¹³ < M₅₀₀[M_⊙] < 1.44 × 10¹⁵. In particular, we identify 448 objects with M₅₀₀ > 10¹⁴ M_⊙, from which we can therefore expect ~10⁵ candidate pairs of individually identified cluster-size haloes.

Motivated by the properties of observed pairs of galaxy clusters, for which a variety of observations have been dedicated to investigate their properties, by the stage of interaction, and by the possible detection of filamentary structures (bridges) that physically connect them, we further restricted the search to the pairs that are close in projection. We find that only 535(135) pairs appear closer than 10(5) Mpc in at least one projection (considering the three Cartesian axes x, y, and z as possible l.o.s. directions). By inspecting the three-dimensional spatial distribution of all the cluster-size haloes in the simulated volume, we nevertheless ind that only for 146, 88, and 40 unique pairs is the three-dimensional distance between the two haloes smaller than 20, 15, and 10 megaparsecs, respectively (see Table A.1). The number of close pairs decreases further when we impose additional constraints on the mass ratio between the two haloes or on the initial mass-range selection of the clusters. In particular, for haloes with masses in the range 1.5-3.5 × 10¹⁴ M_⊙, only 25 candidate pairs are closer than d_3D < 20 Mpc, of which ~30% are closer than d_3D < 10 Mpc.

These results for Box2/hr at z= 0.07 indicate that we expect roughly ~400 pairs per Gpc³ with cluster size members (M₅₀₀ > 10¹⁴M_⊙) separated by a physical distance smaller than 10 Mpc. In most of the cases (~70%), the mass ratio of the member clusters does not exceed 2, and in 20%, the members have comparable masses (M1/M2 ~ 1). In order to confirm this, we also investigated larger cosmological volumes (i.e. boxes with 0.9 h⁻¹ cGpc and 2.7 h⁻¹ cGpc per side) at similar resolution, which are part of the Magneticum simulation set. In the local Universe (z ~ 0), we consistently find ~300 cluster pairs per Gpc³ volume with a three-dimensional separation smaller than 10Mpc and member masses lower than M₅₀₀ > 10¹⁴M_⊙. By restricting the sample to member clusters with $1.5<M_{500}^{1,2}[{10}^{14}{M}_{\odot }]<3.5$, we expect to find ~50 close (d_3D < 10Mpc) pairs per Gpc³, 50% of which with a mass ratio M1/M2 ≲ 1.2-1.3.

In order to compare this with the specific case of A3391/95, we chose one of the seven close pairs at z= 0.07 (member separation of d_3D < 10 Mpc and mass ratio 1 < M1/M2 ≲ 1.2) for deeper investigation as a study case in the simulations that resembles the observed A3391-A3395 system fairly well. The specific selection was made based on a visual inspection of each candidate, taking the projected distances between the member clusters in the three main projections of the simulation box into consideration as well (i.e. we restricted the comparison to those along the x, y, and z Cartesian axes for simplicity). We remark nonetheless that each of the candidates with similar member masses, mass ratio, and physical separation might have been chosen in principle by inspecting all the possible l.o.s. and choosing the projection that resembles the observed system best.

3.1 A3391-A3395: Simulation analogue

The close galaxy cluster pair that we selected at z = 0.07 as a case study comprises two haloes with masses $M_{500}^1=1.96\times {10}^{14}{M}_{\odot }(GC1)$ (GC1) and $M_{500}^2=1.67\times {10}^{14}{M}_{\odot }\text{(GC2)}$, respectively, namely with a mass ratio of ~1.17. The two pair members are physically separated by d_3D = 4.54 Mpc and the projected distance is 2.6, 4.5, and 3.8 megaparsecs in the xy, xz, and yz projection planes, respectively. The physical distance between their centres approximately corresponds to three times the sum of their R₅₀₀ radii, $R_{500}^1=862\text{kpc}$ and $R_{500}^2=819$ kpc, respectively. This is equivalent to $d_{3\text{D}}~2\times (R_{200}^1+R_{200}^2)$ considering R₂₀₀ ≈ 1.5 × R₅₀₀. The top-hat virial radii (R_vir) are similar, both approximately 50% larger than R₂₀₀ (i.e. $R_{\text{vir}}^1=1.88$ Mpc and $R_{\text{vir}}^2=1.87\text{Mpc)}$ Mpc), and do not yet overlap. However, the spheres approximating 3 × R₂₀₀ are already overlapping physically and therefore also in projection.

We selected this pair, in particular, in the reference xy projection (i.e. with the l.o.s. along the z-axis of the cosmological box), motivated by the similarities with the A3391-A3395 system, as shown in Fig. 1. In the central panel of the figure, we report the eROSITA PV observation of the ~15 deg² region around A3391/95 (Reiprich et al. 2021) and label all the main systems. In the right panel of Fig. 1 we report the gas density map of the simulated analogue at z= 0.07, comprising a cubic volume of 20 h⁻¹ cMpc per side. Compared to the observed system, the GC1 and GC2 clusters have a broadly similar orientation and are surrounded by other systems aligned in a filamentary structure from top to bottom over 20 h⁻¹ cMpc in projection (approximately 26 physical Mpc at this redshift). Figure 1 left shows a simulation zoom onto the central region with 10 h⁻¹ cMpc per side, which approximately covers the same size as the eROSITA field. We label the two main clusters and their sizes, an infalling group-size halo (see Sect. 4.4), and mark in blue the interconnecting gas bridge. On the plane of the sky, they appear to be at a relatively similar distance (~3 Mpc) compared to the observed A3391/95, in which the projected 3 × R₂₀₀ extent of the two main systems also overlaps.

Compared to the simulations, the mass of the observed A3391/95 system reported by Piffaretti et al. (2011) is similar, although ~25% more massive overall, with M₅₀₀ = 2.16 × 10¹⁴ M_⊙ and M₅₀₀ = 2.4 × 10¹⁴ M⊙ for A3391 and A3395, respectively (Piffaretti et al. 2011; Alvarez et al. 2018). Similarly, we mainly considered the A3391/95 system as a cluster pair for the purpose of comparing them to the simulations, treating A3395 as a whole although it is a double-peaked merging cluster (Reiprich & Böhringer 2002; Reiprich et al. 2021). We note nonetheless that the southern GC2 cluster in the simulations shows a broad similarity with A3395. It hosts two massive substructures that recently merged (see further discussion in Sect. 5).

In Fig. 2 we show the evolution of the system in its local environment from z = 1 to z = 0.07 (from left to right). The maps show the gas surface density (upper row) and ideal X-ray emission in the [0.5-2] keV band (lower row) at three representative redshifts z = 1.04, z = 0.47, and z = 0.07, projected along the z-axis and centred on the pair centre of mass. In order to visualise the local environment around the system, these maps comprise a comoving volume of (20h⁻¹ cMpc)³. Overdense and higheremission regions are marked with lighter colours in the top panels and with blue to magenta colours in the bottom panels.

The evolution of this region shows how the assembly of structures proceeds along filaments by the accretion and merging of smaller structures and diffuse matter. The X-ray emission (Fig. 2, lower panels) shows that the main systems are aligned in a top-down direction, and the two pair members are enclosed by a common higher-emissivity envelope that traces the X-ray emitting diffuse gas. The left-most panel shows that at z ~ 1, the cores of the main cluster progenitors can already be distinguished, together with the other main haloes in the surroundings. Compared to the gas density maps, this allows us to better determine the principal filamentary structure that is defined by the main massive structures that align in the field and are embedded in the diffuse gaseous component. Minor gas filaments that are visible in the density maps are too faint to significantly emit in the X rays, even when the soft band is considered. The X-ray maps in Fig. 2 indicate a remarkable resemblance to the eROSITA PV observations of the A3391/95 field, where a ~ 15 Mpc gas emission filament has been detected based on eROSITA large FoV (~ 1 deg) and superior soft X-ray effective area. The observed filamentary gas emission extends north and south from the A3391/95 system and spans ~4 degrees in the sky. It connects other extended structures, namely groups and clusters, all at the same redshift (Reiprich et al. 2021).

Fig. 1 Comparison between the simulated galaxy pair candidate at z = 0.07 and the observed A3391/95 system. Left: zoom onto the projected gas density map of the simulated cluster pair. The image is 10 h−1 cMpc per side and projected for 10 h−1 cMpc along the l.o.s. We label the two clusters GC1 and GC2 and an infalling group (group B), and mark their extent. In blue we show the gas in the interconnecting bridge. Middle: eROSITA PV observation of a ~15 deg2 region around the multiple A3391/95 system, with the main clusters and Northern Clump marked (Reiprich et al. 2021). Here, 1 deg ≈ 3.9 Mpc, at the redshift of A3391 and for the cosmology adopted in Reiprich et al. (2021) (for comparison to simulations, we also report the comoving scale). Right: projected gas density map of the simulated cluster pair in its local environment. The map encloses a (20 h−1 cMpc)3 volume that spans ~26 Mpc (physical) in projection for the chosen redshift and cosmology.

3.2 Gas bridge between the pair clusters

The nature of the gas between the two main clusters in the observed A3391/95 system is debated. Following the common terminology, the gas filament between clusters in close pairs is typically referred to as a "bridge". Depending on the interaction stage of the two systems, this could be either stripped gas from the outer cluster atmospheres or actual filament gas.

Distinguishing pure filament gas from outer cluster atmospheres is very challenging per se in observations as well as in simulations. Therefore it is difficult to define the exact boundaries of the bridge. In our simulations, we defined the interconnecting bridge as a three-dimensional cylinder-like volume in between the two main clusters. For simplicity, we defined the bridge main axis as the three-dimensional line connecting the centres of GC1 and GC2. For the radius of the cylinder, we assumed ~660kpc. Along the main axis, the bridge was limited to the region outside of R₂₀₀ of each cluster. With the physical separation of the two clusters, this corresponds to ~2 Mpc along the spine direction. We verified that either taking a right cylinder (with flat ends) or modifying the main axis direction to account for asymmetries in the mass distribution of the clusters does not affect our conclusions.

We primarily considered the volume within a radius of 500 h⁻¹ ckpc from the spine, namely ~660kpc for the chosen redshift and cosmology, to focus on the core of the bridge and investigate the origin and properties of the gas. This radius approximately corresponds to ~0.8 × R₅₀₀ for both clusters. We adopted a similar definition as in Brüggen et al. (2021), who investigated the bridge of thermal gas between A3395 and A3391 that was recently observed with eROSITA, to search for signatures of radio synchrotron emission. They also defined the bridge as a cylindrical volume with a radius of 600 kpc and a length of 1.3 Mpc (see Sect. 5 for further discussion).

A schematic view of the system and the bridge, rotated to visualise the maximum physical distance between GC1 and GC2, is reported in Fig. 3. The bridge gas is marked in blue, as in Fig. 1. The gas mass in the bridge thus defined is roughly ~10% of the gas mass comprised within R₂₀₀ of either cluster in the pair.

Fig. 2 Projected maps of gas density and ideal X-ray emission in the [0.5-2] keV band. Each map is 20 h−1 cMpc per side, integrated for 20 h−1 cMpc along the l.o.s. (z-axis) and centred on the pair centre of mass at z = 0.07. The squares (10 h−1 cMpc per side) approximate the size of the eROSITA observation.

4 Results

4.1 Origin of the pair system

In Fig. 4 we show the evolution of the system at four redshifts between z ~ 0.47 and z = 0.07. The gas density maps, projected along the three major axes of the simulations (z, y, and x from left to right), comprise the local environment within a comoving volume of (20 h⁻¹ cMpc)³ centred on the pair centre of mass. In addition to the pair members GC1 and GC2, we also mark the position and R₂₀₀ extent of the most massive groups identified at z = 0.07 in the region (see Sect. 4.4). The maps allow following the assembly of the knot region in the last ~4 gigayears and show the complicated three-dimensional geometry connecting the multiple haloes and filaments around it.

Filamentary structures in the gas distribution are observed in all maps, with a complex three-dimensional geometry. Several gaseous filaments can be noticed even without significant group-size haloes located in them, such as in the regions left and right of the main halo chain in the z-projection (left panels). These can also be recognised in the other two projections. These minor filaments, however, are expected to be faint in the X rays, as shown in the lower panels in Fig. 2. The redshift evolution in Fig. 4 shows that the whole region, filaments and haloes, collapses to form the central cosmic web knot in which the cluster pair is finally located.

When we focus on the pair system in particular, simulations allow us to explicitly investigate the origin of the gas between the two member clusters. The evolution of the halo progenitors plotted in Fig. 4 shows that the pair members GC1 and GC2 have not yet merged by z = 0.07. Although R₂₀₀ of GC1 and GC2 almost overlap in the reference projection (see Fig. 1, left), we recall that the two clusters are indeed physically separated by a larger distance (i.e. by two times the sum of their R₂₀₀ radii, as described in Sect. 3.1), which might also be the case for the observed A3391/95 system (Tittley & Henriksen 2001; Reiprich et al. 2021), although it is difficult to assess the true physical separation from redshift arguments alone. In the simulations, we can explicitely visualise this from the additional projections along the y - and x -axes in Fig. 4 (last row, central and right panels).

We aim at assessing whether the gas in the bridge can be stripped gas from the outer atmosphere of the two clusters or is true ilament-like gas. To this end, we selected the gas that resides within GC1, GC2, and in the interconnecting bridge at z = 0. 07 and tracked it back in time.

Figure 5 shows the time evolution of the gas distances from the closest cluster progenitor, between redshift z = 0.47 and z = 0.07. For the two GC1 and GC2 clusters, we considered all the gas particles enclosed within their R₅₀₀, which are marked with dashed red and solid yellow lines, respectively. The gas bridge is selected in the region between GC1 and GC2, as deined in Sect. 3.2, and is marked with dotted blue lines.

The bulk (90-95%) of the gas in the bridge at z = 0.07 was never inside R₂₀₀ of either cluster progenitor, and the large majority (≳80%) was always beyond 2 × R₂₀₀ down to z ~ 0.25. Since R_vir ~ 1.5 × R₂₀₀, ≳80% of the bridge gas was essentially never within the virial radius (R_vir) of the main progenitors of the clusters. At lower redshifts, when the pair system approaches the final configuration, the gas bridge is increasingly confined between R₂₀₀ and 2 × R₂₀₀, where it is finally identified. From our analysis, we conclude that the diffuse gas in the bridge in this case is mostly ilament-like and comprises only a minor fraction of gas that was mixed in from the outer atmosphere of the two main clusters.

Compared to the bridge gas, Fig. 5 shows that the gas selected within R₅₀₀ in both clusters has been more smoothly accreted during their formation history. Below z ~ 0.25 it almost entirely resides already within 2 × R₂₀₀. Both GC1 and GC2 undergo significant mergers, which occur at z ≲ 0.5 for GC1 and more recently, at z ~ 0.2, for GC2. The latter, in particular, is clearly shown in the top right panel of Fig. 5, where the two clear components in the gas distance distribution of GC2 (solid yellow line) approach and merge between z ~ 0.25 and z ~ 0.14 (also visible in the gas density maps in Fig. 4). In the case of GC2, a significant fraction of the gas selected within R₅₀₀ at z = 0.07 comes from the merger, as ~55(34)% of it is still beyond R₂₀₀ at z ~ 0.25(0.18). In contrast, most (70%) of the gas that finally resides within R₅₀₀ of GC1 is already within its R₂₀₀ by z = 0.25, whereas only ~5% is farther than 2 × R₂₀₀. From this, we therefore also conclude that the gas within the R₅₀₀ of GC1 does not come from the northern group B, which finally enters the atmosphere (R_vir) of GC1 at z ~ 0.16 for the first time. Compared to the observed A3391/95 system, this is consistent with a picture in which the ICM of the Northern Clump, located at 3 × R₂₀₀ relative to A3391, has not yet been stripped and mixed with the ICM of A3391.

Figure 6 shows the spatial origin of the gas particles of the three components through their trajectories towards their final positions at z= 0.07. In the background, we report the projected gas density maps in the (20 h⁻¹ cMpc)³ comoving volume centred on the pair at z = 0.07 (same as Fig. 4 bottom row). The lines mark the trajectories of the gas in the different components, from their position at z ~ 2 till z= 0.07, all in comoving units. The colours distinguish the gas in GC1, GC2, and in the bridge between them (red, yellow, and blue, respectively; as in Fig. 5), with a colour gradient from early times (dark) until z = 0.07 (light). In order to facilitate visualisation, we only show the trajectories for 1000 randomly selected gas particles that sample each gas component. From the visual inspection of the trajectories, we note that the gas in the bridge (blue lines) collapses from directions almost orthogonal to the main accretion directions of the gas selected from the clusters (red and yellow lines).

Fig. 3 Schematic view of the simulated system and bridge, rotated to maximise the physical separation between GC1 and GC2.

4.2 Thermal and chemical properties of the diffuse gas

The geometrical configuration presented by this system, with a large physical separation relative to the R₂₀₀ of the member clusters and a mostly independent origin of the majority of gas in the bridge, makes it an optimal target for studying the thermo-chemical properties of the warm-hot gas populating cosmic filaments in comparison to ICM in clusters. Observationally, indications of the presence of colder gas in the interconnecting bridge of the A3391/95 system have been in fact found by eROSITA (Reiprich et al. 2021) as well. This encourages a deeper dedicated spectroscopic analysis of the warm gas in the bridge (Ota et al., in prep.).

In the simulated system, GC1 and GC2 are small-size clusters with mass-weighted average temperatures of $T_{\text{mw,500}}^1=2.60\text{keV}$ and $T_{\text{mw,500}}^2=2.93\text{keV}$, respectively. As often done in simulation analysis for a better comparison to X-ray observations, we can estimate a spectroscopic-like ICM temperature (T_s1), that is, an average temperature weighted by ρ²T^−0.75 (Mazzotta et al. 2004) instead of gas mass. In this case, we obtain $T_{\text{sl,500}}^1=2.03\text{keV}$ and $T_{\text{sl,500}}^2=2.64\text{keV}$ for GC1 and GC2, respectively.

In Fig. 7 we show a zoom-in mass-weighted temperature map of the system at z = 0.07, projected along the z-axis, in which we mark the two main clusters GC1 and GC2 and their R₅₀₀ radii. The map comprises a region of 7 h⁻¹ cMpc per side, centred on the cluster pair centre of mass. From the map we note that a radial gradient in temperature is present in both systems with higher temperatures (≳3.5 keV) in the cores. For GC2, we also note an asymmetry in the innermost region, with a more extended high-temperature core. This is related to the two main substructures, that recently merged (~3 Gyr ago; see also Fig. 5 and the related discussion). Within R₅₀₀, we note that the ICM in GC2 is overall hotter than that in GC1, although its M₅₀₀ is somewhat lower. In the region between the two clusters, the colour code shows colder (bluer) gas, which can be linked to the bridge connecting the two clusters (see also Reiprich et al. 2021). The figure further shows the colder region in the top right corner, which corresponds to the smaller group-size halo that approaches R₂₀₀ of the GC1 cluster (diamond; group B).

At z = 0.07, this sub-halo is found to be already gravitationally bound to GC1, and the merging started at z ~ 0.16.

In Fig. 8 we report the density-temperature phase-space diagram of the gas particles enclosed in the (7 h⁻¹ cMpc)³ region that we showed in Fig. 7. The gas density is here reported in terms of overdensity with respect to the mean background baryon density of the simulation, δ = ρ/ <ρ_gas >. The vertical line in the diagram corresponds to an overdensity of δ = 100, and the horizontal lines mark the temperature thresholds T = 10⁵ K and T = 10⁷ K. In this way, interesting gas phases can easily be separated depending on the gas temperature and density (see also Cen & Ostriker 2006; Ursino et al. 2010; Martizzi et al. 2019). We define as "hot" gas the phase with T > 10⁷ K (roughly corresponding to T > 1 keV). This is typically the gas within deep potential wells, such as those of galaxy clusters or groups, as well as in shocks. The warm gas is instead defined as the gas with 10⁵ < T[K] < 10⁷. In terms of density, the diffuse gas that is expected to fill the filaments and bridges interconnecting clusters in the large-scale structure is typically characterised by an overdensity of ~100 or lower. In particular, we refer here to the so-called WHIM as the gas with 10⁵ < T[K] < 10⁷ and δ < 100.

The colour code in Fig. 8 shows the gas Fe abundance, weighted by the gas mass in each temperature-density bin. The phase diagram shows that the hot phase and the dense and cold gas are typically characterised by the highest abundance level, with Fe abundances at about the solar value or higher. The low-density gas is characterised by a lower Fe abundance, Z_Fe,⊙ ~ 0.01-0.1 (see also the following sections), especially when the warm phase (10⁵ < T[K] < 10⁷) is considered.

We overplot contour lines (cyan curves) to mark the regions where ~25 and ~75% (inner and outer curves, respectively) of the total selected gas mass is located. They show that the majority of the gas mass (~75%) in the pair and its surroundings is in the warm and hot phase and at overdensities between a few and ~2500. The gas in smaller haloes (low temperatures and high densities) and the hot low-density gas associated with feedback processes only contribute ~10% of the gas mass in the whole region.

For comparison, in Fig. 9 we show instead the density-temperature diagram for the gas included in one of the two pair members (upper panel) and in the bridge (lower panel). As an example, we selected the GC1 cluster and show the gas enclosed within its R200 radius in Fig. 9a. The phase-space diagram is therefore only populated at overdensities higher than ~100. As indicated by the external contour level, 75-80% of the GC1 gas mass within R₂₀₀ has T ≳ 10⁷ K and median overdensities δ ~ 1000, as marked by the mass-weighted values reported in the diagram (dot-dashed lines). The temperature-density distribution also shows that the hot phase (T ≳ 10⁷ K) of the gas shows a gradient in the Fe abundance: at lower densities, the typical metallicity is also lower. This essentially corresponds to the X-ray emitting ICM, for which the radial abundance profile is indeed expected to be a decreasing function of the cluster-centric distance (see Mernier et al. 2018; Biffi et al. 2018b, for recent observational and numerical reviews), with higher values in the central regions (corresponding to higher overdensities) and lower values in the outskirts (at lower overdensities). The cluster potential well also comprises high-density cold gas, which is typically found in (or close to) star-forming regions, or within substructures, for example, galaxies. Typically, this is consequently characterised by high values of the metallicity because it has been more easily polluted by stellar sources. Although we only show the case of GC1 as an example for the pair clusters, we verified that the density-temperature diagram of the gas in GC2 presents very similar features.

Figure 9b shows that the majority of the gas in the core of the pair bridge (i.e. closer than <500 h⁻¹ ckpc from its spine and outside the R₂₀₀ of either cluster; Sect. 3.2) occupies the upper envelope of the distribution, as marked by the mass contours. Specifically, we find that in the simulated system, ~75% of the gas mass in the bridge is in the warm-hot phase, with a typical temperature of ~1 keV, a median overdensity of ~100 (marked by the dot-dashed lines), and an homogeneous iron abundance of Z_Fe ~ 0.1 Z_Fe.,⊙. The colder phase of the WHIM, with T ~ 10⁵-10⁶ K and δ ≲ 100, is not significant. This is consistent with recent simulation results by Galárraga-Espinosa et al. (2021), for instance, who reported hot gas within 1 Mpc from the spine of short filaments. At larger distances from the filament spine, they predicted the colder WHIM phases to dominate. By enlarging the bridge radius up to 1 h⁻¹ cMpc (approximately 1.3 Mpc in physical units), we also confirm a shift of the gas mean overdensity towards lower values δ ≲ 100. Finally, Fig. 9b shows that a minor fraction of the bridge gas mass is characterised by high values of overdensity and iron abundance, and low temperatures (T ≲ 10⁵ K). This halo gas is associated with galaxies, either comprised within small substructures or stripped from galaxies at earlier times. A small but non-negligible fraction of halo gas in short filaments was also reported by Galárraga-Espinosa et al. (2021).

Overall, in this system, the bridge gas closely reflects the chemo-energetic properties of the WHIM. The average temperature, density, and metallicity are clearly lower than the ICM within the two clusters. Combined with the analysis of its spatial origin (Sect. 4.1), this further supports the conclusion that this is filament-like gas and not ICM stripped by the interaction between the pair clusters.

Fig. 4 Evolution of the main haloes in the pair region in the three main projection axes (z, y, and x, from left to right). The maps show the projected gas density within a (20 h−1 cMpc)3 comoving volume centred on the final (at z = 0.07) pair centre of mass for different redshifts between z = 0.47 and z = 0.07 (from top to bottom). We mark GC1 and GC2 (red cross and yellow asterisk, respectively), the two northern groups A and B (green square and red diamond), and the additional main group-size haloes in the region (white triangles) and their progenitors. For each halo, circles approximate their R200 extent, except for group B at z ≲ 0.16, when it is identified as a substructure gravitationally bound to GC1. For GC1 and GC2, the dot-dashed circles at z = 0.07 correspond to 3 × R200.

Fig. 5 Evolution between z ~ 0.47 and 0.07 of the gas distance with respect to the closest progenitor of the two pair members in units of its R200 at the given redshift. Gas is selected at z = 0.07 within the bridge (dotted blue line) or within R500 of each cluster (dashed red line for GC1 and solid yellow line for GC2).

Fig. 6 Gas density maps in the three main projection axes (z, y, and x simulation axis, from left to right) at z = 0.07. Each map, centred on the cluster pair, encloses a cubic volume of 20 h−1 cMpc per side. We overplot in red, yellow, and blue the spatial trajectories of the gas particles that are finally selected within R500 of GC1 and GC2 and in the bridge at z = 0.07. Each trajectory shows a colour gradient with time, from dark at z ~ 2 to light at z= 0.07.

Fig. 7 Projected mass-weighted temperature map of the galaxy cluster pair at z = 0.07. The map is centred on the system centre of mass and is 7 h−1 cMpc per side, and the integration along the l.o.s. is also performed over 7 h−1 cMpc. The dotted white circles indicate the R500 radii of GC1 and GC2. The position of group B (diamond) is also marked.

Fig. 8 Density-temperature phase-space diagram of the gas in the (7 h−1 cMpc)3 region centred on the cluster pair at z = 0. 07, colour-coded by gas Fe abundance (in solar units with respect to Anders & Grevesse 1989). We overplot the contours (cyan lines) enclosing ~25 (inner) and ~75% (outer) of the selected gas mass.

Fig. 9 Same as Fig. 8, but only for the gas within R200 for one of the member clusters, namely GC1 (upper panel) and for the gas in the bridge (lower panel) between the pair clusters, at z = 0.07. Dot-dashed lines mark the average temperature and overdensity values (mass-weighted estimate) of the selected gas.

Metallicity distribution

Figure 10 shows the iron abundance distribution for the gas in the (7 h⁻¹ cMpc)³ volume around the pair at z = 0.07 for different thermal phases and spatial selections. In particular, we consider six cases: all the gas in the hot phase (T > 10⁷ K; black), the warm gas (10⁵ < T[K] < 10⁷ ; solid grey), the WHIM (10⁵ < T [K] < 10⁷ and δ < 100; dot-dashed grey), the X-ray-emitting ICM in GC1 and GC2 (with T > 0.3 keV, and located inside R₂₀₀ ; red and yellow respectively), and the gas selected in the bridge between them (blue; defined as in Sect. 4.1). In the figure we report the number distribution of iron abundance (upper panel) and oxygen-over-iron abundance ratio (lower panel), normalised to the total number of Fe- or O/Fe-rich particles in each gas subsample.

Considering the iron abundance (Fig. 10, top), we note that a broad double-peak feature is common to all curves, with differences in the gradient of the low-abundance tail and the relative amplitude of the peaks depending on the gas selection considered. The two peaks roughly correspond to abundances of Z_Fe ~ 10⁻³-10⁻²Z_Fe,⊙ and to solar values (~Z_Fe,⊙), respectively. For the bridge gas, the low-metallicity peak is further shifted to ~10⁻⁴ Z_Fe.⊙, with a more extended valley between the two peaks. Despite the differences, we note that the mass-weighted abundance of each sample is always ~10⁻¹ Z_Fe,⊙, namely between 0.1 Z_Fe,⊙ for the warmer phases and 0.2 Z_Fe,⊙ for the ICM/hot gas For the ICM in GC1 and GC2, the global estimate within R200 is dominated by the outer lower-density regions, where the bulk of the cluster mass resides (see also mass contours in Fig. 9). Compared to other selections, the whole hot-phase gas and the ICM residing within clusters typically show a similarly steeper low-metallicity tail and a more pronounced peak at low metallicity relative to the solar metallicity one. Overall, a large fraction of the ICM is significantly enriched, with abundances ≳10⁻² Z_Fe,⊙. Compared to the hotter gas phases, for the gas in the warm phase and the WHIM, the low-metallicity tail is shallower and has a higher normalisation (consistently with previous studies, e.g. Cen & Ostriker 2006), with a shallower peak at abundances ~10⁻³Z_Fe,⊙. This broad low-abundance component extends to even lower values, ~10⁻⁴Z_Fe,⊙, for the gas in the core of the pair bridge. Essentially, we can divide the different gas selections into two broad categories based on the trends shown in Fig. 10: the hotter gas within the main clusters, and the warmer, lower-density gas phase.

The warm diffuse gas outside the main clusters and in the bridge is typically poorly enriched. It still always comprises a highly enriched component originating from or still associated with galaxy-size haloes. This causes the distribution peak centred on solar abundances.

A broadly similar distribution is found for the oxygen abundance, with the main difference that both peaks are more pronounced and the peak at lower abundances is broader and shifted towards lower values. This can be investigated via the O/Fe abundance ratio, shown in the lower panel of Fig. 10. The oxygen-to-iron ratio distribution shows in all cases a very similar peak around solar values, originated by the corresponding solar abundance peaks present in both iron and oxygen distributions. Small differences are still observed between hot cluster gas and warmer bridge plasma in the left tail of the distribution. This is broadly flat for the warmer gas phases, whereas it presents a secondary peak for the hot ICM and a decrease towards Z_OFe ≲ 0.1 Z_OFe,⊙. The mass-weighted average O/Fe ratio is three times lower in the warm phase than within clusters, indicating that a larger fraction of warmer gas has a very low O/Fe abundance ratio. This is consistent with the fact that the vast majority of warm diffuse gas is expected to be located far from dense star-forming sites that are primarily responsible for the oxygen enrichment.

Fig. 10 Distributions of Fe abundance and O/Fe abundance ratio for the gas selected at z = 0.07 in different spatial regions and thermal phases. For each gas selection, the corresponding distribution is normalised to the total number of enriched (i.e. ZFe > 0 or ZOFe > 0) gas particles.

4.3 Evolution of the chemical properties of the gas

With the possibility of tracking the gas particles selected at z = 0. 07 back in time, we investigated the evolution of the chemical properties of the gas depending on its final location or thermal phase. This is shown in Fig. 11 for the gas selected at z = 0.07 within GC1 (top) and in the pair bridge (bottom), which is representative of the hotter and warmer phases, respectively. In the four panels of the figure we report the redshift evolution between z = 4.3 (purple curves) and z = 0.07 (red curves) of the iron abundance (left panels) and O/Fe ratio distributions (right panels) of these two gas subsamples.

From the evolution of the iron distribution, two general trends can be inferred: the gas that is finally selected in the hot phase and in clusters has undergone a stronger evolution in time than the gas that is finally located in the warm, low-density gas, such as the WHIM or the gas selected within the pair bridge. When we compare the two left panels in Fig. 11, we note that the shape of the iron abundance distribution changes in time in the GC1 case, whereas the gas that is finally selected in the bridge shows little variation, especially in the low-metallicity tail. The mass-weighted abundance increases with time in both gas subsamples, however, indicating that the mass fraction of highly enriched gas increases in both cases at low redshift. The solar abundance component noted in Fig. 10 seems to be present already from z ≲ 2-3 (blue and cyan lines). This can also be observed from the evolution of the O/Fe abundance ratio distribution in the right panels of Fig. 11, where a peak around solar values is always present below z ≲ 3 and is more prominent at high redshift, where it is driven by SNII oxygen enrichment. The latter also causes the higher values of the O/Fe abundance ratio that are reached at high redshifts, where the right tail is higher in normalisation and more extended. This confirms that part of the gas has already been enriched at early times (especially with oxygen) and survives in the diffuse form until z = 0.07 (see also Biffi et al. 2018a). Residual late-time enrichment drives the change in the O/Fe distribution at lower values, where a central peak around Z_OFe ~ 0.2 Z_OFe,⊙ builds up at low redshift, more prominently in the hot ICM within GC1 than in the gas ending up in the bridge. The gas selected within the bridge in contrast shows a shallower distribution at central O/Fe values, although some increase is also visible, and the peak at solar values still decreases in time relative to the lowest abundance ratios. This can still be attributed to late iron pollution due to long-lived SNIa.

While Fig. 11 shows the differences in the distribution shape at different redshifts, we also verified that despite preserving the general shape features, the normalisation of the curves increases with time when the distributions are computed in terms of gas mass (see Appendix C), and normalised at each redshift to the total mass of the selected gas. The mass fraction of iron- and oxygen-rich gas increases with time. This picture is summarised by the trends shown in Fig. 12, where we quantify the redshift evolution of the star-forming (i.e. with instantaneous SFR > 0) and enriched (Z_Fe > 0) gas fractions with respect to the total gas mass for each gas selection. The figures shows the enriched gas fraction increases with redshift in all cases we analysed, with ~90(50)% of the hot (warm) gas finally enriched by z = 0.07. Over the whole redshift range, 0.07 ≲ z ≲ 4, the gas in the hot phase and within clusters at z = 0.07 was systematically more enriched by a factor of ~1.7 compared to the gas that is finally in the warmer phases. The gas that will end up in the bridge shares a chemical evolution similar to the diffuse, warmer phases down to z ~ 0.5. At lower redshifts, however, the mass fraction of the (poorly) enriched component increases more significantly in the bridge gas than in the warm WHIM case.

In Fig. 12 we also show the time evolution of the highly enriched gas mass fraction, which we define as the mass fraction of the selected gas particles that at any given redshift have an iron abundance higher than the mass-weighted average value computed at z = 0.07 (i.e. with $Z_{\text{Fe}}>Z_{\text{Fe,mw}}^{z=0.07})$). Similarly, this also grows in time. This highly enriched component, nevertheless, was already present at high redshift in all selections, and typically increases by a factor of ~3 in mass between z ~ 2 and z = 0.07. Only for the ICM selected at z = 0.07 within the two clusters it increases more prominently (by a factor of ≳4). Although not shown here, an overall similar but shallower trend is observed for oxygen, with less pronounced differences in the various gas selections. The curves shown in the bottom inset of Fig. 12 show that the gas selected in the different cases at z = 0. 07 has always undergone a peak of star formation activity around z ~ 2. Regardless of where the gas finally resides, the mass fraction of tracked gas that was star forming is about a few percent at most at the peak of star formation. At redshifts z ≳ 3 in contrast, the star-forming gas adds up to less than 1% of the total mass of the selected gas, with slightly higher values for the gas that ends up within structures and negligible or zero star formation in the other cases. This trend is also observed at z ≲ 1, with no star-forming gas within clusters at z = 0.07, by definition.

Fig. 11 Redshift evolution of the iron abundance (left) and O/Fe abundance ratio (right) distributions for the ICM in one of the clusters (GC1; top) and the bridge gas (bottom), selected at z = 0.07 and tracked back in time. Vertical lines mark the MW average abundance at any corresponding redshift. The bridge histogram line for z = 4.23 is thinner to mark the very low statistics of metal-enriched gas particles (only ≲ 100, i.e. ~0.6% of the total).

Fig. 12 Evolution with redshift of enriched and star-forming mass fractions for the gas in different spatial regions and thermal phases, selected at z = 0. 07 and tracked back in time. Upper inset: mass fraction of Ferich gas (ZFe > 0; diamonds) and highly enriched gas ($(Z_{\text{Fe}}>Z_{\text{Fe,mw}}^{z=0.07};\text{triangles})$ triangles). Lower inset: redshift evolution of the star-forming (SFR > 0) gas mass fraction.

4.4 Additional structures in the local environment

Several other extended sources the sizes of groups and clusters of galaxies have been detected in the eROSITA PV observation of the A3391/95 field, some of them at the same redshift as the pair system (Reiprich et al. 2021; Ramos-Ceja et al., in prep.). Similarly, in the local environment surrounding the simulated cluster pair, we find additional massive structures, all comprised within the considered volume of 20 h⁻¹ cMpc per side (about 26 Mpc for the redshift and cosmology we considered). In the simulations, we identify four group-size haloes with masses in the range M₅₀₀ ~ 1.4-9.6 × 10¹³ M_⊙ in particular. Additionally, we considered group B, a structure with the size of a small galaxy group, which entered the atmosphere of GC1 at about z = 0.16 and reached its outskirts by z= 0.07 (when it was located at $d_{\text{3D}}~1.5\times R_{200}^1)$). At z ≲ 0.16, the halo is no longer identified as an independent group, but rather as a gravitationally bound substructure of the GC1 cluster by the SUBFIND substructure-finding algorithm. The position and R₂₀₀ extent of these five systems is also reported in the projected maps of Fig. 4 at different redshifts z ≲ 0.5. For the reference projection at z = 0.07 (bottom left panel in Fig. 4), the haloes are all labelled.

In Fig. 13 we instead report the redshift evolution of the total mass enclosed by R₅₀₀ (M₅₀₀) for the five groups and the main clusters GC1 and GC2 as well. All the haloes identified at z = 0.07 are tracked back in time, and the main progenitor (i.e. the most massive progenitor of each halo) is considered at every redshift up to z ~ 1. The figure shows an increases in the M₅₀₀ masses of GC1 and GC2 after the merging events around redshift z ~ 0.5-0.4 and z ~ 0.2 (as shown in Fig. 5).

Motivated by the eROSITA observation of the A3391/95 system, where some clumps that appear to be infalling towards the system have been discovered in the Northern and Southern Filaments (Reiprich et al. 2021), we concentrated in particular on two systems of these neighbouring structures to search for signatures of motion towards the pair. Considering the most massive groups (h159, group A, and group B; Fig. 13), we focused on the northern-most group A and the smaller group B, marked with the green square and red diamond in Fig. 4, respectively. Although it is massive (M₅₀₀ ~ 6.7 × 10¹³ M_⊙), we instead excluded the more complex group h159, which constitutes another close pair together with halo h9503. They are separated by d_3D ~ 3.5 h⁻¹ cMpc and have a mass ratio of almost 7:1. Through closer inspection, we also note that group h159 is a multiple system itself, comprising three massive gravitationally bound substructures of similar size (M_tot ~ 1.4–4.9 × 10¹³ M_⊙), all within its boundaries (≲1.5 × R₂₀₀). In contrast, groups A and B are both isolated single structures, but represent two extreme phases during the infall towards the GC1-GC2 system. Group A is still far away from the pair location, and group B already merges with one of the member clusters. By inspecting the evolution of the group B main progenitor, we find that its mass as a galaxy group is about M₅₀₀ ~ 3 × 10¹³ M_⊙ (more specifically, varying between M₅₀₀ ~ 2.5 ⊙ 10¹³ M_⊙ and M₅₀₀ ~ 3.6 × 10¹³ M_⊙) at redshift 1 ≳ z ≳ 0.16, that is, before it enters the GC1 atmosphere (see Fig. 13). Group A is an isolated, more massive group with M₅₀₀ ~ 5 × 10¹³ M_⊙ and is located at a distance of 7.7 h⁻¹ cMpc (i.e. ~10 physical megaparsecs) from the pair centre of mass, at z = 0.07. Figure 13 shows that groups A and B both have a rather smooth assembly history, without indications of significant merger events.

In the eROSITA A3391/95 field, the Northern Clump (MCXC J0621.7-5242) cluster located in the Northern Filament in particular shows several pieces of evidence that suggest a motion towards the A3391 cluster and a possible interaction with the filament gas in which it is embedded. These comprise, for instance, an emission enhancement towards the south, a tail of emission extending towards the north, and a pair of bent asymmetric radio lobes of its central wide-angle tails (WAT) radio galaxy as well (Brüggen et al. 2021; Veronica et al. 2022). We therefore investigated and compared the signatures of the movement of groups A and B towards the galaxy cluster pair at redshift z ≲ 1. This is quantified in Fig. 14 (upper insets) by the cosine of the angle (cos θ) between the group velocity and the direction pointing towards the final position of the pair at z = 0.07, where the cosmic web knot at the centre of the whole (20 h⁻¹ cMpc)³ region is forming. The figure shows the evolution of cos θ with redshift from z ~ 1 until z = 0.07. For comparison, we investigated the three directions defined by the positions at z = 0. 07 of the pair centre of mass (in blue), GC1 (in red), and GC2 (in yellow) separately. The figure shows that the cosine values for the gas and bulk velocity typically lie between -0.7 and -0.95 below z ≲ 1. This indicates that the velocity component along the infall direction is indeed dominant, and both groups are moving towards the cluster pair.

Small differences are nevertheless found in the two cases, with the smaller group B moving more specifically towards GC1 (cos θ is significantly lower for the red curve) and eventually merging with it. The value of cos θ for the direction towards the final position of GC2 (yellow curves) starts to significantly increase at z ≲ 0.5, marking the divergence of the group B trajectory. Instead, given the larger distance from the knot region, we find no significant differences for the three directions for group A, especially for redshift z < 0.5.

When we consider all the gas within R500 (dotted lines) instead of the total bound gas (dashed lines), then a more significant difference with respect to the halo bulk (DM dominated; solid lines) velocity is noted. The average gas velocity within R500 can be influenced by smooth accretion as well as by the infall of substructures onto the group. For group A, these accretion events are visible for instance around redshifts z ~ 0.75, 0.45, and 0.2, where the increase in | cos θ| is accompanied by a subsequent increase of the subhalo mass (Fig. 14, lower inset). In the case of group B, the accretion of the group itself onto the GC1 cluster is characterised by a corresponding decrease in the group total gas(DM) bound mass by a factor of ~3(1.6) between redshifts z ~ 0.16 and z = 0.07. The decrease in bound mass is the consequence of a stripping phenomenon that group B undergoes while penetrating the atmosphere of the more massive GC1 cluster.

The radial component of the velocity along the infall direction is reported in Fig. 15 across the redshift range 0.07 < z < 1. The evolution is shown as a function of radial distance from the pair, with larger radial distances corresponding to earlier times. As in Fig. 14, we investigated the three infall directions separately, and we focused in particular on the velocity of the bound components, with the bulk velocity essentially dominated by the dark matter. Comparing the two panels in the Fig. 15, we note different trends for the two haloes. Group A shows an oscillating behaviour of both gas and bulk velocities, with the oscillations amplified especially in the gas component. This translates into a relative difference between gas and bulk radial velocities that varies between -10% and +10% throughout the redshift-distance range we inspected. Since the feature is preserved in all the three directions considered (always defined with respect to the final positions of GC1, GC2 and pair centre of mass), it is indeed associated with the group itself. This reflects a sloshing trend of the gas component in the xy plane, where the movement along the filament mostly occurs. It is likely due to accretion processes of smaller substructures and diffuse gas from the filaments. This is further confirmed by the centre shift between gas centre of mass and halo centre, whose x and y components visibly present opposite oscillating trends (see Fig. B.1). While some similar signatures for sloshing are also shown at early times, the smaller group B presents a monotonic increase in relative difference between the gas and bulk radial velocities while approaching closer than ~2.6 h⁻¹ cMpc from GC1, namely below z ≲ 0.25. Even though the reference positions of the pair and its members considered in Fig. 15 are the final positions at z = 0.07, we can compute the actual distance between the centre of group B and the GC1 progenitor at z ~ 0.25, which corresponds to $d_{\text{3D}}~3.5h^{-1}\text{cMpc}\text{~}\text{4}\times R_{200}^1$. On the xy plane, the group B is found at a projected distance of $~3\times R_{200}^1$ at z ~ 0.25. Even though group B does not appear to be located within a main filament, this is the typical distance from clusters at which simulations predict that changes in the gas properties along filaments connected to them are observed, such as an increase in the gas temperature and a steeper radial gas density profile (i.e. at three to four times the cluster virial radius; e.g. Dolag et al. 2006). The configuration of group B roughly ~3 Gyr ago and the whole picture shown since then, that is, between z ~ 0.25 and z ~ 0.16, is consistent with the properties of the observed Northern Clump in the eROSITA A3391/95 field (Veronica et al. 2022), which is located at a similar projected distance of about three virial radii from the A3391 cluster (Reiprich et al. 2021; Veronica et al. 2022).

Fig. 13 Redshift evolution of the total mass M500 for all the clusters and groups identified at z = 0.07 in the (20 h−1 cMpc)3 volume around the cluster pair. At each redshift up to z ~ 1, we consider the main progenitor of the given halo identified at z = 0.07.

Fig. 14 Infall motion and mass evolution of groups A (left panel) and B (right panel) towards the cluster pair. Upper panels: redshift evolution of the alignment between radial direction towards the pair (centre of mass, blue; GC1, red; GC2; yellow) and subhalo velocity (as in the legend), quantified by the cosine of the angle between them. For Group B (right), the gas < R500 at z = 0.16 is tracked to lower redshifts after entering the GC1 virial boundary. Lower panels: redshift evolution of the subhalo DM (solid lines) and gas (dashed lines) bound mass.

Fig. 15 Radial bulk (solid lines) and bound-gas (dashed lines) velocities (upper panels) and their relative absolute difference (lower panels) for groups A and B (left and right, respectively) as a function of the radial distance from the system (centre of mass in blue, GC1 in red, and GC2 in yellow).

5 Discussion

We investigated the properties of a close pair of galaxy clusters showing several similarities with the observed A3391/95 system. Because the simulated analogue was extracted from a cosmological hydrodynamical simulation and not from an ad hoc simulation, some difference with respect to the observations is naturally present.

In our investigation, we identified the A3391/95 analogue by searching among close pairs of clusters. We note, however, that the A3395 cluster in the observed system is rather a double merging cluster itself (A3395n/A3395s; Reiprich & Böhringer 2002; Reiprich et al. 2021), even though in the literature it is sometimes treated as a whole but double-peaked system (Alvarez et al. 2018). The selected simulated pair is instead composed of two clusters, although the GC2 member also comprises two main substructures. These are two massive group-size haloes that underwent a major merging around 3-3.5 Gyr ago (z ~ 0.30.25, see Sect. 4.1) with a mass ratio of 1:1.6. At the time of merging, the two systems had masses of $M_{500}^1=7\times {10}^{13}{M}_{\odot }$ and $M_{500}^2=1.1\times {10}^{14}{M}_{\odot }$, respectively. Two emission peaks are still visible in the innermost regions of GC2 at z= 0.07, for instance in the bottom right panel of Fig. 2. Compared to A3395s and A3395n, the two structures in the simulated GC2 cluster are at a later stage of interaction. Furthermore, despite this broad similarity, we recall that the overall mass of the system is at least 20% lower than that of A3391/95, which is consistent with the smaller components we find in GC2 compared to A3395s and A3395n.

Because the masses of the simulated clusters are lower, we can also expect a smaller number of prominent filaments connected to them, as shown in Fig. 2. Numerical studies predict that most of the close pairs of clusters with separations of 5 h⁻¹ cMpc or less are typically connected by straight filaments, and less massive clusters are typically connected to fewer filaments (one or two instead of up to four or five) (Colberg et al. 2005).

In general, some of the structures in the observed A3391/95 system have no direct match in the simulations. This is the case, for instance, of the ESO 161-IG 006 galaxy group located in the bridge region. In the eROSITA observations, part of the X-ray emission in the bridge between A3391 and A3395 is associated with the ESO 161-IG 006 group, although it is not dominated by it. Reiprich et al. (2021) showed that the X-ray diffuse emission in the bridge qualitatively spans a horizontal scale of ~3 Mpc, which is larger than the region occupied by the ESO 161-IG 006 group. The definition of the bridge is somewhat uncertain, however. Brüggen et al. (2021) modelled the bridge in A3391/95 as a cylindrical volume with 600 kpc radius and 1.3 Mpc length, excluding the group emission. A similar slightly larger size (i.e. ~800kpc and 1.6Mpc for radius and length, respectively) was also assumed by Sugawara et al. (2017) in their modelling. We assumed a relatively similar radius (~660kpc) and a length of about 2 Mpc, which corresponds to the region outside the R₂₀₀ of both GC1 and GC2, given their physical separation. In the simulated analogue, none of the galaxy groups in the area lie between GC1 and GC2. This implies that the gas selected in the interconnecting region is mostly warm and diffuse, except for a minor fraction of colder gas located in galaxy-size haloes (see Sect. 4.2). The transition between outer cluster haloes and pure inter-cluster gas is difficult to determine, however. In our definition, gas and subhaloes in the outermost periphery of the clusters (R₂₀₀ < r < R_vir) lie within the geometrical boundaries of the bridge (see Sect. 3.2). The presence of clumps in the far outskirts of clusters has recently been investigated in cosmological simulations by Angelinelli et al. (2021), who reported that clumps around (and beyond) ~3 × R₅₀₀ can indeed be used to locate filaments connected to clusters and to study the thermodynamics of diffuse baryons before they are processed by the interaction with the ICM of the host cluster. A moderate (≲ 10% in mass) presence of dense halo gas associated with galaxies in short filaments was also found by Galárraga-Espinosa et al. (2021), who investigated the properties of filaments in different cosmological simulations, including the Magneticum set.

From the analysis done on the selected pair and given its geometrical configuration, we nonetheless conclude that the majority of the bridge gas has not been stripped from the member clusters due to their interaction. The two clusters GC1 and GC2 are separated by a physical distance that is ~3 × ($(R_{500}^1+R_{500}^2)$), and their virial boundaries (R_vir) do not overlap. The GC1 and GC2 spheres of influence (~3 × R₂₀₀) instead already overlap, and the two clusters move towards each other, as confirmed by the evolution shown in Fig. 4. We can expect that they will finally merge. Observationally, it is more difficult to determine the physical distance of the A3391/95 components along the l.o.s. In this respect, the indications of some cold gas component in the interconnecting bridge (Reiprich et al. 2021) can still be consistent with a scenario similar to the simulated case. That is, the filament has been forming in time during the LSS assembly and the two clusters can be in a pre-merging phase, in which no stripping or direct interaction between the cluster atmospheres has occurred so far. A longer interconnecting bridge of pure filament gas along the l.o.s. in case of a larger physical separation could be observation-ally confirmed by the presence of colder (T < 1 keV) WHIM gas in a detailed spectral analysis of the bridge region.

An early interaction phase of the two pair members was instead found in Planck Collaboration Int. VIII (2013) for a similar simulated system resembling A399-A401, with a significant part of the bridge gas stripped from the two clusters. The main difference with respect to the A3391/95-like system considered here was the three-dimensional distance of the two clusters relative to their size. In the case reported by Planck Collaboration Int. VIII (2013), the two clusters (more massive than GC1 and GC2) were separated by a smaller physical distance relative to their size, and the two virial radii overlapped in all projections.

Despite the several similarities between the simulated system and the observed A3391/95, we finally note that none of the groups found in the surroundings of the selected pair at z = 0.07 exactly matches the properties of the observed Northern Clump. Nonetheless, by selecting two interesting candidates, we were able to investigate the expected signatures of infall towards the overdense cosmic web knot in which the pair is located. Groups A and B served as extreme limits, both with similarities and differences with respect to the real Northern Clump. Group A is an isolated massive group moving along with the main northern filament towards the pair system, like the Northern Clump, but is located at a too large distance to experience the direct influence of any of the main clusters. Group B instead is a smaller group outside of the main north-south filament and is already merging with one of the clusters. However, group B showed a configuration similar to that of the Northern Clump about ~3 Gyr ago (z ~ 0.25), when it was at a distance of ~3 × R200 from GC1. We used this example to speculate that the features presented by the Northern Clump in the eROSITA observations are consistent with a radial motion of the system towards A3391 and the pair (Veronica et al. 2022). From the evolution of group B between z ~ 0.25 and the final snapshot (z = 0.07), we observe a slowing down of the gas compared to DM, and a stripping of the gas in the later merging stages (z ≲ 0.16) when group B finally enters the virial boundary of GC1.

6 Summary and conclusions

We presented results from the Magneticum cosmological hydrodynamical simulations on a case study resembling the binary cluster system in the observed eROSITA A3391/95 field. By investigating a cosmological comoving volume of (352 h⁻¹ cMpc)³ at high resolution, we considered all possible close pairs of cluster-size haloes and selected a candidate at z = 0. 07 presenting several global similarities with the observed A3391/95 system. Specifically, by restricting ourselves to all clusters with M₅₀₀ > 10¹⁴M_⊙ in the local Universe, we expect to find 300-400 pairs with a physical separation smaller than 10 Mpc per (Gpc)³. This was confirmed when larger cosmo-logical volumes of the Magneticum suite were inspected. For our study, we focused in particular on the large-scale region of 20 h⁻¹ cMpc per side around the pair, approximating a physical size of ~26Mpc for the given redshift and cosmology. The selected simulated system is located in a node of the cosmic web that is is assembling, with filaments building up in time and substructures moving along with them towards the cluster pair. The main clusters in the pair are slightly less massive than the members of the observed A3391/95 system, althought the mass ratio (M1/M2 ~ 1.2) and configuration are somewhat similar. The projected distance between the pair clusters is indeed similar (2.6 Mpc in the reference projection), and the regions enclosed by 3 × R₂₀₀ overlap. The simulated system is one of the about seven close pairs found at z=0.07, with member masses of $1.5<M_{500}^{1,2}[{10}^{14}{M}_{\odot }]<3.5$ and a three-dimensional separation of d_3D ~ 4.5 Mpc. Typically, our simulations predict that ~25 close (d_3D ≲ 10 Mpc) pairs per (Gpc)³ are found whose members have masses in this range and a mass ratio of M1/M2 ≲ 1.2.

Based on the encouraging similarity to the A3391/95 system in the eROSITA field, we then explored the properties and origin of the pair and of the diffuse gas selected at z = 0.07 in different thermal phases and spatial regions. We summarise the main results we obtained below.

• The gas located in the pair bridge originates from distances of a few to ten comoving megaparsecs away with respect to its final position. Furthermore, its accretion trajectories are almost perpendicular to the filamentary directions through which the gas found within the main clusters has been accreted, marking a distinctly different origin.

• Most (~90%) of the gas in the bridge was never inside the R200 radius of either cluster progenitor, and a substantial fraction (≳80%) was always beyond 2 × R₂₀₀ at z ≳ 0.25. The simulation study thus predicts that the majority of the bridge gas has not been stripped from the member clusters despite their closeness.

• The bridge gas is characterised by a typical temperature of ~1 keV and by an overdensity of about δ ≲ 100 with respect to the mean baryon density of the Universe. Lower densities are expected at larger distances from the bridge spine (≳ 1 Mpc) or in longer filaments, where the colder WHIM phase is expected to become more significant (Galárraga-Espinosa et al. 2021);

• In the whole pair region, the enrichment level of the warm diffuse gas and WHIM is fairly homogeneous and lower than for the hot ICM in the clusters (by a factor of ~2 for the iron abundance) or the cold-dense gas in galaxies and star-forming regions. Nonetheless, a highly enriched component is present at z ≲ 2 in all gas selections, and is mostly related to early enrichment within galaxies and star formation sites. This component, roughly ~ 10-20% of the gas mass depending on the specific selection, also persists in the bridge, within small clumps (galaxies), or in the diffuse form after being stripped from those haloes.

• Similarly to the A3391/95 field observed by eROSITA, we also find additional group-size objects in the surroundings of the pair in the simulation (see also Ramos-Ceja et al., in prep.), mostly aligned along a prominent north-south filament structure spanning about ~26Mpc (20 h⁻¹ cMpc) in projection (Figs. 2-4). We identify in particular five groups in the mass range M₅₀₀ ~ 1.4-9.6 × 10¹³ M_⊙ at z = 0.07.

• The additional systems, along with filaments, trace the assembly of the LSS around the overdense node in which the pair is finally located. By inspecting the properties and trajectories of two group-size haloes in particular, groups A and B, we find clear signatures of their motion towards the pair, with strong alignment of the halo bulk (gas and dark matter) velocity with the infall direction towards the pair.

• The farther and isolated group A is characterised by a visible sloshing of the gas component with respect to the DM component, while no significant mergers characterise its mass accretion history, except for accretion of gas and smaller clumps (galaxies) through the filament. This is marked by the oscillating evolution of gas and DM infall velocity components as well as by the centre shift between halo centre and gas centre of mass. - The smaller group B is identified as a substructure that is gravitationally bound to the pair cluster GC1 at redshift z ≲ 0.16 and is located in its outskirts at z = 0.07. Tracing the trajectory of the group towards the pair, we find evidence of the influence of the GC1 cluster. At a projected distance of about $3\times R_{200}^1$ from GC1, at about z ~ 0.25, the gas radial velocity starts to systematically decrease compared to the DM velocity. While entering the outer atmosphere of GC1 (z ≲ 0.16), it undergoes gas stripping, as marked by the decrease in bound gas mass by a factor ~3. The configuration corresponding to 3 Gyr ago (z ~ 0.25) shares several similarities with the case of the real Northern Clump relative to A3391, supporting the picture of its actual infalling motion towards the A3391/95 system (see Veronica et al. 2022).

The eROSITA superior soft response and large FoV has allowed the first direct detection of the faint X-ray emission from the diffuse gas in the bridge and filaments connected to A3391/95, drawing a consistent picture of the LSS that compares well to theoretical predictions. Our conclusions from simulations further support that the geometrical configuration of the A3391/95 system, its pre-merger phase and physical separation, together with the large-scale structure of the field in which it is located, provide an optimal target for unveiling the elusive warm-hot gas that populates the cosmic web filaments and for characterising its physical properties.

The eROSITA all-sky survey mode will finally enable a more statistical approach by providing a significantly larger sample of candidate multiple cluster systems and gas emission filaments. Combined with and compared to theoretical predictions from simulations, where samples of binary or multiple cluster systems in various configurations can be statistically investigated as well, this will greatly benefit the baryon census in the Universe, especially through the detailed characterisation of the chemo-energetic properties of the diffuse pristine gas.

Acknowledgements

The authors would like to thank the anonymous referee for constructive comments and suggestions that contributed to improve and clarify the presentation of this work. V.B. kindly acknowledges A. Merloni for providing comments on the manuscript and inspiring discussions with U. Maio on cosmic gas chemistry. This research was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - 415510302, and also partially supported by the Excellence Cluster ORIGINS, which is funded by the DFG under Germany's Excellence Strategy - EXC-2094-390783311. Part of this work has been funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - 450861021. The Magneticum Pathfinder simulations have been performed at the Leibniz-Rechenzentrum with CPU time assigned to the projects pr86re and pr83li. K.D. acknowledges support through the COMPLEX project from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program grant agreement ERC-2019-AdG 860744. This work was supported in part by the Fund for the Promotion of Joint International Research, JSPS KAKENHI Grant Number 16KK0101. This work is based on data from eROSITA, the soft X-ray instrument aboard SRG, a joint Russian-German science mission supported by the Russian Space Agency (Roskosmos), in the interests of the Russian Academy of Sciences represented by its Space Research Institute (IKI), and the Deutsches Zentrum für Luft- und Raumfahrt (DLR). The SRG spacecraft was built by Lavochkin Association (NPOL) and its subcontractors, and is operated by NPOL with support from the Max Planck Institute for Extraterrestrial Physics (MPE). The development and construction of the eROSITA X-ray instrument was led by MPE, with contributions from the Dr. Karl Remeis Observatory Bamberg and ECAP (FAU Erlangen-Nuernberg), the University of Hamburg Observatory, the Leibniz Institute for Astrophysics Potsdam (AIP), and the Institute for Astronomy and Astrophysics of the University of Tübingen, with the support of DLR and the Max Planck Society. The Argelander Institute for Astronomy of the University of Bonn and the Ludwig Maximilians Universität Munich also participated in the science preparation for eROSITA.

Appendix A Cluster pairs in the Magneticum simulations

In Table A.1 we report the number of cluster pairs found in the Magneticum simulation (Box2/hr) at z = 0.07 for different separations and mass ratios³. We considered unique pairs of clusters with $M_{500}^{1,2}>{10}^{14}{M}_{\odot }$ (448 haloes, upper half in the table) as well as pairs in which the two haloes have masses in the range $1.5<M_{500}^{1,2}[{10}^{14}{M}_{\odot }]<3.5$ (165 haloes; lower half of the Table).

Appendix B Clumps

In Fig. B.1 we report the centre shift evolution for the groups A (left) and B (right) discussed in Sec. 4.4. The centre-shift is computed as the three-dimensional difference between the centre of the halo (defined as the position of the minimum of the potential well) and the bound-gas centre of mass, and is given in units of the halo R₅₀₀ radius at any given redshift between z ~ 1 and z ~ 0.07. In addition to the modulus of the shift, we also show the three components along the major simulation axes separately (as in the legend). Comparing the two panels, we note that the amplitude of the shift is overall larger in the smaller infalling group B relative to its size. Considering the separate shift components, an oscillating opposite trend in the xy plane is observed in both cases, which indicates a sloshing feature of the gas component relative to the DM.

Fig. B.1 Redshift evolution of the centre shift between the (bound) gas centre of mass and the halo centre in units of R500 for groups A (left) and B (right) analysed in Sec. 4.4. We report the modulus and each separate component, as in the legend.

Appendix C Metallicity distribution

In Fig. C.1 we show the distributions of iron abundance and O/Fe abundance ratio for the hot gas (T > 10⁷K) and WHIM (10⁵ < T [K] < 10⁷ and δ < 100) in the (7 h⁻¹ cMpc)³ region around the selected pair (see Sec. 4.2). Gas particles in the two cases are selected at z = 0.07 and tracked back in time. By comparison with Fig. 11, we note that the hot and WHIM gas present overall similar trends to the hot ICM in clusters and the bridge gas, respectively, as discussed in Secs. 4.2-4.3. Furthermore, the mass fraction distribution of the iron abundance is reported (left-most panels) to demonstrate that the normalisation of the distributions for a given gas selection indeed increases with time, indicating as expected that the mass fraction of the enriched gas grows from high to low redshift (as in Fig. 12).

Fig. C.1 Iron abundance and O/Fe ratio distribution for the hot gas and WHIM (upper and lower rows, respectively) for different redshifts as in the legend. The left panels show the mass fraction distribution of the hot-gas iron abundance for comparison. At each redshift, the mass fractions are normalised to the total gas mass of the selected subsample.

References

All Tables

All Figures

Fig. 1 Comparison between the simulated galaxy pair candidate at z = 0.07 and the observed A3391/95 system. Left: zoom onto the projected gas density map of the simulated cluster pair. The image is 10 h−1 cMpc per side and projected for 10 h−1 cMpc along the l.o.s. We label the two clusters GC1 and GC2 and an infalling group (group B), and mark their extent. In blue we show the gas in the interconnecting bridge. Middle: eROSITA PV observation of a ~15 deg2 region around the multiple A3391/95 system, with the main clusters and Northern Clump marked (Reiprich et al. 2021). Here, 1 deg ≈ 3.9 Mpc, at the redshift of A3391 and for the cosmology adopted in Reiprich et al. (2021) (for comparison to simulations, we also report the comoving scale). Right: projected gas density map of the simulated cluster pair in its local environment. The map encloses a (20 h−1 cMpc)3 volume that spans ~26 Mpc (physical) in projection for the chosen redshift and cosmology.

In the text

	Fig. 2 Projected maps of gas density and ideal X-ray emission in the [0.5-2] keV band. Each map is 20 h−1 cMpc per side, integrated for 20 h−1 cMpc along the l.o.s. (z-axis) and centred on the pair centre of mass at z = 0.07. The squares (10 h−1 cMpc per side) approximate the size of the eROSITA observation.
In the text

	Fig. 3 Schematic view of the simulated system and bridge, rotated to maximise the physical separation between GC1 and GC2.
In the text

Fig. 4 Evolution of the main haloes in the pair region in the three main projection axes (z, y, and x, from left to right). The maps show the projected gas density within a (20 h−1 cMpc)3 comoving volume centred on the final (at z = 0.07) pair centre of mass for different redshifts between z = 0.47 and z = 0.07 (from top to bottom). We mark GC1 and GC2 (red cross and yellow asterisk, respectively), the two northern groups A and B (green square and red diamond), and the additional main group-size haloes in the region (white triangles) and their progenitors. For each halo, circles approximate their R200 extent, except for group B at z ≲ 0.16, when it is identified as a substructure gravitationally bound to GC1. For GC1 and GC2, the dot-dashed circles at z = 0.07 correspond to 3 × R200.

In the text

	Fig. 5 Evolution between z ~ 0.47 and 0.07 of the gas distance with respect to the closest progenitor of the two pair members in units of its R200 at the given redshift. Gas is selected at z = 0.07 within the bridge (dotted blue line) or within R500 of each cluster (dashed red line for GC1 and solid yellow line for GC2).
In the text

Fig. 6 Gas density maps in the three main projection axes (z, y, and x simulation axis, from left to right) at z = 0.07. Each map, centred on the cluster pair, encloses a cubic volume of 20 h−1 cMpc per side. We overplot in red, yellow, and blue the spatial trajectories of the gas particles that are finally selected within R500 of GC1 and GC2 and in the bridge at z = 0.07. Each trajectory shows a colour gradient with time, from dark at z ~ 2 to light at z= 0.07.

In the text

	Fig. 7 Projected mass-weighted temperature map of the galaxy cluster pair at z = 0.07. The map is centred on the system centre of mass and is 7 h−1 cMpc per side, and the integration along the l.o.s. is also performed over 7 h−1 cMpc. The dotted white circles indicate the R500 radii of GC1 and GC2. The position of group B (diamond) is also marked.
In the text

	Fig. 8 Density-temperature phase-space diagram of the gas in the (7 h−1 cMpc)3 region centred on the cluster pair at z = 0. 07, colour-coded by gas Fe abundance (in solar units with respect to Anders & Grevesse 1989). We overplot the contours (cyan lines) enclosing ~25 (inner) and ~75% (outer) of the selected gas mass.
In the text

	Fig. 9 Same as Fig. 8, but only for the gas within R200 for one of the member clusters, namely GC1 (upper panel) and for the gas in the bridge (lower panel) between the pair clusters, at z = 0.07. Dot-dashed lines mark the average temperature and overdensity values (mass-weighted estimate) of the selected gas.
In the text

	Fig. 10 Distributions of Fe abundance and O/Fe abundance ratio for the gas selected at z = 0.07 in different spatial regions and thermal phases. For each gas selection, the corresponding distribution is normalised to the total number of enriched (i.e. ZFe > 0 or ZOFe > 0) gas particles.
In the text

Fig. 11 Redshift evolution of the iron abundance (left) and O/Fe abundance ratio (right) distributions for the ICM in one of the clusters (GC1; top) and the bridge gas (bottom), selected at z = 0.07 and tracked back in time. Vertical lines mark the MW average abundance at any corresponding redshift. The bridge histogram line for z = 4.23 is thinner to mark the very low statistics of metal-enriched gas particles (only ≲ 100, i.e. ~0.6% of the total).

In the text

Fig. 12 Evolution with redshift of enriched and star-forming mass fractions for the gas in different spatial regions and thermal phases, selected at z = 0. 07 and tracked back in time. Upper inset: mass fraction of Ferich gas (ZFe > 0; diamonds) and highly enriched gas ($(Z_{\text{Fe}}>Z_{\text{Fe,mw}}^{z=0.07};\text{triangles})$ triangles). Lower inset: redshift evolution of the star-forming (SFR > 0) gas mass fraction.

In the text

	Fig. 13 Redshift evolution of the total mass M500 for all the clusters and groups identified at z = 0.07 in the (20 h−1 cMpc)3 volume around the cluster pair. At each redshift up to z ~ 1, we consider the main progenitor of the given halo identified at z = 0.07.
In the text

Fig. 14 Infall motion and mass evolution of groups A (left panel) and B (right panel) towards the cluster pair. Upper panels: redshift evolution of the alignment between radial direction towards the pair (centre of mass, blue; GC1, red; GC2; yellow) and subhalo velocity (as in the legend), quantified by the cosine of the angle between them. For Group B (right), the gas < R500 at z = 0.16 is tracked to lower redshifts after entering the GC1 virial boundary. Lower panels: redshift evolution of the subhalo DM (solid lines) and gas (dashed lines) bound mass.

In the text

	Fig. 15 Radial bulk (solid lines) and bound-gas (dashed lines) velocities (upper panels) and their relative absolute difference (lower panels) for groups A and B (left and right, respectively) as a function of the radial distance from the system (centre of mass in blue, GC1 in red, and GC2 in yellow).
In the text

	Fig. B.1 Redshift evolution of the centre shift between the (bound) gas centre of mass and the halo centre in units of R500 for groups A (left) and B (right) analysed in Sec. 4.4. We report the modulus and each separate component, as in the legend.
In the text

	Fig. C.1 Iron abundance and O/Fe ratio distribution for the hot gas and WHIM (upper and lower rows, respectively) for different redshifts as in the legend. The left panels show the mass fraction distribution of the hot-gas iron abundance for comparison. At each redshift, the mass fractions are normalised to the total gas mass of the selected subsample.
In the text