The Next Generation Virgo Cluster Survey (NGVS). III. A Catalog of Surface Brightness Fluctuation Distances and the Three-dimensional Distribution of Galaxies in the Virgo Cluster

The surface brightness fluctuation (SBF) method is a robust and efficient way of measuring distances to galaxies containing evolved stellar populations. Although many recent applications of the method have used space-based imaging, SBF remains a powerful technique for ground-based telescopes. Deep, wide-field imaging surveys with subarcsecond seeing enable SBF measurements for numerous nearby galaxies. Using a preliminary calibration, Cantiello et al. presented SBF distances for 89 bright, mainly early-type galaxies observed in the Next Generation Virgo Cluster Survey. Here we present a refined calibration and SBF distances for 278 galaxies extending several magnitudes fainter than in previous work. The derived distances have uncertainties of 5%–12% depending on the properties of the individual galaxies, and our sample is more than 3 times larger than any previous SBF study of this region. Virgo has a famously complex structure with numerous subclusters, clouds, and groups; we associate individual galaxies with the various substructures and map their three-dimensional spatial distribution. Curiously, subcluster A, centered around M87, appears to have two peaks in distance: the main peak at ∼16.5 Mpc, and a smaller one at ∼19.4 Mpc. Subclusters B and C have distances of ∼15.8 Mpc. The W and W ′ groups form a filament-like structure, extending more than 15 Mpc behind the cluster with a commensurate velocity increase of ∼1000 km s−1 along its length. These measurements are a valuable resource for future studies of the relationship between galaxy properties and local environment within a dynamic and evolving region.


INTRODUCTION
The nearby Virgo cluster of galaxies, along with the groups and filaments that surround and feed into it (e.g., Tully 1982) offers a close-up view of the ongoing assembly of a massive cluster and the transformation of galaxies as they fall into the mix.The properties of galaxies are influenced by both the large-scale environment (e.g., the familiar morphology-density and color-density relations; Dressler 1980;Blanton et al. 2005;Bamford et al. 2009;Wang et al. 2018) and the proximity of neighbors on smaller scales (e.g., Park & Hwang 2009;Ellison et al. 2009Ellison et al. , 2010)).Virgo is actively accreting both isolated galaxies and small groups.The accreted galaxies may undergo dramatic changes in morphology and starformation activity as they encounter the cluster environment and their host groups are subsequently disrupted (Paudel et al. 2013;Boselli et al. 2014;Lisker et al. 2018;Benavides et al. 2020).
Many past studies have divided the galaxies in and around the Virgo cluster into distinct substructures.The first detailed investigation was probably that of de Vaucouleurs (1961), who defined the E, S, W, and X clouds, along with the tighter S , Wa, Wb, and W groups. Cloud E referred to the main body of the Virgo early-type galaxies, while S was the more dispersed distribution of spirals, and X referred to Virgo's Southern Extension, which was further split into two components; additional sub-groupings were added by de Vaucouleurs & de Vaucouluers (1973).These studies relied on sky positions, velocities, morphologies, and some very crude distance estimates based on properties such as angular diameters and apparent magnitudes.
A revised nomenclature, describing the same general structures but in finer detail, was introduced by Binggeli et al. (1987) based on the Virgo Cluster Catalog (VCC; Binggeli et al. 1985).These authors used a combination of velocity and morphological characteristics to exclude obvious background galaxies, then grouped the remaining galaxies based on sky position into six substructures: cluster A, centered on the cD M 87; cluster B, centered on M 49, the brightest cluster galaxy (BCG); cluster C, centered on the giant elliptical M 601 ; the M cloud, identified as a distinct velocity component by Ftaclas et al. (1984); the W cloud; and the W group (with the latter two adopted from de Vaucouleurs 1961).The VCC did not include the Southern Extension.
In this paper, we refer to the A, B, and C groupings as subclusters.Binggeli et al. (1987) described subcluster C as being within the elongated halo of A, with a continuous bridge of early-type galaxies linking the two concentrations.These substructures were further refined by Binggeli et al. (1993), using more extensive radial velocity data.From a combination of the luminosity functions, radial velocities, and some 21 cm line-width distance estimates for spirals, these authors argued that the M, W, and W substructures were significantly more distant than Virgo's three main subclusters.
To probe the three-dimensional structure of Virgo in more detail, Gavazzi et al. (1999) used both fundamental plane (FP) and Tully-Fisher (TF) distances for over 130 early-type and spiral galaxies.Surprisingly, they found that subcluster B, associated with M 49, was about 9 Mpc (one magnitude) more distant than subcluster A, associated with M 87 (they did not distinguish between A and C).Since the two structures have nearly the same mean velocity, the interpretation was that B is falling towards A from behind, at about 760 km s −1 .Interestingly, the distance to M 49 itself was much closer to the mean of subcluster A than to subcluster B, while the reverse was true for M 87.However, the authors es-timated their individual distance errors to be 0.35 and 0.45 mag, or 17% and 21%, for the TF and FP distances, respectively, making it difficult to assign individual galaxies to specific substructures.They further concluded that the M and W clouds were at roughly twice the distance of subcluster A. Several other subgroups that they identified (i.e., clouds N, S and E) were located close to the mean distance of A, but with significantly different mean velocities; this could suggest that they are the bound remnants of groups being accreted by the main cluster.This work supported the general view of Virgo as a cluster still in the process of assembly.
A number of studies have explored the relationship between galaxy properties and the local environment within the extended Virgo cluster region.For instance, Boselli & Gavazzi (2006) analyzed a sample of 868 UVselected galaxies in the Virgo region having velocities below 3500 km s −1 to study the role of environment in the cessation of star formation and the transformation from late-to early-type morphology within high-density regions.In addition to UV and optical photometry, the authors made use of extensive mid-infrared and Hi radio data.They then studied the distribution of redsequence, blue star-forming, and transitional "green valley" galaxies as a function of the estimated local galaxy density.In particular, they examined the UV-optical color-magnitude diagrams of subclusters A, B, and C, as well as the M, W, and W groups, the low-velocity cloud (LVC), and the surrounding field.As expected, red sequence galaxies predominate within the A and C regions, but also within the W and W groups. Subcluster B shows more of a mix of red and blue galaxies, while the M group, LVC, and surrounding field are dominated by blue galaxies.More generally, the red fraction correlates with estimated local galaxy density.However, since reliable distance information was unavailable for most of their sample, the authors made the assignments of galaxies to substructures based on sky position and velocity, and they used only projected densities in their analysis, rather than volume densities.Thus, projection effects were a significant limitation in their analysis.
Going further in eliminating contamination to disentangle the various environmental factors in such a complex region would require precise information on the relative distances of the galaxies.Precise and reliable distances make it possible to estimate the local volume density of galaxies, the depth of a galaxy within the cluster well, and the likelihood of its association with any of the substructures that have been previously identified in wide-field surveys.The virial radius of the Virgo cluster (combining the A, B, C subclusters) is approximately 1.6 Mpc (Ferrarese et al. 2012), or about one-tenth of the distance (Tonry et al. 2001;Cantiello et al. 2018a) to the cluster itself.This means that resolving the internal structure of Virgo would require a distance precision of ∼ 5%.More compact systems at similar distances (e.g., Fornax), or larger but more distant clusters, would require correspondingly higher degrees of precision.For example, any attempt to resolve the line-of-sight structure of the Coma cluster at ∼ 100 Mpc would be wholly impractical.
But even for Virgo, measuring distances of sufficient precision has been a major challenge.Cepheids are the premier extragalactic distant indicator, but Cepheid distances are observationally expensive and cannot be used for the early-type galaxies that dominate the densest parts of the cluster and its associated groups.The tip of the red giant branch (TRGB) method is similarly expensive and requires Hubble or JWST image quality; for this reason, relatively few TRGB distances are available in Virgo (see Bird et al. (2010), Lee & Jang (2017), Mihos et al. (2022) and the Appendix of Blakeslee et al. (2021)).Type Ia supernovae are too rare for mapping individual structures in detail, but the planetary nebula luminosity function (PNLF) method applied using wide-area integral field spectroscopy (Spriggs et al. 2021) presents an intriguing new possibility; the globular cluster luminosity function has also been used to measure distances within Virgo (e.g., Whitmore et al. 1995).Thus far, however, the most detailed work on the three-dimensional structure of Virgo has been carried out using the SBF method; for a recent review of this technique, see Cantiello & Blakeslee (2023).
The first attempt to resolve the depth of the Virgo cluster with SBF distances was that of Tonry et al. (1990) soon after the method was introduced (Tonry & Schneider 1988).These authors showed definitively that NGC 4365 was more distant than the other Virgo giant ellipticals and suggested that it was a member of the W cloud (now it is recognized as belonging to W ; Mei et al. 2007;Cantiello et al. 2018a).However, the stellar population dependence of the absolute I-band SBF magnitude M I (Tonry 1991;Tonry et al. 1997) was not recognized at the time, and the other depth effects reported by Tonry et al. (1990) are attributable to bluer galaxies having brighter M I values.
A decade later, West & Blakeslee (2000) used data from the ground-based SBF Survey of Galaxy Distances (Tonry et al. 1997(Tonry et al. , 2001)), fully corrected for stellar population effects, to constrain the shape and orientation of the main body of Virgo which runs along the axis joining subclusters A and C.This study concluded that the fourteen brightest ellipticals in this central region follow a roughly co-linear distribution in three dimensions along an axis connecting to a filament that extends in the direction of Abell 1367.Neilsen & Tsvetanov (2000) observed a similar trend using SBF distances measured from Hubble Wide Field Planetary Camera 2 observations of a dozen bright Virgo galaxies along the A-C subcluster axis; they also confirmed that NGC 4365 (located south of this axis) is significantly more distant.West & Blakeslee (2000) did not consider subcluster B or the W/W groups in their analysis, but they showed that NGC 4168, which is now known to be a member of the M group (Cantiello et al. 2018a), lay along an extension of the same "principal axis".However, the study was limited by the relatively small number of sample galaxies and the ∼ 10% errors for these ground-based distance measurements.Mei et al. (2007) examined the detailed structure of the Virgo cluster using SBF distances for 84 galaxies observed in the ACS Virgo Cluster Survey (ACSVCS; Côté et al. 2004).Thus, their sample size was much larger, and the mean random distance errors of the space-based measurements were more than a factor of two smaller, than for the data set analyzed by West & Blakeslee (2000).Five of the galaxies (including NGC 4365) were found to be members of the W group, about 6 Mpc beyond the main cluster.Excluding these objects, the rms depth of the remaining 79 galaxies, including members of the A, B, and C subclusters, was estimated to be just σ d = 0.6 ± 0.1 Mpc.Analysis of the three-dimensional positions of these galaxies showed that they defined a triaxial distribution with axis ratios 1:0.7:0.5 inclined about 30 • from the line of sight.SBF distances measured from ACSVCS data for several additional Virgo galaxies, and a slightly revised SBF calibration, were presented by Blakeslee et al. (2009).Despite the enormous advance over earlier studies of Virgo's structure, the targeted ACSVCS covered only a small fraction of the cluster area, missing many galaxies suitable for SBF distance measurements; for instance, the ACSVCS sample did not include any galaxies in the M or W clouds. Ideally, one would like high-quality imaging over the entire cluster.
The Next Generation Virgo Cluster Survey (NGVS Ferrarese et al. 2012) is a deep, multi-band imaging survey with the MegaCam instrument on the Canada-France-Hawaii Telescope (CFHT) covering an area of 104 square degrees and extending out to the virial radii of the A and B subclusters.The NGVS has already produced numerous results on the properties, varieties, and distributions of Virgo galaxies and their globular cluster systems (e.g., Durrell et al. 2014;Ferrarese et al. 2016;Powalka et al. 2016;Sánchez-Janssen et al. 2016, 2019;Roediger et al. 2017;Liu et al. 2020;Lim et al. 2020), as well as other objects along the line of sight (e.g., Chen et al. 2013;Raichoor et al. 2014;Lokhorst et al. 2016;Fantin et al. 2017).The best seeing conditions during the acquisition of NGVS imaging were reserved for the i band, and this makes it ideal for SBF measurements.
As a first step to mapping the three-dimensional structure of Virgo from NGVS, Cantiello et al. (2018a) published SBF measurements for 89 VCC galaxies with total B magnitude B T 13 mag and provided a preliminary calibrations for M i as a function of various colors and color combinations that could be constructed from the NGVS u * , g, i, z passbands.The broadest baseline (u * −z) color provided the best overall calibration with the lowest scatter.In a few cases where u * photometry was not available, a calibration that combined both (g−i) and (g−z) was used.The sample included five galaxies in the M and W clouds, both of which lie at roughly 30 Mpc but are widely separated on the sky, and a similar number in the W group at ∼ 22.5 Mpc.
That study demonstrated the potential of the method in mapping the structure of the Virgo cluster and its surrounding groups using NGVS data.However, a more extensive sample encompassing fainter galaxies is desirable due to the significantly larger number of available targets, enabling a more detailed examination of the cluster's internal structure.Expanding to fainter galaxies necessitates refining the calibration, particularly towards the blue colors typical of low-mass galaxies.A significant advantage of NGVS lies in its ability to employ a single dataset comprehensively and consistently, covering a range of masses (and colors) from bright galaxies reaching M B ∼ −23 mag, like M 87 or M 60, to much fainter ones with M B ∼ −12 mag.This paper presents SBF measurements and calibrated distances for about 280 galaxies observed as part of NGVS.This is the largest and most comprehensive sample of SBF measurements in Virgo to date, extending to fainter and significantly more diverse galaxies than those analyzed by Cantiello et al. (2018a).The following section, §2The Galaxy Samplesection.2, summarizes the sample of galaxies that we study, while §3SBF measurementssection.3 describes the SBF measurements and §4Analysissection.4 presents the refined SBF calibration based on this larger sample and the resulting distances.We find a median relative distance uncertainty of 7.5% for bright galaxies with B T < 13 mag, and ∼ 10% for fainter galaxies (B T > 13 mag).The final section, §5Results and perspectivessection.5, discusses the three-dimensional distribution of the sample galaxies and the effect of internal stellar population gradients on the SBF magnitudes, with a comparison to theoretical models.We conclude with some thoughts on the prospects for going beyond the current work.

