Comparing PyMorph and SDSS photometry. II. The differences are more than semantics and are not dominated by intracluster light

TL;DR

This paper compares PyMorph and SDSS photometry, revealing that differences are mainly due to model fitting choices and not intracluster light, impacting galaxy mass estimates and related cosmological models.

Contribution

It demonstrates that photometric differences are driven by model fitting, not intracluster light, and discusses implications for galaxy mass functions and halo-galaxy relations.

Findings

SDSS underestimates galaxy brightness for luminous galaxies

Model choice impacts photometric measurements more than intracluster light

Previous mass functions were biased due to SDSS underestimates

Abstract

The Sloan Digital Sky Survey pipeline photometry underestimates the brightnesses of the most luminous galaxies. This is mainly because (i) the SDSS overestimates the sky background and (ii) single or two-component Sersic-based models better fit the surface brightness profile of galaxies, especially at high luminosities, than does the de Vaucouleurs model used by the SDSS pipeline. We use the PyMorph photometric reductions to isolate effect (ii) and show that it is the same in the full sample as in small group environments, and for satellites in the most massive clusters as well. None of these are expected to be significantly affected by intracluster light (ICL). We only see an additional effect for centrals in the most massive halos, but we argue that even this is not dominated by ICL. Hence, for the vast majority of galaxies, the differences between PyMorph and SDSS pipeline photometry cannot be ascribed to the semantics of whether or not one includes the ICL when describing the stellar mass of massive galaxies. Rather, they likely reflect differences in star formation or assembly histories. Failure to account for the SDSS underestimate has significantly biased most previous estimates of the SDSS luminosity and stellar mass functions, and therefore Halo Model estimates of the z ~ 0.1 relation between the mass of a halo and that of the galaxy at its center. We also show that when one studies correlations, at fixed group mass, with a quantity which was not used to define the groups, then selection effects appear. We show why such effects arise, and should not be mistaken for physical effects.

Figures

Equations

$$\mathtt{Model}-\mathtt{SerExp}=(\mathtt{deV}-\mathtt{SerExp})-(\mathtt{deV}-\mathtt{Model}),$$

$$n(m,H)=n(m){\int }_H\mathrm{d}hp(h|m),$$

$$\displaystyle\bigl{\langle}s-m|m,H\bigr{\rangle}$$

$$\bigl{\langle}s-m|m,H\bigr{\rangle}=\int{\rm d}s\,p(s|m)\,(s-m).$$

$$\displaystyle\bigl{\langle}s-m|s,H\bigr{\rangle}$$

$$=\frac{n(s)}{n(s,H)}{\int }_H\mathrm{d}hp(h|s)\int \mathrm{d}mp(m|s,h)(s-m)$$

$$n(s,H)={\int }_H\mathrm{d}h\int \mathrm{d}mn(m,s,h)=n(s){\int }_H\mathrm{d}hp(h|s).$$

$$\bigl{\langle}s-m|s,H\bigr{\rangle}=s-\langle m|s\rangle-\frac{n(s)\,C}{n(s,H)}\int_{H}{\rm d}h\,p(h|s)\,[h-\langle h|s\rangle].$$

Full text

1 Introduction

The observed magnitudes reported by the Sloan Digital Sky Survey are underestimated especially at the highest luminosities (e.g. Bernardi et al. 2007; von der Linden et al. 2007; Bernardi et al. 2013; Meert et al. 2015; D’Souza et al. 2015; Bernardi et al. 2017). In a companion paper [Fischer et al. 2017], we show that this is due to a combination of sky background and model fitting effects: (i) the SDSS overestimates the sky background (Blanton et al. 2005; Bernardi et al. 2007; Hyde & Bernardi 2009; Blanton et al. 2011; Meert et al. 2015), and (ii) single- or two-component Sersic model based estimates of galaxy luminosities are more reliable than estimates based on single exponential, single de Vaucouleurs, or a linear combination of the best fitting exponential and de Vaucouleurs models (the so-called cModel magnitudes) used by the SDSS pipeline (e.g. Bernardi et al. 2010; Bernardi et al. 2014). Fischer et al. show that PyMorph sky estimates are fainter than those of the SDSS DR7 or DR9 pipelines, but are in excellent agreement with the estimates of Blanton et al. (2011). The PyMorph Sersic-based estimates [Meert et al. 2013, Meert et al. 2015, Meert et al. 2016] are more reliable and return more light than do SDSS estimates, and the difference is most pronounced at the highest luminosities. This can have a dramatic impact on the estimated $z\sim 0.1$ luminosity and stellar mass functions [Bernardi et al. 2013, Bernardi et al. 2016, Bernardi et al. 2017], and hence on Halo Model [Cooray & Sheth 2002] based estimates of the relation between stellar and halo mass (Kravtsov et al. 2014; Shankar et al. 2014). They also impact models for the formation of the massive galaxies because they affect estimates of the mass scale on which galaxy scaling relations show curvature [Bernardi et al. 2011].

The most luminous galaxies reside in or at the centers of clusters. Since clusters are known to possess intercluster light (hereafter ICL), it is natural to ask if the PyMorph Sersic-based photometry, with its improved sky estimates, is brighter than the SDSS measurements (i.e. Model or cModel magnitudes) because it includes more of this ICL component. Previous work, based on stacked images of central galaxies of massive clusters, suggests that the ICL only contributes substantially to the surface brightness profile on scales larger than 50 kpc [Zibetti et al. 2005], where the surface brightness has dropped below about 27 mags arcsec*-2*. A stacking analysis of LRGs suggests this scale may be even larger for lower mass groups [Tal & van Dokkum 2011]. For the vast majority of central galaxies in the SDSS, this corresponds to a scale where the surface brightness has dropped to less than 1% of the sky brightness, making it extremely difficult to detect in individual images. While this means that it is very unlikely that the ICL contributes substantially to the PyMorph reductions, in what follows, we provide additional evidence that the PyMorph-SDSS difference is unlikely to be due to ICL.

