Estimating Sky Level

Detection and surface photometry of faint objects rely heavily on accurately estimating the underlying sky background flux. From scattered light originating from astronomical objects such as the Sun, the Moon, the Milky way, stars, and galaxies, to light pollution from the ground, there are many sources that contribute to the night sky surface brightness (Roach & Gordon 1973). Therefore, all ground-based telescopes encounter the challenge of estimating and subtracting the night sky surface brightness (μ_sky ≃ 21 mag arcsec⁻² at typical "dark" locations) in order to access the far smaller flux levels characteristic of faint galaxy halos. Not only ground-based telescopes, but also space-based telescopes must tackle the issue of sky subtraction. The sky surface brightness measured by the Hubble Space Telescope (HST) is 1–2 mag arcsec⁻² fainter than the sky surface brightness measured by ground-based telescopes (Trujillo & Fliri 2016), but contaminating flux sources remain: above the atmosphere zodiacal light, airglow from the Solar wind, and excitation of residual propellant gas from spacecraft all contribute to the sky background.

The proper sky level can be different for detection of objects than it is for the optimal measurement of photometry. This is due to the fact that faint, unresolved and undetected objects underlie the object for which photometry is desired. This is true for stars as well as galaxies. Unbiased sky estimation has been attempted widely in the literature: prominent examples include FOCAS (Tyson & Jarvis 1979), DAOPHOT (Stetson 1987), SExtractor (Bertin & Arnouts 1996), SDSS Photo (Lupton et al. 2002), GALFIT (Peng et al. 2010), PyMorph (Virkam et al. 2010; Bernardi et al. 2017; Fischer et al. 2017), the LSST Data Management Stack (LSST Stack hereafter, Bosch et al. 2018), Huang et al. (2018), and Jenness (2015). The problem of using biased sky background around detected objects is typically encountered on scales that are large compared with the point-spread function (PSF), where the pixel counts from the object become indistinguishable from the sky pixel counts. Traditionally, detection of low surface brightness galaxies has relied on background sky estimation precision of one part in 10,000. Extreme dwarf galaxies in the Local Group have mean surface brightnesses as faint as ∼32 mag arcsec⁻² (e.g., McConnachie et al. 2009; Homma et al. 2016). Thus, accurate surface brightness measurements would require a sky unbiased at a level of ∼34 mag arcsec⁻², or about 6 parts-per-million (ppm) of the typical R-band sky level.

Current and upcoming surveys such as the Dark Energy Survey (DES; Flaugher 2005), the Hyper Suprime-Cam (HSC) survey (Aihara et al. 2018), and the Large Synoptic Survey Telescope (LSST; LSST Science Collaboration et al. 2009) are likely to reveal new aspects of galaxies as low surface brightness (LSB) objects (Ivezic et al. 2008; Robertson et al. 2017). The discovery space is large: LSB features can exist on scales of arcseconds to many arcminutes, spanning the majority of faint galaxies at high redshift to more nearby LSB galaxies. Tidal tails have already been detected at surface brightness levels of ∼30 mag arcsec⁻² (e.g., Martínez-Delgado et al. 2010; van Dokkum et al. 2014) and surely exist at lower levels. A relatively unexplored area is the ultra-low surface brightness morphology and tails over a wide range of angular scales at levels of 31–32 mag arcsec⁻². Discoveries are likely at even fainter levels of surface brightness still, which may become accessible in upcoming deep field observations such as the LSST Deep Drilling Fields (hereafter LSST DDFs), which are expected to achieve a coadded 5σ depth of ∼29 magnitude in the r-band filter (LSST Science Collaboration et al. 2009).

Proper sky background estimation is an important tool for studying galaxy formation and evolution, and even more important when estimating galaxy types based on surface brightness profiles. However, it has been challenging to calculate the correct value of sky background. Many automated photometry programs estimate a biased sky background (Bosch et al. 2018; Huang et al. 2018; Jenness 2015). Previous techniques typically mask detected objects and use the remaining pixel values to estimate the background level; however, sky estimates are generally biased high because pixels in the outskirts of detected galaxies survive the masking processes, and undetected low signal-to-noise sources contaminate the background. Accurate sky estimation is thus a prerequisite for photometric studies of faint objects (e.g., LSB galaxies or low-level features around galaxies).

As mentioned above, at these low levels of surface brightness a sufficiently accurate model of camera scattered light must be used for each exposure. Indeed on a wide range of angular scales the sky surface brightness will be dominated by scattered light from bright stars. Such modeling is beyond the scope of this paper, instead we focus on the challenge of sky bias introduced by the detected object's faint outer halo and by the high density of undetected galaxies. Thus, suppressing sky bias from these two known effects is a necessary but not sufficient condition for ultra-low surface brightness photometry (LSST Science Collaboration et al. 2009).

Sky background estimation directly impacts astronomical object detection. Detection of both stars and galaxies requires accurate characterization of the sky background. As an example, photometry of faint stars whose surface brightnesses are very near that of the sky is described in Stetson (1987). Typically, detection algorithms begin by marking a collection of CCD pixels as belonging to an astrophysical source if they are above some threshold, usually after convolution with a spatial filter optimized for some angular scale. Calculating the flux due to this source (and crucially, this source alone) requires that we quantify the flux those CCD pixels would have in the absence of contaminants, e.g., the sky background. In virtually all sky surveys the flux from unresolved, undetected faint galaxies form a component of this background sky. However, their number and luminosity distribution is known statistically from existing deep surveys. Compilations of deep imaging data provide the number of galaxies as a function of magnitude up to m ≃ 30 for various astronomical filters. For example, Metcalfe et al. (2001) showed that the galaxy count slope in R-band is d(logN)/dm_R ∼ 0.37 for 20 ≲ m_R ≲ 26 and becomes shallower for 26 ≲ m_R ≲ 30. This complete galaxy number count was achieved by compiling a number of observations from both ground-based and space telescopes. Data obtained by a single survey rarely satisfy both a large field of view and a very deep image depth. In deep imaging covering a sufficiently large area (to avoid sample variance) one can statistically expect the same galaxy number counts from any observation at that wavelength. We may thus adopt the well measured mean number of faint undetected galaxies which are responsible for biased sky estimates if the imaging is sufficiently deep.