THE GALAXY SAMPLE
As with the first SBF distance paper in the NGVS series (Cantiello et al. 2018a, NGVS Paper XVIII, hereafter referred to as Paper I), the present study is based on CFHT/MegaCam imaging data from the NGVS.Full details on the NGVS survey observations and image processing are presented in Ferrarese et al. (2012).As a brief summary, the NGVS exploits the capabilities of CFHT/MegaCam to reach 5-σ limiting magnitudes for point sources of 26.3, 26.6, 25.8, and 24.8 mag in the stacked u * , g, i, and z images, respectively, across the 104 deg 2 area of the survey.These limits are well beyond the turnover magnitude for the globular cluster (GC) luminosity function (e.g., Durrell et al. 2014), even when the GCs are superimposed on a bright galaxy background.Moreover, the median FWHM for each NGVS band (u * giz) is better than 0. 9, with a median FWHM of 0. 54 in the i band.These properties make the NGVS images extremely well suited to an SBF analysis, as most of the potentially contaminating sources can be identified and removed.
For this second NGVS paper on SBF distances, we have analysed ∼ 300 galaxies brighter than B T ≈ 19 mag (median B T = 14.5 mag) in the Virgo Cluster Catalogue (VCC) of Binggeli et al. (1985).This includes approximately half a dozen faint galaxies (B T ∼ 18 mag) that were not included in the VCC, but are available in the NGVS catalog of Ferrarese et al. (2020).Although we imposed no strict selection on morphology, we did choose galaxies that exhibit some degree of morphological regularity.Such regularity is an essential requirement for conducting a meaningful analysis of fluctuation amplitudes (see, e.g., Cantiello & Blakeslee 2023).A future paper in the NGVS series will introduce a customized morphological classification system for NGVS galaxies, including the ∼ 300 galaxies that make up our current sample.For the time being, we simply note that the sample is dominated by galaxies with "earlytype" morphologies.For example, 184 of the 215 galaxies having quality classification codes of q1 or q2 (see the following section) are classified as E or ES galaxies (usually corresponding to the VCC classifications of E, S0, dE or dS0).The remaining systems have later type morphologies but usually with the presence of a prominent bulge or bar.We point out that this requirement for some degree of morphological regularity could introduce environment-related selection effects in our sample, since galaxies within dense groups and clusters tend to have less star formation and a smoother appearance.Table 1The Galaxy Sampletable.1 lists the main properties of the galaxies included in the present work.For each galaxy, we list the VCC number from Binggeli et al. (1985) (Col. 1) its alternate and IAU name (Cols.2, 3); J2000 celestial coordinates in degrees (Cols.4, 5); total B T magnitude from the VCC (Col.6); galactic extinction (Schlegel et al. 1998, Col. 7); heliocentric velocity v h from the NASA Extragalactic Database (Col.8); morphological T type (Col.9) and its error (given in parenthesis) from Hyperleda (Makarov et al. 2014).
The SBF analysis for galaxies in the present sample is carried out as detailed in Paper I and briefly described in the following sections.