Our logic is simple: Since the ICL should be centered on the cluster center, and is expected to be fainter at larger cluster-centric distances, it is reasonable to suppose that it will affect the photometry of the central galaxy more than the satellites. The ICL is also expected to be fainter around central galaxies in less massive groups. Therefore, we study if the difference between PyMorph and SDSS magnitudes depends on whether the galaxy is a central or a satellite, as well as if it depends on cluster mass. Section 2 describes our galaxy sample and shows why the ICL is unlikely to play a major role in the images of individual galaxies. Section 3 describes the group catalogs we use for separating centrals from satellites. Section 4 presents the main results of our analysis and Section 5 summarizes our findings. While the main text presents the results using the Sersic-Exponential photometry, Appendix A shows a similar analysis for the single-component Sersic fits since these are much more common in the literature. Finally, Appendix B discusses a pernicious selection effect which arises when working with group catalogs, and which we were careful to avoid.

When necessary, we assume a spatially flat background cosmology with parameters $(\Omega_{m},\Omega_{\Lambda})=(0.3,0.7)$, and a Hubble constant at the present time of $H_{0}=70$ km s*-1Mpc-1*.

2 Comparison of PyMorph and SDSS luminosities: The main galaxy sample

2.1 The parent galaxy sample

The analysis which follows is based on the Main Galaxy sample of the Ninth Data Release of the Sloan Digital Sky Survey (hereafter SDSS DR9; Aihara et al. 2011); the sample is limited to a Petrosian $r$-band apparent magnitude of $m_{r{\rm Pet}}\leq 17.7$. (Our conclusions are unchanged if we use DR7 values instead of DR9.) To about $\sim 660,000$ of these (those at $z\leq 0.25$), Huertas Company et al. (2011) have assigned Bayesian Automated Morphological Classification weights which represent the probabilities that the galaxy is Elliptical, S0, Sab, or Scd. The vast majority of the most luminous galaxies are Es or S0s, so, in what follows, we define the E+S0 weight to equal $p$(E+S0) $\equiv p$(E) + $p$(S0).

For every DR9 galaxy, the SDSS pipeline provides Model magnitudes, which are based on separately fitting an exponential and a de Vaucouleurs form to the surface brightness profile and choosing the value returned by the model which fits best. Hence, for essentially all E+S0s, the Model magnitude is that from the de Vaucouleurs fit. (Our conclusions are unchanged if we use cModel magnitudes).

We will compare these Model magnitudes with single-component Sersic and two-component Sersic-Exponential fits (hereafter Ser and SerExp) to these same objects provided by Meert et al. (2015). These Ser and SerExp fits are returned by the PyMorph algorithm which is described and tested in Meert et al. (2013, 2015, 2016) and used by Bernardi et al. (2013, 2014, 2016, 2017). Although those tests were based on DR7 objects, Fischer et al. (2017) show that the difference between DR7 and DR9 is negligible for PyMorph. I.e., the parameters provided by Meert et al. (2015) for DR7 can also be used for DR9. The two-component SerExp fits are the most accurate (i.e. least biased) of the PyMorph outputs [Meert et al. 2013, Bernardi et al. 2014]. For this reason, we use these in the main text. However, single-component fits are much more common in the literature; we show our analysis of Ser photometry in Appendix A.

We also present our results as a function of stellar mass $M_{*}$. For a given luminosity (e.g., Model, or SerExp), we obtain $M_{*}$ by multiplying $L$ by the dust-free $M_{*}/L$ estimates of Mendel et al. (2014) assuming a Chabrier (2003) IMF. Bernardi et al. (2017) describe a number of other reasonable choices of $M_{*}/L$. The results which follow are robust to changes in this choice.

2.2 Example profiles and fits

Before showing results of a statistical analysis, it is useful to look at a few representative images. This helps build intuition for what it is that we will be averaging. Our aim here is to directly address what previous work based on stacked images imply for our analysis of individual images. We are particularly interested in the effects of the ICL. If ICL is defined as being an excess above the best-fitting deV profile, then it is only expected to contribute significantly on scales larger than 50 kpc [Zibetti et al. 2005]. On the other hand, if ICL is a departure from the best-fitting Ser profile, then the relevant scale may be about a factor of two larger [Tal & van Dokkum 2011].

Figures 1–3 show the surface brightness profiles of three luminous ($M_{r}\leq-23$) galaxies in the SDSS main galaxy sample, which are at a range of distances ($z=0.06$, 0.15 and 0.25). As our primary interest here is in the role of the ICL, all three objects are central galaxies in the group catalogs we use extensively later in this paper. In all cases, the top panels show the image, the middle panels show the best 2D PyMorph SerExp fit and residuals from it, and the bottom panels show 1D angular averages of the profile, the fit and residuals (the latter are computed as fit-data). The 2D residuals show the image including sky whereas those in the bottom the sky has been subtracted. These bottom panels are for illustration only, as the fitting was done in 2D with more care taken to remove light associated with the other objects in the field. The legends on either side provide a wealth of information about the parameters of the fit (see Meert et al. 2015 for details). In all cases, the image cutouts are labeled in units of arcsecs, and the color scale is in mags/arcsec2.

