© The Authors 2025

1. Introduction

An accurate determination of distance is essential for deriving almost all the key properties of cosmic structures and objects, such as mass, luminosity, and size. Furthermore, resolving the 3D structure of galaxies within clusters requires distance measurements for individual targets with a precision better than the ratio between the depth of the cluster and its mean distance, typically corresponding to uncertainties of less than ∼0.1 mag in distance modulus (5% in distance) in the local Universe (Mei et al. 2007; Blakeslee et al. 2009; Cantiello et al. 2018; Rodríguez-Beltrán et al. 2024). Only a few distance indicators can achieve such precision. Among them, one that offers both the broad applicability necessary to resolve the depth of galaxy clusters out to the Hubble flow and the capability to be applied to a significant fraction of galaxies within the same environment is the surface brightness fluctuation method (SBF). For a detailed and updated description of the method, we refer to the review from Cantiello & Blakeslee (2023).

Since the first definition of the SBF technique by Tonry & Schneider (1988), a major observational effort has been made to improve the accuracy and precision of the method both as a distance indicator (e.g., Tonry et al. 1997, 2001; Cantiello et al. 2007, 2018, 2024; Mei et al. 2007; Blakeslee et al. 2009; Jensen et al. 2021) and as a diagnostic for the properties of the unresolved stellar populations in galaxies (e.g., Buzzoni 1993; Worthey 1994; Cantiello et al. 2003; Raimondo et al. 2005). The precision of the method depends on the quality of the measurements (image quality, background, extinction, and photometric uncertainty) and the wavelength-dependent scatter that arises from variations in the stellar populations of the galaxies themselves. Overall, in the optical passbands, the errors of the measurements can achieve statistical uncertainties on the SBF magnitude of ≲0.1 mag for ground-based observations (e.g., Cantiello et al. 2024) and ≲0.05 mag for space-based observations (e.g., Blakeslee et al. 2009). The SBF in the near-IR has some advantage over optical bands, as the amplitude of the fluctuations is stronger, the impact from undetected sources is less dramatic than in optical bands, and the effect of residual dust is negligible in most cases. Recent studies based on HSC/WFC3 data achieved SBF measurement precisions of ∼0.08 mag in the best cases (Jensen et al. 2021), not including the systematic uncertainty in the zero point calibration.

These characteristics make SBF a promising distance indicator for the forthcoming generation of wide-area imaging surveys such as Rubin, Euclid, and Roman, which will image many thousands of galaxies suitable for SBF measurement. In particular, the Legacy Survey of Space and Time (LSST; Ivezić et al. 2019) of the Vera C. Rubin Observatory will scan about 18 000 deg² of the southern sky, reaching an i-band depth of ≈26.5 for the 10 year co-added images with a median full width at half maximum (FWHM) of $\sim 0\stackrel{″}{.}6$, which will allow us, in principle, to estimate distances with SBF up to ∼50 Mpc, and possibly farther for selected objects (Cantiello & Blakeslee 2023). This represents an important opportunity to constrain the 3D structure of the local Universe. Considering the large amount of data that will be produced by Rubin and other deep surveys, it is essential for the SBF analysis to be carried out in a fast, robust, and automatic way. However, the existing SBF measurement procedures require extensive user intervention.

In this work, we present the first results from a new pipeline for the measurement of the SBF signal, specifically developed to meet the requirements of speed and automation. We test the procedure on the Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP) survey data, a precursor to LSST, and with the recent measurements from the Next Generation Virgo Cluster Survey (NGVS Ferrarese et al. 2012; Cantiello et al. 2024). The paper is organized as follows. In Sect. 2 we describe the data and the HSC-SSP survey. In Sects. 3 and 4 we present the SBF measurement procedure and then compare our measurements to previous results and to predictions from stellar population models. Our conclusions are presented in Sect. 5.

2. Data and sample selection

This study is mainly based on imaging data from the second data release (PDR2) of the HSC-SSP. The primary motivation for using HSC data lies in its intrinsic similarities with LSST. Both HSC-SSP and LSST are designed to conduct deep, wide-field optical surveys using similar photometric bands (grizy for HSC and ugrizy for LSST) and, most importantly, the data reduction pipelines developed for LSST are built upon the same framework as the HSC science pipelines (Ivezić et al. 2019). The accuracy of the SBF signal also depends on the data reduction strategies adopted (Cantiello et al. 2005; Mei et al. 2005). Thus, HSC-SSP data serve as an exceptionally valuable resource for the optimization and the automation of the SBF measurements procedure for the LSST dataset.

Full details on the HSC-SSP survey strategy and image processing are provided by Aihara et al. (2019). We briefly summarize the specifics of the survey. The Subaru Telescope is a wide-field optical-infrared imaging telescope with a 8.2 meter aperture, located at the Mauna Kea Observatory in Hawaii. The telescope is equipped with the 870 megapixel Hyper Suprime-Cam (Miyazaki et al. 2018), which covers a 1.5 deg² field-of-view with a pixel scale of $0\stackrel{″}{.}17$. The HSC-SSP multiband imaging survey was designed to map a fraction of the sky in three layers: Wide (1400 deg²), Deep (27 deg²) and Ultradeep (3.5 deg²). For this work, we chose the PDR2 Wide layer, which covers ∼1100 deg², reaching a 5σ limiting magnitude for the co-added images of 26.6, 26.2, 26.2, 25.3 and 24.5 mag in the g, r, i, z, and y bands, respectively, with a median i band FWHM of $0\stackrel{″}{.}58$. We used the stacked images from the PDR2 instead of the more recent and spatially complete PDR3 because the background subtraction strategy adopted for the PDR3 images can severely oversubtract the galaxy halo for bright and extended galaxies (Aihara et al. 2022), potentially biasing the SBF measurement, which relies on an accurate model of the galaxy light distribution. The images obtained from the HSC-SSP archive were calibrated and normalized with identical zero-point in all bands (m_zp = 27). The magnitude of the sources was derived as m = −2.5log₁₀(f)+m_zp, where f is the flux measured on the image.

From the PDR2 sample, we identified a subset of galaxies with existing SBF measurements in the literature, focusing on large galaxy catalogs (Tonry et al. 2001; Blakeslee et al. 2009; Jensen et al. 2015). At the time the cross-match of the PDR2 archive with such catalogs was performed, we found five matches for galaxies observed both in the i and g band, all of which are from Tonry et al. (2001, hereafter Ton01). We chose the i and g bands for coherence with the existing literature, using the i-band for SBF measurements, and the combined bands to get the g − i color required for estimating the absolute SBF magnitude (see Sect. 3.6). Ton01 presented SBF distances to ∼300 galaxies and 20 nearby groups. The five selected galaxies are distributed across three distinct sky regions, shown in Fig. 1. We further expanded our sample by visually inspecting the ≈1.5 × 1.5 deg area around the galaxies to identify additional candidates potentially associated with the reference galaxies and suitable for SBF analysis. In total, we identified 31 sources, two of which are relatively bright early-type galaxies while the remaining are fainter dwarf-like galaxies, all with no previous distance available from the literature. The majority of the dwarf galaxies identified are also observed in the MATLAS survey (Duc et al. 2015; Habas et al. 2020; Poulain et al. 2021). The expanded sample offers the opportunity to assess the performance of our SBF measurement code on fainter objects. While we did not impose strict morphological selection criteria, we prioritized galaxies with smooth morphology. For this reason, our sample is dominated by galaxies with early-type morphologies. The remaining systems are later types, but usually with the presence of a prominent bulge or a region showing a smooth profile. This regularity is key to observe an effective SBF signal for distance estimates, as early-type galaxies with smooth profiles minimize artificial residuals arising from the modeling procedure (see Sect. 3.2) and limit the contribution of young stellar populations that could bias the SBF calibration. (e.g., Moresco et al. 2022). Table A.1 presents the properties of the galaxies included in our sample.

