Comparing PyMorph and SDSS photometry. I. Background sky and model fitting effects

TL;DR

This paper compares galaxy photometry from PyMorph and SDSS, highlighting how background sky estimation and model choice significantly affect luminosity measurements, especially for luminous galaxies.

Contribution

It demonstrates that PyMorph provides more accurate sky estimates and luminosities than SDSS, emphasizing the importance of background and model fitting choices in galaxy photometry.

Findings

PyMorph sky estimates are fainter than SDSS, aligning with Blanton et al. (2011).

SDSS sky biases luminosities by over 0.1 mag for large galaxies.

PyMorph fits yield more light than deVaucouleurs but less than Sersic models for E+S0 galaxies.

Abstract

A number of recent estimates of the total luminosities of galaxies in the SDSS are significantly larger than those reported by the SDSS pipeline. This is because of a combination of three effects: one is simply a matter of defining the scale out to which one integrates the fit when defining the total luminosity, and amounts on average to < 0.1 mags even for the most luminous galaxies. The other two are less trivial and tend to be larger; they are due to differences in how the background sky is estimated and what model is fit to the surface brightness profile. We 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). Using the SDSS sky biases luminosities by more than a few tenths of a magnitude for objects with half-light radii > 7 arcseconds. In the SDSS main galaxy sample these are typically luminous galaxies, so they are not necessarily nearby. This bias becomes worse when allowing the model more freedom to fit the surface brightness profile. When PyMorph sky values are used, then two component Sersic-Exponential fits to E+S0s return more light than single component deVaucouleurs fits (up to ~0.2 mag), but less light than single Sersic fits (0.1 mag). Finally, we show that PyMorph fits of Meert et al. (2015) to DR7 data remain valid for DR9 images. Our findings show that, especially at large luminosities, these PyMorph estimates should be preferred to the SDSS pipeline values.

Equations

$$L_{\mathrm{trunc}}=L_{\infty }\frac{\gamma (2n,b_n{\rho }_{\mathrm{trunc}}^{1/n})}{\Gamma (2n)},\mathrm{where}{\rho }_{\mathrm{trunc}}\equiv \frac{{\theta }_{\mathrm{trunc}}}{\sqrt{ab}},$$

Full text

1 Introduction

There is substantial interest in quantifying the luminosity and stellar mass functions in the local universe (Bernardi et al. 2017a and references therein). The Sloan Digital Sky Survey (hereafter SDSS), which surveyed about a quarter of the sky to a median redshift of about $z\sim 0.1$, is the benchmark database for such studies. Recently Meert et al. (2015, 2016) have made available a re-analysis of the galaxies in the SDSS DR7 release (Abazajian et al. 2009). Their analysis determines photometric parameters, such as luminosity, half-light radius, a measure of the steepness or central concentration of the profile, etc., by fitting a number of different models to the surface brightness profile: a single component deVaucouleurs profile, a single component Sersic profile, and a two component Sersic bulge plus exponential disk profile (hereafter deV, Ser and SerExp). The fitting algorithm is called PyMorph (Vikram et al. 2010; Meert et al. 2013, 2015, 2016; Bernardi et al. 2014). The PyMorph catalog yields substantially more light at high luminosities (Bernardi et al. 2013, 2016a,b, and Figure 1 below) than previous work based on SDSS pipeline photometry. The differences impact Halo Model (Cooray & Sheth 2002) based interpretations of the relationship between galaxies and dark matter halos at $z\sim 0.1$ (e.g. Shankar et al. 2014). Pinning down this relationship locally is crucial for studies of how this relationship evolves.

In addition, as first identified by Bernardi et al. (2011) at the high mass (luminosity) end there is a special mass (luminosity) scale: $2\times 10^{11}M_{\odot}$ (which corresponds to an $r-$band luminosity scale of $\sim-22.5$ mag). Various scaling relations change slope at this scale, and this is thought to be related to a change in the assembly histories – e.g. minor versus major dry mergers. It is also the mass (or luminosity) scale where the stellar mass (or luminosity) function starts to drop exponentially. For all these reasons, identifying and accounting for all possible biases so as to have reliable photometric estimates at these luminosity and mass scales is important. Here we address the reasons for the differences between PyMorph and the SDSS, and show that PyMorph should be used, especially at large luminosities.