Figure 1 shows a galaxy which was selected because of its large apparent brightness ($m_{r}\sim 14.3$), so that the dynamic range between the half light radius $R_{hl}$ and the scale $R_{1\%{\rm sky}}$ on which the profile has dropped to 1% of sky, is large. The SerExp absolute magnitude of this object is $\sim-23$ and $R_{hl}\approx 13$kpc. For this object, $R_{1\%{\rm sky}}\approx 50$ kpc; note that there is no obvious feature in the profile shape on scales smaller than this. Indeed, a single Ser fit to this object has $n_{\rm Ser}=5$ and is not very different from the best fitting SerExp profile shown by the solid curve in the bottom left panel. The dashed curve in the same panel shows the bulge component of the SerExp fit (a single Sersic profile with $n\approx 4$, meaning that this bulge component has a deV profile). Beyond about 10 arcsec from the center – i.e., beyond about $R_{hl}$ – the second component (dotted curve) is clearly necessary. Note in particular that this second component, which describes light in excess of a deV profile, is necessary on scales which are much smaller than the 50 kpc associated with the ICL.

Figure 2 shows a more luminous ($M_{r}=-23.65$) and distant ($z=0.15$) central galaxy with $m_{r}\approx 16$. For this object too, the bulge component is not too different from a deV profile, but the second component is necessary even on scales as small as 15 kpc. A single Ser profile with $n=5.2$ also provides a good fit on all scales smaller than $R_{1\%{\rm sky}}$.

Figure 3 shows another object of similar luminosity, but at even higher redshift, so its apparent magnitude is significantly fainter ($m_{r}\sim 17.4$). Comparison with the bottom left panels of the previous figures shows the reduced dynamic range which makes it more difficult to detect departures from a single Ser fit which, in this case, has $n_{\rm Ser}=7.9$. As for the previous two examples, the bulge is closer to deV, and the need for a second component is already evident on scales of order 20 kpc.

Thus, to the extent that these galaxies are representative of all central galaxies, these figures make two points. First, departures from a deV profile are almost always detected with high significance (see also e.g. Gonzalez et al. 2005; Bernardi et al. 2007). Moreover, a pure deV profile becomes a poor fit on much smaller scales than is reasonable to associate with the ICL. Second, if the ICL is associated with departures from a single Ser (rather than deV) fit, then if this occurs, it is at surface brightnesses which are too faint to be seen in individual images. Therefore, the departures from a pure deV profile, which PyMorph detects in its Ser or SerExp fits and which are the subject of this paper, likely reflect differences in star formation or assembly histories, rather than ICL.

Having illustrated that differences from deV photometry are common, we now turn to a statistical analysis of these differences.

2.3 Comparison of SDSS and PyMorph pipelines

Figure 4 shows a comparison of SDSS Model and PyMorph SerExp magnitudes; a similar analysis using Ser magnitudes is shown in Appendix A. Here, as in most of the figures which follow, objects have been weighted by $p$(E+S0). This ensures that we are working with a sample for which Model = de Vaucouleurs and removes the question of how the morphological mix affects the PyMorph-SDSS comparison.

Symbols with error bars show the median and the error on it (only bins with more than 50 objects are shown), and dashed curves show the region which encloses 68% of the sample in each absolute magnitude bin. The SDSS Model magnitudes are increasingly fainter as luminosity increases. Some of the differences in Figure 4 are simply due to the fact that SDSS magnitudes are based on integrating the fitted profile to approximately $7.5\times$ the semi-major axis $a_{e}$, whereas PyMorph does not truncate. For this reason, we have shown the result of truncating the PyMorph fits as well: this makes them fainter (by less than 0.1 mags on average), but other than this shift, the overall trends with luminosity are unchanged. Therefore, truncation is not the primary reason why SDSS is fainter (see Sections 2.3 and 2.4 in Fischer et al. 2017 for more discussion of truncation).

Fischer et al. (2017) show that, once truncation has been accounted for, there are two effects which contribute to biasing SDSS low (rather than biasing PyMorph high). These are due to differences in how the background sky is estimated, and what model is fit to the surface brightness profile. Regarding the first effect, there is now general consensus that the SDSS treatment of the sky is flawed. This affects nearby galaxies (as emphasized by Blanton et al. 2011) but also high luminosity galaxies which have relatively large angular sizes and tend to be in crowded fields (Fischer et al. 2017 and references therein). PyMorph attributes less light to the sky than does the SDSS; as a result, it assigns more light to the galaxy than does the SDSS. Fischer et al. (2017) also show that PyMorph sky estimates are in excellent agreement with those of Blanton et al. (2011). About half of the bias at high luminosities arises from fitting different models (the second effect). At the high luminosity end, most galaxies are E+S0s, for which the SDSS Model magnitudes are essentially deV magnitudes, and deV fits return less light than SerExp or Ser fits.

3 Galaxy groups in the SDSS

Our goal is to check if the differences between PyMorph and Model photometry shown in Figure 4 depend on whether a galaxy is a central or a satellite. To achieve our goal, we use two group catalogs in which centrals and satellites have been identified. One of these is from Yang et al. (2007) (the DR7 version, hereafter Yang+), and the other is the redMaPPer sample of Rykoff et al. (2014). Whereas the former identifies groups spanning a wide range of masses, the latter only identifies the very most massive clusters (the redMapper algorithm is not well-suited for identifying lower mass/richness groups).

3.1 Description of group catalogs

The redMapper sample is drawn from the SDSS DR8 footprint which covers $\sim 10,000$ deg2. Groups are identified on the basis of angular positions and color. Of the $1.7\times 10^{6}$ objects in $2.6\times 10^{4}$ groups in the redMaPPer sample, Meert et al. only provide PyMorph reductions for $1.8\times 10^{4}$: these are the subset which have spectroscopic information and were in the SDSS DR7 $\sim 7,700$ deg2 footprint. These objects are in about 3,400 clusters, of which about 2,400 are at $z\leq 0.25$, and 1,200 are at $z\leq 0.2$.

