The 2.4 $μ$m Galaxy Luminosity Function as Measured Using WISE. III. Measurement Results
S. E. Lake E. L. Wright R. J. Assef T. H. Jarrett S. Petty S. A. Stanford D. Stern C. -W. Tsai



0000-0002-4528-7637]S. E. Lake
Physics and Astronomy Department, University of California, Los Angeles, CA 90095-1547

0000-0001-5058-1593]E. L. Wright
Physics and Astronomy Department, University of California, Los Angeles, CA 90095-1547

0000-0002-9508-3667]R. J. Assef
Núcleo de Astronomía de la Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejército Libertador 441, Santiago, Chile


Astronomy Department University of Cape Town Private Bag X3 Rondebosch 7701 Republic of South Africa

0000-0003-0624-3276]S. Petty
NorthWest Research Associates
4118 148th Ave NE
Redmond, WA 98052-5164


Institute of Geophysics and Planetary Physics, Lawrence Livermore National
Laboratory, Livermore CA 94551
Department of Physics, University of California, Davis, CA 95616



Jet Propulsion Laboratory, California Institute of
Technology, 4800 Oak Grove Dr., Pasadena, CA 91109

0000-0002-9390-9672]C.-W. Tsai
Jet Propulsion Laboratory, California Institute of
Technology, 4800 Oak Grove Dr., Pasadena, CA 91109
Physics and Astronomy Department, University of California, Los Angeles, CA 90095-1547




S. E. Lake
lake@physics.ucla.edu


The  satellite surveyed the entire sky multiple times in four infrared wavelengths <cit.>. 
The unprecedented combination of coverage area and depth gives us the opportunity to measure the luminosity function of galaxies, one of the fundamental quantities in the study of them, at 2.4 to an unparalleled level of formal statistical accuracy in the near infrared. 
The big advantage of measuring luminosity functions at wavelengths in the window ≈ 2 to 3.5 is that it correlates more closely to the total stellar mass in galaxies than others. 
In this paper we report on the parameters for the 2.4 luminosity function of galaxies obtained from applying the spectroluminosity functional based methods defined in <cit.> to the data sets described in <cit.> using the mean and covariance of 2.4 normalized spectral energy distributions (SEDs) from <cit.>. 
In terms of single Schechter function parameters evaluated at the present epoch, the combined result is: ϕ_⋆ = 5.8 ± [0.3_stat,  0.4_sys] × 10^-3Mpc^-3, L_⋆ = 6.4 ± [0.1_stat,  0.3_sys] × 10^10  L_2.4 ⊙ (M_⋆ = -21.67 ± [0.02_stat,  0.05_sys]AB mag), and α = -1.050 ± [0.004_stat,  0.04_sys]. 
The high statistical accuracy comes from combining public redshift surveys with the wide coverage from , and the unevenness in statistical accuracy is a result of our efforts to work around biases of uncertain origin that affect resolved and marginally resolved galaxies. 
With further refinements, the techniques applied in this work promise to advance the study of the spectral energy distribution of the universe.








§ INTRODUCTION

The luminosity function (LF) is one of the most basic statistical properties measured for any class of objects in astronomy. 
The fundamental nature of the LF means that it has been measured for galaxies many times, in many different bandpasses <cit.>. 

The release of the AllWISE catalog generated from the data gathered by the  satellite, described in <cit.> and <cit.>, marks the availability of 3.4 (W1) and 4.6 (W2) photometric data that is better than 95% complete over the vast majority of the sky down to 44 and 88 (19.79 and 19.04 AB mag), respectively. 
This new data set presents the opportunity to utilize the large number of public redshift surveys, in tandem with a small -selected survey of our own, to measure the near-IR luminosity function of galaxies at 2.4 to unprecedented accuracy. 
The advantage of measuring the luminosity function in this range of wavelengths is that fluxes suffer from minimal dust extinction in both the target galaxy and the Milky Way, according to dust extinction models like the one from <cit.>. 
Further, near infrared light traces the target galaxy's stellar mass in evolved stars more faithfully than optical wavelengths <cit.>, as long as the contribution of thermally pulsing asymptotic giant branch (TP-AGB) stars can be correctly accounted for in the population synthesis models <cit.>. 
We measure the luminosity function at 2.4, in particular, because it is the wavelength directly observed by W1 for galaxies at the median redshift, z=0.38, of galaxies with F_W1 > 80 in <cit.>.

This new opportunity also presents new challenges.
First, even limiting the redshift surveys to those that are publicly available and that are, primarily, selected by flux at a single wavelength meant that there were a lot of details that needed to be addressed in the characterization and selection process for the six surveys used here. 
Because of this, the primary characterization is done separately in a companion paper <cit.>. 
Second, we limited the surveys to be well above the sensitivity limits of the AllWISE data set, a minimum flux of 80 in W1 (19.14 AB mag), in order to minimize the additional incompleteness from the cross-match and to match the properties of <cit.> as closely as possible. 
This means that the surveys have flux limits at two wavelengths. 
The two flux limits, combined with the wide range of redshifts included in the analysis, 0.01 < z ≤ 1.0, meant that the existing luminosity function measurement tools were not adequate to the challenge of analyzing the collected data set. 
For this reason, a new estimator based on analyzing the likelihood of a galaxy's entire spectral energy distribution (SED) is derived in a companion paper <cit.>. 
The likelihood estimate used in this work is based on the mean and covariance of SEDs as measured in <cit.>. 

The structure of this paper is as follows: Section <ref> summarizes properties of the data used to measure the 2.4 luminosity function described in LW18II, Section <ref> summarizes the estimators used in this paper to perform the measurements, Section <ref> contains the results of the analysis (including comparisons to other measured LFs), Section <ref> places the results of this work in context and outlines possible improvements to the methods used, and Section <ref> contains the conclusions drawn from the analysis.