Fig. 1. Spatial distribution of the galaxies in the three fields we studied. The color represents their velocity in respect to the CMB frame. Small circles represent dwarf galaxies, while bright galaxies are represented with bigger size and are labeled. The box size is 1.3 × 1 deg and is the same for all the panels.

In addition to the comparison with the Ton01 measures, we also ran our new pipeline on data with more recent ground-based SBF measurements. As is demonstrated in Cantiello et al. (2018), the roughly thirty-year-old SBF survey by Tonry & collaborators contains some outliers, especially for the most distant targets in their sample. For this comparison, we included additional data from the Next Generation Virgo Cluster Survey (NGVS, Ferrarese et al. 2012), which were previously studied in Cantiello et al. (2024). The NGVS survey obtained ugriz band imaging data of 104 sq. degrees in the Virgo cluster using CFHT MegaCam data, with i-band observations optimized for the SBF measurements. Here, we run our pipeline on a selection of 20 galaxies from the NGVS with different total magnitudes and colors, to verify the results of our procedure. The data are the same used in Cantiello et al. (2024), and were obtained from the NGVS archives, now public; the images are normalized and calibrated such that the zero-point is m_zp = 30 mag in all bands.

3. SBF measurements with FAST-SBF

To measure the SBF magnitudes, we followed an approach similar to the well established procedures already developed in previous works (e.g., Blakeslee et al. 2009, 2021; Cantiello et al. 2018, 2024; Jensen et al. 2021; Mei et al. 2005), automated as much as possible to efficiently handle large datasets. We made an effort to use publicly available Python packages for this purpose. Below, we provide an overview of the technique, adding some practical notes, and detail the changes and advancements introduced. Figure 2 shows the flow chart of the adopted procedure. One of the main purposes of the present work is to show the progress in the automation and Python porting of the code, named FLEXIBLE AUTOMATED SELF-CONTAINED TOOL FOR SBF (FAST-SBF), and establish its consistency with previous literature results. The pipeline is currently under active development, and certain components will be modified and improved prior to its public release¹. The procedure was developed to work for any optical to near-IR passband. The band-dependent parameters (e.g., GCs and background galaxies luminosity function fitting, extinction correction, instrumental zero-point, etc.) are user-defined, and are provided in a configuration file.

Fig. 2. Flow chart of the main process of the FAST-SBF pipeline.

3.1. Sky background

The first step of our procedure entails accurately determining the sky background level to be subtracted from the image of the galaxy cutout². Any uncertainty in the sky level directly impacts the SBF signal and the galaxy color measurements, as well as the associated errors. Typically, the sky uncertainty contributes ≤0.02 mag of the total error on the SBF, corresponding to less than 1% on the distance (Cantiello & Blakeslee 2023). The survey data used here have already been sky-subtracted in the image-processing stage. However, a residual sky value might be present and needs to be estimated and subtracted.

The tests and results presented here rely on wide-field imaging data; therefore, the sky can generally be determined far from the galaxy center. We first estimate the median flux in circular annuli concentric to the galaxy, extending outward. The sky background value is then determined as the sigma-clipped median value of the radial flux profile, evaluated in a region away from the galaxy core, where the radial flux profile flattens out. The uncertainty on the background, σ_sky, is derived from the median absolute deviation (MAD) estimated within the same region (for a normal distribution, σ = 1.483 × MAD, where σ is the standard deviation). Figure 3 shows an example of the approach adopted for NGC 5831, one of the five galaxies with HSC-SSP data and in the Ton01 sample. This method may be less reliable when the galaxy covers most of the field. The present version of the code is tailored for wide-field surveys such as the HSC-SSP and NGVS, and hence estimating the sky does not generally represent a problem. Future versions of the code, already in testing stage, will address the issues arising from observations where the galaxy covers entirely the detector (such as massive elliptical galaxies in the local Universe imaged with the Hubble or James Webb space telescopes).

Fig. 3. Example of the procedure to estimate the i-band background of NGC 5831. The solid black line shows the radial profile of the flux in the input image, azimuthally averaged in circular bins of 50 pixels. The shaded gray area is the region where the background is estimated. The dashed red line represents the background estimated with the method described in the text.

3.2. Galaxy modeling and subtraction

The SBF signal is defined as the ratio between the intrinsic variance of the stellar light distribution and the mean surface brightness within a given region of the galaxy. Thus, a smooth isophotal model of the galaxy surface brightness after sky subtraction is needed; the resulting model frame corresponds to the first moment of the light distribution. In order to perform the modeling of the galaxy, we first derived the photometry of compact and extended sources using SExtractor (Bertin & Arnouts 1996). The photometric catalog was used to mask the brightest sources in the frame, adopting the appropriate source shape³, magnitude cuts (in the same magnitude scale of the images), and masking radius scaled to source magnitude.

The background subtracted masked frame was then used to derive a 2D model of the galaxy surface brightness, with the ELLIPSE routine (Jedrzejewski 1987) from the Photutils⁴ package. This routine offers reliable results with minimal user intervention, but it is computationally time-consuming for extended galaxies (> 10⁵ pixels). Moreover, a relevant limitation occurs when modeling isophotes heavily contaminated by globular clusters, bright companions, background sources, dust patches, or bright stars. Most of these sources are automatically masked using the compact and extended source catalogs. However, some objects (such as extended dust patches) or features (e.g., bright stellar halos or diffraction spikes) may evade the preliminary masking stage. To address this, we incorporated an optional manual masking process, based on an interactive graphical interface.

Usually, the model-subtracted image shows residual patterns from the model subtraction stage, especially in cases where central dust lanes, shells, bars or spiral sub-structures are present in the target galaxy. As will be discussed later (Sect. 3.5), the presence of large-scale residual patterns is not a limitation for the procedure. Nevertheless, the large-scale residuals still present in the frame after subtracting the galaxy model are removed using a background map obtained with SExtractor⁵. In contrast, small-scale residual structures, such as localized rings, shells, or dust clumps, can introduce biases in the measurement and are therefore masked during the analysis (see Sect. 3.5). Hereafter, we refer to the image with the sky, galaxy model, and large-scale residuals subtracted as the residual frame. The top panels of Fig. 4 show the procedure described in this section for the case of NGC 5813. The images show the i-band galaxy cutout, the residual, and residual masked image after the automatic masking of the detected sources (left to right). The right panel also indicates the annulus used for the SBF measurement.

