The SRG/eROSITA All-Sky Survey View of the Virgo Cluster

Context. As the closest galaxy cluster, the Virgo Cluster is an exemplary environment for the study of large-scale filamentary structure and physical e ff ects that are present in cluster outskirts but absent from the more easily studied inner regions. Aims. Here, we present an exploration of the SRG / eROSITA data from five all-sky surveys. Methods. eROSITA allows us to resolve the entire Virgo cluster and its outskirts on scales between 1 kpc and 3 Mpc, covering a total area on the sky of about 25 ◦ by 25 ◦ . We utilize image manipulation techniques and surface brightness profiles to search for extended emission, surface brightness edges


Introduction
The Virgo Cluster is the closest (∼ 16.1 Mpc, Tonry et al. 2001) galaxy cluster and one of the best studied.In the X-ray band, in which the strongest emission processes are those of thermal bremsstrahlung and line emission in the hot intracluster medium (ICM), the interior of Virgo and its brightest cluster galaxy (BCG) M87 have been the subject of much previous study (Böhringer et al. 1994(Böhringer et al. , 1995;;Schindler et al. 1998;Churazov et al. 2001;Young et al. 2002;Forman et al. 2005;Werner et al. 2006;Forman et al. 2007;Simionescu et al. 2008;Million et al. 2011;Arévalo et al. 2016;Gatuzz et al. 2022Gatuzz et al. , 2023)).Other works have characterized portions of the outskirts of Virgo (Urban et al. 2011;Simionescu et al. 2017;Mirakhor & Walker 2021).The outskirts of clusters, which Reiprich et al. (2013) defines as between r 500 and 3r 200 , can be a location of physical effects that are no longer observable in the more relaxed central regions.These include breakdown of hydrostatic, thermal, equipartition, and ionization equilibrium states, as well as structure formation effects, such as clumpy gas distribution or accre-tion shocks (see Reiprich et al. 2013, Walker et al. 2019 for a review of these processes in cluster outskirts).Limited by high particle backgrounds, however, previous X-ray instruments have been largely unable to observe significant emission in the outskirts, or else have had too narrow a field of view (FOV) to cover the entirety of the closest clusters.
The extended Roentgen Survey and Imaging Telescope Array (eROSITA) was launched aboard the Spektrum Roentgen Gamma (SRG) mission (Sunyaev et al. 2021) in 2019 with a main science goal of detecting the hot ICM of galaxy clusters and groups to enhance the study of cosmic structure evolution (Predehl et al. 2021).eROSITA is expected to detect 10 5 galaxy clusters and groups through its all-sky survey (Merloni et al. 2012;Liu et al. 2022;Bulbul et al. 2022).eROSITA's combination of soft-band sensitivity and a wide FOV makes it the ideal instrument for the study of faint emission in cluster outskirts that may have previously gone undetected.Other current generation X-ray instruments, such as Chandra and XMM-Newton, have a relatively small FOV and are thus inefficient in mapping large volumes of the Universe; only through mosaics have they been able to study out to and beyond r 200 .The only imaging Xray all-sky survey to precede eROSITA is ROSAT, which was performed over six months in 1990.By the completion of the eight all-sky surveys planned for eROSITA, its survey results will be 25 times more sensitive than ROSAT in the soft band (0.2-2.3 keV, Merloni et al. 2012;Predehl et al. 2021).
The Virgo Cluster's distance, at which 1 arcminute corresponds to 4.65 kpc, provides the opportunity to study the ICM on smaller scales than is possible with any other cluster.Recent studies with other X-ray instruments have taken advantage of Virgo's distance, modest physical size, and dynamism to better analyze its outskirts.These include a study by Urban et al. (2011), which uses a series of thirteen XMM-Newton pointings out to the virial radius in the north to characterize the outskirts with spectral analysis; observations with XMM-Newton to study gas clumping by Mirakhor & Walker (2021); and a Suzaku mosaic along four arms out to the virial radius, which Simionescu et al. (2017) uses to analyze the thermodynamics of the outskirts.All of these studies were limited to analysis along a single radial line to r 200 due to the narrow FOV of the instruments used.
Here, we present an eROSITA imaging analysis of the Virgo Cluster produced from the cumulative data of its first four complete surveys and partial fifth survey.This is the first comprehensive X-ray imaging analysis of the cluster since ROSAT and it reveals significantly more detail.This article is organized as follows: in Sect.2, we describe the data reduction steps and analysis strategy.In Sect.3, we present what the imaging analysis and surface brightness analysis reveal about the large-scale structure of the cluster.We describe the results of deprojection analysis in Sect. 4. Section 5 describes a region to the southwest of M49 where eROSITA soft-band X-ray emission is found to be coincident with LOFAR radio emission.We conclude with a summary in Sect.6.
Throughout this work, we assume a ΛCDM cosmology with Ω m = 0.27, Ω Λ = 0.73, and H 0 = 70 km/s/Mpc.The distance to the Virgo Cluster is taken to be that of its central galaxy, M87, at 16.1 Mpc (Tonry et al. 2001), or z ≈ 0.00377.For consistency with other works (Urban et al. 2011;Simionescu et al. 2017), we consider r 200 , the radius within which the mean density is 200 times the critical density of the Universe, to be the virial radius.The virial radius of the Virgo Cluster is taken to be 1.08 Mpc (Urban et al. 2011), or a projected radius of ∼ 3.9 • .Other radii (i.e., r 500 ) are calculated with the relationships given in Reiprich et al. (2013).Abundances are given relative to the solar abundance of Asplund et al. (2009).Unless otherwise stated, all uncertainties are given at a 68 percent confidence level.