In contrast, Yang+ identify $6\times 10^{5}$ galaxies in about $4.7\times 10^{5}$ groups at $z\leq 0.2$ in the SDSS DR7 footprint. Most of these are much less massive groups, of course. Yang+ defined their catalog using a complex iterative procedure that makes use of observed angular positions, redshifts and Model photometry. They estimate the halo masses $M_{\rm Halo}$ of their groups using two simple proxies for halo mass. One is based on summing up the Model luminosities of the central+satellite galaxies. As a result, there is a tight correlation between $M_{\rm Halo}$ and Model which will be important below (see also Appendix B). The other is based on summing their stellar masses, which they estimate by multiplying the Model luminosities by a simple color-based estimate of $M_{*}/L$ (see their equation 2). We use their luminosity based estimate, but the results which follow are robust to changes in this choice.

The Yang+ catalog only extends out to $z=0.2$, so our first step was to identify objects which appear in the redMaPPer catalog as well. There are 13,253 such objects of which 8,081 have $p$(E+S0)$\geq 0.7$. Figure 5 shows the distribution of these objects in the $M_{\rm Halo}-$Model magnitude plane.

Red and green symbols show objects Yang+ classify as centrals and satellites, respectively. Notice that the centrals define a rather tight $M_{\rm Halo}$-Model magnitude correlation which has a sharp boundary. This is because of the Yang+ definition of $M_{\rm Halo}$, and will be important later (see discussion related to Figures 7 and 8, and Appendix B). More important in the present context are the brown symbols, which show the centrals in redMaPPer (the redMaPPer satellites are all the points which are not brown). Below $M_{\rm Halo}\sim 10^{14}M_{\odot}$ only a small fraction of the Yang+ centrals are also centrals in redMaPPer (i.e., redMaPPer labels most Yang+ centrals as satellites – remember that this is the subsample of galaxies which is common to both group catalogs). If Yang+ are correct, then redMaPPer is wrongly linking together objects which are really in separate halos – this will compromise efforts to use redMaPPer to address assembly bias like effects. On the other hand, if redMaPPer is correct, then the Yang+ misclassifications will make centrals and satellites seem more similar than they really are. This disagreement between the two groups complicates any attempt to draw unambiguous conclusions about the central-satellite difference in lower mass halos.

In view of the central-satellite disagreement, in what follows, we work primarily with the Yang+ objects which Yang+ classified as being in halos more massive than $10^{14}M_{\odot}$. While Figure 5 suggests this is a reasonable choice, things are not completely straightforward, since only about one in three of the Yang+ centrals is in the redMaPPer catalog. We show this in Figure 6, where $N$ is the total number of objects in each catalog (i.e. not the small subset which were common to both catalogs) weighted by $p$(E+S0). (Since we only have BAC weights at $z\leq 0.25$, we set $p$(E+S0)=1 for all objects at $z>0.25$, since this is very likely to be realistic. We have checked that making cuts on concentration index or color smooths out the small plateau at $z=0.25$ in the grey histogram, but makes essentially no difference to the others.) The bottom panel only includes galaxies at $z\leq 0.25$. Remember that the Yang+ catalog only includes objects a $z<0.2$. Restricting the full sample (grey histogram) to galaxies with $z\leq 0.2$ gives the same distribution as for the full Yang+ catalog (black histogram). For completeness, the bottom panel also shows the subsample of the centrals in the redMaPPer catalog at $z<0.2$ (brown dotted line). (We show the number of objects $N$, rather than the comoving number density, since our goal is to show how many objects contribute to each bin in the Figures which follow.)

Note that, in the bottom panel, there are $3\times$ more centrals in Yang+ than in redMaPPer, even when restricting Yang+ to $M_{\rm Halo}\geq 10^{14}M_{\odot}$. This remains true if we use the Yang+SerExp catalog (described in Section 3.2) and show the counts as a function of SerExp magnitudes instead (the main difference is a shift to brighter magnitudes). We assume that the reason for this difference (i.e. the factor of $3\times$) is not the central-satellite designation, but the group richness: presumably many of the Yang+ objects are centrals of lower richness groups which failed the redMaPPer richness cut. In partial support of this, we note that the number of redMaPPer satellites is nearly the same as Yang+, despite having $3\times$ fewer centrals; this is consistent with them being more massive. In what follows, we use all the redMaPPer objects within $z\leq 0.25$. We do not include higher $z$ objects because we do not have E+S0 classifications above $z=0.25$. (We have checked, but do not show here, that including the red and/or high-concentration objects from $z\geq 0.25$ makes no difference to our results, other than to improve the statistical significance.) Moreover, Yang+ only extends to $z=0.2$, and we did not want questions of evolution being different in the two samples to complicate our results.

Finally, the magenta curves show the distribution of $\sim 5500$ E+S0 central galaxies in the MaxBCG catalog of Rykoff et al. (2012) restricted to $z\leq 0.25$ and which have PyMorph reductions from Meert et al. (2015). At the highest luminosities, their comoving density is greater than redMaPPer but smaller than the $\geq 10^{14}M_{\odot}$ Yang+ centrals, so we expect them to be intermediate in mass as well. This is our primary reason for including this sample.

3.2 Combining PyMorph with Yang+

We remarked in the previous subsection that Model photometry plays an important role in the Yang+ catalog. This raises the question of how the catalog is modified if we use PyMorph photometry instead. Although the appropriate thing to do is to re-run the algorithm, this is well-beyond the scope of the current study. Instead, we have performed the following simple but reasonable procedure.