There are expected to be three main culprits. An important step in the determination of the amount of light we receive from an object is the estimation of the amount of light which is contributed by the background sky. Over-estimating the contribution from the sky will lead to an underestimate of the size and total light, and perhaps a decrease in the estimate of how centrally concentrated the object is. Bernardi et al. (2007) (see also, e.g., SDSS DR7 documentation) noted that the SDSS pipeline reductions underestimated the sky, especially in crowded fields. In the years since, the SDSS has revised its pipelines (see the DR9, Ahn et al. 2012, and subsequent data releases).

In addition, a number of other analyses have also provided improved estimates (Simard et al. 2011, Blanton et al. 2011, Meert et al. 2015, 2016). One of the main goals of the present work is to compare different estimates of the sky in the SDSS footprint, and to quantify the impact this has on the estimated sizes, shapes and luminosities of galaxies. Blanton et al. (2011) argue that the SDSS values can be biased by as much as a magnitude for nearby objects with large angular size (half-light radius $\geq 40$ arcseconds). However, because the bias is really associated with having a large angular size, the bias can still be significant (a few tenths of a magnitude) for large objects (half-light radius $\geq 7$ arcseconds) whether or not they are nearby. There is a tight correlation between luminosity and physical size, so even though the majority of luminous galaxies in the SDSS main galaxy sample tend to be more distant ($z\sim 0.2$) they still have relatively large angular sizes ($\geq 7$ arcseconds).

In addition to the sky, two other effects contribute to differences between SDSS pipeline and more recent estimates of galaxy luminosities and sizes. One is trivial: when reporting the total light in an image, the SDSS only integrates the surface brightness profile out to about $\sim 7\times$ the half-light radius. Others, such as PyMorph (Meert et al. 2013), do not truncate. This amounts to a small systematic difference of order 0.05 mags for deV profiles, but can be larger for Ser profiles (e.g. Kelvin et al. 2012). The second effect is more interesting: it is the fact that the luminosity and size estimates depend on the model which is fitted to the image. In what follows, we will be careful to distinguish between these three effects. E.g., it is not obvious if models which have more freedom to better fit the image will end up predicting more light or less.

There is another potential observational systematic: the deblending of overlapping galaxies. However, this is resolved in Meert et al. (2015), who discuss how PyMorph handles nearby neighbours, as well as polluted fits (those that could not be deblended). Their rate of occurence is sub-percent, and PyMorph provides a flag identifying them, so it was simple to exclude them from the analysis which follows.

The present study is timely because the Meert et al. analyses are based on SDSS DR7 images. However, significant changes to the SDSS imaging pipeline were implemented in DR9, and remain in place in subsequent data releases. These are described on the SDSS website: www.sdss.org. Therefore, after defining the sample we work with in Section 2.1, our first step is to compare PyMorph analyses of the DR7 and DR9 images. This is the subject of Section 2.2. Section 2.3 and Section 2.4 quantify the effects of truncation. Section 2.4 also highlights the fact that, because the most massive objects may be a different population having different profile shapes it is important to specify the choice of regression, i.e. the average magnitude difference may depend on the luminosity being used as x-axis. Section 3 compares sky estimates from the SDSS DR7 and DR9 pipelines with determinations from Blanton et al. (2011), Simard et al. (2011), and Meert et al. (2015, 2016) (hereafter B11, Simard11 and PyMorph DR7, respectively). Section 4 shows how the choice of model to fit affects the estimated total light. A final section summarizes.

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 SDSS and PyMorph

