The Massive End of the Stellar Mass Function
MMon. Not. R. Astron. Soc. , 1–11 (2002) Printed 10 October 2018 (MN L A TEX style file v2.2)
The Massive End of the Stellar Mass Function
Richard D’Souza (cid:63) , Simona Vegetti , Guinevere Kauffmann Max Plank Institute for Astrophysics, Munich, Germany
Accepted 1988 December 15. Received 1988 December 14; in original form 1988 October 11
ABSTRACT
We derive average flux corrections to the
Model magnitudes of the Sloan Digital SkySurvey (SDSS) galaxies by stacking together mosaics of similar galaxies in bins of stel-lar mass and concentration. Extra flux is detected in the outer low surface brightnesspart of the galaxies, leading to corrections ranging from 0.05 to 0.32 mag for the high-est stellar mass galaxies. We apply these corrections to the MPA-JHU (Max-PlanckInstitute for Astrophysics - John Hopkins University) stellar masses for a completesample of half a million galaxies from the SDSS survey to derive a corrected galaxystellar mass function at z = 0 . in the stellar mass range . < log( M ∗ /M (cid:12) ) < . .We find that the flux corrections and the use of the MPA-JHU stellar masses have asignificant impact on the massive end of the stellar mass function, making the slopesignificantly shallower than that estimated by Li & White (2009), but steeper thanderived by Bernardi et al. (2013). This corresponds to a mean comoving stellar massdensity of galaxies with stellar masses log( M ∗ /M (cid:12) ) (cid:62) . that is a factor of 3.36larger than the estimate by Li & White (2009), but is 43% smaller than reported byBernardi et al. (2013). Key words:
Galaxy Formation – Stellar haloes
The stellar mass function of galaxies is a basic probe ofgalaxy formation and evolution enabled by large redshiftsurveys. In recent years, major advances have been made bylarge redshift surveys, such as the 2dF Galaxy Redshift Sur-vey and the Sloan Digital Sky Survey (SDSS), in estimatingthe stellar mass function in the low-redshift Universe (Coleet al. 2001; Bell et al. 2003; Blanton et al. 2003). For exam-ple, Li & White (2009) have used a uniform sample of almosthalf a million galaxies from SDSS DR7 to derive the stellarmass function at z = 0 . . This has been complemented bythe effort of the Galaxy and Mass Assembly Survey (GAMABaldry et al. 2012), which has accurately constrained thefaint end slope of the stellar mass function down to stellarmasses ∼ M (cid:12) .The calculation of the stellar mass function hinges onthe proper determination of the stellar mass of a galaxy,which in turn depends critically on the estimation of its to-tal flux in a given pass-band. Systematic differences in theestimation of the stellar mass of a galaxy may arise fromdifferent choices of the initial mass function (IMF) and thestellar mass-to-light ratio (M/L), as well as from differentestimations of the galaxy total flux. Determining the fluxaccurately for a large number of galaxies in an all-sky survey (cid:63) E-mail address: [email protected] (RDS) is a challenging task. In particular, quantifying the flux inthe outer low surface brightness (LSB) regions of a galaxyhas proven to be difficult and is still subject of much de-bate (Bernardi et al. 2013; Simard et al. 2011). These un-certainties mean that the slope at the massive end of themass function is not very well determined. This has signifi-cant implications for several astrophysical problems, includ-ing halo occupation models, the mean baryon fraction inthe Universe, X-ray and Sunavey-Zeldovich studies of highmass galaxies, and understanding the evolution of massivegalaxies to high redshifts.Different approaches have been employed by SDSS inits photometric pipeline (
PHOTO ) to estimate the total fluxof a galaxy. In addition to SDSS
Petrosian magnitudes,two dimensional models (e.g. exponential or de Vaucouleurs)have been used to model the surface brightness distributionof galaxies (SDSS
Model magnitudes). Further improvementhas been provided by SDSS cModel magnitudes, for whichfluxes are estimated as a linear combination of an exponen-tial and a de Vaucouleurs model. In recent years, severalstudies have tried to fit Sérsic and multi-component modelsto the surface brightness distribution (Simard et al. 2011;Lackner & Gunn 2012; Bernardi et al. 2013).Each of these approaches provides a progressively betterestimate of the total flux of a galaxy, but they all suffer fromthe same intrinsic drawback, namely that the models arefits to the central, high signal-to-noise ratio (SNR) regions c (cid:13) a r X i v : . [ a s t r o - ph . GA ] S e p R. D’Souza, S.Vegetti, G. Kauffmann of the galaxy and assumptions are required about the outerlower SNR (beyond µ r ∼
27 mag arcsec − ) part of the galaxyprofile. Additionally, the total flux estimated through modelfitting can be biased in a number of ways.The biggest source of systematic bias in the flux deter-mination is related to the estimation of the sky background,especially for large nearby objects or those located in denseenvironments (von der Linden 2007; Bernardi et al. 2007).In principle, this can be overcome by considering extremelylarge fields of view. For example, considerable progress hasbeen achieved by Blanton et al. (2011) by fitting the maskedsky background for each SDSS scan with a smooth continu-ous function.However, even with improvements to the sky back-ground algorithm, one is still limited by the depth of thesurvey. The relatively short exposure time of SDSS (53.9secs) limits the accuracy of the background determinationand subtraction. This in turn limits ones ability to distin-guish between the flux of the outer stellar halo and the skybackground, leading to an over- or under-estimation of thetotal flux of a galaxy. In particular, multi-component modelfitting of the main galaxy can lead to biased results. Thismay explain why recent attempts to trace the low SNR LSBpart of a galaxy through fitting multi-component models tosingle SDSS photometric images have yielded divergent re-sults (Simard et al. 2011; Bernardi et al. 2013; Meert et al.2015).Other sources of systematic error in determining theflux of a particular object are the procedures employed fordeblending and masking, as well as the radial extent of themodels used for the surface brightness fitting. Finally, in ad-dition to photometry, several other effects have a consider-able impact on the massive end of the stellar mass function,such as evolutionary corrections and fiber collisions (i.e. thefraction of galaxies not targeted for spectroscopy due to thefact that fibres cannot be positioned closer together than 55arcseconds on the SDSS plug plates).An alternative but viable approach to fitting models toindividual images of galaxies, is to stack images of similargalaxies to quantify the average total amount of extra lightin the outer parts (Tal & van Dokkum 2011; D’Souza etal. 2014). By stacking galaxies as a function of their stellarmass and galaxy-type, D’Souza et al. (2014) have reacheda depth of µ r ∼
32 mag arcsec − . The increased depth ofgalaxy stacks helps to reliably constrain the total amountof light especially in the LSB component. In addition, thebackground for stacked galaxies can be determined moreaccurately. This then provides a direct handle on the cor-rections to the Model magnitudes as a function of the stellarmass and galaxy type.In this paper, we attempt to derive flux corrections tothe
Model magnitudes and re-derive the galaxy stellar massfunction at redshift z = 0 . using MPA-JHU (Max-PlanckInstitute for Astrophysics & John Hopkins University) stel-lar masses (Kauffmann et al. 2003) and the sample of Li &White (2009). We estimate corrections to the Model magni-tudes by stacking volume-limited samples in bins of stellarmass, concentration and model type. We also consider var-ious effects that may systematically bias the stellar massfunction.In Section 2, we define the samples used for derivingthe corrections as well as the full sample used to derive the stellar mass function. In Section 3, we derive the flux cor-rections to the
Model magnitudes. In Sections 4 and 5, wederive the galaxy stellar mass function and the luminosityfunction respectively. In Sections 6 and 7, we summarise anddiscuss our results. Throughout this paper, we assume a flat Λ CDM cosmology, Ω m = 0 . and Ω Λ = 0 . . We furtherassume a Hubble parameter h = 0 . for the calculation ofphysical distance scales wherever necessary. Following Li & White (2009), we select SDSS spectroscopicgalaxies from the NYU-VAGC (New York University - ValueAdded Catalogue) catalogue (Blanton et al. 2005) withredshifts in the range . (cid:54) z (cid:54) . and Petrosian r -band magnitudes in the range (cid:54) m rPet (cid:54) . . Thisgives us a total of 533442 galaxies, which are ideal for largescale structure studies. We further pruned the sample to523476 galaxies by retaining only those galaxies with a validMPA-JHU stellar mass. We estimate the “effective” surveyarea to be deg (2.0084 steradians), by taking into ac-count the incompletness and the masked-out regions (due tobright stars) of the survey.For the stellar mass function, we use the stellar massesprovided in the DR7 version of the MPA-JHU catalogue ,which assumes a universal Chabrier initial stellar mass func-tion ( Chabrier 2003).To derive the luminosity function, we use the r-bandabsolute Model magnitude ( M r . ), corrected for evolutionand K-corrected to its value at z = 0 . according to thefollowing equation: M = m − DM ( z ) − K ( z ; z ) + Q e ( z − z ) , (1)where M is the absolute magnitude, DM ( z ) is the dis-tance modulus at redshift z , m the apparent magnitude, K ( z ; z ) is the K -corrections relative to a passband blue-shifted by z and the luminosity e -correction is parametrisedlinearly by Q e . The K -corrections were calculated using thecode Kcorrect v4_3 (Blanton & Roweis 2007). In general,we assume a uniform luminosity evolutionary correction of Q r = 1 . as derived by (Blanton et al. 2003). To derive the corrections to the
Model magnitudes, we stackvolume-limited sub-samples of isolated galaxies defined fromthe parent sample in various ranges of stellar mass, concen-tration ( R /R ) and redshift (See Table 1). In each sub-sample, galaxies that were better fit by an exponential ( Exp )or a de Vaucouleurs ( deV ) model by the SDSS pipeline (de-fined by comparing the likelihood values of the model fitsfrom the SDSS
PhotoObjAll database) were stacked sepa-rately.We select isolated galaxies by requiring that there areno brighter companions in the spectroscopic sample within Available at http://sdss.physics.nyu.edu/vagc/ . We also include the three survey strips in the Southern Cap. (cid:13) , 1–11 tellar Mass Function log M ∗ / ( h − M fl ) z
33% 37.6%42.8%49.5%60.6% 76.9% 93.4%
Figure 1.
Volume-limited samples used for stacking: Shaded con-tours show the distribution of the parent sample galaxies in theplane of stellar mass versus redshift. The seven coloured boxesindicate the redshift limits of the seven stellar mass sub-samples.These sub-samples are further divided by concentration (C) (notshown in the figure). The numbers in the coloured boxes indi-cate the fraction of centrals in these volume limited stellar masssub-samples. R (cid:54) (where R is the projected comoving separation)and | δz | < − .In Figure 1, the redshift limits of the various stellar masssub-samples are shown projected along the plane of stellarmass versus redshift of the parent NYU-VAGC sample. Thefraction of centrals in each stellar mass sub-sample are alsoshown. MODEL
MAGNITUDES
In this section, we derive corrections to the original SDSS
Model magnitudes derived from the standard DR7 photo-metric pipeline ( photo v5_4 ). We first demonstrate that theoriginal SDSS
Model magnitudes are biased, in agreementwith other studies (Simard et al. 2011; Bernardi et al. 2013).We then proceed to derive corrections to the
Model magni-tudes.
Model magnitudes
The original
Model are affected by two sources of systematicbias related to over-subtraction of the sky background (vonder Linden 2007; Bernardi et al. 2007) and to simplisticchoices of the model for the surface brightness profile of thegalaxy (exponential or de Vaucouleurs). In this section, weallow for more complex models to describe the light profilesand we also allow the sky background to vary. We fit 2D axisymmetric models (single Sérsic and dou-ble Sérsic models along with a constant background) usingthe Bayesian analysis described by D’Souza et al. (2014) toindividual postage-stamp cutouts of the highest stellar massand high-concentration galaxies ( . < log( M ∗ /M (cid:12) ) < . , . < C < . ) in the redshift range . < z < . (covered by the sample G4 above - 38 galaxies) and . Model mag-nitudes reported by the SDSS photo v5_4 pipeline. In Figure2, we plot a histogram of M model − M fit for each galaxy. Thedistribution is broad with a standard deviation of 0.25 mag-nitudes and is positively skewed. The median is shifted byis 0.03 magnitudes and the mean by 0.08 magnitudes. Thelarge spread in the histogram arises from a degeneracy be-tween the best-fit model and the level of sky background.The results in Figure 2 demonstrate that shallow single-exposure SDSS images are insufficient to accurately quantifythe total amount of light in massive early-type galaxies tobetter than 0.25 mag. We also note that estimates of thetotal flux from a single SDSS image will also be affected bythe deblending and masking algorithm .Because of the limitations in estimating total fluxesfrom single SDSS images, we have chosen to correct Model magnitudes using stacked images , where the increasedsignal-to-noise ratio better constrains both the model andthe level of sky background. In order to derive the flux corrections to the Model mag-nitudes, we used the sky-subtracted SDSS Data Release 9images to create mosaics in the g , r and i bands centredon each galaxy in the sub-samples defined in Section 2.2.The mosaics extend out to radii of 0.6 - 1 Mpc dependingon the stellar mass and redshift range. We follow the stack-ing procedure outlined in D’Souza et al. (in preparation)and similar to that used by D’Souza et al. (2014). In short,each mosaic was deblended, masked, corrected for galacticextinction (Schlegel et al. 1998), transformed to the highestredshift in that respective bin, rotated so that the majoraxis of each galaxy is aligned, and then stacked using thetruncated-mean algorithm. The g - and the i - band imageswere only used to create the final mask along with the r -bandimages. Conservative masking was used. The final stackingwas done using the masked and transformed r -band images. Measuring the total integrated flux of a galaxy stack byfitting a model to its light distribution misses a fair amount In this paper, we follow the deblending and masking techniqueoutlined in D’Souza et al. (in preparation). For the truncated-mean stack, we removed 5% of the extrememinimum and maximum values for each pixel.c (cid:13)000 Model mag-nitudes, we used the sky-subtracted SDSS Data Release 9images to create mosaics in the g , r and i bands centredon each galaxy in the sub-samples defined in Section 2.2.The mosaics extend out to radii of 0.6 - 1 Mpc dependingon the stellar mass and redshift range. We follow the stack-ing procedure outlined in D’Souza et al. (in preparation)and similar to that used by D’Souza et al. (2014). In short,each mosaic was deblended, masked, corrected for galacticextinction (Schlegel et al. 1998), transformed to the highestredshift in that respective bin, rotated so that the majoraxis of each galaxy is aligned, and then stacked using thetruncated-mean algorithm. The g - and the i - band imageswere only used to create the final mask along with the r -bandimages. Conservative masking was used. The final stackingwas done using the masked and transformed r -band images. Measuring the total integrated flux of a galaxy stack byfitting a model to its light distribution misses a fair amount In this paper, we follow the deblending and masking techniqueoutlined in D’Souza et al. (in preparation). For the truncated-mean stack, we removed 5% of the extrememinimum and maximum values for each pixel.c (cid:13)000 , 1–11 R. D’Souza, S.Vegetti, G. Kauffmann Table 1. Volume-limited samples of isolated galaxies selected by stellar mass from the NYU-VAGC sample for the purpose of stacking.Sample Stellar mass Concentration Redshift N gal Exp N gal deV A1 . < log( M ∗ /M (cid:12) ) < . 89 1 . < C < . . < z < . 797 117A2 . < log( M ∗ /M (cid:12) ) < . 89 2 . < C < . . < z < . 66 501B1 . < log( M ∗ /M (cid:12) ) < . 09 1 . < C < . . < z < . . < log( M ∗ /M (cid:12) ) < . 09 1 . < C < . . < z < . 83 1111C1 . < log( M ∗ /M (cid:12) ) < . 29 1 . < C < . . < z < . 638 15C2 . < log( M ∗ /M (cid:12) ) < . 29 2 . < C < . . < z < . 752 308C3 . < log( M ∗ /M (cid:12) ) < . 29 2 . < C < . . < z < . 121 1499C4 . < log( M ∗ /M (cid:12) ) < . 29 2 . < C < . . < z < . . < log( M ∗ /M (cid:12) ) < . 49 1 . < C < . . < z < . 342 38D2 . < log( M ∗ /M (cid:12) ) < . 49 2 . < C < . . < z < . 534 535D3 . < log( M ∗ /M (cid:12) ) < . 49 2 . < C < . . < z < . 89 1468D4 . < log( M ∗ /M (cid:12) ) < . 49 2 . < C < . . < z < . . < log( M ∗ /M (cid:12) ) < . 69 1 . < C < . . < z < . 239 74E2 . < log( M ∗ /M (cid:12) ) < . 69 2 . < C < . . < z < . 555 1093E3 . < log( M ∗ /M (cid:12) ) < . 69 2 . < C < . . < z < . 72 1981E4 . < log( M ∗ /M (cid:12) ) < . 69 2 . < C < . . < z < . - 1867F1 . < log( M ∗ /M (cid:12) ) < . 09 1 . < C < . . < z < . 199 264F2 . < log( M ∗ /M (cid:12) ) < . 09 2 . < C < . . < z < . 303 1510F3 . < log( M ∗ /M (cid:12) ) < . 09 2 . < C < . . < z < . 76 2919F4 . < log( M ∗ /M (cid:12) ) < . 09 2 . < C < . . < z < . . < log( M ∗ /M (cid:12) ) < . 69 1 . < C < . . < z < . . < log( M ∗ /M (cid:12) ) < . 69 2 . < C < . . < z < . 11 220G3 . < log( M ∗ /M (cid:12) ) < . 69 2 . < C < . . < z < . 15 792G4 . < log( M ∗ /M (cid:12) ) < . 69 2 . < C < . . < z < . of light due to the inability of the model to reproduce thebulge/disk component of the inner part of the galaxy. Forexample, a double Sérsic model can miss upto 0.22 mag nearthe centre of a galaxy stack. On the other hand, “isophotal”magnitudes are unable to measure the LSB features of thegalaxy stack, especially in the low S/N regime.In order to measure the total integrated light in eachstack, we consider, therefore, a hybrid approach betweena “model” and an “isophotal” magnitude. In particular, wefirst fit the large mosaics using two-dimensional axisymmet-ric double Sérsic models with a flat background using theBayesian analysis described by D’Souza et al. (2014). Duringthis fit, the inner component is truncated outside a radiusequal to R e , while the outer component extends out to in-finity. Then, to the flux derived by the double Sérsic model,we add the total residual flux ( = Data − Model ) within a cir-cular aperture of limiting radius R lim , defined as the radiusat which the residual flux is maximised. The advantage ofthis hybird approach is that the double Sérsic model mea-sures the slow decline of flux into the low S/N regime, whilethe sum of the residual and model fluxes reproduces thetotal flux in the high S/N part of a galaxy stack.In the next subsections, we explore the different factorsthat may bias our measurements of the total flux of thestack. The sky residuals in the individual SDSS DR9 images are re-sponsible for some small amount of residual sky backgroundin the final stacks ( < × − nanomaggies per pixel). Wequantify this residual by adding a flat background compo-nent as a free parameter of our models (Section 3.3). As showlater in Section 3.3.2, this bias is minimal ( (cid:54) . mag).For the individual DR9 images, Blanton et al. (2011)quantified the spread in the sky background residuals to be σ ∼ . × − nanomaggies per pixel. The median biasin the r -band magnitudes was estimated to be at the most0.1 mag independent of R . In addition, higher stellar massgalaxies are found predominately at higher redshifts in theSDSS spectroscopic sample, limiting the bias in the flux ofindividual images caused by faulty sky-subtraction. In order to test how the choice of model may affect the cor-rections to the total luminosity using the “hybrid” magni-tudes, we fit each of our galaxy stack using two-dimensionalexponential, de Vaucouleurs, Sérsic and double Sérsic mod-els. Each model also includes a flat sky residual. By com-paring the evidences generated from the Bayesian fitting,we find that the double Sérsic models are preferred by more c (cid:13) , 1–11 tellar Mass Function M Model − M Fit N u m b e r C o un t s Figure 2. Bias in SDSS Model magnitudes: A histogram of thedifference between the flux derived from our best fit models ofhigh stellar mass high-concentration galaxies and the Model mag-nitudes from DR7 for galaxies with . < log( M ∗ /M (cid:12) ) < . , . < C < . , . < z < . and . < z < . .The green solid line indicates the fit to the skewed-normal distri-bution. than 10- σ over the other models in all cases. The de Vau-couleur model gives the highest estimate of the total amountof light, followed by the single Sérsic, the double Sérsic andthe exponential model respectively. Calculating the magni-tudes in the “hybrid” manner as described above gives verylittle difference in the total flux derived from different mod-els. Each model also yields different estimates of the residualbackground in the stacks. However, determining the back-ground level independently and keeping it fixed during thefitting process does not alter our estimates of the extra light(at the 0.01 mag level). This is due to the fact that the re-sults of the fitting are driven primarily by the inner highSNR part of the galaxy stack.We conclude that the combination of the depth of ourstacked image and our “hybrid” magnitudes enables us toaccurately constrain the total flux in the galaxy stack. Ourouter models are not truncated, but instead extend out toinfinity. The difference between models which are truncatedat R e and models which instead extend out to infinity isat most 0.05 mag. For each stellar mass, concentration range and model fit type(exponential or de Vaucouleurs), we measure the averageextra flux correction to the Model magnitudes as the differ-ence between the total integrated light in the stack and themedian Model flux of the galaxies in the stack. The median Model magnitude was calculated by taking the median of the individual fluxes of galaxies in the stack. We find that themedian Model magnitude is on average higher than the mean Model magnitudes. We use a two-dimensional-interpolationscheme to calculate the average extra light as a continuousfunction of stellar mass and concentration for each modeltype. These are shown in Figure 3. As can be seen, there isan extra light contribution from those galaxies which werefit by an exponential model both for high concentrationsand for high stellar masses. The extra light correction fromthose galaxies fit by a de Vaucouleurs model comes predomi-nately from the massive, high concentration galaxies. On theother hand, the de Vaucouleurs model often over-estimatesthe flux of a galaxy for low concentration massive galaxies.We note that the large width of the stellar mass bins forthe highest stellar mass galaxies may influence the correctionderived in the stacking procedure. To account for this, wedivide our highest mass sample, G4 ( . < log( M ∗ /M (cid:12) ) < . , . < C < . , . < z < . ), into smaller massbins of size 0.1 dex. We find that the relative correctionsranges from 0.23 to 0.31 mag, gradually increasing from thelowest to the highest stellar mass bin (see Figure 4). Themean correction derived by stacking the entire sample G4 is0.29 mag.For galaxies outside the the mass limits defined in Sec-tion 2.2, we extrapolate assuming the same mass correctionsof the nearest defined mass bin. In particular, at the highmass end, there are 116 galaxies with stellar masses largerthan log( M ∗ /M (cid:12) ) > . , the highest stellar mass bin usedabove. For these galaxies, we assume the corrections to bethe same as found for the highest stellar mass bin ( 0.31mag).Assuming a constant M/L for each galaxy, we calculatethe extra mass for each galaxy in our main sample given itsstellar mass, concentration and model type (by comparingthe likelihoods of the Model fits from the SDSS database)as: log M ∗ + δM ∗ M ∗ = − ∆Mag / . (2) We estimate the abundance of galaxies as a function of theirstellar mass using the /V max method outlined by Li &White (2009). In combination with the depth and the largespectroscopic sample of SDSS, the /V max method providesan unbiased estimate of the stellar mass function and its nor-malisation.In Section 4.2, we demonstrate that the /V max estimator is unbiased against large scale structure at stellarmasses of log( M ∗ /M (cid:12) ) (cid:62) . , which is the regime studiedin this work. We limit ourselves to this regime since as esti-mated by Figure 4 of Baldry et al. (2008), all galaxies abovestellar masses of log( M ∗ /M (cid:12) ) (cid:62) . , will be detected irre-spective of their central surface brightness. Moreover, ourflux corrections begin from log( M ∗ /M (cid:12) ) (cid:62) . upwards.For each observed galaxy i , we define the quantity z max,i to be the maximum redshift at which the observedgalaxy would satisfy the apparent magnitude limit of oursample m r Pet (cid:54) . . Evolutionary and K-corrections are c (cid:13)000 Model magnitudes. We use a two-dimensional-interpolationscheme to calculate the average extra light as a continuousfunction of stellar mass and concentration for each modeltype. These are shown in Figure 3. As can be seen, there isan extra light contribution from those galaxies which werefit by an exponential model both for high concentrationsand for high stellar masses. The extra light correction fromthose galaxies fit by a de Vaucouleurs model comes predomi-nately from the massive, high concentration galaxies. On theother hand, the de Vaucouleurs model often over-estimatesthe flux of a galaxy for low concentration massive galaxies.We note that the large width of the stellar mass bins forthe highest stellar mass galaxies may influence the correctionderived in the stacking procedure. To account for this, wedivide our highest mass sample, G4 ( . < log( M ∗ /M (cid:12) ) < . , . < C < . , . < z < . ), into smaller massbins of size 0.1 dex. We find that the relative correctionsranges from 0.23 to 0.31 mag, gradually increasing from thelowest to the highest stellar mass bin (see Figure 4). Themean correction derived by stacking the entire sample G4 is0.29 mag.For galaxies outside the the mass limits defined in Sec-tion 2.2, we extrapolate assuming the same mass correctionsof the nearest defined mass bin. In particular, at the highmass end, there are 116 galaxies with stellar masses largerthan log( M ∗ /M (cid:12) ) > . , the highest stellar mass bin usedabove. For these galaxies, we assume the corrections to bethe same as found for the highest stellar mass bin ( 0.31mag).Assuming a constant M/L for each galaxy, we calculatethe extra mass for each galaxy in our main sample given itsstellar mass, concentration and model type (by comparingthe likelihoods of the Model fits from the SDSS database)as: log M ∗ + δM ∗ M ∗ = − ∆Mag / . (2) We estimate the abundance of galaxies as a function of theirstellar mass using the /V max method outlined by Li &White (2009). In combination with the depth and the largespectroscopic sample of SDSS, the /V max method providesan unbiased estimate of the stellar mass function and its nor-malisation.In Section 4.2, we demonstrate that the /V max estimator is unbiased against large scale structure at stellarmasses of log( M ∗ /M (cid:12) ) (cid:62) . , which is the regime studiedin this work. We limit ourselves to this regime since as esti-mated by Figure 4 of Baldry et al. (2008), all galaxies abovestellar masses of log( M ∗ /M (cid:12) ) (cid:62) . , will be detected irre-spective of their central surface brightness. Moreover, ourflux corrections begin from log( M ∗ /M (cid:12) ) (cid:62) . upwards.For each observed galaxy i , we define the quantity z max,i to be the maximum redshift at which the observedgalaxy would satisfy the apparent magnitude limit of oursample m r Pet (cid:54) . . Evolutionary and K-corrections are c (cid:13)000 , 1–11 R. D’Souza, S.Vegetti, G. Kauffmann log ( M ∗ /h − ) C o n c Exp log ( M ∗ /h − ) C o n c deV Figure 3. The flux corrections ( ∆Mag ) as a function of stellar mass and concentration using an interpolation scheme for exponential( Exp ) and de Vaucouleurs ( DeV ) fit galaxies. log ( M ∗ /h − ) E x t r a L i g h t ( ∆ M a g ) Figure 4. The flux corrections ( ∆Mag ) as a function of stellarmass for galaxies in the sample G4 (blue circles). We also indicatethe uncertainty in the corrections by showing the flux correctionsderived from the mean Model flux of the stacks (red circles). included when calculating z max,i . Hence, z max,i is the mini-mum of the upper limit of the redshift slice and the solutionof the equation: M i = m Faintr Pet − DM ( z max ) − K ( z max ) + Q e ( z max − z i ) (3) Similarly, we also define z min,i as the minimum redshiftat which the galaxy would be present in our sample. Hence, z min,i is the maximum of the lower limit of the redshift sliceand the solution to the equation: M i = m Brightr Pet − DM ( z min ) − K ( z min ) + Q e ( z min − z i ) (4)This then allows us to calculate V max,i for the galaxyin question as the total co-moving volume of the survey be-tween z min,i and z max,i . The stellar mass function can bethen estimated as: Ψ( M ∗ )∆ M ∗ = (cid:88) i ( f norm coll,i V max,i ) − (5)where f norm coll,i is the normalised fiber collision factor de-fined below, and the sum extends over all sample galaxieswith stellar masses in the range M ∗ ± . δM ∗ . The error barsare estimated by taking into consideration both Poissonianand bootstrapping errors, as well as errors due to cosmicvariance (See 4.2).We calculate the stellar mass function in the total red-shift range . (cid:54) z (cid:54) . as well as in three redshift slices: . (cid:54) z (cid:54) . , . (cid:54) z (cid:54) . and . (cid:54) z (cid:54) . . /V max Estimator In this work, we estimate the abundance of galaxies usingthe /V max method. Given the large effective surface area(nearly deg ) and the depth of spectroscopic sample,the /V max method will be invariant to large-scale structureup to a limiting stellar mass. To test this, we divide our sam-ple into three independent but contiguous parts (Sample A,Sample B and Sample C split by right ascension), and calcu-late the standard deviation in the stellar mass function as a c (cid:13) , 1–11 tellar Mass Function -4-3-2-1 l og ( Ψ / h M p c − l og M − ∗ ) Full Sample Sample A Sample BSample C8.5 9.0 9.5 10.0 10.5 11.0 log ( M ∗ /h − ) % E rr o r Figure 5. Top panel: The stellar mass function for the completesample as well as the three independent smaller samples A, B andC, using the MPA-JHU stellar masses for the full redshift range . (cid:54) z (cid:54) . . Bottom Panel: The standard deviation in thestellar mass function as a function of stellar mass. function of stellar mass. In the bottom panel of Figure 5, thestandard deviation is plotted as a function of stellar mass.The difference in the estimates of the stellar mass functiondue to the /V max method from the three independent sam-ples is less than 10% for stellar masses log ( M ∗ /M (cid:12) ) > . .The standard deviation gives us also a handle on the er-rors in our estimates of the stellar mass function due to thecosmic variance. In calculating the stellar mass function, various systematicand random effects combine to affect the final result. Wediscuss each of these effects in turn in the following subsec-tions: The first source of systematic bias comes from the estimationof the stellar mass of individual galaxies. In this work, we usethe MPA-JHU stellar masses to calculate the stellar massfunction. This involves a change of flux (from Petrosian to Model magnitudes) and M/L ratio (from NYU-VAGC toMPA-JHU) relative to Li & White (2009).We find that the use of NYU-VAGC stellar massesbased on the Model magnitudes rather than the Petrosian magnitudes introduces a shift beyond the knee of the stellarmass function towards a shallower slope at the higher massend. This shift is then further increased when we switch toMPA-JHU stellar masses based on the Model magnitudes.The slope of the massive end of the mass function is shal-lower than that obtained by shifting the mass function de-rived from the NYU-VAGC stellar masses by 0.1 dex (Seeappendix of Li & White 2009: ∆ log M ∗ = 0 . ). At a stellarmass of log M ∗ ∼ . M (cid:12) , this accounts for an increase inthe stellar mass function by a total of 1.24 dex (a 0.57 dexincrease due to the change from Petrosian to Model mag-nitudes and a 0.67 dex increase due to the change from theNYU-VAGC to MPA-JHU M/L ratios).Our assumed M/L ratio affects our estimation of thestellar mass function. The use of the MPA-JHU stellar M/Lratios makes the slope at the massive end shallower than theNYU-VAGC M/L ratios. We note that the MPA-JHU M/Lratios are derived from models that include the possibility ofcomplex star formation histories, whereas the NYU-VAGCassumes that red galaxies can be described by single stellarpopulations. Analysis of spectra of massive galaxies in theBOSS survey by Chen et al. (2013) indicates that the starformation histories of the most massive galaxies are charac-terised by episodic star formation histories.The extra flux derived from the photometry of stackedgalaxies introduces a further shift, making the slope at themassive end of the stellar mass function even shallower. Wefind that this shift of the stellar mass function is indepen-dent of whether we apply the corrections only to the centralgalaxies, or to all the galaxies in the sample. Although asmall difference is found at the knee of the mass function,both results are consistent with each other within the errorbars. The second source of systematic bias is caused by fiber col-lisions. The NYU-VAGC catalogue lists the spectroscopiccompleteness f sp of each galaxy, defined as the fraction ofphotometrically defined target galaxies in the subarea forwhich usable spectra are obtained. The NYU-VAGC cata-logue calculates the average completeness for each of thesesubareas by taking into consideration overlapping plates. Inthe jargon of the NYU-VAGC catalogue, these subareas arecalled sectors . f sp contains information about the missinggalaxies due to lack of fibers in dense regions, missing galax-ies due to spectroscopic failures, and missing galaxies due tofiber collisions. The average f sp for the sample defined aboveis 0.9146. However, f sp assumes that all galaxies with mea-sured spectra are randomly distributed within a sector, andhence cannot account for specific differences between highand low density regions in the same sector. In particular, dueto fiber collisions, certain galaxies (e.g. satellite galaxies oflarge clusters found at high redshifts) will be preferentiallymissed.To account for fiber collisions, we define the fiber colli-sion f coll,i for each galaxy, as the fraction of photometricallydefined target galaxies that fall within a area of 55” in ra- c (cid:13)000 The second source of systematic bias is caused by fiber col-lisions. The NYU-VAGC catalogue lists the spectroscopiccompleteness f sp of each galaxy, defined as the fraction ofphotometrically defined target galaxies in the subarea forwhich usable spectra are obtained. The NYU-VAGC cata-logue calculates the average completeness for each of thesesubareas by taking into consideration overlapping plates. Inthe jargon of the NYU-VAGC catalogue, these subareas arecalled sectors . f sp contains information about the missinggalaxies due to lack of fibers in dense regions, missing galax-ies due to spectroscopic failures, and missing galaxies due tofiber collisions. The average f sp for the sample defined aboveis 0.9146. However, f sp assumes that all galaxies with mea-sured spectra are randomly distributed within a sector, andhence cannot account for specific differences between highand low density regions in the same sector. In particular, dueto fiber collisions, certain galaxies (e.g. satellite galaxies oflarge clusters found at high redshifts) will be preferentiallymissed.To account for fiber collisions, we define the fiber colli-sion f coll,i for each galaxy, as the fraction of photometricallydefined target galaxies that fall within a area of 55” in ra- c (cid:13)000 , 1–11 R. D’Souza, S.Vegetti, G. Kauffmann dius. f coll,i takes the values between 0.111 and 1.0 (that is,8 closest neighbours and no neighbours respectively). Theaverage of f coll,i over our whole sample is 0.93819. We nor-malise f coll,i such that it’s average value is the same as thatof f sp . f norm coll,i = fac ∗ f coll,i , where fac is defined as < f sp > / < f coll > and takes the value 0.9749. The nor-malised fiber collision f norm coll now has the general averageproperties of f sp , but can better account for fiber collisions.Weighting by normalised fiber collision maintains thenormalisation of the stellar mass function at the low massend and increases the mass function up to 22% at the highmass end. Li & White (2009) did not include fiber collisionsin their derivation and we can only reproduce their stellarmass function by using Petrosian magnitudes and by ne-glecting the effect of f sp . The third main source of systematic error is related to theassumption about the passive evolution of galaxies both intheir number density and luminosity. In order to constructa stellar mass function from a large redshift range (( . (cid:54) z (cid:54) . ), we would need account for the passive evolution ofgalaxies using a so-called evolutionary correction. Assumingsuch a uniform evolutionary correction is problematic, sincegalaxy evolution is a function of galaxy type and cannotbe described by a simple linear model. For example, star-forming galaxies will evolve more slowly in luminosity thanearly-type galaxies.In order to quantify the effects on the stellar massfunction related to the assumptions about galaxy evolution,we consider two approaches. In the first approach, we as-sume a uniform evolutionary correction ( Q r = 1 . ), whichwould represent an upper limit for the evolution of early-type galaxies with high stellar masses and stellar popula-tions that evolve passively with time (i.e. in the absence ofany mergers). In the second approach, we derive the stel-lar mass function without evolution in three redshift slices: . (cid:54) z (cid:54) . , . (cid:54) z (cid:54) . and . (cid:54) z (cid:54) . .In Figure 6, we plot the stellar mass function derivedusing the MPA-JHU stellar masses, including a uniform evo-lutionary correction, accounting for fiber collisions and forthe additional stellar mass corrections due to the extra lightat large radii (red solid curve). In addition, we also indicatethe mass function calculated in the three redshift slices men-tioned above, without evolution. As seen from Figure 6, theevolutionary correction has only a small effect on the stellarmass function ( ∼ at the massive end). This is relatedto the fact that the luminosity evolution is implicitly foldedinto the derivation of the M/L ratio. Another source of systematic bias is related to binningthe data in calculating the mass function via the /V max method. In particular, this introduces further uncertaintyat the massive end of the mass function due to a combina-tion of the low number statistics and the steep slope of themass function over this mass range. In order to quantify thisuncertainty, we recalculate the mass function with differentvalues for the bin sizes, from 0.05 dex to 0.4 dex. In particu-lar, larger bin sizes tends to bias the slope at the high mass log ( M ∗ /h − ) l og ( Ψ / h M p c − l og M − ∗ ) MPA-JHU Evol+Photo Corr.MPA-JHU 0.001 <= z < 0.15MPA-JHU 0.15 <= z < 0.30MPA-JHU 0.30 <= z < 0.50 Figure 6. The effect of evolutionary corrections on the stellarmass function: The red solid line shows the stellar mass func-tion calculated using MPA-JHU stellar masses corrected for miss-ing light from photometry of stacked galaxies, corrected for fibercollisions and with uniform evolutionary correction Q = 1 . inthe redshift range . (cid:54) z (cid:54) . . Also plotted are the stellarmass function without evolutionary corrections in redshift slices . (cid:54) z (cid:54) . , . (cid:54) z (cid:54) . and . (cid:54) z (cid:54) . dashed blue,green and violet lines. end of the mass function towards shallower values. Reduc-ing the bin size increases the steepness of the slope until asaturation limit of about 0.1 dex. The variation caused bychanges in the bin size around the saturation limit is withinthe uncertainties derived by bootstrapping and within thePoissonian errors. Hence, we calculate the stellar mass func-tion in bins of 0.1 dex. Another source of systematic bias in the stellar mass func-tion is caused by the random errors in the flux and M/L ra-tios of individual galaxies. Such an “Eddington” bias causesthe stellar mass function to be higher in the low-numberdensity part because of scattering from the lower stellarmasses (higher number density). This becomes particularlyacute because of the steepness of the stellar mass functionat higher stellar masses.To correct for this bias, we assume a parametrized formfor the stellar mass function. We convolve this function witha distribution of the uncertainties in the stellar mass. Wethen fit this convolved function to the binned values ofthe stellar mass function calculated from the data using amaximum-likelihood method. The best fit paramteric func-tion is thus our true stellar mass function corrected for theEddington bias. For the parametric function, we assume a c (cid:13) , 1–11 tellar Mass Function Table 2. Parameters of a double Schechter function fit to thestellar mass function of SDSS galaxies. Φ ∗ α log M ∗ ( h Mpc − log M − ) ( h − M (cid:12) )0.008579 -1.082 10.6150.000355 -1.120 10.995 double Schechter function, given by Ψ M d M = (cid:20) Ψ ∗ M ∗ e − M/M ∗ (cid:18) MM ∗ (cid:19) α +Ψ ∗ M ∗ e − M/M ∗ (cid:18) MM ∗ (cid:19) α (cid:21) d M , (6)where Ψ M d M is the number density of galaxies between M and M + d M . This provides a much better fit to the datarelative to a single Schechter. We further assume that theuncertainties in the stellar mass are distributed normally in log ( M ∗ /M (cid:12) ) .To estimate the uncertainties in the stellar mass, wefirst estimate the M/L uncertainties as a function of stel-lar mass from the MPA-JHU database. We find that theaverage uncertainty ∆ log ( M/L ) ranges from 0.08 to 0.1as a function of stellar mass. We then estimate the averageuncertainty in the Model magnitude as a function of stel-lar mass. We find that the average uncertainty in the Model magnitude is ∼ . mag across the stellar mass range con-sidered. Hence the M/L uncertainty is much larger than theflux uncertainty.We find that correcting the stellar mass function for theEddington bias reduces it at the high mass end by as muchas 0.48 dex. In Figure 7, we present our final estimate of the stellar massfunction corrected for missing flux, fiber collisions, evolutionand Eddington bias with that of the original Li & White(2009) in red and the Bernardi et al. (2013) (Sersic-Expfits) stellar mass function in green.We provide a parametric representation of the stellarmass function for stellar masses greater than log( M ∗ /M (cid:12) ) (cid:62) . . The parameters of the double Schechter function arelisted in Table 2. An integration of our stellar mass func-tion for stellar masses greater than log( M ∗ /M (cid:12) ) (cid:62) . givesthe mean comoving stellar mass density of the low redshiftuniverse as φ ∗ = 3 . ± . h Mpc − . This amounts toa 35% increase in the mean comoving stellar mass densitycontributed from the same stellar mass range for the Li &White (2009) stellar mass function. In particular, focussingon the high stellar mass end: the mean comoving stellar massdensity of galaxies with stellar masses log( M ∗ /M (cid:12) ) (cid:62) . is a factor of 3.36 larger than the estimate by Li & White(2009), but is 43% smaller than reported by Bernardi et al.(2013). log ( M ∗ /h − ) l og ( Ψ / h M p c − l og M − ∗ ) Li & White 2009 MPA-JHU Evol+Photo Corr.Bernardi et al. (2013) SersicExp Figure 7. The stellar mass function: The blue solid line showsthe stellar mass function calculated using MPA-JHU stellarmasses, corrected for missing flux, fiber collisions, evolution andEddington bias. The red dot-dashed line shows the original Li& White 2009 stellar mass function calculated using the NYU-VAGC stellar masses based on Petrosian magnitudes. Also shownare the Bernardi et al. (2013) stellar mass function values (green)based on Sersic-Exp fits to individual galaxies. The dashed ver-tical line indicates our lowest stellar mass limit above which thestellar mass function is not affected by surface brightness com-pletness issues. Similar to the galaxy stellar mass function, we also calculatethe galaxy luminosity function using the /V max method.However, more careful attention needs to be paid to the evo-lutionary corrections which affects the luminosity functionnot only via the derivation of V max , but also via the calcula-tion of a galaxy luminosity via equation 1. We calculate theluminosity function using two approaches: in redshift slices( . (cid:54) z (cid:54) . , . (cid:54) z (cid:54) . and . (cid:54) z (cid:54) . )without evolution and using a uniform evolutionary correc-tion of Q r = 1 . . In Figure 8, we present the results of M . r band luminosity function considering Model magni-tudes with photometric corrections from stacking, fiber col-lisions and evolutionary corrections in bins of 0.25 dex. Wealso indicate the luminosity funtion without evolution cor-rections in three redshift slices. A comparision of our resultswith those of Bernardi et al. (2013) would require a morecareful treatment of luminosity evolution which is beyondthe scope of this paper. In this paper, we have shown that stacking similar galax-ies together in volume-limited stellar mass and concentra-tion bins allows one to derive average flux corrections to theSDSS Model magnitudes. In particular, we find that these c (cid:13)000 Model magnitudes. In particular, we find that these c (cid:13)000 , 1–11 R. D’Souza, S.Vegetti, G. Kauffmann M . r − log h l og ( Φ / h M p c − m a g − ) Evolution Corrected LF (z~0.1)Bernardi et al. (2013) SerExp (z~0.1)LF: 0.001 <= z < 0.15LF: 0.15 <= z < 0.3LF: 0.3 <= z < 0.5 Figure 8. The luminosity function: The M . r luminosity func-tion calculated with photometric corrections, fiber collisions andflux uncertainty in three redshift slices and assuming an uniformevolutionary correction of Q r = 1 . . We also show the corresponding luminosity function fromBernardi et al. (2013) Sersic-exponential fits. corrections range from 0.02 to 0.31 magnitude, dependingon the stellar mass and concentration of the galaxy.We apply these corrections to the Model fluxes andre-derive the stellar mass function using MPA-JHU stellarmasses, accounting for galaxy evolution corrections and fibercollisions. We find that the slope of the massive end of thestellar mass function is shallower than reported by Li &White (2009), but much steeper than derived by Bernardiet al. (2013).The biggest change in the slope at the massive end ofthe mass function comes from our adoption of the MPA-JHUstellar masses (as much as a 1.24 dex increase at log M ∗ ∼ . M (cid:12) with respect to Li & White 2009). This involves anincrease of 0.57 dex and 0.67 dex due to the changes in fluxand M/L ratio respectively. The second major contributoris the bias caused by the uncertainty in M/L ratio and fluxmeasurements of individual galaxies which accounts for adecrease of ∼ . dex in the mass function at the massiveend. Fiber collisions contributes to an increase of nearly 22%at the massive end. Galaxy evolution corrections accountsfor a decrease of maximum 10% at the massive end of themass function.We also derive the r -band galaxy luminosity functionand obtain similar results. In particular, the biggest sourceof systematic uncertainty in the galaxy luminosity functionis related to the model assumed for the galaxy evolution cor-rection. In this Paper, we use the evolution correction valuesderived by Blanton et al. (2003), which serves as an upperlimit for galaxies at the bright end of the galaxy luminosityfunction. The flux corrections to the SDSS Model magnitude and theirrespective uncertainties derived in this work by stacking mo-saics of similar galaxies in volume limited stellar mass andconcentration bins are consistent with those presented bySimard et al. (2011). We find no evidence for the need oflarge flux corrections of the order of 0.5 magnitudes as pro-posed by Bernardi et al. (2013).Our results are also consistent with extremely deepimaging of nearby early-type galaxies, obtained with theMegaCam camera on the Canada-France-Hawaii Telescopewhich indicate that outer LSB light contributes 5 to 16 per-cent to a galaxy’s total luminosity (Duc et al. 2015). Stack-ing results for luminous red galaxies (average redshift of z ∼ . ) from Tal & van Dokkum (2011) also indicatethat typical SDSS-depth images miss about 20 percent ofthe total stellar light.A number of systematic differences could contribute tothe discrepancy between our results and those by Bernardiet al. (2013). In the limit of low SNR, the determination ofthe sky background level can influence the measured flux ofa galaxy derived from fitting models to the surface bright-ness distribution. The depth of an image limits ones abilityto distinguish between the flux of the outer LSB featuresof the galaxy and the sky background, especially for largestellar mass galaxies at higher redshifts. The use of multi-component models aggravates this problem.The simultaneous estimation of the model parametersand the sky background level may be prone to system-atic bias, since these are often degenerate with each other.Bernardi et al. (2013) use the PyMorph algorithm (based on GALFIT ), which estimates the galaxy flux based on model fit-ting along with a simultaneous estimation of the sky back-ground. Meert et al. (2013) and Meert et al. (2015) havealready highlighted the effect of a bias in the sky subtrac-tion on the total flux of a galaxy. On the other hand, SDSS Photo pipeline estimates the Model magnitudes by first in-dependently estimating and subtracting the local sky back-ground. A similar procedure is followed by Simard et al.(2011). In this work, we use the background subtracted im-ages provided with SDSS DR9 to derive the flux corrections.In addition, the depth of our stacked images allows us to ac-curately determine the residual sky background.Estimating the total flux of a galaxy is dependent onthe exact procedure used for deblending and masking ( seeBlanton et al. 2011 and Simard et al. 2011). In particular,the amount of masking employed has a substantial effect onthe amount of flux that is derived for a specific galaxy. Inthis Paper, we use the conservative masking described byD’Souza et al. (2015), which involves using multiple runs of SExtractor (Bertin & Arnouts 1996).Guo et al. (2010) calculated the stellar mass functionusing the NYU-VAGC stellar M/L ratios and Model mag-nitudes using the methodology of Li & White (2009). Thestellar mass function derived here has a large shift and shal-lower slope than Guo et al. (2010), owing primarily to the useof the MPA-JHU stellar masses and the flux corrections tothe Model magnitudes. The results of our work will affect themajority of recent halo occupation and abundance matchingstudies (e.g. Moster et al. 2013) that use the measurementsof the stellar mass function from Guo et al. (2010). c (cid:13) , 1–11 tellar Mass Function Finally, we comment that the majority of studies of theevolution of the massive end of the stellar mass functionhave found suprisingly little change out to z ∼ (Marastonet al 2013; Moustakas et al 2013; Fritz et al 2014). Theco-moving number density of galaxies with stellar massesgreater than M (cid:12) has apparently remained constant overthe past 9 Gyr, calling into question the late build-up ofthese systems through mergers and accretion. Our work hasshown that a significant fraction of the mass of these systemsmay be “hiding” in low surface brightness outer componentsthat are systematically missed by conventional photometricextraction software. Accurately quantifying the evolution ofthe stellar mass in these halos will be an important challengefor next generation deep imaging surveys. ACKNOWLEDGEMENTS Funding for SDSS-III has been provided by the Alfred P.Sloan Foundation, the Participating Institutions, the Na-tional Science Foundation, and the U.S. Department ofEnergy Office of Science. The SDSS-III web site is .SDSS-III is managed by the Astrophysical ResearchConsortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, theBrazilian Participation Group, Brookhaven National Lab-oratory, University of Cambridge, Carnegie Mellon Uni-versity, University of Florida, the French ParticipationGroup, the German Participation Group, Harvard Univer-sity, the Instituto de Astrofisica de Canarias, the MichiganState/Notre Dame/JINA Participation Group, Johns Hop-kins University, Lawrence Berkeley National Laboratory,Max Planck Institute for Astrophysics, Max Planck Insti-tute for Extraterrestrial Physics, New Mexico State Univer-sity, New York University, Ohio State University, Pennsyl-vania State University, University of Portsmouth, PrincetonUniversity, the Spanish Participation Group, University ofTokyo, University of Utah, Vanderbilt University, Universityof Virginia, University of Washington, and Yale University. REFERENCES Baldry I. K, Driver S. P, Loveday J. et al., 2012, MNRAS,421, 1621Baldry I. K., Glazebrook K., Driver S. P., 2008, MNRAS,388, 945Bell E. F., McIntosh D. H., Katz N., & Weinberg M. D.2003, ApJS, 149, 289Bernardi M., Hyde J. B., Sheth R. K., Miller C. J. et al.2007, AJ, 133, 1741Bernardi M., Meert A., Sheth R. K., Vikram et al. 2013,MNRAS, 436, 697Bertin E., & Arnouts S. 1996, AAPS, 117, 393Bertin E., Mellier Y., Radovich M., et al. 2002, ASPC,281, 228Blanton M. R., Hogg D. W., Bahcall N. A., et al. 2003,ApJ, 592, 819Blanton M. R., Schlegel, D. J., Strauss, M. J., et al. 2005,ApJ, 129, 2562 Blanton M. R., Eisenstein D., Hogg D.W., Schlegel D.J.,Brinkmann, J. 2005, Apj, 629, 143Blanton M. R., Roweis S. 2007, AJ, 133, 734Blanton M. R., Kazin E., Muna D., Weaver B. A., & Price-Whelan, A. 2011, AJ, 142, 31Chen Y., Kauffmann G., Heckman T. M., Tremonti C., etal. 2013, MNRAS, 429, 2643Cole S., Noberg P., et al. 2001, MNRAS, 326, 1, 255D’Souza R., Kauffmann G., Wang J., Vegetti S. 2014, MN-RASD’Souza R., Kauffmann G., Vegetti S. 2015, (in prepara-tion)Duc P., Cuillandre J., Karabal E., Cappellari M., 2015,MNRAS, 446, 120Fritz A., et al., 2014, A&A, 563, A92Guo Q., White S., Li C., Boylan-Kolchin M., 2010, MN-RAS, 404, 1111Kauffmann G., Heckman T. M., White S. D. M., et al.2003, MNRAS, 341, 33Kauffmann G., White S. D. M., Heckman T. M., et al.2004, MNRAS, 353, 713Kelvin L.S., Driver S.P., Robotham A.S.G., Hill D. T. etal. 2012, MNRAS, 421, 1007Lackner C. N., & Gunn J. E. 2012, MNRAS, 421, 2277Li C., White S. D. M. 2009, MNRAS, 398, 4, 2177Maraston C., et al., 2013, MNRAS, 435, 2764Meert A., Vikram V., Bernardi M., 2013, MNRAS, 433,1344Meert A., Vikram V., Bernardi M., 2015, MNRAS, 446,3943Mendel, J. T., Simard L., Palmer M., Ellison S.L., PattonD.R., 2014, ApJS, 210,3Moustakas J., et al., 2013, ApJ, 767, 50Schlegel D. J., Finkbeiner D. P., & Davis M. 1998, ApJ,500, 525Simard L., Mendel J. T., Patton D. R., Ellison S. L., Mc-Connachie A. W., 2011, ApJS, 196,11Tal T., & van Dokkum P. G. 2011, ApJ, 731, 89Von Der Linden A., Best P.N., Kauffmann G., WhiteS.D.M., 2007, MNRAS, 379, 867 c (cid:13)000