The cosmology used in this paper is based on the WMAP 9 year ΛCDM cosmology <cit.>[<http://lambda.gsfc.nasa.gov/product/map/dr5/params/lcdm_wmap9.cfm>], with flatness imposed, yielding: Ω_M = 0.2793, Ω_Λ = 1 - Ω_M, redshift of recombination z_recom = 1088.16, and H_0 = 70 ^-1^-1 (giving Hubble time t_H = H_0^-1 = 13.97Gyr, and  Hubble distance D_H =  c t_H = 4.283 Gpc). All magnitudes are in the AB magnitude system, unless otherwise specified. 
In cases where the source data was in Vega magnitudes and a magnitude zero point was provided in the documentation, they were used for conversion to AB (2MASS[<http://www.ipac.caltech.edu/2mass/releases/allsky/faq.html#jansky>] and AllWISE[<http://wise2.ipac.caltech.edu/docs/release/allsky/expsup/sec4_4h.html#WISEZMA>]). 
For the surveys without obviously documented zero points (NDWFS[<http://www.noao.edu/noao/noaodeep/>], SDWFS[<http://irsa.ipac.caltech.edu/data/SPITZER/SDWFS/>]) we converted Vega magnitudes to AB magnitudes using those provided in <cit.>. 
When computing bandpass solar luminosities we utilized the 2000 ASTM Standard Extraterrestrial Spectrum Reference E-490-00[<http://rredc.nrel.gov/solar/spectra/am0/>]. 
For our standard bandpass, W1 at z=0.38, we get an absolute magnitude of M_2.4 ⊙ = 5.337AB mag, L_2.4 ⊙ = 3.344 × 10^-8^2.



§ DATA SUMMARY

The full description of the data sets used in this work can be found in LW18II.  

Table 1 of LW18II summarizes the properties of the spectroscopic surveys used here, and we refer interested readers there.





























As described in LW17I, the parametric model used here requires models for the relationship between noise and flux. 
Because the flux selection cuts are all relatively bright a single power law of the form

    (σ_F/F)^2 = (F/A)^B

was adequate. 
The model also requires an estimate of the mean, ⟨τ⟩, and standard deviation, σ_τ, of the dust obscuration for each selection filter in the field of the redshift survey. 
Table <ref> contains the parameters measured for each survey, with separate values for the optical selection filter and W1. 
The calculation of the incompleteness due to flux variety also uses the parameters in the table to calculate the selection function in redshift-luminosity space, S(L,z).

l|cccccccc
	
	
	
	Noise Models Used for Parametric Fits and Completeness
	 Survey     
		A_opt     B_opt     
		A_W1     B_W1    
		 ⟨τ_opt⟩     σ_τ opt     
		 ⟨τ_W1⟩     σ_τ W1 

		     
		             
		            
		  10^2     10^2     
		  10^3     10^3 
	6dFGS     43410     -0.9068     5.469× 10^-6     -0.2774     2.559     3.142     12.69     15.58

		SDSS     203.7     -1.174     0.9552     -0.5334     0     0     6.173     6.107

		GAMA     66.15     -1.419     6.793     -0.6317     0     0     5.914     1.801

		AGES     37.12     -1.472     1.570     -0.5651     2.002     0.4561     2.015     0.4591

		/DEIMOS     89.52     -0.457     675.6     -1.016     23.91     23.89     15.8     15.79

		zCOSMOS     102.1     -1.379     139.1     -0.8576     3.134     0.2279     3.141     0.2283

	Noise model parameters fit to each survey separately. 
	The `opt'  parameters pertain to the optical filter used for selection in the survey, 
	and the `W1' parameters are the parameters for the W1 filter. 
	The A and B parameters relate to a power law fit of signal-to-noise ratio to flux (see Equation <ref>), and ⟨τ⟩ and σ_τ are the mean and standard deviation of the optical depth from foreground dust obscuration, averaged over targets in the survey. 
	When the dust obscuration parameters are zero, target selection was performed on fluxes after extinction correction instead of before. 
	All quantities given to four significant figures, regardless of the statistical or systematic uncertainties in the quantities. 
	


There is one additional plot set not in LW18II that is necessary for understanding the choices made about how to analyze the data, found in Figure <ref>. 
It shows that w1rchi2, the reduced χ^2 of fitting a point spread function to the source in the W1 images <cit.>[<http://wise2.ipac.caltech.edu/docs/release/allwise/expsup/sec5_3bii.html#review_stationary_model>],  behaves about as one would expect for galaxies with z<1 where, even if galaxies in the past were the same size as today, the increasing angular diameter distance gives them, overall, a smaller radius on the sky. 
The vertical line at z=0.2 is placed there after trial and error as a dividing line between samples that are significantly contaminated by resolved and marginally resolved objects (z ≤ 0.2), and those that are sufficiently point-like to render resolution concerns moot (z > 0.2). 
In this work we describe the set of all galaxies that come from below z=0.2 as the low z sample, and those from above the line as the high z sample. 



As a check on the impact the new analysis techniques from LW17I had on the results, we have also analyzed high and low redshift subsamples that had completeness above 98% of the maximum value for the survey, as defined by Equation 21 of LW17I, and as shown by the light blue contours on the luminosity-redshift plots in LW18II. 
These samples, called the `trim' simples, have significantly lower numbers of galaxies with the accompanying increase in statistical uncertainty, but also reduced systematic uncertainty from the constancy of the selection function. 

The number of sources each survey makes to the combined samples, as well as their overall sizes, can be found in Table <ref>.

l|rrrrrr|r
	
	
	
	Combined Samples Sizes
	Sample     N_6dFGS     N_SDSS    
		N_GAMA     N_AGES     
		N_WDa     N_zCb     N_tot 
	Low z     27,071     450,731     28,619     2,000     36     133     508,590 

		High z     0     26,556     15,872     3,741     171     1,091     47,431 

		Low z Trim     15,891     106,003     9,513     1,212     25     116     132,760 

		High z Trim     0     21     125     452     17     647     1,262

	Number of sources contributed by each survey to the combined samples named in the first column. 
	The last column contains the total number of sources in each combined sample.
	a
	bzCOSMOS
	




§ ANALYSIS METHODS AND MODELS SUMMARY

The analyses in this work are centered around maximum likelihood estimates of the quantities related to the LF. 
As is customary in most works on the LF, we analyze the data using both binned and parametric estimators. 
The full description and derivation of the estimators used here is in LW17I. 
Two standard binned estimators were adapted for this work: 1/V_max from <cit.> as modified by <cit.>; and N_obs/N_mdl from <cit.>. 
The adaptation was to handle situations where the selection function, S(L,z), varies significantly over the bin in luminosity-redshift space. Note that the selection function, S, is here not regarded as a particular function but as an assessment of the probability that a galaxy with the properties specified in the arguments would be selected for inclusion in the data set.


For the 1/V_max estimator, the varying completeness corrected version of the estimator, with log-spaced bins in luminosity, is given by:

    Φ(L_i,z_i) = 1/S(L_i,z_i) L_i Δln L_i∑_j1/Δ V_j,

where the redshift-luminosity bin is labeled by the index i, Φ(L_i,z_i) is the constant estimate of the luminosity function for that bin, the sum is over sources that fall into the bin, and Δ V_j is the volume available to the source to still be in both the bin and the selection criteria of the survey. 

The varying completeness corrected version of N_obs/N_mdl is, formally, much closer to the form of the original estimator from <cit.>:

    Φ(L_i,z_i)     = Φ^mdl(L_i,z_i) ·N_i/⟨ N^mdl_i⟩
        = Φ^mdl(L_i,z_i)/∫_binΦ^mdl(L,z)  S(L,z)ḶṾ· N_i,

where Φ^mdl(L,z) is an approximate (`model') luminosity function that brings the estimator closer to evaluating the luminosity function at the center of the bin, (z_i, L_i), as long as Φ^mdl(L,z) is closer to the true luminosity function than the implicitly assumed Φ^mdl(L,z) = constant of an uncorrected estimator. 
For these purposes, this work assumes Φ^mdl(L,z) is a Schechter function that has a faint end slope α = -1 and M_⋆ = -22 (near the peak of the luminosity histograms in LW18II).

The parametric estimator used in this work is based on the spectro-luminosity functional, Ψ. 
The final estimator of the likelihood of the data is a little complicated, requiring a few nested equations to express. 
The outermost equation is given by:

    ln(ℒ)     = ∑_galaxiesln (S(F_sel, F_0, x⃗) (F_sel, F_0|L_0) Φ(L_0, z)) 
        - ∫ S(F_sel, F_0, x⃗) (F_sel, F_0|L_0) Φ(L_0, z) F̣_selF̣_0 Ḷ_0 Ṿ_c,

where the selection function, S, calculates the probability that a source with given observed fluxes used for primary and secondary target selection, F_sel and F_0 (respectively), and at spatial position x⃗ is selected for inclusion in the survey. 
Here F_sel differs by survey (see the band column of Table 1 in LWII), and F_0 is always the W1 flux because the source needs to be well detected there to use F_W1 in the calculation of L_2.4.
This flux selection function is assumed to take only two values: 0 for excluded regions, and and s for sources in selected regions, where s (0 < s < 1) is the overall completeness of the survey. 
Basically, its primary purpose is to set the limits of the integration. 

 is the likelihood for a galaxy to have a particular SED given it has 2.4 luminosity L_0, projected down from the full function space to the fluxes F_sel and F_0. 
In this work we approximate the full  as a Gaussian and incorporate an uncertainty model for the fluxes and dust extinction to get:

    (F_sel, F_0|L_0)     = ^τ_1 + τ_2/√((2π)^2 det(σ))·(4π D_L(z)^2/(1+z)L_0)^2 ·exp(-1/2∑_i,j=2^2 [ℓ_i - μ_i][σ^-1]_ij·[ℓ_j - μ_j]), and
    ℓ_i    = F_i 4π D_L(z)^2 ^τ_i/(1+z)L_0.

τ_i is the optical depth due to dust in the Milky Way present for flux in filter i, when the selection fluxes weren't extinction corrected, 0 when they were. 
D_L(z) is the luminosity distance to the source. The F_i are one of F_sel or F_0 and are the measured fluxes. 
μ_i is the mean value ℓ_i takes, as predicted from the mean SED. 
σ is a covariance matrix that takes the following form (no summation):

    σ_ij = Σ_ij + (δ_ij A_i F_i^B_i + σ_τ  iσ_τ  j) μ_i μ_j,

where A and B are model parameters fit using ordinary least squares in log-space of σ_F^2 F^-2 to F, σ_τ  i is the standard deviation of the optical depth present for fluxes in channel i over the survey targets, and Σ_ij is the covariance of SEDs projected down to apply to the space spanned by F_sel and F_0, similar to the calculations described in <cit.>. 
The optical depths are calculated from the dust extinction model of <cit.> using the E(B-V) dust maps from <cit.>.

Φ(L_0,z) is a Schechter luminosity function, originally defined in <cit.>,

    Φ(L, z)     = ϕ_⋆/L_⋆(L/L_⋆)^α^-L/L_⋆.

The parameterization for evolution in ϕ_⋆ and L_⋆ are:

    ϕ_⋆    = ϕ_0 ^-R_ϕ t_L(z), and
    
    	L_⋆    = L_0 ^-R_L t_L(z)(1 - t_L(z)/t_0)^n_0,

where t_L(z) is the lookback time of a source at redshift z, R_ϕ is the specific evolution rate in ϕ and is assumed to be constant, R_L is the luminosity evolution rate, t_0 is the time of first light (set here to t_0 = t_L(z_recom) since our data cannot meaningfully constrain it),  and n_0 is the power law index for the early time increase in L_⋆.

As explained in LW17I, the actual fitting was done using derived parameters that proved to be more statistically orthogonal than the ones given above. 
In terms of the parameters from Equation <ref>, they are:

    R_n     = R_ϕ - min(1+α, 0) [R_L + n_0/t_0], and
    κ_⋆    = ϕ_0 L_0^3/2/Ω_sky 4π^1/2Γ(α + 5/2),

which are named the specific rate of change in galaxy number density at z=0, and the normalization to the source flux counts distribution. 
If the luminosity function were a static Schechter function in a static Euclidean universe of infinite radius then Ṇ/F̣Ω̣ = κ_⋆ F^-5/2.



 §.§ Error Analysis Details

The statistical uncertainty in the binned estimators is assumed to be fundamentally Poisson combined with propagation of errors that assumes all factors other than the number of galaxies in a bin are constants, as described in LW17I. 
For the parametric estimator in Equation <ref>, its complicated form makes performing a second order expansion about the maximum likelihood parameters, as is done in frequentist statistics, to find the uncertainty in those parameters tedious and error prone. 
Further, when there is any region of the parameter space where the likelihood becomes unusually flat, as happens frequently when selection functions are involved, the Taylor expansion must be carried out to higher order to get an estimate in the parameters' uncertainties. 
These conditions make the work done here an ideal case for the application of a Bayesian analysis using Markov Chain Monte Carlo (MCMC) to estimate the uncertainties in the parameters. 
The software tool used to perform this error analysis is the Python package known as emcee version 2.1.0[<http://dan.iel.fm/emcee/current/>], described in <cit.>, an implementation of the stretch-move algorithm proposed in <cit.>.

When doing any MCMC analysis the algorithm needs to run for a number of steps before it starts to provide an accurate and uncorrelated sample of the posterior. 
This process is called `burn in.' 
If the initial point is far from the mode of the posterior, then this first step is dominated by a pseudo-random walk toward that mode, effectively making it an inefficient optimization algorithm. 
This process is even less computationally efficient for emcee because it uses an entire ensemble of independent walkers. 
To short circuit this process, we started the walkers in a ball around an extremum found using the optimize package of SciPy. 
When fitting models as complex as the ones used here, there is the added downside that there are frequently multiple local extrema. 
This factor causes a tradeoff in the usage of emcee and how spread out the initial positions of the walkers are: if the spread is large then the odds of finding a better minimum than the current guess goes up, but walkers will also get stuck in local minima that are uninteresting outliers; if the spread is small, then all of the walkers are characterizing the minimum of interest, but it takes much longer to find any possible lower minima.

In principle, it would be possible to design an algorithm that was mostly emcee, but that periodically trimmed outliers and started new walkers at large distances to search for possible new minima. 
In practice, it is easier to break the process up into different steps. 
In minimum finding mode, the initial spread of the walkers is large, centered on the minimum found by a comparatively efficient algorithm, and it restarts any time a new minimum is found that is lower than the initial one. 
The software then switches to minimum characterizing mode where the initial spread is small, and after a number of burn in steps that are discarded the final sample is produced. 
If a new minimum is found during minimum characterizing mode, then the minimum characterizing process begins again from the beginning, but the switch is not made back to minimum finding mode under the assumption that the improvement from doing so is marginal.

With any Bayesian analysis the prior must be described. 
Table <ref> contains the explicit descriptions of the ranges of the parameters, and the priors assumed on the parameters, for the LF. 
In most cases the priors are flat in the given parameter, and the majority of the remainder are flat in the logarithm of the parameters. 
The exceptions to this are the faint end slopes, α, and the initial luminosity index, n_0. 
Because α=-2 gives an unphysical infinite background radiation, we chose to impose a mildly informative prior with a beta distribution shape that excluded the end points of the allowed interval. 
With n_0 the buildup of L_⋆ should neither be discontinuously fast, so n_0 > 0, nor should it be much slower than the integral of a linear accretion rate, so n_0 < 10. 
Further, the data used here does not constrain the buildup of L_⋆ directly, so an informative prior that constrains the value of n_0 is used. 
Explicitly, it is expected that their initial luminosity is proportional to the gas accretion rate, and should therefore be reasonably near linear, so we impose n_0 = 1 ± 0.2 with a log-normal distribution.

cccccc
	
	0.5
	LF Parameter Priors
	 Param     Min     Max     Units     Prior     Notes 
	κ_⋆     10^-3     10^2     Jy^3/2^-1     κ_⋆^-1     a 

		R_n0     -50     50     t_H^-1     flat     b

		L_0     10^9     10^13     L_2.4 ⊙     L_0^-1     c  

		R_L     -50     50     t_H^-1     flat     

		α     -2     0     —     -α(2+α)     

		n_0     10^-15     10     —     n_0^-1e^-25 ln(n_0)^2 / 2     d

	Parameter (Param) ranges (from Min to Max) and priors used in fitting the luminosity function. 
	Priors are given in unnormalized form.
	a Normalization to the static Euclidean number counts. 
	Its relationship to standard Schechter function parameters can be found in Equation <ref>.
	b Specific rate of change of the galaxy number density at z=0. 
	Its relationship to standard Schechter function parameters can be found in Equation <ref>.
	c L_2.4 ⊙ = 33.44 ^2 ↔ M_2.4 ⊙ = 5.337 AB absolute.  
	d This prior has mean 1 and standard deviation 1/5. 
	




§ LUMINOSITY FUNCTIONS

Running emcee produces a sequence of model parameters that are distributed as though they were samples taken from the posterior distribution. 
The sequence of parameters is, therefore, the result of the analysis from which all other results are derived. 
The chains are published with this paper and at <www.figshare.com> under the digital object identifier (DOI) https://dx.doi.org/10.6084/m9.figshare.410962510.6084/m9.figshare.4109625 in https://www.gnu.org/software/gzip/gzipped http://irsa.ipac.caltech.edu/applications/DDGEN/Doc/ipac_tbl.htmlIPAC table format. 
An few example lines from one of the chains can be found in Table <ref>.

rrcccccc
	
	
	0.97
	
	Example Lines from MCMC Chains
	  StepNum     WalkerNum    
			ln_KappaStar     R_n      ln_Lstar      R_L     alpha     ln_n0 

		—     —     
			ln(^3/2^-1)      ^-1     ln(^2)     ^-1     —     — 
	0           0          2.0414       -0.056289      8.7190     -0.13463       -1.5071       -0.30514 
   
		0           1          1.8404        -0.013686      8.2423     -0.19569       -1.2721       -0.03764  

		0           2          2.1309        0.224666      7.6358      -0.29812      -0.9931      -0.01672   

		0           3          2.5518        0.095079      8.6992      -0.22876      -1.4027       0.41558   

		0           4          1.8430        -0.007250      8.2185      -0.21643       -1.3473       0.03394   

	Example lines from one of the chains produced by emcee in the tables under DOI https://dx.doi.org/10.6084/m9.figshare.410962510.6084/m9.figshare.4109625. 
	Floating point values truncated here for brevity, but not in the downloadable tables. 
	StepNum is the zero indexed step number that the ensemble was at in the chain, and WalkerNum is the number of the walker which was at the position defined by the row for that step. 
	ln_KappaStar is the natural logarithm of κ_⋆ in ^3/2^-1 (see Equation <ref>). 
	R_n is the specific rate of change of the galaxy number density evaluated at the present time in ^-1 (see Equation <ref>). 
	ln_Lstar is the natural logarithm of L_⋆ in ^2 evaluated at z=0 ( M_⋆ = -2.5/ln 10ln[L_⋆/3631 4π 10^-10^2 ]≈ -1.086ln[L_⋆/^2 ] -13.352  for AB absolute mag). 
	R_L is the long term decay constant in L_⋆ in ^-1. 
	alpha is the faint end slope of the luminosity function. ln_n0 is the natural logarithm of the early time power law index of the evolution of L_⋆ (see Equation <ref>) 
	


The content of that file set is as follows: there is an individual chain for each individual spectroscopic survey (with file names matching the names of the survey), one for each merged subsample described in Table <ref> with matching names, and two chains, labeled `High z Prior' and `High z Trim Prior', where the High z samples were analyzed with the mean and standard deviation in α from the matching Low z samples used as additional Gaussian priors in the analysis. 
The reason for only using the faint end slope as a prior is discussed in Section <ref>. 
Of all of the analysis chains produced, High z Prior is the canonical one for this work based on a subjective evaluation of statistical accuracy and bias. 
A comparison of the N_obs / N_mdl binned estimator using all of the data to the LF with the mean parameters from the High z Prior chain can be found in Figure <ref>. 
The solid red lines, of varying opacity and brightness, are the result of evaluating the High z Prior LF model at five equally spaced redshifts from 0.01 to 1 (0.01, 0.258, 0.505, 0.743, 1). 
The three 3.6 LFs from <cit.>, color corrected to 2.4 using the mean SED from <cit.>, are plotted at z=0.38 using black dashed (`all'), red dotted (`early'), and blue dash-dotted (`late') lines.

Two features stick out most prominently in Figure <ref>. 
First, the parametric estimator places the value of ϕ_⋆ a factor of about 1.7 higher than the binned estimator, though this is within the statistical uncertainty in ϕ_⋆ for Hi z Prior. 
Second, the falloff at the bright end appears to be better described as a power law than exponential. 
The former is likely caused by the fact that the data content of Figure <ref> is drawn from all the samples, without restriction on redshift, while the plotted LF is of the mean parameters from the High z Prior chain. 
The cause of the latter is uncertain. 
Averaging the evolving fit LF over different redshifts using the observed redshift distribution did not produce the observed power law shape. 
Further, limiting the plot to data in the High z sample does not alter this feature significantly, either, so it is unlikely to be an artifact related to photometry of resolved sources. 
The explanation that, qualitatively, seems most likely to cause the feature is the presence of active galactic nuclei (AGN) in the sample, which are observed to have a power law falloff to the LF on the bright end <cit.>. 
While this means that the fit LF doesn't match all of the details of the real LF, this was an expected consequence of using a single Schechter function for the LF and not separate LFs for different galaxy types. 
The global properties of the LF, like predictions of galaxy count and luminosity density, should still correspond with the observable values in the same way that a Gaussian fit to a data set will reproduce the mean and standard deviation, even if the data are not Gaussian distributed.

Detailed analyses of the different posterior chains is done in Subsection <ref>, including an examination of the necessity of splitting the combined samples by redshift, and evidence for luminosity uncertainty contributing to the softening of the high luminosity falloff. 
Comparisons with the results of other measurements of the luminosity function is done in Subsection <ref>.





 §.§ Internal Comparisons

Because the parametric estimator used in this work is new, it is important to analyze multiple data set that all have different selection criteria and compare the results to see if the systematic biases have been correctly managed. 
The mean parameters from each of the posterior chains can be found in Table <ref>, and parameters derived from those mean parameters are found in Table <ref>. 

All of the uncertainties given in the tables are purely statistical uncertainties derived from the Bayesian posterior of the data. 
They do not include sources of error that are, for the purposes of this work, systematic, including: the accuracy of the cosmological parameters (≈ 3%), the peculiar velocity the Milky Way , the accuracy of the completeness assessments of the different surveys , selection effects not modeled, the accuracy of the AllWISE W1 photometric zero point <cit.>, and the accuracy of the numerical integration algorithms used . 
Cosmic variance, the additional variability of the data set induced by the tendency of galaxies to cluster more than random chance, was estimated to be about 5% for SDSS in <cit.>, and <cit.> found GAMA to be 15% under-dense with respect to SDSS. Based on this, we assign cosmic variances of: 15%, 5%, 15%, 20%, 20%, and 20% to the 6dFGS, SDSS, GAMA, AGES, , and zCOSMOS surveys, respectively, with 4% for the combined analyses.
The combination of these effects implies 5.3% systematic uncertainty in the determination of L_⋆, 5.2% in κ_⋆, 7.3% in ϕ_⋆, an assumed -fold per Hubble time in the evolution rate parameters, and an assumed 4% in α, for the combined samples.

All of the parameters on a given line in Tables <ref> and <ref> are correlated, to greater or lesser degrees. 
Including tables or plots of the correlation among the parameters would take up a prohibitive amount of room, so this work only contains a single example of the covariance matrix among the primary parameters constructed from the High z Prior chain in Table <ref>. 
It should be noted that the full correlation among parameters is not necessarily encapsulated by a covariance matrix, particularly the correlation between ln n_0 and R_n. 
In other words, an examination of the pairwise distribution of parameters in each Markov chain does not always show the elliptical structure that would suggest they are well characterized by a Gaussian distribution, and hence by only a mean vector and covariance matrix, particularly when the sample size is small.


The biggest trend in Table <ref> is in R_n, with the lower depth surveys consistent with the number density of galaxies currently declining and the higher depth surveys with the opposite. 
While it is possible that this is a real feature of the data, whether from a turnover in the comoving number density of galaxies or cosmic variance from the Milky way existing in an low density region, selection biases must first be ruled out. 
A preliminary analysis of the SDSS subset with a W1 flux maximum of 2  (15.65AB mag) suggested that this effect was being driven by the galaxies with high apparent flux, and not the faint galaxies that dominate number counts for samples with minimum luminosity significantly lower than L_⋆. 
Because faint galaxies dominate the number density of galaxies, it is unlikely that a trend driven by the presence of bright objects in the sample is a real phenomenon, neither cosmic trend nor cosmic variance.

As the error bars on R_n show, this parameter is poorly constrained by the data, so a small bias in a correlated parameter, even a weekly correlated one, can drive a big change in R_n. 
Since bright galaxies are also more likely to be resolved or marginally resolved, we decided to work around the problem by dividing the combined samples at z=0.2 where the number of galaxies with high w1rchi2 fell to a level low enough to be negligible (see Figure <ref>). 
The advantage of a redshift split instead of a flux cut is that it does not increase the impact of the systematic uncertainties inherent in the the final estimator used. 
The only parameter from the untrimmed low redshift sample analysis that is plausibly not affected by a bias affecting resolved sources is the faint end slope. 
Given the luminosity range covered by the Low z samples makes it the strongest constraint on α available, we use the mean and standard deviation of α from the corresponding Low z samples as Gaussian priors on the High z Prior samples. 

In order to estimate the possible systematic impact the Trim samples were constructed. 
The Trim samples are identical to their corresponding combined samples, but each survey is limited to the region in luminoisty-redshift space where the completeness is at least 98% of its maximum value; the regions enclosed by the faint blue lines in the L-z plots of LW18II. 
This substantially reduces the size and depth of the sample, so it is more vulnerable to cosmic and statistical variance. 
The Low z and Low z Trim chains agree within the statistical uncertainties, with the exception of L_⋆. 
Even though the effects on the High z samples was more dramatic, this is to be expected given the reduced effective depth and loss of low luminosity sources, especially in the AGES and  samples. 
When the power law part of the LF is not directly sampled, where the LF is linear in a log-log plot, the information about α is encoded in the higher order moments of what was observed, increasing sensitivity to statistical and cosmic variance fluctuations. 

lcccccc 
	
	
	
	
	Luminosity Function Bayesian Mean Parameters
	 Survey      κ_⋆     R_na      L_⋆a      R_L     α     n_0 

		—      ^3/2^-1     t_H^-1     10^10 L_2.4 ⊙ h^-2     t_H^-1     —     — 
	6dFGS     2.94 ± 0.07     -4.9 ± 0.8     3.15 ± 0.07     -3.4 ± 0.3     -0.91 ± 0.02     0.9 ± 0.2

		SDSS     4.06 ± 0.03     -3.55 ± 0.09     3.34 ± 0.03     -2.4 ± 0.2     -0.957 ± 0.005     1.0 ± 0.2

		GAMA     4.0 ± 0.1     0.2 ± 0.2     2.67 ± 0.09     -12.7 ± 0.8     -0.95 ± 0.02     8.3 ± 0.6

		AGES     5.4 ± 0.3     1.8 ± 0.3     2.2 ± 0.2     -5.0 ± 0.7     -0.57 ± 0.06     1.8 ± 0.4

		/DEIMOS     8 ± 2     1 ± 2     5 ± 2     -3 ± 1     -1.1 ± 0.3     1.0 ± 0.2

		zCOSMOS     4.4 ± 0.5     1.8 ± 0.8     2.6 ± 0.5     -4.5 ± 0.5     -0.9 ± 0.1     0.9 ± 0.2

		Low z     3.11 ± 0.02     -7.47 ± 0.09     3.41 ± 0.02     -1.5 ± 0.2     -1.059 ± 0.004     0.8 ± 0.2

		High z     5.0 ± 0.2     0.6 ± 0.2     2.66 ± 0.05     -2.8 ± 0.1     -0.68 ± 0.03     0.56 ± 0.08

		High z Prior     5.7 ± 0.2     0.1 ± 0.2     3.12 ± 0.05     -2.6 ± 0.1     -1.050 ± 0.004     0.50 ± 0.07

		Low z Trim     3.24 ± 0.03     -7.7 ± 0.2     3.62 ± 0.04     -0.4 ± 0.2     -0.972 ± 0.008     0.9 ± 0.2

		High z Trim     9 ± 1     -0.8 ± 0.7     3.9 ± 0.4     -4.8 ± 0.5     -1.93 ± 0.04     1.0 ± 0.2

		Hi z Trim Prior     5.0 ± 0.6     1.9 ± 0.7     2.9 ± 0.2     -4.3 ± 0.4     -0.935 ± 0.008     1.0 ± 0.2

	Mean parameters from the Bayesian posterior functions. 
	The top half of the table is broken down by survey, and the bottom half is one of the combined analyses of all data sets. 
	κ_⋆ is the Euclidean flux counts normalization (see Equation <ref>). 
	R_n is the specific rate of change of the numeric density of galaxies (see Equation <ref>). 
	R_L is the long time decay constant in L_⋆, and n_0 is the initial luminosity index (see Equation <ref>). 
	For κ_⋆, L_⋆, and n_0 the means are geometric means, in keeping with the values published in the posterior chains. 
	a Parameter evaluated at z=0. 
	


lccccccc 
	
	
	
	
	Luminosity Function Derived Parameters
	  Survey      ϕ_⋆a     R_ϕ     M_⋆a - 5 log_10 h     z_⋆     ρ_L2.4a     R_ρa     z_ρ 

		—      10^-2h^3 ^-3     t_H^-1     AB mag     —     10^8 L_2.4 ⊙^-3      t_H^-1     — 
	6dFGS     0.86 ± 0.03     -5.2 ± 0.7     -20.91 ± 0.03     1.7 ± 0.3     1.81 ± 0.04     -7.6 ± 0.6     4.1 ± 0.6

		SDSS     1.09 ± 0.01     -3.61 ± 0.09     -20.973 ± 0.008     1.0 ± 0.1     2.49 ± 0.01     -4.93 ± 0.07     2.8 ± 0.4

		GAMA     1.51 ± 0.07     0.0 ± 0.2     -20.73 ± 0.04     0.43 ± 0.01     2.76 ± 0.05     -4.2 ± 0.2     0.43 ± 0.02

		AGES     2.4 ± 0.2     0.4 ± 0.2     -20.53 ± 0.09     1.3 ± 0.2     3.4 ± 0.2     -2.8 ± 0.4     1.1 ± 0.1

		/DEIMOS     1.2 ± 0.8     1 ± 2     -21.4 ± 0.5     1.3 ± 0.6     5 ± 1     -1 ± 1     0.5 ± 0.7

		zCOSMOS     1.8 ± 0.5     1.5 ± 0.5     -20.7 ± 0.2     2.3 ± 0.3     3.1 ± 0.3     -2.0 ± 0.4     1.5 ± 0.3

		Low z     0.814 ± 0.009     -7.43 ± 0.09     -20.994 ± 0.008     0.65 ± 0.09     2.013 ± 0.009     -8.09 ± 0.06     4.8 ± 0.6

		High z     1.81 ± 0.08     -0.1 ± 0.2     -20.73 ± 0.02     2.4 ± 0.2     3.0 ± 0.1     -2.4 ± 0.2     2.5 ± 0.3

		High z Prior     1.69 ± 0.08     0.2 ± 0.2     -20.90 ± 0.02     2.5 ± 0.3     3.8 ± 0.1     -1.9 ± 0.2     2.3 ± 0.3

		Low z Trim     0.77 ± 0.01     -7.7 ± 0.2     -21.06 ± 0.01     -0.7 ± 0.2     1.92 ± 0.02     -7.1 ± 0.2     4.2 ± 0.6

		High z Trim     1.1 ± 0.3     3 ± 1     -21.2 ± 0.1     2.4 ± 0.4     50 ± 40     -1.1 ± 0.6     0.8 ± 0.5

		Hi z Trim Prior     1.7 ± 0.3     1.7 ± 0.7     -20.80 ± 0.08     2.1 ± 0.3     3.3 ± 0.4     -1.6 ± 0.5     1.2 ± 0.3

	Bayesian mean values related to the luminosity function calculated from the posterior chains separately from the parameters in Table <ref>, with h=0.7. 
	The means of ϕ_⋆ and ρ_L2.4 are geometric means. 
	R_ϕ is the specific rate of change of ϕ_⋆ (assumed constant for all redshifts, see Equations 31 of LW17I). 
	z_⋆ is the redshift at which the model predicts L_⋆ will peak (Equation 32 of LW17I). 
	ρ_2.4 is the present day 2.4 luminosity density of galaxies (see Equation 41 of LW17I), R_ρ is its specific rate of change (Equation 42 of LW17I), and z_ρ is the redshift at which the model predicts j to have peaked (Equation 45 of LW17I). 
	a Parameter evaluated at z=0. 
	


ll|llllll
	
	0.95
	Hi z Sample with Low z Prior Luminosity Function Bayesian Parameter Posterior Covariance
	 Parameter      σ     δlnκ_⋆      δ R_n t_H      δln L_⋆(0)     δ R_L t_H     δα     δln n_0 
	δlnκ_⋆     0.03029     1.000     0.9023     -0.5915     -0.1510     -0.04423     -0.2330

		δ R_n t_H     0.2040     0.9023     1.000     -0.8257     -0.4389     0.06990     -0.02677

		δln L_⋆(0)     0.01524     -0.5915     -0.8257     1.000     0.7285     -0.1243     -0.2852

		δ R_L t_H     0.1210     -0.1510     -0.4389     0.7285     1.000     -0.009374     -0.8508

		δα     0.004266     -0.04423     0.06990     -0.1243     -0.009374     1.000     -0.009308

		δln n_0     0.1431     -0.2330     -0.02677     -0.2852     -0.8508     -0.009308     1.0000 

	The σ column contains the standard deviation of the parameters, and the remaining rows and columns comprise the matrix of correlation coefficients between the parameters, considered pairwise. 
	All values given to four significant figures.
	


The numerical comparisons discussed above provide a nice overview of the behavior of the different sets, but no work is complete without graphical comparisons of the fit models to binned estimators for the data. 
The comparisons of the unbinned evolving model ML fits to the entire samples with two binned estimators, 1/V_max and N^obs / N^mdl, can be found in Figures <ref>, <ref>, and <ref>. 
The first notable feature of the grid of plots in Figure <ref> is that the N^obs/N^mdl agrees better with the fit luminosity functions than the 1/V_max estimator. 
This is to be expected since N^obs/N^mdl is closer to a maximum likelihood estimator, and the fits were done using unbinned maximum likelihood. 
In particular, the N^obs/N^mdl outperforms 1/V_max where the approximation that the selection function is constant with redshift is a bad one. 
The difference is particularly stark for panel b, SDSS, where the shallow optical selection makes SED variability particularly relevant. 
The next obvious feature is the disagreement of the fitted faint end slope with the binned ones in panel d, AGES. 
There is a local maximum in the likelihood with a faint end slope closer to α = -1, but it is not a global maximum. 
There are two reasons this fails to be a global maximum: the low faint end completeness of AGES giving the less luminous galaxy bins bigger error bars, and a fluctuation in the data. 
These facts are more apparent in an examination of the first row of Figure <ref>, where the apparent disagreement vanishes.

The final feature of note in the luminosity function plots is the upturns at the faint ends of panels b (SDSS) and c (GAMA). 
It is likely the same feature that caused <cit.> and <cit.> to use a double luminosity function to fit the data. 
While a double LF would provide better agreement to the data, it is unclear without a deeper examination of the data the extent to which the additional LF is modeling a fundamental feature of the universe (for example, the split between red and blue galaxies) or a cosmic variance fluctuation in the data. 
One example of an even bigger fluctuation can be seen in panel c of Figure 12 of LW18II. 
There is a significant over-density in the Sloan data near its peak at around z=0.75. 
The over-density that causes that bump goes by the name the Sloan Great Wall, discovered in <cit.>, and it is a good example of how large cosmic variance can get.

It is also profitable to compare the observed redshift histograms against the predictions based on the luminosity function and selection function; plots containing such comparisons can be found in Figure <ref>. 
The power of this comparison is that, unlike the binned/unbinned LF comparisons in the earlier figures, the data in this figure need not be limited to sources with measured  fluxes. 
Thus, the comparison between the black histogram and darker lines is not entirely one of a fit with the data it was fit to, but the extrapolation of the LF and selection function with new data. 
The first noteworthy feature of the plots in Figure <ref> is that the High z Prior LF (blue dashed lines) provide more accurate extrapolations, overall, than the individual survey fits (solid red lines). 

That said, there are some features of the extrapolation in Figure <ref> that need explaining. First, the model over-predicts the number of galaxies observed at low z in Panel a, 6dFGS.
The likely dominant culprit there is the bias against unresolved galaxies induced by the use of the 2MASS Extended Source Catalog (2MASS XSC) to produce the target list for 6dFGS, as evidenced by the dearth of sources with w1rchi2 < 3 relative to other surveys in Figure <ref>.
Next, the high redshift tails of of Panels a–c also disagree with the model.
Confusingly, the model under-predicts the number of galaxies in the high z tail of a, and over-predicts the tails of b and c.
Because the high-z tail of the histograms is controlled primarily by the way the selection function limits sources to those with L>L_⋆, this is the same as saying that: one, there are more L>L_⋆ galaxies than the model predicts; and two, that those galaxies are redder than average in the optical for the r selection of SDSS and GAMA to remove them.
This appears consistent with the textbook level knowledge that the galaxies found in clusters are redder and larger than the smaller and bluer field galaxies.
Put simply, this is probably a limitation of the single SED single LF model used in this paper.
We have not performed any quantitative investigations into whether this explanation is sufficient because to do so would be to use multiple luminosity functions with their own mean SEDs to build a more accurate spectro-luminosity functional and then fit that to the data, which is beyond the scope of this work.

We have, however, investigated a number of factors that should contribute to the disagreement, but all of them are either the wrong magnitude or incapable of explaining both of the high z tails: 


  * any  related selection effects (the extrapolated graphs remove all -based selection criteria),

  * contamination of the aperture photometry in the target selection catalogs causing sources that are too faint to be included (number of sources is too small, and makes the under-predictions worse), and

  * the spread of predictions consistent with the uncertainty in the model parameters (the 1-σ band is thinner than the lines in the high z regions). 
 
 























 §.§ External Comparisons

In order to make the Schechter parameters in other papers comparable to the ones measured in this one it was necessary to use the mean SED from this work to color correct their values of L_⋆ and, wherever possible, use the evolution measured in the other papers to bring them all to redshifts 0, 0.38, and 1.5. 
The literature parameters can be found in Table <ref> alongside the mean High z Prior parameters from this work. 
Further, the external works often have measurements of the LF in multiple filters; when that is the case, the observation filter with wavelength closest to the W1's 3.4 was used for the primary parameter comparisons. 
On the whole, this work's estimate of L_⋆ is lower than the literature and its estimate for ϕ_⋆ is higher, but not radically so, especially compared to the spread among the literature values. 
The values for L_⋆ at high redshift (z=1.5) have a much larger spread, making this epoch ripe for studies based on deeper imaging surveys.

The spread in measured specific evolution rates, shown in Table <ref>, is considerably larger than the primary parameters. 
The uncertainties are not included in the table, but they're generally more than 0.1 t_H^-1 and less than 1 t_H^-1. 
The most directly comparable values, the ones at ^0.1z, [3.6], and 2.4, are all largely consistent. 
Most importantly, the specific rate of change in the density of galaxies, R_n, should be the same for all the different surveys in all bandpasses. 
That the spread is so large is likely attributable to a combination of selection effects, the inadequacy of a single Schechter function to describe all types of galaxies in all bandpasses, and the fact that R_n is one of the least well constrained parameters by the data.

The final comparison is a graphical one of the models for the evolution of L_⋆ and ρ_L2.4 from this work to an empirical models from <cit.> and <cit.> in this work's Figure <ref>. 
The model from the review in <cit.> is for the evolution of the cosmic star formation rate density, ψ (Equation 15 there), and the the model from <cit.> is a 1-σ variability band in luminosity density evolution for the  K filter (λ≈ 2.2) scaled using the mean SED (L_K / L_ν(2.4) = 1.199). 
The data used to produce the empirical model from <cit.> was, essentially, scaled luminosity densities measured in the far ultraviolet and far infrared. 
Most importantly, the measurements used to produce the empirical model for ψ spans the range of redshifts from 0 to 8, so their estimate would place the peak of the luminosity density somewhere between redshifts 1.3 and 2.5 based on data. 
The data on ρ_K evolution from <cit.> is primarily below redshift 2, and based on statistically interpolating measured values from the literature, and constructing the band from an ensemble of such interpolations. 
Considering that all of the data used in this work has a redshift below 1, and most of that less than 0.5, the crude evolution model used here does surprisingly well at locating the epoch where L_⋆ and ρ_L2.4 peak. 
Figure <ref> also shows that the peaks in L_⋆ and ρ_L are, possibly, too broad, but this is not surprising given the simplicity of the model and limits of the data used.





lccc|ccc|ccc
	
	
	
	Comparison with Other IR Schechter LF Measurements
	Paper     Obs Band     α     10^-8 ρ_L2.4(0)     ϕ_⋆(0)     ϕ_⋆(0.38)     ϕ_⋆(1.5)     L_⋆(0)     L_⋆(0.38)     L_⋆(1.5) 

	—     —     —     h L_2.4 ⊙^-3          10^-2h^3 ^-3              10^10 L_2.4 ⊙ h^-2     
	<cit.>     K     -1.2±0.2     7±6     1.2±0.8     —     —     5±2     —     — 

		<cit.>     K_s     -1.09±0.06     7.5±0.9     1.2±0.1     —     —     4.5±0.2     —     — 

		<cit.>b     K_s     -1.2 ± 0.1     2±3     0.1±0.1     0.2     0.4     19±9     17     14 

		<cit.>     K_s     -0.86 ± 0.04     16±1     4.2±0.2     —     —     4.1±0.2     —     — 

		<cit.>     c1a     -0.92 ± 0.04     5± 5     1±1     2     2     4.4±0.4     0.28     0.25 

		<cit.>b     K     -0.99±0.04     5.0±0.8     0.65±0.09     0.80     0.63     7.8±0.7     7.8     10.0 

		<cit.>b     K     -1.1±0.2     12±3     1.5±0.6     1.1     0.64     8±3     8     11 

		<cit.>     c1a     -1.12 ± 0.16     3.9±0.9     1.08±0.03     —     —     3.3±0.6     5.1     18 

		<cit.>     K     -0.81±0.04     5.7±0.3     1.66±0.08     —     —     3.7±0.1     5.3     14.9 
 
		<cit.>     z     -1.07±0.02     5.1±0.8     1.3±0.2     1.1     0.63     3.9±0.2     7.0     41 

		<cit.>     K     -1.16±0.04     7±2     0.7±0.1      —     —     9±2     —     —   

		This work     W1     -1.050± 0.004     5.4 ± 0.2     1.69 ± 0.08     1.6     1.5     3.12 ± 0.05     5.6     9.9 

	Paper is the work from which the measurements came, in chronological order. 
	`Obs Band' is the observation band used to calculate the luminosity function. 
	ρ_L2.4(0) is the 2.4 luminosity density at z=0, ϕ_⋆(z) is the LF normalization at redshift z, and L_⋆(z) is the value of L_⋆ at redshift z. 
	Uncertainties from this work do not include systematic uncertainties, and so are underestimated. 
	The conversion factors used are:    and 
	aSpitzer/IRAC channel 1, λ=3.6.
	bParameters here are based on linear interpolation/extrapolation.
	



lcc|cc|cc
	
	0.68
	
	Evolution Rate Comparisons
	Paper     z_0      Band     ∂ln L_⋆ / ∂ t     R_ϕ     R_n     R_ρ 

	            t_H^-1     t_H^-1     t_H^-1     t_H^-1 
	<cit.>     0.1     ^0.1u     -3.4     -2.6     -2.3     -5.9 

		        ^0.1g     -1.6     -0.3     -0.1     -1.9 

		        ^0.1r     -1.3     -0.1     -0.2     -1.4 

		        ^0.1i     -1.3     -0.5     -0.5     -1.8 

		        ^0.1z     -0.6     -1.8     -1.9     -2.4 

		<cit.>     0.25     [3.6]     -0.8     0a     -0.1     -0.8 

		        [4.5]     -0.7     0a     0.0     -0.7 

		<cit.>     0.3     ^0.1r     -1.0     0.4     0.4     -0.5  

		<cit.>     0.13     ^0.1u     -4.8     6.5     6.1     1.8 

		        ^0.1g     -2.2     1.2     0.9     -1.1 

		        ^0.1r     -0.5     -1.4     -1.5     -1.9 

		        ^0.1i     -1.1     0.0     -0.1     -1.2 

		        ^0.1z     -1.3     0.4     0.3     -0.9 

		<cit.>     0.2     ^0.1r     -0.7     -0.7     -0.9     -1.4 

		This work     0.38     2.4     -1.8     0.2     0.1     -1.7 

	Paper is the work from which the measurements came, in chronological order. 
	z_0 is the mean or median redshift of the data in the work, and is the redshift at which the parameters were evaluated to make this table. 
	Band is the passband in which the luminosity function was measured. 
	∂ln L_⋆ / ∂ t is the specific rate of change of L_⋆, measured in units of inverse Hubble times. 
	R_ϕ is the specific rate of change of ϕ_⋆ (see Equations 31 of LW17I), R_n is the specific rate of change in the number density of galaxies (see Equation <ref>), and R_ρ is the specific rate of change in the 2.4 luminosity density (see Equation 42 of LW17I) 
	aThis parameter was set to this value, not measured.
	






§ DISCUSSION
 
The consistency of the results in this work with the literature shows that, while the goal of increasing the statistical accuracy of the LF measurements has been met, that increased accuracy has not, yet, uncovered any new facets of galaxy evolution.
The internal comparisons show that there is still room for improvement in the techniques used here. 
In particular, closer attention paid to ensuring that the full flux of resolved galaxies is measured without contamination from foreground stars will permit a further large jump in sample size. 
Improvement in the performance of numerical integration of arbitrary high dimensional Gaussian functions over rectangular regions would make spectro-luminosity functional, Ψ, based techniques able to analyze samples constructed using the sort of complicated color selection done for DEEP2 <cit.>, the SDSS luminous red galaxy sample <cit.>, and the SDSS quasar sample <cit.>. 

Improvements in the form of Ψ that could allow it to fit the data more closely need a little more thought to the overall approach of the analysis before implementation. 
The primary approach used in the measurement of LFs is to classify galaxies into types and measure separate LFs for each type. 
It is not obvious, but the improvements suggested in LW17I, where the total Ψ is written as a sum over components, falls into the classification category. 
The way that the sum over components with Ψ becomes a classification scheme, effectively, is because the fluxes measured for each source will place it closest to one mean SED. 
If we define the square distance to the SED as the exponent in  then the majority of sources will most strongly affect the the LF parameters that correspond to the SED to which they are closest by that distance measure – with a few in a boundary region dividing their influence among the terms. 
Thus, especially if the mean SEDs are fixed, the model effectively classifies the sources by which term in the total spectro-luminosity functional the source has the greatest impact on. 
The more fluxes per source brought to bear in the analysis, the sharper that divide between source types is. 

The less common approach would be to, instead of classifying galaxies, analyzing their composition. 
Think of it as the difference between deciding whether a galaxy `is a' versus how much the galaxy `has a'. 
The advantage of the composition approach is that it addresses an ambiguity not dealt with directly by the model for estimating Ψ developed here: how many galaxies does each source represent? 
We know from images of local galaxies that the moderate to large galaxies have smaller galaxies in their halos (for example: the Large and Small Magellanic Clouds). 
As distance increases, the light from any satellite galaxies must, inevitably, be merged into the light from their primaries. 
The hallmark of a compositional approach is that it does not just assign probabilities to a source being in different classes, it divides the source's luminosity among them. 
So, for example, it could be possible to talk about dividing the Ψ into massive stars (say O, B, and A), intermediate stars (F and G), light stars (K, M, and lighter), stellar and supernova remnants, nebular emission, and AGN. 
The primary challenge would be to figure out what form the base luminosity function, Φ, of each of these should take. 
While it would seem that an increased richness of photometric data would also be required in order to analyze the composition of each galaxy, it isn't necessary to get a detailed analysis of every galaxy to get an accurate picture of the average composition of galaxies that is encoded in Ψ.

One question that is answerable by improving the techniques developed here, and adding data for galaxies in the redshift range of 1 to 3 is: which peaked first, L_⋆ or ρ_L? 
The reason this question is of interest is because it encodes information about the comparative rate of star formation versus galaxy mergers; the former boosts both quantities, the latter only boosts L_⋆. 
We would, for example, expect L_⋆ to peak first in a universe where star formation continued in small galaxies after it and collisions slowed down for large ones. 
The converse would be the case if star formation was a relatively short epoch that cut off in most galaxies and the majority of individual galaxy luminosity growth was through accretion. 
The data  and simple evolution models used here cannot answer this question, though they hint that L_⋆ peaked first.

Another question that will require even greater care to answer is: when did (or will) the number density of galaxies peak? 
The question is equivalent to asking when R_n=0. 
Though that is an equation that can, in principle, be solved using the parameters in this work, the spread in values produced are so large as to make the answer meaningless. 
It would take an evolutionary model with greater physical fidelity combined with more data to provide an answer worth examining.



§ CONCLUSION

The combination of the six different redshift surveys described in LW18II, made possible by the analysis techniques derived in LW17I, has produced measurements of the Schechter LF parameters that are comparable to the literature in terms of statistical precision for ϕ_⋆, and a marked improvement for L_⋆ and α. 
The parameters describing the evolution of L_⋆ and ϕ_⋆, R_L and R_ϕ, are less well constrained, but still comparable to the literature. 
Improving the constraints will require refinements in the process from end to end.

The photometry of bright, resolved and marginally resolved, galaxies requires improvements that bring them in line with the quality of our photometry for point source. 
In the ideal case, a photometric survey with sufficient sensitivity to resolve the wings of the Airy profile produced by even the largest galaxies, and the software tools needed to remove foreground contamination from stars, would guarantee that close enough to all of the light in the galaxy has been directly observed to measure accurate luminosities and colors for each whole galaxy. 

The accuracy of the spectro-luminosity functional, Ψ[L_ν](z), can be improved in the straightforward way: writing it as a sum of spectro-luminosity functionals, each with its own LF and Gaussian . 
It may also be possible to improve the performance of Ψ by a rethinking how the LF is defined – instead of classifying galaxies into mutually exclusive categories, split into constituent parts with their own luminosity. 
The upside of a such an approach is that it naturally handles cases where multiple unresolved galaxies are contained in the same object.

Perhaps more important than increasing the fidelity of Ψ is adding a model of the effective radius of galaxies so that surface brightness limits on galaxy selection can be modeled. 
Likewise, re-deriving an approximation of the estimator for the likelihood of observing an entire catalog so that it includes the effect of galaxy environment. 
Adding the effect of environment, for example using the two point function (ξ(r)), in the likelihood of the catalog is the most natural way to introduce cosmic variance to the process.

Finally, measuring fluxes and redshifts for even fainter galaxies using instruments like the Multi-Object Spectrometer for Infra-Red Exploration (MOSFIRE) on the Keck II telescope, and the James Webb Space Telescope (JWST), will allow for explorations of the evolution of the faint end slope of galaxies, better constrain the evolution of ϕ_⋆ and L_⋆ with more high redshift data, and even, potentially, find a downturn in the LF at faint luminosity where galaxies and star clusters overlap as gas accreting gravitationally bound systems.



We would like to thank the  team.

This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, and NEOWISE, which is a project of the Jet Propulsion Laboratory/California Institute of Technology. WISE and NEOWISE are funded by the National Aeronautics and Space Administration.

We would like to thank the SDSS team.

Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. The SDSS-III web site is http://www.sdss3.org/.

SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofisica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University.

We would like to thank the GAMA team.

GAMA is a joint European-Australasian project based around a spectroscopic campaign using the Anglo-Australian Telescope. The GAMA input catalogue is based on data taken from the Sloan Digital Sky Survey and the UKIRT Infrared Deep Sky Survey. Complementary imaging of the GAMA regions is being obtained by a number of independent survey programs including GALEX MIS, VST KiDS, VISTA VIKING, WISE, Herschel-ATLAS, GMRT and ASKAP providing UV to radio coverage. GAMA is funded by the STFC (UK), the ARC (Australia), the AAO, and the participating institutions. The GAMA website is http://www.gama-survey.org/ .

Based on observations made with ESO Telescopes at the La Silla or Paranal Observatories under programme ID 175.A-0839.

We would like to thank the 2MASS team.

This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.

We would like to thank the MAST team.

Some/all of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX13AC07G and by other grants and contracts.

We would like to thank the NDWFS team.

This work made use of images and/or data products provided by the NOAO Deep Wide-Field Survey (Jannuzi and Dey 1999; Jannuzi et al. 2005; Dey et al. 2005), which is supported by the National Optical Astronomy Observatory (NOAO). NOAO is operated by AURA, Inc., under a cooperative agreement with the National Science Foundation.

We would like to thank the IPAC team.

This research has made use of the NASA/ IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

We would like to thank the GALEX team.

Based on observations made with the NASA Galaxy Evolution Explorer. 
GALEX is operated for NASA by the California Institute of Technology under NASA contract NAS5-98034.

We would also like to thank the teams behind 6dFGS, AGES, zCOSMOS, SDWFS, and COSMOS.

RJA was supported by FONDECYT grant number 1151408.


