Fig. 4. SBF analysis images and plots for NGC 5813. Top row: i-band image, residual, and residual masked image (left to right). The black circles in the third panel show the inner and outer radii of the annulus adopted for the SBF measurements. Middle row: Example of the PSF selection process. In the left panel are shown the stars pre-selected based on photometric criteria. Stars highlighted with red borders represent a subsample of those subsequently chosen for the fitting. In the right panel, the radial profile of the full sample of PSF selected is shown. The red line shows the median profile, while thin gray lines are azimuthal averaged profiles of single PSFs. The PSFs are normalized to have a total energy of one. Bottom left panel: Fit luminosity function of external sources. The green circles represent the binned observational data. The solid green curve is the best fit to the data. The two components of the luminosity function, the background galaxies and GCLF, are shown with dotted blue and dashed red curves, respectively. Bottom middle panel: Azimuthal average of the residual image power spectrum using one of the good PSFs. The gray dots represent the observational data. The solid blue curve is the fit obtained according to the procedure described in the text. The dashed red lines represent the range of wavenumbers used for the fit. The dashed green curve shows the PSF power spectrum. The horizontal dotted green line represents the fit white noise value P1. Bottom right panel: Estimated P0 (gray dots) at different kstart; the vertical dotted black lines indicate the most stable region in which P0 is calculated, while the horizontal one is the final value of P0 for the specific PSF.

3.3. Point spread function selection

In the ideal case, without atmospheric or instrumental blurring and in the absence of external sources of fluctuation beyond the stellar counts, the SBF signal would appear indistinguishable from the instrumental white noise. In practice, the convolution with the point-spread function (PSF) introduces correlations between adjacent pixels. The SBF signal is measured in Fourier space, where the power spectra of the signals convolved in the physical space (PSF and stellar variance) are multiplied by each other (see Sect. 3.5). For this reason, a key ingredient for accurate SBF measurements is the image PSF (see Blakeslee et al. 1999; Moresco et al. 2022; Cantiello & Blakeslee 2023). Owing to the wide-field images available, in this work we characterize the PSF relying on local stars.

We compiled a preliminary library of PSF reference stars by extracting bright and isolated point sources from the regions around the target galaxies obtained by running SExtractor on the residual frame. For each field, we selected compact objects with magnitudes ranging from half a magnitude fainter than the saturation limit (m_i ∼ 18 mag for all the fields considered), down to two magnitudes fainter.At the distance of our sample and the given image resolution, globular clusters (GCs) appear as point sources and are therefore indistinguishable from stars. However, the magnitude range of the GCs host in our sample galaxies (see Sect. 3.4) is significantly fainter than the range adopted for the PSF selection. Therefore, when constructing the PSF sample, we do not expect any confusion between GCs and stars. The magnitude range is defined in the parameter configuration file; hence, it can be extended for fields lacking isolate, bright, and compact PSF candidates. To minimize contamination from the galaxy light, we selected stars located at a distance of r_PSF = 10R_e effective radii from the target. To identify isolated stars, only sources at least 10 FWHM from any other detected source were kept. The middle row in Fig. 4 shows an example of the PSF selection process around NGC 5813. Some of the selected stars exhibit a residual contamination from undetected faint sources or residual image artifacts. To refine the selection we adopted a two-steps strategy. First we performed an azimuthal average of the normalized flux for each candidate, then the profiles were combined using a sigma-clipping algorithm to build a reference PSF (solid red line in Fig. 4, middle right panel). The selection of best PSF candidates was made by comparing the azimuthal profile of each individual star to the reference PSF, through a chi-square minimization. Second, we further cleaned the PSF cutout by replacing all outlier pixels beyond a given distance from the PSF core (20 × FWHM for the data analyzed here) with a pixel value extracted randomly from the distribution along the same radial bin. Finally, the background was estimated around each individual PSF, and is subtracted from the PSF cutout.

In this work we selected an average of ten good PSFs for each galaxy. A subsample of the good PSFs identified for the case of NGC 5813 is shown in Fig. 4, highlighted with red borders (middle, left panel). The PSF selection and cleaning procedure can be tailored by changing the isolation criteria, the magnitude ranges, SExtractor stellarity index, etc. The entire procedure of PSF selection is completely automatic and computationally quite efficient (≈10 seconds to run on each galaxy in our reference sample). To obtain the amplitude of the fluctuation signal, the power spectrum analysis is performed independently as many times as the number of good PSF identified (see Sect. 3.5).

3.4. Residual power from globular clusters and background galaxies

In addition to the stellar counts fluctuations that generate the SBF signal, the total fluctuation signal is also affected by contamination from compact and extended sources, mainly GCs and background galaxies (e.g., Blakeslee & Tonry 1995; Tonry et al. 1997). The first step to reduce the contamination from the non-stellar fluctuations is to mask all detected GCs and background sources. We used the photometric catalog derived running SExtractor on the residual frame (Sect. 3.3) to mask all such contaminating sources. Sources fainter than the detection limit may still represent a source of contamination to the total fluctuations measured. The procedure we adopted with FAST-SBF to correct for the contamination from undetected sources is a python version of the well-established strategy adopted in previous works (e.g., Blakeslee & Tonry 1995; Cantiello et al. 2005; Jensen et al. 2015). This step entails modeling the GC luminosity function (GCLF) as a Gaussian distribution and the background galaxies luminosity function as a power law, using the photometric catalog derived in the previous step. The parameters adopted for fitting both components to the luminosity function are user-defined. In particular, we assumed a Gaussian width in the range σ_GCLF = 0.9 − 1.4 mag for the dwarf/bright galaxies (Villegas et al. 2010) and an i-band GCLF peak magnitude of M_i = −8.0 mag (Harris et al. 2002; Durrell et al. 2014). The i-band background galaxy luminosity function is fit assuming a power law with slope of γ = 0.34 (Bernstein & Jarvis 2002; Cantiello et al. 2005).

The residual power contribution, P_r, arising from faint undetected GCs and background galaxies, is then computed by integrating the composite LF beyond the photometric completeness limit (Blakeslee & Tonry 1995; Cantiello et al. 2007). This residual power was subtracted from the total fluctuation power, P₀ (Sect. 3.5), to isolate the stellar SBF signal: P_f = P₀ − P_r.

Previous studies have shown that in optical passbands the GCLF can be reliably constrained with observations reaching > 0.5 mag fainter than the peak of the GCLF (Cantiello & Blakeslee 2023). At such a depth, the residual contribution from unmasked, contaminant sources to the SBF signal is reduced to < 10%, and it becomes negligible if the completeness limit extends > 1 mag beyond the GCLF peak. For the bright galaxies in this study, we have an average contamination of P_r/P₀ ∼ 16%. GC contamination does not represent an issue for the dwarf galaxies inspected here, because their GC populations are typically negligible (e.g., Marleau et al. 2024, and references therein). Thus, for the dwarfs, we neglected the contribution from GCs to P_r, and only accounted for the contribution from background galaxies in the correction of the SBF signal.

A common ingredient in SBF measurements pipelines is the modified weight map used for the (second, in our case) SExtractor run. This is required to account for the high-amplitude fluctuations close to the central regions of the galaxy where the high surface brightness causes the SExtractor to detect the SBF star-count fluctuations as if they were GCs, which can then end up masked and bias P_r. To prevent false detections of the SBF signal as spurious sources, especially near the core of the brightest galaxies, the variance map used by SExtractor includes the galaxy model profile scaled by a constant factor, κ, to balance noise suppression and source detection (see e.g., Jordán et al. 2004). Through empirical calibration, we determined κ ∼ 0.02 for the NGVS sample and κ ∼ 0.002 for the HSC-SSP dataset. Future developments of FAST-SBF will include a procedure for determining κ.

