Testing the Surface Brightness Fluctuation Method on Dwarf Galaxies in the COSMOS Field

Dwarf galaxies are important tracers of small-scale cosmological structure, yet much of our knowledge about these systems comes from the limited sample of dwarf galaxies within the Local Group. To make a comprehensive inventory of dwarf populations in the local Universe, we require effective methods for deriving distance estimates for large numbers of faint, low surface brightness objects. Here we test the surface brightness fluctuation (SBF) method, traditionally applied to brighter early-type galaxies, on a sample of 20 nearby dwarf galaxies detected in the COSMOS field. These objects are partially resolved in HST ACS images, and have confirmed redshift distances in the range 17-130 Mpc. We discuss the many model choices required in applying the SBF method, and explore how these affect the final distance estimates. Amongst other variations on the method, when applying the SBF method, we alter the standard equation to include a term accounting for the power spectrum of the background, greatly improving our results. For the most robust modelling choices, we find a roughly Gaussian SBF signal that correlates linearly with distance out to distances of 50-100 Mpc, but with only a fraction of the power expected. At larger distances, there is excess power relative to that predicted, probably from undetected point sources. Overall, obtaining accurate SBF distances to faint, irregular galaxies remains challenging, but may yet prove possible with the inclusion of more information about galaxy properties and point source populations, and the use of more advanced techniques.


INTRODUCTION
The Lambda-Cold Dark Matter (ΛCDM) model underpins modern cosmology, and has proven a great success in predicting the observed abundance and distribution of galaxies on large scales (Mo et al. 2010).On small scales, galaxy formation is strongly modulated by baryonic effects, which clearly suppress the abundance of dwarf galaxies through internal and external feedback processes, although the exact details are unclear and problematic (Bullock & Boylan-Kolchin 2017).These feedback processes are multiple, and will often depend strongly on local environment, so large statistical samples of dwarfs from a range of different environments are needed to tease them apart.
Unfortunately, a significant part of our observational understanding of dwarf galaxies is limited to the Local Group and its immediate surroundings, where dwarf galaxies can be detected to faint magnitudes and low surface brightness limits, using a variety of techniques (McConnachie 2012).Within the local volume out to 10 Mpc, ongoing efforts are building up extensive samples of satellites (see e.g.Carlsten et al. 2022;Nashimoto et al. 2022, for full lists of references), but most are in systems similar to the Local Group in mass.Our census of dwarfs in more and in less massive systems is much less complete, despite the important information these environments provide about the underlying efficiency of galaxy formation (e.g.Sales et al. 2013;Grossauer et al. 2015;Müller & Jerjen 2020;Garling et al. 2021;Roberts et al. 2021, Xi & Taylor in prep.).Frustratingly, many of these systems are easily detectable in current wide-area ground based surveys, or will be detected in forthcoming spacebased surveys, and their average abundance is also increasingly well constrained (Wang & White 2012;Speller & Taylor 2014;Nierenberg et al. 2016;Tanaka et al. 2018;Xi et al. 2018;Roberts et al. 2021;Wang et al. 2021;Wu et al. 2022).The challenge is to identify exactly which of the huge numbers of objects detected in these surveys are in fact local, without resorting to complete spectroscopic surveys (e.g.Geha et al. 2017), which are both expensive and challenging, particularly for early-type, low surface brightness dwarfs.This has prompted renewed interest in extending distance estimates based on imaging to fainter populations (Xi et al. 2018;Carlsten et al. 2019;Polzin et al. 2021;Kim & Lee 2021;Greco et al. 2021;Carlsten et al. 2022).Xi et al. (2018) identified a set of galaxies in the COSMOS field that appear to be partially resolved into point sources and/or surface brightness fluctuations.The objects in the sample are confirmed to be local by their photometric or spectroscopic redshifts.Polzin et al. (2021) applied the surface brightness fluctuation (SBF - Tonry & Schneider 1988) method on the nearest object in the sample (COS-MOS ID 549719), and obtained a distance estimate in agreement with the redshift distance.This was an exciting but somewhat surprising result, as the SBF method is normally used on older, brighter ellipti-cal galaxies, whereas they also found the galaxy to have been recently star forming.A natural next step is to test whether this method can produce accurate distance estimates for the larger sample of partially resolved objects compiled by Xi et al. (2018).
In this paper, we apply the SBF method to a sample of 20 objects from the Xi et al. (2018) serendipitous catalogue, selected to have spectroscopic redshifts, and redshift-based distances of less than ∼130 Mpc.We compare the objects' SBF and proper distances, to explore the limits of the reliability of SBF distance estimates to galaxies of this type at these distances.
The outline of the paper is as follows.In Section 2, we describe our sample and the observational data used.In Section 3, we summarize each of the steps in the SBF method.In Section 4, we present a fiducial SBF distance estimate for each of the galaxies, following the SBF procedure outlined in Section 3. In Section 5, we test variations on the fiducial method, and determine their impact on the final results.In Section 6, we conclude by discussing the limitations of using the SBF method on dwarf galaxies at distances of 20-130 Mpc, and consider how the method might be improved in future work.

Galaxy Sample
The effective range of the SBF technique depends on the depth of the imaging used, but for bright galaxies imaged with HST, reaches 100-130 Mpc (Blakeslee et al. 2021;Moresco et al. 2022).For faint, low surface brightness galaxies, the range is unclear, but Xi et al. (2018) found a reasonable correspondence between general optical morphology (i.e.visual indications that systems were partially resolved) and distance for dwarf galaxies with  + < 20 mag, out to distances of ∼ 200 Mpc.To test the effective range of the technique for fainter, less regular galaxies, we considered the 34 objects in their serendipitous catalogue with spectroscopic redshifts that put them within ∼130 Mpc.
Of these, 14 galaxies were excluded from our final sample; 5 because they were clearly intrinsically bright and/or star-forming galaxies, 2 because they had an extremely irregular morphology that would make model fitting difficult (see 3.1.2),and 7 because they were too small, faint, or low surface brightness.The COSMOS2015 IDs, heliocentric redshift distances (d = / 0 , assuming  0 = 70 km s −1 Mpc −1 ), and coordinates of the remaining 20 objects are listed in Table 1, along with the cutout image sizes used in our subsequent analysis.