The analysis which follows is based on the SDSS DR7 and DR9 Main Galaxy samples. For these galaxies, the SDSS provides a number of photometric parameters on its website: www.sdss.org. We are most interested in the total magnitudes and half-light radii, the best SDSS pipeline estimates of which are based on fitting exponential or deVaucouleurs profiles to the sky subtracted image. Model magnitudes simply choose the better of the two fits, whereas cModel magnitudes use a linear combination of the two best fits (a $\chi^{2}$-like goodness of fit metric is minimized to set the relative amplitudes of the components). Thus, although they are the result of fitting two profile shapes, cModel magnitudes are not really two-component fits.

In contrast to the SDSS cModel photometry, the best PyMorph SerExp photometry is based on true two-component fits – a Sersic bulge with an exponential second component – in which the sky, assumed to be constant across the image, is also fit simultaneously (e.g. Meert et al. 2015). These fits were made using the DR7 release.

2.1 Motivation

We begin with a comparison of what are considered to be the best SDSS and PyMorph photometry: cModel and SerExp magnitudes. Figure 1 shows that the two are in good agreement, except at the bright end, where PyMorph is substantially brighter. The bottom panel shows the result of replacing cModel with Model magnitudes. Except for an offset at low and intermediate luminosities, both panels show similar trends. The similarity observed at the bright end is expected because the vast majority of the most luminous galaxies are E+S0s, so Model = deV and cModel $\approx$ deV.

The main goal of the present study is to determine which of the three culprits mentioned in the Introduction are responsible for the offsets in Figure 1. In particular, it may be that the agreement between cModel and SerExp at faint and intermediate luminosities is fortuitous.

Figure 1 was made using DR7 galaxies. However, between DR7 and DR9, a number of parts of the SDSS pipeline were changed. The most important change is the SDSS sky estimate, but how flux calibration is done, and so on, also changed (see Aihara et al. 2011 and Ahn et al. 2012 for details). Therefore, our first step is to determine if the changes from DR7 to DR9 matter. Figure 2 shows a similar comparison as in Figure 1, but now using DR9 values. To make this figure we ran PyMorph on a subset of $10^{4}$ DR9 galaxies. The chosen objects are the same as those used by Meert et al. (2013) when developing and testing PyMorph. The distribution of the measured parameters of this subset reproduces the distribution of all the observed galaxies in the SDSS DR7 main galaxy sample (see their Figure 1). Comparison of Figure 2 with Figure 1 shows little difference: the discrepancy between SDSS and PyMorph which was known to exist in DR7 persists in DR9.

2.2 Comparison of SDSS DR7 and DR9

We now consider if the best fitting PyMorph parameters have changed between DR7 and DR9. Since PyMorph fits for the sky itself – it does not use the SDSS value – we expect the change to the SDSS sky estimate to have little impact on the PyMorph fits. Figure 3 shows that this is indeed the case: the apparent magnitudes, sizes, and Sersic indices for PyMorph Ser fits are essentially unchanged. Figure 4 shows that this is also true for PyMorph SerExp fits; because these are two-component fits, the bottom panel shows bulge/total ratios rather than Sersic indices. Both figures show that, although there is scatter between the DR7 and DR9 values, it is similar to the statistical uncertainty on the parameters (Meert et al. 2013). It should be noted that there is larger scatter for SerExp than for Ser, because there are more free parameters and hence more potential degeneracies. Therefore, Figures 3 and 4 indicate that the PyMorph parameters of Meert et al. can be used essentially without modification even for DR9. (There are, of course, other studies for which the difference between DR7 and DR9 or DR13 recalibrations do matter.)

2.3 Effect of truncation

In what follows, we would like to compare the luminosity estimates of PyMorph and the SDSS. Both report values based on fitted models; however, whereas PyMorph integrates the fitted profile to infinity, the SDSS does not. If a two-dimensional Sersic profile with semi-major axis $a$ and axis ratio $b/a$ is truncated along a line of constant surface brightness, then

[TABLE]

$\gamma(m,x)$ is the incomplete gamma function, $\gamma(m,\infty)=\Gamma(m)$, and $b_{n}$ is defined by requiring $\gamma(2n,b_{n})=\Gamma(2n)/2$. E.g., $b_{n}\approx 7.669$ when $n=4$.

