© The Authors 2024

1. Introduction

The concept of surface brightness fluctuations (SBFs) was first introduced by Tonry & Schneider (1988) and Tonry et al. (1990). Since then, SBFs have been extensively studied as a powerful tool for understanding the properties of galaxies and their environments (e.g. Jensen et al. 2003). Traditionally, SBFs are used to determine extragalactic distances (e.g. Blakeslee et al. 2010; Cantiello et al. 2018). However, besides distance indicators, SBFs have shown potential to constrain stellar populations in galaxies. Stellar population analysis is generally performed by comparing the mean¹ ‘standard’ luminosity of a given population with stellar population synthesis models. In this sense, SBFs arise as a complement to obtain stellar properties (Buzzoni 1993; Worthey 1994; Raimondo et al. 2004, 2007; Cerviño 2013; Vazdekis et al. 2020; Rodríguez-Beltrán et al. 2021), among others.

Surface brightness fluctuations refer to the variation in the light across the surface of a galaxy, which arises from fluctuations in the distribution of stars among different pixels (Tonry & Schneider 1988). Surface brightness fluctuations are calculated by subtracting a mean reference image, which is the modelled surface brightness of the galaxy, correspondent with the average luminosity of the stellar population at each pixel, and then measuring the local variance of the light. From a theoretical point of view, SBFs are defined as the ratio of the variance and the mean of the luminosity distribution of individual stars. It can be demonstrated that this ratio is independent of the number of stars when a stellar population is considered (Cerviño & Luridiana 2006; Cerviño et al. 2008). In this sense, SBFs are the consequence of the pixel-to-pixel variations in the sampling of the luminosity distribution function of the stellar population, that is, pixels have a different luminosity even with a similar evolutionary status and number of stars.

Surface brightness fluctuations have been used as a tool for studying the properties of galaxies across a wide range of masses and types. However, there are several limitations to their use that should be considered. One major limitation is the dependence of SBF measurements on factors such as the quality of the image data, the signal-to-noise ratio (S/N), the point spread function (PSF), or the brightness, the size of the object and its photometry, among others (Jensen et al. 2015; Moresco et al. 2022; Cantiello & Blakeslee 2023). Additionally, the presence of background galaxies, globular clusters (GCs), foreground stars, or other sources of contamination must be masked, as they interfere with the quality of the SBF measurement.

In this context, it is necessary to establish a way of retrieving an uncertainty from the measured SBF. For instance, in Jensen et al. (2015), the uncertainty of the SBF was measured from several sources: the standard deviation of the SBF fitting, the PSF adjustments, the background variability, and the subtracted mean galaxy model. Since Jensen et al. (2015) focussed on the use of SBFs to obtain distances, they also included other additional sources of uncertainty related to the colour and distance calibrations. The total statistical uncertainty given in Jensen et al. (2015) reaches 0.1 magnitudes. Among other methods, previous authors have also provided SBF uncertainties associated with the stochasticity in the residual signal from unmasked sources, an estimate through Monte Carlo simulations while slightly varying the galaxy mask and the fitted frequencies, applying bootstrap resampling, considering the PSF mismatch, or a combination of the above (Blakeslee et al. 2001; Cohen et al. 2018; Carlsten et al. 2019). To the best of our knowledge, most authors employ the deviation of the fitting to measure the uncertainty. Besides the uncertainties inherent to the observation, such as an unknown PSF or sky background, we aim to study another source of uncertainty: the variability of a system that is intrinsically stochastic. This is how the same SBF value could coincide with different pixel distributions of the light. In the current paper, we investigate the accuracy of the SBF measurement under controlled conditions and propose a way of evaluating the precision of the observations.

Among the applications of the SBF, multiple authors have obtained SBFs from dwarf or diffuse galaxies in recent years, for instance Kim & Lee (2021), Greco et al. (2021), and Jerjen et al. (1998, 2000, 2004). As low-mass systems, dwarf galaxies are thought to be the building blocks of larger galaxies (Grebel 2001; Tosi 2003). Understanding their properties is crucial for constraining models of galaxy formation and evolution. On account of this, SBFs offer a powerful tool for studying such objects, as they are sensitive to variations in the underlying stellar population. Additionally, SBFs have been used to measure the distances of dwarf galaxies, which are notoriously difficult to determine using other methods. Nevertheless, sources with a low number of stars per pixel might not be a reliable representation of the SBF stellar population (Cerviño & Luridiana 2006; Cerviño et al. 2008; Cerviño 2013). In this work we discuss the limitations of using SBFs to study faint sources and dwarf galaxies.

Aside from studying the applicability of SBFs on dwarf galaxies, we are interested in evaluating the detection of SBF gradients in massive galaxies. Different formation histories can produce very different radial gradients in galaxies (Sánchez-Blázquez et al. 2007; Martín-Navarro et al. 2018). The spatial distribution of SBFs across a galaxy can provide additional information about galaxy structure and evolution. Specifically, the gradient of SBFs changing as a function of the radius can reveal important information about the underlying stellar population, the presence of a substructure, and the history of galaxy interactions and mergers (Cantiello et al. 2007). As was put forward in Rodríguez-Beltrán et al. (2021), a combination of mean and SBF colours is able to constrain composite stellar populations and, so, predict galaxy formation models. Several authors have studied SBF gradients by applying annular masks over different regions of the galaxy, and they have been able to retrieve the SBF gradient; this includes Cantiello et al. (2005, 2007), Sodemann & Thomsen (1995a,b) and Jensen et al. (2015). However, those authors addressed the precariousness of the observations or did not conclusively detect SBF gradients (Jensen et al. 1996). Hence, investigating the possibility and limitations of the SBF gradient measurement is a necessary task.

Having presented the state of the art on the topic, we introduce the main goals of this work: (1) To analyse the SBF derivation methodology, addressing the approximations taken and proposing improvements for its estimation; (2) To provide a measure of the uncertainty of SBF estimations due to the variability of a system that is intrinsically stochastic; (3) To evaluate the reliability of SBF retrieval under a wide range of conditions (varying parameters such as the brightness, the fluctuation, the mask applied, the PSF, the sky background, the size of the galaxy, etc.); (4) To analyse the possibility of SBF gradient detection and, if present, to address the influence when measuring the whole galaxy; and (5) to propose a recommended procedure to retrieve the SBF. We used mock galaxies as an ideal laboratory in which to perform such experiments, free from the inherent challenges associated with actual observations.

This paper is organised as follows. In Sect. 2 we present the galaxy data we use as a reference, the modelling of our mock galaxy, the SBF derivation, and its consequent fitting. In Sect. 3 we display the results of retrieving the SBF while varying a wide range of parameters. In Sect. 4 we address the SBF derivation procedure, and we warn readers about SBF biased measurements due to low flux level pixels (as in dwarf galaxies) or other sources of offset and the possibility of tracing SBF gradients. In Sect. 5 we summarise the contents of this work. In Sect. 6 we provide our computational implementation in the form of a straightforward code that estimates the uncertainty of the measured SBF. As a set of conclusions, in Sect. 7 we give recommendations and ideal conditions for a proper SBF retrieval. Finally, in Appendix A we give a more detailed description of the SBF derivation and in Appendix B we summarise the notation used in the paper.

2. Methodology

2.1. Reference galaxy data

In order to study the reliability when obtaining SBFs, we have created a synthetic galaxy image using NGC 4649 (VCC 1978) as reference. The apparent i-band magnitude of the galaxy is m_i = 9.493 ± 0.001 mag (AB) from Sloan Digital Survey (Abazajian et al. 2009), its effective radius is R_eff = 82 arcsec (van der Marel 1991) and its velocity dispersion is σ = 330.5 ± 4.6 km s⁻¹ (Davies et al. 1987). The rest of data used to define the galaxy are derived from the work of Cantiello et al. (2018) and CFHT/MegaCam imaging data from the NGVS survey (Ferrarese et al. 2012). The MegaCam general specifications state the plate scale at the centre of the field is 0.187 arcsec/pixel, so the effective radius is R_eff ≈ 438 pixels. According to Cantiello et al. (2018) the SBF magnitude of the galaxy is $\overline{m_i}=30.64\pm 0.07$ mag (AB) and its distance, derived from the SBF magnitude, is D = 16.7 ± 0.6 Mpc. From now on, we abandon the i-band notation, which is only taken as an initial reference value. The exposure time can be retrieved from Ferrarese et al. (2012) as a summed stacking of five exposures of 411 s each, leading to a total of t_exp = 2055 s. We obtain the sky background from the header of the stacked NGVS+3+0.I2 image², where NGC 4649 is located. It is specified that the minimum and the maximum sky counts found among the 5 exposures are 1725 and 1938 counts, respectively. Therefore, we take an average value for the sky background of N_Sky = 5 × 1831 = 9155 counts associated with each pixel, when all the images are added.

Aside from the observational values we have gathered, we need other additional parameters to model a synthetic galaxy, such as a Sersic index and a PSF. The observed values of the Sersic index found in the literature for NGC 4649 range approximately from n = 3.4 ± 0.5 (Vika et al. 2013) to $n=5.{36}_{-0.32}^{+0.38}$ (Kormendy et al. 2009). So, we choose an intermediate Sersic index value of n = 4, typical of an elliptical galaxy (De Vaucouleurs 1953; Caon et al. 1993). We show in Sect. 3.3 that the selected Sersic index does not drastically affect the SBF measurement. On the other hand, we generate a Gaussian PSF with a standard deviation such that 3 × σ_PSF = 4 px (σ_PSF = 1.33 px), centred in a square frame of size 2 × 3 × σ_PSF + 1 px. In comparison, the full width at half maximum (FWHM) of the observed galaxy is FWHM ≈ 0.″55, according to Cantiello et al. (2018). With the MegaCam pixel scale, this corresponds to FWHM ≈ 2.94 px. If we assume a Gaussian PSF for the NGVS observation, we find that σ_PSF ≈ FWHM/2.355 ≈ 1.24 px. This is similar to our assumption, as 3 × σ_PSF = 3 × 1.24 = 3.72 ≈ 4 px. Finally, the galaxy is centred in a squared n_pix  ×  n_pix image, where n_pix = 1605 pixel wide.

Using the total magnitude we calculate the number of counts in every pixel of the galaxy. In order to do so, first we integrate the total light of the galaxy as

$$\begin{array}{c}\hfill L=2\pi {\displaystyle {\int }_0^{\infty }}rI(r)\mathrm{d}r,\end{array}$$(1)

where we describe the light profile of the galaxy I(r) as a Sersic profile:

$$\begin{array}{c}\hfill I(r)=I_{\mathrm{eff}}{e}^{-b_n[{\left(\frac{r}{R_{\mathrm{eff}}}\right)}^{1/n}-1]}.\end{array}$$(2)

Here, b_n ≈ 2n − 1/3 (Ciotti & Bertin 1999) and I_eff is the intensity per unit area at the effective radius. Then, integrating Eq. (1) returns:

$$\begin{array}{c}\hfill L=I_{\mathrm{eff}}2\pi n{e}^{b_n}{b}_n^{-2n}{R}_{\mathrm{eff}}^2\mathrm{\Gamma }(2n),\end{array}$$(3)

with Γ being the mathematical gamma function. If we apply the negative logarithm multiplied by 2.5 we find the enclosed magnitude profile as:

$$\begin{array}{c}\hfill m={\mu }_{\mathrm{eff}}-2.5{log}_{10}(R_{\mathrm{eff}}^2)-2.5{log}_{10}\left(2\pi n{e}^{b_n}{b}_n^{-2n}\mathrm{\Gamma }(2n)\right),\end{array}$$(4)

where μ_eff is the surface brightness at the effective radius. Using our reference data, we obtain μ_eff = 22.45 mag/arcsec² or 26.09 mag per pixel.

Finally, we transform the surface brightness at the effective radius (μ_eff) to counts, N_Reff, according to MegaCam general specifications³:

$$\begin{array}{c}\hfill N_{\mathrm{Reff}}=t_{\exp }{10}^{({\mu }_{\mathrm{eff}}-\mathrm{P}\mathrm{H}\mathrm{O}\mathrm{T}\mathrm{\_}\mathrm{C}0)/-2.5},\end{array}$$(5)

where t_exp is the exposure time and PHOT_C0 = 25.743 is the nominal camera zero point defined by ELIXIR-LSB software used in the image reduction process. Then, according to Eq. (5) we get N_Reff = 1494.96 counts at the effective radius. Even though the count number returned is not an integer, we note that this value does not represent individual counts from the galaxy, but an estimation obtained from the magnitude⁴.

Subsequently, the number of counts in each pixel of the galaxy is then obtained from the Sersic profile presented in Eq. (2). The steepness of the Sersic profile might overestimate the number of counts in the innermost region of the galaxy, so the experiments performed in this work do not consider pixels in the centre.

On the other hand, the count number associated with the SBF, $\overline{N}$, is obtained similarly to that of N_Reff, in order to keep a coherent procedure. We introduce $\overline{m_i}$ in Eq. (5) as:

$$\begin{array}{c}\hfill \overline{N}=t_{\exp }{10}^{(\overline{m}-\mathrm{P}\mathrm{H}\mathrm{O}\mathrm{T}\mathrm{\_}\mathrm{C}0)/-2.5}.\end{array}$$(6)

We find an SBF value of ${\overline{N}}_{\mathrm{ref}}=22.59$ counts associated with each pixel. Here, we want to address a detail of the nomenclature applied throughout this document. $\overline{N}$ is a general way of addressing the number of counts associated with the SBF. If $\overline{N}$ has a subscript the connotations are different: ‘input’, refers to SBF values used for building a mock galaxy, it can be replaced either for ‘real’ (if the value is actually known, as for example, in a simulation) or for ‘obs’ (if the value is measured from an observation). Subscript ‘ref’ alludes to the reference value of 22.59 counts. Subscript ‘fit’ is used if the value is the result from the fitting.

The parameters presented in this section (shown in Table 1) are the values chosen for our reference image in most of the experiments of this work, unless otherwise stated.

2.2. Mock galaxy creation

Having specified the effective radius (R_eff), the image size (n_pix × n_pix), the number of counts at the effective radius (N_Reff), the number of counts associated with the fluctuation ($\overline{N}$) in each pixel, the sky background counts (N_Sky), the Sersic index (n) and the PSF model (σ_PSF), consequently, the synthetic galaxy is created as follows:

1. We create a two-dimensional Sersic image (Gal_mean(x, y)) enclosed in a square of n_pix × n_pix size, assuming that the galaxy can be well described by a Sersic model (Eq. (2)), with index (n), an effective radius (R_eff) in pixels and an intensity per unit area expressed in counts (N_Reff). So, the count value for each pixel of the mean galaxy model is denoted as N(Gal_mean(x, y)). To simplify notation, from now on we do not write the (x, y) dependence of the images. Additionally, in the notation used throughout this work we define the number of counts ‘N’ of any magnitude ‘𝒳’ as N(𝒳) or N_𝒳. The magnitude ‘𝒳’ can refer to an image or a location within an image. For instance, N(Gal_mean) is the number of counts of the mean model image (Gal_mean) for each (x, y) pixel, N_Reff is the number of counts at the effective radius or N_sky is the number of counts of the sky (which is the same for every pixel).