The process to mask sources and estimate P_r is fully automated beyond some iteration required to determine κ. This step is quite efficient in terms of computational time and for the convergence of the fitting algorithms. The parameter settings in the detection algorithm and luminosity function fitting are strictly telescope- and band-dependent, and thus are flexible and user-defined.

3.5. Fluctuation magnitudes

The previous steps produced the residual masked frame, a list of reference PSF stars and a table containing the input for estimating the residual fluctuation due to undetected sources, P_r. The next step is to identify the region of the galaxy to be used for the measurement of the SBF signal.

For bright galaxies, the fluctuation signal is strongest in the central regions, but the fraction of useful pixels in this region can be relatively small (because of poor modeling of the steep brightness profile gradients, the presence of contaminating dust, core saturation, etc.) and the P_r correction is largest due to the central concentration of GCs. In the outer galaxy regions, and for fainter galaxies, the S/N ratio of the fluctuation signal is typically lower, making the measurement more vulnerable to background variations and increased noise levels. However, the region available is, typically, wider and the contamination from GCs and background galaxies is lower. There is therefore a region between the two extremes where the SBF signal is optimally measured.

In our procedure, the SBF signal was measured within a single annulus, individually defined for each galaxy. At present, the inner and outer radii of the measurement annulus are still not well automated, although they are relatively stable for galaxies with similar magnitudes and morphologies. The outer boundary of the annulus, r_out, was chosen such that the ratio between the rms of the galaxy flux and the sky+detector background rms is S/N(r_out)≥5, for bright galaxies⁶. The inner boundary, r_in, was set to exclude regions potentially affected by contamination, such as dust lanes or residual modeling artifacts. Additionally, the optional interactive graphic interface was run again at this stage to manually mask any remaining contamination left out in the previous steps.

The Fourier spatial power spectrum of the annular residual-masked frame was then derived using the Numpy (Harris et al. 2020) package Fast Fourier Transform⁷ and the resulting 2D image was azimuthally averaged to obtain the profile, P(k), as a function of the wavenumber, k. The fluctuation amplitude, P₀, was estimated by fitting the equation

$$\begin{array}{c}\hfill P(k)=P_0\times E(k)+P_1.\end{array}$$(1)

where E(k) is the (azimuthally averaged) power spectrum of the PSF (normalized to have a total energy of one) convolved with the mask (Tonry et al. 1990), and P₁ is the instrumental white noise component arising from photon shot- and readout-noise, independent of the wavenumber, k. The convolution was performed using the Numpy package convolve⁸, while for the fit of P(k) to E(k) we used the linear regression method with the Scipy (Virtanen et al. 2020) package curve_fit⁹. The parameters, P₀ and P₁, were then obtained as outputs of this fitting procedure. The power, P₀, was converted to apparent SBF magnitude, $\overline{m}$, using the relation

$$\begin{array}{c}\hfill \overline{m}=-2.5{log}_{10}(P_0-P_{\mathrm{r}})+m_{\mathrm{zp}},\end{array}$$(2)

where m_zp is the photometric zero point for the image. Finally, the band-dependent foreground extinction correction was applied, using A_i = 1.684 × E(B − V) (Schlafly & Finkbeiner 2011).

The fitting procedure for Eq. (1) can be replaced by a robust fitting algorithm based on the “random sample consensus” method (RANSAC, Fischler & Bolles 1981). In normal cases, the RANSAC routine produces results equivalent to the default fitting method, and is more time-consuming than the default method. However, this option is useful in cases where extreme outliers are present in the P(k) or E(k) radial profiles.

The power spectrum, P(k), for the residual image can be distorted by large-scale background patterns from the galaxy model fitting procedure, which often adds power at low wavenumbers. Similarly, the flat noise power spectrum (P₁) can be tilted when the pixel interpolation kernel used in combining images produces correlations between pixels. We therefore excluded both low and high wavenumbers when fitting for P₀ Eq. (1) (see Moresco et al. 2022, and references therein). The region of valid k numbers used for the fitting, k_ini − k_end, is highlighted with vertical red lines in Fig. 4 (mid-lower panel).

Here, we summarize our method of identifying the optimal k wavenumber interval to estimate P₀, which is basically the same as the method adopted by, for example, Cantiello et al. (2024). For a given E(k), corresponding to a specific PSF, a series of fits to Eq. (1) was performed by incrementally increasing the starting wavenumber, k_start, up to the user-defined upper limit, k_end. Each iteration yields a corresponding estimate of P₀. The optimal wavenumber range is defined as the interval where P₀ exhibits the highest stability (lower rms), as indicated by the vertical dotted lines in Fig. 4 (bottom row, middle panel). The lower right panel of the figure shows the fit P₀ value as a function of the starting wavenumber k_start used for the fit. A key advantage of this approach is its efficiency and automation. The procedure is then iteratively applied to each selected PSF, and the final P₀ value is determined as the 3σ-clipped median of the individual estimates derived for each PSF.

3.6. Color measurement

The SBF signal is by definition sensitive to the properties of the stellar populations within the galaxy. Integrated colors have been used as a proxy for these characteristics to calibrate the absolute fluctuation amplitude (Tonry et al. 1990). Thus, an accurate measurement of the color is essential for reliable determination of $\overline{M}$ and, consequently, of the distance (Tonry et al. 1997; Blakeslee et al. 2001; Jensen et al. 2015; Cantiello et al. 2018). Our procedure takes full advantage of the multiband HSC-SSP and NGVS survey observations to measure integrated colors in the same annuli used for the SBF, and propagating the same masks (obtained from both bands) optimized during the fluctuation analysis stage. Total fluxes in both bands were estimated using this mask to get the galaxy color. The MW foreground extinction correction was applied using A_g = 3.237 × E(B − V) and A_i = 1.684 × E(B − V) (Schlafly & Finkbeiner 2011).

To quantify the color uncertainty, we propagated the contribution from the mean sky level, the instrumental readout noise and the sky background error. For each band, the uncertainty on the total flux, σ_fTot², was computed as

$$\begin{array}{c}\hfill {\sigma }_{f_{\mathrm{Tot}}}^2=(f_{\mathrm{Tot}}+N_{\mathrm{p}}\times Sky)+{(N_{\mathrm{p}}\times Sk{y}_{\mathrm{err}})}^2+N_{\mathrm{p}}\times 2RO{N}^2,\end{array}$$(3)

where f_Tot is the band-dependent total galaxy flux in the annulus, N_p is the number of unmasked pixel in the annulus, RON is the readout noise of the instrument, Sky is the mean sky level of the image, and Sky_err is the uncertainty on the mean sky level. The error on the (g − i) color was then derived as

$$\begin{array}{c}\hfill {\sigma }_{(\mathrm{g}-\mathrm{i})}=2.5/log(10)\times \sqrt{\frac{{\sigma }_{{\mathrm{f}}_{\mathrm{Tot},\mathrm{g}}}^2}{f_{\mathrm{Tot},\mathrm{g}}^2}+\frac{{\sigma }_{{\mathrm{f}}_{\mathrm{T}}\mathrm{ot},\mathrm{i}}^2}{f_{\mathrm{Tot},\mathrm{i}}^2}}·\end{array}$$(4)