The ratio $L_{\rm trunc}/L_{\infty}$ clearly depends on $n$. Notice that if $\theta_{\rm trunc}$ is a multiple of $\sqrt{ab}$, then, at fixed $n$, the correction is the same for all axis-ratios. For example, in their work with the GAMA survey, Kelvin et al. (2012) set $\theta_{\rm trunc}=10\sqrt{ab}$. For reasons which will become clear shortly, the grey solid curve in Figure 5 shows $\theta_{\rm trunc}=7.5\sqrt{ab}$. This shows that, when $n=4$, the correction is 0.07 mags. Unfortunately, the SDSS truncation is more complicated: the SDSS website says that it truncates with a function which drops from unity to zero between $7$ and $8\times$ the half-light radius. However, in the database, the quantity which is called $r_{e}$ is the semi-major axis $a$, rather than $\sqrt{ab}$. In addition, the actual form of this truncation has never been published. As we show below, we are able to reproduce the SDSS values if we use a sharp truncation radius of $7.5a$ making $\rho_{\rm trunc}^{\rm SDSS}\approx 7.5a/\sqrt{ab}=7.5\,\sqrt{a/b}$. (In particular, $7.5a$ works substantially better than $7.5\sqrt{ab}$.) Hence, at fixed $n$, $L_{\rm trunc}/L_{\infty}$ is a monotonic function of $b/a$: since $0\leq b/a\leq 1$, the correction is maximal when $b/a=1$ and $L_{\rm trunc}\to L_{\infty}$ as $b/a\to 0$. Thus, at fixed $n$, there is a range of corrections which depends on the distribution of $b/a$. Since our goal is to compare with the SDSS, the black solid line in Figure 5 shows the median of $2.5\log_{10}(L_{\rm trunc}/L_{\infty})$ as a function $n$, and the scatter around this median (black dashed lines), for the PyMorph Sersic reductions of SDSS E+S0 galaxies when $\theta_{\rm trunc}=7.5a$. This shows that, when $n=4$, the correction is $\sim 0.05$ mags, but when $n=8$, then the median correction is $\sim 0.16$ mags. (For later type galaxies $n$ is smaller so the correction is smaller; the blue dot shows the correction if $n=1$ and one truncates at $3.5\times$ the half light radius.)

In what follows, we will be careful to indicate if the reported magnitudes were based on truncation or not. However, the half-light radii we report are always those which include $L_{\infty}/2$; we never use the scale associated with $L_{\rm trunc}/2$.

2.4 Choice of regression and truncation

We remarked in the Introduction that the most massive galaxies appear to be a structurally different population. So it should not be surprising if their surface brightness profiles are also different in some way. If these are objects for which SDSS and PyMorph photometry is particularly different, then plots versus PyMorph may look rather different from plots versus Model, for the same reason that, in a Gaussian mixture model, plots of $y$ versus $x$ can look very different from plots of $x$ vs $y$. Figure 6 shows that something like this happens in the SDSS data: The differences between Ser and Model magnitudes increase at the bright end, but they look much larger when shown as a function of Ser rather than Model magnitudes. Figure 7 shows that the same is true of SerExp magnitudes.

The differences are reduced slightly if one uses truncated Ser or SerExp magnitudes, since this reduces the analog of $m$ in the example above, but it does not change the fact that the choice of $x$-axis matters. While truncation matters, the net effect of truncation is about half of what one would naively have expected from Figure 5. This is because the correction depends on $n$, but the $n$-$L$ correlation is weak. Although large $L$ have larger $n$, so truncation matters more for large $L$, there is substantial scatter around the mean $n$ which reduces the net effect. This is also why, in practice, it matters little ($\leq 0.01$ mags) whether one truncates using $7.5a$ or $7.5\sqrt{ab}$. Of course, truncation matters even less for the SerExp fits.

3 Comparison of sky estimates