Imaging and Spectroscopy
Each of the galaxies in our final sample has multi-band photometry from the COSMOS 2015 catalogue (Laigle et al. 2016), as well as a spectroscopic redshift given in Xi et al. (2018).We used these catalogue redshifts in all cases except for 549719, where Polzin et al. (2021) measured a redshift of 1222 ±64 km s −1 , and 213165, whose catalogue redshift places the galaxy too close.For this galaxy, a corrected redshift was obtained from the NASA Extragalactic Database (NED) 1 .
The HST images used in the SBF analysis were HST/ACS F814Wband mosaics (Koekemoer et al. 2007;Massey et al. 2010)   from the NASA/IPAC IRSA COSMOS Cutouts Service, while Subaru  + , IA464, and IA484-band images (Taniguchi et al. 2015) from the same source were used to estimate the  − colour of the galaxies.

METHOD
If the light emitted by a galaxy were produced entirely by stellar point sources of fixed luminosity, we would expect the number  of these sources falling within a single pixel of the image to scale with the physical area subtended by the pixel, and thus to vary as the square of the geometric distance to the galaxy.Poisson fluctuations in this number from pixel to pixel, across a region of uniform brightness, should vary as  1/2 , that is as .For more realistic stellar populations with a range of intrinsic luminosity, the amplitude of fluctuations should depend on a weighted moment of the luminosity function, and thus on the age and metallicity of the population, and the resulting distribution would no longer be purely Poissonian (Cerviño & Luridiana 2006;Cerviño et al. 2008).Given an absolute calibration for a similar stellar population, however, one should still be able to infer the distance from the amplitude of pixel-to-pixel variations, relative to the mean surface brightness.Based on this idea, Tonry & Schneider (1988) first proposed the Surface Brightness Fluctuation (SBF) method for determining absolute distances to nearby galaxies.
In practice, the SBF method requires a number of distinct steps, including masking individual point sources and foreground/background objects, modelling and subtracting and/or dividing by the average light distribution across the galaxy, measuring the spatial power spectrum of the residual map, the background sky around the object, and the point spread function (PSF) of the image, fitting the overall power spectrum to these components, and adjusting the final result, cast as an 'SBF magnitude', for an age and metallicity-dependent zero-point.We discuss these steps, and the possible choices to make at each, below.

Masking and Modelling
An important first step in the SBF method is to mask any potential contaminants in the image that could give rise to extra fluctuations unrelated to discreteness noise in the underlying stellar population.Examples include foreground stars, background galaxies, globular clusters, and cosmetic artifacts in the image.Once the image is masked, we can create a smooth model of the galaxy's light distribution, to be used in subsequent steps.Masking and modelling are interconnected and happen simultaneously, so we will discuss them both in this section.
We use three methods for masking the image to identify bad regions and exclude them from subsequent analysis: manual masking, masking with automatic point source detection, and masking using the galaxy image pixel histogram.Each of the three methods builds on and incorporates the previous one(s).

Initial Manual Masking
We first masked the images manually, by inspecting them visually, noting which areas of the image contained bad pixels or contaminating objects, and masking them as simple circular or rectangular regions.The goal of this process was to remove the most obvious contaminating sources from each image, so that the subsequent steps were more efficient.

Sérsic Model
Given our initial masking, the next step was to model the smooth light distribution of the galaxy.To do this, we initially created elliptical Sérsic models of each galaxy, using imfit (Erwin 2015), with our manual mask indicating regions to ignore, and the absolute value of the most negative pixel in the image taken to be the value of the previously subtracted constant sky background.We used reasonable guesses for the initial parameters of the fit in each case, including  = 0.5,  (  ) = 0.01 counts/s, and   = 100 pixels, but found these had little impact on the final fitted parameters (see table A1), as did the constraints we put on parameters.To determine errors on the fit, we used imfit's bootstrap resampling command, with 200 iterations.Given the Sérsic Model fit, we restricted our subsequent SBF analysis to an elliptical region centred on the galaxy with a semi-major axis of 2  , and the same axis ratio as the Sérsic Model.

Bicubic Spline Model
We also created a second, non-parametric model for the smooth light distribution of the galaxy using bicubic spline interpolation (Jordán et al. 2004).First, we made a lower resolution version of the masked image, with each 40x40 pixel area on the original image corresponding to a single pixel on the lower resolution image.Then, we fit a bicubic spline function to the logarithm of the lower resolution version of the image using scipy.interpolate.griddata.Finally, we created the model by calling the bicubic spline function at the original resolution grid spacing, and taking the reverse logarithm of the resulting image.
To account for the mask affecting adjacent pixels during the process of lowering the resolution and interpolating between pixels, we also applied the spline fitting process to the manual mask.We masked all pixels in the resulting 'mask model' that had values less than 0.9, in effect expanding the edges of the manual mask.

Automatic Point Source Mask
Given our manual mask and smooth model for the light distribution in each galaxy, we then used the photutils.detection.IRAFStarFinder method from the photutils package (Bradley et al. 2022) to automatically detect any point sources in the images that had been missed when masking manually, choosing a threshold of 0.018 counts/s, a background of 0, and a PSF FWHM of 3 pixels.For each galaxy, we masked each detected source using a circular area with a radius of 6 pixels.For efficiency, we only searched for point sources within the 2  region of interest defined above.Point source properties are described further in Appendix B.

Smoothed Image Model
We then created a third, non-parametric model for the smooth light distribution in each galaxy, by simply convolving the masked image with a Gaussian kernel with a 5 pixel radius.This empirical model has the benefit of more accurately capturing any irregular structure in the galaxies, and being more customized to each individual galaxy.