Data reduction and analysis
Observations of the Virgo Cluster from the first five eROSITA all-sky surveys (eRASS:5) were used in this analysis, although the fifth survey only covers approximately half of the area (see Merloni et al. 2024 for all-sky survey information, sky tile definition, etc.).A comprehensive list of the 62 utilized sky tiles, which span an approximate area of 25 • by 25 • , similar to the area covered by >2500 XMM-Newton EPIC pointings, can be found in Appendix A. The processing configuration c020, extended Science Analysis Software (eSASS, Brunner et al. (2022)) version eSASSusers_211214, and Heasoft version 6.29 were used throughout the data reduction.

Data preparation
Upon retrieving the data, the first action was to merge the eRASS:5 event lists of all seven1 telescope modules (TMs) across the full energy range 0.2-10 keV for each tile using the eSASS task evtool.In this work, we use the notation that TM 0 is the combination of all seven TMs.The parameters pattern=15 and flag=0xe00fff30 were applied to select single, double, triple, and quadruple patterns and remove bad pixels and strongly vignetted corners of CCDs.The eSASS task flaregti was then run with parameters source_size = 150, which is in units of arcseconds and dictates the diameter of the extraction area for dynamic threshold calculation, pimin=5000, which is in units of electron volts and controls the lower energy range used, and gridsize=26, which determines the number of grid points per dimension for the dynamic threshold calculation.This task produces updated good time intervals (GTIs) that exclude times during which there is flare contamination.In addition to removing several short flares, flaregti aided in the removal of a flare that spanned the entire area of analyzed sky during eRASS 2.
Following flare filtering, the eSASS task radec2xy was applied to recenter all tiles around a common center, R.A. 187.7042 • , Dec. +12.3911• .Updated GTIs were applied using evtool with the parameters gti="FLAREGTI" and rebin=124, which results in a pixel size of 6.2 arcseconds.This choice allowed the entire sky area to fit into the maximum eSASS image size.The combined event file was then split into its seven individual TMs for further image processing, with the ultimate goal of a cleaned and corrected TM 0 image.The energy band for the analysis was restricted to 0.3-2.0keV for TMs with on-chip filter (TMs 1, 2, 3, 4, and 6; combined, they are referred to as TM 8) and 0.8-2.0keV for TMs without on-chip filter (TMs 5 and 7; referred to as TM 9).The reason for this differing lower limit is the optical light leak contamination of TM 9, which is further detailed in Predehl et al. (2021).Both vignetted and nonvignetted exposure maps were created for every image using the task expmap.

Image corrections
In order to better isolate the ICM emission visible in the eROSITA image, it is necessary to account for the particle background, spatially variable absorption, and variable coverage of the field by the different TMs.

Particle-induced background
The particle-induced background (PIB) subtraction steps are based on knowledge of the eROSITA Filter-Wheel-Closed (FWC) data (Freyberg et al. 2020;Yeung et al. 2023) and were performed as described in Reiprich et al. (2021).Important TMspecific values for this calculation can be found in Table A.1.This process results in PIB maps for each TM, which are then added to create a TM 0 PIB map that can be subtracted from the photon image.

Relative absorption correction
The next step was to account for the spatially variable absorption, since N HI varies across the large sky area and can poten-  2016).Its values span (0.71−5.95)×10 20 cm −2 .Thus, in order to calculate the N HI correction factor across the entire image, spectra for 525 N HI values in this range were simulated to quantify the effect of varying absorption on the X-ray spectra and overall count rates.The simulations were conducted in xspec version 12.12.0(Arnaud 1996).The fakeit command was used with the model where apec 1 accounts for the X-ray emission from the unabsorbed Local Hot Bubble (LHB), tbabs corresponds to the absorption along the line-of-sight, apec 2 accounts for the X-ray emission of the Milky Way Halo (MWH), and pow is the power law component originating from unresolved active galactic nuclei (AGN).The temperatures, metal abundances, and redshifts of the LHB and MWH are fixed, along with the power law photon index.Due to the origin within the galaxy, both metal abundances are set to 1 Z ⊙ , while both redshifts are set to zero.All values can be found in Table A.2.The simulations were run for TMs with and without on-chip filter.This required the use of their separate response files and soft energy bands.The result of the simulated spectra was count rates for every N HI value.From these, a ratio of the count rate for a given N HI to the count rate of the reference N HI can be calculated.The reference N HI is the median N HI across the total sky area and in this instance was 3.33 × 10 20 cm −2 .The correction factor map is then acquired by distributing the ratio corresponding to a given N HI value to its area on the sky.The correction maps for TM 8 and TM 9 are applied by multiplying them by the exposure maps.
The variation of the correction factors is 18.6% for TM 8 and 7% for TM 9. Since absorption is weaker at higher energies, the smaller variation of TM 9 is expected.

Exposure correction
Since TM 8 and TM 9 feature different soft responses and different energy bands to account for the light leak, it is necessary to apply a correction factor when combining all TMs and observations.This process, referred to as exposure correction and described in larger detail in Reiprich et al. (2021), involves using the TM 9-to-TM 8 ratio of the N HI -corrected, PIB-subtracted count rates.The correction factor for this analysis is 0.392.Exposure maps prior to and following absorption and exposure correction can be seen in Fig The final image is found by adding the individual TM event files together to create the TM 0 photon image, subtracting the PIB map, and dividing by the corrected exposure map.It is shown in Fig. 1.Several prominent groups, galaxies, and clusters are labeled.