Our goal is to compare PyMorph and SDSS sky estimates. However, when fitting a model to the observed galaxy image, PyMorph fits for the sky – assumed to have constant surface brightness across the image – simultaneously. Therefore, it is possible that the fitted sky varies when the model which is fitted to the galaxy surface brightness profile changes. This would make comparisons with the SDSS sky estimate depend on the fitted model. Fortunately, Figure 8 shows that the estimated sky is essentially the same whatever the fitted model. (We have plotted versus truncated magnitudes. Of course, the y-axis legend does not specify truncated because the sky estimates do not depend on (i.e. are the same) whether or not we truncate.)

This has two consequences. First, when comparing the PyMorph sky with other estimates, we do not need to specify if it is the deV sky, the Ser sky, or the SerExp sky, since, for the present purposes, they are all the same. We exploit this fact in Section 3.1. Second, the similarity in sky values indicates that differences between PyMorph models are not driven by the sky. We use this fact in Section 4.

We are now ready to compare background sky estimates with those from PyMorph, for which we use the SerExp sky value. Figure 9 compares background sky estimates from SDSS DR7, Simard11, and PyMorph DR7 SerExp. The PyMorph sky is faintest and SDSS brightest, with the Simard11 sky lying closer to the SDSS at the faint and intermediate luminosities and in between at the bright end. The differences from PyMorph are particularly large for objects with large angular sizes or luminosities.

While it is tempting to conclude that Simard11 is the most prudent choice because it lies between the other two, Figure 10 shows that the PyMorph sky estimate is in excellent agreement with that of B11. In contrast to the previous figure, this one uses DR9 images, for which Simard11 values are not available.

The PyMorph and B11 sky values were determined in very different ways. Those of B11 are based on fitting the masked background sky for each SDSS scan with a smooth continuous function across the sky. (In Figure 10, we used the B11 sky value measured at the center of the galaxy image since the variation of the sky value on the scale of a galaxy is very small.) In contrast, the PyMorph sky is determined on an object-by-object basis. Therefore, the agreement between the two is nontrivial, and strongly suggests that these two estimates are to be preferred over the others. Note also that the scatter around the median is symmetric, whereas in the comparison with SDSS it is not.

B11 argue that their sky estimates represent a substantial improvement over the standard SDSS catalog results and should form the basis of any analysis of nearby galaxies using the SDSS imaging data. Figure 10 shows that, in fact, this is not restricted to nearby galaxies: E.g., for all galaxies with apparent sizes $\geq 7$ arcseconds, the SDSS sky is biased (left panel) (we quantify its effect on photometric parameters in Section 3.1). In the SDSS main galaxy sample these tend to be galaxies with large luminosities (right panel) which are typically in crowded fields. The agreement between B11 and PyMorph in both panels suggests that, in contrast to the SDSS, PyMorph is unbiased for large luminous galaxies.

3.1 Sky-related biases when fitting deVaucouleurs profiles to E+S0s

Having determined that the PyMorph/B11 sky is to be preferred, we now consider how the choice of sky biases the inferred parameters. We begin with a study of the only case in which a direct comparison (i.e. same model fit) with SDSS is possible: fitting a deVaucouleurs profile to images of E+S0 galaxies. For morphological type, we use the Bayesian automated classifications of Huertas-Company et al. (2011): each galaxy is assigned weights which represent the probabilities that it is Elliptical, S0, Sab or Scd. We restrict to E+S0s since a deVaucouleurs profile is known to not fit other morphological types well, and we do not wish to confound the question of sky-related biases with biases arising from fitting a bad model.

The top panel of Figure 11 shows the difference between the SDSS DR9 and PyMorph estimates of the total (truncated) magnitude. The black curve shows results for the DR9 E+S0 subset while the gray curve shows the results for the larger ($\sim 60\times$) PyMorph DR7 E+S0 sample. These curves show that SDSS is fainter, and this difference increases for the largest (left) and most luminous (right) galaxies. This is a consequence of three effects: (i) the SDSS sky is brighter, so galaxies with large angular radii tend to have their sizes reduced by a bigger factor, as a result of which less light is assigned to the galaxy; (ii) the total magnitude is computed by integrating the surface brightness profile, and our model of how the SDSS truncates this integral (equation 1 and related discussion) may not be accurate; (iii) the SDSS and PyMorph fitting routines are systematically different.