SBF MEASUREMENTS
For measuring the SBF amplitudes in the sample galaxies, we adopted the same procedures developed and described in Paper I, which are based on well established methods (Blakeslee et al. 2001(Blakeslee et al. , 2009(Blakeslee et al. , 2010;;Cantiello et al. 2005Cantiello et al. , 2007aCantiello et al. ,b, 2011Cantiello et al. , 2013;;Mei et al. 2005a,b;Moresco et al. 2022;Cantiello & Blakeslee 2023).A few minor changes were made, as detailed below.
In broad terms, the SBF distance measurement entails the following steps: a. First, we subtract the local sky background, mask bright stars and neighboring galaxies, then model the 2-D galaxy surface brightness distribution as a series of elliptical isophotes and subtract this model from the image.Large-scale model residuals are then fitted with a bicupic spline and subtracted to produce a clean residual image.In terms of the SBF definition, which involves the ratio of the second to the first moment of the luminosity function of the observed stellar population (Tonry & Schneider 1988), this process generates a model image representing the first moment of the light distribution and a residual image related to the second moment of the light distribution.
b.We then detect all point-like and extended sources (foreground stars, GCs in the target galaxy, faint background galaxies) down to a fixed signal-tonoise threshold using an automated photometry program.These sources, along with other contributors to non-SBF variance (e.g., visible dust, brighter satellite galaxies, tidal features, regions of poor model residuals), are then combined in a global mask.This masking stage, along with the previous step of galaxy modeling, subtraction, and a low-order fit to the background, follows an iterative procedure.
c. Next we obtain one or more accurate point spread function (PSF) templates for the image.This is a crucial step because stellar fluctuations are convolved with the PSF.In the Fourier domain, these fluctuations multiply the Fourier transform of the PSF (convolved with the window function of the mask; see item d below).Thus, accurate measurement of the SBF amplitude requires a robust determination of the PSF.To ensure the required reliability, we generate several PSF templates from stars in the field, normalized to unit flux.Generally, for NGVS, we find that empirical PSF templates provide more robust results (i.e., a better match to the azimuthally averaged powerspectrum) compared to model PSFs.

d.
We then analyze the power spectrum of the masked residual frame, normalized to the square root of the galaxy model.After azimuthally averaging the power spectrum, P (k), the total fluctuation amplitude corresponds to the P 0 coefficient in the fitted equation: where E(k) is the "expectation power spectrum," calculated by convolving the power spectrum of the normalized PSF with the power spectrum of the mask window function, and the P 1 term is independent of wavenumber k in the absence of correlated noise.
e.The fitted P 0 includes contributions from all astrophysical sources of fluctuation, both the stellar SBF that we seek to measure and contamination from faint sources below the detection limit.The presence of dust could potentially impact the power spectrum and, consequently, the fitted P 0 .However, the u * and g NGVS images make it possible identify the dusty regions in color maps, so that the SBF measurements can be confined to clean areas.
Because of their high frequency and radial concentration in early-type galaxies, GCs below the direct detection threshold are the main source of contamination in the P 0 measurement.The surface density of background galaxies, being lower and relatively uniform, has a less dramatic impact and is more easily constrained.The "residual power" P r is estimated from sources fainter than the detection limit (which varies with galactocentric radius) by extrapolating the fitted com-  Note-This table is available in its entirety in a machine-readable form in the online journal.A portion is shown here for guidance regarding its form and content.
Table 2. SBF, colors and distances for the 278 galaxies in our sample (Excerpt) Note-This table is available in its entirety in a machine-readable form in the online journal.A portion is shown here for guidance regarding its form and content.All magnitudes are AB mag.
Figure 1.SBF analysis images and plots for one of our sample galaxies (VCC 1146).Starting from upper left: i-band image, residual, and residual masked image (first to third panel).The black annuli in the second and third panels show the inner and outer radii of the region adopted for SBF measurements.Lower left panel: Fitted luminosity function of external sources.Filled green circles mark observational data, with a downturn at faint magnitudes due to incompleteness; the solid green line represents the best fit to the data, after accounting for incompleteness; the two components of the total luminosity function, i.e. the background galaxies and the GC luminosity function, are shown with blue dotted and red dashed curves, respectively.Lower central panel: azimuthal average of the residual image power spectrum (gray dots) and the fit obtained according to the procedure described in text (solid black curve).Lower right panel: the fitted P0 value as a function of the lowest wavenumber kstart used for the fit (the highest wavenumber used in all the fits is chosen as k end = 450).The vertical lines at 80 and 160 show a factor-of-two range in kstart over which there is little variation in the fit, and we adopt the median P0 (indicated by the dotted horizontal line) over this range; see text for details.
bined GC and background galaxy luminosity function, multiplied by the square of the source flux, from the radially-dependent detection limit down to zero flux (see Figure 1figure.1,lower left panel).The residual power P r is subtracted from P 0 to obtain intrinsic stellar fluctuations, typically denoted as P f = P 0 −P r from which the SBF magnitude is derived as m = −2.5 log P f + m zp , where the zeropoint term is m zp = 30 mag in the NGVS images.
f. Adopting a value for M from either an empirical or theoretical SBF calibration (generally based on galaxy color; see Blakeslee 2012, for a discussion), we finally obtain the distance modulus (m − M ).
Figure 1figure.1 illustrates the basic steps in the SBF measurement.The top panels show the image of the target galaxy (VCC 1146 in this example), the residual frame (step a in the list above) after galaxy subtraction, and the residual frame after masking the regions contaminated by foreground/background sources and image defects (step b).The lower panels show the fitted background luminosity function model (left panel; step e), the power spectrum of the residual frame compared to the scaled PSF power spectrum (middle panel; step d), and the fitted P 0 value as a function of the starting wavenumber k start for the fit to the power spectrum.
The upturn in the power spectrum at low wavenumber k (lower middle panel) is due to the imperfect subtraction of large-scale features in the residual frames (step a).The downturn at high k occurs because of the correlation of noise in adjacent pixels caused by the sinc-like interpolation kernel during image stacking (see Tonry et al. 1990;Cantiello et al. 2005;Mei et al. 2005a).For these reasons, the lowest (k < k start ) and highest (k > k end ) wavenumbers are excluded from the power spectrum fits.The downturn at high k from the noise correlation is straightforward to determine, and here we conservatively stop each fit at a maximum wavenumber of k end = 450.However, the fitted P 0 is more sensitive to excess power at low k, and this depends to some degree on the quality of the isophotal fit.As an example, each dot in the lower-right panel of Figure 1figure.1 shows the (logarithm of) the P 0 value obtained when starting the fit at k start in the range 0 < k start < 450.When using k start 80, the fitted P 0 values show an upward trend with decreasing k start , a consequence of the excess power due to imperfect galaxy model residuals on large scales.
Compared to Paper I, we have made a small change in the procedure for determining P 0 .In the previous work, we examined the log P 0 versus k start plots (lower right panel in Figure 1figure.1)to visually iden-tify the best values of k start for stable P 0 determinations, and averaged over a representative range.For the present measurements, we follow the approach described by Blakeslee et al. (1997).First we identify the value of k start where the log P 0 versus k start relation stops declining; call this k s,1 .We then perform a series of fits with k start in the range k s,1 ≤ k start ≤ k s,2 , where k s,2 2 k s,1 .For example, in Figure 1figure.1,we use (k s,1 , k s,2 ) = (80, 160), typical values for the galaxies we analyze.We then take the median of the fitted P 0 values from among all these fits, and estimate the error by combining the P 0 scatter in quadrature with median fit uncertainty.This approach, as compared to basing the range purely on visual inspection, results in a more homogeneous set of P 0 values for the sample galaxies.However, in practice, it negligibly changes the m i measurements presented in Paper I, with a median difference in the final m i values ∆(m i,new −m i,Paper I ) < 0.01 mag.
Another change with respect to Paper I is that we now present a quality classification code, or flag, for each SBF measurement: q1 for excellent quality; q2 for good quality; q3 for low quality.These flags are assigned based on our confidence in the result of the full SBF analysis, considering the severity of the isophotal model residuals, possible effects of bright nearby objects, the degree of masking, and the regularity of the power spectrum (which can be affected by the masking, chip defects, isophotal irregularity, etc.).
We also measure the integrated colors of the galaxies in the same annuli used for the SBF measurements.Table 2The Galaxy Sampletable.2collects the results of the SBF analysis, and the distances derived as detailed in the forthcoming section § 4Analysissection.4.In particular, we report the (u * −z), (g−i) and (g−z) colors and their uncertainties (Cols.2-4)2 ; SBF magnitude (Col.5); the preferred distance modulus and the distance in Mpc (Cols.6 and 7); the mean distance using all four SBF versus color calibrations derived here (Col.8); the area and the median radius of the annulus adopted for each galaxy in our sample (Cols.9 and 10); and the quality flag for the SBF measurement (Col. 11).
In all the figures that follow, the colors and magnitudes reported in Table 2The Galaxy Sampletable.2have been extinction corrected using the values from Schlegel et al. (1998) 3 .Figure 2figure.2plots the SBF magnitudes for the sample of galaxies reported in Table 2The Galaxy Sampletable.2versus various colors.In the panels of the figure, different shades of blue refer to different quality classification codes: i.e., dark-filled symbols for q1, light-filled symbols for q2 and open circles for q3.In all, we measured SBF distances for 278 galaxies4 .

ANALYSIS
To derive distances from SBF measurements, an accurate calibration of M i is required (step f in the previous section).In this section, we describe the calibration procedure adopted for the new sample of galaxies, discuss the differences with respect to the calibrations given in Paper I (and other relations available in the literature) and compare the empirical calibrations to those predicted by stellar population models.

Calibrating absolute SBF magnitudes
As shown in Figure 2figure.2, the relation between m i and the galaxy integrated colors (with m i becoming fainter as the color gets redder) appears even with no correction for the presence of bright galaxies in the background of Virgo, including members of the more distant substructures of the cluster, like the W or W clouds (see section § 5.1Distances and Substructuresubsection.5.1).To model the SBF versus color relation, we adopted a procedure similar to the approach of Blakeslee et al. (2009, ACSFCS-V hereafter) for the analysis of SBF measurements from HST/ACS data of Fornax and Virgo cluster galaxies.The procedure takes into account the cluster depth, as well as the mean uncertainties on the measured SBF, and rejects the outliers with a σ-clipping iteration.
In more detail, using the data in Table 2The Galaxy Sampletable.2 and plotted in Figure 2figure.2,we first selected the targets flagged as q1 or q2 and with uncertainty ∆m i ≤ 0.25 mag, to avoid targets with less reliable measurements.The selected sample contained N sel = 212 galaxies.We then fitted the m i versus color relations with a clipping algorithm to reject the outliers (mostly galaxies in the W and W clouds).For the fitting, we used the astropy (Astropy Collaboration et al. 2022) class F ittingW ithOutlierRemoval, which performs a fit according to the chosen model and clips the outliers at each iteration, until no outliers remain.We chose to use a polynomial model function and tested degrees between one and five; for the clipping, we used σ MAD , a robust estimate of the scatter derived from the median absolute deviation (MAD; σ MAD = 1.48 MAD), with a 3σ rejection parameter.
We adopted a third-order polynomial because, on average for the four colors, it gave a slightly better χ 2 ν ∼ 1.0 and lower scatter compared to lower and higher degree polynomials.This is consistent with previous studies.For example, in ACSFCS-V, a third degree polynomial proved the most adequate description of the z 850 versus (g−z) relation.
After just three iterations, the clipping and fitting procedure converged in all four colors considered.The numbers of non-rejected galaxies used in the fitting procedure were N f it = 193, 196, 192, 195 for the (u * −z), (u * −i), (g−z), and (g−i) relations, respectively.The reduced χ 2 ν in all cases was ∼ 1.0, with rms scatter of σ 0 ∼0.18 mag for the (u * −z) and (g−z) relations, and σ 0 ∼0.20 mag for the other two colors.The results of the polynomial clipping and fitting procedure are shown in Figure 3figure.3.
We adopt a two-step procedure to determine the absolute zero points of the M i versus color relations.Similar to the approaches used by ACSFCS-V and Mei et al. (2007), we start by assuming the Cepheid-calibrated mean SBF distance modulus of Virgo from Tonry et al. (2001), modified as described in Blakeslee et al. (2002) and then reduced by 0.023 mag for consistency with the revised LMC distance of Pietrzyński et al. (2019) giving a fiducial distance modulus of 31.067mag, with a systematic uncertainty of 0.09 mag (see Blakeslee et al. (2021) for further discussion).We subtract this value from the zero points of the apparent m i versus color relations to get a first estimate of the zero-point calibration.
Next, to ensure the best consistency of our distances with previously published high-quality SBF distance measurements, we refine the zero points of the fitted SBF-color relations by comparing the first-step calibrated distances with measurements for the same individual galaxies tabulated in ACSFCS-V (again adjusted for the revised LMC distance modulus).In total, we find 83 galaxies in common between our sample and with ACSFCS-V for the (u * −z) and (u * −i) calibrations, and 85 for (g−z) and (g−i).The small difference in sample size is due to the fact that some NGVS targets lack u * -band data.Using these sets of galaxies in common between NGVS and ACSFCS-V, we do final shifts in the zero points of the calibration relations so that the median difference between the ACS and NGVS SBF distance moduli is equal to zero for each color calibration.Table 3Calibrating absolute SBF magnitudesequation.4.3 provides the results of the cubic polynomial calibration equations for each color in the form: where x ≡ col − col ref , col is the color index used for the calibration, and col ref is the reference color, taken as the median color (rounded to the nearest 0.05 magnitude bin) in that index for the galaxies in common with ACSFCS-V.The absolute zero points, M i (x=0) = a 3 , are set by the comparison to ACSFCS-V, as described above.
Thanks to the large color range of the galaxies in our sample, we have characterized the intrinsic scatter for both the blue and red side of the SBF versus color relations, using col ref as blue/red separation limit.In Table 3Calibrating absolute SBF magnitudesequation.4.3 we report some other results from the fitting procedure for each color: the scatter σ 0 of the fitted curve, the reduced χ 2 ν , the intrinsic cosmic scatter in the SBF calibration σ cos , together with σ cos,blue and σ cos,red , the cosmic scatter derived for the galaxies bluer/redder than col ref .To determine the cosmic scatter 5 , σ cos , for each relation, we first estimated the depth of Virgo, σ est , adopting the same approximations used in ACSFCS-V (their eq.3): 5 The cosmic scatter definition we adopt here is the irreducible scatter in SBF magnitudes at fixed galaxy properties (specifically, the color, in our case, as we calibrate M i against the galaxy color).By definition, this represents the difference in M i one might measure from ideal observations for the SBF signal in two different galaxies with the same color, due to the intrinsic difference in the luminosity function of their stellar populations (see Tonry et al. 2000) based on the positional scatter (σ R.A. and σ Dec. ) of the sample of Virgo cluster galaxies from the VCC.After removing galaxies with B T > 18 mag or with v h > 3000 km s −1 , we find σ est = 0.12 mag.
We then estimated the cosmic scatter in the m i -color relations as where σ 0 is the observed scatter reported in Table 3Calibrating absolute SBF magnitudesequation.4.3, and σ err is the median measurement error on m i .We find σ cos ∼ 0.07, 0.06 mag for (u * −z), (g−z), respectively, and slightly larger values for (u * −i) and (g−i).Values of σ cos ∼ < 0.08 mag are consistent with the upper limits of previous estimates (Tonry et al. 2000;Blakeslee et al. 2009).Larger scatters may be explained by the wider color and magnitude ranges of our sample compared to previous studies.For example, the ACSVCS-V (g−z) calibration sample includes a total of 14 galaxies in Virgo and Fornax with (g−z) < 1.1 mag; in our sample we have ∼ 200 galaxies bluer than this.Using the same fitting procedures above, and further limiting the sample to class q1 and bright class q2 galaxies (N sel ∼ 120, N f it ∼ 110) to better match with the magnitudes and color ranges of previous SBF calibration studies, we find σ cos,red ∼ < 0.06 mag for the (g−z) and (g−i) calibrations, and σ cos,red ∼ 0.08 for (u * −z) and (u * −i).The (g−z) result agrees well with that of ACSFCS-V.Such a conservative selection, though, narrows the color range by ∼ 20% on the blue side, greatly reducing the usefulness of the M i calibrations for dwarf galaxies.
The increased SBF scatter at blue colors is expected due to the more complex age and metallicity variations in the stellar populations of low-mass, blue galaxies compared to bright, massive ones (e.g.Carlsten et al. 2019;Greco et al. 2021;Kim & Lee 2021;Moresco et al. 2022).Any contribution of very young stars especially increases the scatter of the SBF magnitude at a given color; most often, the scatter is observed to increase in moving from massive red ellipticals to bluer, low-mass dwarfs.Interestingly, in the case of (u * −z), the scatter is actually smaller on the blue side, suggesting that this broadest baseline color is more effective in characterizing the stellar populations of blue galaxies.
Combining the m i values and the colors in Table 2The Galaxy Sampletable.2,with the equations in Table 3Calibrating absolute SBF magnitudesequation.4.3, four different estimates of distance can be obtained for each galaxy in our sample.As in Paper I, we adopt the (u * −z) calibration for reference6 , using the mean of the (g−i) and (g−z) distances for the targets with no u *band data.The reference distances and distance moduli are reported in Table 2The Galaxy Sampletable.2,indicated with the ref label.In the table, we also provide the distance d mean obtained by averaging the estimates from all four calibrations (or two, in cases where no u * photometry is available).
Finally, the uncertainties reported on the distance moduli in Table 2The Galaxy Sampletable.2are the square sum of the measurement errors ∆m i , the propagation of the color uncertainty on the cubic polynomial, and the cosmic scatter σ cos (adopting σ cos = σ cos,red for galaxies with color larger than col ref , and σ cos = σ cos,blue for galaxies with colors col < col ref ).
While extensive literature exists on the distances of galaxies in our sample, a fair comparison would necessitate rescaling all available distances to a common reference calibration.Moreover, given the wealth of literature on Virgo, the number of available distance estimates likely surpasses a thousand measurements with diverse and often poorly characterized uncertainties.A comprehensive comparison, encompassing numerous galaxies with diverse distance estimates and standardizing to the same zero point, is beyond the scope of this work.
Nevertheless, we recognize the importance of comparing our SBF distances with those available from the TRGB method.The preference for TRGB distances is due to their reliability and potential for providing a Cepheid-independent calibration of the SBF method, useful for comparison to the Cepheid-based distance ladder.Adopting the TRGB distances for M 87 (VCC 1316, NGC 4486) and M 60 (VCC 1226, NGC 4472) from Bird et al. (2010) and Lee & Jang (2017), homogenized as in Blakeslee et al. (2021), we find a very close agreement between our distance moduli in Table 2The Galaxy Sampletable.2 and literature TRGB value for M 87, with ∆(SBF − T RGB) = 0.03 ± 0.13 mag, and a 1-σ match of ∆(SBF − T RGB) = 0.13 ± 0.13 mag for M 60.

Comparison with Paper I and previous calibrations
Paper I presented the detailed SBF measurement procedure and derived a calibration for the SBF versus color relation, using a sample of 32 galaxies with accurate distances from ACSFCS-V7 .Due to the relatively small color range, it was appropriate in that case to fit the color dependence of M i with a linear relation.In addition to the different fitting scheme used in Paper I, we adopt here a revised LMC distance and a richer sample of galaxies to set the zeropoint: 85 galaxies for (g−z) and (g−i), and 83 for (u * −z) and (u * −i).
Figure 4figure.4shows the difference in distance moduli for galaxies in common between ACSVCS-V and the present set of NGVS measurements, as well as the earlier Paper I values, for all four of the color calibrations derived.Note the few calibrators in Paper I toward the blue edge of the sample at (u * −z)∼ 2.2 (the red sym- As a further comparison, Figure 5figure.5 shows the two sets of calibrations (gray for the present relations; blue for Paper I) along with the difference between the two (red curve in the figure).The largest difference in the predicted M i between the relations and, therefore, in the derived distances, is observed over distinct ranges for each color -for the (g−i) and (g−z) relations toward the very red end of the color range, where the upturn in the third degree polynomial fit diverges sensibly from the linear trend; toward intermediate colors with respect to the range of validity of Paper I (the yellow shaded regions in the figure, at (u * −z)∼ 2.5 mag) for (u * −z) and (u * −i), and at the very blue tail of the (u * −i) and, most notably, the (g−i) relations.
The preference of the third degree polynomial over a linear relation at red colors cannot be robustly constrained with the present data, as the difference between the two relations starts to diverge toward colors more red than the range spanned by our sample.As already explained, our preference for the third degree polynomial is based on the χ 2 ν values and the scatter of the fits.Moreover, even for the very red objects in our sample, over the color range where both relations are valid, the separation between the fitted relations is consistent within our estimated scatter, with a higher divergence observed in the (g−i) color.
However, for the blue side of the relations, especially for the (g−i) calibration, the linear extrapolation from Paper I (blue solid lines beyond the yellow shaded areas in Fig. 5figure.5)would imply overly large mean distances for blue galaxies.Hence, a measurable difference between a linear and a third-degree fit emerges.For the specific case of the (g−i) relations, we direct the reader to Section 4.3Comparison with other blue calibrationssubsection.4.3.
In conclusion, the combination of factors reported in the above paragraphs further support our choice for the degree of the fitting polynomial.We emphasize that the distances presented here supersede those presented in Paper I.
The Carlsten et al. calibration refers to the same instrument and passbands used in this work (MegaPrime/MegaCam), hence no further passband transformation is required to compare it with our data.The Kim & Lee relation is instead derived from Subaru/Hyper Suprime-Cam data, which use a photometric system slightly different from that of CFHT.We use the equations derived by Kim & Lee from PAR-SEC isochrones (Bressan et al. 2012) and, like these authors, assume that the i-band magnitudes of the two systems are almost identical, while for the g-band a correction is required.To this aim, we invert the equation given in section 5 of Kim & Lee and derive (g−i) HSC = 1.045(g−i)CF HT + 0.006.Hence, the revised M i versus (g−i) in the CFHT passbands is: The relations are plotted in Figure 6figure.6,along with our best-fit equation.In the figure, our (g−i) calibration is illustrated with an error region rms = 0.15 mag, obtained from the quadrature sum of σ cos (g−i) and the median measurement error in m i .In spite of the differences between the three equations -for example, in the approach used for data analysis (e.g., the way the P r correction is handled), in the standard candle used to anchor the calibration (TRGB for for the literature relations; Cepheids for our relation), or in the typical targets used -the relations overlap within their rms scatter over the range of validity in color.
In Figure 6figure.6we also show, with a blue line, the extrapolation of the Paper I M i versus (g−i) relation derived for colors (g−i) ∼ > 0.8 mag.The plot clearly evidences the systematically too bright SBF amplitudes extrapolated for (g−i)≤ 0.7 mag, independent of the calibration with which they are compared.As discussed in the previous section, this evidence supports our preference for a higher degree polynomial over the wide color range inspected in this work, compared to the linear relation adopted in previous studies, like Paper I.

SBF calibrations vs simple stellar population models
As discussed in the previous section, obtaining a distance from m requires an estimate of M .This typically involves an empirical calibration that parameterizes the correlation between M and the properties of the unresolved stellar population in the galaxy via its integrated color.However, theoretical predictions of M from stellar population synthesis models also represent a viable alternative for deriving absolute SBF versus color relations.
There have been many efforts in the literature to predict absolute SBF magnitudes using stellar population models (e.g., Worthey 1993;Blakeslee et al. 2001;Raimondo et al. 2005;Cook et al. 2020;Chung et al. 2020;Greco et al. 2021).SBF magnitudes are heavily weighted toward the brightest evolutionary stages, and theoretical calibrations of the SBF method from different sets of models can differ by several tenths of a magnitude or more, particularly in near-infrared passbands where luminous stars in fast evolutionary stagesespecially the thermally pulsating asymptotic giant branch (e.g., Raimondo 2009)-are still inadequately modeled.Nevertheless, despite the disagreement among model predictions (or, arguably, because of it), SBF magnitudes, gradients, and distance-independent SBF "fluctuation colors" (Cantiello et al. 2005(Cantiello et al. , 2007b;;Jensen et al. 2015;Rodríguez-Beltrán et al. 2021) hold great promise as tools for scrutinizing the properties of integrated stellar populations within galaxies, complementary to more traditional spectrophotometric indicators.
In this section, we focus on the role of model SBF predictions in illuminating the physical properties underlying the empirical calibration equations derived in the previous section.With some exceptions (e.g., Biscardi et al. 2008), the vast majority of published SBF distances rely on empirical calibrations (Tonry et al. 2001;Blakeslee et al. 2001Blakeslee et al. , 2009;;Cantiello et al. 2018b;Jensen et al. 2003Jensen et al. , 2021)).In this section, we present a comparison of our empirical SBF-color relations with predictions from the Teramo-SPoT simple stellar population (SSP) models (Cantiello et al. 2003;Raimondo et al. 2005;Raimondo 2009).
The SPoT SBF models are based on a stellar population synthesis code first described and tested observationally by Brocato et al. (1999Brocato et al. ( , 2000)).The code is optimized to reproduce both the observed distribution of stars in the color-magnitude diagrams and the integrated properties (color and spectral indices, energy distribution, etc.) of star clusters and galaxies.A detailed description of the SPoT models and ingredients can be found in literature cited above.The updated color and SBF models used here, and in Paper I, are reported in Table 4SBF calibrations    Interestingly, the empirical calibration for blue galaxies with M i −1 mag matches well with the set of SSP models with [Fe/H] = −0.35dex for the full age range from 1 to 14 Gyr.However, we do not expect that the dominant stellar component in the bluest sample galaxies is, in general, as young as 1 Gyr, since this this would likely be associated with an irregular galaxy morphology and/or the presence of dust that would complicate analysis of the SBF signal.As noted above, target galaxies were selected partly on the basis of having regular, "clean" morphologies.Thus, these blue galaxies are likely to be somewhat older and more metal poor than the bluest models shown here.Based on previous SPoT simulations, SSP models for older ages and [Fe/H] < −0.5 dex are also expected to cover the blue/bright side of the curves in Figure 7figure.7 (see, e.g., Cantiello et al. 2003).New SSP models covering a wider range of [Fe/H] values and new photometric systems (like JWST, Euclid) are currently being developed by our team.
We emphasize that the comparisons shown here are based on SSP models derived from stellar evolution theory (e.g., stellar tracks, initial mass function, stellar atmospheres) combined with Monte-Carlo simulations.This demonstrates the robustness of the SBF technique as a distance indicator, in that it can be calibrated using purely theoretical arguments.It also highlights the usefulness of SBF measurements for studying the properties of unresolved stellar populations.We further discuss SBF as a stellar population tracer in §5.3SBF versus color gradientssubsection.5.3.

RESULTS AND PERSPECTIVES
In this section, we use our new catalogue of SBF distances to examine the three-dimensional structure of the Virgo cluster.We include a discussion of possible SBF gradients in our program galaxies, and conclude with a summary of this work as well as some thoughts on areas for future research.

Distances and Substructure
Using the reference distances, d ref , reported in Table 2The Galaxy Sampletable.2,we examine the distribution of the galaxies in our sample, both spatially and in the velocity-distance diagram.We restrict the analysis to only those galaxies with q1 and q2 quality codes.
The left panel of Figure 8figure.8shows position on the sky for our sample galaxies, with the outline of several Virgo cluster substructures identified in Boselli et al. (2014).In the figure, the half-virial radii for the A and B subclusters are from McLaughlin (1999) and Ferrarese et al. (2012), while Boselli et al. defined the other radii based on constraints to the angular distance from the overdensity peak.
For each substructure, we selected candidate member galaxies using sky coordinates, distances, and recession velocities with the following membership criteria.For all substructures, except the M cloud and LVC, we: i) selected all galaxies projected within the substructure area, adopting the central positions and radii reported in Table 5Distances and Substructuresubsection.5.1 (Cols.1-4); ii) for the selected sources, we derived the median and the standard deviation (from the median absolute deviation) of the distances and the recession velocities reported in Table 2The Galaxy Sampletable.2;iii) rejected from the selected sample all sources that are more than 3-σ outliers with respect to the mean distance and velocity.For the distant M cloud (where we can identify only a single member in our sample), in addition to angular projection proximity, we required d ref > 25 Mpc.For membership in the LVC, we required projection proximity and v h < 800 km s −1 (Boselli et al. 2014).
The substructures and their associated galaxies are color coded in the left panel of Figure 8figure.8,as labeled.Sample galaxies that fall outside the boundaries of all substructures ("Virgo field") are represented by light-gray circles.The right panel of the figure shows a Hubble diagram of the substructures including Virgo field galaxies, using the same color coding as in the left panel.As in Paper I, there is no obvious distinction in the Hubble diagram among the subclusters around M 87, M 49 and M 60, although both the W and W clouds do appear distinct from the main cluster.
The transitional region between the main Virgo cluster and the W /W clouds is better appreciated in the three-dimensional plots shown in Figure 9figure.9.In this figure, we see a sequence of about a dozen galaxies lying along a 'filament' from the core of subcluster B toward the W cloud, and covering the W region.This feature may be a projection resulting from sampling effects of the so-called W-M sheet (see Kim et al. 2016;Castignani et al. 2022, and references therein).
The main properties of the substructures inspected here are summarized in Table 5Distances and Substructuresubsection.5.1, which contains the substructure identifier, coordinates, radius from Boselli et al. (2014), the number of galaxies in our sample that we identify as members (N str ), median velocity, median distance, and estimated rms depth after correcting for expected scatter from distance errors (with the bootstrap uncertainty given in parentheses).Although we formally calculate a depth of ∼ 4 Mpc for the W cloud, we note that this structure is prone to contamination from both the richer foreground subcluster B and the background W cloud.It may be that a more useful physical definition of W would be the giant elliptical NGC 4365 at ∼ 23 Mpc and its close satellites.In that case, the group would be a much more compact knot within the filament-like W cloud extending out to beyond 30 Mpc.As previously mentioned, the SBF analysis favors galaxies with smooth surface brightness profiles.This preference may result in early-type galaxies, with smoother profiles, being measured across a wider range of distances than later type galaxies.Additionally, smooth early-type galaxies tend to reside in the denser regions of a cluster (e.g., Dressler 1980).The combination of these factors implies that, at a fixed image depth and quality, more reliable SBF distances would be measured toward the densest substructures.Hence, the properties of the structures reported in Table 5Distances and Substructuresubsection.5.1, especially the depth, should be interpreted in light of this bias.
Table 5Distances and Substructuresubsection.5.1 also reports the number of galaxies detected in the NGVS survey around the projected area of each subcluster region using the catalogs from Ferrarese et al. (2020) and Ferrarese et al. (2024), adopting the same redshift properties used to identify galaxies in the substructure, as well as the numbers without any redshift constraint (these numbers are given in parentheses in the N str column).While it would be possible to incorporate additional Virgo galaxy distances from other methods (e.g., Gavazzi et al. 1999), as noted in Section 4.1Calibrating absolute SBF magnitudessubsection.4.1, the combined heterogeneity of literature samples and the typically much larger scatter from distances derived using alternative indicators (e.g., both the Tully-Fisher and fundamental plane methods have distance uncertainties ∼ 20%, larger than the sizes of most of the structures we have discussed), as compared to the homogeneity of the NGVS sample and the small intrinsic scatter of SBF distances (Sec.4.1Calibrating absolute SBF magnitudessubsection.4.1), discourage us from combining our measurements with literature data for this analysis.