3.7. Uncertainties and quality flag

Each step described above has an associated error, which then propagates in the following stages of the procedure. We treated these uncertainties as independent and, to be as consistent and conservative as possible, we added them in quadrature. These include:

• Power spectrum fitting. The statistical uncertainty of the fitting procedure for Eq. (1) reflects the quality of the residual image power spectrum. Typically this is a negligible error term, of the order of 0.005 mag on $\mathrm{\Delta }\overline{m}$;

• PSF templates. As the measure of $\overline{m}$ is performed on different individual PSF stars (defined as “good PSFs”, Sect. 3.3), the statistical difference between the PSFs has to be taken into account. We used the σ_MAD (defined in Sect. 3.1) of the P₀ distribution for the selected PSF references to estimate the uncertainties. This amounts for $\mathrm{\Delta }\overline{m}\sim 0.02$ mag, on average;

• Residual power from undetected sources. As is usually done in the relevant literature, we adopted an empirical uncertainty of 25% of the value of P_r (Cantiello et al. 2007; Blakeslee et al. 2009; Jensen et al. 2021, 2025). In this work, this contribution is of the order of $\mathrm{\Delta }\overline{m}\sim 0.04$ mag;

• Background subtraction. This uncertainty impacts both the normalization factor in the SBF amplitude (Tonry & Schneider 1988) and the color used in the $\overline{M}$ calibration. To estimate the effect of background uncertainty, we offset the background by ±1σ_sky (see Sect. 3.1) and repeated the entire analysis with the new sky value. Here, the average contribution of the background uncertainty is on the order of $\mathrm{\Delta }\overline{m}\sim 0.02$ mag. The typical impact is larger for fainter galaxies as it impacts on average for $\mathrm{\Delta }\overline{m}\sim 0.003$ mag on bright galaxies, and $\mathrm{\Delta }\overline{m}\sim 0.03$ mag for faint ones;

• Extinction correction. We adopted a relative error of 10% on the extinction (e.g., Schlafly & Finkbeiner 2011).

In addition to propagated photometric errors, we associated a quality flag with each galaxy. To define the confidence of the SBF analysis, we assigned a quality classification code to each $\overline{m}$ presented here, similar to the approach adopted in Cantiello et al. (2024). The classification is assigned considering the quality of the isophotal model residuals, possible residual contamination effects of bright nearby objects, the degree of masking, the regularity of the power spectrum fitting stage, and the stability of the SBF signal measurement to small changes in the input parameters. We classify the sources as q1 for excellent quality, q2 for good quality, and q3 for poor quality. All bright galaxies are classified as q1, while the dwarf galaxies range between q2 and q3. In the present version of FAST-SBF, the combination of these factors (e.g., quality of the residuals) is not necessarily reflected in the global error budget. In the future we plan to automatize the procedure to define q_i using quantitative tests on the effects described above.

3.8. Distance determination

For SBF, to estimate the distance moduli of the measured galaxies, an estimate of the absolute SBF magnitude $\overline{M_{\mathrm{i}}}$ is needed. For this purpose we adopted the calibration equation from Cantiello et al. (2024) for the (g − i) color:

$$\begin{array}{c}\hfill \overline{M_{\mathrm{i}}}=a_0{x}^3+a_1{x}^2+a_{2}x+a_3,\end{array}$$(5)

where a₀ = 10.085 ± 6.832, a₁ = 5.20 ± 2.042, a₂ = 2.81 ± 0.192 and a₃ = −1.129 ± 0.022, and x = (g − i)−0.90 mag represents the difference between the measured integrated color and the reference color. This calibration was obtained in the CFHT magnitude system, which we transformed to HSC passbands using the equations derived by Kim & Lee (2021) from PARSEC isochrones (Bressan et al. 2012) and, similarly to these authors, assume that the i-band magnitudes of the two systems are almost identical, while for the g-band a correction is required. Inverting the equation given in Sect. 5 of Kim & Lee (2021), we derive

$$\begin{array}{c}\hfill {(g-i)}_{\mathrm{CFHT}}=0.956\times [{(g-i)}_{\mathrm{HSC}}-0.006]\end{array}$$(6)

with an rms = 0.002 mag from the PARSEC isochrones (see Sect. 5 of Kim & Lee 2021). Table A.2 reports the results of the SBF analysis for the selected galaxies.

4. Analysis

The primary object of this study is to introduce a new, self-contained SBF pipeline, developed to be flexible and automated, and to study the prospects for measuring SBF in the LSST by using the HSC precursor data. In this section we first assess the robustness of the new pipeline on the NGVS sample, then estimate the distance to the galaxies observed by the HSC-SSP and compare the results to literature and stellar population models.

4.1. Comparison with independent NGVS measurements

We compared our measurements for the SBF signal, color, and distance modulus with the results recently obtained by Cantiello et al. (2024). The comparison was carried out for 20 galaxies in the Virgo cluster observed in the NGVS observational campaign. The galaxies were selected to cover a wide range of magnitude and color. Furthermore, since our goal is to demonstrate the consistency between the results of the two different analysis procedures, we limited the comparison to galaxies labeled with code q1 in Cantiello et al. (2024). This was done to avoid discrepancies arising from the specific galaxy or field properties, rather than from differences in the algorithms themselves. Figure 5 shows the difference in $\overline{m}$ and distance modulus as a function of color (left panels) and luminosity (right panels) between the two sets of measurements. The galaxies selected for this comparison span a range of color 0.7 ≤ (g − i)≤1.2 and magnitude 9 ≤ m_i ≤ 15 mag. The SBF measurements are consistent within the errors. Even in the few cases in which the discrepancy appears larger (e.g., VCC 0538), we find excellent agreement with previous measurements when the comparison is expressed in terms of distance modulus. The difference can be attributed to minor discrepancies in the masking and the annulus size used for the measurements. These factors introduce small offsets in the $\overline{m}$ estimates but, since it also affects the color measurement, $\overline{M}$ is affected too. Hence, the final distance modulus estimates remain well within 1σ consistency. In conclusion, the tests performed on the NGVS sample show an excellent agreement between the previous results (based on a pipeline which required heavy human intervention at all stages) and the estimates derived with the new FAST-SBF pipeline, where some modules have been successfully automated, with a net gain in speed and flexibility.

Fig. 5. Comparison with the measurements based on NGVS data. The left panels show the SBF magnitude (upper panel) and distance modulus (lower panel) differences between our measurements and Cantiello et al. (2024), as a function of the (g − i) color. The right panels show the same comparison but as a function of the galaxy total i-band magnitude.

4.2. Measurements on the HSC-SSP sample