X-ray images
An RGB image can be seen in Fig. 2.This was created following the steps described in Sect.2, except that different energy ranges were used.The colors and corresponding bands are: red for 0.3 − 0.6 keV, green for 0.6 − 1.0 keV, and blue for 1.0 − 2.3 keV.We note that TM 9 was not used to create the red band, and was only used in the range 0.8 − 1.0 keV for the green band.In the image, the diffuse ICM emission around M87 is visible in a hazy white, while many of the galaxies in Virgo, most obviously M86 and NGC 4636, appear green-blue.This implies that these regions have a harder spectrum than the surrounding medium, which could result from a higher temperature.Outside of the virial radius, extended white-green emission is also visible surrounding and to the southwest of M49.The point sources are likely a combination of foreground stars and background AGN.The variety of point source colors implies a plethora of spectral shapes.The northern eROSITA bubble (Predehl et al. 2020), also known as the North Polar Spur, is bluer towards the Galactic center (outside the analyzed area to the southeast), and then fades to green and then red closer to the center of Virgo, indicating that the highest energies are closest to the Galactic center.This observation is in agreement with the assertion in Predehl et al. (2020) that the bubbles are likely due to energy injections from the Galactic center.
We can infer from the color of the eROSITA bubbles that the soft Virgo Cluster emission in the southeast is likely contaminated by bubble emission.Comparison of the surface brightness in the lower left and upper right of the image revealed that the bubble emission enhances the foreground by a factor of about two.In the blue band, however, there is almost no contribution from the eROSITA bubbles within ∼ 3×r 200 , and the bubbles are absent inside r 200 .This inspired the surface brightness analysis of the cluster in the 1.0 − 2.3 keV band, described in Sect.3.2.
As the focus of this study is the ICM emission, wavelet filtering was used for point source detection and enhancement of diffuse emission.The steps presented in Pacaud et al. (2006) and Reiprich et al. (2021) were followed to apply the wavelet filtering algorithm.Different parameters were set to manipulate the scales and thresholds depending on whether the goal of the filtering run was to detect point sources or emphasize significant ICM emission.For point source detection, the signal-to-noise threshold was lower than for ICM detection, a smaller background cell size was used, the maximum scale was lowered, and fewer iterations were required.The 0.3 − 2.0 keV wavelet filtered image can be seen in Fig. 3, while the "blue band" 1.0 − 2.3 keV wavelet filtered image is in Fig. B.1.Point source detection was performed with the source detection tool SExtractor (Bertin & Arnouts 1996), with some modifications made after manual inspection.Additionally, background clusters were identified via the MCXC catalog (Piffaretti et al. 2011).Our approach to point source detection also identifies small extended sources if they are above the signal-to-noise threshold of four; therefore, we have also excised additional background clusters that are not included in the MCXC catalog with this approach.Point sources and background clusters were then masked, with the radius iden-tified via SExtractor and out to r 500 , respectively, for all subsequent imaging manipulation and surface brightness analysis steps.
Adaptive smoothing produces an image where the signalto-noise ratio (S/N) at each pixel is approximately equivalent, meaning that fainter areas become more smoothed than brighter areas in an image.In this work, the task asmooth from the XMM-Newton Science Analysis System (xmmsas, SOC 2022) with parameter smoothstyle='adaptive' was used to apply adaptive smoothing to the PIB-subtracted, absorption-and exposure-corrected image with point sources masked.An S/N of ten was chosen, along with minimum and maximum values of zero and 25 respectively for the normalized Gaussian convolvers.Masked point source areas were then effectively filled in with the signal level that surrounds them.The point source excised, adaptively smoothed image is shown in Fig. 4.
Gradient filtering has previously been used to examine features in clusters and simulations (Sanders et al. 2016;Walker et al. 2016).The Gaussian gradient magnitude (GGM) filter was applied following Sanders et al. (2016) using scipy (Virtanen et al. 2020) and σ = 1, 2, 4, 8, 16, 32 detector pixels on an image with cosmetically removed point sources.A sample of the images that best show edges is in Fig. 5.The σ = 4 map shows the well-known M87 arms, while the σ = 8 map reveals a spiral edge feature to the northwest.On the largest scales, an extension to the southwest is prominent.
Looking at the wavelet filtered and adaptively smoothed images in Figs. 3 and 4, faint X-ray emission is visible out to ∼ 4 − 5 • from M87, beyond the virial radius, with the most conspicuous emission to the south out to M49 and beyond.An exception to this is the southwestern side, where there is a noticeable emission drop-off prior to the virial radius, a feature which is especially clear in the wavelet filtered image.To the south within the virial radius, the ICM curves in an extension to the southwest in a spiral shape.Further south, the emission halo around M49 appears more extended to its southern edge than to its northern, and is joined by an apparent extension, or tail, which curves off of the halo to the southwest in the direction of NGC 4261.To the north of the cluster, the emission appears more irregular than in the other directions.The emission seems to clump into a number of extensions stretching away from the cluster, with noticeable examples to the northeast and in the radial direction beyond M86.This finger-like tail protruding from M86, which appears out to a distance of 80 arcminutes or 370 kpc from the center of the galaxy, has been observed previously with other instruments.Using Chandra data, Randall et al. (2008) asserted that the tail is related to ram pressure stripping as M86 falls into the Virgo Cluster and interacts with the ICM.To the southeast, the eROSITA bubble emission dominates all images, contaminating the ICM emission of Virgo.Overall, the cluster structure does not appear spherical or symmetrical.