A Closer Look at Cluster A
In examining the spatial distribution of the galaxies in our sample, we noticed an apparent gap in distance near 18 Mpc.The feature is most pronounced in the distribution of distances for galaxies grouped into the A subcluster. Figure 10figure.10shows the histograms for all of the identified structures with more than one member.The bimodality in the histogram of subcluster A is clear; we used the Gaussian Mixture Modeling (GMM) code of Muratov & Gnedin (2010) to quantify this impression.The GMM code compares the goodness of fit for single and double Gaussian models of the unbinned distribution, then uses bootstrap resampling to assess the uncertainties.In the case of subcluster A, it finds that the distances prefer a double Gaussian model with ∼ 99.9% confidence.For the preferred double Gaussian model, the two peaks occur at 16.5 ± 0.1 and 19.4 ± 0.2 Mpc, with dispersions (uncorrected for measurement error) of 1.0 ± 0.1 and 0.5 ± 0.1 Mpc, respectively.According to GMM, the nearer peak comprises 82% of the sample.Based on the appearance of the histogram, the strong preference for two Gaussians is not surprising, but there is, of course, no reason to expect that distances should follow a Gaussian distribution.
A more meaningful question is: how significant is the gap in the distance distribution?In an attempt to address this, GMM calculates the separation of the fitted peaks in units of the combined dispersion, and assesses the significance of the gap based on this statistic.Using this approach, GMM finds the gap is significant with 90% confidence.However, there is again an underlying assumption of Gaussianity in this analysis.Finally, the software package provided by Muratov & Gnedin (2010) also includes code to calculate the "dip" statistic of Hartigan & Hartigan (1985), which provides a measure of the significance of multimodality based on the maximum difference of the sorted empirical distribution with respect to the nonparametric unimodal distribution that minimizes this maximum difference.According to this test, we cannot rule out a unimodal distance distribution for subcluster A at even 1σ significance.Muratov & Gnedin (2010) find a similar result for the distribution of metallicities in Galactic globular clusters, and conclude the dip test is "less powerful" than GMM.However, it is robust against the assumption of Gaussianity.
In any case, the distance distribution of galaxies in subcluster A is certainly suggestive of bimodality.We performed tests with alternative distance calibrations for M i , such as using (g−i) instead of (u * −z) for the stellar population dependence, and adopting a linear rather than a polynomial relation.However, the apparent bimodality remains, with similarly high significance 99.8% from GMM, and low significance 50% from the dip test.Thus, the double-peaked structure is not the result of the adopted calibration relation.
The most luminous galaxy in the secondary distance peak is M 86 (VCC 881) at 19.1 ± 1.1 Mpc and with a velocity v h = −224 km s −1 .While there are a few other galaxies (VCC 828, 997, and 1414) with distances of 19 to 20 Mpc, v h 500 km s −1 , and projected locations within 0.5 Mpc of M 86, there is no obvious spatial concentration of the galaxies in the secondary peak around this bright galaxy.Figure 11figure.11compares the distances, magnitudes, sky locations, velocities, and colors of galaxies in the primary and secondary distance peaks of subcluster A, where we use 18 Mpc as the dividing line.M 86 is marked with a square box in the figure.Most of the galaxies in the secondary peak are dwarfs  5Distances and Substructuresubsection.5.1).The structures are color coded and labelled, with the size of the symbol scaling with galaxy magnitude.Note that the circles provide only rough indications of the structures, and in some cases represent upper values.Field galaxies, not associated with any of the substructure analysed, are indicated by the gray circles.Right: Hubble diagram for these same galaxies.Symbols and color coding are the same as in the left panel.
with g magnitude m g > 15 mag, but overall the galaxies are not distinguished from the rest of the subcluster A galaxies by any parameter other than distance.We conclude that the apparent bimodality in the distances within subcluster A is intriguing, but needs confirmation and further exploration with an even larger sample of high-quality distances.