To remove the latter two effects, the yellow curve shows the difference between forcing PyMorph to use the SDSS sky values when fitting and the original PyMorph value. Since both estimates are from PyMorph DR9, effects (ii) and (iii) have been removed, so the yellow curves differ from zero entirely because of the differences in sky values (the SDSS sky is brighter). Moreover, the fact that these yellow curves agree with the previous black ones to better than 0.01 mags strongly suggests that we have modelled the SDSS truncation algorithm correctly: PyMorphSDSSsky,deV,trunc is a good proxy for SDSSdeV. (As an aside, this means that the good agreement at magnitudes fainter than $\sim-22$ mag in the top panels of Figures 1 and 2 is fortuitous, at least where the contribution from E+S0s is significant.)

Figure 10 shows that while the SDSS sky is brighter than PyMorph, the B11 sky is in excellent agreement across the entire population. The bottom panels of Figure 11 show that if PyMorph is forced to use the B11 sky estimate rather than its own (in practice, this means PyMorph is made to fit the B11 sky-subtracted image provided on the SDSS website, while forcing its own additional sky estimate to be zero across the image), then the median difference in magnitude is negligible. Notice that the scatter around the median is less than 0.03 mags; this level of agreement is remarkable. Comparison of the top and bottom panels shows that the sky can introduce biases of order 0.1 mags or more for the most luminous objects when fitting deVaucouleurs profiles.

3.2 Sky-related biases in Ser and SerExp fits

We now consider sky-related biases when fitting other models. Figure 12 shows results for PyMorph SerExp fits to all galaxies as the restriction to E+S0s is no longer necessary. The yellow curves in the different panels show that the brighter SDSS sky biases the estimated SerExp magnitude fainter, and this bias is most severe for the largest (top left) and/or most luminous (top right) galaxies; it also biases the half-light radii and B/T values to smaller values (bottom left and right, respectively). For PyMorph SerExp fits the biases arising from the SDSS sky are significantly larger than when fitting deVaucouleurs profiles. However, there are no such biases associated with the B11 sky values (red curves). On the other hand, although the scatter is similar to that in the bottom left panel of Figure 11, degeneracies between the fitted SerExp parameters and the fitted sky contribute to increased scatter at high luminosities.

Figure 13 shows a similar analysis of Sersic rather than SerExp fits. In this case, there is only one component, so the bottom right panel shows the Sersic index $n$ rather than B/T. Again, the SDSS sky biases the estimated Ser magnitude fainter, and this bias is most severe for the most luminous (top left) and/or largest (top right) galaxies; note that now the bias can be as large as 0.4 mags – substantially larger than when fitting deVaucouleurs profiles. The SDSS sky also biases the half-light radii and Sersic indices to smaller values (bottom left and right, respectively). While there are no such biases associated with the B11 sky values, there are hints of a small bias at the largest angular sizes and luminosities. Since there are fewer free parameters compared to SerExp, and therefore fewer degeneracies, we would expect the scatter around the zero-median to be smaller. This is indeed the case for intermediate and low luminosity galaxies which usually have a Sersic index $n<4$. A Sersic fit with a higher $n$ is more sensitive to differences in the background sky (Meert et al. 2013). Thus, the larger scatter observed at large sizes and/or luminosities is due to the fact that the most luminous galaxies usually have $n\geq 4$.

The results of this subsection have an interesting connection to recent work. D’Souza et al. (2015) state that image stacking is essential for recovering unbiased estimates of the total light. Their stacks were of DR9 sky subtracted images, meaning that they assumed the B11 sky estimate was correct. The results in the bottom halves of each panel in Figures 11–13 were based on analyses of individual images. Since no stacking was performed when fitting, the lack of bias between the full PyMorph values and those when the sky is fixed to that of B11 shows that stacking is not a prerequisite for obtaining unbiased results.