2. We replicate the fluctuation of the stellar population luminosity of every pixel as a random Gaussian distribution (Rand_Gauss)⁵ with mean N(Gal_mean) and variance⁶ ${\sigma }_{\mathrm{fluc}}^2=\overline{N}·N({\mathrm{Gal}}_{\mathrm{mean}})$, according to each pixel count value. This step returns an image of the modelled galaxy with its fluctuations Gal_mean, fluc. Also, we can describe this fluctuation as the addition of a random Gaussian distribution (Gal_fluc) with mean value in zero and the same variance⁷ to the mean galaxy image Gal_mean:

$$\begin{array}{cc}& {\mathrm{Gal}}_{\mathrm{mean},\mathrm{fluc}}={\mathrm{Rand}}_{\mathrm{Gauss}}[{\mathrm{Gal}}_{\mathrm{mean}},{\sigma }_{\mathrm{fluc}}^2]\hfill \\ \hfill & ={\mathrm{Gal}}_{\mathrm{mean}}+{\mathrm{Rand}}_{\mathrm{Gauss}}[0,{\sigma }_{\mathrm{fluc}}^2]={\mathrm{Gal}}_{\mathrm{mean}}+{\mathrm{Gal}}_{\mathrm{fluc}}.\hfill \end{array}$$(7)

At this point, the presence of globular clusters or background sources could be added, although we will not consider them in our experiments. For instance, GCs are influenced by a combination of factors, including galaxy properties (e.g. brighter galaxies tending to host more GCs), the observational setup or factors such as the PSF of the image. We acknowledge that we vary parameters such as the galaxy brightness or the PSF in the current work (Sect. 3), so the fluctuation contribution due to GCs may be relevant. However, in our analysis, we neglect the GCs impact, as our goal is to determine when the SBF fitting is reliable under ideal conditions. If the SBF retrieval is not trustworthy without the GCs and the background sources, it certainly will not be with them.

3. We sum a flat image of size n_pix × n_pix with a value of N_Sky in every pixel as the sky background (Sky),

$$\begin{array}{c}\hfill {\mathrm{Gal}}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n},\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{c},\mathrm{s}\mathrm{k}\mathrm{y}}={\mathrm{Gal}}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n},\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{c}}+\mathrm{Sky}.\end{array}$$(8)

4. We create a PSF as a two-dimensional Gaussian⁸ with a standard deviation of σ_PSF, centred in a square frame of size 2 × 3σ_PSF + 1 px. Then, Gal_{mean, fluc, sky} is convolved⁹ with the PSF:

$$\begin{array}{c}\hfill {\mathrm{Gal}}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n},\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{c},\mathrm{s}\mathrm{k}\mathrm{y},\mathrm{P}\mathrm{S}\mathrm{F}}={\mathrm{Gal}}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n},\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{c},\mathrm{s}\mathrm{k}\mathrm{y}}\otimes \mathrm{PSF}.\end{array}$$(9)

5. Finally, we imitate the instrumental noise. To do so, we vary every pixel count number with a random Poisson distribution (Rand_Poi)¹⁰ centred around each pixel value, that is, N(Gal_{mean, fluc, sky, PSF}). Thus, our mock galaxy is built as:

$$\begin{array}{c}\hfill {\mathrm{Gal}}_{\mathrm{mock}}={\mathrm{Rand}}_{\mathrm{Poi}}\left[{\mathrm{Gal}}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n},\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{c},\mathrm{s}\mathrm{k}\mathrm{y},\mathrm{P}\mathrm{S}\mathrm{F}}\right].\end{array}$$(10)

Moreover, we can express the addition of instrumental noise as Gal_mock = Gal_{mean, fluc, sky, PSF} + R. Here, R depicts the instrumental noise, which would be calculated separately as:

$$\begin{array}{c}\hfill R={\mathrm{Rand}}_{\mathrm{Poi}}\left[{\mathrm{Gal}}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n},\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{c},\mathrm{s}\mathrm{k}\mathrm{y},\mathrm{P}\mathrm{S}\mathrm{F}}\right]-{\mathrm{Gal}}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n},\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{c},\mathrm{s}\mathrm{k}\mathrm{y},\mathrm{P}\mathrm{S}\mathrm{F}}.\end{array}$$(11)

In this work we make use of Poisson noise to represent any source of variance that is not convolved with the PSF and we call it generically ‘instrumental noise’.

In brief, the unfolded expression for the mock galaxy is:

$$\begin{array}{c}\hfill {\mathrm{Gal}}_{\mathrm{mock}}={\mathrm{Rand}}_{\mathrm{Poi}}[({\mathrm{Rand}}_{\mathrm{Gauss}}[{\mathrm{Gal}}_{\mathrm{mean}},{\sigma }_{\mathrm{fluc}}^2]+\mathrm{Sky})\otimes \mathrm{PSF}].\end{array}$$(12)

As an example, in the left panel of Fig. 1 we show the image of a mock galaxy (Gal_mock), where the input data are taken from Sect. 2.1. In the right panel of the same figure, we show its associated radial profiles for the mean model (Gal_mean), the final synthetic galaxy (Gal_mock) and the background sky with the instrumental noise (Sky + R).

Fig. 1. Representation of our models computed with the data presented in Sect. 2.1. Left: image of the mock galaxy model (Galmock). Right: radial profiles of the mean model (Galmean) as a blue solid line, the mock model (Galmock) as an orange dotted line, and the sky background with the instrumental noise as a solid green line. The effective radius Reff is shown as a red circumference (left) and a vertical line (right).

2.3. SBF derivation

For measuring the SBF amplitude in the mock galaxy we start by rewriting Eq. (12) as:

$$\begin{array}{c}\hfill {\mathrm{Gal}}_{\mathrm{mock}}=({\mathrm{Gal}}_{\mathrm{mean}}+{\mathrm{Gal}}_{\mathrm{fluc}}+\mathrm{Sky})\otimes \mathrm{PSF}+R.\end{array}$$(13)

Again, Gal_fluc represents the Gaussian fluctuation around the mean value (Gal_mean) due to the stellar population luminosity variation, as introduced in Eq. (7). The term R is the instrumental (Poisson) noise introduced in Eqs. (10) and (11). In this case, the PSF, which is contained in a square of side 2 × 3σ_PSF + 1 px, is now re-inscribed in the centre of a blank-template of the same size as the rest of images. This is done by convolving the PSF with a n_pix × n_pix image of zero values for all the pixels except the central one, which is assigned a value of one. Also, we note that Sky ⊗ PSF = Sky, since Sky is a constant image. We emphasise that the experiments of this study are conducted under ideal conditions, and we treat the sky as known and flat. In observations, however, the sky value has its own uncertainty and is not necessarily spatially or temporally invariant.

In order to arrive to the SBF term, first we subtract the sky background and the mean model (smoothed with the PSF) to the observed galaxy:

$$\begin{array}{c}\hfill {\mathrm{Gal}}_{\mathrm{mock}}-\mathrm{Sky}-{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}={\mathrm{Gal}}_{\mathrm{fluc}}\otimes \mathrm{PSF}+R.\end{array}$$(14)

Next, we normalise by the square root of the mean model convolved with the PSF. This step is necessary to arrive at the SBF definition, that is, the stellar luminosity distribution variance divided by its mean (it will be squared later on, in Eq. (16)).

$$\begin{array}{c}\hfill \frac{{\mathrm{Gal}}_{\mathrm{mock}}-\mathrm{Sky}-{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}=\frac{{\mathrm{Gal}}_{\mathrm{fluc}}\otimes \mathrm{PSF}+\mathrm{R}}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}.\end{array}$$(15)

We denote as Gal_mock fluc the left side of Eq. (15). The right-hand side of the equation contains two noise components with a null mean value: the first one is the population fluctuation (the SBF), which is convolved with the PSF; the second one is the instrumental noise, both varying from pixel to pixel. The fluctuation contribution Gal_fluc and the PSF can be disentangled in the Fourier space. We do so applying the power spectrum to Eq. (15), this is PS(f) = |𝔉(f)|² = 𝔉(f)⋅𝔉(f)^†, where † denotes the complex conjugate.

$$\begin{array}{c}\hfill PS({\mathrm{Gal}}_{\mathrm{m}\mathrm{o}\mathrm{c}\mathrm{k}\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{c}})=PS\left(\frac{{\mathrm{Gal}}_{\mathrm{fluc}}\otimes \mathrm{PSF}+\mathrm{R}}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right).\end{array}$$(16)

As the power spectrum involves a squared Fourier transform, the summed terms of Eq. (16) are developed as the square of complex numbers¹¹ (PS(f + g) = |𝔉(f)+𝔉(g)|² = |𝔉(f)|² + 𝔉(f)^† ⋅ 𝔉(g)+𝔉(f)⋅𝔉(g)^† + |𝔉(g)|² = PS(f)+𝔉(f)^† ⋅ 𝔉(g)+𝔉(f)⋅𝔉(g)^† + PS(g)). Then, by applying the convolution theorem, 𝔉(f ⊗ g) = 𝔉(f)⋅𝔉(g), we separate Gal_fluc and the PSF contribution.

$$\begin{array}{cc}& PS({\mathrm{Gal}}_{\mathrm{m}\mathrm{o}\mathrm{c}\mathrm{k}\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{c}})=PS\left(\frac{{\mathrm{Gal}}_{\mathrm{fluc}}}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right)·PS(\mathrm{PSF})\hfill \\ \hfill & +PS\left(\frac{R}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right)\hfill \\ \hfill & +\mathfrak{F}{\left(\frac{{\mathrm{Gal}}_{\mathrm{fluc}}}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right)}^{†}·\mathfrak{F}{\left(\mathrm{PSF}\right)}^{†}·\mathfrak{F}\left(\frac{R}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right)\hfill \\ \hfill & +\mathfrak{F}\left(\frac{{\mathrm{Gal}}_{\mathrm{fluc}}}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right)·\mathfrak{F}\left(\mathrm{PSF}\right)·\mathfrak{F}{\left(\frac{R}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right)}^{†}.\hfill \end{array}$$(17)

Finally, we apply an azimuthal average to these power spectra, which in the following we denote with the subindex ‘r’ to indicate a radial profile. Applying an azimuthal average reduces the dimension of the 2D Fourier space images (n_pix × n_pix) of Eq. (17) into 1D radial profiles, each with length n_pix/2, dependent on the frequency (k) in px⁻¹ units. Note that, as far as only sums and scalar multiplications are involved, we can consider the total azimuthal average as a sum of the azimuthal averages of the different components. Thus, the SBF value to recover ($\overline{N}$) appears as a constant value corresponding to the $PS{({\mathrm{Gal}}_{\mathrm{fluc}}/\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}})}_{\mathrm{r}}$ term:

$$\begin{array}{cc}& PS{({\mathrm{Gal}}_{\mathrm{m}\mathrm{o}\mathrm{c}\mathrm{k}\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{c}})}_{\mathrm{r}}=\overline{N}·PS{(\mathrm{PSF})}_{\mathrm{r}}+PS{\left(\frac{R}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right)}_{\mathrm{r}}\hfill \\ \hfill & +\sqrt{{\overline{N}}^{†}}·{(\mathfrak{F}{\left(\mathrm{PSF}\right)}^{†}·\mathfrak{F}\left(\frac{R}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right))}_{\mathrm{r}}\hfill \\ \hfill & +\sqrt{\overline{N}}·{(\mathfrak{F}\left(\mathrm{PSF}\right)·\mathfrak{F}{\left(\frac{R}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right)}^{†})}_{\mathrm{r}}.\hfill \end{array}$$(18)

The two fluctuation terms considered here, in essence, the stellar population luminosity variability ($PS{({\mathrm{Gal}}_{\mathrm{fluc}}/\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}})}_{\mathrm{r}}$) and the instrumental noise ($PS{(R/\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}})}_{\mathrm{r}}$), are constant on average with respect to the frequency; as the example of Fig. 2 shows, using the data of our reference galaxy. In this figure, both radial profiles are flat and, therefore, neither of their associated images have any structure. Thus, the fluctuation term ($PS{({\mathrm{Gal}}_{\mathrm{fluc}}/\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}})}_{\mathrm{r}}$) can be represented by a constant value, which in this case matches the input SBF (${\overline{N}}_{\mathrm{ref}}$).

Fig. 2. Azimuthally averaged power spectrum of both fluctuation contributions: the stellar population luminosity variability ($PS{({\mathrm{Gal}}_{\mathrm{fluc}}/\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}})}_{\mathrm{r}}$) is in black and the instrumental noise ($PS{(R_{\mathrm{approx}}/\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}})}_{\mathrm{r}}$) is in grey. In this case we consider a mask, introduced as in Eq. (21), with the same data used in Fig. 4. Note that the input SBF value (${\overline{N}}_{\mathrm{ref}}=22.59$ counts) coincides with the average value of the fluctuation term.

At this point, most of the literature assumes neglecting the crossed term that appears due to the square modulus of the power spectrum. In Fig. 3 (again, using the data of our reference galaxy) we demonstrate that the power spectrum of the imaginary component is null and the real component is three orders of magnitude lower than the rest of the terms. Additionally, numerical calculations show that the influence of the real part of the crossed term is negligible (see Sect. 4.1, Table 2). Consequently, Eq. (18) is rewritten and fitted as:

$$\begin{array}{c}\hfill PS{({\mathrm{Gal}}_{\mathrm{m}\mathrm{o}\mathrm{c}\mathrm{k}\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{c}})}_{\mathrm{r}}\approx \overline{N}·PS{(\mathrm{PSF})}_{\mathrm{r}}+PS{\left(\frac{R}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right)}_{\mathrm{r}}.\end{array}$$(19)

Fig. 3. Real and imaginary components of the crossed term presented in Eq. (18). Using the same data as in Fig. 4.

In Appendix A we show this mathematical development with more detail, also including the presence of a mask.

2.4. Applying a mask

Commonly, during the SBF measurement, masks are applied in order to cover areas of the image that could hinder obtaining the SBF. For instance, in the centre of the galaxy either the number of counts could saturate the detector or the profile obtained from a theoretical model could overestimate the light, as in our case. Moreover, the mean model (Gal_mean) subtracted from the galaxy image is not representative in the pixels with the highest count value (this corresponds to the central pixels in our work), as the mean will always return lower values than the maximum. Thus, in this work all the experiments with mock masked galaxies are performed, at least, 4 pixels away from the centre. On the other hand, in the external regions of the galaxy there might be not enough light to obtain the SBF, either because the observed galaxy is too faint or because the theoretical model undervalues the profile. Additionally, there might be other saturated or intrusive elements that should be hidden, such as foreground stars, globular clusters, other galaxies, cosmic rays, etc.

These masks are applied by multiplying zero-one (False-True) n_pix × n_pix images to our synthetic galaxy image (Gal_mock ⋅ Mask = Gal_mock mask). In so doing, Eq. (13) would look:

$$\begin{array}{c}\hfill {\mathrm{Gal}}_{\mathrm{mock}}·\mathrm{Mask}=(({\mathrm{Gal}}_{\mathrm{mean}}+{\mathrm{Gal}}_{\mathrm{fluc}}+\mathrm{Sky})\otimes \mathrm{PSF}+R)·\mathrm{Mask}.\end{array}$$(20)