SBF versus color gradients
In addition to the main annulus used for our SBF analysis, for every galaxy in our sample, we measured SBF and colors in three to five smaller and distinct circular concentric annuli, to analyze the radial behaviour of integrated stellar population properties using combined SBF and color gradients.In examining these radial gradients, we limited our sample to galaxies that have either q1 or q2 quality flags, and for which the SBF and color gradients are monotonic.
The left panels of Figure 12figure.12show SBF amplitudes versus color data for the galaxies where the uncertainty on the fitted slope, ∆M i /∆color, is smaller than one half of the expected pure age gradient, as explained below.Using the SPoT SSP models in Table 4SBF calibrations vs simple stellar population modelstable.4,we estimated the M i versus color gradient as follows: we first fix the metallicity, [Fe/H], and, for the entire set of ages available, derive the slope of the M i -color relation for each of the three [Fe/H] available.The dashed lines in the upper left panel of Figure 12figure.12indi-cate sequences in age at fixed [Fe/H].For the two colors considered, the mean slopes derived are: Similarly, we derived the slopes of M i versus color at fixed ages.In the lower left panel of Figure 12figure.12,the dashed lines connect the predictions from the SPoT SSP models for coeval (fixed age) populations with differing values of [Fe/H].The median slopes of the fitted linear equations are: Rather than the exact values obtained from SSP models, what we highlight here is the evidence that SBF versus color trends driven by [Fe/H]variations at fixed age are predicted to have a slope two times steeper than those at fixed metallicity [Fe/H].In the left panels of Figure 12figure.12,we illustrate the direction of age gradients (indicated by the red arrow labeled "Fixed [Fe/H]") and metallicity gradients (shown by the blue arrow labeled "Fixed Age").   ) is marked by a black square, while three galaxies near it on the sky and with similar velocities are marked by blue circles.Apart from these few galaxies, there is no obvious grouping in position and velocity of the secondary peak galaxies around M 86.The two distance subgroups fall along the same color-color relation, so the putative bimodality does not appear to be a consequence of stellar population differences.
In the right panels of Figure 12figure.12,we plot the histograms of the SBF-color slopes derived for all class q1 and q2 galaxies with monotonic radial gradients in SBF magnitude and color, without the limitation on maximum error of the slope (adopted in the left panels for sake of clarity).In the panels, we represent with green histograms the sample of ∼40 targets thus selected and the mean slopes derived from SSP models, with the same color coding as in the left panels.
From the histograms of gradients for both colors, we find that the median for the observed data lies between the two fiducial lines, with a slight preference for smaller slopes (i.e., age-driven, fixed [Fe/H]) compared to larger slopes (fixed Age, [Fe/H]-driven).Additionally, we assessed the gradient differences between bright and faint galaxies, using a threshold at B T = 13 mag.Fainter galaxies tend to exhibit slightly larger gradients, although the statistical significance is weak.
Of course, most galaxies likely have both age and [Fe/H] gradients to some degree.Nevertheless, it is interesting to note that in the lower-left panel of Figure 12figure.12,most of the multi-annuli data (thick broken lines with different colors) have a ∆M i /∆(u * −i) < 2, indicative of significant age gradients.The biggest exception is VCC 1720 (dark-red line, at (u * −i)∼ 2.3 mag and M i ∼ −0.9 mag), which appears to align with the model sequences at fixed age, with a ∆M i /∆(u * −i) = 3.1 ± 0.1.For the same galaxy, we find ∆M i /∆(g−z) = 5.9 ± 0.5, more than twice as steep as the median slope for the sample in the upper-right panel.