As is discussed in Sect. 1, the SBF method is one of the few standard candles that allows measuring distances to dwarf galaxies in the local Universe when single stars are not resolved. From the perspective of the forthcoming start of the LSST observations, it is realistic to expect that distances to thousands of dwarfs (or more) will be within reach. HSC-SSP data serve as a valuable test bed for such expectations, since the data reduction pipelines developed for LSST are built upon the same software framework as the HSC-SSP science pipelines. Here, we discuss the results obtained on the sample of five bright galaxies with contemporary HSC-SSP g and i band data and available in the Ton01 catalog (Sect. 2). In addition to these reference targets, we extended our measurements to a collection of visually identified dwarf galaxies with brightness profiles appropriate for SBF analysis, brighter than m_B ∼ 18 mag and within two degrees of the bright galaxies. In all the figures that follow, the colors and magnitudes reported in the table have been extinction corrected.

In Fig. 6 we show separately the measured SBF magnitudes as a function of galaxy (g − i) color in the three fields. To use the calibration equation given in Eq. (5), we converted the HSC (g − i) colors to the NGVS system using Eq. (6). In the figure, the color bar represents the velocity of individual galaxies relative to the CMB frame. The dashed black line indicates the calibration relation from Cantiello et al. (2024), shifted according to the mean distance modulus of the group, which we derived from our distances using only the confirmed group members (Kourkchi & Tully 2017). Below we discuss the three fields separately.

Fig. 6. SBF magnitude as a function of color. Left, middle and right panels show the three different fields analyzed in this work, with the same order as in Fig. 1). The color bar represents the velocity of the sources, vCMB. The various symbol shapes represent the assigned quality flag for the SBF measurement: circles for galaxies classified as q1, triangles for q2, and “X” for q3. Different sizes distinguish bright from dwarf galaxies. The dashed black lines represent the calibration formula from Cantiello et al. (2024) shifted as discussed in the text. The shaded gray area is the mean cosmic scatter from the relation σcos, as estimated in Cantiello et al. (2024).

• Virgo III Cloud – Three of the bright galaxies in our sample, NGC 5813, NGC 5831 and NGC 5839, are members of NGC 5846 group in the Virgo III cloud (Aaronson et al. 1982). Kourkchi & Tully (2017) report that the three galaxies belong to a larger group comprising 84 members. We identified a sample of 20 dwarf galaxies which all appear, together with the bright galaxies, to lie on the expected $\overline{m}$-color relation for galaxies at the same adopted group distance. This result, combined with the recession velocities, suggest that the measured galaxies are gravitationally bound to the same group. With the exception of PGC 1199471, all the identified galaxies are reported as group members in Kourkchi & Tully (2017). Considering only the three bright galaxies, we estimate a group distance modulus of m − M = 32.09 ± 0.12 corresponding to a distance of D = 26.2 ± 1.5 Mpc. We ignore the dwarfs to avoid any bias arising from an incomplete selection of dwarf galaxies. When the sample of dwarfs is included, we find m − M = 32.08 ± 0.21 (or 26.1 ± 2.6 Mpc).

We observe the majority of dwarf galaxies in this field (left panel in the Fig. 6) to align with the calibration relation, with an average scatter of 0.21 mag. The observed rms with respect to the mean relation corresponds to an approximate depth of 0.15 mag, namely ∼1.8 Mpc depth of the group, after subtracting the expected cosmic scatter (σ_cos = 0.09, Cantiello et al. 2024) and the median $\overline{m}$ uncertainty. We find all sources, with the exception of PGC 1230503, MATLAS 1944, MATLAS 1945, and MATLAS 2005, to be located in projection around 1.5 Mpc from the nearest reference elliptical galaxy. The remaining four dwarf galaxies, instead, are found at the projected distances of ≈3.5 Mpc from the main group. The information derived is of paramount interest for membership confirmation for such loose, low-surface-brightness galaxies where spectroscopic membership typically requires prohibitive integration times, even with 8 m-class telescopes.

To test the usability of the SBF signal when the signal is low but still measurable for dwarf low-surface-brightness galaxies, we compared the power spectrum of one of the average dwarfs in our sample (namely PGC 1165764) to the power spectrum of an empty field and a background galaxy. Figure 7 shows, from left to right, the power spectra derived over an empty region, for a background galaxy, and of PGC 1165764, analyzed using the same configuration parameters. The background galaxy is PGC 1174699, at v_CMB = 12 522 km/s, and M_g ≈ −20 mag, with a distance from the Hubble law of ≈170 Mpc, assuming H₀ ∼ 73.5 km/s/Mpc (Jensen et al. 2025). The fluctuation signal measured for the dwarf, as observed by the relative amplitude of the PSF power spectrum (dash-dotted green lines in the panels), compared to the background galaxy and the empty field, is evidence of a stellar count fluctuation signal being detected, and therefore of its relative closeness. For a quantitative comparison, the ratio between the fluctuation signal and the white noise component is P₀/P₁ ≈ 20 for the background galaxy PGC 1174699, in comparison with P₀/P₁ ≈ 300 for the dwarf galaxy PGC 1165764. For PGC 1174699, given its large distance, the small signal measured is attributed entirely to its undetected GC population (Blakeslee & Tonry 1995). Moreover, the test on the empty field indicates that the amplitude measured for the dwarf is not a bias or an artificial effect from, for example, data processing.

Fig. 7. Comparison of power spectra measured in an empty region (left panel), for a background galaxy (PGC 1174699, middle panel), and one of the dwarf galaxies in our sample (PGC 1165764, right panel). The colors, points and line styles are the same as in Fig. 4 (middle low panel).

Before comparing our distances for the single galaxies with previous estimates, we have to consider that the Ton01 measurements were performed on a collection of data acquired in the 1990’s with 2.4 to 4 m aperture telescopes and using smaller format CCDs (Tonry et al. 1997). As already discussed in Cantiello et al. (2018), in a few cases the distances from the first successful ground-based SBF survey could be significantly off. The lower resolution could, for example, result in a poor characterization of the GCLF, possibly leading to issues in the estimates of P_r.

Our distance modulus for NGC 5813 is ≈0.6 mag smaller (∼30% closer in distance) than Ton01. The two estimates differ by ≈2.8σ. The Ton01 quality flags for the galaxy (PD and Q parameters) are slightly below average. We find m − M = 31.92 ± 0.12 mag (namely D = 24.2 ± 1.3 Mpc), which places NGC 5813 closer to the mean group distance (26.7 ± 1.1 Mpc) with respect to Ton01. Figure 4 shows the i-band image of the galaxy-subtracted frame, with the annulus used for the SBF measurement. The galaxy has a large population of GCs, concentrated in the central regions, making the constraint of P_r crucial. Moreover, we also observe the presence of dust lanes around the galaxy core, which we masked before performing the SBF measurement. The poorer resolution of the older Ton01 data (see also Fig. 11 of Cantiello et al. 2018) could at least in part explain the observed difference. In the case of the HSC-SSP imaging, the photometric completeness drops around ≈26 mag, which is ≈1.5 mag fainter than the i-band turnover of the expected peak of the GCLF (see Fig. 4), and is ideal for precise SBF measurement (Cantiello et al. 2007; Blakeslee et al. 2009).