Applying a similar procedure as in Sect. 2.3, we find that Eq. (19) with a mask is¹²:

$$\begin{array}{cc}& PS{({\mathrm{Gal}}_{\mathrm{m}\mathrm{o}\mathrm{c}\mathrm{k}\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{c}\mathrm{m}\mathrm{a}\mathrm{s}\mathrm{k}})}_{\mathrm{r}}\approx \hfill \\ \hfill & \overline{N}·{(PS(\mathrm{PSF})\otimes PS(\mathrm{Mask}))}_{\mathrm{r}}+PS{\left(\frac{R·\mathrm{Mask}}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right)}_{\mathrm{r}}.\hfill \end{array}$$(21)

This expression is commonly found in the literature as P(k) = P₀ ⋅ E(k)+P₁ (e.g. Tonry et al. 1990). Here, P(k) corresponds to the fluctuation frame term, P₀ is the average flux from the fluctuations to fit and E(k) is the expectation power spectrum (this is, the convolution of the PSF power spectrum and the mask power spectrum). P₁ is the constant instrumental noise component (or any other sources of variance that are not convolved with the PSF).

The derivation of the power spectrum is mathematically not completely rigorous, for example, by neglecting the crossed term or altering the operational order of the fluctuation, the PSF and the mask. Nonetheless, the effect of these approximations on the power spectrum is minimal (Jensen et al. 1998), as demonstrated in Sect. 4.1 and in the Appendix A.

2.5. Measuring SBFs from images

Our goal consists in fitting the fluctuation value from Eq. (21) and comparing the measured ${\overline{N}}_{\mathrm{fit}}$ with the known input value ${\overline{N}}_{\mathrm{real}}$. In the procedure of deriving an SBF from an observation there are some considerations to take into account: on the one hand, for calculating Gal_{mock fluc mask} it is necessary to use a model for Gal_mean, which serves as a reference to obtain the fluctuations. When using real observations, such a mean model is obtained by smoothing the image, by an isophote fitting of the galaxy, with Sersic profiles or other methods (e.g. Pahre et al. 1999; Cantiello et al. 2018; Carlsten et al. 2019). The sky value is often obtained from an empty region of the observed image and the PSF can be derived from the profile of isolated stars, always considering the particularities of each observation. However, for our experiments, we already know the mean model behind our mock galaxy, the sky count number and the shape of the PSF. So, we can subtract them directly and recover an accurate version of the fluctuation image. Here, we reiterate that we do not account for fluctuations coming from other sources than the galaxy and the instrumental noise. That is, in our ideal image we neglect globular clusters, background galaxies, foreground stars, etc.

2.5.1. Modelling the instrumental noise

In our work, we already know the non-correlated noise added to the image (R), as shown in Eq. (11). However, for each fitting we assume as unknown any fluctuation present on the image, that is, neither the fluctuation of the stellar population luminosity nor the instrumental noise. Therefore, R needs to be calculated differently, so, we propose modelling the instrumental noise with a slight difference with respect to Eq. (11), this is, without the term Gal_fluc:

$$\begin{array}{c}\hfill R_{\mathrm{approx}}={\mathrm{Rand}}_{\mathrm{Poi}}[{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}+\mathrm{Sky}]-({\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}+\mathrm{Sky}).\end{array}$$(22)

With this equation we derive a map for R_approx(x, y). This approach could be used both for mock galaxies and for some observations, as far as the mean galaxy, PSF and sky background are known. We note that this work considers Poisson noise, but this approach should be adapted to the nuances of each observation.

Modelling independently the instrumental noise leaves Eq. (21) with only the SBF left to fit. This reduces the uncertainty, as we demonstrate in Sect. 4.1 (see Table 2). This approximation is reliable if the Gal_fluc contribution is small enough compared to the rest of the terms. In the current work this is backed up by the third criterion presented in Eq. (25), in which we require the contribution of galaxy counts to be 10 times larger than the SBF counts.

2.5.2. SBF fitting example

In this section we create mock galaxy images using Eq. (12) based on the data from Sect. 2.1. Then we attempt to recover its SBF value following Eq. (21), where we already know Gal_mean, the PSF, the Sky and R_approx images. We fit¹³ the fluctuation value (${\overline{N}}_{\mathrm{fit}}$) and we compare it with the actual input value (${\overline{N}}_{\mathrm{real}}$) of Gal_mock.

As an example, we take the reference galaxy created in Fig. 1 (with size of n_pix ≈ 1605 px) and attempt to fit its SBF value imitating the mask used in Cantiello et al. (2018). The results are shown in Fig. 4. In the inset panel we show the resulting image of calculating Gal_{mock fluc mask} with a centred annular mask of radii r₁ ≈ 80 px (with N(Gal_mean(r₁)) ≈ 21 132 counts) and r₂ ≈ 332 px (with N(Gal_mean(r₂)) ≈ 2492 counts). We choose these values based on Fig. 1 and Table 2 (⟨Rad⟩ column) of Cantiello et al. (2018), with an inner radius of r₁ = 15 arcsec and an external radius of r₂ = r₁ + ⟨Rad⟩ = 15 + 47.2 arcsec. As our work considers an ideal laboratory for the mock galaxies, no contaminants are present, so we only assume an annular mask. In a real observation, any other light source needs to be taken into account and covered, the Mask term must be a combination of all of these contributions.

Fig. 4. Rescaled (multiplied by PS(PSF)r/(PS(PSF)⊗PS(Mask))r) logarithm of the radial power spectrum profiles: the masked instrumental noise term ($PS{(R_{\mathrm{approx}}·\mathrm{Mask}/\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}})}_{\mathrm{r}}=PS{(R_{\mathrm{approx}}^{\mathrm{norm}})}_{\mathrm{r}}$) is the cyan line; the observed masked fluctuation (PS(Galmock fluc mask)r) is the blue line; the right part of Eq. (21) with the real input value of the SBF (${\overline{N}}_{\mathrm{real}}$) is the red line, and with the fitted fluctuation (${\overline{N}}_{\mathrm{fit}}$) is the dashed green line. The selected range of frequencies where the fitting is performed happens between kfit, i = 75 px−1 and kfit, f = 400 px−1, marked with a pale yellow vertical region. In this example, the input SBF value is ${\overline{N}}_{\mathrm{real}}={\overline{N}}_{\mathrm{ref}}=22.59$ and the fitted SBF is ${\overline{N}}_{\mathrm{fit}}=22.47$. The embed image is the observed fluctuation Galmock fluc mask with a circular mask of r1 ≈ 80 px and r2 ≈ 332 px applied.

In the main panel of Fig. 4 we show the logarithm of the radial power spectrum obtained for the different components of Eq. (21): we show with a cyan line the power spectrum of the normalised (by $\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}{)}_{\mathrm{r}}$) and masked instrumental noise ($PS{(R_{\mathrm{approx}}^{\mathrm{norm}})}_{\mathrm{r}}=PS{(R_{\mathrm{approx}}·\mathrm{Mask}/\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}})}_{\mathrm{r}}$); with a blue line we show the power spectrum of the (‘observed’) mock galaxy fluctuation (PS(Gal_{mock fluc mask})_r); with a red solid line and a green dashed line we show the right part of Eq. (21) for the real input value of the SBF (${\overline{N}}_{\mathrm{real}}$) and for the fitted fluctuation (${\overline{N}}_{\mathrm{fit}}$), respectively. As commonly done in the literature, each one of these power spectra has been rescaled by multiplying with the ratio of the PSF loss due to the mask, that is, PS(PSF)_r/(PS(PSF)⊗PS(Mask))_r. In such way, the y-axis of different SBF figures can be compared independently of the mask used.

The selected range of frequencies where the fitting is performed (between k_fit, i = 75 px⁻¹ to k_fit, f = 400 px⁻¹) is marked with a pale-yellow vertical region. All the examples presented in this work fit the fluctuation between these frequencies (exceptions are mentioned when necessary). After numerous tests, in this range we can ensure a proper following of Gal_mock fluc shape, without entering too much into the noisy or flat, non-informative power spectrum frequency intervals (at low and high frequencies, respectively). These k_fit values are selected for an image size of n_pix = 1605 px and a given PSF of 3 × σ_PSF = 4 px, when these parameters are changed, it is necessary to adjust the interval of fitting.

Our modelled masked instrumental noise term ($PS{(R_{\mathrm{approx}}·\mathrm{Mask}/\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}})}_{\mathrm{r}}=PS{(R_{\mathrm{approx}}^{\mathrm{norm}})}_{\mathrm{r}}$, or commonly P₁ in the literature) is displayed to show the contrast between it and the mock masked fluctuation (PS(Gal_{mock fluc mask})_r, as a blue line). In reality, PS(Gal_{mock fluc mask})_r would correspond to the observed fluctuation. A significant contrast between the two is key for a reliable measurement. Then, PS(Gal_{mock fluc mask})_r is meant to be compared with the fitted fluctuation (the right part of Eq. (21), with $\overline{N}={\overline{N}}_{\mathrm{fit}}$, as a green dashed line). This would correspond to comparing an observed fluctuation (commonly P(k) in the literature) and its fitting (P₀ ⋅ E(k)), respectively. We also plot input ‘real’ introduced fluctuation (the right part of Eq. (21), with $\overline{N}={\overline{N}}_{\mathrm{real}}$, as a red line), which is meant to be compared to the green line associated with ${\overline{N}}_{\mathrm{fit}}$. Even if both are almost identical, we show them because they represent our methodology for evaluating the input against the output values (see next section).

Although not shown in Fig. 4, the PSF is responsible for the shape of the PS. If the PSF is narrower in the physical plane the shape of PS(Gal_{mock fluc mask})_r increases its frequency width. It is worth mentioning other features found when changing parameters in Eq. (21), such as the size of the image, the mask or the sky background. The number of pixels (n_pix) employed for the SBF measurement fixes the final frequency k = n_pix/2, constraining the width of the PS and, therefore, the range of frequencies where the fitting is worth. For instance, an image with a larger n_pix × n_pix presents a larger range of frequencies where the fitting can be performed. For a given mask and brightness, an image with a larger n_pix × n_pix value provides a lower power spectrum. In a different sense, applying a mask reduces the value of the power spectrum of the SBF and the uncorrelated noise; but the contrast between both remains constant. However, the larger the number of masked pixels the less information available for the measurement. And finally, the sky background is responsible for the difference between the value of $PS{(R_{\mathrm{approx}}^{\mathrm{norm}})}_{\mathrm{r}}$ and the expected value of PS(Gal_{mock fluc mask})_r at k = 0 (corresponding to the $\overline{N}$ value). A high sky count value enlarges the effect of R, making such contrast lower and increasing the noise effect when fitting.

In Fig. 4, our known input fluctuation was ${\overline{N}}_{\mathrm{real}}={\overline{N}}_{\mathrm{ref}}=22.59$ counts or ${\overline{m}}_{\mathrm{ref}}=30.64$ mag, while the fitted result is ${\overline{N}}_{\mathrm{fit}}=22.47$ counts or ${\overline{m}}_{\mathrm{fit}}=30.646$ mag. Additionally, we find clear similarities in the shape of the power spectrum when comparing the results of Fig. 4 with the results of Fig. 1 in Cantiello et al. (2018) for the galaxy selected in Sect. 2.1.

2.5.3. Reliability of the SBF measurement

We evaluate the reliability of the SBF estimate with two parameters, the relative error and the relative standard deviation of the fitting:

$$\begin{array}{cc}& {\epsilon }_{\mathrm{rel}}=\frac{|{\overline{N}}_{\mathrm{real}}-{\overline{N}}_{\mathrm{fit}}|}{{\overline{N}}_{\mathrm{real}}},\hfill \end{array}$$(23)

$$\begin{array}{cc}& {\sigma }_{\mathrm{fit}}=\frac{{\sigma }_{\mathrm{cov}}}{{\overline{N}}_{\mathrm{fit}}},\hfill \end{array}$$(24)

where σ_cov is the standard deviation in counts returned from the computational fitting¹³.

The relative standard deviation (σ_fit) estimates the quality of the least squares fitting, while the relative error (ε_rel) refers to the accuracy of the ${\overline{N}}_{\mathrm{fit}}$ result. Thus, the experiment performed in Fig. 4 returns a relative error of ε_rel = 0.53% and a relative standard deviation of σ_fit = 0.18%.

2.6. Criteria

For this work we take as a valid SBF estimation those measurements with a σ_fit and a ε_rel lower than a 10%. In addition, we only consider non-masked pixels of the galaxy (Gal_mock mask) with values larger than 10 times the input SBF in counts. We check if every pixel fulfils this criterion (stated in the third line of the equation below), otherwise the measurement is not performed. This condition assures a Gaussian probability distribution of the integrated light among the pixels, as approximated from the galaxy modelling of Tonry & Schneider (1988), Cerviño et al. (2008). For lower count values, the Gaussian condition is not assured (Cerviño & Luridiana 2006), so the traditional SBF modelling is not necessarily physically correct. This is discussed further in Sect. 4.2.

In summary, our adopted criteria are:

$$\begin{array}{c}\hfill \mathrm{Criteria}\{\begin{array}{cc}& {\epsilon }_{\mathrm{rel}}<10\%,\hfill \\ \hfill & {\sigma }_{\mathrm{fit}}<10\%,\hfill \\ \hfill & N({\mathrm{Gal}}_{\mathrm{m}\mathrm{o}\mathrm{c}\mathrm{k}\mathrm{m}\mathrm{a}\mathrm{s}\mathrm{k}}(x,y))-N(\mathrm{S}k{y}_{\mathrm{mask}}(x,y))\ge 10·\overline{N}.\hfill \end{array}\end{array}$$(25)

Throughout this work, we take as known the luminosity distribution in each pixel, on the contrary, in real observations the SBF value is unknown a priori. Therefore, in observations we can calculate σ_fit, but we are not able to calculate ε_rel, as ${\overline{N}}_{\mathrm{real}}$ is unknown. This is one reason why the modelling presented in this work is a very useful tool: we select an input value of ${\overline{N}}_{\mathrm{real}}$ with which we can foresee if an observation will return a reliable result. Moreover, the relative standard deviation of the fitting σ_fit is not representative of the accuracy of the returned SBF, but only the quality of the fit. Thus, for the purposes of this work ε_rel serves as a guide for obtaining a trustworthy SBF under different conditions.

In observations, other consequence of being oblivious to the real SBF value is that the third criterion cannot be guaranteed before performing the measurement. However, we suggest checking the condition $N({\mathrm{Gal}}_{\mathrm{mock}\mathrm{mask}})-N({\mathrm{Sky}}_{\mathrm{mask}})\ge 10\times {\overline{N}}_{\mathrm{obs}}$ on every non-masked pixel after the fitting, to assure the reliability of the observed SBF. We are aware that the parameters we vary are not necessarily independent (e.g. they depend on distance or stellar population properties), but we consider here that these parameters are unrelated and selected ad hoc, as a way to explore several scenarios.

Note that for the relative error and the relative standard deviation of the fitting, the sky is implicitly embedded. As we show in Sect. 3.3, the uncertainty increases for large sky background values, due to a smaller contrast between P₀ and P₁. The third condition of Eq. (25) already includes the sky contribution to the image. The latter criterion should only consider the brightness of the galaxy stellar population itself, without any other contribution.

Besides the criteria of Eq. (25), in real observations the sky itself has its own uncertainty and an offset in its estimation can distort the SBF measurement. Similarly, a bad modelling of the mean reference galaxy or the PSF lead to an offset in the SBF. These cases are quoted in more detail in Sect. 4.3.

3. Results

In order to evaluate the SBF estimation procedure we create mock galaxies for a wide range of conditions. We compare the ${\overline{N}}_{\mathrm{fit}}$ obtained from the mock galaxies with our known input ${\overline{N}}_{\mathrm{real}}$ and we study the uncertainty of the fitting, as explained in Sects. 2.5.2 and 2.5.3. In this way we can explore the parameter space looking for those galaxies where the retrieved SBF is reliable based on the criteria shown in Sect. 2.6. The galaxy parameters we choose to vary are: the magnitude of the galaxy, the magnitude of the fluctuation, the effective radius, the PSF, the Sersic index and the exposure time. For these variations of the galaxy parameters we study the fitting results when applying masks of different sizes. As in an observational setup, for each iteration the image remains unchanged when different masks are applied. We emphasise that the mean galaxy model, the sky and the PSF are assumed to be known in the fitting process. Therefore, the uncertainties addressed in this section are representative of the variability of the system due to the intrinsic stochasticity of the stellar population luminosity distribution and the instrumental noise.

With these experiments, we aim to evaluate the uncertainty of the SBF via Monte Carlo simulations. To do so, we create n_sim = 50 mock galaxies each time a parameter is varied¹⁴ except the mask. For each mock galaxy we compute ε_rel and σ_fit. From the resulting distribution of those parameters we choose the higher 90% percentile (ε_rel, 90%, σ_fit, 90%), while also satisfying the criteria of Eq. (25). In a distribution of 50 simulations, this corresponds to the 45th highest value. This provides a conservative estimate of how unfavourable the obtained SBF would be. Additionally, we calculate the relative 90% width of the ${\overline{N}}_{\mathrm{fit}}$ distribution results. This is calculated as the subtraction of the 95% and 5% percentile values, divided by the mean value of the ${\overline{N}}_{\mathrm{fit}}$ distribution $⟨{\overline{N}}_{\mathrm{fit}}⟩$. In our case those 5% and 95% percentile values are the 2^nd lowest and the 48^th highest ${\overline{N}}_{\mathrm{fit}}$ values found after sorting the 50 simulations results. Thus, Δ_90% represents the width covered by the 90% of the distribution of ${\overline{N}}_{\mathrm{fit}}$ values, giving a measurement of the precision of the fitting results:

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{90\%}=\frac{{\overline{N}}_{\mathrm{fit},95\mathrm{th}\%}-{\overline{N}}_{\mathrm{fit},5\mathrm{th}\%}}{⟨{\overline{N}}_{\mathrm{fit}}⟩}.\end{array}$$(26)

In summary, the key parameters for the analysis of this work are: ε_rel, 90%, σ_fit, 90% and Δ_90%, as introduced above. All the following figures are calculated with the previous procedure, making use of the distribution of results obtained from the mock galaxy simulations, except for those in Sect. 3.3, where further explanations are provided. In Sect. 5 we summarise the situation with a flowchart of the procedure explained in this section.

3.1. Parameter space: Masks

The calculation of SBFs is performed using masks of different sizes. Having masked point sources or bad pixels, the SBF is then commonly measured within an annulus (of properly selected width and eccentricity). This annular mask permits hiding regions of the galaxy that negatively affect the SBF measurement, as explained in Sect. 2.4. Also, in order to study the possible detection of SBF gradients we require annuli of different width and radius. Thus, we study the SBF fitting for all the varying combinations of annular rings defined by an internal and an external radius (r₁, r₂), with circular shape. We move from 4 pixels (avoiding the overestimation of the centre) to R_eff (as representative of a region of the galaxy with proper signal-to-noise), with a step of Δr = 20 px. The rest of parameters are taken as described in Sect. 2.1. We perform 50 Monte Carlo simulations of our mock galaxy, then, each one of these realisations is analysed for every (r₁, r₂) pair.

In Fig. 5 we show the colour maps of the 90% percentile values of the relative error (ε_rel, 90%) in panel (A), 90% percentile values of the relative standard deviation (σ_fit, 90%) in panel (B), the relative 90% width of the ${\overline{N}}_{\mathrm{fit}}$ distribution (Δ_90%) in panel (C) and the number of pixels within a mask (n_mask pix) in panel (D), all of them dependent on the masks size determined by (r₁, r₂). The contour lines found over the map in panels (A, B, C) correspond to the number of pixels shown in panel (D). We mark as a black circled-cross the reference galaxy presented in Sect. 2.1 with the mask applied in its source reference (r₁ ≈ 80 px, r₂ ≈ 332 px; Cantiello et al. 2018).

Fig. 5. Colour maps for the reference galaxy when applying different masks (r1, r2): the 90% percentile of the relative error in panel (A), up to εrel, 90% ≥ 10%; the 90% percentile of the relative standard deviation in panel (B), up to σfit, 90% ≤ 2%; the relative 90% width of the ${\overline{N}}_{\mathrm{fit}}$ distribution in panel (C), up to Δ90% ≤ 10%; and the number of pixels within each one of the masks (nmask pix) in panel (D). The contour lines found over the map in panels (A), (B) and (C) correspond to the number of pixels shown in panel (D). The black circled-cross signals the reference galaxy presented in Sect. 2.1 with the mask applied in its source reference Cantiello et al. (2018); r1 ≈ 80 px, r2 ≈ 332 px).

Comparing panels (A, B, C) with (D) we find a relation between the number of pixels of the mask and ε_rel, 90%, σ_fit, 90% and Δ_90%. This behaviour is especially clear for the relative error. Comparing those panels with (D) there are larger errors in regions of low numbers of pixels, that is, in very thin annuli or in small masks. For instance, any mask with 60 000 pixels or fewer will most likely produce an error higher than ε_rel, 90% > 9 %. On the other hand, panel (B) shows that the relative standard deviation of the fitting fulfils the criterion for every mask, with values lower than σ_fit, 90% < 2 %. In panel (C) we find Δ_90% < 10% for masks with a number of pixels approximately larger than 120000. Among these three parameters (panels A, B, C), we find that ε_rel, 90% is about a factor 10 larger than σ_fit, 90%, then, Δ_90% is about 2–4 times larger than ε_rel, 90%.

From the results of this section we can draw some general notions:

– First, the relative error (ε_rel, 90%), the relative standard deviation of the fitting (σ_fit, 90%) and the relative 90% width of the ${\overline{N}}_{\mathrm{fit}}$ distribution (Δ_90%) are tightly related to the number of pixels within the mask (n_mask pix).

– Second, the relative error (ε_rel, 90%) estimates the accuracy of the result and is a more restrictive constraint than the relative standard deviation (σ_fit, 90%). The relative standard deviation just measures the quality of the least squares fit, but is not representative of how close the fitted SBF (${\overline{N}}_{\mathrm{fit}}$) is to the real SBF value (${\overline{N}}_{\mathrm{real}}$).

– Third, the relative 90% width of the ${\overline{N}}_{\mathrm{fit}}$ distribution (Δ_90%) is the most restrictive parameter. Although, it does not provide information about the accuracy for finding a ${\overline{N}}_{\mathrm{fit}}$ value similar to ${\overline{N}}_{\mathrm{real}}$, it provides a much more conservative uncertainty estimate than ε_rel, 90% and σ_fit, 90%.

Since we are interested in the accuracy in the SBF estimate, henceforth we use as a proxy the number of pixels (n_mask pix) against the relative error (ε_rel, 90%) in our experiments. Moreover, studying ε_rel, 90% takes advantage of creating mock galaxies, as in this work. The other parameters, σ_fit, 90% and Δ_90%, are presented only in certain cases of interest. Additionally, we would like to propose that in real observations it is always possible to build mock galaxies as in Sect. 2.2 using the retrieved (observationally fitted) SBF, then calculate σ_fit, 90% and Δ_90% with respect to our initially fitted ${\overline{N}}_{\mathrm{obs}}$.

3.2. Parameter space: Galaxy brightness and SBF magnitude

Once our reference galaxy has been studied through different mask sizes, we study the parameter space when varying the brightness of the galaxy and its fluctuation contribution. We use as a reference the galaxy described in Sect. 2.1 and we vary both its apparent magnitude from m = 3 mag to m = 14 mag, in intervals of Δ(m) = 0.5 mag, and the SBF magnitude from $\overline{m}=25$ mag to $\overline{m}=37$ mag, with $\mathrm{\Delta }(\overline{m})=0.5$ mag. We recognise that some of these values might be unrealistic, but we keep such ranges for illustrative purposes and exploring the parameter space. As justified in the previous section, for these experiments we study the behaviour of the relative error only.

Again, we perform 50 galaxy simulations for each case and we obtain the 90% percentile of the distribution of ε_rel, 90%. In Fig. 6 we show the relative error (ε_rel, 90%) for each ($m,\overline{m}$) combination. The SBF estimate is obtained for three centred annular masks: we fix r₁ = 4 px, then we select r₂ = R_eff for panel (A), r₂ = R_eff/2 (B) and r₂ = R_eff/3 (C). Also, in Fig. 6 we mark with a black circled-cross the reference galaxy based on Sect. 2.1 data ($m=9.493,\overline{m}=30.64$).

Fig. 6. Colour map of the relative error (εrel, 90%) obtained after fitting the SBF for mock galaxies of different brightness and SBF magnitudes ($m,\overline{m}$). The fitting was performed for three centred annular masks, with fixed r1 = 4 px, then r2 = Reff for panel (A), r2 = Reff/2 (B), and r2 = Reff/3 (C). The black circled cross signals the reference galaxy presented in Sect. 2.1 ($m=9.493,\overline{m}=30.64$).

All cases where the error is higher than the criterion, ε_rel, 90% ≥ 10, are shown with the same red colour as ε_rel, 90% = 10. In general, high brightness and low fluctuation magnitudes, that is, increasing the count number of both, returns lower relative errors. When increasing $\overline{m}$ the relative error ascends up to ε_rel, 90% = 10. We find errors higher than our criterion for $\overline{m}\approx 34-37$ mag, depending on the mask applied. When the fluctuation luminosity is larger than the galaxy luminosity itself, our mock galaxy is no longer physically realistic. This happens for lower fluctuation magnitudes than $\overline{m}\approx 24$.

The empty region found in the right and bottom-right corner represents galaxies where the third criterion of Eq. (25), $N({\mathrm{Gal}}_{\mathrm{mock}\mathrm{mask}})-N({\mathrm{Sky}}_{\mathrm{mask}})\ge 10\times \overline{N}$, is not fulfilled, therefore the SBF is not obtained. In our scenario, total galaxy magnitudes higher than m ≈ 10.5 − 13.5 (depending on the mask) do not fulfil the condition and neither do cases where $\overline{m}-m\gtrsim 19$. Such a region is larger for larger masks, as we are considering regions of the galaxy with lower brightness. We observe that using smaller masks increases the relative error. As we have shown in Sect. 3.1, the larger the number of pixels within the mask, the lower the relative error.

It is worth highlighting the differences when studying a galaxy through the integrated surface luminosity, or when studying it through surface brightness fluctuations. The first case depends on the luminosity profile, so a few pixels with high luminous stars can dominate the flux. In the second case, the SBF strongly depends on the number of pixels with a given fluctuation, since the SBF is defined as the luminosity normalised by variance. Additionally, the SBF range of possible values is more restricted, unlike the large variation of the integrated surface luminosity. The integrated surface luminosity and the SBF complement each other and a combination of both provides more constrained information about the galaxy (Rodríguez-Beltrán et al. 2021).

3.3. Varying other parameters

After studying general SBF measurements related to the brightness of the galaxy and the size of the masks applied (their number of pixels), we proceed to analyse other parameters such as: the PSF size (σ_PSF), the exposure time (t_exp), the sky number of counts (N_Sky), the Sersic index (n) and the effective radius (R_eff). Again, we use as a reference the galaxy described in Sect. 2.1, then, we vary a certain parameter while fixing the rest. In the figures of this subsection we show the reference case as black dashed line. For every case we show the relative error (ε_rel, 90%) against the number of pixels within a mask (n_mask pix). In Sect. 3.1 we demonstrated how the uncertainty is similar in masks with the same number of pixels. Therefore, we make use n_mask pix as a proxy for our analysis, with the procedure we explain below.

First, we select which parameters are held constant based on our reference galaxy and which one is varied. Second, for this chosen configuration we perform simulations for every mask combination (r₁, r₂), as we did in Fig. 5. Third, the results are grouped into subsets with similar number of pixels. The subset size is determined by the maximum number of pixels, which corresponds to a mask with r₁ = 4 and r₂ = R_eff, divided in 10 partitions. These subsets are analogous to each one of the regions displayed with different colours in panel (D) of Fig. 5 (although this figure is divided in 20 subsets, instead of 10). Fourth, we repeat the process of the previous two points (second and third) using 50 Monte Carlo simulations. Fifth, from this group of simulations we combine every subset with similar number of pixels. This is, every subset is now fifty times larger. Sixth, from each one of these subsets we calculate the 90th percentile of the relative error. With this method we ensure that the percentile is derived from masks with a similar number of pixels, and that it is applied simultaneously to the 50 simulations.

We note that, as expected, in all our experiments the relative error generally decreases when the number of pixels increases. The curves do not descend smoothly due to the finite number of simulations that we have performed and their random nature. However, even our limited number of simulations is sufficient to address the general behaviour shown in the plots.

3.3.1. PSF size

In Fig. 7 we show how the relative error changes when varying the PSF width, using 3 × σ_PSF = 2, 3, 4, 6 and 8 px. The power spectrum shrinks in the frequency domain (k-axis) when the PSF is wider in the physical domain (px-axis). In order for the comparison to be fair, we adjust the range of frequencies used for the fitting with respect to the point spread function of reference (σ_PSF, ref) from Sect. 2.1 and the fitting frequencies of reference (k_fit, ref) from Sect. 2.5.2. The new fitting frequencies are calculated as k_fit = k_fit, ref/(σ_PSF/σ_PSF, ref).

Fig. 7. Upper 90% relative error against the number of pixels within the mask when varying the PSF pixel size: 3 × σPSF = 2 is in blue, 3 × σPSF = 3 in cyan, 3 × σPSF = 4 in dashed black, 3 × σPSF = 6 in lime yellow, and 3 × σPSF = 8 in orange.

The ε_rel, 90% worsens as the physical width of the PSF increases: the power spectrum is narrower, the contribution of the noise is enlarged and there is less relevant information to fit. Hence, in the physical plane the smaller the PSF the better, nonetheless, we remind that the presence of the PSF is required to disentangle the instrumental noise from the spatially correlated SBF signal.

3.3.2. Exposure time