Histogram Mask
Finally, we constructed a third, more sophisticated mask based on a pixel histogram of the galaxy's normalized residual image (NRI; see Section 3.2 for details on the construction of this image), masked with the automatic detection mask.The pixel histogram of the residual image (Fig. 1) includes a Gaussian component at low counts/s, but in some cases also an excess tail at high counts/s.We expect intrinsic surface brightness fluctuations to be Gaussian in the limit where large numbers of stars contribute significantly to the light in each pixel.Since a Gaussian component appears in the NRI distributions for all the objects in the sample, whereas only a few show a significant high-residual tail, for this masking method we assume the Gaussian component corresponds to intrinsic surface brightness fluctuations, whereas the rest of the distribution is produced by other effects.We mask any pixels with residuals in the range where the NRI histogram greatly exceeds the Gaussian component, as follows.
We create the histogram mask by estimating the peak of the distribution to be the bin with the most pixels,   , and then fitting a normal distribution to a region around and below the peak, below a threshold value: (1) This ensures that the fit is not affected by the non-Gaussian tail.If, due to noise, the fitted value of the peak is above  threshold , we iterate until it drops below this value.
Given our fit to the Gaussian component of the pixel distribution, we then mask any pixels beyond the point where the fit intercepts the line  = 1, i.e., where we would only expect one pixel in the bin with this pixel value.This produces the most aggressive of our three masks.
We note that the galaxies in our sample have on the order of 10 3  ⊙ per pixel.From (Tonry & Schneider 1988), and the results of (Cerviño & Luridiana 2006;Cerviño et al. 2008), this seems too low to expect fully Gaussian fluctuations, so it may be that the high-residual tails seen in the NRI for 6-7 of our galaxies are genuinely due to stellar discreteness noise, and the corresponding pixels should be included in the SBF calculation.On the other hand, the fact they generally appear in objects with larger stellar surface mass density makes this seem less likely.As discussed in appendix B , these galaxies also have the largest number of bright point sources, so residuals in the point source subtraction may explain the presence of the tails.

Computing the Power Spectrum
The next step in the SBF method is to calculate the normalized residual image, and determine what fraction of the pixel-to-pixel variation in this image is due to SBFs.In practice, this is usually done in Fourier space, since for a raw CCD image, the power spectrum of the image should include an 'on-sky' component that is convolved with the PSF, and a white noise component from the readout electronics (Blakeslee et al. 1999).
After masking and modelling the galaxy, we create a normalized residual image (NRI).The NRI is calculated as: where we first subtract the model from the image, so that only the fluctuations remain, and then divide by the square root of the model, to normalize the scale of the fluctuations to the expected value.
To compute the power spectrum of the NRI, we create an elliptical aperture around the galaxy and combine this with the mask, to ensure that any measured SBF variance is coming from the galaxy rather than from other sources.
Since SBFs occur in the real image of the galaxy on the sky, they will be convolved with the telescope's PSF.Taking the Fourier transform of the image converts this convolution to a multiplication, so we can simply divide the power spectrum of the NRI, minus any white noise component due to readout noise, by the power spectrum of the PSF, to estimate the amplitude of the SBFs.
Two complications arise.First, COSMOS mosaic images have been combined from multiple raw exposures, using drizzling with an interpolation kernel for geometric corrections, so the white noise component is not actually 'white', but has a slight dependence on wavenumber (as discussed further in section 5.3; see also Mei et al. 2005).To correct for this, we estimated the power spectrum of the background from blank regions in the image around the galaxy, and subtracted this empirical form when fitting for the SBFs.
A second complication is masking.This has a multiplicative effect on the image in real space, and thus corresponds to a convolution in frequency space.To account for this, we multiply the SBF component of the power spectrum by a Fourier transform of the PSF, convolved with a Fourier transform of the full mask (taken to be the product of the manual mask, the automatic point source mask, the histogram-based mask, and the elliptical aperture described above).All Fourier transforms were computed in 2D using the Python function numpy.fft.fft2;azimuthal averages were then taken for each bin in wavenumber.
To compute the PSF, we used eight non-saturated stars from across the COSMOS field.We took circular areas of radius 25 pixels around each star and masked the rest of the image.We then computed the 2D power spectrum of each star and took the average of all eight power spectra as the mean power spectrum of the PSF.Fig. 2 shows the resulting 1D power spectra of the eight individual stars, as well as the mean PSF.
To compute the power spectrum of the background, we used four nearly empty areas (i.e. with only minimal masking required) in the COSMOS field, ranging in size from 15-20 ′′ (500-667 pixels) on a side.We modelled each area using imfit's FlatSky function, which only has one parameter, the surface brightness of the sky, to calculate a normalized residual image for each.Calculating the power spectrum of each and normalizing it to have an integral of 1, we then interpolated the resulting background spectrum to have the same sampling in  as each individual galaxy image.Finally, we took the mean of the four normalized and interpolated spectra to get the power spectrum of the background for each galaxy.Fig. 3 shows the resulting normalized 1D power spectrum of the background.

Noise Estimate
We have found experimentally that the shape and amplitude of the NRI power spectrum change slightly with the choice of aperture size and position, so to avoid creating any bias from a particular choice of these parameters, we created 50 different power spectra for each galaxy, randomly selecting the size of the semi-major axis  and the position of the centre ( 0 ,  0 ) from uniform distributions:   ≤  < 2  , and  0 − 5 ≤  <  0 + 5,  0 − 5 ≤  <  0 + 5, where   is the effective radius from imfit and  0 and  0 are the coordinates of the aperture centre in pixels.We held the ellipticity and angle constant, using the best-fit imfit values for these parameters.

Estimating the SBF Variance
The SBF variance is usually calculated by fitting the azimuthallyaveraged power spectrum of the masked NRI, (), to the equation (Tonry et al. 1990): where  0 , the quantity of interest, is the SBF variance,  1 is the white noise variance, and  is the spatial frequency in units of inverse pixels (e.g.Greco et al. 2021). () is known as the (azimuthallyaveraged) expectation power spectrum, and is a 1D average of the two-dimensional convolution of the power spectra of the PSF and the mask (Greco et al. 2021): We used astropy.convolution.convolve_fftto compute the expectation power spectrum, given the power spectra of the PSF and the mask, and normalized it to have an integral of 1 before fitting.
As mentioned above, however, we found that the white-noise component of the spectrum actually has a dependence on  (see Section 5.3).To account for this, we modified equation (3) to include the power spectrum of the background variance: where  () and () are the expectation power spectrum and the power spectrum of the background respectively, computed as described above.We fit to equation (5) using scipy.optimize.curve_fit,which uses a non-linear least squares algorithm.Given the model from equation ( 5), it returns the best fitting  0 and  1 values and the covariance for these.Note that when fitting, we use units of inverse pixels for  for convenience, i.e. we divide the standard wavenumber  by the number of pixels on one side of the image.
In fitting the expectation power spectrum, we limit the range of wavenumbers considered.Very low wavenumbers (large spatial scales) are affected by the smoothing used to make the empirical model of the smooth light distribution, while high wavenumbers are compromised by correlated noise between pixels (Greco et al. 2021) and/or errors in our background variance estimation.As discussed in Section 5.5, we tested many wavenumber ranges for the entire sample, finding the range 0.1-0.3pix −1 to be optimal.561851, 401988, 733922, 686606, 213165, 259971, 709026, 918161, 458976, 331749, 880547, 380820, 589205, 279307, 377112, 824852.