The idea of correct sky estimation over all angular scales is actually an ill-posed problem. The proper background sky for barely resolved galaxies at high redshift is, in principle, quite different from the correct sky level for large angular scale LSB features. Indeed, the flux from barely resolved galaxies sits on top of the fainter, larger, angular scale flux associated with arcminute scale LSB extragalactic features, which in turn sits on top of the starlight reflected by Galactic cirrus, the zodiacal light, the night sky surface brightness caused by atmosphere emission, and scattered light from bright objects in the camera and the atmosphere. Thus there could be a separate sky estimate appropriate for each of the different morphological classes of LSB objects. To make the problem tractable a multi-component sky model must be built.

The sky model, in principle, can be built using knowledge of the camera and telescope system, locations of bright objects, observational data, and statistical summaries of faint galaxy counts from ultra deep images like the HST. The first step is detecting all objects above a position-variable local sky estimate and masking them. The remaining pixels still contain flux from both undetected galaxies and the faint outer isophotes of the masked detected galaxies which, if left uncorrected, gives an over-estimate of the sky level around compact objects. Because of this, fitting the remaining "sky" pixels with a Gaussian profile, as if it were pure Poisson noise, is incorrect; the distribution of remaining pixels would follow a Gaussian if the pixels contain only the true sky. However, the real distribution has a tail of positive pixels due to the two contributors mentioned above. While 3σ clip and/or one-sided Gaussian fitting improves the estimate, these approaches are arbitrary and lead to a small positive sky background bias (Robertson et al. 2017).

Using the known statistical faint galaxy counts beyond the detection limit together with growing masks around detected objects by a defined amount scaled by total flux help significantly in making these corrections. Indeed, both biases must be removed if the sky is to be correctly estimated at the sub-percent level. This entire process is recursive on every angular scale where there are important sky components. In this paper we focus on the more tractable task of estimating the sky level in the generic case of the extragalactic sky superposed on a slowly varying foreground, thus focusing on the ∼few arcsecond scales associated with typical faint galaxies. To explore the effect of modified masking and accounting for undetected galaxies in a controlled way, we develop a set of simulated images with known inputs and properties.

We begin by describing our image simulations in  Section 2. In  Section 3 we outline the methodology for creating detection and measurement catalogs with SExtractor, where we detail the software specific settings used in this analysis.  Section 4 continues with a discussion of two widely used sky estimators and the techniques they employ. In  Sections 4.3 and 4.4, we present our new sky estimation technique which deals with biases current sky estimators suffer. In Section 5 we examine the importance of accurate sky background estimation in the fitting of galaxy surface brightness profiles. We close this work in Section 6 with a discussion of the difficulties of correctly estimating sky backgrounds, the effectiveness of our algorithm to overcome them, and prospects for future directions.

We use GALSIM (Rowe et al. 2015) to generate galaxy images and sky background. Galaxy and PSF parameters are chosen to be similar to those observed in the R-band imaging data of the Deep Lens Survey (DLS; Wittman et al. 2002). In all, three types of simulations are performed in which we vary the sky level, sky noise, galaxy placement, and magnitude distribution. Each simulation is 8000 × 8000 pixels in area. In Table 1 we list the parameters used in our simulations.

In order to isolate the effect of detection masks from that of undetected faint galaxies, in our first simulation (hereafter referred to as the "uniform distribution simulation") 10,000 galaxies are evenly spaced in the image in a grid pattern, and are not surrounded by any other galaxies inside the sky analysis areas. This ensures the galaxies are well separated and do not contaminate their neighbors with stray flux. As discussed earlier, extended features in detected galaxies often "bleed" beyond the detection mask, contaminating the sky background estimate. As we will show in further detail in Section 4.1, galaxies in the uniform distribution simulation are all detected, meaning there will be no sky background contribution from unresolved or undetected sources. Having completely detected all simulated galaxies, well localized sources will enable us to directly test the impact of growing detection masks, independent of any effects from undetected sources in more realistic images. The number count as a function of magnitude in R band follows a power law with a shallow slope of d(logN)/dm_R = 0.1 over the magnitude range of 19 < m_R < 25 for galaxies in this "uniform distribution" simulation. We choose a shallower slope for the number counts to ensure that our simulated population includes a balanced mix of bright and faint galaxies in our sample of 10,000. A slope of d(logN)/dm_R = 0.1 yields a sample that contains adequate bright galaxies to test our algorithms while minimizing problems due to flux overlaps that may occur from two neighboring bright galaxies on our simulation grid.

The second simulation (hereafter referred to as the "random distribution simulation") contains 793,116 galaxies which are randomly placed in position over the image. As a consequence of their more realistic placement, galaxies may contribute flux to neighboring profiles. Importantly, the random distribution simulation includes a much larger number of galaxies, with an n(m) distribution that extends to a much fainter magnitude limit of m_R = 29. Magnitudes of galaxies in the random distribution simulation follow a power law with a slope of d(logN)/dm_R =0.4. The fainter galaxies with m_R > 25.5 are mostly not detected by the detection algorithms tested in this paper. We call the flux from these undetected galaxies the extragalactic background (EBL). The EBL will contaminate both the estimate of the sky background, and the flux of nearby detected galaxies. These undetected galaxies will be important in Section 4.3. Taken together, the two simulations enable us to test the effect of mask growth in isolation in the uniform distribution simulation, and the joint effects of mask growth and unresolved galaxies in the random distribution simulation.

In Figure 1 we show a small portion of each of the first two simulations generated by GALSIM. The simulated images are intended to reflect the observing conditions in the DLS, with 20 co-added 900 s R band exposures on a 4 m telescope. We constructed 20 simulated images with different random seeds and co-added the images for each simulation. In doing so, we increase the signal-to-noise ratio and limiting magnitude. The simulated sky in both the uniform and random distribution simulations properly emulates conditions of the DLS. The sky surface brightness is μ_sky = 21 mag arcsec⁻² (3240 ADU pixel⁻¹), which corresponds to a single exposure time of 900 s in DLS R band with Poisson noise. We also assume that the sky surface brightness is spatially flat in the uniform and random distribution simulations. The root-mean-square value of sky background in the co-added image is reduced by a factor of the square root of the number of co-added images and becomes σ_sky = 12.73 ADU pixel⁻¹.