For each Yang+ group, we assume that although changing the photometry may change the central-satellite designation within a group, it will not change group membership. We then replace the Model photometry with PyMorph values for each group member, and sum these new values to obtain a new estimate of the total group luminosity. We rank-order this quantity. Since this new ordering is different from that based on Model photometry, we reassign halo masses based on this new rank ordering. We also define the most luminous galaxy in a group as the central. Doing this separately for SerExp and Ser photometry gives us Yang+SerExp and Yang+Ser catalogs. For both catalogs, the set of halo masses associated with the original Yang+ catalog is unchanged (by construction), but the mapping between $M_{\rm Halo}$ and central galaxy luminosity is modified.

The main effect to a plot like Figure 5 is to shift the points slightly to brighter magnitudes – we explore the consequences of this for Halo Model like analyses elsewhere – but the motivation for splitting the sample at $M_{\rm Halo}\sim 10^{14}M_{\odot}$ remains. Since the shifts are small, we have not included a plot showing, e.g., $M_{\rm Halo}$ versus SerExp magnitude in the Yang+SerExp catalog. However, defining the Yang+SerExp and Yang+Ser catalogs is important for what follows (see also Appendix B).

4 Comparison of centrals and satellites over a range of galaxy and halo masses

4.1 Centrals and satellites in poor groups

Figure 7 shows the difference between SerExp and Model photometry, as a function of Model magnitude. (The SerExp magnitudes are based on truncating the fits at $7.5a_{e}$ as described previously.) Grey symbols show the full E+S0 sample; red and green curves show the subset of galaxies identified as being centrals (red) and satellites (green) in halos having masses between $10^{13}M_{\odot}$ and $10^{14}M_{\odot}$ in the original Yang+ catalog. They are both very similar to the grey curve, which shows the median trend for all E+S0 galaxies (from top panel of Figure 4). Hence, either the Yang+ central/satellite designations are completely random, or the difference between Model and SerExp magnitudes does not depend on whether a galaxy is a central or a satellite. If the latter, then either these groups are too low mass to have a significant ICL component, or what ICL is present does not play an important role in the PyMorph-SDSS difference. In either case, the PyMorph-SDSS difference, at least for galaxies in groups less massive than $10^{14}M_{\odot}$, is real – it is not just semantics.

Although Figure 4 indicates that the difference depends strongly on which quantity is used as reference, we show in Appendix B that simply plotting versus SerExp magnitudes (instead of Model) is biased by a selection effect which arises because $M_{\rm Halo}$ in the Yang+ group catalog is tightly correlated with Model magnitude (c.f. Figure 5), but less so with SerExp. For this reason, Figure 8 shows the difference in the Yang+SerExp catalog, where $M_{\rm Halo}$ is correlated with SerEXp rather than Model. (For a given group, the central-satellite designation in Yang+SerExp may be different from that in Yang+. Defining the ‘central’ as the ‘brightest’ is important: using the Yang+ designation here produces a noticable bias, see Figure 15.) Note in particular that although the median difference (grey) is larger than in Figure 7, centrals and satellites are still remarkably similar. Hence, the agreement of the red and green lines with the grey one in Figures 7 and Figure 8 strongly suggests that centrals and satellites are similar, and that the PyMorph-SDSS difference is not just semantics.

In this context it is interesting to note that Luminous Red Galaxies (LRGs) are believed to populate groups having masses of a few times $10^{13}M_{\odot}$. This is similar to the mass scale we are considering here. For LRGs at $z\sim 0.34$, a stacking analysis shows that the characterisic scale where the ICL becomes apparent is $\sim 100$ kpc [Tal & van Dokkum 2011], where the surface brightness is $28$ mags arcsec*-2*. Accounting for surface brightness dimming between $z\sim 0.15$ and $z=0.34$ would make this $27.4$ mags arcsec*-2*. This is substantially fainter than the $\sim 26$ mag/arcsec2 which corresponds to $\sim 1\%$ of the sky value in the individual $r$-band SDSS images of the objects which contribute to Figures 7 and 8. I.e., the individual images will show no sign of the ICL, so it is very unlikely that PyMorph’s Sersic-based fits are sensitive to it. And indeed, Figures 1 and 2 show that the PyMorph fits are rather good within $R_{1\%{\rm sky}}$. Stated differently, the ICL component in the LRG stacks shows up as excess light compared to a Sersic profile (with $n=5.5$) on scales larger than $\sim 8R_{e}$, but Figures 1 and 2 show that the individual images typically have $R_{1\%{\rm sky}}\approx 4R_{e}$, which is why they show no sign of the ICL.

Finally, note that a de Vaucouleurs profile is not a good fit to the regions within $\sim 8R_{e}$, neither for LRGs nor for the galaxies we are considering here. Thus, if one thinks of a Sersic profile as describing light in excess of a de Vaucouleurs profile, then this excess light is not ICL. The ICL shows up as an additional departure from a Sersic profile beyond $8R_{e}$. This is consistent with our assertion above: at least for galaxies in groups less massive than $10^{14}M_{\odot}$, the PyMorph-SDSS difference is not due to the ICL.

Now, the PyMorph-SDSS difference for the full sample of galaxies (grey curve) is similar to that for these centrals and satellites. Since the vast majority of these are in lower mass halos, ICL effects do not play a role for the vast majority of objects which contribute to the grey curve. These are the objects which Bernardi et al. (2013, 2016, 2017) used to estimate the $z\sim 0.1$ luminosity and stellar mass functions. Hence, the difference between these mass functions and those based on SDSS pipeline photometry are real: it is incorrect to attribute it to the semantics of whether or not one includes the ICL when estimating the light from a galaxy.

4.2 Centrals and satellites in massive clusters

