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

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 $i^{+}<20$ mag, out to distances of $\sim$ 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 $\sim$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 = $cz/H_{0}$, assuming $H_{0}=70{\rm\,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.

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 $\pm 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)¹¹1https://ned.ipac.caltech.edu.

The HST images used in the SBF analysis were HST/ACS F814W-band mosaics (Koekemoer et al., 2007; Massey et al., 2010) obtained from the NASA/IPAC IRSA COSMOS Cutouts Service, while Subaru $i^{+}$, IA464, and IA484-band images (Taniguchi et al., 2015) from the same source were used to estimate the $g-i$ colour of the galaxies.

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).

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$\arcsec$ \textcolorblack(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 $k$ 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.

The SBF variance is usually calculated by fitting the azimuthally-averaged power spectrum of the masked NRI, $P(k)$, to the equation \textcolorblack(Tonry et al., 1990):

where $P_{0}$, the quantity of interest, is the SBF variance, $P_{1}$ is the white noise variance, and $k$ is the spatial frequency in units of inverse pixels (e.g. Greco et al., 2021). $E(k)$ is known as the (azimuthally-averaged) 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_fft to 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 $k$ (see Section 5.3). To account for this, we modified equation (3) to include the power spectrum of the background variance:

where $E(k)$ and $B(k)$ 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 $P_{0}$ and $P_{1}$ values and the covariance for these. \textcolorblackNote that when fitting, we use units of inverse pixels for $k$ for convenience, i.e. we divide the standard wavenumber $k$ 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.3\,\text{pix}^{-1}$ to be optimal.

Each of the 50 SBF runs for a single galaxy returns a value of $P_{0}$, which needs to be divided by the number of unmasked pixels, $N$, as $P_{0}$ is the sum of pixel variances (Tonry & Schneider, 1988; Greco et al., 2021). We computed the mean and standard deviation of the 50 $P_{0}/N$ values rather than the mean of the SBF distances to estimate the uncertainty because some of the mean $P_{0}/N$ values are negative, mapping to an infinite distance. \textcolorblackThis can happen when the data is noisy and the background term dominates equation (5). We calculated the apparent fluctuation magnitudes from:

where $m_{\rm 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 $M_{\rm SBF}$. We used the calibration developed by Carlsten et al. (2019), that was calculated using tip of the red-giant branch distances for low surface-brightness dwarf galaxies. $M_{\rm 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 $M_{\rm SBF}$ on colour (Worthey, 1994; Carlsten et al., 2019; Polzin et al., 2021).

We tested two separate methods of computing the $g-i$ 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 $g-i$ value using the COSMOS 2015 catalogue’s 3$\arcsec$ \textcolorblack(corresponding to 100 pixels in the F814W imaging) aperture magnitudes in the Subaru bands $i^{+}$, 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 $g-i$ estimation but using the Subaru $i^{+}$, IA464, and IA484 images directly. For each galaxy, we created an approximate g-band image using equation (7), and used this to make $g-i$ 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 $g-i$ colour, and the standard deviation as its error.

To use the Carlsten et al. (2019) calibration, we converted the galaxy’s $g-i$ colour to a F475W-F814W colour by solving \textcolorblackequation (3) in Carlsten et al. (2019):

where $x$ is the F475W-F814W colour. Next, we used \textcolorblackequation (4) in Carlsten et al. (2019):

to calculate the SBF absolute magnitude in the $i$-band. Since our imaging is in the F814W band, however, this was converted to this band using \textcolorblackequation (2) in Carlsten et al. (2019):

Finally, the SBF distance was calculated using:

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 \colorblackGaussian 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 \colorblackGaussian component and centre around 0 as expected\colorblack, assuming fluctuations are in the large-$N$ 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.

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 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.

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.3 $\text{pix}^{-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. \textcolorblackThis approach is a significant improvement on the assumption of white noise, as it gives a much better fit to the data.

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 $P_{0}/N$ 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 $P_{0}/N$ 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 $P_{0}/N$ 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.

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\,\text{pix}^{-1}$ to $0.20\,\text{pix}^{-1}$ and upper bounds ranging from $0.20\,\text{pix}^{-1}$ to $0.50\,\text{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 $\sim 0.10\,\text{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 $P_{0}/N$ against the lower bound, using a fixed upper bound of $0.30\,\text{pix}^{-1}$ in Fig. 16, and against the upper bound, using a fixed lower bound of $0.10\,\text{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\,\text{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\,\text{pix}^{-1}$ for most galaxies. From this, we infer that the SBF signal dominates at low values of $k$, 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\,\text{pix}^{-1}$ to capture as much of the SBF signal as possible while excluding the dip in the power spectrum below $0.1\,\text{pix}^{-1}$.

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 $g-i$ images within $1r_{e}$. 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\textcolorblack, 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.