X-ray surface brightness analysis
To acquire the surface brightness profiles, PyProffit (Eckert et al. 2020) was used.For the full azimuthal profile, and subsequent profiles divided into eighths, 656 concentric annuli centered at the X-ray surface brightness peak, (187.705• , +12.391 • ), and stretching to an outermost radius of 700 arcminutes or ∼ 3r 200 were utilized.Their radial bin distribution is logarithmic and has a minimum bin size of 26 arcseconds, approximately the same as the survey average point spread function (PSF) of eROSITA (Predehl et al. 2021).Fig. 2: RGB map of the Virgo Cluster out to three times the virial radius of the cluster (red for 0.3 − 0.6 keV, green for 0.6 − 1.0 keV, and blue for 1.0 − 2.3 keV) is shown.The surface brightness for all channels has a dynamic range from 10 −5 to 1.8 × 10 −3 counts/sec/arcmin 2 .The map was rescaled by a factor of 5 in each dimension, but no smoothing was applied.
To estimate the magnitude of the cosmic X-ray background (CXB), two box regions of dimension 1.8 • by 3.1 • were extracted from outside 3r 200 to the northwest to avoid cluster emission and contribution from the northern eROSITA bubble.The average value for 0.3-2.0keV was found to be 2.8 × 10 −3 cts/s/arcmin 2 , and for 1.0−2.3keV, 1.4×10 −3 cts/s/arcmin 2 .
A single β-model (Cavaliere & Fusco-Femiano 1976;Voit 2005): was fit to the 0.3 − 2.0 keV azimuthal profile out to 300 arcminutes, where R is the projected distance from cluster center, S X (0) is the normalization factor, r c is the core radius, β determines the slope of the surface brightness profile, and K is the background value.The best-fit r c and β can be found in Table 1.This simple model was chosen to allow for comparison with previous ROSAT publications, which also made use of the single β-model and comprehensive Virgo data.The best-fit slope of 0.43 ± 0.02 is in good agreement with ROSAT (Nulsen & Böhringer 1995), which cites β = 0.46, but the core radius of 47 ± 5 arcseconds is about half the ROSAT value of 115 arcseconds.This difference seems mainly due to the poor fit of the single β-model to the innermost ∼ 2 arcminutes of the cluster combined with eROSITA's much better PSF as compared to ROSAT; additional factors could be different treatment of point sources, background clusters, and the central AGN between this work and Nulsen & Böhringer (1995).
Although the spherical β-model allows us to directly compare to ROSAT results, we can also achieve a better fit to cluster emission by allowing for ellipticity in the model.This has the same form except that R becomes where ϵ is the eccentricity and the coordinates are defined as x e = (x − x 0 ) cos(θ) + (y − y 0 ) sin(θ) and y e = (y − y 0 ) cos(θ) − (x − x 0 ) sin(θ), where θ is the ellipse's angle of rotation, and x 0 and y 0 are the central coordinates.The best-fit elliptical model out to r 200 (also in Table 1) was used to create a flattened image of the cluster center, meaning the surface brightness image was divided by the 2D model image.These data products can be seen in Fig. 6.In the flattened image, M87's characteristic arms are visible, as are several surface brightness edges.One is a weak shock identified in Forman et al. (2007) with Chandra.Two have been identified as candidate cold fronts in previous works and are explored further in the following section.

Eight sectors
To better study features in any given direction, eight 0.3−2.0keV surface brightness profiles were constructed.The choice of eight sectors was motivated by the decision not to bisect M86 or M49, the two largest subgroups, and to achieve a balance between information in more directions vs. to larger radii.Unlike, for example, ten sectors, eight sectors allows significant measurements to the outskirts of Virgo (beyond r 500 ) while offering more directional information than, for example, six sectors.The regions are shown in Fig. 7 and named for the eight half-wind compass points (north-northwest, east-northeast, etc.).The resulting surface brightness profile for each sector is seen in Fig. 8.In this figure, each plot has a single profile for every sector, which is plotted in the color it corresponds to in Fig. 7.The other seven sectors for any given profile are plotted in gray and the azimuthal average in black to offer a better idea of how the profiles compare.Profiles for the eight sectors were also created in the blue band of the RGB image, 1.0 − 2.3 keV, which is approximately free from eROSITA bubble contribution.These can be seen in Fig. B.7, and are referred to as the "blue" profiles below.
The WNW sector contains a bump in the profile at 10 arcminutes corresponding to a surface brightness excess, and is followed by a drop-off at 20 arcminutes.A prominent peak in the profile is visible at ∼ 70 arcminutes or 330 kpc, corresponding to the X-ray halo around M86.The profile then gradually flattens out into the outskirts past r 200 , where it shows a low level of surface brightness as compared to the other sectors.In the blue band, M86 is still prominent.
In the NNW sector, which does not contain any of the brightest Virgo galaxies or groups, the profile displays a number of interesting features.Between 10 and 20 arcminutes, there is an excess of brightness followed by a discontinuity in the profile.This behavior is often attributed to shocks or cold fronts, but to discern between the two, temperature and density profiles of the region are necessary (Markevitch & Vikhlinin 2007).This particular feature was previously studied with XMM-Newton , Chandra, and Suzaku data (Urban et al. 2011;Werner et al. 2016;Simionescu et al. 2017) and identified as a cold front.Another excess as compared to the other sectors is visible between 50 and 100 arcminutes.This could be associated with the apparent extension of ICM emission in that direction.Past r 200 , this sector shows among the lowest overall brightness when compared to the other sectors.This is expected because NNW is the furthest from the eROSITA bubbles.In the blue band, all of the same features are visible.
In the NNE, a similar discontinuity as in NNW is visible at the same distance, ∼ 20 arcminutes, indicating that the cold front extends along the entire northern direction.There is a bump in the profile just past r 500 that could be due to the finger-like extension of the ICM observed in that direction.The outskirts past r 200 are at a similarly low level as in NNW with a similar shape in the full and blue profiles, indicating that there may not be substantial contribution by the eROSITA bubble in this region.
To the ENE, contribution from the eROSITA bubble becomes apparent both in the images and in the surface brightness profile, where the brightness increases from r 200 to 3r 200 .The blue profile instead flattens out in this range.
The ESE sector shows a slow increase from r 200 to 3r 200 , a segment in which it registers the highest surface brightness of all sectors due to the eROSITA bubble.The blue profile for this region is not completely flat, but instead also begins to increase around 400 arcminutes, implying some contribution from the eROSITA bubble in the 1.0 − 2.3 keV band to the ESE.However, the blue surface brightness in this range is still a factor three below the 0.3 − 2.0 keV surface brightness when normalized to the same central value, implying that most of the contribution from the eROSITA bubble is in the 0.3 − 1.0 keV band.
The SSE sector is brighter on average than other sectors around 40 − 50 arcminutes, and then drops off gradually.The eROSITA bubble clearly also affects the outer regions of this sector, which increases between r 200 and 3r 200 and shows a prominent bump just after r 200 which becomes even more prominent in the blue band.This bump is prominent even following the removal of extended sources, or clumps, before production of the surface brightness profiles.
The SSW sector contains M49 and its extended ICM tail, features that appear in the profile in the form of a bump between 200 and 300 arcminutes.A discontinuity, again possibly indicative of either a shock or cold front (Markevitch & Vikhlinin 2007), is present in this profile at the same location as in the SSE around 40 − 50 arcminutes, except that it appears more pronounced in this sector.
Finally, to the WSW, the profile records the lowest values of all segments between 10−20 arcminutes, indicating that the cluster emission is not stretched along this direction, which matches our imaging analysis.There is a minor bump and drop around 50 arcminutes, perhaps indicative of the same feature present in the SSW and SSE.
Overall, if we compare the blue and full surface brightness profiles normalized to the same central values, all of the same features observed in the full profile are present in the blue band.The scaled difference to the median is the difference between each sector and the median divided by the uncertainty ("in sigmas").This confirms that these features are not due to the eROSITA bubble contribution.