We now consider galaxies in more massive groups. Figure 9 shows the difference between SerExp and Model photometry, as a function of Model (left) and SerExp magnitude (right). The grey curves in the two panels are for the full sample and are taken from Figure 4 as before. Red and green symbols and curves show the subset of galaxies which Yang+ identify as being centrals and satellites in halos more massive than $10^{14}M_{\odot}$ – these are shown only on the left panel due to the selection effect just mentioned above (and discussed in Appendix B). In the right panel, the red and green symbols and curves show the corresponding subset of galaxies in our Yang+SerExp sample. (I.e., they are in halos more massive than $10^{14}M_{\odot}$, but because the assignment of halos masses – and in some cases the central satellite designation – has changed, they are not exactly the same objects as in the panel on the left.) Brown and blue symbols and curves in both panels show centrals and satellites in the redMaPPer catalog, and magenta shows the MaxBCG centrals.

In both panels, the satellites are in good agreement with the grey: satellites are similar to the average over the full population, even at $M_{r}<-23$. In the panel on the right, centrals and satellites are also in good agreement (only ther redMapper centrals are slightly offset), even though the median differences between SDSS and SerExp are larger compared to the panel on the left. Only in the panel on the left do the centrals tend to be slightly more luminous, with the difference increasing with group mass (recall that the redMaPPer groups are more massive than Yang+). At $M_{r}<-23$, where the grey curve indicates the average difference is 0.15 mag, the redMaPPer centrals show an additional 0.08 mag difference. Clearly, the centrals of the most massive clusters, which are less than 30% of the rarest most luminous objects (Figure 6), are different from the vast majority of the galaxy population.

However, even for these centrals, it is not obvious that one can attribute the additional offset entirely to ICL-related effects for the reasons given in Section 2.2. As Figures 1–3 illustrate (Figure 3 is a redMaPPeR central), the PyMorph fits are usually accurate out to $\sim 1\%$ of sky. If there is a difference beyond this scale, it is in the sense that the ICL will be yet another addition to the SerExp estimate.

4.3 Differences in sky and profile shape

When discussing Figure 4 we noted that the median trend for the full E+S0 sample (labelled ‘All’) is due to two effects: one is due to differences in the PyMorph and SDSS sky estimates, and the other to the increased freedom which the two component SerExp profile has compared to Model, which, for E+S0s, is basically a single de Vaucouleurs profile. This raises the question of which effect dominates the central-satellite differences we see in Figure 9?

To address this, the lower set of curves in Figure 10 compare the sky estimates in the PyMorph deV, Ser, and SerExp fits. These show no significant differences from zero. This is in marked contrast to the top set of curves showing a large asymmetric scatter. These indicate that the SDSS DR7 sky estimate is brighter than PyMorph’s for the most luminous objects. Moreover, the SDSS overestimate is even brighter for centrals than for satellites. The former trend was reported by Fischer et al. (2017, and references therein) who show that the SDSS sky estimate is brighter than both PyMorph and Blanton et al. (2011), while the latter two are in excellent agreement with one another. Fischer et al. argued that it was a consequence of the fact that the SDSS does not use a large enough region from which to estimate the background sky. The issue with the SDSS sky will only be exacerbated in clusters, so the fact that it is worse for centrals than satellites is not unexpected. (The comparison here uses the SDSS DR7 rather than DR9 sky values, since the PyMorph photometry is based on the DR7 rather than DR9 flux calibrated images. However, Section 3 and Figures 9 and 10 in Fischer et al. 2017 show that, despite the changes from DR7 to DR9, the bias in the SDSS DR7 sky estimates is also present for the DR9 reductions.)

To see if this can account for all the differences shown in Figure 9 it is useful to write

[TABLE]

where deV is the result of using PyMorph to fit a (truncated) de Vaucouleurs profile to the image. The first term (i.e. deV-SerExp) isolates the effect of fitting different models whereas the second term (i.e. ${\tt deV}-{\tt Model}$) is entirely due to the difference between the PyMorph and SDSS sky estimates.

Figure 11 shows these two contributions for the various samples shown in Figure 9. Comparison with the left hand panel of Figure 9 shows that the sky accounts for a little less than half the total difference, and it affects the redMaPPer centrals the most. At $M_{r}\sim-23$, the choice of profile shape alone accounts for about 0.08 mags when averaged over all E+S0s, with the average for redMaPPer centrals being larger by an additional 0.08 mags. However, in the panel on the right, the central-satellite difference is entirely due to differences in sky.

Figure 11 makes one more point. The grey curve in the left panel shows that the ${\tt deV}-{\tt SerExp}$ difference is not monotonic with luminosity. Section 4 in Fischer et al. (2017) discusses this in more detail and argues that it suggests there are two populations at high luminosities. However, their analysis did not provide a physically motivated model for the two populations. Our analysis suggests that the two populations may correspond to centrals and satellites.

Figure 12 shows that none of our conclusions are changed if we plot the magnitude difference versus the corresponding stellar mass estimate rather than luminosity. And Figures 13 and 14 show that they are unchanged if we replace SerExp magnitudes with Ser. (For all these figures, we use Yang+ for the panels on the left, and Yang+Ser for the panels on the right.) In all cases, once the bias in the SDSS sky has been accounted for, the luminosity or stellar mass dependence of the SDSS-PyMorph difference is the same for all satellites, whatever the mass of the cluster they inhabit. Although for centrals of the most massive clusters, there is an additional effect which increases as cluster mass increases, this additional effect is substantially smaller when shown as a function of Sersic-based rather than SDSS Model magnitudes.

To explore the mass dependence further, we divided the Yang+ and Yang+SerExp samples into objects having $14\leq\log_{10}(M_{\rm halo}/M_{\odot})\leq 14.3$ and $\log_{10}(M_{\rm halo}/M_{\odot})\geq 14.3$. We saw a small difference between the two, which is consistent with the mass trend we infer from the comparison with redMaPPeR. Similarly, when redMaPPer is split into two subsamples at richness $\lambda=32$, the one with higher richness shows a slightly larger difference. In both cases, the smaller sample size means the statistical significance of the finer mass trend is smaller, so we have not shown these additional trends here.