In Fig. 8 we show ε_rel, 90% versus n_mask pix while changing the exposure time t_exp = 411, 822, 1233, 2055 and 4110 s. This is equivalent to changing the number of counts of N_Reff, $\overline{N}$ and N_Sky. We find that very short exposure times lead to larger relative errors, as the number of counts of the image is not high enough compared to the instrumental noise source.

Fig. 8. Upper 90% relative error against the number of pixels within the mask when varying the exposure time (in seconds) of the mock galaxy: texp = 411 is in magenta, texp = 822 in blue, texp = 1233 in cyan, texp = 2055 in dashed black, and texp = 4110 in lime yellow.

3.3.3. Sky background

In Fig. 9 we show ε_rel, 90% versus n_mask pix while changing the sky background of the image (N_Sky), which is similar to varying the S/N of the image. We show the results for a lower and higher sky count with respect to Sect. 2.1: Sky = 2288, 4575, 9150, 18 300, 36 600 and 73 200 counts.

Fig. 9. Upper 90% relative error against the number of pixels within the mask when varying the sky background count number: NSky = 2288 is in blue, 4575 in cyan, 9150 in dashed black, 18 300 in lime yellow, 36 600 in orange, and 73 200 in brown.

For a given galaxy flux (Gal_mean, fluc), we find that decreasing the sky background reduces the relative error, while increasing the sky worsens the SBF retrieval, up to a limit where the criteria of this work are not fulfilled (around Sky ≈ 120 000 counts). These results are similar to those found in Fig. 6, where instead of varying the sky, we vary the luminosity of the galaxy (both m and $\overline{m}$, for a fixed sky value). The instrumental noise increases with higher sky values. If there is not enough contrast between the correlated noise (the SBF with the PSF or, traditionally, P₀) and the uncorrelated noise (instrumental noise or P₁) the fitting will worsen. This is different from a bad evaluation of the observed sky, introducing an offset, which also leads to high relative errors (see Sect. 4.3).

3.3.4. Sersic index

In Fig. 10 we test changing the Sersic index n = 1, 2, 3, 4, 6 and 10, which is equivalent to varying the steepness of the light profile. We find that, for our reference galaxy at least, the Sersic index is not a dominant factor when evaluating the SBF fitting.

Fig. 10. Upper 90% relative error against the number of pixels within the mask when varying the Sersic profile of the mock galaxy: n = 1 is in magenta, n = 2 in blue, n = 3 in cyan, n = 4 in dashed black, n = 6 in lime yellow, and n = 10 in orange.

3.3.5. Effective radius

In Fig. 11 we vary the size of our mock galaxy by changing its effective radius R_eff. We take the effective radius presented in Sect. 2.1 (R_eff, ref) and present the fitting for galaxies with a size proportional to this radius (0.5R_eff, ref,  0.75R_eff, ref,  R_eff, ref,  1.25R_eff, ref,  1.5R_eff, ref). The range of each line is limited by the maximum number of pixels of each effective radius, as we only study masks from r₁ = 4 px to r₂ = R_eff.

Fig. 11. Upper 90% relative error against the number of pixels within the mask when varying the effective radius of the mock galaxy: 0.5 × Reff, ref is in blue, 0.75 × Reff, ref in cyan, Reff, ref in dashed black, 1.25 × Reff, ref in lime yellow, and 1.5 × Reff, ref in orange.

As expected, the number of pixels limits the reliability of the fitting, for instance the case of 0.5R_eff, ref does not reach ε_rel, 90% ≲ 4%, while R_eff, ref and 1.5R_eff, ref do. Aside from limiting the region to study, larger effective radii only increase slightly the relative error when sharing the same number of pixels.

4. Discussion

The results presented in Sect. 3 together with the diverse literature previously discussed in Sect. 1 are aimed at addressing the reliability and applicability of the SBF retrieval. In this section we discuss the limitations of the SBF computations, such as some of the approximations taken in the procedure, calculating biased measurements and the possibility of measuring robust SBF gradients.

4.1. Reliability of SBF derivation

In Sects. 2.3, 2.4 and Appendix A we describe the mathematical development for deriving the SBF. Some approximations must be taken to derive the final expression, Eq. (21). Hence, in this section we analyse the influence of the crossed term, the azimuthal average, the operational order when applying the convolution theorem and the fitting of the instrumental noise.

First, after applying the power spectrum (Eq. (16) or Eq. (A.3) with a mask), a crossed term appears (Eq. (17) or Eq. (A.4)) which, to our knowledge, is not discussed in the literature. This crossed term appears to be negligible. In Fig. 3 we test this by calculating the contribution of the crossed term. We found that the crossed term has a null imaginary component and a real component three orders of magnitude lower than the radial power spectrum of the galaxy. We presently do not have an interpretation for the non-zero real component of the crossed term. We demonstrate this contribution is insignificant by fitting the SBF of the galaxy presented in Fig. 1 considering the crossed term as in Eq. (A.7) (the galaxy is masked as in Fig. 4). The results are shown in Table 2, where we present the upper 90% percentile relative error, the 90% percentile of the relative standard deviation and the relative 90% width of the ${\overline{N}}_{\mathrm{fit}}$ distribution, all from 50 simulations of the galaxy. In column 1 (named ‘1 Param.’) we present the results for the fitting obtained as in Fig. 4 and the rest of the work, that is, using Eq. (21). In Col. 2 (named ‘C.T.’) we present the results for the fitting when considering the crossed term. Both columns are equal, finding again that the crossed term does not change the results of the fitting of ${\overline{N}}_{\mathrm{fit}}$ (at least for this case, up to the fifth decimal digit).

Second, another source of uncertainty appears in Eq. (18) (or Eq. (A.6) when masked), where an azimuthal average is applied to the power spectrum of the images. Each point in the profile found after performing this average has a scatter, since this image is noisy and does not necessarily present radial symmetry. Therefore, the resulting distribution of values for a fixed frequency k will follow, in general, a non-Gaussian asymmetric distribution. This is shown in the top panel of Fig. 12 where, for each k value, we display in pairs, the lower 16% and higher 84% percentiles and the lower 32% and higher 68% percentiles of the scatter associated with the distribution. Note that the azimuthal average is close to the 68% percentile instead of the 50% percentile, as it should happen in a symmetrical distribution. This means that, if the SBF extraction procedure makes use of an azimuthal mode or median instead of the azimuthal average, it could lead to distorted SBF results. In addition, we have an estimate of the standard deviation of the averaged value at each k, σ_az(k), which is shown in the bottom panel of Fig. 12. In this manner, we can perform a weighted fitting¹⁵ by considering this σ_az(k). The results obtained when performing the weighted fitting are found in column 3 (named σ_az) of Table 2, showing an increase in both the relative error and the relative standard deviation with respect to our standard way of fitting. This indicates that including the σ_az(k) margins of the radial profile is a more conservative way of fitting. On the other hand, the relative 90% width of the ${\overline{N}}_{\mathrm{fit}}$ distribution (Δ_90%) is slightly lower than the original fitting. This shows less dispersion (better precision) in the fitting results when applying the σ_az(k) margins.

Fig. 12. Power spectrum (not rescaled) similar to the one of Fig. 4, using the same data and mask. Top panel: blue line as the radial power spectrum of the observed fluctuation, PS(Galmock fluc mask)r, with a scatter associated with the distribution of values found at each radius (lower 16% and higher 84% percentiles in pale-blue and lower 32% and higher 68% in pale-cyan); cyan line as the radial power spectrum of the masked instrumental noise term, $PS(R_{\mathrm{approx}}\times \mathrm{Mask}/{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}_{\mathrm{r}}=PS{(R_{\mathrm{approx}}^{\mathrm{norm}})}_{\mathrm{r}}$; the range of frequencies where the fitting is performed (kfit, i = 75 px−1, kfit, f = 400 px−1) is marked with a pale-yellow vertical region. Bottom panel: standard deviation associated with the azimuthal average (σaz).

Third, when tackling the fitting procedure, we contemplate different options. The results presented in this work consider a single fitting of the SBF as an unknown parameter, because we model the Poisson noise as in Eq. (22). Most authors fit both the SBF and the instrumental noise term simultaneously, this is, with the traditional nomenclature, fitting together P₀ and P₁ (e.g. Pahre et al. 1999; Mitzkus 2017). Other authors have already studied how to model the noise and the sky (e.g. for optimised drizzling algorithms, Mei et al. 2005), although it was not directly applied in the SBF fitting. In column 4 (named ‘2 Param.’) of Table 2 we present the results of the dual fitting of the SBF and PS(R^norm)_r. We find an increase in the 90% percentile of the relative error, the 90% percentile of the relative standard deviation and the relative 90% width of the ${\overline{N}}_{\mathrm{fit}}$ distribution, with respect to our fitting. This result shows how modelling the instrumental noise helps in measuring a more precise SBF, instead of fitting both parameters.

Fourth, we review a certain step in the SBF derivation procedure that, to the best of our knowledge, was not previously assessed. Tonry et al. (1990) defines the expectation power spectrum (E(k)) as the convolution of PS(PSF)⊗PS(Mask), which is scaled by the SBF value. Using the expressions of this work, the calculation is performed in the following order: ${\left|\mathfrak{F}\left(\frac{{\mathrm{Gal}}_{\mathrm{fluc}}}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right)\right|}^2\times ({\left|\mathfrak{F}\left(\mathrm{PSF}\right)\right|}^2\otimes {\left|\mathfrak{F}\left(\mathrm{Mask}\right)\right|}^2)$. Instead, the rigorous order according to Eq. (A.5), consists of first multiplying the fluctuation term by the PSF and, then, convolving the result with the mask, this is, $({\left|\mathfrak{F}\left(\frac{{\mathrm{Gal}}_{\mathrm{fluc}}}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right)\right|}^2\times {\left|\mathfrak{F}\left(\mathrm{PSF}\right)\right|}^2)\otimes {\left|\mathfrak{F}\left(\mathrm{Mask}\right)\right|}^2$. We note that the order followed conventionally by the literature is correct only if the fluctuation term, $PS({\mathrm{Gal}}_{\mathrm{fluc}}/\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}})$, is constant. And so it appears to be, at least for the experiments of this work, as we demonstrated in Sect. 2.3, with Fig. 2. In a galaxy with an SBF gradient (as presented in Sect. 4.4) we find the same constant behaviour. To this extent, and given our ideal experiments, the above approximation is valid.

In summary, the fitting approach presented in this work appears to be a proper estimation for the SBF measurement. We note that considering the uncertainties when performing the azimuthal average shows a more conservative estimation for ε_rel, 90%) and σ_fit, 90%, but improves the precision of the fitting. Modelling the instrumental noise (R), instead of fitting it, reduces the uncertainty of the calculation (see Table 2). We suggest modelling the instrumental noise as in Eq. (22), which requires knowledge of the mean galaxy value, the PSF and the sky background, as well as considering Poisson noise. In this regard we encourage adapting, if necessary, the equation for each observation or studying other procedures, such as the one presented in Mei et al. (2005) for correlated noise. The rest of the approximations taken during the SBF derivation are negligible (i.e. the crossed term and the order after applying the convolution theorem).

4.2. Biased SBF measurements due to low flux levels

This section is intended to point out an important caution to consider: examining a low flux source using the traditional SBF extraction methodology could potentially introduce a bias into the result. Here, we explore how this bias appears, how to avoid or mitigate it, as well as the conditions necessary for applying confidently the standard SBF derivation.

To begin with, we analyse how this bias can appear in the procedure. A reference ‘mean image’ (Gal_mean(x, y)) always can be obtained by different methods, such as applying a mean filter, smoothing the image, by an isophote fitting, with Sersic profiles or other approaches. Usually, this image is used to obtain an SBF measurement by subtracting it to the original image, dividing the result by its square root and obtaining the power spectrum of the resulting image. However, having an adequate Gal_mean(x, y) model image of the mean brightness profile is not a sufficient condition: it is required that this is a proper representation of the mean of the stellar population luminosity distribution¹⁶ in each pixel (⟨ℒ(x,y)⟩)¹⁷ and, consequently, obtaining the variance (${\sigma }_{\text{fluc}}^2(x,y)$) of such a distribution. We recall that the mean and variance of the population luminosity distribution scale linearly with the number of stars in each pixel. This number of stars must cancel out, as a way for all the pixels to be equivalent when measuring the fluctuation (Tonry & Schneider 1988; Cerviño & Luridiana 2006; Cerviño et al. 2008; Cerviño 2013). Otherwise, it cannot be used neither for stellar population analysis nor for distance calculations.

To illustrate this, let us consider a scenario in which we have made a biased estimate of the mean of the population luminosity distribution along the galaxy pixels, ⟨ℒ(x,y)⟩_real, by an amount of ±L_offset(x, y), which could vary depending on the (x, y) positions (see Cerviño et al. 2008). That is:

$$\begin{array}{c}\hfill {⟨\mathcal{L}(x,y)⟩}_{\mathrm{biased}}={⟨\mathcal{L}(x,y)⟩}_{\mathrm{real}}\pm L_{\mathrm{offset}}(x,y).\end{array}$$(27)

A simple calculation¹⁸ shows that, if we use ⟨ℒ(x,y)⟩_biased as the reference value at a position (x, y), its associated variance is:

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{fluc},\mathrm{biased}}^2(x,y)={\sigma }_{\mathrm{fluc}}^2(x,y)+L_{\mathrm{offset}}{(x,y)}^2,\end{array}$$(28)

and, therefore, we obtain a biased fluctuation at that position:

$$\begin{array}{c}\hfill {\mathrm{SBF}}_{\mathrm{biased}}(x,y)=\frac{{\sigma }_{\mathrm{fluc},\mathrm{biased}}^2(x,y)}{{⟨\mathcal{L}(x,y)⟩}_{\mathrm{biased}}}=\frac{{\sigma }_{\mathrm{fluc}}^2(x,y)+L_{\mathrm{offset}}{(x,y)}^2}{{⟨\mathcal{L}(x,y)⟩}_{\mathrm{real}}\pm L_{\mathrm{offset}}(x,y)},\end{array}$$(29)

which deviates from the actual definition of the SBF (i.e. ${\mathrm{SBF}}_{\mathrm{real}}(x,\mathit{y})={\sigma }_{\mathrm{fluc}}^2(x,\mathit{y})/{⟨\mathcal{L}⟩}_{\mathrm{real}}(x,\mathit{y})$), with no straightforward method of eliminating this offset by factorisation. The situation is worse when several pixels are taken into consideration, since the fluctuations measured along the pixels cannot be compared in a common framework, and it is not assured that there is the required independence of the SBF with the number of stars per pixel.

To assure a correct mean value of the population luminosity distribution it is required that all possible evolutionary phases are well-sampled along the pixels. This is achieved if there is a large enough number of stars distributed in these pixels, about a total of ∼10⁷ stars, depending on the band (Cerviño & Luridiana 2004, 2006). In practice, the galaxy model Gal_mean(x, y) is commonly obtained from an observed image by considering the flux of nearby pixels and making some kind of average for each collection¹⁹ of them. Therefore, it is mandatory that every collection has enough stars to represent robustly the whole stellar population. Additionally, the pixels of each collection should have similar characteristics, as number of stars and stellar populations properties.

