¿Cómo creo una función de masa estelar de galaxia?

Su enfoque es completamente correcto, solo tenga en cuenta tres cosas:

Distribución logarítmica

Primero, dado que la distribución de masas es de naturaleza logarítmica (como la mayoría de las cosas), asegúrese de agruparlas logarítmicamente. De lo contrario, sobremuestreará (submuestreará) los contenedores en el extremo de baja (alta) masa.

Densidades de movimiento

En segundo lugar, para poder comparar funciones de masa en diferentes corrimientos al rojo, se utiliza el volumen comomóvil en lugar del volumen físico, de modo que se excluye la expansión del Universo. Los dos están relacionados como$V_\mathrm{com} = V_\mathrm{phys}(1+z)^3$.

Maldito seas, pequeño$h$!

Finalmente, los observadores y los modeladores tienden a usar una definición ligeramente diferente de la unidad de volumen. Mientras que los observadores suelen utilizar$\mathrm{Mpc}$para las distancias, y por lo tanto$\mathrm{Mpc}^{-3}$para densidades numéricas, si sus galaxias provienen de una simulación cosmológica donde los parámetros cosmológicos pueden ajustarse a voluntad, es costumbre factorizar la constante de Hubble$H_0$. En las simulaciones, las masas y las distancias se miden en$h^{-1}M_\odot$y$h^{-1}\mathrm{Mpc}$, respectivamente, por lo que las densidades numéricas se miden en$h^3 \mathrm{Mpc}^{-3}$. Aquí$h\equiv H_0\,/\,(100\,\mathrm{km}\,\mathrm{s}^{-1}\,\mathrm{Mpc}^{-1})$.

Esto es probablemente una reminiscencia de la época en que la constante de Hubble era bastante incierta. Hoy en día, en mi opinión, no hay necesidad de hacer esto, pero como todo el mundo lo hace, es difícil ir contra la corriente. Para una discusión de este tema, ver Croton (2013) .

código pitón

Dado que ha etiquetado la pregunta con python , escribí este pequeño fragmento que debería hacer el trabajo (elegí al azar$10^5\,\mathrm{Mpc}$como el volumen de su encuesta; tenga en cuenta también que en este ejemplo no factorizo$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')