How do I create a galaxy stellar mass function?

Question

Suppose that I have an array of stellar masses of tens of thousands of galaxies, and that I also already know what the total "survey volume" is (in Mpc^3) within which all of those galaxies are contained. For context, my data is from a simulation of galaxy formation at redshift ~ 0, rather than from an observational simple -- so in principle I don't think I need to worry about incompleteness corrections at the faint (low-mass) end.

Given my array of stellar masses and the total survey volume in Mpc^3, how do I create and plot the "galaxy stellar mass function" that I always see people talking about? Isn't it just basically a histogram of the stellar mass array, and then divide the # of galaxies in every bin by the survey volume (it's almost just like a normalization constant)?

(Of course this question could equally apply to galaxy luminosity functions.)

Answer 1

Your approach is completely correct, just note three things:

Logarithmic distribution

First, since the distribution of masses is logarithmic in nature (as is most other things), be sure to bin them logarithmically. Otherwise you will oversample (undersample) the bins at the low-(high-)mass end.

Comoving densities

Second, to be able to compare mass functions at different redshifts, ones uses the comoving volume rather than the physical volume, such that the expansion of the Universe is factored out. The two are related as $V_\mathrm{com} = V_\mathrm{phys}(1+z)^3$.

Damn you, little $h$!

Finally, observers and modellers tend to use a slightly different definition of the unit volume. Whereas observers usually use $\mathrm{Mpc}$ for distances, and hence $\mathrm{Mpc}^{-3}$ for number densities, if your galaxies come from a cosmological simulation where the cosmological parameters can when tuned at will, it is custom to factor out the Hubble constant $H_0$. In simulations, masses and distances are then measured in $h^{-1}M_\odot$ and $h^{-1}\mathrm{Mpc}$, respectively, so number densities are measured in $h^3 \mathrm{Mpc}^{-3}$. Here $h\equiv H_0\,/\,(100\,\mathrm{km}\,\mathrm{s}^{-1}\,\mathrm{Mpc}^{-1})$.

This is probably a reminiscence of the time when the Hubble constant was rather uncertain. Nowadays, in my opinion there is no need to do this, but since everybody does, it's difficult to go against the stream. For a discussion of this issue, see Croton (2013).

Python code

Since you've tagged the question with python, I wrote this little snippet that should do the work (I randomly chose $10^5\,\mathrm{Mpc}$ as your survey volume; note also that in this example I don't factor out $h$):

import numpy as np
import matplotlib.pyplot as plt

M     = np.loadtxt('Mstar.dat')         #Read stellar masses in Msun
logM  = np.log10(M)                     #Take logarithm
nbins = 10                              #Number of bins to divide data into
V     = 1e5                             #Survey volume in Mpc3
Phi,edg = np.histogram(logM,bins=nbins) #Unnormalized histogram and bin edges
dM    = edg[1] - edg[0]                 #Bin size
Max   = edg[0:-1] + dM/2.               #Mass axis
Phi   = Phi / V / dM                    #Normalize to volume and bin size

plt.clf()
plt.yscale('log')
plt.xlabel(r'$\log(M_\star\,/\,M_\odot)$')
plt.ylabel(r'$\Phi\,/\,\mathrm{dex}^{-1}\,\mathrm{Mpc}^{-3}$')
plt.plot(Max,Phi,ls='steps-post')