Computing the Fluctuation Magnitude and the SBF Distance
Each of the 50 SBF runs for a single galaxy returns a value of  0 , which needs to be divided by the number of unmasked pixels, , as  0 is the sum of pixel variances (Tonry & Schneider 1988;Greco et al. 2021).We computed the mean and standard deviation of the 50  0 / values rather than the mean of the SBF distances to estimate the uncertainty because some of the mean  0 / values are negative, mapping to an infinite distance.This can happen when the data is noisy and the background term dominates equation ( 5).We calculated the apparent fluctuation magnitudes from: where  zpt is the zero point magnitude of the original image.
To convert the apparent fluctuation magnitude to a distance, we need an independent estimate of the calibration magnitude  SBF .We used the calibration developed by Carlsten et al. (2019), that was calculated using tip of the red-giant branch distances for low surfacebrightness dwarf galaxies. SBF is calculated using the colour of the galaxy, because the amplitude of SBF is dependant on the stellar luminosity function of the galaxy, and thus on its stellar population.The age and metallicity of the population can be accounted for using an optical colour, thus the net dependence of  SBF on colour (Worthey 1994;Carlsten et al. 2019;Polzin et al. 2021).
We tested two separate methods of computing the  −  colour needed for the calibration, calculating this either from the COSMOS photometry, or from the images directly.See Section 5.6 for a comparison of the two methods.For the first method, we estimated the  −  value using the COSMOS 2015 catalogue's 3 ′′ (corresponding to 100 pixels in the F814W imaging) aperture magnitudes in the Subaru bands  + , IA464, and IA484.This aperture was selected, as the automatic aperture magnitudes give noisier colours (Laigle et al. 2016).The g-band magnitude is not available, so we estimated it as: The second method involves the same  −  estimation but using the Subaru  + , IA464, and IA484 images directly.For each galaxy, we created an approximate g-band image using equation ( 7), and used this to make  −  image.We then applied to this colour image an elliptical aperture identical to that used on the galaxy.We took the average of all the values within the aperture as the estimated  −  colour, and the standard deviation as its error.
To Finally, the SBF distance was calculated using:

RESULTS
The derived SBF distances were compared to distances based on the redshifts described in Section 2. To correct the latter for the effects of nearby structure, including the Virgo Cluster, the Great Attractor, and the Shapley Supercluster, we also derived flow-corrected proper distances based on the model of (Carrick et al. 2015), by looking up the CMB-frame velocity in NED and then using the on-line tool at http://cosmicflows.iap.fr, with the cosmological parameters (Ω m ,  0 ) = (0.315, 67.5 km s −1 Mpc −1 ).(Note CMB-frame velocities are typically ∼ 300 km s −1 larger than heliocentric velocities in the COSMOS field.)We assigned all flow-corrected distances an error of 4.3 Mpc, corresponding to 300 km s −1 for an  0 value of 70 km s −1 Mpc −1 , to account for peculiar velocities.We use the proper distances given in Table 2, for all subsequent comparisons.Our fiducial SBF distances assume the following choices or parameters: • the smoothed image model to compute the NRI • the histogram mask (incorporating the preceding manual and automatic masks) to remove contaminating sources and bad regions • the realistic background SBF equation (equation ( 5)) for background correction • a  −  colour derived directly from the image • a wavenumber range of 0.1 to 0.3 pix −1 when fitting the power spectrum .
Also, we run the SBF method 50 times on each object, selecting a random aperture each time (see Section 3.2), to estimate the uncertainties in  0 and  1 , and propagate the former through the distance calculations.
The results are shown in Fig. 4, where we plot SBF distances versus proper distances, with the dashed grey line representing a 1:1 correspondence.Error bars are derived from the set of 50 variants on the fiducial aperture described above (±1 in  0 /), and probably underestimate the true errors; these are hard to define precisely, as they depend on the range of choices we are willing to make in the method.in  0 / .Objects with a finite lower limit but an infinite upper limit are plotted as upward-pointing arrows.Note the labels in the legend are sorted by increasing distance, as in Fig. 1.
For galaxies at flow-corrected distances of less than 100 Mpc, there is a clear correspondence between the SBF fluctuation amplitude and actual distance.Fig. 5 shows the mean power spectra for galaxies grouped into four bins by proper distances: 21.3-49.3Mpc (8 objects), 49.3-77.4Mpc (5 objects), 77.4-105.5 Mpc (1 object), and 105.5-133.6Mpc (6 objects).Individual power spectra (within 1  aperture) were normalized by the number of unmasked pixels and averaged in each bin.We see a clear trend of decreasing power with distance, except in the fourth bin, which has a mean power spectrum similar to the third.Since we expect  0 / to scale linearly with distance, we can normalize the results by distance by multiplying the power spectrum for each object by its distance before averaging; once again, Fig. 6 shows that the first three bins have comparable intrinsic power, while the last bin does not.
Finally, Fig. 7 shows residual NRI power spectra, after subtracting the fitted background component and dividing by the fitted SBF/Mask  power spectrum: As before, the results for each galaxy have been averaged in four distance bins.If our modelling and fit are correct, the residuals should be constant with  = 1 (black dashed line).This is appears to be somewhat true in the range  = 0.13 -0.33, but outside this range we see significant positive and negative variations.(We also note that the third bin is very noisy because it contains only one object.)From all this, we conclude that SBF power in our sample scales as expected out to ∼60-100 Mpc, but that beyond that range, some other component that is independent of distance dominates the power spectrum.At smaller distances, there is an additional problem with the zero-point of the relation.All of our sample points lie above the 1:1 line in Fig. 4, indicating that we are systematically underestimating the SBF power in all of them.Fig. 8 illustrates these two problems for individual galaxies; distant galaxies generally have more power than expected, while nearby  galaxies have only a third the power expected.We will explore these problems further in Section 5.

TESTING VARIATIONS ON THE SBF METHOD
There are a number of choices to make when proceeding through the steps in the SBF method.They include the type of smooth model used, the type of mask used, how the background contribution is subtracted, the size of the aperture applied, the range of wavenumber values to fit, and the optical colour to use in the calibration.In this section, we discuss how we made each of these choices, and show how variations on our fiducial choices affect the results.