Zoom In Zoom Out Reset image size

The range of signal-to-noise ratio is 5 ≲ S/N ≲ 300 for 19 ≤ m_R ≤ 26 at DLS depth. Each model galaxy follows a single Sérsic profile with index ranging from 0.5 to 5. The PSF profiles in the DLS are broader than a Gaussian, and are well described by a Moffat profile (Moffat 1969), which is given by

$$\begin{eqnarray}&&\mathrm{PSF}(R)=\displaystyle \frac{\beta -1}{\pi {\alpha }^{2}}{\left[1+{\left(\displaystyle \frac{R}{\alpha }\right)}^{2}\right]}^{-\beta },\end{eqnarray} \tag{ 1 }$$

where α is the scale length, and β is the slope of the profile. To match the DLS observations, we use Moffat profiles of 0.8 < α [arcsec] < 1 and 2.9 < β < 3.2. Under these parameter ranges, Moffat PSFs are randomly distributed in the entire image for both simulations.

To investigate effects of noise, it is informative to simulate deeper data with higher signal to noise for the target galaxies. To this end, for our third image simulation, we regenerate the uniform distribution simulation, this time to LSST DDF depth, in addition to the simulations to DLS depth. Anticipated typical observing conditions for the LSST DDF are as follows: 10,000 co-added 15 s R band exposures on a 6.7 m (effective aperture) telescope. The sky surface brightness of μ_sky = 21 mag arcsec⁻² is estimated as 1000 ADU pixel⁻¹. The root-mean-square value of sky background in the co-added image is σ_sky = 0.1 ADU pixel⁻¹. Simulating 10,000 simulations with varying random seeds requires a substantial amount of computing time. To impart realistic sky noise fluctuations in our LSST DDF depth image, we take the following steps instead: we re-normalize pixel values of the uniform distribution simulation without sky noise in DLS depth to meet the observing condition for the LSST DDF. We subsequently add the Poisson noise of σ_sky = 0.1 ADU pixel⁻¹ to each pixel.

The correct sky level to be used in detection and photometry can differ. The proper sky level for object detection is the sky underlying all objects, bright and faint. This is true even though the faintest objects are generally not detected and form an unresolved extragalactic background. Any detection algorithm should use the true underlying sky after EBL subtraction. However, current algorithms are not sensitive at the levels discussed above. Developing a new detection algorithm which takes full advantage of the high precision sky estimates is beyond the scope of this paper. For the purposes of the inter-comparisons in this work we use SExtractor. SExtractor (Bertin & Arnouts 1996) is an automated catalog builder used to identify and measure various properties of astronomical objects on a CCD image. We run SExtractor using a detection threshold of DETECT_THRESH = 0.5σ where σ is the root-mean-square sky noise in the entire image, a minimum detection area of 6 pixels, and the number of deblending sub-thresholds of DEBLEND_NTHRESH = 10 with a deblending contrast of DEBLEND_MINCONT = 0.0001. The images are filtered through a 5 pixel × 5 pixel Gaussian convolution kernel with FWHM = 3 pixels. A mesh of 80 pixels × 80 pixels is used to estimate sky background for all identified galaxies by SExtractor. A PHOT_FLUXFRAC = 0.5 is used to estimate the half-light radius for each galaxy. The numbers of cataloged galaxies are 10,000 and 111,409 for uniform and random distributions, respectively. As in all current object detection algorithms, SExtractor fails to detect faint galaxies of low signal-to-noise.

In this section, we compare sky background values estimated by various methods on our simulated images. Because we have a priori knowledge of the true underlying sky brightness that was input into our simulations, we can directly assess the accuracy of these estimators by comparing their results with the truth. To estimate sky background around each of our galaxies, we first run SExtractor to construct a detection catalog of galaxies. By running SExtractor, we obtain SExtractor's sky background estimate at the position of each galaxy. The SExtractor catalog is used as input to GALFIT, which provides a second catalog of background estimations. An overview of the background estimation procedures and our parameters used in SExtractor and GALFIT follows in the next two subsections.

Motivated by the strengths and weaknesses observed in these methods, we develop a new scheme for estimating the local sky background around galaxies. Finally, we couple our local sky estimation technique with a global polynomial background model in each of the simulations to obtain accurate sky background estimates with high precision.

4.1. SExtractor Sky Estimation

SExtractor provides a local sky background estimate. In SExtractor this quantity is estimated by performing an iterative 3σ clip of the pixel values within a user-specified mesh grid that covers the image. We use a mesh of 80 pixels × 80 pixels to estimate sky background. Sextractor considers the cell to be "non-crowded" if σ drops by less than 0.2σ per clipping iteration, and crowded otherwise. Based on these two cases the sky background is given by