Cold fronts
Two surface brightness discontinuities are apparent in the eight sector analysis, including one to the north that was previously identified in full by other X-ray instruments, and another to the south that is visible in full for the first time in this work.Both are marked in Figure 4. We extracted surface brightness profiles across these discontinuities (Fig. B.2 and Fig. B.3) in order to better investigate their properties.
At the location of the discontinuity found with eROSITA in the NNE and NNW sectors, XMM-Newton , Suzaku, and Chandra all record the same feature at this distance from cluster center to the north.These other instruments did not all cover the same portion of this feature; Simionescu et al. (2010) 2016) uses a northwestern Chandra pointing to chart the sharpest part of the front.The location of the discontinuity, at 19.1 arcminutes or 90 kpc from cluster center, is identical across instruments.Although spectral analysis is outside the scope of this work, Simionescu et al. (2010) uses spectral analysis to derive the temperature, density, and pressure profiles for this region, finding that the discontinuity is evidence of a cold front rather than a shock.The signature of cold fronts is a discontinuous drop in temperature paired with a jump in density, while shocks show an increase in temperature, density, and pressure (Zinger et al. 2018).The analysis by Werner et al. (2016) finds that it is sharper to the north than the west, which is in agreement with our surface brightness findings.Simionescu et al. (2010) also shows that the cold front is not centered spherically around M87, a conclusion supported by using the GGM and flattened images and experimentation with the eROSITA surface brightness of different regions around the location of the cold front to best determine its location.The cold front location is marked by a series of white arrows in Fig. 4. The front is closer to M87 in the west than the east, matching the spiral morphology of gas sloshing in simulations (Markevitch & Vikhlinin 2007).A single β-model with a discontinuity was fit to the eROSITA data in order to quantify the cold front.This model is found by scaling a standard single β-model by a sigmoid function.This takes the form where a is the relative decrease of the surface brightness at the discontinuity, b is the steepness of the jump, and R 0 is the location of the discontinuity.The values for the best fit can be found in Table 1 and the surface brightness profile across the discontinuity and model are shown in Fig. B.2.We measure a surface brightness decrease at the discontinuity of 41%.
To the south, the feature between 40 − 60 arcminutes (∼220 kpc) from center, with the distance dependent on the sector in which it is found (SSW, SSE, WSW) looks similar to the cold front feature in the north, although less sharp.A surface brightness discontinuity within this range, presumed to be the same feature, is also present in Suzaku data.Simionescu et al. (2017) performs spectral analysis of two discontinuities, one at 233 kpc in the western arm and one at 280 kpc to the south, and Notes.The total is fit to both spherical and elliptical β-models, while the "north CF" and "south CF" profiles are across the cold fronts and fit to a β-model with a discontinuity. (a) The core radius r c is in units of arcseconds.
finds that their temperature profiles are also indicative of a cold front.These two discontinuities are along the single discontinuity that we find, which is not spherical with respect to M87 but is closer to it in the west than the south.Future work should confirm that the discontinuity is consistent with a cold front along its entire extent.We fit Eq. ( 4) to the southern profile and found a less steep and less drastic decrease than the northern cold front.Best-fit values can be seen in Table 1 with the profile shown in Fig. B.3.The off-set cold fronts in both the northern and southern directions could together be evidence of large-scale gas sloshing in the Virgo Cluster, a phenomenon caused by the infall of subclusters to the cluster core that displace the cool gas (Markevitch & Vikhlinin 2007).Similar evidence for this type of gas sloshing in the Virgo Cluster was found in Simionescu et al. (2007); Simionescu et al. (2010); Gatuzz et al. (2022).