Modelling the Smooth Light Distribution
In Section 3.1, we described three different models of the smooth light distribution in the galaxy: a Sérsic profile created using imfit, a bicubic spline model, and a smoothed-image model created by convolving the image with a Gaussian kernel.The model of the light distribution is used when creating the NRI, using equation (2).
Using the parametric imfit model, we found that the NRI contains a lot of large-scale residual structure (Fig. 9).This is because the surface brightness profiles of the dwarf galaxies in our sample are generally too complicated to model with a simple Sérsic profile.Like the imfit model, the bicubic spline model also leaves behind a lot of residual structure in the NRI (Fig. 10).We also found that many of the NRI pixel histograms created using the bicubic spline models (Fig. 11) do not show a clear Gaussian component, and do not centre around 0 as expected.Furthermore, fitting the distribution for the furthest galaxy in the sample, 824852, did not produce a convergent result, and thus we could not obtain a SBF distance for this object.
By comparison, the smoothed-image model leaves much less residual structure in the NRI, as it is more customized to shape and gross features of the galaxy (Fig. 12).In addition, the pixel value distributions of the NRIs created using the smoothed-image model (Fig. 1) show a Gaussian component and centre around 0 as expected, assuming fluctuations are in the large- limit (see Section 3.1.6).Thus, we have adopted the smoothed-image model for our fiducial analysis, although we still use the Sérsic parameters given by imfit to position and scale the apertures around each galaxy.

Masking
As discussed in Section 3.1, there are several ways to create a mask.Each represents a compromise between removing all contaminating sources and bad regions, and preserving genuine SBFs.Of our three masking options, the histogram mask is the most effective because it captures the most contaminating sources.
The initial manual mask removes only the most obvious contaminants, while the automatic point source detection mask removes more point sources, but leaves noticeable oscillations in the power spectrum of the masked NRI (Fig. 13).These appear to be particularly strong in objects with large numbers of bright central point sources.We have determined experimentally that the oscillations originate from the point-source masking itself.Using the same aperture size for all the point sources detected introduces structure on this spatial scale, and thus produces periodic oscillations in the power spectrum, which will degrade the quality of the SBF distance estimate, as discussed further in appendix C.
To remove these features, an even more aggressive approach is required, so we introduce the histogram mask.This mask removes  Power spectra of the masked normalized residual images using different masks.These are created using the smoothed image models, and a 1  aperture around the galaxy.Each subsequent masking method includes the previous method(s).Galaxies are ordered as in Fig. 1. the oscillations seen in the automatic point-source detection mask spectra, and in some cases, covers all the region of the automatic mask, meaning a histogram masking approach could be effective on its own without further point source removal.Thus, we adopt the histogram mask (which includes our previous two masks) in our fiducial method.3) -red dashed curve), and a background determined directly from blank areas of the field (cf.equation ( 5) -blue dashed curve).

Background Subtraction
Examining our NRI power spectra, we found that the power does not become constant at large wavenumbers, but continues to decrease down to the Nyquist frequency.This is an artifact of our images, which have been combined from multiple pointings and exposures, correlating and smoothing pixel noise over a range of scales (Mei et al. 2005;Mitzkus et al. 2018).Thus, rather than fitting the NRI power spectrum using the original SBF equation (equation ( 3)), we estimated the background from blank regions in the COSMOS field, as described above in Section 3.2, and fit using equation (5).
To illustrate the differences between the two approaches, we created an average power spectrum by combining the individual NRI spectra of the 5 nearest galaxies, normalized to the same integrated power, and fit to the average using both background subtraction methods.In place of the expectation power spectrum, we used the power spectrum of the PSF alone, since the expectation power spectrum is only slightly different from the power spectrum of the PSF, and is different for each object since they have different masks.The fit was performed over the wavenumber range 0.1-0.3pix −1 .
Fig. 14 shows the the averaged masked NRI power spectrum (solid black curve), along with the two fits.It is clear that the fit using equation (3) (red dashed curve) does not accurately capture the shape of the power spectrum.Accounting for a realistic background contribution using equation ( 5), however (blue dashed curve), gives a very good fit to the spectrum, justifying its use in our fiducial method.This approach is a significant improvement on the assumption of white noise, as it gives a much better fit to the data.

Dependence on Aperture Size
The choice of aperture size to use when measuring SBFs is unclear a priori.Too large an aperture will dilute the signal with a larger background contribution, while too small an aperture may leave too little area to measure the SBFs.
As explained previously, to account for some of the uncertainty introduced by this choice, we run the SBF measurement process 50 times, choosing a random aperture size each time, taking the final value of  0 / to be the mean of the 50 resulting values, and its error to be the standard deviation of the distribution.Fig. 15 shows the value of  0 / as a function of aperture size for each galaxy.In most cases, the variations seem small and/or random, but in some (e.g.733922), we see indications of a systematic decrease in  0 / with aperture size.This could indicate excess power from point sources or smooth model errors in the centre of the image, and we have used these plots iteratively to test the effectiveness of our masking choices.

Wavenumber Range
As discussed in Section 3.3, both the low and high wavenumber ends of the NRI power spectra are compromised by various factors, and thus need to be excluded from the fitting.We have tested many different wavenumber ranges, with lower bounds ranging from 0.08 pix −1 to 0.20 pix −1 and upper bounds ranging from 0.20 pix −1 to 0.50 pix −1 .We note that using the smoothed image model removes some large-scale power from the NRI, introducing a dip in the power spectrum at values below ∼ 0.10 pix −1 , so we need to set the lower end of the range to at least this value.
To study the effect of wavenumber range, we plotted the values for  0 / against the lower bound, using a fixed upper bound of 0.30 pix −1 in Fig. 16, and against the upper bound, using a fixed lower bound of 0.10 pix −1 in Fig. 17.For each plot, for efficiency, one point represents 25 rather than 50 SBF runs.Otherwise the process is the same as described in Section 3. A flat curve in these plots means the SBF signal is stable against small changes in the range, which suggests it is a good region in which to select a bound.
From Fig. 16, it is clear that the lower bound must be selected carefully as the SBF signal can vary significantly with the lower bound.We estimate that a lower bound in the region of 0.10 to 0.13 pix −1 is the best choice overall, as this region is stable across most of the galaxies.Fig. 17 shows that the upper bound is much less important than the lower bound, as the plot is generally flat past about 0.30 pix −1 for most galaxies.From this, we infer that the  SBF signal dominates at low values of , and that focusing on lower wavenumbers will give the best SBF fit.For the fiducial method, we choose to fit the wavenumber range 0.1 to 0.3 pix −1 to capture as much of the SBF signal as possible while excluding the dip in the power spectrum below 0.1 pix −1 .