NGC 5839 has several redshift-independent distance estimates in the literature that we can compare with. Our distance modulus is ≈0.4 larger than Ton01. The ground-based distance reported by Ton01 is superseded by the more recent estimate by Jensen et al. (2021), based on NIR measurements from HST/WFC3 data. They report (m − M) = 32.32 ± 0.08 mag, while we find (m − M) = 32.17 ± 0.19, with a difference of ≈0.7σ. Another available measurement comes from the Type Ia supernova 2015bp as part of the Pantheon+SH0ES program (Riess et al. 2022). They report (m − M) = 32.38 ± 0.32, within 0.55σ of our measurement. We note that the central region of NGC 5839 is contaminated by the presence of an additional source, most likely a background galaxy, making it difficult to accurately model. Moreover, Jensen et al. (2021) reported a central spiral component observed with HST F110W band, which is barely resolved in HSC i-band imaging. To mitigate these issues, we masked the central region when performing the SBF measurement. Unlike NGC 5813, the sensitivity limit of the field of NGC 5839 drops around ≈24.5 − 25 mag. Because of this, we are able to observe ≈0.5 − 1 mag fainter than the turnover of the GCLF, placing the galaxy at the limiting distance to have a reliable measurement of the SBF signal.

For NGC 5831, Ton01 reported a distance of (m − M) = 32.14 ± 0.17 mag, in perfect match with our measurement of (m − M) = 32.17 ± 0.13 mag. The galaxy is average in terms of Ton01 quality flags, and in the HSC-SSP image there are no visible artifacts or anomalies.

• Virgo II Cloud – In the middle panel of Fig. 6 we show the results for the field of NGC 4753. As reported in Kourkchi & Tully (2017), NGC 4753 is the brightest in a group of eight galaxies, five of which are observed in the HSC-SSP PDR2, and analyzed in this work. Moreover, within 2 deg of NGC 4753, we identified another bright galaxy, NGC 4690, and three fainter galaxies (namely MATLAS 1477, MATLAS 1491 and PGC 135809) which presented characteristics useful for the SBF analysis. We find that all the galaxies reported as group members of NGC 4753 in Kourkchi & Tully (2017) follow very well the NGVS $\overline{M}$-(g − i) calibration (Eq. (1)) shifted to the group distance, which we estimated as m − M = 31.77 ± 0.17 mag (D = 22.6 ± 1.9 Mpc) using all the confirmed members. The group had no previous mean distance measured in literature. As the majority of the galaxies analyzed here are dwarfs with relatively high SBF uncertainties, and given the limited number of members, we lack the precision necessary to constrain the depth of the group. PGC 135809 shows slightly brighter $\overline{m}$ compared with the rest of the group, consistent with it being in the foreground of the group, a finding which is possibly coherent with its lower recession velocity.

For NGC 4753, Ton01 reported a distance of (m − M) = 31.90 ± 0.19 mag, in excellent agreement with our estimate (m − M) = 31.92 ± 0.13 mag (D = 24.3 ± 1.5 Mpc). The galaxy presents diffuse dust lanes extending from the center of the galaxy to the outskirts, which we masked based on the HSC-SSP g-band image.

In the same field, NGC 4690, MATLAS 1477 and MATLAS 1491 exhibit much fainter fluctuations, consistently with their velocity ∼1500 km/s higher than NGC 4753. The combined information about recession velocities and newly measured distances suggest membership to a single group in the background of NGC 4753. Indeed, NGC 4690 and MATLAS 1491 are reported as part of the same group by Kourkchi & Tully (2017). For NGC 4690, there are no previous distance measurements. We derive (m − M) = 32.62 ± 0.11 mag, which places the galaxy at D = 33.4 ± 1.7 Mpc.

• Leo II Cloud – In the last panel we show the results for the field of IC 0745, reported as member of the Leo II cloud (Castignani et al. 2022). The galaxy is reported as isolated in Kourkchi & Tully (2017), and indeed we visually found only two dwarf galaxies within 2 deg of the bright reference. According to our distance estimates, none of the observed dwarfs appears to be associated with IC 0745.

For IC 0745, Ton01 reported a quality flag largely below average. The galaxy shows strong contamination from an intermediate mass irregular galaxy at small galactocentric distance, and morphological irregularities that suggest the presence of dust patches. Comparing the i- and g-band images, we masked all regions potentially affected by dust contamination and restricted the analysis to the outskirts of the galaxy to avoid poorly modeled regions. Ton01 estimated the SBF signal in an inner annulus more centrally located than ours, covering a region possibly affected by dust extinction. They reported a distance modulus of (m − M) = 31.38 ± 0.30 mag, in comparison with our analysis which yields (m − M) = 31.63 ± 0.11 mag (D = 21.2 ± 1.1 Mpc). The two estimates agree within ≈0.8σ, also owning to the relative large errors of the older measurements. For the additional bright galaxy found in the field, IC 0753, we are only complete down to ∼0.5 mag brighter than the turnover of the GCLF rendering the constraint on the luminosity function unreliable. Furthermore, at its distance of ∼100 Mpc, it is similar to the background galaxy discussed previously (PGC 1174699), where the measured $\overline{m}$ is expected to be dominated by the contribution from the undetected GC population of the galaxy, rather than the unresolved stellar component. For these reasons, we chose to exclude IC 0753 from the final analysis.

4.3. Comparison with stellar population models

In this section we compare our measured SBF amplitudes with predictions from the Teramo-SPoT stellar population (SSP) models (Cantiello et al. 2005; Raimondo et al. 2005; Raimondo 2009). The SPot-SBF models are based on a stellar population synthesis code originally described by Brocato et al. (1999). The code is designed to accurately reproduce both the observed distribution of stars in color–magnitude diagrams and the integrated properties of star clusters and galaxies, including colors, spectral indices, and energy distributions. The models adopted are an updated version of the one showed in Cantiello et al. (2024), covering a wider [Fe/H] range to include the intermediate-mass and dwarf galaxies regime.

A comparison between empirical data and SSP model predictions is shown in Fig. 8. To compare the measured fluctuation magnitudes with the theoretical predictions, we shifted the apparent $\overline{m}$, adopting the mean distance modulus of each group. Most of the observed (g − i) galactic colors agree with the sequences predicted by SSP models across all metallicity ranges, for ages between 3 and 14 billion years, with two outliers at the bluest end. Note that SSP models can provide insights into the dominant stellar populations in the galaxy or in the galactic regions where SBFs were measured, while composite population models would be more appropriate for a detailed study of the stellar content of the galaxies (e.g., Vazdekis et al. 2020).

Fig. 8. Comparison of our measurements with models. The diamonds show the SPoT SSP models for different metallicities, given by the color bar. The symbols and sizes of the galaxies are the same as in Fig. 6. The dashed black line represents the calibration equation from Cantiello et al. (2024) shifted as discussed in the text. The shaded gray area is the mean cosmic scatter from the empirical calibration.

All the galaxies of the Virgo III cloud are very well aligned with predictions. We find four reference galaxies (NGC 5813, 5831, 5839, and 4753) and NGC 4690 to be consistent with the set of models around solar metallicity, and stellar populations with ages on the order of 3–14 Gyr. As was expected, the three elliptical galaxies in the Virgo III cloud (NGC 5813, NGC 5839, NGC 5831) cover the region for old stellar population age. NGC 4753 and NGC 4690 instead, are aligned with ages of around 3–6 Gyr, which is not unexpected considering their relatively irregular morphology. IC 0745 shows sign of star-formation in the very center, which can explain its blue color and position in the region of models having sub-solar metallicity and young ages. On the other hand, most of the dwarf galaxies in the sample align with SSP models having [Fe/H] = –0.10 dex and lower, and with ages ranging from 3 to 8 Gyr, in agreement with previous studies suggesting that such systems generally host relatively old (5–10 Gyr) and subsolar-metallicity stellar populations (Kirby et al. 2013; Weisz et al. 2014). In conclusion, the agreement between measures and predictions for the appropriate range of age and metallicity provides further support to the reliability of the results presented.

