<figure><img src="image_1.png" /></figure>

# The high-mass end of the stellar mass function: dependence on stellar population models and agreement between fits to the light profile

M. Bernardi,<fnref n="1" />★ A. Meert,<fnref n="1" /> R. K. Sheth,<fnref n="1" />,<fnref n="2" /> J.-L. Fischer,<fnref n="1" /> M. Huertas-Company,<fnref n="3" /> C. Maraston,<fnref n="4" /> F. Shankar<fnref n="5" /> and V. Vikram<fnref n="1" />

Accepted 2017 January 19. Received 2017 January 17; in original form 2016 April 4

## ABSTRACT

We quantify the systematic effects on the stellar mass function that arise from assumptions about the stellar population, as well as how one fits the light profiles of the most luminous galaxies at $z \sim 0.1$. When comparing results from the literature, we are careful to separate out these effects. Our analysis shows that while systematics in the estimated comoving number density that arise from different treatments of the stellar population remain of the order of $\le 0.5$ dex, systematics in photometry are now about 0.1 dex, in contrast to some recent claims in the literature. Compared to these more recent analyses, previous work based on Sloan Digital Sky Survey pipeline photometry leads to underestimates of $\rho_*(\ge M_*)$ by factors of 3–10 in the mass range $10^{11}-10^{11.6} M_{\odot}$, but up to a factor of 100 at higher stellar masses. This impacts studies that match massive galaxies to dark matter haloes. Although systematics that arise from different treatments of the stellar population remain of the order of $\le 0.5$ dex, our finding that systematics in photometry now amount to only about 0.1 dex in the stellar mass density is a significant improvement with respect to a decade ago. Our results highlight the importance of using the same stellar population and photometric models whenever low- and high-redshift samples are compared.

**Key words:** galaxies: fundamental parameters – galaxies: luminosity function, mass function – galaxies: photometry.

## 1 INTRODUCTION

Our knowledge of galaxy formation and evolution has grown tremendously over the last two decades. This knowledge comes from analysis of observed galaxy images, colours and spectra. However, to quantify evolution, one must compare populations at different redshifts, so it is necessary to correct the observed colours or spectra to a common rest-frame waveband. This is known as applying the $k+e$ correction. In principle, this correction depends on the formation history of the stellar population that contributes to the observed light, so, in principle, it is different for each object. Galaxy samples are now deep enough that $k+e$ corrections are necessary, and sample sizes are large enough that $k+e$ corrections can be a significant source of systematic error. Since determining the appropriate $k+e$ correction for an object boils down to inferring something about its star formation history, $k+e$ correction codes naturally also return an estimate of the stellar mass of the object.

Hence, the conversion from apparent brightness $\ell$ to luminosity $L$ in a fixed rest-frame waveband depends on the $k+e$ correction, so converting from $\ell$ to stellar mass $M_*$ does not require much more work. This is one of the reasons why, over the last decade and a half, the emphasis has shifted to presenting results in terms of stellar mass $M_*$ rather than luminosity $L$. Unfortunately, the conversion from $\ell$ to $L$ or $M_*$ is rather sensitive to how the stellar population was modelled – what we will loosely refer to as the stellar mass-to-light ratio $M_*/L$ in what follows – for which there is still no consensus (PEGASE; Fioc & Rocca-Volmerange 1997, Fioc & Rocca-Volmerange 1999; Bruzual & Charlot 2003; Blanton & Roweis 2007; Conroy, Gunn & White 2009; Maraston et al. 2009).

During the same period, improvements in survey photometry have driven the development of more sophisticated algorithms for estimating the observed flux $\ell$. A few years ago, Bernardi et al. (2010) noted that improving on Sloan Digital Sky Survey (SDSS) pipeline photometry (sky subtraction issues, especially in crowded fields – see e.g. Bernardi et al. 2007; Hyde & Bernardi 2009; Blanton et al. 2011) indicated that the most luminous galaxies at $z \sim 0.1$ were more abundant than expected from the most commonly used

<footnote n="1">¹Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA</footnote>

<footnote n="2">²The Abdus Salam International Center for Theoretical Physics, Strada Costiera 11, I-34151 Trieste, Italy</footnote>

<footnote n="3">³GEPI, Observatoire de Paris, CNRS, Univ. Paris Diderot, Place Jules Janssen, F-92190 Meudon, France</footnote>

<footnote n="4">⁴Institute of Cosmology and Gravitation, Dennis Sciama Building, Burnaby Road, Portsmouth PO1 3FX, UK</footnote>

<footnote n="5">⁵School of Physics and Astronomy, University of Southampton, Southampton SO17 1BJ, UK</footnote>

<footnote>* E-mail: bernardm@sas.upenn.edu</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Comparison of six $\phi(M_*)$ estimates in which $M_*/L$ is fixed to the value used by B13, and the measured apparent brightness $\ell$, from which $L$ is derived, is varied. To compare a large dynamic range, each $M_*\phi(M_*)$ is shown in units of $[M_*\phi(M_*)]_{\text{fid}}$, where the fiducial value uses B13-Meert15 SerExp photometry. All six estimates are in remarkable agreement below $10^{11} M_{\odot}$; above this mass, the *cModel* estimate becomes increasingly smaller than the others.</figcaption></figure>

B13), but $M_*\phi/(M_*\phi)_{\text{fid}}$, where $(M_*\phi)_{\text{fid}}$ uses the SerExp $\ell$ values of B13-Meert15. Three other choices of $\ell$ are also based on PyMorph: one is the B13-Meert15 Sérsic estimate, and the other two show the result of truncating the fitted (SerExp and Sérsic) profiles at $7.5 \times$ the half-light radii when estimating the total magnitude [Fisher et al. (2017) show that this closely approximates the SDSS algorithm]. The two other curves are based on the *Simard11* Sérsic and DR9 *cModel* estimates. (Replacing *cModel* with *Model* magnitudes gives almost identical results.) The trends here are similar to those shown in fig. 1 of B13, except that there, the different $\ell$ estimates were compared, whereas here it is the impact on $\phi(M_*)$ that is shown.

All six estimates are in remarkable agreement below $10^{11} M_{\odot}$. They become increasingly different above this mass. B13 highlighted the fact that these systematic differences give rise to discrepancies in $\phi(M_*)$ determinations of more than 1 dex at $\ge 10^{12} M_{\odot}$. However, they also noted that estimates based on fitting Sérsic or SerExp profiles were in better agreement with one another than with *cModel* magnitudes. What is the real level of disagreement between these other estimates?

We begin with a brief discussion of the fact that the B13-Meert values are based on integrating the fitted profile to infinity, whereas the SDSS truncates. Briefly (see Fischer et al. 2017 for a more detailed discussion), the effect of SDSS-like truncation is largest for single Sérsic fits with large $n$ (~0.16 mag for $n=8$). There is a correlation between $L$ and $n$ – large $L$ have larger $n$ – so truncation matters more for large $L$. However, the $n-L$ correlation is weak; at the largest $L$, the average $n$ is 6 for which the correction is ~0.12 mag. But because there is substantial scatter around the mean $n$, the net effect of truncation is about half than one would naively have expected. Of course, truncation matters even less for the SerExp fits (see figs 6 and 7 in Fischer et al. 2017).

Comparison of the black and grey curves in Fig. 1 shows that truncation only matters above $\log_{10}(M_*/M_{\odot}) \ge 11.5$, where it introduces a systematic difference of 0.1 dex for SerExp, and slightly more for Sérsic. There is general agreement that the truncation scale should be of the order of $7-10 \times$ the half-light radii D'Souza et al. (2015) use the lower value; Kelvin et al. (2012) use the higher

value; the SDSS is in between]. Therefore, the actual uncertainty associated with truncation is smaller than the 0.1 dex shown in Fig. 1. This is much less than the difference with respect to SDSS-*cModel*, which is our main point. Henceforth, we ignore the effects of truncation: all PyMorph values we show are untruncated, since these are the values that are publicly available.

Discussion of the differences between Sérsic or SerExp fits is complicated by the fact that B13-Meert15 and *Simard11* are not in perfect agreement. This is in part because the B13-Meert15 estimates are biased slightly bright (Meert et al. 2013; Bernardi et al. 2014). This difference is mainly driven by the fact that the sky values assumed by *Simard11* are closer to the SDSS pipeline values (which are biased, see Section 2.1) than to PyMorph (see Fischer et al. 2017). The sky background difference leads to a median difference in the fitted Sérsic magnitude of about 0.05–0.3 mag (Meert et al. 2013; Fischer et al. 2017).

Potentially more important is the fact that recent work does not use the *Simard11* photometric reductions ‘as is’. Namely, Thanjavur et al. (2016) advocate using what they call ‘fiducial’ values: they choose the *Simard11* Sérsic value for most objects, and the *Simard11* deVExp value (note that they do not advocate using the *Simard11* SerExp fits at all) for the others. (This choice depends on the quality of the fit to the surface brightness profile.) Although we do not have their ‘fiducial’ values, the solid curve in the top panel of their fig. B1 lies slightly below the dashed curve, indicating that their ‘fiducial’ values produce slightly smaller $M_*\phi(M_*)$ than one obtains if one simply uses the *Simard11* Sérsic values. Note that the difference between these two only matters at large $M_*$, and is approximately the same as the amount by which the blue symbols in our Fig. 1 differ from zero. This is reassuring, since we believe that the B13-Meert15 SerExp $\ell$ values are our least biased estimates. That is to say, results based on what each group believes to be its best estimate (B13-Meert15 SerExp and Thanjavur et al. (2016) ‘fiducial’) are in even better agreement with one another than the 0.1 dex difference between the B13-Meert15 SerExp- and the *Simard11* Sérsic-based results shown in Fig. 1.