Colour Variation
In Section 3.4, we outlined two methods for estimating galaxy colour.Both assume the same filter conversions, but one uses the catalogue magnitudes, while the other uses the images directly.For the second method, Fig. 18 shows pixel histograms of the  −  images within 1  .It is clear from these histograms that in some cases, there is a large spread in the individual pixel colours relative to the mean catalogue colour, and that the latter may not correspond to the mean of the individual pixel colours.For these reasons, we choose to use the pixel-based approach in our fiducial method.The spread in some of the colour distributions (e.g.824852) also suggests we should treat with caution SBF distances derived assuming a single mean colour, as the colour can vary widely across the extent of the object.An improvement on the method might be to account for the changes in colour across the object, which could be explored in future work.

DISCUSSION
Dwarf galaxy populations in the local Universe represent a challenge and an opportunity to improve our understanding of galaxy formation, feedback, and dark matter physics.To identify local dwarfs and place them in their environmental context, distance estimates are essential, yet obtaining redshift-based distances for large numbers of faint, low-surface brightness dwarfs remains challenging.With a new generation of imaging surveys forthcoming, including spacebased surveys with Euclid and Roman, it is worth reconsidering the potential for other methods such as SBF to tackle this problem.
The SBF method has traditionally been applied to relatively bright, high surface-brightness early-type galaxies with relatively simple stellar populations (Moresco et al. 2022), but in recent work, several groups have pushed to apply the technique to fainter systems within ∼ 25 Mpc (Carlsten et al. 2019;Polzin et al. 2021;Kim & Lee 2021;Greco et al. 2021;Carlsten et al. 2022).In this paper, we have tested the possibility of applying the SBF method to a set of 20 dwarf galaxies at distances of 20-150 Mpc, drawn from the HST/ACS imaging of the COSMOS survey.There are several important challenges in getting the method to work for this sample, including the overall faintness and low surface brightness of the galaxies, internal gradients or scatter in stellar populations and/or extinction, irregular morphology, and a correlated background component introduced by mosaicing.
After making careful choices in how the SBF method is applied to the images, we find a Gaussian SBF signal that does correlate linearly with distance, out to distances of 50-100 Mpc.On the other hand, the recovered SBF power is only ∼ 1 /3 of that expected, leading to a corresponding overestimate in distance.
Beyond 100 Mpc, the measured SBF power is generally too high, leading to underestimates in distance (e.g., Blakeslee (1997)).Examining the point source luminosity functions of our most distant galaxies (Appendix B), we infer that this excess power comes from faint point sources undetected in the relatively shallow COSMOS imaging.In principle, deeper observations that better resolve the point source luminosity function would overcome this limitation.
Estimating SBF distances also requires a number of (only partly constrained) choices in the detailed method, including how to model the galaxy surface brightness profile, how to identify and mask contaminating sources, how to estimate and subtract the background contribution to the NRI power spectrum, the colour used in computing the SBF calibration magnitude, and the wavenumber range used in fitting the power spectrum.For the galaxies in our sample, many possible choices in these details lead to systematic variations significantly larger than the simplest errors we have estimated by varying the aperture around each galaxy slightly.Several choices increase the average SBF power for nearby galaxies until it is closer to the expected level, but in doing so produce oscillations in the NRI power spectrum, non-Gaussian features in the NRI pixel distribution, larger uncertainties, and/or reduce the correlation between detected power and actual distance.
Recent work by Polzin et al. (2021) estimated a SBF distance to the nearest galaxy in our sample, 549719.They obtained a result in agreement with the proper distance (  = 21.3Mpc), SBF = 24 ± 3Mpc (their full galaxy measurement, which averages over two regions with slightly different colours).By comparison, with our fiducial method, we obtain a SBF distance of  SBF = 56Mpc with lower and upper limits of 50Mpc and 67Mpc respectively.Thus, our fiducial estimate is incorrect by a factor of almost 3, and inconsistent with the Polzin et al. (2021) result.We can recover their result by modifying our methodology to make it more similar to theirs.If we use a single Sérsic profile to model the smooth light distribution, flat background subtraction assuming the standard SBF equation (Eq.3), a less aggressive mask (using only a manual mask modelled after their figure 2), and a wavenumber range of 0.15 to 0.4 pix −1 , we obtain a SBF distance of 23 ± 5 Mpc, in excellent agreement with their result.These are not optimal choices for the rest of our sample, however; if we apply this modified method to the whole sample, we find all our galaxies have estimated SBF distances in the range of ∼10-30 Mpc, and there is little or no correlation between  SBF and   .
There are a few ways we may imagine making progress on the problem of generalizing the SBF method to a more diverse, distant, low surface-brightness populations of galaxies.With larger samples of partially resolved objects, we might be able to develop a new calibration for our current fiducial method, that corrected for the zeropoint offset seen in Fig. 4. We note that with a different calibration, most objects in this figure with true distances of less than 100 Mpc might have SBF distances consistent to within the (10-15%) statistical errors shown in this plot.
We also have several important pieces of additional information about the galaxies in our sample, including their associated point source luminosity functions (see Appendix B), their SEDs (see Appendix A), and the detailed shape of the pixel histogram (Fig. 1).These diagnostic features may allow corrections to the SBF method for age, metallicity, galaxy type and/or star formation history.More general statistical methods that incorporate the additional information, such as those used in machine learning (ML), might be be applied to this problem, in order to determine the best choices to make in applying the SBF method, and which factors are most important in the analysis (see Tanoglidis et al. 2020;Müller & Schnider 2021, for recent examples of the application of ML to the morphological classification of low surface brightness dwarf galaxies).Such methods might also be used to identify and correct the NRI in distant galaxies where undetected point sources contribute significant power, as discussed in Appendix B. We will explore these possibilities in future work.
Given the imaging data expected from forthcoming surveys (Moresco et al. 2022), extending the SBF method to a broader class of galaxies is a priority.The Euclid Wide Survey (Euclid Collaboration et al. 2022), for example, will cover 15,000 square degrees on the sky, almost 10 4 times the effective area of the COSMOS ACS imaging, albeit with lower resolution (∼ 0.2") and less collecting area.If the 20-30 galaxies detected in the COSMOS field at distances of less then 200 Mpc are indicative, Euclid may see thousands of partiallyresolved dwarf galaxies out to ∼50 Mpc.The Roman Space Telescope (Spergel et al. 2015), with 100 times the field of view of HST and better infrared sensitivity, could produce samples of thousands of objects out to distances of 200 Mpc or beyond.From the ground, the Vera Rubin Observatory (LSST Science Collaboration et al. 2009) will similarly produce large samples of partially resolved galaxies over the distance range of the COSMOS sample (Greco et al. 2021).These samples could provide a trove of information about galaxy formation on the smallest scales, if we can identify these nearby objects and correctly place them in their three-dimensional context.

APPENDIX B: POINT SOURCE LUMINOSITY FUNCTIONS
In addition to SBFs, most of the galaxies in our sample also contain populations of bright point sources, which may be globular clusters, unresolved HII regions, young star clusters, or even individual young stars.These sources contribute additional power to the NRI, unrelated to SBFs, so it is important to identify and mask them.We detected sources using photutils, as described in Section 3.1.The resulting luminosity functions are shown in Fig. B1, both in terms of apparent magnitude (left-hand panel), and in terms of absolute magnitude (right-hand panel, using a distance modulus based on the proper distance).The effective depth of the luminosity functions varies due to differences in the background and smooth galaxy model, but is generally in the range  814 = 26-27.5,judging from the slope of the counts.For the nearer galaxies in our sample, this corresponds to absolute magnitudes −6 to −4, around the bottom of the globular cluster luminosity function, and the beginning of the red giant branch (RGB) whose stars contribute strongly to SBF.Carlsten et al. (2019) have argued that masking to this depth should eliminate essentially all globular clusters expected in dwarf galaxies.On the other hand, Kim & Lee (2021) recently tested the effect of point source masking to greater depths on a sample of nearby dwarfs, observed with Hyper Suprime-Cam.Masking all point sources down to a depth of   = −3.5 to −5, they still find significant variations in the SBF absolute magnitude over this range, at least for the bluer part of their sample ( −  < 0.5).Thus, for the nearest galaxies in our sample, it seems we have just enough depth to account for point source contamination, and may fall a bit short of this for blue galaxies.
On the other hand, for the most distant galaxies located at 100 Mpc or more, the point source depth of the COSMOS imaging corresponds to absolute magnitudes −10 to −8, suggesting these are bright globular clusters or HII regions, and many fainter point sources remain undetected.We note that the 3-4 most distant galaxies are also those for which we measure significantly more power in the NRI (see Fig. 8).It seems likely that this power ("excess", relative to the power measured in the rest of the sample, though the resulting total power is close to the SBF expectation) comes from undetected fainter point sources in these objects.
While we cannot detect and mask out these fainter sources without deeper data, close examination of the apparent luminosity functions suggests we may be able to identify galaxies where undetected point sources are a problem.For the nearest galaxies in the sample, we detect many point sources, and thus can remove much of the associated power from the NRI with masking.For some more distant galaxies, the only detected point sources are very faint, and thus any undetected sources will contribute relatively little power.The difficult cases seem to be those where we detect 5-10 fairly bright point sources, but few fainter ones (e.g.galaxies 824852, 37712, 279307, and 589205).This suggests that with further consideration of the shape of the point source luminosity function, together with colours, radial distributions, and/or general information about the galaxy (magnitude, surface brightness, concentration, etc.), we might be able to distinguish nearby galaxies from more distant ones, e.g. using a machine-learning approach trained on a larger sample, and correct for undetected point sources statistically (e.g.Tonry & Schneider 1988;Mei et al. 2005).

APPENDIX C: EFFECT OF ALTERNATE MASKING STRATEGIES
To determine the origin of the oscillations in the masked NRI power spectra, we created a mask for a single automatically-detected point source (one circle with a 6 pixel radius).We then compared the power spectrum of this single-source mask with the power spectrum of their full mask, as well as the power spectrum of their masked NRI (Fig. C1).We see clearly that the oscillations in the power spectrum of the single-source mask correspond closely to the oscillations in the other power spectra.Thus, we conclude that these oscillations originate from the automatic point source mask.While our method includes this point source mask in the expectation power spectrum, we suspect that noise in the power spectrum around the oscillations leads to poor fitted values for the automatic point-source mask (cf.C4).
We tested two alternate versions of the automatic mask, to try to remove the oscillations in the spectra.First, we simply smoothed the automatic mask using the same Gaussian kernel used for the smoothed image model (see section 3.1.5).This was not effective, as it did not fully mask the point sources, i.e. the masked areas around these objects were no longer set to 0. A second method was to multiply this smoothed mask by the original mask, effectively smoothing the edges of each masked region.
The first smoothed mask was effective at removing the oscillations in the power spectra, while the smoothed edge mask still showed some oscillations (Fig. C3).The smoothed mask produced SBF distance results closer to the proper distance, and no longer consistently overestimates the distance of nearby objects (see Fig. C4).We consider this result is misleading, however, as the extra power available to decrease the SBF distance estimate is clearly related to point sources, and not genuine SBFs.Given the limitations of these versions of the automatic mask, for our fiducial method we chose to use the histogram mask.

APPENDIX D: VARIATIONS ON THE BICUBIC SPLINE MODEL
When creating the bicubic spline models, we tested different scales for the low-resolution image, binning together 20x20, 25x25, 30x30, 35x35, and 40x40 pixel areas.Using the same choices as our fiducial method, but the bicubic spline model for the smooth light distribution, we obtain the results given in Fig. D1.The results are generally similar for all the resolutions considered.The nearest galaxies tend to be clustered at similar SBF distances, and the results are less correlated with proper distance compared to Fig. 4.
As described in Section 3.1.3,we also tested using an expanded manual mask to remove pixels contaminated by the mask when producing the lower resolution image.Fig. D2 shows the results for various resolutions, with this additional masking.The additional masking reduces the uncertainties in some cases, but generally the results are similar to those in Fig. D1.Since the SBF distances show less correlation with proper distance, relative to our fiducial method, and the NRI histograms are less Gaussian, as described in Section 5.1, we have not adopted the bicubic spline modelling method in our main analysis.

Figure 1 .
Figure 1.Pixel value distributions in the Normalized Residual Images (NRIs -coloured histograms) created using the smoothed image model.Note that in this and most subsequent plots, the sample is ordered by distance from top left to bottom right.The thick points show the bins used to fit the Gaussian component, while the black curves and vertical lines show the final fit and fitted peak value respectively.The histogram mask excludes pixels with values in excess of the Gaussian expectation (i.e. the tail to the right of the Gaussian component).Galaxies are ordered from left to right as follows: 424575, 549719, 677414, 260583,  561851, 401988, 733922, 686606, 213165, 259971, 709026, 918161, 458976, 331749, 880547, 380820, 589205, 279307, 377112, 824852.

Figure 2 .
Figure 2. Power spectra of the eight non-saturated stars used to estimate the PSF (thin coloured curves), as well as the power spectrum of the mean PSF (thick black curve).

Figure 3 .
Figure 3. Power spectra of the four empty regions used to estimate the power spectrum of the background (thin coloured curves), as well as the mean power spectrum of the background (thick black curve).

Figure 4 .
Figure 4. SBF distance  SBF versus proper distance   .Objects with upper and lower limits are plotted with vertical error bars corresponding to ±1 in  0 / .Objects with a finite lower limit but an infinite upper limit are plotted as upward-pointing arrows.Note the labels in the legend are sorted by increasing distance, as in Fig. 1.

Figure 6 .
Figure 6.Average NRI power spectra for each of the four bins in proper distance, with the individual power spectra multiplied by their distance before averaging.

Figure 7 .
Figure 7. Residuals of the average NRI power spectra for each of the four bins in proper distance, after subtracting the fitted background component and dividing by the fitted SBF component.The dashed line indicates the theoretical expectation,  = 1.

Figure 8 .
Figure 8. Measured SBF power  0 / versus expected power, as a function of expected power.

Figure 9 .
Figure 9. Normalized residual images created using the imfit models.Each image has a 2  aperture applied around it, to focus on the residuals left behind after fitting.Galaxies are ordered as in Fig. 1.

Figure 10 .
Figure 10.Normalized residual images created using the bicubic spline models.Each image has a 2  aperture applied around it.Galaxies are ordered as in Fig. 1.

Figure 11 .
Figure 11.Pixel value distributions in the NRIs created using the bicubic spline model.Galaxies are ordered as in Fig. 1.

Figure 12 .
Figure 12.Normalized residual images created using the smoothed image models.Each image has a 2  aperture applied around it.Galaxies are ordered as in Fig. 1.

Figure
Figure13.Power spectra of the masked normalized residual images using different masks.These are created using the smoothed image models, and a 1  aperture around the galaxy.Each subsequent masking method includes the previous method(s).Galaxies are ordered as in Fig.1.

Figure 14 .
Figure 14.Average power spectrum of the masked normalized residual image for the 5 nearest objects (black curve), along with fits assuming a flat white-noise background (equation (3) -red dashed curve), and a background determined directly from blank areas of the field (cf.equation (5) -blue dashed curve).

Figure 15 .
Figure 15.Dependence of  0 / on aperture size.The values of  0 / were computed using the smoothed-image model, with the histogram mask, and a fitting range  = 0.1 to 0.3 pix −1 .Galaxies are ordered as in Fig. 1.

Figure 16 .
Figure 16.The normalized SBF signal plotted against the lower  bound, given a fixed upper bound of 0.3 pix −1 .The error bars represent the standard deviation on the values obtained using 25 randomly selected apertures.Galaxies are ordered as in Fig. 1.

Figure 17 .
Figure 17.The normalized SBF signal plotted against the upper  bound, given a fixed lower bound of 0.1 pix −1 .The error bars are as in Fig. 16.Galaxies are ordered as in Fig. 1.

Figure 18 .
Figure 18.Pixel histograms for the  −  images.The dotted vertical line represents the catalogue value and the dashed line represents the mean of the histogram values.The mean value and its corresponding standard deviation are what we choose to represent the  −  value of the galaxy.Galaxies are ordered as in Fig. 1.

Figure A1 .
Figure A1.SEDs for each galaxy in the sample, normalized to the value at 16311Å.Galaxies are ordered as in Fig. 1.

Figure B1 .
Figure B1.Luminosity functions for each galaxy in the sample.Left: in terms of apparent magnitude.Right: in terms of absolute magnitude, calculated using the proper distance.

Figure C1 .
Figure C1.Power spectrum of the automatic point source mask, for each of the subset of the sample with the strongest oscillations.The power spectrum of a single point-source mask (solid pink curve) is shown for comparison.

Figure C2 .
Figure C2.Power spectrum of the NRI, for each of the subset of the sample with the strongest oscillations.The power spectrum of a single point-source mask (solid pink curve) is shown for comparison.

Figure C3 .
Figure C3.NRI power spectra derived using the normal automatic mask, smoothed edge automatic mask, and smoothed automatic mask.Note that the systems with the largest oscillations in the automatic mask are also those with the largest extra-Gaussian components in their pixel distribution (see Fig.1).Galaxies are ordered as in Fig.1.

Figure C4 .
Figure C4.SBF distance estimates derived using the normal automatic mask (left panel), the smoothed-edge automatic mask (centre panel), and the smoothed automatic mask (right panel).A vertical dashed line indicates an object for which SBF were not detected successfully.

Figure D1 .
Figure D1.SBF distance estimates derived using the bicubic spline model.From top to bottom, the pixel grid sizes are 20x20, 25x25, 30x30, 35x35, and 40x40 corresponding to a single pixel on the lower resolution image.

Figure D2 .
Figure D2.SBF distance estimates derived using the bicubic spline model and the expanded manual mask.From top to bottom, the pixel grid sizes are as in Fig. D1.

Table 1 .
COSMOS galaxy ID, redshift distance, coordinates, and image size

Table 2 .
COSMOS galaxy ID, redshift distance, and estimated proper distance

Table A1 .
Galaxy properties.Position angle  is given in degrees counter clockwise from the positive y-axis.

Table A2 .
SBF results.SBF magnitudes are given in the F814W band.