$$\begin{eqnarray}\mathrm{sky}=\left\{\begin{array}{ll}\mathrm{Mean} & {\sigma }_{i}-{\sigma }_{f}\leqslant 0.2{\sigma }_{i}\\ 2.5\times \mathrm{Median}-1.5\times \mathrm{Mean} & \mathrm{otherwise}\end{array}\right.,\end{eqnarray} \tag{ 2 }$$

where σ_i and σ_f are the standard deviations of the pixel values in a mesh before and after the 3σ clip, respectively.

4.2. GALFIT Sky Estimation

GALFIT (ver. 3; Peng et al. 2010) is a two-dimensional model fitter designed to model multiple categories of astronomical objects. In the course of measuring a model galaxy GALFIT estimates sky background, which is subtracted in order to find the best-fit parameters for a functional model on a CCD image. GALFIT estimates the sky background at the object's centroid as follows:

$$\begin{eqnarray}\mathrm{sky}({x}_{0},{y}_{0}) & = & \mathrm{sky}({x}_{c},{y}_{c})+({x}_{0}-{x}_{c})\displaystyle \frac{d\mathrm{sky}}{{dx}}\\ & & +({y}_{0}-{y}_{c})\displaystyle \frac{d\mathrm{sky}}{{dy}},\end{eqnarray} \tag{ 3 }$$

where (x₀, y₀) is the centroid of the object in pixel coordinates, (x_c, y_c) is the center of an image cutout and dsky/dx and dsky/dy are gradients of sky background in x and y directions, respectively.

Care must be taken in choosing the background estimation parameters for GALFIT, particularly when choosing an image cutout size. If the size is too small, the image cutout does not include enough sky pixels to make an accurate estimate. However, too large of an image cutout not only requires expensive computational resources, but also results in inaccurate estimation of the sky background due to the increasing number of undetected galaxies (Barden et al. 2012; Vikram et al. 2010). In our study, we adaptively choose the width w, and height h, of the image cutout centered on each galaxy's position. It is crucial that the cutout does not truncate the faint tails of galaxies. To ensure this is the case, we adopt a scheme to conservatively estimate the radius at which the galaxy profile reaches a surface brightness of μ = 30 mag arcsec⁻², R₃₀. We assume all galaxies are described by an n = 4 Sérsic profile and use the half-light radius as measured by SExtractor as the profile's half-light radius. After the mock profile is constructed, R₃₀ can be readily calculated for each galaxy. The width w, and height h, of the image cutout are then defined as

$$\begin{eqnarray}&&w={f}_{\mathrm{img}}{R}_{30}(| \cos \theta | +(1-e)| \sin \theta | ),\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&h={f}_{\mathrm{img}}{R}_{30}(| \sin \theta | +(1-e)| \cos \theta | ),\end{eqnarray} \tag{ 5 }$$

where θ is THETA_IMAGE, e is ELLIPTICITY, and f_img is a free parameter to set the optimal size of an image cutout. We empirically determine that a value of f_img = 2 results in a sufficient number of background sky pixels to determine an accurate estimate of the sky background. We use this same image cutout for galaxy surface brightness profile fitting in Section 5.

4.3. New Sky Estimation Technique

Sky estimation methods which use a sample of local image CCD pixels to estimate the background level at the position of a galaxy can suffer a high bias from two factors: flux from the outer tails of galaxy profiles which extend beyond their respective masks, and flux from undetected (and hence completely unmasked) faint galaxies that reside in the image pixels used to estimate the sky. To deal with these two effects, we develop a new sky estimation technique. This technique employs a two-filter estimator, one spatial and one statistical, to minimize the contribution from pixels coming from the unmasked outskirts of galaxy profiles and from unidentified objects. Our method consists of three high level steps and ultimately yields an estimation of the sky level at the positions of detected galaxies. A flow chart describing the overall process of our sky estimation is shown below in Figure 4.

In the spatial filtering step, we create an updated object mask for each cataloged source. These new masks more effectively exclude flux from the extended tails of galaxy profiles, which previously contaminated the pixels used to estimate the sky background. The procedure for creating new masks was calibrated on the uniform distribution simulation, where galaxies are laid down on a regular grid, with no neighboring galaxies inside the cutout area. The masks are generated by first creating a mock one-dimensional Sérsic profile for each source. Pixels are then masked out if their positions satisfy

$$\begin{eqnarray}&&{{\rm{C}}}_{\mathrm{xx}}{(x-{x}_{c})}^{2}+{{\rm{C}}}_{\mathrm{xy}}(x-{x}_{c})(y-{y}_{c})\\ &&\quad +\,{{\rm{C}}}_{\mathrm{yy}}{(y-{y}_{c})}^{2}\lt {({R}_{30}/{a}_{\mathrm{rms}})}^{2},\end{eqnarray} \tag{ 6 }$$

where a_rms is the 2nd moment along the semimajor axis (A_IMAGE), (C_xx, C_xy, C_yy) is the object ellipse parameter (CXX_IMAGE, CXY_IMAGE, and CYY_IMAGE, respectively) measured by SExtractor, (x_c, y_c) is the centroid of a galaxy, and R₃₀ is the cutout radius (see Section 4.2). We find that using a fainter surface brightness for our masks does not significantly change the remaining pixel statistics. Additionally, the measured SExtractor ellipse parameters are used to assign orientation angles and ellipticities to the masks. We stress that these masks are not meant to perfectly model the profiles of galaxies, but rather effectively mask out their flux. Using an n = 4 Sérsic parameter is sufficient to mask out galaxies which are best described by 0.5 ≤ n ≤ 4 Sérsic profiles.

We then proceed to the statistical filtering step to estimate and subtract the flux contribution from undetected faint galaxies which make up the EBL. To model the EBL the number counts of undetected galaxies in the entire field and their flux must be considered. In the uniform distribution simulation, all galaxies are detected by SExtractor, so accounting for undetected EBL galaxies is unnecessary. In the random distribution simulation, however, SExtractor fails to detect and catalog some galaxies at a true magnitude fainter than m_R ∼ 25.5, as the signal-to-noise for these galaxies approaches the user-set detection threshold. In Figure 2 we show the galaxy number counts for the random distribution simulation. The blue curve shows the histogram of true magnitudes for the objects in the random distribution simulation, which is an excellent match to the the input number-counts slope that was used to generate the mock galaxies, shown in black. The red histogram shows the actual number of detected objects, and the filled sea-green histogram indicates the number of objects not detected by SExtractor, the difference of the red and blue histograms. We exclude a small number of galaxies with observed magnitude m_R < 25.5 when modeling the EBL; residuals in the observed number counts compared to the input number counts for m_R < 25.5 galaxies are due to small measurement errors and the effects of blending. The total number of undetected galaxies with m_R > 25.5, which make up the EBL, is 679,753. Since the number of EBL galaxies is large, for simplicity we assume that the galaxies are uniformly distributed in the field.

Zoom In Zoom Out Reset image size

We use the residual in the observed number counts of galaxies as a function of magnitude compared to their expected value to calculate the total flux of all undetected galaxies. Because we assume EBL galaxies are uniformly distributed across the field, we also assume the total EBL flux is uniformly distributed as well. As a result, we obtain an estimate for the EBL flux per pixel for our simulated data, μ_EBL = 0.898 ADU pixel⁻¹, by simply dividing the total EBL flux by the number of unmasked pixels in the image. This "pedestal" level of flux is then subtracted from each unmasked image pixel to mitigate the effects of the EBL in background estimation. To examine whether the method is sensitive to the exact cutoff in the simulated faint galaxies, we test our method for EBL estimation on an additional image simulation where galaxies are generated up to m_R = 31 with a simple power-law of d(logN)/dm_R = 0.4. We find a background consistent with the added flux from the 29 ≤ m_R ≤ 31 galaxies: the EBL in this case is μ_EBL = 1.670 ADU pixel⁻¹ while the median pixel value is μ_median = 1.593 ADU pixel⁻¹. Thus, the method is not sensitive to the EBL faint end cutoff beyond 30 mag arcsec⁻². Statistical galaxy counts as a function of magnitude are now complete to m_R ≃ 30 (Metcalfe et al. 2001). The slope of the galaxy count at m_R ≃ 29 becomes so shallow that the EBL from the galaxies with m_R > 29 decreases rapidly (Tyson 1995). Because of this, there is little difference in the EBL estimates even though we simulate galaxies following the real galaxy counts. It is therefore safe to simulate galaxies with m_R < 29.

Lastly, we measure local sky background estimates for each galaxy. We select an image cutout centered on the centroid of each target galaxy. The initial width and height of the images are 15 R_e, where R_e is FLUX_RADIUS with PHOT_FLUXFRAC = 0.5 in our SExtractor catalog. If the number of unmasked pixels is less than 4000 or the width (or height) is less than 80 pixels, we iterate by increasing the width and height with an increment of 10 pixels. Once the number of unmasked pixels residing in the image cutout is greater than 4000, the mean of the unmasked pixels is calculated and used as an estimate of the local sky value for the center of the image cutout.

In Figure 3 we compare the performance of background estimation by SExtractor and our method for a particular galaxy. It is known that SExtractor tends to overestimate sky background (Häussler et al. 2007), and the estimate becomes worse when the number of unidentified objects increases. Although the 3σ clip by SExtractor described in Section 4.1 removes excessively bright pixels from consideration when estimating the sky background, the remaining pixels are still contaminated by the flux from undetected galaxies and insufficiently masked galaxies. This can be seen directly in the top right panel in Figure 3, which shows the distribution of unmasked pixel values when using SExtractor's detection masks. While the distribution is well described by a Gaussian, the mean value of the distribution is 1.23 ADU pixel⁻¹ above the true sky value. This bias results from the excess number of pixels in the bright wing of the distribution. In the bottom right panel of Figure 3, we show the distribution of unmasked pixel values after using our new masking scheme. The bias from extended galaxy profiles, which otherwise survives SExtractor's masking procedure, is mitigated with our new masks, as the mean value is μ = 0.910 ADU pixel⁻¹ above the true sky value. However, undetected galaxies continue to pollute the pixels used for sky estimation with excess flux, even when new masks are used. To deal with flux contamination from undetected galaxies, we must also subtract the calculated EBL flux per pixel to obtain an accurate estimation of the sky background (μ = 0.01 ADU pixel⁻¹). The precision of the estimate can be increased by utilizing many such samples of sky over a much larger area and requiring smoothness. For this we assume the sky underlying all galaxies varies on scales much larger than galaxy scales.

Zoom In Zoom Out Reset image size

4.4. Global Background Model Coupled with New Local Sky Estimator

In this sub-section, we describe the procedure to create a global background model for entire simulated images. The background model is created by computing the average background value in semi-local uniformly spaced subsections of the CCD image, and subsequently fitting a smooth two-dimensional polynomial to the average background values. Once created, the global background model can be evaluated at any point in the image to produce a sky estimation (Bosch et al. 2018). Note, this is in contrast to the estimation techniques detailed above. Aforementioned techniques are local estimators of the sky background: using pixels from the surrounding ∼1 arcminute diameter of a galaxy to construct an estimate of the sky brightness in the immediate vicinity of each galaxy. The details specific to creating the global background model (subdividing the CCD, computing the average background values, and polynomial fitting) are discussed below.

We begin by sub-dividing the image into 2500 evenly spaced, equally sized image subsections, where each subsection is 160 × 160 pixels in size (∼40 square arcseconds), and image subsections do not overlap. The centers of these subsections define a 50 × 50 point spatial grid. In each image subsection, an iterative 3σ clip mean and variance are computed on all pixels which do not correspond to detected objects (i.e., the pixels not included in the corresponding SExtractor segmentation maps). The sigma-clipped mean of each image subsection is assigned to its corresponding grid point. The average position of the non-masked pixels in each image subsection is used to place the points for the spatial grid in their respective image subsections. Subsequently, a 6th-order two-dimensional Chebyshev polynomial is fit to the spatial grid. Chebyshev polynomials are more robust to over fitting than spline interpolation, as they are not strictly required to pass through the grid points obtained from the 3σ clip. Each grid point is inverse-variance-weighted in the fit, so image subsections where many pixels are masked have a reduced impact on the fidelity of the fit (Bosch et al. 2018). The fitted polynomial model may be used as a background model which may be evaluated at any location in the image to predict the local sky value.

Several factors must be considered when choosing the image subsection size, as this will determine spatial scales the background model is sensitive to. If the image subsections are too small, they will include too few surviving unmasked pixels from which to calculate a mean. Additionally, if the subsections are smaller than the spatial scales of galaxy profiles, extended features in galaxy profiles risk being subtracted out. However, all other variations in the sky must happen at spatial scales lower than that of the sky model if they are to be fitted and removed. Because we have restricted our focus to the general case of an extragalactic sky superposed on a slowly varying foreground sky level for this work, we have chosen image subsection sizes that yield a background model sensitive to spatial variations on the order of 40 square arcseconds. Investigations in varying the subsection size showed the model is insensitive to varying the subsection size between 25 square arcseconds and 1 square arcminute. Using subsections smaller than 25 square arcseconds created subsections which were completely masked and had no usable pixels from which to estimate statistics, and subsections larger than 1 square arcminute are sufficiently large to avoid over subtracting extended tails of the galaxies simulated here. A summary flow chart of our methods is shown in Figure 4.

Zoom In Zoom Out Reset image size

The appropriate polynomial order for the background model must also be chosen with care, and will depend on a variety of factors. The background in data taken with real cameras on telescopes (in contrast to our idealized simulated data here) will contain contributions due to the optics, like scattered light and diffuse ghosts. Artifacts in detectors, such as tree rings, can give rise to variations across the CCD itself. The spatial scales over which these contributions occur will vary from instrument to instrument and telescope to telescope. This is also true for the astrophysical background itself, as the scales over which it varies can depend on the field of view. The model must vary on spatial frequencies at least as small as those discussed immediately above. If the model varies on scales finer than this, then the model will be susceptible to over-fitting noise, or fitting and subtracting extended features in galaxy profiles. As a result, the optimal polynomial order in our scheme will depend largely on the data being considered.

Fitting a polynomial background is a somewhat heuristic procedure when the true underlying sky model is not explicitly known, as is the case for most real observational data. As a consequence, the order of the polynomial used in the background model can be somewhat ad hoc. In Bosch et al. (2018), the authors use the same polynomial background model discussed above, and find a 6th order polynomial is well suited for modeling the background on the 4 K × 2 K CCDs used on the Subaru Hyper Suprime-Cam over the appropriate scales. Because our simulated data is meant to emulate DLS data, we directly turn to DLS data to investigate the appropriate polynomial order. We find that while such a model can suffer over subtraction in the immediate vicinity of bright field stars, it otherwise leaves galaxy profiles intact. Our idealized simulations include only galaxies, freeing us from this potential issue. We do not include a spatially varying component of the sky background in our idealized simulation, but for consistency with the methods discussed above, we employ a 6th order polynomial fit in our simulations.

In Figure 5, we show the result of this background estimator when used on our DLS depth random distribution simulation. In the blue histogram, we show the distribution of background estimations at the centroids of detected galaxies. The mean of this distribution is 1.115 ± 0.019 ADU pixel⁻¹ above the true sky value. We argue this overestimation is due to flux in the extended, unmasked outskirts of galaxy profiles, and undetected galaxies. We attempt to remove this bias using techniques described in the previous subsection. The overestimation in sky background can be partially ameliorated by recreating the background polynomial model, where we use our new masks in lieu of the SExtractor masks. Additionally, a simple mean of the unmasked pixels is used to estimate the sky in each image sub-section, instead of the 3σ clip. In the green curve in Figure 5, we show the distribution of background estimations at the centroids of detected galaxies using this technique. The mean sky value, 0.892 ± 0.013 ADU pixel⁻¹, while an improvement to the 3σ clip method, still suffers from an overestimation. As before, we can remove the persistent bias by subtracting the EBL flux in addition to dealing with previously unmasked extended galaxy profiles. In the red curve in Figure 5, we subtract the EBL flux level, estimated from the residual between the true and measured galaxy number counts, from each unmasked pixel in the random simulation. We then repeat the previous procedure to create the green distribution. The resulting distribution has a mean of −0.006 ± 0.013 ADU pixel⁻¹. As discussed in  Section 4.3, the biases from extended galaxy profiles and undetected galaxies must be dealt with in order to create accurate estimations of the sky background.

Zoom In Zoom Out Reset image size

Note, in its implementation in  Section 4.3, our sky estimator initially takes an image cutout centered on a galaxy, and grows the image cutout by 10 pixel increments in width and height as needed to ensure at least 4000 unmasked pixels reside in the image cutout. For the global background models discussed in this sub-section, however, we do not grow the 160 × 160 pixel sub-image when determining the mean and variance pixel value. This is done to ensure that an equally spaced grid is used in the polynomial fit, and that the statistics computed for each sub-image are representative of the pixels in that sub-image alone.

By combining the new sky estimation technique advocated here with a global polynomial background model, we can benefit from "the best of both worlds." The polynomial background model captures and smooths over spatial fluctuations on the ∼40 arcsecond scale. This ultimately yields a low variance in the distribution of predicted background values at the centers of galaxies. The combined model also benefits from the accuracy of our new sky estimation technique by accounting for flux contributions from undetected sources, and by excluding flux from tails of galaxy profiles from detected sources by growing our detection masks.

4.5. Comparing Different Sky Estimators

We compare sky background estimates for detected galaxy images in our DLS-depth simulated images. In the uniform distribution simulation, all galaxies are detected and cataloged; consequently, we use our estimators on all galaxies in this simulation. In the random distribution simulation, however, we only consider galaxies which SExtractor detected and for which GALFIT is able to converge on a sky estimation. The numbers of galaxies that meet the criteria above are 10,000 and 98,149 for the uniform and random simulations, respectively. In Figure 6, we compare sky estimations by SExtractor, GALFIT, Gaussian fit + new masks (with and without EBL correction), and polynomial background model with EBL correction and new masks. Note that the true sky background is 3240 ADU pixel⁻¹ and this pedestal has been subtracted in the uniform and random distribution simulations, although the shot noise from this sky level is included.

Zoom In Zoom Out Reset image size

The ensemble of SExtractor local background estimations has a mean of 0.567 ± 0.310 ADU pixel⁻¹ and 1.767 ± 0.198 ADU pixel⁻¹ for the uniform and random distribution simulations, respectively. As we discussed in Section 2, each image cutout in the uniform distribution simulation is completely isolated from the flux of neighboring galaxies, allowing us to test background estimation independent of blending and crowding effects. Even so, SExtractor overestimates the sky background for simulated galaxies in the uniform distribution. This is due to the flux residing in the extended galaxy profiles SExtractor fails to mask. If the measurement is done in a crowded region with a number of undetected galaxies (the random distribution simulation), the overestimation is compounded by the flux of these undetected galaxies.

GALFIT background estimates are better than SExtractor background estimates; GALFIT estimates average local backgrounds of 0.142 ± 0.331 ADU pixel⁻¹ and 0.990 ± 0.393 ADU pixel⁻¹ for the uniform and random distribution simulations, respectively. The distributions of GALFIT background estimates have noticeable extended tails, and comparatively larger scatters than other methods. This is likely due to GALFIT's sensitivity to noise. Nevertheless, the peak values in the pixel histograms as shown in Figure 6 are close to the true value in the uniform distribution simulation and the EBL in the random distribution simulation.

Estimation of sky background in our hybrid method has the highest precision of all in terms of mean and standard deviation of histogram, and is immune to the environments we have tested. Before using the polynomial fit, average local backgrounds are 0.057 ± 0.152 ADU pixel⁻¹ and −0.001 ± 0.199 ADU pixel⁻¹ for the uniform and random distribution simulations, respectively. As we discussed in Section 4.3, we can reduce the noise in sky estimation by applying the polynomial model to the sky estimates obtained from fitting a Gaussian profile to the pixel distribution in each image cutout. The resulting sky backgrounds (after EBL correction for the random distribution simulation) are 0.071 ± 0.014 ADU pixel⁻¹ and −0.006 ± 0.012 ADU pixel⁻¹ for the uniform and random distribution simulations, respectively. By doing so, we reduce the uncertainty of sky estimation by a factor of 10 or more; most sky estimates lie within ±0.0004% (0.015 ADU pixel⁻¹). Below we investigate the precision of recovery of surface brightness profiles of simulated galaxies.

In this section we focus on the effects of sky estimation accuracy on the apparent galaxy morphology. Since Edwin Hubble's first study of galaxy classification using their appearance (Hubble 1926), connections between galaxy morphology, shape, and color have provided insight into galaxy formation and evolution. Further methods for galaxy classification have been developed using the one-dimensional radially averaged profile of galaxy surface brightness (de Vaucouleurs 1948; Sérsic 1963). Among several functional forms for the profiles, the Sérsic profile is one of the most popular. The Sérsic profile is a fitting function that describes the surface brightness profile (the intensity of light as a function of distance from the center) given by

$$\begin{eqnarray}&&{\rm{\Sigma }}(R)={{\rm{\Sigma }}}_{e}\exp \left[-{b}_{n}\left\{{\left(\displaystyle \frac{R}{{R}_{e}}\right)}^{1/n}-1\right\}\right],\end{eqnarray} \tag{ 7 }$$

where Σ_e is the effective intensity (the surface brightness at the effective radius), R_e is the effective radius which encloses half of the total light (half-light radius, hereafter), b_n is the concentration which is defined so that half of the total light is inside the half-light radius for a given Sérsic index (Graham & Driver 2005), and n is the Sérsic index which describes the shape of profile and is correlated with galaxy surface brightness morphology.

As an illustrative example of the impact different background estimation techniques can have on source measurement, we measure the surface brightness profiles of model galaxies after subtracting the sky background using the three techniques discussed in Section 4: SExtractor, GALFIT, and our new method. We measure the Sérsic index in all simulations using GALFIT. As galaxy profile fits are sensitive to flux from nearby objects that are blended with the galaxy of interest, we restrict this study to the uniform distribution simulations where blending is not an issue to isolate the sky subtraction effect from the blending effect. We do so for both DLS and LSST DDF depths. To no surprise, we find that the shape of the faint outer surface brightness tail is affected by sky level mis-estimation.

The GALFIT output parameters are very sensitive to the parameters in the GALFIT start file. For this test, we use parameters derived from SExtractor. The initial parameters are determined as follows: the size of image cutout is the same as used for sky estimation (see Section 4). The total magnitude is given by MAG_AUTO; the half-light radius is given using FLUX_RADIUS with PHOT_FLUXFRAC = 0.5. The axis ratio b/a and the position angle are derived by taking 1—ELLIPTICITY and THETA_IMAGE, respectively. A Sérsic index of n = 2.5 is assumed as an initial guess in the fitting. For the PSF convolution, we generate the PSF image using GALSIM with the true parameters after matching the SExtractor's position and true position.

We investigate the effect of sky background on determining galaxy surface brightness profile types. In Figure 7 we compare the measured Sérsic index versus the true value of Sérsic index using different sky estimators for both the DLS-depth and DDF-depth uniform simulations. We compare all galaxies with which GALFIT fits do not fail nor produce problematic parameters. Galaxy surface brightness fits do not converge for fainter galaxies because their sizes become small, and the image cutout becomes noise-dominated. The numbers of galaxies with good fits are 8,428 and 9,784 for DLS and LSST DDF depths, respectively.

Zoom In Zoom Out Reset image size

Overall, our estimation of sky background performs more robustly than the other methods. We note that the scatter of our technique increases for large Sérsic index. Noise fluctuation in the faint, extended tails of large-Sérsic-index galaxies negatively impacts the fidelity of profile fitting. By contrast, the surface brightness profile of low-Sérsic-index galaxies drops rapidly, so galaxy surface brightness profile estimations are less affected by noise in these cases.

We also investigate the dependence of the Sérsic-index estimation on the Sérsic-index and magnitude of galaxies. For a given sky estimate, sky background and Sérsic index are anti-correlated: for sky estimators that tend to overestimate the background, there is a corresponding underestimate of Sérsic-index. This is more clearly seen as the brightness decreases or Sérsic-index increases. This trend is also found in previous survey data like the SDSS where sky background is overestimated (Blanton et al. 2005, 2011). There is an obvious explanation as to why there is an overestimate of Sérsic-index for high-Sérsic-index galaxies: in these cases, the true galaxy profile tends to be truncated at large radii due to the over-subtraction of sky background.

In our DLS-depth uniform simulations, our new sky background estimator results in an unbiased and less-scattered estimation of Sérsic-index. For bright galaxies (i.e., m_R ≤ 24 for our study), there is no discernible trend in over- or under-estimation of galaxy surface brightness morphology, i.e., Sérsic-index. For fainter galaxies, however, uncertainties in the fits increase, and a strong dependence on galaxy surface brightness morphology is seen. Nonetheless, our new estimator still outperforms the other sky estimators in this regime.

The performance of GALFIT runs using different sky estimators are less distinguishable at LSST DDF depth. Unlike DLS depth, the mean values of estimated Sérsic indices match the true values for all GALFIT runs. No bias is found as the magnitude becomes fainter or as the true Sérsic index increases. This is because the sky noise in LSST DDF depth is much smaller than DLS depth. Also, there is no discernible difference between our work (Gaussian fit and Gaussian fit → Poly. fit) and SExtractor. However, the GALFIT runs are still noisier than the other methods. This is due to the fact that the sky value estimated by GALFIT (μ_sky = −0.049 ± 0.225) has larger error than the other methods (μ_sky = 0.031 ± 0.022, 0.011 ± 0.011, and 0.014 ± 0.002 for SExtractor, Gaussian fit, and Gaussian fit → Poly. fit, respectively). We conclude that for the isolated galaxy case signal-to-noise ratio is a major factor in the surface brightness profile estimation. For example, see Taghizadeh-Popp et al. (2015), who show an underestimate of galaxy size near the detection limit at multiple depths. This truncation in size is closely related to the bias in profile estimate. Good signal-to-noise is a necessary but not sufficient condition for accurate profile fitting when faint undetected galaxies are included. As we found before in the DLS depth case, accounting for the undetected galaxies corrects for this sky level bias, leading to more accurate outer profile fits. However, the presence of unresolved background galaxies disturbs the Sérsic-index fit for individual galaxies even at fairly high signal-to-noise, leading to increased bias and scatter for all estimators tested compared to the uniform simulation case, where no confounding galaxies are present nearby. We will explore these effects in more detail in a future paper.

In this paper we present sky background estimation using various publicly available packages, and compare with results using our new method that grows the spatial masks around detected objects, and statistically accounts for flux from undetected faint galaxies. Our algorithm is able to recover the sky level to 4 ppm in our simulated data, an improvement over existing sky background estimation techniques. Our analysis is confined to simulations of the extragalactic sky components with added spatially uniform sky foreground; optics ghosts and scattered light around bright stars are beyond the scope of this paper.

We demonstrate that insufficiently masking the extended features of galaxies can bias sky estimation high. This occurs because flux from these galaxy regions contaminate the pixels used to estimate the sky background. While widely used estimators suffer from this bias, our technique is able to overcome it; by conservatively masking galaxies, our estimator considers pixels which are truly more representative of the true sky background. We also show that successful estimation of the underlying sky background must consider flux contributions from undetected faint background galaxies. To correct this bias, we use knowledge of galaxy counts as a number of magnitude to accurately estimate the flux contribution from these galaxies.

To demonstrate the power of our new technique, we obtain galaxy surface brightness profile fits, via the Sérsic index, using different sky estimators. Previous methods overestimate sky background, resulting in incorrect Sérsic estimates, and all show large scatter. In contrast, our two-filter estimator has the highest precision and is least affected by simulation details such as the brightness, the surface brightness profile, and the galaxy number density.

While it is beyond the scope of this paper, there are additional steps that may be taken to improve our technique, so that real imaging data may fully benefit from it. In our demonstration of this hybrid sky estimation algorithm we have used a simplified galaxy number count distribution, n(m). For a more realistic approach in applying this sky estimator to real data one should use a n(m) slope that becomes shallower beyond m_R ≃ 26 to more realistically represent the observed faint end of magnitude distribution. Additionally, we do not simulate internal reflections and scattered light in the camera, or other sources of sky variation, all of which will have to be adequately modeled for each exposure in real data. Ultimately, the level of accuracy offered by our technique is a necessary but not sufficient condition for exploration of ultra-low surface brightness. Our simulation placed galaxies randomly in the image plane; however, real galaxies are embedded in large scale structures, and cluster with each other. This may lead to a slight overdensity of undetected EBL galaxies in the vicinity of detected objects that is not reflected in our random simulation. Such an excess could still pollute the outer isophotes of the galaxy profile. The large number of undetected galaxies over a wide range of redshifts in projection strongly mitigates this effect.

We emphasize that developing detection algorithms operating at ultra low surface brightness which take full advantage of such high precision sky estimates is beyond the scope of this paper, though such sky precision would be a prerequisite. An added challenge is fitting sky across CCDs in a mosaic. For this, we note that the LSST Project has recently implemented a superior sky estimator which fits the background over an entire exposure.¹ This allows using a larger scale for the background model, including the removal of static structures (such as the average response of the camera to the sky) in the background.

The accuracy of fitting galaxy surface brightness profiles will be improved with multi-wavelength photometry since structural parameters of galaxies vary in different photometric bands, though they are correlated (Häußler et al. 2013). In such a joint fit, our sky estimator can contribute to enhanced profile accuracy, particularly in low signal-to-noise bands. Optimal detection depends on the angular scale of the object: therefore the underlying sky background must be defined on that scale and larger scales. Thus, scale dependent sky models must be developed which include all components of the apparent sky on relevant scales, from telescope optics ghosts to large scale dust.

There are many science drivers which rely on detection of low surface brightness features. Increased precision in the measurement of sky background can be applied to a better understanding of the evolution of galaxy surface brightness morphology (Conselice 2003, 2014). One can also probe the evolution of mass structure, surface brightness profile, and star formation of galaxies as a function of redshift with less bias at low surface brightness. Studies of the dark halo stellar halo connection would be less biased: there is a correlation between dark matter structure and light distribution of galaxies at late cosmic time (Wetzel & Nagai 2015; Somerville et al. 2017; Huang et al. 2017) because galaxies have different star formation histories depending on stellar mass (Qu et al. 2017) and galaxy surface brightness morphology (Wuyts et al. 2011). Finally, it is possible, and even likely, that unexpected discoveries lie at low surface brightness levels. Depending on the angular scale of the object being studied, such applications of sky estimation will require a full multi-component model of the apparent sky. The most demanding application is the unbiased detection and photometry of ultra LSB galaxies of large half-light radius. Using a noise-based non-parametric technique may be a better approach to detect such ultra faint sources. The faint-source detection capability resulting from this method has shown improved results for faint sources relative to the signal-based source detection algorithm employed by SExtractor (Akhlaghi & Ichikawa 2015). For most of these cases, a robust sky estimation with accuracy of a few parts in 10⁵ or better, and at high precision is required.

We thank Perry Gee, Lee Kelvin, Robert Lupton, Chien Peng, Brant Robertson, and Paul Price for helpful discussions. The sky sample fitting algorithm which we use is part of the sky estimator in the 2017 LSST Stack, for which we acknowledge the efforts of Steve Bickerton and Russell Owen. We thank Andrew Bradshaw and Craig Lage for comments on the manuscript. We thank the referee for suggestions that improved the manuscript. Support from DOE grant DE-SC0009999 and NSF/AURA grant N56981C is gratefully acknowledged.