Perspectives
Taking advantage of the depth and image quality of the NGVS dataset, we have presented SBF measurements and distances for 278 galaxies within the 104 deg 2 survey footprint.This sample of distances is several times larger than in any previous SBF study of the Virgo region.These results highlight the great potential of the SBF method for deriving high-quality distances for large sample of galaxies in upcoming wide-field optical and near-IR imaging surveys.
With a 5-σ point-source depth of ∼ 25.8 mag and i-band seeing ≤ 0. 7, the NGVS data are about one magnitude shallower than what is expected from Rubin/LSST.Since the distance to which SBF can be reliably measured is often limited (for bright galaxies) by the ability to remove contaminating sources, the greater depth of the LSST data implies a distance limit potentially several times larger than in this study, for which the maximum measured distance is ∼33 Mpc.Similarly, the Euclid satellite, now in operation, will have an H-band point-source limit of ∼ 24 AB mag, with an (undersampled) PSF of ∼ 0. 3.This combination of angular resolution, depth, sky coverage, and the benefit of minimal impact from dust, will be ideal for measuring SBF distances for large samples of nearby galaxies, complementing those measured in the optical with Rubin/LSST.
On slightly longer timescales, the Nancy Grace Roman Space Telescope will be of enormous interest for SBF measurements.Roman will have the same aperture as the Hubble Space Telescope and similar resolution, but a ∼ 100 times larger field of view and superior sensitivity at infrared wavelengths.With a 5σ point-source depth of 28 mag in 1 hr in the J and H bands, Roman should deliver phenomenal survey depth and breadth, making it the ultimate machine for churning out SBF distances in abundance (e.g., Greco et al. 2021;Blakeslee et al. 2023).Notably, all three of the above-mentioned facilities will observe the sky in multiple passbands, making it possible to derive exquisite calibrations for the SBF distances, as well as to measure SBF colors (i.e., differences of SBF magnitudes in multiple bands), a unique tool for probing stellar populations.
To realize the full benefit of these deep, wide-area, multi-band surveys, a new paradigm for SBF measurements will be essential.Efficient and robust automated procedures must be developed for measuring SBF amplitudes and distances in large numbers of galaxies without requiring extensive human intervention.We are presently developing such automated SBF pipelines, taking advantage of the large, well-tested sample of SBF measurements that we have presented here based on the multi-band NGVS observations.
Finally, in addition to the potential of wide-area surveys for measuring unprecedented numbers of SBF distances across the sky, the large apertures, high nearinfrared sensitivities, and excellent spatial resolution of JWST and future AO-assisted instruments on ELTs, will allow the SBF method to be pushed to unprecedented distances for targeted observations of individual galaxies.Once calibrated from the TRGB method, this will provide a fully independent distance ladder rivaling the more familiar one based on Cepheids and Type Ia supernovae.Altogether, the outlook for the SBF technique as a distance and stellar population indicator is excellent.