We conclude that the analysis in this Section has shown that, for the vast majority of luminous galaxies, the difference between PyMorph and SDSS photometry is not dominated by ICL-like effects. Rather, it reflects real structural differences, and is not just a matter of semantics.

5 Discussion

The SDSS and PyMorph luminosity estimates differ (Figure 4), because the SDSS sky estimate is biased, and because the fitted models differ [Fischer et al. 2017] . Biases in the SDSS sky estimate are worst for central galaxies in the most massive halos (brown curve in Figure 10). In general, sky-related biases account for a little less than half the total difference between the SDSS and PyMorph luminosities (lower set of curves in Figures 11, 12, 13 and 14). The remaining difference (upper set of curves in Figures 11, 12, 13 and 14) is due to fitting different models to the surface brightness profile (i.e. deV, Ser or SerExp).

For the vast majority of galaxies, once biases in the SDSS sky estimate have been accounted for, the SDSS-PyMorph difference does not depend strongly on whether a galaxy is a central or a satellite. The difference, averaged for the full sample of galaxies (grey curve), is the same as that of central or satellite galaxies of less massive groups (red and green curves in Figures 7 and 8, and of satellites in more massive groups (blue curves in Figures 11 and 12). We only see an additional effect for centrals in the most massive halos (which Figure 6 shows are rare), and then only when shown as a function of SDSS rather than Sersic-based photometry. This conclusion is similar whether one uses two-component SerExp photometry or single component Ser fits (compare Figures 11 and 12 with Figures 13 and 14).

The ICL is expected to be fainter at larger cluster-centric distances or around central galaxies in smaller groups (e.g. Tal & van Dokkum 2011). Now, the difference between Sersic-based and SDSS Model magnitudes, when averaged over the full population is the same as for centrals or satellites in smaller groups (Figures 7 and 8). In fact, Figures 11 and 12 show that it is also the same for satellites in more massive groups. Therefore, we conclude that the vast majority of massive galaxies are not well fit by a simple de Vaucouleurs profile: the difference from Model magnitudes is not due to the ICL, but indicates real structural differences. Hence, the difference between PyMorph Sersic-based luminosity and stellar mass functions in Bernardi et al. (2013, 2016 and 2017) with respect to estimates based on SDSS Model magnitudes is real – it is not just semantics.

The issue is slightly more complicated for the small number of galaxies which are centrals of the most massive clusters. Even though the top set of curves in the right hand panels of Figures 11, 12, 13 and 14 are so similar, they are offset brightwards in the left panel (compare brown and grey curves). Hence, there must be differences between these objects and the vast majority of the galaxy population. However, we believe that even these differences – which affect less than 30% of the rarest most luminous objects (Figure 6) – should not be attributed to the ICL entirely. This is because the surface brightness where the ICL becomes apparent is fainter than $\sim 27$ mag/arcsec2 (Zibetti et al. 2005; based not on individual galaxies, but a stacking analysis of $\sim 600$ brightest cluster galaxies). This is fainter than the $\sim 26$ mag/arcsec2 associated with $\sim 1\%$ of the background sky in individual $r-$band SDSS images, so it is very unlikely that the ICL plays a role when fitting to individual images (Figures 1–3 and related discussion).

Moreover, we note that the right hand panels of Figure 12 and especially Figure 14 show a dramatic change in slope around $2\times 10^{11}M_{\odot}$. This is the same mass scale where a number of other scaling relations also change (Bernardi et al. 2011, 2014), and where the galaxy population becomes dominated by slow rotators (Cappellari et al. 2013). These other features are thought to indicate a change in the typical assembly mechanism of the population. While this change may also be related to the build-up of the ICL, the features are not caused by it.

For all these reasons, we believe that PyMorph SerExp photometry represents a significant improvement to SDSS-Model photometry. For the vast majority of massive galaxies, the differences between PyMorph and SDSS are real. In particular, our work shows that previous halo model analyses which used SDSS pipeline photometry when relating the stellar mass of the central galaxy to the dark matter mass of the halo which surrounds it should be redone. Ascribing the difference compared to analyses based on Sersic magnitudes to the semantics of whether or not one includes the ICL when describing the stellar mass of the central galaxy is incorrect.

We also showed that, when studying correlations at fixed group mass, care must be taken to ensure that one does not mistake selection effects for physical effects. The SerExp and Model magnitude difference, when plotted versus SerExp magnitude, shows strong trends even though no such trends were apparent when plotting versus Model magnitude (compare red and green curves in Figures 7 and 15). We argued that this is because the group masses in the Yang+ catalog are tightly correlated with Model magnitudes. A simple model illustrated why selection effects appear when one studies correlations at fixed group mass with a quantity which was not used to define the groups (Appendix B). Indeed, when we reassign halo masses to the Yang+ groups on the basis of SerExp magnitudes, then the selection effect appears when plotting versus Model rather than SerExp magnitudes (compare Figures 8 and 17). This conclusion is general. For example, we have checked (but do not show here) that the selection effect is even stronger when Ser rather than SerExp magnitudes are used. We hope our demonstration of the existence of this sort of pernicious selection effect, our explanation of its cause (Appendix B), and how one should correct for it will prevent future confusion.

Acknowledgements

We would like to thank the referee for a helpful report. MB and RKS are grateful to the ICTP for its hospitality during the summer of 2016.

Appendix A Comparison with single Ser photometry

The main text studied differences between PyMorph SerExp and SDSS pipeline photometry. Here we show the corresponding results if we use (truncated) PyMorph Ser instead of (truncated) SerExp values. Meert et al. (2013) and Bernardi et al. (2014) have argued that the SerExp photometry we present in the main text is more unbiased than the Ser photometry here. However, single Ser fits are more common in the literature, which is why we show them here.

Figure 13 shows that the differences with respect to SDSS sky are similar to when fitting SerExp profiles (compare Figure 11), but the remaining difference is substantially larger than before. In addition, even though the differences are larger, centrals and satellites are similar when shown as a function of PyMorph photometry but different when shown as a function of Model magnitudes (left and right panels, respectively), as was the case in Figure 11. Figure 14 shows that this is also true if we use Ser photometry to determine the stellar masses.

Appendix B Selection effects for analyses at fixed halo mass

We noted in the main text that it was important to present results using the same photometry that was used to define the centrals, satellites, and halo masses in the Yang+ group catalog. Figure 7 shows the result of plotting the PyMorph-SDSS difference versus Model magnitude. Figure 15 shows that the scalings change dramatically when we plot versus SerExp instead. The faintest satellites (green) and brightest centrals (red) are remarkably similar to the average (grey), but the brightest satellites and faintest centrals curve away. We see similar effects if we replace the two-component SerExp with single-component Ser photometry, and if we use $M_{*}=(M_{*}/L)\,L_{\rm Ser}$ or $(M_{*}/L)\,L_{\rm SerExp}$ instead of $(M_{*}/L)\,L_{\rm Model}$. We now show that these offsets are selection effects which result from the fact that the group catalog was constructed using Model magnitudes. For brevity, we only consider the trends for centrals.

Let $m_{i}$ denote the model magnitude of galaxy $i$, $s_{i}$ its SerExp magnitude, and $h_{i}$ the mass of the halo to which it belongs. Let $n(m,H)$ denote the number density of objects which have Model magnitude $m$ and reside in halos with mass $H\equiv(h_{\rm min}\leq h\leq h_{\rm max})$. Then,

[TABLE]

where the $H$ indicates the limited range in halo masses. The measurement in Figure 7 is

[TABLE]

If the scatter between $s$ and $m$ does not depend on halo mass, then we can set $p(s|m,h)=p(s|m)$, and so

[TABLE]

This is the same for all choices of $H$, and is consistent with the fact that all the curves in Figure 7 are the same.

In contrast, the measurement in Figure 15 is

[TABLE]

where

[TABLE]

To gain insight, suppose that $p(m|s,h)$ is Gaussian. Then the mean $\langle m|s,h\rangle$ can always be written as $\langle m|s\rangle+C\,[h-\langle h|s\rangle]$, where $C$ does not depend on either $s$ or $h$. Then

[TABLE]

The second term is zero if $C=0$; this only happens if the $mh$ correlation is entirely due to the $ms$ and $sm$ correlations. If $C\neq 0$, then the second term equals zero only if $H$ allows the full range of halo masses. However, if $H$ selects only a subset of halos, then the integral is non-vanishing, yielding a correction factor which depends on $H$. If $p(s|h)$ were also Gaussian, then the integral above would equal $\sigma_{h|s}\,[\exp(-y_{\rm min}^{2}/2)-\exp(-y_{\rm max}^{2}/2)]/\sqrt{2\pi}$ and the denominator would be $[{\rm erf}(y_{\rm max}/\sqrt{2})-{\rm erf}(y_{\rm min}/\sqrt{2})]/2$, where $y_{\rm max}=(h_{\rm max}-\langle h|s\rangle)/\sigma_{h|s}$ and similarly for $y_{\rm min}$.

To test this, we assumed that the distribution of $m$ is Gaussian with unit variance, that $p(s|m)$ is Gaussian with mean $1.025m$ and rms 0.2, and $p(h|m)$ is Gaussian with mean $-0.4m$ and rms 0.1. While these scalings are not exactly the same as those in the data, they are similar. Note in particular that $h$ is directly related to $m$; it only correlates with $s$ because $s$ too is directly related to $m$. (I.e., the factor $C$ is proportional to the $hm$ correlation.)

We first made Monte-Carlo realizations of this joint distribution of $m$, $s$, and $h$, and checked that measurements of \bigl{\langle}s-m|m,H\bigr{\rangle} do not depend on $H$, in good agreement with equation (4). Figure 16 shows a similar analysis of \bigl{\langle}s-m|s,H\bigr{\rangle}; i.e., when we plot the difference between $s$ and $m$ versus $s$ rather than $m$. The straight black line shows $s-\langle m|s\rangle$ as a function of $s$, and the other three curves show equation (7) with $(h_{\rm min},h_{\rm max})=(-5,-0.5)$, $(-0.25,0.25)$, and $(0.75,1.75)$. We have chosen the different bin centers to illustrate the mass dependence of the offsets, and the different bin widths to show that these offsets appear primarily at the edges of the bins. The error bars show the corresponding measurements of \bigl{\langle}s-m|s,H\bigr{\rangle} in the Monte-Carlo realizations in which \bigl{\langle}s-m|m,H\bigr{\rangle} showed no offsets. The agreement between the measurements and equation (7) indicates that we understand the origin of the offsets from the straight line.

In particular, the fact that we see no offsets when we plot versus $m$ indicates that the offsets here, and in Figure 15, are selection effects – they are not physical. Additional evidence that they are selection effects comes from comparison of Figures 8 and 17. Both use the Yang+SerExp catalog defined in Section 3.2 for which halo masses were assigned to the Yang+ groups on the basis of SerExp rather than Model magnitudes. In this case, plots as a function of SerExp magnitudes were well-behaved (Figure 8), whereas those for Model show offsets (Figure 17). I.e., the selection effect is reversed, as expected based on the analysis above.