In this context, we could analyse two opposite cases: (a) if the pixels of the observation gather low number of stars per pixel, each collection needs to assemble a large enough number of pixels to achieve this statistically meaningful estimate of the mean; or, (b) if every pixel in the observation has a sufficient number of stars, we do not need a necessarily large number of pixels per collection, as in a traditional SBF derivation. Both cases are fully understood when knowing the shape of the population luminosity distribution function, which has been studied in detail in Cerviño & Luridiana (2006) and Cerviño et al. (2008). These works show that the distribution of the integrated luminosity in the collection of pixels follows an L-shape in the extreme case (a) with a single star per pixel, and a Gaussian shape in case (b) with infinite stars per pixel. There is a continuous transition between both cases, dependent on the number of star per pixel.

Remember that in most studies the methodology applied considers a situation as case (b), typical of elliptical galaxies, where a traditional SBF derivation is considered. The mean value of each pixel is obtained with a reasonable small group of neighbouring pixels, each with sufficient number of stars. Here, we warn about the risks of studying an observation of case (a) applying the traditional procedure made for case (b). In doing so, the mean is estimated from neighbouring pixels that might be dominated by a few luminous stars. If those pixels are considered, they lead to an overestimation, while, if absent, they lead to an underestimation (e.g. if pixels with extremely bright stars are masked due to its luminosity excess). In both situations we introduce an offset, such as illustrated in Eq. (29).

To tackle this situation it is necessary to stablish a threshold in which we are confident to be working in scenario (b). In this sense, Cerviño & Luridiana (2004) states how many stars are required for a statistically meaningful distribution to be found. A Gaussian-like regime can be reached with ∼10⁶ stars per pixel for the optical bands (even a larger number for the infrared bands), assuming a simple stellar population and a standard Initial Mass Function (Cerviño & Luridiana 2006). This requirement is similar to the assumption of having about 20 giant stars per pixel in old stellar populations, quoted by Tonry & Schneider (1988), in the visible and infrared bands. In the case of an old stellar population, the SBF flux value is similar to the flux of a single giant star (Tonry & Schneider 1988). Hence, we find the requirement that the SBF should be obtained from pixels with fluxes larger than ∼20 times the SBF flux. Such a criterion can be extended to more complex stellar populations, as inferred from Cerviño & Luridiana (2004, 2006). In Sect. 2.6 we relaxed this criterion to 10 times the SBF flux ($N({\mathrm{Gal}}_{\mathrm{mock}\mathrm{mask}})-N({\mathrm{Sky}}_{\mathrm{mask}})\ge 10\times \overline{N}$, in count numbers). This is a quick and easy test to minimise possible effects of non-Gaussianity in the stellar population distribution²⁰. Finally, aside of having enough stars per pixel, it is desirable to avoid too steep light profiles or abrupt flux changes, as shown in Cerviño et al. (2008).

In conclusion, in this section we highlight the problems derived when obtaining the SBF following the standard procedure, which assumes the condition of high number of star per pixel (case (b)), when the observed source has a low number of stars per pixel (case (a)). The caveats discussed here are particularly relevant when studying external parts of galaxies, faint sources, diffuse or dwarf galaxies (see references of Sect. 1). For instance, we found this scenario in Fig. 6, for low surface brightness and high fluctuation counts. In that figure, there is a blank region where the criterion of ‘studying pixels with, at leats, fluxes of 10 times the SBF’ is not fulfilled. For these cases, even if the estimated errors were low, the mean reference value is potentially incorrect. Then, the chance of obtaining a biased SBF increases in a non-trivial manner, so it is less reliable for the calibration of distances or stellar populations. Therefore, when working with an observation it is advisable to verify whether all the pixels utilised in the measurement meet the criterion of Eq. (25), once the SBF has been estimated.

4.3. Other sources of bias

Beyond the experimental conditions of an SBF measurement, such as an adequate number of counts associated with each pixel or achieving a favourable S/N, three essential parameters must be accurately estimated: the mean galaxy image, the PSF, and the sky background. These terms are responsible for obtaining the fluctuation image (Gal_mock fluc) and modelling the instrumental noise (R_approx). Additionally, the mean image and the PSF normalise Eq. (15) by $\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}$, following the SBF definition. Currently, the SBF measurement procedure is already well established by using several PSF templates, running experiments on the sky background and the galaxy profile until is well determined. Nevertheless, in this section we want to highlight how the modelling of these parameters should not be treated lightly.

On the one hand, even with enough counts in each pixel (as discussed in the previous section), when Gal_mean is not properly calculated, an offset is introduced in the SBF definition, as addressed in Eq. (29). On the other hand, if a uniform sky is assumed in an image with a non-uniform real sky, the estimated value of N_Sky would present an offset with respect to the real one, at least in certain areas of the image. Even if the observed sky is flat, it carries its inherent uncertainty. Finally, the PSF is commonly obtained from one or more foreground stars using different techniques, introducing variability in the PSF characterisation. In addition, although not considered in our following examples, in a real observation, a wrong sky estimation implicitly could lead to an erroneous mean galaxy model (Gal_mean).

To assess how an offset in any of these terms affects the SBF measurement, we propose some example experiments where we perform 50 simulations underestimating or overestimating the values of either Gal_mean, N_Sky or σ_PSF of our reference galaxy from Fig. 4, modifying one parameter while leaving the other two fixed. We choose to introduce offsets of χ_offset = ±3% and ±9%, applied as, for instance, Gal_{mean, offset} = Gal_mean ± χ_offset × Gal_mean. In the case of the PSF, it represents an offset of ±0.04 and ±0.12 px with respect to the reference value of σ_PSF = 1.33 px. For the sky, it represents an offset of ±275 and ±824 counts with respect to the original value of N_Sky = 9155 counts. The results for the mean fitted SBF $⟨{\overline{N}}_{\mathrm{fit}}⟩$ and the 90% percentile of the relative error (ε_rel, 90%) are presented in Table 3. They should be compared to those of the reference case (Fig. 4), which has an input SBF of ${\overline{N}}_{\mathrm{ref}}=22.59$ and, after performing 50 simulations, $⟨{\overline{N}}_{\mathrm{fit}}⟩=22.53$ counts and ε_rel, 90% = 2.88%. In Table 3 we do not show neither the relative standard deviation nor the precision, as they do not change significantly compared to the results found for Fig. 4.

The results presented in Sect. 3 have a variety of different relative errors depending on the circumstances, but their mean value of the fitted SBF are found surrounding ${\overline{N}}_{\mathrm{real}}$. However, when an offset is introduced, the results of $⟨{\overline{N}}_{\mathrm{fit}}⟩$ show a systematic bias with respect to the reference case of Fig. 4, as it is shown in Table 3. We note that for all the cases presented in Table 3 the accuracy worsens with the offset. In this example the most sensitive parameter is the sky, followed by the mean galaxy model and the PSF. Underestimating Gal_mean returns worse accuracies than overestimating it, due to the reduced contrast between the sky background and the galaxy light.

We note that the impact of a biased sky depends on the luminosity of the galaxy (Gal_mean, fluc). A brighter galaxy will not be as highly distorted by an offset in the sky. Similarly, when applying a mask, the relative error will be lower when studying the inner (and more luminous) regions than the outer parts of the galaxy. For instance, the ±275 counts of sky offset introduced in the previous example quantify differently depending on the number of counts at each part of the galaxy: at the effective radius, with N_Reff = 1494.96 counts, it represents an ±18% offset; at the outer radius of the mask (r₂), with N_Reff = 2492 counts, the offset is ±11%; at the inner radius of the mask (r₁) with N_Reff = 21 132 counts it represents a ±1% offset. As already shown in Sect. 3, studying bright galaxy pixels with respect to the sky count value facilitates retrieving a reliable SBF.

To summarise, these results show how an offset due to a non-properly estimated mean model, PSF or sky background, could severely increase the relative error. In particular, the accuracy drastically worsens when the sky is not properly measured. Meanwhile, the relative standard deviation and the precision are not drastically affected. In a real observation, we are not able to obtain the accuracy, but the relative standard deviation and the precision can be used as measurements of the uncertainty. Therefore, it is worth emphasising that achieving a favourable relative standard deviation or precision values might lead into misinterpretations of how close our fitted SBF is from the real fluctuation value. Even if this section aimed at pointing out how the results vary due to a biased modelling, generally such issues can be detected when checking the power spectrum plot used for the fitting (as the one shown in Fig. 4), where any inconsistency should be apparent and serves as a warning that an offset is being introduced.

4.4. SBF gradient detection

As mentioned in the introduction, some authors have investigated SBF gradients in galaxies. Several studies have found such SBF profiles (Cantiello et al. 2005, 2007; Sodemann & Thomsen 1995a,b; Jensen et al. 2015), but often show the radial outline to be relatively small, uncertain, or not conclusively detect the variability under given observational conditions (Jensen et al. 1996). Thus, we study the possibility of SBF gradient detection in a controlled environment, with the intention of establishing under what conditions such gradients (if present) would be found, and what is the optimal procedure to obtain them.

First, we need to determine what is the actual SBF value when the fluctuations are not homogeneous in all regions of a studied mask. We have found that the value recovered after the fitting is a weighted average of all the SBFs present in the masked region. This is, the sum of every SBF value (${\overline{N}}_1,{\overline{N}}_2,...,{\overline{N}}_f$) times the number of pixels where it appears (n_pix, 1, n_pix, 2, ..., n_pix, f), divided by the total number of pixels within the mask (n_mask pix):

$$\begin{array}{c}\hfill {\overline{N}}_{⟨\mathrm{SBF}⟩}=\frac{1}{n_{\mathrm{mask}\mathrm{pix}}}{\displaystyle \sum _{i=1}^f{n}_{\mathrm{pix},i}{\overline{N}}_i.}\end{array}$$(30)

Note that the SBF of an individual pixel is not measurable in observations, as more than one pixel is required for measuring a variance, but in practical terms we do know at each point of our mock galaxy the integrated luminosity distribution associated with the stellar population in that point and, hence, its associated SBF value (see Cerviño & Luridiana 2006).

We study the presence of a variable SBF using the reference galaxy presented in Sect. 2.1 and some complementary literature information. We built a gradient based on the work of Martín-Navarro et al. (2018). From that work, in their Fig. 5 the average stellar population gradient versus the radius is analysed. In particular, we focus on the metallicity profile for galaxies with a velocity dispersion over 300 km s⁻¹, typical of massive galaxies. We obtain a metallicity of [M/H]≈0.26 for log₁₀(r/R_eff) = − 1.1 and [M/H]≈0.06 for log₁₀(r/R_eff) = − 0.3. Additionally, the age is approximately constant with radius, with t ∼ 12.5 Gyr. We can relate these values to SBF magnitudes using the models presented in Vazdekis et al. (2020), their Fig. 9. There, it is shown how a metallicity of [M/H]≈0.26 corresponds to ${\overline{M}}_i=0.1$ mag, while [M/H]≈0.06 corresponds to ${\overline{M}}_i=-0.25$ mag, both for a 13 Gyr. Then, we calculate the apparent SBF magnitude for our reference galaxy, using the distance presented in Sect. 2.1, D = 16.7 Mpc. We find ${\overline{m}}_i=31.21$ mag for log₁₀(r/R_eff) = − 1.1 and ${\overline{m}}_i=30.86$ mag for log₁₀(r/R_eff) = − 0.3. Finally, if we follow the gradient shown in Fig. 5 from Martín-Navarro et al. (2018) for the metallicity, we can approximate a linear profile for ${\overline{m}}_i$ versus log₁₀(r/R_eff). Thus, we can calculate ${\overline{m}}_i(r)$ for any given radius as a linear function and then transform it into counts using Eq. (6). Note that the fluctuation count increases with the radius.

Once a galaxy with an SBF gradient is created, we test the possibility of an actual SBF gradient detection. In Fig. 13 we show the radial profile of the SBF gradient in counts as a black line and its respective profile in magnitudes as a grey line. We mark with cyan vertical regions each of the annular masks used to fit the SBF in different annuli, from 4 pixels up to the effective radius. Each one of these masks share approximately the same number of pixels (n_mask pix ≈ 75 000 px), this way we assure that the uncertainties are comparable. Within an annular mask, the weighted average SBF, obtained with Eq. (30), is associated with a radius that depends on the SBF profile. In Fig. 13, we chose the radius (r_⟨SBF⟩) such that the value of the SBF at that distance from the centre matches the weighted average SBF of the studied annulus, this is SBF(r_⟨SBF⟩) = ⟨SBF⟩. We perform the SBF estimation for 50 galaxy simulations. Then, we present the mean SBF fitting value as red dots, with an error interval defined with the 5th% minimum and the 95th% maximum percentile values (as explained in Eq. (26)) among the 50 iterations of the mock galaxy. We show the results at each mask for two different example mock galaxies retrieved from the 50 simulations, using green and yellow star-markers.

Fig. 13. Variable SBF radial profile shown in counts as a black line and in magnitudes as a grey line. Annular masks are indicated as cyan vertical regions. For each mask, the SBF mean fitting results over 50 simulations are marked as red dots, as well as its error bar. Results at each mask of two different example galaxies (retrieved from the 50 simulations) are shown with green and yellow star markers.

In Table 4 we show in more detail the results for each one of the masks. Column 1 presents the radii (r₁, r₂) that define each annular mask, column 2 shows the exact number of pixels within the mask (n_mask pix), column 3 marks the radius associated with the weighted average SBF of each mask (r _⟨SBF⟩), column 4 shows the weighted average of the SBF in counts (N_⟨SBF⟩), column 5 corresponds to the measured SBF in counts ($⟨{\overline{N}}_{\mathrm{fit}}⟩$, with a lower and upper limit calculated as $⟨{\overline{N}}_{\mathrm{fit}}⟩-{\overline{N}}_{\mathrm{fit},5\mathrm{th}\%}$ and ${\overline{N}}_{\mathrm{fit},95\mathrm{th}\%}-⟨{\overline{N}}_{\mathrm{fit}}⟩$, respectively), column 6 and 7 show the higher 90% percentile of the relative error and the relative standard deviation (ε_rel, 90% and σ_fit, 90%) respectively, column 8 presents the relative 90% width of the ${\overline{N}}_{\mathrm{fit}}$ distribution (Δ_90%), column 9 is the magnitude of the weighted average SBF and column 10 is its corresponding SBF measured magnitude (m_⟨SBF⟩ and ${\overline{m}}_{\mathrm{fit}}$ magnitudes are retrieved from columns 4 and 5 making use of Eq. (6)).