Gas clumping
Based on Zhuravleva et al. (2013), which shows that the median surface brightness in simulated galaxy clusters is robust against gas inhomogeneities, Eckert et al. (2015) proposed a method for recovering unbiased density profiles using the azimuthal median of the surface brightness in radial annuli.They apply the method to 31 clusters observed with ROSAT/PSPC.In this work, we applied the method, which involves making use of a 2D binning algorithm based on Voronoi tessellation (Cappellari & Copin 2003), to eROSITA data.The VorBin python package was used (see Sect. 3.3 of Eckert et al. (2015) for further detail on methods).We computed the median surface brightness profile from the binned image and compared it with the surface brightness profile obtained by averaging.The comparison of the two methods can be seen in Fig. B.4.From these profiles, we computed the emissivity bias, a measure used in Eckert et al. (2015) as a projected 2D surface brightness distribution proxy for the true 3D clumping factor of the gas.The emissivity bias obtained using the median and mean surface brightness profiles is plotted as a function of radius in Fig. 9, out to 1σ detection significance.Our results match the trend found in Eckert et al. (2015); the emissivity bias is low (∼ 1.1) in the inner regions of the cluster, and increases with radius in the outskirts (∼ 1.2).Mirakhor & Walker (2021) studies the clumping in Virgo with XMM-Newton in the northern direction, and also finds mild clumping with a slight increasing trend.They report that their point source detection algorithm was likely able to exclude a significant fraction of clumps; this could well be the case in this work.As in that study, we also find that the exclusion of extended sources (primarily M86, M49, and M58) decreases the clumping factor significantly.After the exclusion of the emission around these three Virgo member galaxies, the measurement is consistent with the findings of Mirakhor & Walker (2021).When the emissivity bias is calculated in the 1.0−2.3keV band, the same trend exists, with the bias increasing to ∼ 1.2 past r 200 , implying that the trend is due to true cluster emission rather than foreground contribution from the eROSITA bubble.
Cosmological simulations predict that gas streams along cosmic web filaments can penetrate the inner regions of clusters, with gas clumping levels higher along filamentary directions (Zinger et al. 2016;Vazza et al. 2013).Simionescu et al. (2017), using the strip of Suzaku pointings, finds that gas density is higher along the north-south axis of Virgo than the east-west.They also find a steep increase in gas mass fraction beyond the virial radius, especially along the north-south axis.These findings support the idea that a cosmic web filament, complete with gas clumping, may run north-south in the Virgo Cluster, an idea also reinforced by the distribution of galaxies in the local Universe (Castignani et al. 2022).Mirakhor & Walker (2021), on the other hand, finds uniform clumping factor measurements in all directions when considering all available XMM-Newton pointings.
To investigate whether eROSITA data indicates significant differences in clumping in varying directions, we also calculated the emissivity bias as a function of radius within the eight sectors identified in Fig. 7. Larger radial bin widths and fewer bins than for the azimuthal assessment were used so as to achieve significant measurements out to ∼ r 200 in every direction.For this calculation, the extended halos around M86, M49, and M58 were removed.
The emissivity bias of each individual sector, shown in corresponding color in Fig. B.5 with the median in black, follows the same trend as the total profile.The values fluctuate around 1.0-1.1 until ∼ 60 arcminutes, after which point they increase to a median value of 1.2 just beyond r 200 .Half of the sectors have single bins that spike to 1.3-1.4;namely, ENE, WSW, ESE, and WNW, all of which span the east-west axis of the cluster.However, no one sector is systematically higher than another.When the median of the four north-south axis sectors and four eastwest axis sectors is compared (see Fig. B.6), the measurements are within uncertainties of one another.That is, there is no apparent excess of clumping in the north-south axis as compared to the east-west out to r 200 .Fig. 10 shows the good agreement among b X measurements in each sector beyond 0.75r 500 .
Overall, the effect of clumping in the Virgo Cluster appears to be mild and broadly consistent in each radial direction.Moreover, with the full and, at the same time, detailed eROSITA view, we are not only able to measure the bias out to ∼r 200 , but are also able to trace back the causes for stronger bias to individual sources.That is, the largest values of the emissivity bias occur at the locations of prominent galaxies and groups within Virgo, such as M86 at ∼ 100 ′ , M58 just beyond it, and M49 at ∼ 300 ′ .

Deprojected gas density and mass
Using PyProffit, deprojection was performed in order to extract gas density and gas mass estimates.We follow Eckert et al. (2020), which gives details of the application of the deprojection to eROSITA simulations.To calculate the conversion factor between count rate and emissivity, we assume an abundance   2017), and the median N H across the analyzed area from Sect. 2. The conversion factor from count rate to emission measure is highly insensitive to temperature and abundance in this case.The multiscale decomposition deprojection was then performed.This allows for the extraction of the gas density profile, and subsequent calculation of the gas mass profile, which is the deprojected gas density integrated over a given volume.
We applied these steps to the mean surface brightness profile.The gas density profile is plotted in Fig. 11, and the gas mass in Fig. 12.
The hydrogen number density at 0.5 arcminutes is estimated to be 0.06 ± 0.01 cm −3 and the density at r 200 = 1.08 Mpc is (4.37 ± 0.51) × 10 −5 cm −3 .The density profile approximately follows a power-law model n H ∝ r −k where k = 1.25 ± 0.05 between 1.5 and 100 arcminutes.The slope is flatter inside 1.5 arcminutes.The overall shape of the density profile matches typical profiles measured by ROSAT in that most clusters have a steepening profile beyond ∼ r 500 (Eckert, D. et al. 2012).
The deprojected mean profile yields gas masses M gas,r<r 500 = (1.05± 0.32) × 10 13 M ⊙ and M gas,r<r 200 = (1.98 ± 0.60) × 10 13 M ⊙ .The gas mass profile has a best-fit power-law model with slope k = −1.75 ± 0.04.Scaling relations from Arnaud et al. (2005) predict a total gravitating mass of M 200 = 1.4 × 10 14 M ⊙ for a cluster temperature of 2.3 keV, while Simionescu et al. (2017) calculates an estimated virial mass of (1.05 ± 0.02) × 10 14 M ⊙ .This places our M gas,r<r 200 gas mass estimate at approximately 14 − 19% of the predicted total mass.Typical estimates of gas fraction in clusters are on the same order as the cosmic baryon fraction, ∼ 15% (Aghanim et al. 2020), though Simionescu et al. (2017) finds increasing gas fraction with radius in Virgo, greater than 20% for some radial directions beyond r 500 .When the deprojection is performed only for portions of the cluster that are not affected by the eROSITA bubble (i.e., northern and western segments), estimated masses are ∼ 10% lower, still within reported uncertainties of the mean profile.
It should be noted that the deprojection process requires the assumption of spherical symmetry, which we have shown in the rest of this work is a poor assumption in the case of Virgo due to its asymmetrical structure outside of the inner region.It is promising, however, that the gas clumping in all directions is shown to be mild and consistent.The non-uniformity of clumpy gas leads to an overestimate of density and therefore of mass; if clumping is mild and uniform, as our analysis shows, the resulting density and mass is less biased.Eckert et al. (2015) finds typical b X values for their cluster sample that broadly agreed with those found here out to r 200 , although with larger values and uncertainties at higher radii.Their clumping factors resulted in a gas mass overestimate in that work of ∼ 10% at r 200 .Therefore, the impact in this case, with lower b X , would be ≤ 10% gas mass bias above the true value due to gas clumping.