Therefore, the most serious systematic difference is driven by the fact that *cModel* magnitudes are much fainter than the others. While B13 argue that this is a significant bias, recent work has questioned this conclusion. For example, D'Souza et al. (2015) do find a bias, but since their analysis of $\phi(M_*)$ resulted in fewer objects at high masses, they conclude that this bias was smaller than reported by B13. However, their conclusions are not based on a direct comparison of the two sets of photometry. As B13 emphasized, when using stellar masses to draw conclusions about photometry, it is important to ensure that differences in $M_*/L$ are not playing a role. The next subsection shows that, once differences in $M_*/L$ have been accounted for, the D'Souza et al. (2015) corrections to SDSS pipeline photometry are in remarkably good agreement with those of B13-Meert15.

## 2.3 Comparison with D'Souza et al. (2015): $\ell$

Recently, D'Souza et al. (2015) have used fits to stacked images to calibrate correction factors that they then apply to the SDSS *Model* magnitudes of individual galaxies. They perform these stacks separately for a number of different bins in stellar mass, where $M_*$ was obtained by combining the MPA-JHU $M_*/L$ values of the objects with the *Model* magnitudes (because they could not stack on the corrected magnitudes, of course).

Fig. 2 shows the median correction factors shown in their fig. 4, which are appropriate for the highest mass bins (redshift
<figure><img src="image_3.png" /></figure>

**Figure 2.** Corrections to SDSS pipeline photometry in the format of fig. 4 of D'Souza et al. (2015). Filled symbols show the median correction applied by D'Souza et al. (2015) to the most massive objects, and black solid line shows that associated with the B13-Meert15 SerExp magnitudes for the same objects, both shown as a function of the $M_*$ estimate associated with MPA-JHU $M_*/L$ values (scaled to a Chabrier IMF) and SDSS Model magnitudes. Black dotted curves show the region that contains 68 per cent of the objects. Since the D'Souza et al. (2015) corrections were computed as a function of MPA-JHU $M_*/L$ values and SDSS Model magnitudes, this is the correct way to do the comparison, and the agreement is remarkable. Red curves show Model-SerExp as a function of the $M_*$ estimate associated with the same MPA-JHU $M_*/L$ values but SerExp magnitudes. When shown this way, the implied corrections are larger, but comparing with the filled symbols, as D'Souza et al. (2015) did, is not justified.

and concentration index) they consider – the ones that are most relevant to our present study. Notice that their values are brighter than the Model magnitudes, in qualitative agreement with B13. To see if the agreement is quantitative, the black solid curve shows the corresponding median values of the differences between the B13-Meert15 SerExp and Model magnitudes for the same objects. The agreement is remarkable.

Why then did D'Souza et al. (2015) come to a different conclusion? The reason is that their correction factor was calibrated in bins of $M_*$ in which Model magnitudes were used to convert $M_*/L$ to $M_*$. However, when they discussed the B13-Meert15 SerExp correction, they appear to be referring to some of the curves in fig. 1 of B13, which are binned in Ser magnitudes. Since Ser is not the same as Model, the objects in each $M_*$ bin are no longer the same, and so their discussion is not based on a direct comparison. To illustrate this point, the red solid line shows the median Model-SerExp correction as a function of SerExp-rather than Model-based $M_*$. Now the correction appears larger, but this is simply a consequence of the well-known fact that there is scatter between SerExp and Model magnitudes (fig. 1 in B13). D'Souza et al. (2015) were incorrectly comparing the 0.3 mag correction associated with the blue symbols here with the 0.5 mag correction factor suggested by the red curve. The correct comparison is with the black curve, and this shows that there is little or no difference.

The agreement is not quite as good as it seems, because D'Souza et al. (2015) obtained their corrections from truncated profiles: the inner component is truncated outside a radius equal to 7 half-light radii, while the outer component extends to infinity. In contrast, B13-Meert15 always integrate to infinity. If we were to account for this difference (see fig. 7 in Fischer et al. 2017), then the D'Souza et al. (2015) correction to Model magnitudes would exceed that of

<figure><img src="image_4.png" /></figure>

**Figure 3.** Same as the previous figure but now for all galaxies in the mass bin.

B13-Meert15 (by about 0.05 mag; see also Fig. 1 here and related discussion). That is, if anything, D'Souza et al. (2015) is brighter than B13-Meert15, not fainter. This is perhaps not surprising; Py-Morph SerExp forces the second component to be an exponential, whereas D'Souza et al. (2015) allow it to be another Sérsic, potentially including more light at large radii.

It is worth noting that, while the agreement between B13-Meert15 and Thanjavur et al. (2016) described in the previous section is perhaps not so surprising – both are based on similar approaches to fitting the surface brightness profiles of individual galaxies, so they differ in the details of how the fit is actually done – the agreement between the B13-Meert15 SerExp and the D'Souza et al. (2015) corrections shown in Fig. 2 really is remarkable. For example, the stacking method must account for light lost to masked pixels, which are different for each member of the stack, and this will become increasingly important (and difficult!) at the high-mass end. This potential systematic is not an issue for approaches based on individual galaxies. Another potential systematic is the sky, but we have already argued that because D'Souza et al. (2015) used the Blanton et al. (2011), rather than the SDSS pipeline, values for the sky, their sky values are in good agreement with those of B13-Meert15.

We now consider the scatter around the median. The dashed lines show that this scatter (the region that encloses 68 per cent of the objects at each $M_*$) is slightly asymmetric, with an rms of approximately ~0.2 mag above $10^{11.3} M_{\odot}$. It is similar to that shown in fig. 2 of D'Souza et al. (2015). This rms is substantially larger than the quoted errors on either Model or SerExp Meert et al. (2013) quote errors of approximately 0.05 mag each; adding in quadrature yields 0.07 mag, so we conclude that most of it is intrinsic (because $0.07^2 \ll 0.2^2$). We return to this point in Section 4.2.

Fig. 3 shows a similar analysis of all galaxies – not just those with high concentration indices. Comparison with the previous panel shows that the median corrections here are smaller than they were for the high-concentration galaxies. (This is consistent with the fact that massive high-concentration galaxies tend to be in crowded fields.) Thus, for the galaxies with the biggest corrections, the median correction is the same as or smaller than made by D'Souza et al. (2015).

Based on the arguments above, we conclude that in fact there are no significant differences between the photometry of D'Souza et al. (2015) and B13. Moreover, the agreement between the solid black line and the big blue symbols in Fig. 2 indicates that cModel is indeed biased by the amount claimed by B13. Therefore, the fact
that cModel is an outlier in Fig. 1 should not feature in discussions of the impact of photometry-related systematics on $\phi(M_*)$. The real level of systematic uncertainty is ~0.1 dex (see the other curves in Fig. 1).

## 3 DEPENDENCE ON STELLAR POPULATION MODELLING

Having shown that systematics with more recent photometry contribute about 0.1 dex uncertainty at high masses (and that SDSS pipeline photometry is substantially fainter, in agreement with B13), we now study the systematics that arise from different assumptions about the stellar populations of the most massive galaxies. Since we would like to separate the effects of stellar population modelling from those due to photometry, we must first account for the fact that not all the groups whose $M_*$ estimates we compare below used the same $\ell$ for estimating $L$. To account for this, we multiply each group's $M_*$ value by $\ell_{\text{fiducial}}/\ell$ where, for $\ell_{\text{fiducial}}$, we use the r-band SerExp values provided by Meert et al. (2015). Note that by using the ratio of the apparent brightnesses, we are making no change to the (stellar population dependent) $k + e$ values assigned by each group. (We do account for the small differences that arise from the fact that not all groups used our fiducial cosmological parameters.)

We have taken a number of stellar mass estimates from the SDSS data base, www.sdss.org/dr12/spectro/galaxy/, which all assumed a Kroupa (2001) initial mass function (IMF). These are as follows:

(i) spectrophotometric-based estimates from the Portsmouth group (Maraston et al. 2013), which assume a passive stellar population (dust-free) that is expected to be appropriate for the most massive galaxies (table stellarMassPassivePort);

(ii) PCA-based estimates from the Wisconsin group (Chen et al. 2012) that are based on the stellar population synthesis (SPS) models of Maraston & Strömbäck (2011) and account for dust [table stellarMassPCAWiscM11 – for the high masses of most interest in our work, these are not very different from the models of Bruzual & Charlot (2003)];

(iii) estimates from the Granada group, which are based on the publicly available Flexible Stellar Population Synthesis (FSPS) code (Conroy et al. 2009), and assume that star formation occurs over a wide range of redshifts (tables stellarMassFSPSGranWide-Dust and stellarMassFSPSGranWideNoDust);