5. Summary

In this work, we have introduced FAST-SBF, a new, flexible Python-based pipeline designed for SBF measurements. The code has been developed for upcoming next generation of large astronomical surveys (particularly LSST, Euclid, and the Roman Space Telescope), with the aim of building a tool that is fully automated and able to handle large datasets. In this work, we tested the procedure on ground-based imaging data from the HSC-SSP survey, an LSST precursor, taking as references five galaxies already studied in the literature. Starting from these, we expanded our analysis to include both bright and dwarf galaxies in the same fields. In general, we obtained accurate measurements for the distances of the galaxies already known in the literature, reporting results typically in good agreement with the previous measures, confirming the robustness of the new pipeline presented here. Moreover, we extended the analysis to dwarf galaxies lacking prior distance estimates, and used their consistency with the empirical SBF–color relation to assess group membership.

Further validation was carried out using measures from the more recent NGVS program, where we compared our results with the estimates from Cantiello et al. (2024) for 20 galaxies. Again, we found excellent consistency across SBF magnitudes, colors, and derived distances. The measurements also show good agreement with predictions from stellar population synthesis models, particularly in the age and metallicity range expected for early-type galaxies and dwarfs in the local Universe.

Overall, the results demonstrate the reliability of the new procedure and its potential for enabling precise, automated distance measurements across a wide variety of galaxy types. The procedure has also allowed us to assess the group membership of low-luminosity dwarf galaxies.

The current version of FAST-SBF is more automated and portable than previous implementations of SBF measurement pipelines. However, the code is not yet ready for public release. At this stage, obtaining reliable results requires specific training to understand the caveats and to correctly handle those modules that demand particular attention (e.g., contamination from residual sources and proper power-spectrum fitting). While the code can produce fluctuation magnitudes for any target, the scientific validity of these measurements depends critically on both a thorough understanding of the data and the correct execution of the analysis steps. Researchers interested in using FAST-SBF are encouraged to contact the second author (at michele.cantiello@inaf.it) to discuss access to the code and receive the necessary training for its appropriate use.

Acknowledgments

GR acknowledges support from the INAF “Astrofisica Fondamentale” Mini-grant 2024 n.17/RSN1 (PI G. Riccio). MC acknowledges support from ASI–INAF grant no. 2024-10-HH.0 (WP8420), the ESO Scientific Visitor Programme, and INAF GO-grant no. 12/2024 (P.I. M. Cantiello). RH and MC acknowledge support from the project “INAF-EDGE” (Large Grant 12-2022, P.I. L. Hunt). NH acknowledges support from the Polish National Science Centre grant 2023/50/A/ST9/00579. MM acknowledges support from the ESO Studentship program. This research has made use of the NASA/IPAC Extra- galactic Database (NED (https://ned.ipac.caltech.edu/)), which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology. Software: Numpy (Harris et al. 2020), Scipy (Virtanen et al. 2020), Source Extractor (Bertin & Arnouts 1996). We are grateful to Adriano Pietrinferni and Santi Cassisi for providing us tables with key points of their evolutionary tracks published in Hidalgo et al. (2018), and Pietrinferni et al. (2021).

References

Appendix A: Data and SBF measurement tables

Here we report the properties of the HSC-SSP galaxies (Table A.1) and the results of the SBF analysis (Table A.2).

All Tables

All Figures

	Fig. 1. Spatial distribution of the galaxies in the three fields we studied. The color represents their velocity in respect to the CMB frame. Small circles represent dwarf galaxies, while bright galaxies are represented with bigger size and are labeled. The box size is 1.3 × 1 deg and is the same for all the panels.
In the text

	Fig. 2. Flow chart of the main process of the FAST-SBF pipeline.
In the text

	Fig. 3. Example of the procedure to estimate the i-band background of NGC 5831. The solid black line shows the radial profile of the flux in the input image, azimuthally averaged in circular bins of 50 pixels. The shaded gray area is the region where the background is estimated. The dashed red line represents the background estimated with the method described in the text.
In the text

Fig. 4. SBF analysis images and plots for NGC 5813. Top row: i-band image, residual, and residual masked image (left to right). The black circles in the third panel show the inner and outer radii of the annulus adopted for the SBF measurements. Middle row: Example of the PSF selection process. In the left panel are shown the stars pre-selected based on photometric criteria. Stars highlighted with red borders represent a subsample of those subsequently chosen for the fitting. In the right panel, the radial profile of the full sample of PSF selected is shown. The red line shows the median profile, while thin gray lines are azimuthal averaged profiles of single PSFs. The PSFs are normalized to have a total energy of one. Bottom left panel: Fit luminosity function of external sources. The green circles represent the binned observational data. The solid green curve is the best fit to the data. The two components of the luminosity function, the background galaxies and GCLF, are shown with dotted blue and dashed red curves, respectively. Bottom middle panel: Azimuthal average of the residual image power spectrum using one of the good PSFs. The gray dots represent the observational data. The solid blue curve is the fit obtained according to the procedure described in the text. The dashed red lines represent the range of wavenumbers used for the fit. The dashed green curve shows the PSF power spectrum. The horizontal dotted green line represents the fit white noise value P1. Bottom right panel: Estimated P0 (gray dots) at different kstart; the vertical dotted black lines indicate the most stable region in which P0 is calculated, while the horizontal one is the final value of P0 for the specific PSF.

In the text

	Fig. 5. Comparison with the measurements based on NGVS data. The left panels show the SBF magnitude (upper panel) and distance modulus (lower panel) differences between our measurements and Cantiello et al. (2024), as a function of the (g − i) color. The right panels show the same comparison but as a function of the galaxy total i-band magnitude.
In the text

Fig. 6. SBF magnitude as a function of color. Left, middle and right panels show the three different fields analyzed in this work, with the same order as in Fig. 1). The color bar represents the velocity of the sources, vCMB. The various symbol shapes represent the assigned quality flag for the SBF measurement: circles for galaxies classified as q1, triangles for q2, and “X” for q3. Different sizes distinguish bright from dwarf galaxies. The dashed black lines represent the calibration formula from Cantiello et al. (2024) shifted as discussed in the text. The shaded gray area is the mean cosmic scatter from the relation σcos, as estimated in Cantiello et al. (2024).

In the text

	Fig. 7. Comparison of power spectra measured in an empty region (left panel), for a background galaxy (PGC 1174699, middle panel), and one of the dwarf galaxies in our sample (PGC 1165764, right panel). The colors, points and line styles are the same as in Fig. 4 (middle low panel).
In the text

Fig. 8. Comparison of our measurements with models. The diamonds show the SPoT SSP models for different metallicities, given by the color bar. The symbols and sizes of the galaxies are the same as in Fig. 6. The dashed black line represents the calibration equation from Cantiello et al. (2024) shifted as discussed in the text. The shaded gray area is the mean cosmic scatter from the empirical calibration.

In the text