M49 region
A bright extended region, that is, to our knowledge, first studied in X-ray here with eROSITA, is located to the southwest of the M49 group, bent at a 50 • angle below the R.A. axis.It spans a projected length of 320 kpc (1.1 • ), is 140 kpc (0.5 • ) wide, and is approximately 1.3 Mpc from the center of the Virgo Cluster, outside r 200 .A surface brightness profile across the extension shows the morphology as singly peaked (in contrast to doublepeaked stripping tails, see e.g., Ge et al. 2023).To determine the significance of this enhancement, the surface brightness and propagated error was extracted from a number of 20 ′ ×20 ′ square regions along the extension (so-called source boxes) and just outside of it (so-called background boxes).A mean value for the source boxes and for the background boxes was calculated.The error was found with where σ src is the standard error of the mean surface brightness of the source boxes and σ bkg is the standard error of the mean surface brightness of the background boxes.The relative difference between the extension and background is found to be 41.75%, or a significance of 3.15σ.We searched the projected area around this feature for clues as to its origin.It is coincident with the so-called W'-cloud, made up of elliptical galaxy NGC 4365 and several member galaxies, which has been previously studied in the optical and infrared (de Vaucouleurs 1961;Binggeli et al. 1987;Mei et al. 2007).This group is located at ∼ 23 Mpc (Mei et al. 2007) and is proposed to be falling into Virgo from behind.NGC 4365 and fellow group member NGC 4342 are coincident and thus likely responsible for the peaks in X-ray emission within the extension.The extended X-ray emission could be the W' intra-group medium, extended due to its fall into the cluster.Nearby, to the northeast of the extension, the M49 group is likely falling into the cluster from the south (Su et al. 2019).Su et al. ( 2019) finds an M49 stripped tail of smaller dimensions (70 kpc long by 10 kpc wide), which they argue is evidence of the brightest group galaxy (BGG) moving relative to the M49 group gas.This smaller tail associated with the BGG is also visible with eROSITA; the longer exten-Fig.13: Image displaying the region to the southwest of M49 in radio with LOFAR (left) and in X-ray with eROSITA (right).Contours of the radio image have been overlaid on the X-ray image for easier comparison, and images are displayed with the same frame.The white line in the upper left is 100 kpc in length.sion matches the morphology and direction of the smaller one.The longer extension could be a stripped tail of intra-group gas, as opposed to the smaller stripped tail of intergalactic gas.
A diffuse feature with a similar extent, location, morphology, and orientation as the X-ray feature was also recently identified in the radio band with LOw-Frequency ARray High-Band Antenna (LOFAR HBA) as part of ViCTORIA (Virgo Cluster multiTelescope Observations in Radio of Interacting galaxies and AGN) project (Edler et al. 2023).The radio extension has a length of ∼ 1 • , about 280 kpc, and shares an orientation with the radio tails around M49, although no clear connection can be seen between structures.Figure 13 shows a cut-out of the region in radio, and in X-ray with radio contours overlaid, showing the similar structures of X-ray and radio emission.
Based on its morphology and shared location, we suggest that the X-ray extension could be causally connected to the radio.The feature could be related to the in-falling of diffuse gas associated with the W' cloud group or stripped gas of the M49 group and its interaction with the ICM.The radio feature could possibly originate in an accretion shock or due to turbulence in the hot plasma (van Weeren et al. 2019).An alternative explanation for the radio feature, which does not, however, explain the similar morphology in the radio and X-rays, is that the radio emission is due to remnant AGN plasma related to a previous active phase of NGC 4365 (Spasic et al. 2024).A follow-up on its X-ray and radio spectral properties is required to gain a better understanding of this feature.