The results presented in Fig. 13 and Table 4 serve as a test for the detection of SBF gradients. With $⟨{\overline{N}}_{\mathrm{fit}}⟩$ we recover the profile correctly, with relative errors lower than ε_rel, 90% < 7.5% and relative standard deviations lower than σ_fit, 90% < 1%. The relative 90% width of the ${\overline{N}}_{\mathrm{fit}}$ distribution shows a range of 9.8%< Δ_90% < 15.9% precision. All three parameters are similar among the different masks, as they approximately share the same number of pixels. The mean value of the 50 solutions of the fitting follows its corresponding weighted average SBF, from $⟨{\overline{N}}_{\mathrm{fit}}⟩=15.8$ counts in the inner mask up to $⟨{\overline{N}}_{\mathrm{fit}}⟩=20.4$ counts in the external mask. Additionally, notice how the position of the radius associated with the weighted average SBF (r_⟨SBF⟩) tends to move towards the region with the larger number of pixels. On this account, the variation of $⟨{\overline{N}}_{\mathrm{fit}}⟩$ through the consecutive masks could be detectable, however, if we consider the uncertainty associated with the Δ_90% width of every annulus, only the results between the first mask and the rest of masks do not overlap. Any other mask comparison, besides the first annulus, has sufficiently large uncertainties for the results not to be distinguishable. For a clearer understanding of the behaviour, we presented in Fig. 13 those two particular simulated galaxies (green and yellow stars) at each mask. In both cases the fitting solutions are found arbitrarily above and below the radial profile of the SBF gradient. This indicates how tracing the gradient with a sole observation could lead to a correct representation of the actual SBF profile. However, it also could lead to misinterpretations if, for instance, most of the mask results fall below the real SBF profile. Among our simulations, we have found cases that both underestimate and overestimate the SBF slope.

In actual observations we do not have 50 fittings from which we extract a mean value, hence, a sole simulation will trace the gradient of the galaxy but with varying points around the mean. So, even working under ideal conditions (as in this work) presents difficulties when comparing the SBF from radius to radius. At least, we can calculate the precision of the results in each annuli, as we advanced in Sect. 3 and we explain in Sect. 6, also illustrated with Fig. 14. To do so, we propose constructing a simulated set of Monte Carlo galaxies, as described in Sect. 2.2, but considering as input the observed SBF of the different annuli, ${\overline{N}}_{\mathrm{obs}}$. Subsequently, we can determine σ_fit, 90% and Δ_90% with respect to the fitted ${\overline{N}}_{\mathrm{obs}}$, following the procedure explained in Sect. 6 and illustrated with Fig. 14. In any case, steeper SBF gradients and better observations would ease constraining the radial profile of the fluctuation.

Fig. 14. Flowchart of the methodology followed throughout this work and applied in the code provided21.

Additionally, we replicate the fitting performed in Fig. 4 (between r₁ = 80 px and r₂ = 332 px), which would be considered as the global SBF of the galaxy, but keeping the SBF gradient of this section, instead of a constant SBF. We have that $\overline{N}(r_1=80)=15.46$ and $\overline{N}(r_2=332)=19.84$ counts. Within this mask, the weighted average of the SBF (from Eq. (30)) and the fitted SBF return the same value with $⟨{\overline{N}}_{\mathrm{fit}}⟩=18.49{\pm }_{0.39}^{0.44}(={\overline{N}}_{⟨\mathrm{SBF}⟩})$ counts and a relative error of ε_rel, 90% = 2.29%. In magnitudes, this is ${\overline{m}}_{\mathrm{fit}}=30.86{\pm }_{0.026}^{0.023}(={\overline{m}}_{⟨\mathrm{SBF}⟩})$ mag. This shows that measuring the SBF of a galaxy as a whole may not reflect the variability of the fluctuation in different regions, which could lead to larger scatter in the calibration. Thus, it is necessary to take into account the potential non-homogeneity of the SBF distribution when constraining stellar populations. When studying the SBF profile with different annular masks (as in Fig. 13), we can validate the partial SBFs using the area covered by each one and comparing with the SBF of the galaxy as a whole, using Eq. (30). We remind the reader that this also occurs in observations through the integrated light of the galaxy region.

5. Summary

In this paper, we have investigated the capabilities and limitations of the surface brightness fluctuation retrieval. To do so we have created a mock galaxy based on the data of NGC 4649 (VCC 1978) from Cantiello et al. (2018). Then, we present the SBF derivation in detail and our fitting procedure. Creating a mock galaxy allows us to study a measurement under ideal laboratory conditions and determine the fidelity of the fitting in relation to the introduced SBF. When performing the fitting we consider the mean galaxy model, the sky background and the PSF as known quantities. Consequently, the errors studied reflect the variability of a system with inherent stochasticity due to the stellar population luminosity fluctuations and the instrumental noise. We define two parameters to evaluate the measurement: the relative error (ε_rel), as a measure of the accuracy between the real SBF (${\overline{N}}_{\mathrm{real}}$) of the mock galaxy and the fitted SBF we obtain (${\overline{N}}_{\mathrm{fit}}$), and the relative standard deviation obtained from the least square fitting (σ_fit). We assume as a valid SBF measurement those where ε_rel < 10% and σ_fit < 10%. Additionally, we require all non-masked pixels of the galaxy to have values larger than 10 times the input SBF, to assure a Gaussian probability distribution of the integrated light among pixels (Tonry & Schneider 1988; Cerviño & Luridiana 2006) and be able to apply the traditional SBF derivation procedure.

We evaluate the SBF estimation under a wide range of conditions. To do so, we create a set of 50 mock galaxies each time a parameter is varied and obtain the distribution of SBF measurements for each set. From this distribution of results we calculate the 90% percentile of the relative error and relative standard deviation (ε_rel, 90% and σ_fit, 90%). Additionally, we calculate the relative 90% width of the ${\overline{N}}_{\mathrm{fit}}$ distribution (Δ_90%), to measure the precision of the fitting results.

We study the SBF fitting for different mask sizes, using all possible combinations of annular rings defined by an internal and external radius. The results show how the relative error (ε_rel, 90%), the relative standard deviation of the fitting (σ_fit, 90%) and the relative 90% width of the ${\overline{N}}_{\mathrm{fit}}$ distribution (Δ_90%) are tightly related to the number of pixels within the mask (n_mask pix). Therefore, the larger the number of pixels within the mask, the better the reliability of the measurement. The precision (Δ_90%) is the most restrictive parameter, followed by the accuracy (ε_rel, 90%) and then by the relative standard deviation (σ_fit, 90%). We focus our study on the relative error, as it is representative of the accuracy between the real SBF value (${\overline{N}}_{\mathrm{real}}$) and the fitted (${\overline{N}}_{\mathrm{fit}}$).

We vary the brightness of the galaxy and its fluctuation contribution. We find that the best SBF measurements results are found for both high brightness and high fluctuation. When the fluctuation contribution is too low, the relative error does not fulfil the 10% uncertainty threshold condition. When the brightness of the galaxy is too low we find that the pixels of the galaxy have lower count number than 10 times the input SBF. Under a traditional SBF extraction procedure this last condition is a first order requirement to obtain a physically reliable representation of the stellar population through the SBF. Otherwise, the estimate of the mean galaxy image can be biased and, so the SBF measurement (Cerviño & Luridiana 2004, 2006; Cerviño et al. 2008; Cerviño 2013).

We also varied other parameters in search of ideal conditions for the SBF measurement. We find that the relative error shrinks with a small PSF up to the size of the pixel. A low sky value and a long exposure time improve the results. For the experiments performed in this work, varying the Sersic index shows no difference in the 90% percentile of relative error (ε_rel, 90%). Varying the size of the galaxy, that is, the effective radius does not significantly change the result either, but does limit the maximum number of pixels in which we can perform the measurement.

Aside from varying the parameters of the galaxy, we study the reliability of the SBF derivation. We have reviewed its mathematical development and verified that any approximations taken do not drastically influence the measurement of the SBF. Additionally, before the SBF fitting, we suggest modelling the instrumental noise in order to improve the accuracy of the SBF as the sole unknown parameter (see Table 2). Of course, the modelling of the instrumental noise must adapt to the particular nuances of any observation (such as considering Poisson noise or any other). Applying Eq. (22) for modelling the noise requires proper approximations of the sky, the mean model and the PSF. Moreover, a poor estimation of these three parameters leads to biased SBF measurements. Also, we find that including the azimuthal margins found after creating the radial profile provides more conservative results, but improves the returned precision.

Finally, we study the possibility of SBF gradient detection by including a variable SBF in our mock galaxy. We propose an expression for measuring a region with a varying SBF, which corresponds to a pixel-weighted average of all the SBFs present in the masked region. In an actual observation the SBF value in each pixel is unknown, but a similar concept applies: an averaged fluctuation calculated with the different SBFs weighted by the area they occupy. If the galaxy fluctuation changes through different regions, measuring the SBF of the whole galaxy returns its weighted average. Consequently, it is crucial to take into account the heterogeneity of the SBF distribution when constraining stellar populations. This reaffirms previous knowledge, as the variability of the stellar population is already considered when studying the integrated colour of an observation. With our mock galaxies we attempt to recover the introduced SBF profile by performing the fitting in concentric annuli, using different masks. When looking for an SBF gradient, we suggest keeping the same number of pixels in every mask used, to ensure the consistency in the results. In our experiments, we can trace the introduced SBF profile using the mean fitted value of each annulus returned from our 50 simulations, allowing us to study the viability and uncertainty of an SBF gradient through the modelling. However, in an actual observation, there are no 50 versions of the observed galaxy from which to retrieve a mean value, so the result is limited to the sole observation performed in each annulus and its associated dispersion does not ensure a reliable gradient is recovered. The solutions will follow the real profile with varying points around the real value. In this case, at the very least we can calculate the precision of the results, as we explain in Sect. 6, illustrated in Fig. 14. In any case, working under ideal conditions is already challenging when comparing SBF from one radius to another. Better observations and a more pronounced SBF gradient will help to reveal the underlying fluctuation radial profile.

6. SBF uncertainty estimator: A code application

We provide a simple code that estimates the uncertainty of an SBF retrieval, (SBF_uncertainty_estimator²¹), due to the variability of a system that is intrinsically stochastic. It implements the procedure described in this document.

The aim of this code is twofold: (i) preliminary study the suitability of a forthcoming SBF measurement or (ii) estimate the uncertainty of an SBF observation. We explain below the steps that our code performs and the nuances in each case, (i) or (ii). Between both cases, the meaning of the introduced fluctuation (${\overline{N}}_{\mathrm{input}}$) has different connotations. The flowchart shown in Fig. 14 also illustrates the procedure:

1. Create a Monte Carlo set of simulated mock galaxies following the methodology described in this work.

– Introduce representative parameters of the target galaxy, including an initial estimate of its SBF, which can be taken from model predictions. Note that although the actual SBF of the galaxy is unknown, this estimate is needed to predict the quality of the SBF derivation. In the current work this value corresponds to ${\overline{N}}_{\mathrm{real}}$.

– Introduce the fitted SBF value (${\overline{N}}_{\mathrm{obs}}$) obtained from the observation, along with the rest of parameters associated with the galaxy.

2. Fit the SBF for each one of the mock galaxies of the simulated set, therefore retrieving a distribution of fitted SBFs from the simulations $[{\overline{N}}_{\mathrm{fit},1},...,{\overline{N}}_{\mathrm{fit},{\mathrm{n}}_{\mathrm{sim}}}]$. Also, we obtain the mean value of the fitted SBFs, $⟨{\overline{N}}_{\mathrm{fit}}⟩$.

– This distribution is conceptually equivalent to those presented throughout this work.

– This distribution would be a “second generation” of fitted mock SBFs with respect to the original observed SBF.

3. Calculate the uncertainties associated with the variability of a stochastic system (ε_rel, 90%, σ_fit, 90%, Δ_90%) from the distribution of fitted mock SBFs. Check for low flux pixels ($N({\mathrm{Gal}}_{\mathrm{mock}\mathrm{mask}})-N({\mathrm{Sky}}_{\mathrm{mask}})\ge 10\times {\overline{N}}_{\mathrm{input}}$) where the modelling might not be correct (see Eq. (25) and Sect. 4.2) and warn the user if necessary.

– With the obtained uncertainties we can predict if our object of study will be a reliable source for the SBF study. The code performs an ideal experiment, so, if the result is untrustworthy according to these simulations, then the situation in a real observation could be more problematic.

– The quality of the fitting (σ_fit, 90%) and the precision (Δ_90%) will be useful for assigning such error bars to the observed SBF. We recall that precision is a more restrictive parameter. Given that the input SBF is ${\overline{N}}_{\mathrm{obs}}$ and the real SBF is unknown, the relative error does not measure the accuracy of the true SBF value of the galaxy. It is a metric for comparing the observed SBF with those derived from the ‘second generation’ Monte Carlo simulations.