Figure 3 .
Figure 3. Illustration of the fitting and 3σ clipping procedures for the different colors.The 212 targets with code q1 and q2 measurements are used.The open circles show the sources rejected after the clipping iteration.The polynomials fitted to data are shown as dotted green (linear fit), dashed red (second degree) and solid gray (third degree) curves.

Figure 4 .
Figure 4. Distance modulus differences between the present NGVS and the ACSFCS-V distances (∆(m−M ) = (m−M )NGV S − (m−M )ACSF CS−V ), used to set the zeropoint of the calibration equations.Red three-pointed stars mark the sample used in Paper I, while the yellow shaded area highlights the color range of the data used in Paper I. Filled blue dots mark the new calibration sample.The vertical lines refer to the range of validity of the revised calibration equations presented in this work.The color interval for the present sample extends toward bluer color than the sample of calibrators because of the fitting approach adopted (see Section 4.1Calibrating absolute SBF magnitudessubsection.4.1).

Figure 5 .
Figure 5.A comparison of the calibration equations derived in this work (gray curves and area) and the ones presented in Paper I (blue solid lines and shaded area, extrapolated outside the yellow shaded area).The yellow shaded regions and vertical dashed blue lines are as in Figure 4figure.4.The red curve shows the difference between the two calibration equations (Paper I minus present).The green horizontal line is the zero-difference level, reported for sake of clarity.
vs simple stellar population modelstable.4.The comparison between empirical data and SSP model predictions is shown in Figure 7figure.7.The observed galaxy colors generally match the sequences expected from the SSP models in the metallicity range −0.35 ≤ [Fe/H] ≤ +0.4 dex, for ages between 1 and 14 Gyr.Each panel in the figure shows one of the SBF