(iv) estimates from the MPA-JHU group (table galSpecExtra), which are based on the galSpec tool (building on the work of Kauffmann et al. 2003; Brinchmann et al. 2004; Tremonti et al. 2004).

There is no consensus on whether or not massive galaxies are dusty. The Portsmouth-Passive estimates assume they are not. This is based on fig. 8 of Thomas et al. (2013) who found very little dust in massive galaxies. On the other side, the MPA-JHU and Wisconsin estimates assume they are. The GranadaWide group allows for both possibilities.

We have also included older estimates from Blanton & Roweis (2007) and a particularly simple estimate used by Bernardi et al. (2010) and B13, in which $M_*/L$ is a linear function of colour [see equation 6 in Bernardi et al. (2010), which is motivated by Bell et al. (2003) who used the PEGASE models], primarily for comparison with previous work. Finally, we include more recent estimates from Mendel et al. (2014, their tables 3 and 5); they provide dusty and dust-free $M_*/L$ estimates, and do not allow for multiple bursts.

When computing stellar masses, there is also no consensus on whether or not one should allow for star formation histories with multiple bursts. The MPA-JHU estimates assume bursty histories,

whereas the Portsmouth-Passive and Mendel et al. (2014) models do not. Mendel et al. (2014) exclude bursts because Gallazzi & Bell (2009) showed that including bursty SPS models can lead to a systematic underestimate of $M/L$ by up to 0.1 dex for galaxies whose star formation histories are generally smooth. This is likely to be the case for the most massive galaxies.

As we show in Section 4.3.1, dust and burstiness contribute significantly to the error budget at the high-mass end of the stellar mass function (e.g. Moustakas et al. 2013; Mendel et al. 2014).

There is also no consensus on the shape of the IMF. We have (arbitrarily) chosen to correct all $M_*$ values to a Chabrier IMF (Chabrier 2003). No correction is necessary for the Blanton–Roweis, B13 and Mendel14 estimates. However, the Portsmouth, Wisconsin, Granada and MPA-JHU estimates assumed a Kroupa IMF (Kroupa 2001); to convert to Chabrier, we simply decrease the quoted (Kroupa) $M_*$ values by 0.05 dex (see e.g. table 2 in Bernardi et al. 2010).

Appendix A compares the $M_*/L$ estimates from these different analyses. It shows that, above $10^{11} M_\odot$, these differences are of the order of 0.05 dex or larger, provided we compare dusty or dust-free models separately (but not dusty with dust-free). However, there is no consensus on the effects of dust. For Mendel14, allowing for dust decreases the inferred $M_*/L$ by 0.05 dex or more (Fig. A1), but for WiscM11 allowing for dust increases it by about 0.06 dex (fig. 13 in Chen et al. 2012). The impact of dust on the GranadaWide $M_*/L$ values is more complicated (see Appendix A for details). The next subsection shows how these differences impact $\phi(M_*)$.

### 3.1 Effect of stellar population model

Fig. 4 compares a number of estimates of $\phi(M_*)$ that are based on the publicly available $M_*$ values listed earlier. Below $10^{10.5} M_\odot$, most estimates of $\phi(M_*)$ differ by of the order of 0.1 dex (we argue later that this should not be used to argue that $M_*/L$ in low-mass galaxies is well understood!). The Portsmouth-Passive estimate should be compared only at higher masses $M_* > 10^{11} M_\odot$ because a passive stellar population is not appropriate at smaller masses. Using their star formation model (table stellarMassStarforming-Port) lowers the value of $\phi(M_*)$ (we have not shown this in the figure because it is not a good choice for the high masses of most interest in our work).

At larger masses, the differences can be much larger. To show this more clearly, Fig. 5 presents these curves normalized by $\phi_{\text{fid}}$, for which we use the one associated with the B13-Meert15 SerExp photometry and the $M_*/L$ value of Mendel14. The results are bracketed by the curves associated with the Blanton–Roweis (lowest) and B13 (uppermost) $M_*/L$ values. Whereas the Blanton–Roweis values are known to be inappropriate for the most massive galaxies (see discussion in Bernardi et al. 2010), the B13 prescription is slightly ad hoc (e.g. the $k + e$ corrections and $M_*/L$ values are not derived self-consistently). Nevertheless, even for the better motivated models, these are not small offsets ($\ge 0.5$ dex or more at $10^{12} M_\odot$), especially in view of the fact that these massive galaxies are expected to have the simplest (single-burst) stellar populations.

### 3.2 Effect of dust

Some of these differences arise from making different assumptions about dust in these galaxies. Fig. 6 compares $\phi(M_*)$ associated with the dusty and dust-free $M_*/L$ values of Mendel et al. (2014); the difference is entirely a consequence of the differences shown in Fig. A1. The dust-free models have larger $M_*/L$, so they result in larger $\phi(M_*)$; at $10^{12} M_\odot$, the difference is about 0.4 dex. [For ease
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Comparison of a number of different determinations of $\phi(M_*)$ as labelled. In all cases, the $M_*$ values were scaled to a Chabrier IMF and B13-Meert15 SerExp photometry, so the differences are entirely due to the stellar population modelling, fitting and assumptions about dust in the galaxies. The Portsmouth-Passive estimate should be compared only at higher masses $M_* > 10^{11} M_{\odot}$ because a passive stellar population is not appropriate at smaller masses.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Same as the previous figure, but now shown in units of $\phi_{\text{fid}}$ for which we have used the cyan curve in the previous figure (i.e. SerExp photometry from B13-Meert15 and $M_*/L$ from dust-free Mendel14). At $M_* > 10^{11} M_{\odot}$ the B13-based estimate is the largest, and the Blanton-Roweis is smallest. The Wisconsin, GranadaWide and Mendel14 determinations are within about 0.1 dex of one another; they lie above the Portsmouth-Passive and MPA-JHU values. The Portsmouth-Passive estimate should be compared only at higher masses $M_* > 10^{11} M_{\odot}$ because a passive stellar population is not appropriate at smaller masses.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Comparison of $\phi(M_*)$ associated with the dust-free and dusty $M_*/L$ values of Mendel et al. (2014). Only in this figure, $M_*$ is computed using Simard11 Sérsic photometry as this eases comparison with fig. B1 in Thanjavur et al. (2016).</figcaption></figure>

of comparison with fig. B1 in Thanjavur et al. (2016), in this figure $M_*$ is computed using Simard11 Sérsic photometry.]

To see if this level of discrepancy is the same in other models, we have compared dusty and dust-free estimates of $\phi(M_*)$ in the GranadaWide and Wisconsin models as well. Fig. 7 shows the results: solid curves represent analyses that assume that the galaxies

are dust-free; dotted curves assume that the galaxies are dusty. [In fact, only dusty WiscM11 $M_*$ values are available. However, fig. 13 of Chen et al. (2012) shows that assuming that the galaxies are dust-free decreases the WiscM11 $M_*/L$ values by about 0.06 dex. So, the WiscM11-NoDust $\phi(M_*)$ curves in Fig. 7 were approximated by simply shifting their $M_*$ values downwards by this amount.]

We noted before that there is not even consensus on the sign of the effect, and this is reflected in the $\phi(M_*)$ estimates. In the Wisconsin models, allowing for dust increases the inferred $M_*/L$
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Same as Fig. 5, but now comparing $\phi(M_*)$ estimates in which galaxies are assumed to be dusty (dotted) or dust-free (solid). Since only dusty WiscM11 $M_*$ values are available, the solid curve for this case is based on making the $M_*/L$ values 0.06 dex smaller (see the text for details).</figcaption></figure>

and hence $\phi(M_*)$; for Mendel14, accounting for dust reduces $\phi(M_*)$ by a similar amount. As a result, the dust-free Mendel14 estimate happens to be reasonably like the WiscM11 dusty estimate, and vice versa. On the other hand, dust has little effect on $\phi(M_*)$ of the GranadaWide models (the effect on $M_*/L$ is actually more complicated; see Appendix A).

Before moving on, it is worth noting that the difference between the B13 $\phi(M_*)$ and the dust-free Mendel14 curve is smaller than the ~0.3 dex reported in Thanjavur et al. (2016). This is because their analysis was based on the dusty Mendel14 models – the dotted cyan curve in Fig. 7. This shows that the differences they report are largely due to differences in $M_*/L$, which are largely due to different assumptions about dust in these galaxies.

Until we know whether dusty or dust-free models are correct, we thought it prudent to show the level of systematics on $\phi(M_*)$ if we compare the two sets of models separately. This is done in Fig. 8. The systematic differences on $\phi(M_*)$ between dust-free models are slightly smaller than for dusty models: $\le 0.3$ dex compared to ~0.4 dex at $10^{12} M_{\odot}$.

### 3.3 Effect of measurement errors

So far we have studied how the measured $\phi(M_*)$ depends on $M_*/L$ and $\ell$. However, the estimated $M_*$ values are noisy, with statistical (not systematic!) uncertainties of ~0.1 dex. While these matter little at the lower masses where $\phi(M_*)$ is relatively flat, they modify $\phi(M_*)$ substantially at the high-mass end of most interest here. If the intrinsic $M_*\phi_{int}(M_*)$ is approximately flat at low masses and falls as $\exp[-(M/M_*)^\beta]$ at high masses, then errors that are Gaussian in $\log_{10}(M_*)$ make