Conclusion
In this work, data collected during the four complete all-sky surveys and partial fifth survey of eROSITA is used to study the Virgo Cluster in an area spanning approximately 25 • by 25 • .Data reduction and preparation include flare removal and corrections for PIB, variation of N HI across the sky area, and exposure differences between TM 8 and TM 9.
Different image manipulation techniques allow for the exploration of larger scale structure.The RGB image reveals extended white-green emission around prominent galaxies and groups, as well as a gradient within the northern eROSITA bubble, which is mainly projected onto Virgo Cluster emission in the red and green bands (0.3-1.0 keV).The bubble fades from blue-green to green (0.6 − 1.0 keV) to red (0.3 − 0.6 keV) as the distance grows from the Galactic center.Wavelet filtering and adaptive smoothing reveal faint ICM emission in several directions out to 4 − 5 • from M87 with a larger radius to the north and south, and a notable drop-off in emission to the west prior to the virial radius.The northern region of the Virgo Cluster contains a number of ICM extensions stretching away from the cluster, while a spirallike extension can be seen to the south.
Surface brightness profiles were created for the full azimuth and eighths of the overall cluster in the 0.3 − 2.0 keV and 1.0 − 2.3 keV bands.The latter band choice arose from the apparent lack of contribution to the Virgo Cluster emission by the northern eROSITA bubble in the blue emission of the RGB image.To the north, in the NNW and NNE subsectors, a prominent discontinuity at ∼ 19 arcminutes marks a 41% decrease in surface brightness.This has been identified as a cold front in other studies that used spectral analysis.To the south, most prominent in the SSW subsector and less sharp in SSE and WSW, another discontinuity in the profile is present at ∼ 48 arcminutes, though with just a 24% decrease.This was also previously identified as a cold front at locations along it to the south and west; however, this is the first time the full extent of the front can be seen.The off-set cold fronts to the north and south of the cluster support a picture of large-scale gas sloshing.The 0.3 − 2.0 keV surface brightness profile to the east is dominated by the eROSITA bubble past r 200 , while the 1.0 − 2.3 keV profile shows only a minor contribution from the bubble in this region.The surface brightness profiles to the west match the visual evidence of the images, displaying few notable features except for an excess at the location of the ICM surrounding M86.All reported features (with the exception of the eROSITA bubble) are present and similarly significant in both the 0.3 − 2.0 keV and 1.0 − 2.3 keV bands.
We make use of Voronoi tessellation to bin the data before extracting a median surface brightness profile.One would expect the median surface brightness profile to be robust against gas clumping down to the scale of the typical Voronoi cell size, unlike the mean surface brightness profile.The comparison of the two can be quantified by the emissivity bias, b X , which is larger in the case of gas clumping.We find a low emissivity bias, about 1.1, in the inner regions, growing to about 1.2 in the outskirts of the Virgo Cluster.Peaks in the emissivity bias can be traced back to M86, M58, and M49, with an increasing trend still apparent following their removal.Emissivity bias was also calculated for each of eight sectors introduced for surface brightness analysis.Although half of the sectors show single bin peaks to values of 1.3-1.4,no single sector has consistently higher gas clumping.The emissivity biases in the outskirts in each direction, and consequently along the north-south and east-west axes, are consistent.The effect of gas clumping from this measure appears to be mild and radially uniform.This detailed measurement is the first of its kind for the full azimuth of Virgo, and offers a high resolution view of clumping that would be difficult to achieve in more distant clusters.
A deprojection analysis was performed in order to determine gas density and gas mass profiles.The result was a virial gas mass of M gas,r<r 200 = (1.98 ± 0.70) × 10 13 M ⊙ , which is on the order of 14-19% of total mass estimates for Virgo, largely consistent with the expected gas mass fraction.This reported gas mass does not incorporate potential overestimates due to contribution from the eROSITA bubbles and due to gas clumping, which could each contribute a ≲ 10% inflation.
An X-ray extension south of M49, outside the virial radius of Virgo, has a length of ∼ 320 kpc.The surface brightness of the extension has a relative difference to the background of 41.75% and a lower-bound significance of 3.15σ.These values were determined using extracted surface brightness values along and just outside of the extension.The X-ray emission is coincident with radio emission observed by LOFAR.The aligned morphology to the M49 BGG stripped tail, shared location with the W' cloud galaxies behind Virgo, and coincident radio emission could together be indication of an accretion process such as a shock or turbulence.
Fig. 1: eROSITA TM 0 particle-induced background subtracted, exposure corrected, absorption corrected count rate image in the 0.3 − 2.0 keV energy band.The image is displayed in logarithmic scaling with a Gaussian smoothing of radius 18 pixels, or 112 arcseconds.The white bar in the lower left spans 1 degree, or 279 kpc.Prominent groups and galaxies within Virgo are labeled in white, while non-members of the cluster are labeled in yellow.

Fig. 3 :
Fig. 3: Wavelet filtered eRASS:5 image of the Virgo Cluster in the 0.3 − 2.0 keV energy band.The r 200 is shown in white.Features mentioned in the text are labeled.

Fig. 6 :
Fig. 6: Series of eROSITA 0.3-2.0keV images centered on M87 with half-size 65 arcminutes.Left: Surface brightness image smoothed with a sigma of 2 pixels, or 12 arcseconds.Point sources are present, but were removed for the fitting of the model.Middle: Best-fit β-model.Right: Flattened image, meaning the image divided by the model, smoothed with a sigma of 5 pixels, or 31 arcseconds.

Fig. 8 :
Fig. 8: Eight surface brightness profiles in the 0.3 − 2.0 keV energy band, one for each sector identified in Fig. 7.The other seven sectors for any given profile are plotted in gray for comparison, with the azimuthal average shown in black.The CXB measured to the northwest beyond 3r 200 is shown in navy blue.A selection of characteristic radii are shown as dash-dotted vertical brown lines.The scaled difference to the median is the difference between each sector and the median divided by the uncertainty ("in sigmas").
utilizes an XMM-Newton mosaic of the inner 200 kpc, Urban et al. (2011) studies the cluster in a strip of pointings directly north from M87 with XMM-Newton , Simionescu et al. (2017) uses a strip of Suzaku pointings to the north, and Werner et al. (

Fig. 9 :
Fig. 9: Emissivity bias, b X , as a function of distance from cluster center.The top panel is the result after removing point sources and background clusters.The bottom panel is after additional removal of the extended gas halos around M86, M58, and M49 to show the upward trend still exists after removal of prominent sources.

Fig. 10 :
Fig. 10: Azimuthal variation in the gas clumping in the radial range 0.75r 500 to 1.2r 200 for all eight sectors.The black line represents the average value.The clumping factor measurements are consistent, with no apparent north-south enhancement.

Fig
Fig. B.7: Eight surface brightness profiles in the 1.0 − 2.3 keV (blue) energy band, one for each sector identified in Fig. 7.The 0.3 − 2.0 keV band for any given profile is plotted in black for comparison.The 0.3 − 2.0 keV CXB measured by eROSITA in the northwest of the cluster is shown as a dotted black line on each plot.The blue band CXB estimate is taken from the same regions as the 0.3 − 2.0 keV CXB, and shown as a dashed navy blue line.A selection of characteristic radii are shown as dash-dotted vertical brown lines.

Table 1 :
Best-fit parameters to surface brightness profiles.