In this context, it is interesting to compare the difference between PyMorph SerExp and SDSS Model magnitudes. Bernardi et al. (2017a) have already shown that the median difference is the same as what D’Souza et al. find from their stacking analyses (see their Figure 2). But they left open the question of the scatter around the median. Since their work used SerExp magnitudes in which PyMorph also fit for the sky, it is possible that some of the scatter is reduced when PyMorph is forced to use the B11 sky. The top panel of Figure 14 shows that this is not the case: the differences between PyMorph SerExp and SDSS Model magnitudes when PyMorph fits its own sky and when the sky is fixed to that of B11 are very similar, not just in the median but also the scatter around it. (Because we are using truncated magnitudes, comparing PyMorph to SDSS quantities, the offset from zero is due to differences in sky and fitted-model only.) This strongly suggests that the scatter reflects true differences between SerExp and Model (i.e. deV) models; it is not dominated by degeneracies arising from fitting the sky simultaneously. To remove trends which arise from morphology, the bottom panel shows a similar analysis for the subset of galaxies classified as E+S0s. While the trends differ especially at low luminosities – where non-E+S0s begin to dominate in the top panel – it is still true that changing from PyMorph to B11 sky values makes little difference.

4 Dependence on fitted model

Having shown the large biases associated with the SDSS, we now turn exclusively to PyMorph values. Recall that the PyMorph sky values are essentially the same for all fitted models (Figure 8), so that comparison of different PyMorph fits show how the luminosity and size depend on the functional form assumed for the surface brightness profile. Also, when comparing results from deVaucouleurs profiles we show results for E+S0s only to avoid the issue of biases which arise from using a functional form which is known to provide a poor fit.

Figure 15 shows that SerExp fits to E+S0s return more light than deV fits especially at large luminosities (up to $\sim 0.2$ mag); when shown as a function of SerExp luminosity, the difference is largest for the most luminous galaxies. (This analysis was done using the DR7 E+S0s, for which PyMorph reductions are available, since it is much larger ($\sim 60\times$) than the subset of DR9 galaxies on which PyMorph was rerun. The result from the DR9 subset is noisier, but otherwise very similar, so we have not included a separate figure showing it.)

The difference due to fitting different models is similar in amplitude to that in Figure 11, which was due to differences in the estimated sky. However, the dependence on choice of regression is more dramatic here than in Figures 6 and 7 because there the effects of the sky somewhat compensated for the difference in profiles. By using only PyMorph quantities here, the sky effects have been removed.

Figure 16 shows a similar comparison, but now between SerExp and Ser fits to DR7 E+S0s. Clearly, SerExp is about 0.1 mags fainter and 10% smaller across the E+S0 population. Finally, Figure 17 compares SerExp and Ser fits to the full DR7 population. At high luminosities, this figure is very similar to the previous one, because most high luminosity galaxies are E+S0s. However, there are small differences at low luminosities. These indicate that Sersic luminosities and sizes of non-E+S0s must be fainter and smaller than the corresponding SerExp values.