$$ \frac{\phi_{\text{obs}}}{\phi_{\text{int}}} - 1 \approx \frac{(\beta \ln(10) \sigma_{\log_{10} M})^2}{2} \left( \frac{M_*}{M_{*-\text{char}}} \right)^{\beta} \times \left[ \left( \frac{M_*}{M_{*-\text{char}}} \right)^{\beta} - 1 \right] \quad (1) $$

(equations 10 and 11 of Bernardi et al. 2010).

Notice that errors matter little below a characteristic mass $M_{*-\text{char}}$; this is not surprising, since the flatness of $M_*\phi_{int}(M_*)$ below $M_{*-\text{char}}$ means that errors move objects between mass bins, but because all bins have the same height, there is no net change to $M_*\phi_{int}(M_*)$.

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Same as the previous figure, but now comparing dust-free models (top) and dusty models (bottom) separately.</figcaption></figure>

Fits that constrain $\beta = 1$ tend to return $M_* \sim 10^{11} M_{\odot}$ (e.g. D'Souza et al. 2015; Thanjavur et al. 2016), so errors only matter at higher masses. (For example, at $10^{12} M_{\odot}$, errors of the order of $\sigma_{\log_{10} M} = 0.1$ dex make $\log_{10}(\phi_{\text{obs}}/\phi_{\text{int}} - 1) \approx 0.4$ dex.) This means that the 0.1 dex agreement between various determinations of $\phi(M_*)$ at lower masses should not be used to argue that $M_*/L$ in low-mass galaxies is well understood (e.g. Thanjavur et al. 2016). As Appendix A shows, this agreement is hiding large systematic differences between different groups.

Fig. 9 shows the measured (solid) and intrinsic (dashed) stellar mass functions for the dusty and dust-free models of Mendel14 (the reference photometry is SerExp of B13-Meert15, so comparing the solid curves here with those in Fig. 6 shows that changing from Simard11 Sérsic photometry to B13-Meert15 SerExp is small). The dashed lines show the ‘error-corrected’ measurements obtained using a modified $V_{\text{max}}^{-1}$ method that accounts for errors (see Sheth 2017 for details). We assumed measurement errors of 0.1 dex in all cases. Clearly, measurement errors matter little below $10^{11} M_{\odot}$, but enhance $\phi(M_*)$ by more than 0.3 dex at ~$10^{12} M_{\odot}$, consistent with the analysis above. Only at these rather large masses are the (statistical) measurement errors comparable to the systematic differences that arise from different treatments of the stellar population. Note that this is of the order of the mass scale where the finite size of the survey volume is also beginning to matter.
<figure><img src="image_11.png" /><figcaption>**Figure 11.** Comparison of $\phi(M_*)$ estimates based on dusty (solid blue) and dust-free (solid cyan) $M_*/L$ ratios of Mendel et al. (2014), and of B13 (solid black); all three curves used B13-Meert15 SerExp photometry. Dashed black curve represents B13 $M_*/L$ and B13-Meert15 Sérsic photometry. The ~0.15 dex difference between the solid cyan and black curves – and not that between the blue and dashed black curves – is a better indicator of the true level of systematic difference between B13-Meert15 and Thanjavur et al. (2016) at high $M_*$.</figcaption></figure>

The lowest (magenta solid) curve in Fig. 12 shows $\phi(M_*)$ reported by Li & White (2009). We only show it because it appears in D'Souza et al. (2015). Li & White used Petrosian luminosities (known even prior to 2009 to underestimate $\ell$ of massive galaxies) and $M_*/L$ ratios from Blanton & Roweis (2007, a choice Blanton and Roweis themselves cautioned against for the most massive galaxies; also see fig. 22 in Bernardi et al. 2010 and related discussion), so it is not surprising that it is the lowest. Replacing Petrosian with cModel DR7 magnitudes yields the dotted grey curve; it too lies substantially below all the others in the figure. Exchanging the Blanton–Roweis $M_*/L$ for MPA-JHU values increases the cModel curve slightly more (dotted green curve), but the difference compared to the Blanton–Roweis curve (magenta) is more than 1 dex at $M_* \ge 10^{12} M_\odot$.

The black solid uppermost curve comes from SerExp luminosities and the colour-based prescription for $M_*/L$ of B13. Using the MPA-JHU $M_*/L$ values instead – as this was the choice made by D'Souza et al. (2015) – yields the curve labelled B13-Meert15 (solid green). This green curve is in rather good agreement with the solid orange curve, which is our rendition of D'Souza et al.'s (2015) determination. Since the only difference between the black and green curves is $M_*/L$, the agreement between the green and orange curves shows that when the comparison is done using the same $M_*/L$, then D'Souza et al. (2015) are in good agreement with B13. This is consistent with our finding (Fig. 2) that the median corrections to SDSS pipeline photometry were similar.

Strictly speaking, D'Souza et al. (2015) used a slightly different background cosmology ($\Omega_m = 0.25$, $\Omega_\Lambda = 0.75$). We expect this to cause negligible difference because, at $z \sim 0.15$, the luminosity distance in their cosmology is larger than ours by a factor of 717.2/713.5, so their luminosities are about 1 per cent brighter than ours. We checked this effect by scaling the B13 sample to D'Souza et al.'s (2015) cosmology and recomputing $\phi(M_*)$, finding it to be nearly indistinguishable from the original. Hence, we expect this would also be true if we instead scaled D'Souza et al. (2015) to our cosmology. Therefore, we have shown the determination of $\phi(M_*)$

by D'Souza et al. (2015) without applying any changes to their cosmology.

However, D'Souza et al. (2015) actually only report the intrinsic, error-corrected dashed orange curve (see their table 2), so we have had to do some work to produce the error-broadened solid orange curve. Each correction is small, so it is easier to see their effects on $\phi(M_*)$ when normalized by a fiducial value, which we take to be that for MPA-JHU $M_*/L$ and B13-Meert15 SerExp $\ell$ values (solid green line). The steps we took are as follows.

(i) The dashed orange curve in Fig. 13 shows D'Souza et al.'s (2015) intrinsic (error-corrected) curve. (All the orange curves in this figure were scaled to a Chabrier IMF by applying the 0.05 dex offset between the Kroupa and Chabrier IMFs.)

(ii) The dotted orange curve shows the result of accounting for errors in the $M_*/L$ determination by using the 0.1 dex they report in our equation (1).

The difference between the dotted orange curve and the solid black curve is the discrepancy with B13 that they report. The difference in Fig. 13 is actually larger than that reported in fig. 7 of D'Souza et al. (2015) because they ignored the offset between the Kroupa and Chabrier IMFs. While this is of the order of 0.3 dex at $10^{11.5} M_\odot$, the appropriate comparison is with the fiducial curve, for which the differences are larger than 0.3 dex only above $10^{12} M_\odot$.

(iii) There is another systematic difference between the D'Souza et al. (2015) photometry and the more traditional estimates of B13 (and Simard11). For a given class of objects, D'Souza et al. (2015) use a single correction factor (based on their stacked images). This does not include the fact that there is scatter around the median correction. The B13 curves do include this scatter, so a fair comparison requires that we account for its effects. In our discussion of Fig. 2, we argued that most of the measured scatter is intrinsic – there is no reason why the profile of each galaxy in the stack should be exactly the same – so we have accounted for an additional 0.2/2.5 dex correction to the orange dotted curve using our equation (1). This produces the solid orange curve.

(iv) Finally, contrary to what is written in the caption of fig. 4 of D'Souza et al. (2015), they actually used the mean instead of median correction. The median is slightly larger – using it instead yields the triple-dot-dashed line in Fig. 13. This is a small effect, which is why, in Fig. 12, we treat the orange solid curve as the final corrected value.

Clearly, once all of these small corrections have been made, the discrepancy with respect to B13-Meert SerExp is within about 0.1 dex (when using the same $M_*/L$). This is much smaller than D'Souza et al. (2015) originally claimed was due to photometry.

Before moving on, we think that it is important to note that below about $10^{10.5} M_\odot$ D'Souza et al. (2015) lie about 0.15 dex above the fiducial value, which is based on B13-Meert15 SerExp photometry. We have argued that this photometry is likely very close to the ‘fiducial’ photometry of Thanjavur et al. (2016). Once differences in $M_*/L$ have been accounted for, the fiducial curve is similar to that of Thanjavur et al. (2016): in fact, Fig. 11 suggests that it is already slightly above the Thanjavur et al. (2016) determination, so D'Souza et al. (2015) would be even higher. This offset is also present in their analysis of the luminosity function: fig. 8 in D'Souza et al. (2015) shows that they lie about 0.2 dex above the $\phi(L)$ determination of B13. One cannot appeal to $M_*/L$ differences to explain the offset in $\phi(L)$, so it must be due to other effects. Since B13 were careful to establish that, at the faint end, their luminosity function was in excellent agreement with the $\phi(L)$ estimates of several other
<figure><img src="image_12.png" /><figcaption>**Figure 12.** Comparison of a number of different determinations of $\phi(M_*)$ as labelled. The curves of most interest here are the ones based on MPA-JHU $M_*/L$ (green lines), so differences between them are entirely due to how $\ell$ was estimated. The orange curves are the determination of $\phi(M_*)$ by D'Souza et al. (2015) obtained using the MPA-JHU $M_*/L$ values – measured (solid) and intrinsic (dashed). The B13 SerExp estimate (black solid line) and the two Blanton–Roweis-based curves are shown as reference.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 13.** Similar to the previous figure, but now shown in units of $\phi_{\text{fid}}$ (solid green line) for which we combined the MPA-JHU $M_*/L$ ratios with the B13-Meert15 SerExp apparent brightnesses. The D'Souza et al. (2015) and B13 analyses agree, to within 0.1 dex, when the comparison is done using the same $M_*/L$.</figcaption></figure>

groups, the vertical offset between D'Souza et al. (2015) and B13 is surprising. The reason for this is an open question.

## 4.3 Comparison with Moustakas et al. (2013) and summary

We have described how a variety of systematics affect determination of the stellar mass function. To put these results in perspective, Fig. 14 shows the state of the art before and after B13. Fig. 15 shows the same curves normalized by the counts associated with

the Mendel14 dust-free $M_*/L$ ratios and the B13-Meert15 SerExp apparent brightnesses (solid cyan curve).

The two bottom-most estimates are based on SDSS pipeline determinations of $\ell$. The magenta solid curve shows the Li & White (2009) estimate; we have already noted that it was biased low because of inappropriate choices for both photometry and stellar population modelling. The brown crosses show the estimate of Moustakas et al. (2013) that uses slightly better (but still SDSS pipeline) cModel photometry $\ell$ and a different stellar population model (about which more in Section 4.3.1).

The other curves are all based on more recent determinations of $\ell$, and they all result in more high-mass objects. The solid orange and green curves show the result of using the same dusty and bursty MPA-JHU $M_*/L$ estimate, but with more sophisticated estimates of the apparent brightness [the stacking analysis of D'Souza et al. (2015) and SerExp values of B13-Meert15, respectively]. Both curves predict substantially more stellar mass than Moustakas et al. (2013) and Li & White: the increase is of the order of 1 dex and 2 dex at $10^{12} M_{\odot}$, respectively. The small differences between these two curves (orange and green) show that photometry now contributes only of the order of 0.1 dex to the systematics error budget. That we are now discussing 0.1 dex systematics rather than 1 or 2 dex represents real progress in treatments of the photometry of the most massive galaxies.

The other curves show results based on alternative treatments of the stellar populations. The cyan and blue curves show the systematic effects of dust in the single-burst Mendel14 models: the cyan curves assume that galaxies are dust-free, suggesting that, within the context of the same modelling and fitting framework, dust leads to differences up to 0.3 dex at $10^{12} M_{\odot}$. However, there is no agreement even on the sign of the systematic. The blue (Mendel14)
<figure><img src="image_14.png" /><figcaption>**Figure 14.** Representative $\phi(M_*)$ estimates before (Li & White 2009; Moustakas et al. 2013) and after B13 as labelled. All curves and symbols only show error-broadened values.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 15.** Similar to the previous figure, but now shown in units of $\phi_{\text{fid}}$ (solid cyan line) for which we combined the Mendel14 dust-free $M_*/L$ ratios with the B13-Meert15 SerExp apparent brightnesses.</figcaption></figure>

and red (WiscM11) curves are both dusty models – so differences between frameworks lead to $\ge 0.3$ dex differences already at $M_* \ge 10^{11.5} M_{\odot}$ – but the dust-free version of the red curve would bring it down, close to the blue curve, rather than increase it. Thus, we conclude that uncertainties in $M_*/L$ are at least 0.3 dex; Fig. 15 in particular shows that they are substantially larger than photometry-related systematics.

The uppermost curve shows the estimate of B13. Some of this can be traced to the fact that its calibration of $M_*/L$ (from colour) is based on the PEGASE stellar population models, and these tend to produce the largest $\phi(M_*)$ estimates (see the right-hand panel of fig. 19 in Moustakas et al. 2013). In this context, it is interesting to revisit the offset between Moustakas et al. (2013) and the others.

### 4.3.1 Effects of bursts

We know that a systematic difference between the Moustakas et al. (2013) determination of $\phi(M_*)$ and that of Mendel et al. (2014) is due to $\ell$. Another is due to the use of bursty (or not) models when determining $M_*/L$. The left-hand panel of fig. 19 in Moustakas et al. (2013) shows that fitting with single-burst models would return $M_*/L$ values that are larger by at least 0.1 dex. Mendel et al. (2014) report a similar difference, and go on to argue that single-burst models are adequate for the vast majority of the galaxy population. As this is even more likely to be true of the most massive galaxies, they do not consider bursty models at all. Therefore, we now explore how the Moustakas et al. (2013) $\phi(M_*)$ determination compares with the others once we have accounted for both $\ell$ and burstiness. We do not have their $M_*$ determinations on an object-by-object basis, so we have transformed their $\phi(M_*)$ estimate in two steps as follows.

The brown symbols in Fig. 16 show their fiducial determination. It happens to be in good agreement with the dotted green curve. This is a coincidence, because, although it too uses the same cModel $\ell$, the $M_*/L$ determination is different (FSPS versus MPA-JHU, both bursty). But we can use this to motivate a simple correction for the effect of $\ell$. Namely, the solid green curve shows how $\phi(M_*)$ changes when the cModel photometry is replaced with SerExp of B13-Meert15. Therefore, we shift each of Moustakas’ brown crosses in the $M_*$ direction by the same amount that the corresponding solid and dotted green curves differ. The magenta diamonds show the result; since the brown crosses were similar to the dotted curve, these diamonds are similar to the solid curve.

To remove the effect of bursts on their $M_*$ estimates, we add an additional 0.1 dex (the value suggested by their fig. 19) to their $M_*$ values: this yields the red triangles. These are now rather close to the determination of B13. Accounting for the difference between the PEGASE models (on which the B13 values were calibrated) and the
<figure><img src="image_16.png" /><figcaption>**Figure 16.** Effect on the $\phi(M_*)$ of Moustakas et al. (2013) when $\ell$ is increased from cModel (brown crosses) to B13-Meert15 SerExp (magenta diamonds), and then when bursty models are not included when determining $M_*/L$ (red triangles). Shifting from FSPS to PEGASE models yields an estimate that lies slightly above that of B13 (black solid curve).</figcaption></figure>

FSPS models they used would shift the red triangles by a further 0.05 dex to the right. We conclude that, once differences in $\ell$ and $M_*/L$ have been removed, Moustakas et al. (2013) are in good agreement with the others. But we emphasize that failure to account for the differences in $\ell$ (at least) will lead to a severe underestimate of the true counts at large $M_*$.

## 4.4 Cumulative counts

Fig. 17 shows the cumulative stellar mass density for a selection of our stellar mass density $\rho_*(\ge M_*)$ determinations that summarize our findings. The bottom-most magenta solid curve shows the Li & White (2009) estimate. The two dotted curves just above it show the result of using the most appropriate SDSS pipeline photometry available at the time (grey dotted), as well as a better estimate of $M_*/L$ that was also available (green dotted). These shifts alone account for more than 1 dex at $10^{12} M_{\odot}$. Fig. 14 shows that the dotted green curve will also provide a good description of the results of Moustakas et al. (2013).

The solid orange and green curves show the result of using the same dusty MPA-JHU $M_*/L$ estimate, but with more sophisticated estimates of the apparent brightness than the SDSS photometry (the stacking-based estimates of D'Souza et al. (2015) and SerExp fits of B13-Meert15, respectively). Note that both curves predict substantially more stellar mass than Moustakas et al. (2013) and Li & White: the increase is of the order of 1 dex and 2 dex at $10^{12} M_{\odot}$, respectively. The small differences between the orange and green curves show that photometry now contributes only of the order of 0.1 dex to the systematics error budget.

The reasonable agreement between the blue curve and the orange and green ones is interesting. Whereas the orange and green curves

assumed bursty star formation histories, the blue curve did not. If the single-burst assumption is more realistic at high masses, then we should shift the orange and green curves to the right by at least 0.1 dex (see the discussion in the previous section). This would decrease the difference with respect to the blue curve. While this is a coincidence, because they are based on different stellar population models (MPA-JHU and FSPS), it may indicate that consensus on this issue would significantly reduce systematics on $\phi(M_*)$.

Unfortunately, this is not the full story. The other solid curves show results based on other treatments of the stellar populations. The cyan, blue and red curves correspond to the various assumptions about dust discussed in Fig. 15.

Whereas all these differences are systematic, there is also a statistical effect associated with measurement errors, which broaden the measured curve relative to its intrinsic value. Comparison of the solid and dashed green curves shows the impact of 0.1 dex statistical measurement errors. Applying this difference to the other curves is a reasonable way of quantifying the effect of measurement errors on the other $\phi(M_*)$ determinations.

## 5 CONCLUSIONS

When comparing different determinations of $\phi(M_*)$, it is important to separate effects that are due to differences in estimating the apparent brightness from those that arise from modelling the stellar population: what we referred to as $\ell$ and $M_*/L$. If the algorithm for assigning $M_*/L$ is the same, then different estimates of $\ell$ appear to give rise to order of 1 dex differences in $\phi(M_*)$ at $10^{12} M_{\odot}$ (Fig. 1). However, much of this is driven by the discrepancy between SDSS pipeline photometry and more recent work based on Sérsic or SerExp fits. We showed that these recent analyses agree, to
<figure><img src="image_17.png" /><figcaption>**Figure 17.** Cumulative stellar mass density for the curves shown in Figs 14 and 16. Solid and dashed green curves, which assume MPA-JHU $M_*/L$ values, show the observed error-broadened and intrinsic error-corrected estimates from B13-Meert15. The other curves only show error-broadened values. The D'Souza et al. (2015) (solid orange) and B13 (solid green) analyses agree, to within 0.1 dex, up to $10^{11.6}$ M$_\odot$; both lie more than a factor of 5 above the estimates of Li & White (2009). At larger masses, while the top set of curves differ from one another by a factor of 2 at most (due to different assumptions in the star formation histories or/and the presence of and effects of dust), they both lie about two orders of magnitude above the Li & White estimate.</figcaption></figure>

within 0.1 dex, up to ~$10^{12}$ M$_\odot$ (or higher). Recent claims of larger discrepancies in the literature are primarily due to $M_*/L$ differences.

Once differences in $M_*/L$ have been accounted for, the ‘fiducial’ $\phi(M_*)$ estimate of Thanjavur et al. (2016) is in good agreement with that based on SerExp photometry of B13 and Meert et al. (2015, Fig. 11). This implies that the ‘fiducial’ photometry of Thanjavur et al. (2016) – which is based on the work of Simard11 – is in good agreement with the B13-Meert15 SerExp photometry.

In addition, we showed that the median corrections applied to the SDSS pipeline photometry by D'Souza et al. (2015) are essentially the same as those advocated by B13 (Fig. 2). Comparison of the associated $\phi(M_*)$ is more complicated, because there is substantial scatter around this median (≥0.15 mag), most of which is intrinsic (cf. Fig. 2 and related discussion). This matters most at the luminous end of $\phi(L)$ and hence the high-mass end of $\phi(M_*)$. Object-by-object analyses like B13-Meert15 include this effect trivially, but stacking-based analyses [like D'Souza et al. (2015)] cannot. Adding the effect of this scatter to the D'Souza et al. (2015) $\phi(M_*)$ estimate brings it, too, to within 0.1 dex of the B13-Meert15 SerExp-based estimate (Fig. 13).

Thus, all three groups agree that if the same $M_*/L$ algorithm is used, then the stellar mass density in objects more massive than $10^{11.3}$ M$_\odot$ is at least five times larger than estimates based on SDSS pipeline photometry (Fig. 17), as first noted by B13. That we are now discussing 0.1 dex differences between groups, rather than factors of 5, represents substantial recent progress in the photometry of massive galaxies.

When photometric parameters are held fixed, differences between stellar population treatments result in $\phi(M_*)$ estimates that can differ by more than a factor of 3 at $10^{11.3-10^{12}}$ M$_\odot$ (Figs 4, 5, 14

and 15). These differences arise because of differing assumptions about the presence of and effects of dust (Figs 6 and 8) and of whether or not the star formation histories were bursty (Fig. 16). Hence, it is systematics in modelling the stellar population, and not the photometry, which now dominate the $\phi(M_*)$ error budget. For example, there is currently not even agreement on whether allowing for dust should increase or decrease the estimated $M_*/L$ (Appendix A).

The shape of $\phi(M_*)$ means that errors – statistical or systematic – mainly matter at the massive end (equation 1). This means that there may be general agreement on $\phi(M_*)$ at masses below ~$10^{10.5}$ M$_\odot$ (Fig. 5), even though there is substantial disagreement about $M_*/L$ (Appendix A). Therefore, the 0.1 dex agreement between various determinations of $\phi(M_*)$ at lower masses should not be used to argue that $M_*/L$ in low-mass galaxies is well understood! At higher masses, statistical measurement errors on the value of $M_*$ can have a big impact: 0.1 dex errors give rise to 0.3 dex differences in $\phi(M_*)$ above $10^{12}$ M$_\odot$ (Figs 10 and 17). As Shankar et al. (2014) have emphasized, these differences matter greatly for matching galaxies to dark matter haloes, and currently limit what we can learn about galaxy formation from this match.

While we do not claim that one particular stellar population model is better than another – Table 1 provides our results in tabular form for a variety of such models (for the photometry, we use the SerExp magnitudes from Meert et al. 2015) – we hope that our summary of the current state of the field will spur work towards reducing these systematic differences. For example, we have focused on the systematics in determining $\phi(M_*)$. However, the clustering of these objects also depends on $M_*$, especially at the high masses of interest here. Therefore, it may be that halo model interpretations
<table><caption>**Table 1.** Table of observed (i.e. error-broadened) and intrinsic stellar mass functions $\Phi \equiv \log_{10}[\ln(10)M_*\phi(M_*)]$, in units of Mpc$^{-3}$ dex$^{-1}$, as a function of $M_* \equiv \log_{10}(M_*/M_\odot)$, for a range of stellar population models with SerExp photometry from Meert et al. (2015). M14$_{d-f}$ and M14$_d$ are the values associated with the dust-free and dusty estimate of $M_*$ from Mendel et al. (2014), respectively, and MPA and WM11 are the MPA-JHU and the WiscM11 (Chen et al. 2012) values associated with the $M_*$ estimates from tables galspecExtra and stellarMassPCAWiscM11 available at www.sdss.org/dr12/spectro/galaxy.</caption><thead><tr><th>M<sub>*</sub></th><th>&Phi;<sup>M14<sub>d-f</sub></sup><sub>Obs</sub></th><th>&Phi;<sup>M14<sub>d-f</sub></sup><sub>Int</sub></th><th>E<sup>M14<sub>d-f</sub></sup></th><th>&Phi;<sup>M14<sub>d</sub></sup><sub>Obs</sub></th><th>&Phi;<sup>M14<sub>d</sub></sup><sub>Int</sub></th><th>E<sup>M14<sub>d</sub></sup></th><th>&Phi;<sup>MPA</sup><sub>Obs</sub></th><th>&Phi;<sup>MPA</sup><sub>Int</sub></th><th>E<sup>MPA</sup></th><th>&Phi;<sup>WM11</sup><sub>Obs</sub></th><th>&Phi;<sup>WM11</sup><sub>Int</sub></th><th>E<sup>WM11</sup></th></tr></thead><tbody><tr><td>9.05</td><td>-2.051</td><td>-2.069</td><td>0.022</td><td>-2.066</td><td>-2.088</td><td>0.019</td><td>-2.097</td><td>-2.103</td><td>0.018</td><td>-2.088</td><td>-2.104</td><td>0.021</td></tr><tr><td>9.15</td><td>-2.075</td><td>-2.080</td><td>0.020</td><td>-2.101</td><td>-2.095</td><td>0.017</td><td>-2.117</td><td>-2.124</td><td>0.017</td><td>-2.101</td><td>-2.112</td><td>0.019</td></tr><tr><td>9.25</td><td>-2.092</td><td>-2.091</td><td>0.018</td><td>-2.095</td><td>-2.108</td><td>0.016</td><td>-2.138</td><td>-2.148</td><td>0.015</td><td>-2.136</td><td>-2.122</td><td>0.017</td></tr><tr><td>9.35</td><td>-2.089</td><td>-2.106</td><td>0.016</td><td>-2.124</td><td>-2.138</td><td>0.014</td><td>-2.152</td><td>-2.170</td><td>0.014</td><td>-2.118</td><td>-2.136</td><td>0.016</td></tr><tr><td>9.45</td><td>-2.118</td><td>-2.139</td><td>0.015</td><td>-2.154</td><td>-2.163</td><td>0.013</td><td>-2.185</td><td>-2.207</td><td>0.013</td><td>-2.141</td><td>-2.164</td><td>0.014</td></tr><tr><td>9.55</td><td>-2.159</td><td>-2.167</td><td>0.014</td><td>-2.175</td><td>-2.188</td><td>0.012</td><td>-2.236</td><td>-2.237</td><td>0.012</td><td>-2.187</td><td>-2.206</td><td>0.013</td></tr><tr><td>9.65</td><td>-2.183</td><td>-2.196</td><td>0.012</td><td>-2.198</td><td>-2.221</td><td>0.011</td><td>-2.246</td><td>-2.252</td><td>0.011</td><td>-2.229</td><td>-2.236</td><td>0.012</td></tr><tr><td>9.75</td><td>-2.208</td><td>-2.212</td><td>0.011</td><td>-2.242</td><td>-2.253</td><td>0.010</td><td>-2.253</td><td>-2.261</td><td>0.010</td><td>-2.241</td><td>-2.252</td><td>0.011</td></tr><tr><td>9.85</td><td>-2.213</td><td>-2.223</td><td>0.010</td><td>-2.263</td><td>-2.265</td><td>0.010</td><td>-2.267</td><td>-2.273</td><td>0.009</td><td>-2.262</td><td>-2.257</td><td>0.010</td></tr><tr><td>9.95</td><td>-2.231</td><td>-2.239</td><td>0.010</td><td>-2.269</td><td>-2.279</td><td>0.009</td><td>-2.281</td><td>-2.280</td><td>0.008</td><td>-2.252</td><td>-2.274</td><td>0.009</td></tr><tr><td>10.05</td><td>-2.252</td><td>-2.246</td><td>0.009</td><td>-2.290</td><td>-2.286</td><td>0.008</td><td>-2.279</td><td>-2.283</td><td>0.008</td><td>-2.296</td><td>-2.285</td><td>0.009</td></tr><tr><td>10.15</td><td>-2.240</td><td>-2.244</td><td>0.008</td><td>-2.285</td><td>-2.285</td><td>0.008</td><td>-2.286</td><td>-2.292</td><td>0.007</td><td>-2.282</td><td>-2.287</td><td>0.008</td></tr><tr><td>10.25</td><td>-2.252</td><td>-2.252</td><td>0.007</td><td>-2.286</td><td>-2.291</td><td>0.007</td><td>-2.293</td><td>-2.300</td><td>0.006</td><td>-2.287</td><td>-2.293</td><td>0.007</td></tr><tr><td>10.35</td><td>-2.250</td><td>-2.265</td><td>0.006</td><td>-2.292</td><td>-2.298</td><td>0.006</td><td>-2.307</td><td>-2.323</td><td>0.006</td><td>-2.298</td><td>-2.294</td><td>0.007</td></tr><tr><td>10.45</td><td>-2.274</td><td>-2.295</td><td>0.006</td><td>-2.304</td><td>-2.321</td><td>0.006</td><td>-2.338</td><td>-2.361</td><td>0.006</td><td>-2.291</td><td>-2.308</td><td>0.006</td></tr><tr><td>10.55</td><td>-2.314</td><td>-2.337</td><td>0.005</td><td>-2.335</td><td>-2.365</td><td>0.005</td><td>-2.382</td><td>-2.409</td><td>0.005</td><td>-2.318</td><td>-2.344</td><td>0.005</td></tr><tr><td>10.65</td><td>-2.357</td><td>-2.392</td><td>0.005</td><td>-2.391</td><td>-2.421</td><td>0.005</td><td>-2.434</td><td>-2.476</td><td>0.005</td><td>-2.370</td><td>-2.390</td><td>0.005</td></tr><tr><td>10.75</td><td>-2.421</td><td>-2.465</td><td>0.005</td><td>-2.449</td><td>-2.497</td><td>0.005</td><td>-2.511</td><td>-2.563</td><td>0.005</td><td>-2.411</td><td>-2.451</td><td>0.005</td></tr><tr><td>10.85</td><td>-2.504</td><td>-2.561</td><td>0.005</td><td>-2.540</td><td>-2.592</td><td>0.005</td><td>-2.608</td><td>-2.664</td><td>0.005</td><td>-2.486</td><td>-2.533</td><td>0.005</td></tr><tr><td>10.95</td><td>-2.611</td><td>-2.675</td><td>0.005</td><td>-2.639</td><td>-2.702</td><td>0.005</td><td>-2.717</td><td>-2.790</td><td>0.005</td><td>-2.576</td><td>-2.632</td><td>0.005</td></tr><tr><td>11.05</td><td>-2.732</td><td>-2.812</td><td>0.005</td><td>-2.758</td><td>-2.844</td><td>0.005</td><td>-2.857</td><td>-2.945</td><td>0.005</td><td>-2.680</td><td>-2.760</td><td>0.005</td></tr><tr><td>11.15</td><td>-2.885</td><td>-2.975</td><td>0.005</td><td>-2.919</td><td>-3.011</td><td>0.005</td><td>-3.022</td><td>-3.127</td><td>0.006</td><td>-2.831</td><td>-2.918</td><td>0.005</td></tr><tr><td>11.25</td><td>-3.055</td><td>-3.159</td><td>0.006</td><td>-3.097</td><td>-3.207</td><td>0.006</td><td>-3.220</td><td>-3.329</td><td>0.006</td><td>-2.996</td><td>-3.100</td><td>0.005</td></tr><tr><td>11.35</td><td>-3.252</td><td>-3.369</td><td>0.007</td><td>-3.304</td><td>-3.427</td><td>0.007</td><td>-3.424</td><td>-3.551</td><td>0.007</td><td>-3.191</td><td>-3.303</td><td>0.006</td></tr><tr><td>11.45</td><td>-3.472</td><td>-3.601</td><td>0.008</td><td>-3.536</td><td>-3.673</td><td>0.008</td><td>-3.660</td><td>-3.795</td><td>0.009</td><td>-3.406</td><td>-3.531</td><td>0.007</td></tr><tr><td>11.55</td><td>-3.715</td><td>-3.841</td><td>0.009</td><td>-3.793</td><td>-3.942</td><td>0.009</td><td>-3.913</td><td>-4.067</td><td>0.011</td><td>-3.642</td><td>-3.772</td><td>0.008</td></tr><tr><td>11.65</td><td>-3.950</td><td>-4.114</td><td>0.011</td><td>-4.066</td><td>-4.249</td><td>0.012</td><td>-4.198</td><td>-4.360</td><td>0.013</td><td>-3.892</td><td>-4.043</td><td>0.010</td></tr><tr><td>11.75</td><td>-4.254</td><td>-4.410</td><td>0.014</td><td>-4.397</td><td>-4.616</td><td>0.015</td><td>-4.507</td><td>-4.681</td><td>0.017</td><td>-4.176</td><td>-4.344</td><td>0.013</td></tr><tr><td>11.85</td><td>-4.554</td><td>-4.732</td><td>0.018</td><td>-4.799</td><td>-5.010</td><td>0.022</td><td>-4.828</td><td>-5.079</td><td>0.023</td><td>-4.483</td><td>-4.673</td><td>0.017</td></tr><tr><td>11.95</td><td>-4.881</td><td>-5.121</td><td>0.024</td><td>-5.202</td><td>-5.457</td><td>0.033</td><td>-5.276</td><td>-5.500</td><td>0.036</td><td>-4.840</td><td>-4.997</td><td>0.023</td></tr><tr><td>12.05</td><td>-5.321</td><td>-5.548</td><td>0.037</td><td>-5.673</td><td>-5.984</td><td>0.052</td><td>-5.696</td><td>-6.073</td><td>0.054</td><td>-5.146</td><td>-5.377</td><td>0.034</td></tr><tr><td>12.15</td><td>-5.742</td><td>-6.039</td><td>0.055</td><td>-6.233</td><td>-6.558</td><td>0.092</td><td>-6.297</td><td>-6.785</td><td>0.102</td><td>-5.539</td><td>-5.900</td><td>0.056</td></tr></tbody></table>

of the dependence of clustering on $M_*$ can be used to help decide between different stellar population assumptions. Until such systematics have been reduced, using the same models and photometry whenever low- and high-redshift samples are compared – as was done by Bernardi et al. (2016) when comparing the CMASS galaxies at $z \sim 0.6$ to this SDSS sample at $z \sim 0.1$ – is essential. The statistical power of large data sets is now sufficiently large that being careless about this will lead to systematically biased conclusions.

## ACKNOWLEDGEMENTS

We are grateful to the staff of the LYTE centre for its hospitality when this work was completed.

## REFERENCES

* Abazajian et al., 2009, ApJS, 182, 543

* Ahn et al., 2012, ApJS, 203, 21

* Baldry I. K. et al., 2012, MNRAS, 421, 621

* Bell E. F., McIntosh D. H., Katz N., Weinberg M. D., 2003, ApJS, 149, 289

* Bernardi M., Hyde J. B., Sheth R. K., Miller C. J., Nichol R. C., 2007, AJ, 133, 1741

* Bernardi M., Shankar F., Hyde J. B., Mei S., Marulli F., Sheth R. K., 2010, MNRAS, 404, 2087

* Bernardi M., Roche N., Shankar F., Sheth R. K., 2011, MNRAS, 412, L6

* Bernardi M., Meert A., Sheth R. K., Vikram, Huertas-Company M., Mei S., Shankar F., 2013, MNRAS, 436, 697 (B13)

* Bernardi M., Meert A., Vikram V., Huertas-Company M., Mei S., Shankar F., Sheth R. K., 2014, MNRAS, 443, 874

* Bernardi M., Meert A., Sheth R. K., Huertas-Company M., Maraston C., Shankar F., Vikram V., 2016, MNRAS, 455, 4122

* Blanton M. R., Roweis S., 2007, AJ, 133, 734

* Blanton M. R., Eisenstein D., Hogg D. W., Schlegel D. J., Brinkmann J., 2005, ApJ, 629, 143

* Blanton M. R., Kazin E., Muna D., Weaver B. A., Price-Whelan A., 2011, AJ, 142, 31

* Brinchmann J., Charlot S., White S. D. M., Tremonti C., Kauffmann G., Heckman T., 2004, MNRAS, 351, 1151

* Bruzual G., Charlot S., 2003, MNRAS, 344

* Chabrier G., 2003, ApJ, 586, 133

* Chen Y. et al., 2012, MNRAS, 421, 314

* Conroy C., Gunn J. E., White M., 2009, ApJ, 699

* Cooray A., Sheth R. K., 2002, Phys. Rep., 372, 1

* D'Souza R., Vegetti S., Kauffmann G. A. M., 2015, MNRAS, 454, 4027

* Fioc M., Rocca-Volmerange B., 1997, A&A, 326, 950

* Fioc M., Rocca-Volmerange B., 1999, preprint (astro-ph/9912179)

* Fischer J.-L., Bernardi M., Meert A., 2017, MNRAS, in press

* Gallazzi A., Bell E. F., 2009, ApJS, 185, 253

* Hyde J. B., Bernardi M., 2009, MNRAS, 394, 1978

* Kauffmann G. et al., 2003, MNRAS, 341, 33

* Kelvin L. S. et al., 2012, MNRAS, 421, 1007

* Kelvin L. S. et al., 2014, MNRAS, 444, 1647
Kroupa P., 2001, MNRAS, 322, 231

Li C., White S. D. M., 2009, MNRAS, 398, 4, 2177

Maraston C., Strömbäck G., 2011, MNRAS, 418, 2785

Maraston C., Strömbäck G., Thomas D., Wake D. A., Nichol R. C., 2009, MNRAS, 394, L107

Maraston C. et al., 2013, MNRAS, 435, 2764

Meert A., Vikram V., Bernardi M., 2013, MNRAS, 433, 1344

Meert A., Vikram V., Bernardi M., 2015, MNRAS, 446, 3943

Meert A., Vikram V., Bernardi M., 2016, MNRAS, 455, 2440

Mendel J. T., Simard L., Palmer M., Ellison S. L., Patton D. R., 2014, ApJS, 210

Moustakas J. et al., 2013, ApJ, 767, 50

Schmidt M., 1968, ApJ, 151, 393

Sérsic J. L., 1963, Bol. Asociación Argentina Astron. La Plata Argentina, 6, 41

Shankar F. et al., 2014, ApJ, 797, L27

Sheth R. K., 2017, MNRAS, in press

Simard L., Mendel J. T., Patton D. R., Ellison S. L., McConnachie A. W., 2011, ApJS, 196, 11 (Simard11)

Taylor E. N. et al., 2015, MNRAS, 446, 2144

Thanjavur K., Simard L., Bluck A. F. L., Mendel T., 2016, MNRAS, 459, 44

Thomas D. et al., 2013, MNRAS, 431, 1383

Tremonti C. A. et al., 2004, ApJ, 613, 898

Weigel A. K., Schawinski K., Bruderer C., 2016, MNRAS, 459, 2150

## APPENDIX A: DEPENDENCE OF $M_*/L$ VALUES ON STELLAR POPULATION MODELLING AND DUST

The main text made the point that it is important to separate the effects of $\ell$ from those of $M_*/L$ on $\phi(M_*)$. This appendix compares the $M_*/L$ values of the different groups, and shows how the inclusion of dust impacts the estimates. In all cases, the photometry is SerExp of B13-Meert15, and the IMF is Chabrier.

We begin with a comparison of the dusty and dust-free $M_*/L$ values from Mendel et al. (2014) in Fig. A1. The solid cyan line shows the median difference for a number of narrow bins in $M_*$. The dashed lines show the range around this median that contains 68 per cent of the objects in each $M_*$ bin. The (large!) differences at low $M_*$ are easy to understand: if dust is not accounted for, red spirals will be assigned larger masses. However, these are the objects of

<figure><img src="image_18.png" /><figcaption>**Figure A1.** Comparison of the dust-free and dusty $M_*/L$ values from Mendel et al. (2014). At the highest masses, the dust-free models return higher masses. Although we only show $10^4$ objects, the solid line shows the median relation for the full sample, and dotted lines show the range within which 68 per cent of the full sample lies.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure A2.** Comparison of the B13 and Mendel et al. (2014) dust-free $M_*/L$ values for a Chabrier IMF. Symbols and line styles same as in Fig. A1.</figcaption></figure>

least relevance to our work. At intermediate masses, there appear to be two populations, one of which is assigned smaller masses when moving from dusty to dust-free models. This is mostly irrelevant for the stellar mass function, since it is relatively flat at masses below about $10^{11} M_\odot$, so shifts in mass make little difference. At the higher masses of most interest to us, these shifts do matter, but in this regime the dust-free models again return higher masses.

The offset at large masses is part of the reason why Thanjavur et al. (2016) reported that their Mendel14-based $\phi(M_*)$ determination was smaller than that of B13. Namely, their determination used dusty models. Fig. A2 shows that the offset between the Mendel14 dust-free $M_*/L$ values and those of B13 is of the order of 0.05 dex at high masses. Combining this with Fig. A1 shows that the difference in $M_*/L$ between B13 and the Mendel14 dusty models is of the order of 0.1 dex. As Fig. 11 in the main text shows, accounting for this offset (i.e. scaling the Mendel14 and B13 $M_*$ estimates to the same $M_*/L$ values) reduces the difference in $\phi(M_*)$ to less than 0.1 dex. This indicates that the ‘fiducial’ photometry of Thanjavur et al. (2016) that is based on Simard11 is very similar to the SerExp values of Meert15 that were used by B13, and that we use as the fiducial photometry in the main text. In other words, the photometry agrees to better than 0.1 dex.

Having shown the systematics associated with only changing the dust model, we now compare $M_*/L$ estimates when the treatment of the stellar population or the method of fitting the model to the data is varied. Because we already know dusty and dust-free models differ, we do not compare the dusty models of one group with dust-free models of another.

We begin with a comparison of the MPA-JHU and Mendel14 models. Fig. 6 in Mendel et al. (2014) suggests that their $M_*/L$ values are in good agreement with the MPA-JHU values. However, they did not account for the difference between their Chabrier and the MPA-JHU Kroupa IMFs. Our Fig. A3 accounts for this difference (as one should), and this shows that the Mendel et al. (2014) values are actually about 0.05 dex larger than the MPA-JHU values. Some of the difference could be related to the inclusion of bursty models (MPA-JHU) or not (Mendel14). Above $10^{11} M_\odot$, the scatter around the median is ~0.05 dex.

Fig. A4 shows a similar comparison between the Wisconsin and Mendel14 values. Comparison with the previous figure shows that the Wisconsin $M_*/L$ values are about 0.1 dex larger than the MPA-JHU values. This is larger than the 0.05 dex scatter between the two
<figure><img src="image_20.png" /><figcaption>**Figure A3.** Comparison of the MPA-JHU and Mendel et al. (2014) $M_*/L$ values for a Chabrier IMF. [Both estimates assume that galaxies are dusty. However, the MPA-JHU includes bursty models while Mendel et al. (2014) do not.] Symbols and line styles same as in Fig. A1.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure A4.** Comparison of the Wisconsin and Mendel et al. (2014) $M_*/L$ values for a Chabrier IMF. (Both estimates assume that galaxies are dusty.) Comparison with the previous figure suggests that the Wisconsin $M_*/L$ values are about 0.1 dex larger than the MPA-JHU values. Symbols and line styles same as in Fig. A1.</figcaption></figure>

models. Finally, the top panel of Fig. A5 shows a similar comparison between the GranadaWide and Mendel14 values. In this case, the scatter is substantially larger than in the previous two figures. The reason for this becomes evident in the bottom panel, which compares the dust-free versions of the two models. The two clouds of points are approximately equally populated in the dust-free models, but not so when dust is included. It is the fact that there appear to be two populations that makes the scatter around the median so much

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /><figcaption>**Figure A5.** Comparison of the GranadaWide and Mendel et al. (2014) $M_*/L$ values for a Chabrier IMF. Symbols and line styles same as in Fig. A1. Top and bottom panels show results when galaxies are assumed to be dusty and dust-free, respectively. In contrast to the previous two figures, there appear to be two populations of galaxies that are approximately equally populated in the dust-free models, but not so when dust is included. As a result, the scatter around the median is much larger here than in the previous two figures.</figcaption></figure>

larger. The appearance of two populations is surprising, since both GranadaWide and Mendel14 are based on the same FSPS library, so the differences can only arise from details in how the library was used, and how the fitting was done.

Finally, for completeness, we compare the Portsmouth-Passive and Mendel14 values (Fig. A6). The disagreement at low masses should be ignored, because the Portsmouth-Passive models are only expected to be realistic at high masses. Above $10^{11.3} M_{\odot}$, they are within 0.05 dex of the Mendel15 et al. values.
<figure><img src="image_24.png" /></figure>

**Figure A6.** Comparison of the Portsmouth-Passive and Mendel et al. (2014) $M_*/L$ values for a Chabrier IMF. (Both estimates assume that galaxies are dust-free.) The disagreement at low masses should be ignored, because the Portsmouth-Passive models are only expected to be realistic at high masses $M_* > 10^{11} M_{\odot}$. Symbols and line styles same as in Fig. A1.

This paper has been typeset from a TeX/LATEX file prepared by the author.