Figure 6 .
Figure 6.A comparison of absolute SBF calibrations for the relation derived in this work (gray solid curve), and the calibrations from Carlsten et al. (2019, green dot-dashed curve) and Kim & Lee (2021, red dashed curve).The shaded areas show the ±1σ regions for each calibration within its range of validity.The solid blue line shows the calibration equation from Paper I extrapolated toward the wider color interval of this paper.

Figure 7 .
Figure 7.Comparison of the observed and model SBF-color relations.The gray lines show the third degree polynomial fits derived as described in the text.Shaded blue circles show the observational dataset used to obtain the calibration itself.The magenta, red and green symbols show the SPoT SSP models for three different metallicities ([Fe/H]=−0.35,0.0 and +0.4,respectively).Ages in the range t=1-14 Gyr are shown with different symbol sizes (smaller for younger ages).

Figure 8 .
Figure8.Left: map of the sky positions of the galaxies in our sample.Only class q1 and q2 galaxies in Table2TheGalaxy Sampletable.2are plotted.The possible members of each substructure are identified through our d ref measurements, recession velocities, and positions and radii (see Table5Distancesand Substructuresubsection.5.1).The structures are color coded and labelled, with the size of the symbol scaling with galaxy magnitude.Note that the circles provide only rough indications of the structures, and in some cases represent upper values.Field galaxies, not associated with any of the substructure analysed, are indicated by the gray circles.Right: Hubble diagram for these same galaxies.Symbols and color coding are the same as in the left panel.

Figure 9 .
Figure 9.The three-dimensional distribution of the sample of galaxies in Figure 8figure.8,shown here from different vantage points.The symbol size scales with galaxy luminosity, and symbol color is coded according to the actual (g−z) color.The locations of the W and W clouds are shown in several panels.Additional representations of the three-dimensional distribution are shown in the Appendix.

Figure 10 .
Figure 10.Histograms of galaxy distances within each of the substructures listed in Table 5Distances and Substructuresubsection.5.1.Subcluster A (top left panel) shows an apparent bimodality in its distance distribution with distinct peaks at 16.5 ± 0.1 and 19.4 ± 0.2 Mpc.

Figure 11 .
Figure 11.A closer look at subcluster A. Plots of distance versus magnitude, right ascension and declination, distance versus velocity, and (g−i) versus (u * −z) color for galaxies in the main d ∼ 16.5 Mpc distance peak (open red circles) and the smaller d ∼ 19.4 Mpc distance peak (filled red circles) of subcluster A, where the separation is taken as 18 Mpc.Other Virgo galaxies in our sample are shown as smaller gray points.The giant elliptical M 86 (VCC 881) is marked by a black square, while three galaxies near it on the sky and with similar velocities are marked by blue circles.Apart from these few galaxies, there is no obvious grouping in position and velocity of the secondary peak galaxies around M 86.The two distance subgroups fall along the same color-color relation, so the putative bimodality does not appear to be a consequence of stellar population differences.
M.C. and G.R. acknowledge INAF for supporting the project Gravitational Wave Astronomy with the first detections of LIGO and Virgo experiments (GRAWITA, PI: Enzo Brocato).M.C. gratefully acknowledges the invaluable professional support provided by Sabrina Ciprietti.This paper is based on observations obtained with Mega-Prime/MegaCam, a joint project of CFHT and CEA/IRFU, at the Canada-France-Hawaii Telescope (CFHT), which is operated by the National Research Council (NRC) of Canada, the Institut National des Sciences de l'Univers of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii.This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency.
made use of the NASA/IPAC Extragalactic Database (NED), which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology. 17 Astropy (Astropy Collaboration et al.   2013, 2018), Source Extractor(Bertin & Arnouts 1996), TOPCAT(Taylor 2005)

Figure 12 .
Figure 12.Left panels: M i versus color diagrams for (g−z) (upper panel) and (u * −i) (lower panel).The gray line and area show the fits of Table 3Calibrating absolute SBF magnitudesequation.4.3 and the scatter.The dashed lines mark the position of SPoT SSP models: connected by [Fe/H] homogeneity (upper panel) or by age homogeneity (lower panel).The thick coloured lines show the SBF-color sequence for a selection of galaxies.The blue (red) arrow indicates the mean slope for SBF versus color derived from SSP models with varying [Fe/H](age) at fixed age (metallicity).The direction of the arrows indicates the increase of Age or [Fe/H].Right panels: histograms of the slopes of the SBF-color relations for the selected sample.The shaded blue and red areas and lines mark the position of the mean slopes derived from the SPoT SSP models for pure [Fe/H] and for pure age variations, respectively.

Figure 13 .
Figure13.Same as in Figure9figure.9, but limited to the sample of galaxies 12.5 ≤ d(M pc) ≤ 22.The distinctive structure of the cluster, elongated along the line of sight, is clearly evident in this representation.

Figure 14 .
Figure14.3D distribution of the galaxies associated with the subclusters, groups, and clouds using the same color-coding adopted in Figure8figure.8.In this view, a filamentary structure appears to connect Cluster B and the W and W clouds.

Figure 15 .
Figure 15.Same as in Figure 14figure.14but limited to the sample of galaxies 12.5 ≤ d(M pc) ≤ 22.

Table 1 .
List of targets with SBF measurements (Excerpt)

Table 3 .
Coefficients and results of the fits M i = a0x 3 + a1x 2 + a2x + a3, where x ≡ col − col ref , and col is the photometric color used for the particular calibration.

Table 4 .
Integrated colors and SBF magnitudes from the Teramo-SPoT Simple Stellar Population models.

Table 5 .
Main characteristics of Virgo cluster substructures.