The differences between the cyan and magenta curves in Figures 15 and 16 at the high luminosity end strongly suggest that the most luminous galaxies have different surface brightness profiles from the bulk of the population. This can be understood as follows. Suppose we have two populations, both of which span the same range of deV. Assume that, for one, Ser$=$deV, but that Ser$=$deV$-m$ for the other (i.e. Ser is $m$ mags brighter). Let $f$ denote the fraction of objects in this second population. Then, a plot of $\Delta M\equiv{\tt deV}-{\tt Ser}$ when shown as a function of deV will look like two horizontal lines, one lying $m$ mags above the other, like this: $=$. The average of $\Delta M$ when shown as a function of deV will equal $fm$. However, when shown as a function of Ser, the second population will be displaced brightwards along the x-axis:   $-$. As a result, where the two populations overlap, the mean $\Delta M$ will still be $fm$, but at the brightest Ser the mean will be $m$. If the second population only spans a limited range of deV, $-=-$, then the average as a function of deV will show curvature, and may not even be monotonic, whereas the other may still be monotonic: $--=$. Alternatively, suppose that when plotted as a function of Ser, $\Delta M$ is made of two populations:    $|$. Then, when plotted as a function deV, this will look like   . Again, the plot versus Ser will be monotonic, whereas that versus deV will not. In practice, this mix of populations means that if one wishes to use the mean of $\Delta M$ as a measure of the difference between deV and Ser, then one must specify which variable was being held fixed (we made a similar point in the context of Figures 6 and 7).

5 Conclusions

In both SDSS DR7 and DR9, PyMorph returns brighter estimates of the total light of a galaxy than either SDSS Model or cModel magnitudes (Figures 1 and 2). While the SDSS values have changed slightly between DR7 and DR9, the PyMorph fits to the DR7 release provided by Meert et al. (2015, 2016) remain accurate for DR9 as well (Figures 3 and 4). Some of the difference with respect to the SDSS arise from the fact that the SDSS value for the total brightness comes from truncating the integral over the surface brightness profile (Figures 5, 6, and 7). We believe we understand the truncation algorithm (Figure 11 and related discussion), and so in all our subsequent comparisons with the SDSS, we have truncated the PyMorph values using a similar algorithm (equation 1) so that truncation plays no further role in the PyMorph-SDSS differences.

The sky estimated by PyMorph is almost completely independent of the model used to fit the galaxy (Figure 8). The PyMorph sky estimates are fainter than those of the SDSS DR7 or DR9 pipelines (Figure 9), but are in excellent agreement with the estimates of B11 (Figure 10). The difference in sky accounts for about half of the discrepancy shown in Figures 1 and 2. In addition, there is an overall offset of about 0.07 mags which comes from the fact that the SDSS value for the total brightness comes from truncating the integral over the surface brightness profile (Figure 7). The remainder arises from fitting different models.

Use of the SDSS sky biases luminosities and half-light radii to lower values; in the main SDSS galaxy sample these biases are significant (a few tenths of a magnitude) at large luminosities: they matter not just for nearby galaxies. The biases become even worse when allowing the model more freedom to fit the surface brightness profile (Figures 11–13).

When PyMorph sky values are used, the SerExp fits to E+S0s return more light than deV fits especially at large luminosities (up to $\sim 0.2$ mag), but less light than Ser fits (Figure 16). For non-E+S0s, which are dominant towards lower luminosities, Sersic luminosities and sizes are slightly fainter and smaller than SerExp (Figure 17).

Our findings show that, especially at large luminosities, SDSS pipeline values should not be used: PyMorph estimates are much more reliable. Of these, Meert et al. (2013) and Bernardi et al. (2014) have already shown that the SerExp values are to be preferred. The PyMorph SerExp values are also consistent with results obtained via the stacking analysis of D’Souza et al. (2015) (Figure14; see also Figure 2 in Bernardi et al. 2017a). This is reassuring because the two analyses are very different. However, this does raise the question of why SerExp is better than SDSS pipeline photometry. E.g., since the largest discrepancies occur at high luminosities, and the most luminous galaxies are preferentially found in clusters, is it possible that the SerExp fits are different because the second component is actually fitting intercluster light? Bernardi et al. (2017b) show that for the vast majority of massive galaxies this is almost certainly not the main reason for the difference.

The assembly history of a galaxy is expected to leave an imprint on its surface brightness profile. Indeed, we find significant evidence that the surface brightness profiles of the most luminous galaxies suggest that they are a distinct population (Figures 6, 7, 15 and 16). Therefore, we hope our results will inform studies of the assembly histories of the most massive galaxies.

Acknowledgements

We thank R.K. Sheth and V. Vikram for helpful discussions, and the referee for a helpful and competent report.