As in this work, the code assumes a Sersic profile, a fluctuation of the stellar population luminosity as a random Gaussian distribution, a flat sky, a 2-D Gaussian PSF and a random Poisson distribution for the instrumental noise. No GC nor background sources are considered. The code inputs are: the number of counts of the galaxy at the effective radius (N_Reff), the number of counts associated with the SBF (${\overline{N}}_{\mathrm{input}}$, either ${\overline{N}}_{\mathrm{real}}$ or ${\overline{N}}_{\mathrm{obs}}$), the sky counts (N_Sky), the size of the point spread function in pixels (σ_PSF), the size of the image in pixels (n_pix), the effective radius in pixels (R_eff), the Sersic index (n) and the applied annular mask (r₁,r₂) in pixels too. The code also asks for the range of frequencies (k_fit, i, k_fit, f) where the fitting will be performed. We recommend checking the azimuthally averaged power spectrum image to ensure that the experiment has been appropriately set and that the selected fitting range covers the frequencies of interest. Finally, the code asks for the number of Monte Carlo simulations to be done, where the randomness comes from the stellar population luminosity fluctuation and the instrumental noise. The input parameters of our reference galaxy (Sect. 2.1) can be found in Table 1. The particular values of the reference galaxy shown in this that table, together with an annular mask of (r₁ = 80; r₂ = 332) px and a range of fitting frequencies of (k_fit, i = 75, k_fit, f = 400) px⁻¹, should replicate Fig. 4 as a validation test. The code returns estimations for the uncertainties: the accuracy (the 90% percentile of the relative error, ε_rel, 90%), the fitting quality (the 90% percentile of the relative standard deviation of the fitting, σ_fit, 90% and the precision (the relative 90% width of the ${\overline{N}}_{\mathrm{fit}}$ distribution, Δ_90%). Additionally, it provides the mean value of the distribution of code fitted SBFs, $⟨{\overline{N}}_{\mathrm{fit}}⟩$, as well as a warning if our third criterion from Eq. (25) is not fulfilled, indicating a possible bias in the derived SBF.

We acknowledge that the programme is a very rough estimator of the SBF uncertainty. Nevertheless, as an ideal laboratory, it aims at providing a conservative benchmark for assessing the uncertainty, whilst real observations will likely yield less favourable results. We encourage users to adapt and modify the code in a way that accomplishes their science cases.

7. Recommendations for the SBF measurement

Having summarised the work, we now propose a set of suggestions for future SBF calculations:

– The recommended conditions for a proper SBF retrieval are as follows. As we find a strong reliance on the number of pixels within the mask, when applying a mask, the largest possible number of non-masked pixels should be kept (as long as they are adequate pixels for the measurement). Aside from the typical SBF measurements in old and metal-rich elliptical galaxies, we recommend selecting bright systems (avoiding faint sources) with high fluctuations (associated, for instance, with young or metal-poor populations, or both, Vazdekis et al. 2020). One should prepare an observation with a small PSF, a low sky background, and a long exposure time. Having the same number of pixels within the mask, the size and light profile of the galaxy do not appear to affect the result drastically, although the size of the galaxy limits the number of pixels with enough light to study.

– Before performing an observation, we can foresee the potential reliability of its SBF measurement, considering representative parameters of the target galaxy and applying the modelling described in this work. As explained in Sect. 6, case (i), this provides a preliminary estimation of the uncertainty due to the stochasticity of the system under ideal conditions. If this initial test is faulty, the observation retrieves even worse errors.

– When fitting the SBF, we strongly recommend modelling the instrumental noise (P₁ in the traditional notation of the SBF literature), in order to leave the SBF as a sole parameter to fit. It can be done using a similar procedure to the one of Mei et al. (2005) or as we proposed in Eq. (22): one can begin by modelling the sky, the mean galaxy, and the PSF for the observation; then, the sky can be added to the mean model convolved with the PSF (in this work referred as Gal_{mean, fluc, sky}); next, the Poissonian randomisation can be introduced to this composite image; lastly, the necessary noise image can be obtained by subtracting the randomised Gal_{mean, fluc, sky} from its original self. We note that to obtain the required noise image it is necessary to adapt this approach to the particularities of each observation. Additionally, a trustworthy estimation of the mean model, the PSF, and, especially, the sky are vital for the SBF derivation. An offset in these parameters leads to a bias in the SBF and high relative errors. In this situation the relative standard deviation of the fitting or the precision are not drastically influenced by an offset. Therefore, in an observation, we can find adequate results for these two uncertainties, while being unaware of how close the fitted SBF is to the real SBF value, that is, the accuracy and a possible bias. This effect is particularly sensitive to small offsets in the sky. A way of preventing this issue is carefully checking the image of the azimuthally averaged power spectrum used for the fitting, in which any inconsistency should be noticeable and provides an indication about any possible offset.

– When performing the SBF fitting in Fourier space, it is necessary to select a pair of limits (k_fit, i, k_fit, f) where the fitting is performed. Among different observations, changes in the size of the analysed image, the scale plate, or in the PSF require this range of frequencies to be adapted. At the moment, there is no universally applicable way of choosing the boundaries of the fitting, but one should manually avoid the lowest frequencies, which have larger uncertainties, and avoid the highest frequencies since the power spectrum profile becomes constant and only provides redundant information about the instrumental noise.

– In a traditional SBF derivation procedure, pixels with a flux similar or lower than ten times the SBF ($N({\mathrm{Gal}}_{\mathrm{mock}\mathrm{mask}})-N({\mathrm{Sky}}_{\mathrm{mask}})\ge 10\times \overline{N}$) do not ensure a correct modelling when retrieving the fluctuation. After estimating the SBF, we suggest checking if all the pixels used in the measurement fulfil this condition. Deriving SBFs from faint sources, such as dwarf galaxies or low brightness regions of the galaxy, might be unreliable. This is because the integrated light among pixels does not follow a Gaussian distribution and the estimate of the mean galaxy image can be biased, and hence the SBF measurement as well (Cerviño & Luridiana 2006; Cerviño et al. 2008). In every pixel of the reference ‘mean image’ (retrieved from the observation as an average of nearby pixels), its value must be a proper representation of the stellar population through the mean of the luminosity distribution.

– Besides the traditional estimation of the SBF uncertainty on an observation, we also suggest computing the error presented in the current work. Specifically, we assessed the inherent stochastic nature of the system due to fluctuations in the luminosity of the stellar population and instrumental noise. To our knowledge, this aspect is not estimated in previous literature. In this case, the actual ${\overline{N}}_{\mathrm{real}}$ is unknown, and therefore we cannot evaluate the accuracy (ε_rel, 90%). However, we can estimate the relative standard deviation of the fitting and the precision of the measurement. Following what is detailed in this paper, we can simulate a Monte Carlo set of mock galaxies (taking ${\overline{N}}_{\mathrm{input}}={\overline{N}}_{\mathrm{obs}}$), and then perform the fitting to calculate σ_fit, 90% and Δ_90% with respect to the observationally measured ${\overline{N}}_{\mathrm{obs}}$. This is illustratively explained in the flowchart of Fig. 14 and its description is provided in Sect. 6, case (ii). In an observation, the intrinsic uncertainty due to the undetermined PSF, the sky estimate, the calibration, or any other source of variability can be complemented with our proposed σ_fit, 90% and Δ_90%. We remind the reader that the relative standard deviation of the fitting (σ_fit, 90%) is a less restrictive parameter than the precision of the fitting (Δ_90%), so we recommend selecting Δ_90% as an uncertainty indicator.

In conclusion, retrieving reliable SBFs requires careful treatment of the measurement and the uncertainty. In order to exploit the potential of the SBFs fully, new studies with high-quality data and proper observational conditions are badly needed.

Acknowledgments

Special thanks to Javier Sánchez Sierras, Isaac Alonso Asensio, Martín Manuel Gómez Míguez and Ignacio Trujillo Cabrera for their invaluable advise and useful discussions. AV and PRB acknowledge support from grant PID2019-107427GB-C32 from the Spanish Ministry of Science, Innovation and Universities MCIU. Also, AV acknowledges support from grant PID2021-123313NA-I00 from the Spanish Ministry of Science, Innovation and Universities MCIU. This work has also been supported through the IAC project TRACES, which is partially supported through the state budget and the regional budget of the Consejería de Economía, Industria, Comercio y Conocimiento of the Canary Islands Autonomous Community. MC has been funded by the project PID2019-107408GB-C41 and PID2022-136598NB-C33 by the Spanish Ministry of Science and Innovation/State Agency of Research MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe”. This research has made use of the SIMBAD data base, operated at CDS, Strasbourg, France. This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency. Based on observations obtained with MegaPrime/MegaCam, a joint project of CFHT and CEA/IRFU, at the Canada-France-Hawaii Telescope (CFHT) which is operated by the National Research Council (NRC) of Canada, the Institut National des Sciences de l’Univers of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. This work is based in part on data products produced at Terapix available at the Canadian Astronomy Data Centre as part of the Canada-France-Hawaii Telescope Legacy Survey, a collaborative project of NRC and CNRS. This paper made use of the IAC Supercomputing facility HTCondor (http://research.cs.wisc.edu/htcondor/), partly financed by the Ministry of Economy and Competitiveness with FEDER funds, code IACA13-3E-2493. We acknowledge the help from the members of the SIE team of IAC.

References

Appendix A: SBF derivation considering a mask

In this appendix we show the derivation of the SBF considering a mask, which is a generalisation of the simpler case without any mask (see Sect. 2.3). We start from the definition of the mock galaxy, as in Eq. (13). Then, multiplying Eq. (13) with a mask image we obtain Eq. (20), which reads:  Gal_mock ⋅ Mask = ((Gal_mean+Gal_fluc+Sky)⊗PSF+R) ⋅ Mask, where Sky ⊗ PSF = Sky, as Sky is a constant image. In order to have all the images with the same size, the PSF image is resized by convolving it with a n_pix × n_pix image, where all the pixels are zero value but the central one, with a value of one.

Next, we rearrange the terms in the previous equation so we can isolate the fluctuation and noise terms:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill ({\mathrm{Gal}}_{\mathrm{mock}}-\mathrm{Sky}-({\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}))·\mathrm{Mask}=\\ \hfill ({\mathrm{Gal}}_{\mathrm{fluc}}\otimes \mathrm{PSF})·\mathrm{Mask}+R·\mathrm{Mask}.\end{array}\end{array}$$(A.1)

Now, we divide by the square root of the mean model, $\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}$. Note that the SBF is obtained from the study of the distribution of the square of the fluctuations around the mean (μ) divided by the mean, that is, (x − μ)²/μ. So we build an image of $(x-\mu )/\sqrt{\mu }$ which will be squared when the power spectrum is applied:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \frac{({\mathrm{Gal}}_{\mathrm{mock}}-\mathrm{Sky}-({\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}))·\mathrm{Mask}}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}=\\ \hfill \frac{({\mathrm{Gal}}_{\mathrm{fluc}}\otimes \mathrm{PSF})·\mathrm{Mask}+R·\mathrm{Mask}}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}.\end{array}\end{array}$$(A.2)

We acknowledge that the SBF definition (the ratio of the variance and the mean images) does not directly mention the PSF. However, for consistency, the definition implies that the term in the denominator must correspond to the one in the numerator, Gal_mean ⊗ PSF, instead of just dividing by Gal_mean.

We denote as Gal_{mock fluc mask} the first term of Eq. (A.2). Then, in order to disentangle the convolved PSF and the fluctuation contribution (Gal_fluc) we shift to the Fourier space. We do so by applying the power spectrum (PS(f) = |𝔉(f)|²):

$$\begin{array}{cc}& {\left|\mathfrak{F}\left({\mathrm{Gal}}_{\mathrm{m}\mathrm{o}\mathrm{c}\mathrm{k}\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{c}\mathrm{m}\mathrm{a}\mathrm{s}\mathrm{k}}\right)\right|}^2=\hfill \\ \hfill & {\left|\mathfrak{F}\left(\frac{({\mathrm{Gal}}_{\mathrm{fluc}}\otimes \mathrm{PSF})·\mathrm{Mask}+R·\mathrm{Mask}}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right)\right|}^2.\hfill \end{array}$$(A.3)

Developing the previous equation by expanding the right term as the square of complex numbers (see step between Eqs. 16 and 17), we find:

$$\begin{array}{cc}& {\left|\mathfrak{F}\left({\mathrm{Gal}}_{\mathrm{m}\mathrm{o}\mathrm{c}\mathrm{k}\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{c}\mathrm{m}\mathrm{a}\mathrm{s}\mathrm{k}}\right)\right|}^2={\left|\mathfrak{F}\left(\frac{({\mathrm{Gal}}_{\mathrm{fluc}}\otimes \mathrm{PSF})·\mathrm{Mask}}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right)\right|}^2+\hfill \\ \hfill & {\left|\mathfrak{F}\left(\frac{R·\mathrm{Mask}}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right)\right|}^2+\hfill \\ \hfill & \mathfrak{F}{\left(\frac{({\mathrm{Gal}}_{\mathrm{fluc}}\otimes \mathrm{PSF})·\mathrm{Mask}}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right)}^{†}·\mathfrak{F}\left(\frac{R·\mathrm{Mask}}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right)+\hfill \\ \hfill & \mathfrak{F}\left(\frac{({\mathrm{Gal}}_{\mathrm{fluc}}\otimes \mathrm{PSF})·\mathrm{Mask}}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right)·\mathfrak{F}{\left(\frac{R·\mathrm{Mask}}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right)}^{†}.\hfill \end{array}$$(A.4)

The convolution theorem, 𝔉(f ⊗ g) = 𝔉(f)⋅𝔉(g) and 𝔉(f ⋅ g) = 𝔉(f)⊗𝔉(g), allows Gal_fluc, the PSF contribution, and the mask to be separated. Thus, 𝔉((f⊗g)⋅m) transforms into (𝔉(f)⋅𝔉(g))⊗𝔉(m). So, we can rewrite Eq. (A.4) as:

$$\begin{array}{cc}& {\left|\mathfrak{F}\left({\mathrm{Gal}}_{\mathrm{m}\mathrm{o}\mathrm{c}\mathrm{k}\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{c}\mathrm{m}\mathrm{a}\mathrm{s}\mathrm{k}}\right)\right|}^2=\hfill \\ \hfill & ({\left|\mathfrak{F}\left(\frac{{\mathrm{Gal}}_{\mathrm{fluc}}}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right)\right|}^2·{\left|\mathfrak{F}\left(\mathrm{PSF}\right)\right|}^2)\otimes {\left|\mathfrak{F}\left(\mathrm{Mask}\right)\right|}^2\hfill \\ \hfill & +{\left|\mathfrak{F}\left(\frac{R·\mathrm{Mask}}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right)\right|}^2\hfill \\ \hfill & +(\mathfrak{F}{\left(\frac{{\mathrm{Gal}}_{\mathrm{fluc}}}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right)}^{†}·\mathfrak{F}{\left(\mathrm{PSF}\right)}^{†})\otimes \hfill \\ \hfill & \mathfrak{F}{\left(\mathrm{Mask}\right)}^{†}·\mathfrak{F}\left(\frac{R·\mathrm{Mask}}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right)\hfill \\ \hfill & +(\mathfrak{F}\left(\frac{{\mathrm{Gal}}_{\mathrm{fluc}}}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right)·\mathfrak{F}\left(\mathrm{PSF}\right))\otimes \hfill \\ \hfill & \mathfrak{F}\left(\mathrm{Mask}\right)·\mathfrak{F}{\left(\frac{R·\mathrm{Mask}}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right)}^{†}.\hfill \end{array}$$(A.5)

For the first right term of Eq. (A.5) we must calculate first the product $(|\mathfrak{F}({\mathrm{Gal}}_{\mathrm{fluc}}/\sqrt{{\mathrm{Gal}}_{\mathrm{mean}})\otimes \mathrm{PSF}}){|}^2·|\mathfrak{F}(\mathrm{PSF}){|}^2)$, and then convolve the result with |𝔉(Mask)|². If the fluctuation term is constant, it can exit the operation as a multiplicative factor. In Sect. 2.3 we show that $(|\mathfrak{F}({\mathrm{Gal}}_{\mathrm{fluc}}/\sqrt{{\mathrm{Gal}}_{\mathrm{mean}})\otimes \mathrm{PSF}}){|}^2$ can be represented by a constant number, at least in the simulations of this work. Therefore, |𝔉(PSF)|² ⊗ |𝔉(Mask)|² can be operated first, as traditionally stated in the literature. The same idea applies to the Fourier transforms in the crossed term. Eventually, we rewrite the power spectrum as ’PS’ and an azimuthal average is applied to Eq. (A.5) (note that the azimuthal average of sums is equal to the sum of azimuthal averages). Hence, radial profiles are obtained, denoted with the subindex ’r’:

$$\begin{array}{cc}& PS{({\mathrm{Gal}}_{\mathrm{m}\mathrm{o}\mathrm{c}\mathrm{k}\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{c}\mathrm{m}\mathrm{a}\mathrm{s}\mathrm{k}})}_{\mathrm{r}}=\overline{N}·{(PS(\mathrm{PSF})\otimes PS(\mathrm{Mask}))}_{\mathrm{r}}+\hfill \\ \hfill & PS{\left(\frac{R·\mathrm{Mask}}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right)}_{\mathrm{r}}+\hfill \\ \hfill & \sqrt{{\overline{N}}^{†}}·{(\mathfrak{F}{\left(\mathrm{PSF}\right)}^{†}\otimes \mathfrak{F}{\left(\mathrm{Mask}\right)}^{†}·\mathfrak{F}\left(\frac{R·\mathrm{Mask}}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right))}_{\mathrm{r}}+\hfill \\ \hfill & \sqrt{\overline{N}}·{(\mathfrak{F}\left(\mathrm{PSF}\right)\otimes \mathfrak{F}\left(\mathrm{Mask}\right)·\mathfrak{F}{\left(\frac{R·\mathrm{Mask}}{\sqrt{{\mathrm{Gal}}_{\mathrm{mean}}\otimes \mathrm{PSF}}}\right)}^{†})}_{\mathrm{r}}.\hfill \end{array}$$(A.6)