We present the bivariate neutral atomic hydrogen (HI)---stellar mass function (HISMF) (phi(M_HI, M_*)) for massive (log M_*/M_sun > 10) galaxies derived from a sample of 480 local (0.025 < z < 0.050) galaxies observed in HI at Arecibo as part of the GALEX Arecibo SDSS Survey (GASS). We fit six different models to the HISMF and find that a Schechter function that extends down to a 1% HI gas fraction, with an additional fractional contribution below that limit, is the best parametrization of the HISMF. We calculate Omega_{HI, M_* >10^10} and find that massive galaxies contribute 41% of the HI density in the local universe. In addition to the binned HISMF we derive a continuous bivariate fit, which reveals that the Schechter parameters only vary weakly with stellar mass: M_HI^*, the characteristic HI mass, scales as M_*^0.39, alpha, the slope of the HISMF at moderate HI masses, scales as M_*^0.07, and f, the fraction of galaxies with HI gas fraction greater than 1%, scales as M_*^-0.24. The variation of f with stellar mass should be a strong constraint for numerical simulations. To understand the physical mechanisms that produce the shape of the HISMF we redefine the parameters of the Schechter function as explicit functions of stellar mass and star formation rate to produce a trivariate fit. This analysis reveals strong trends with SFR. While M_HI^* varies weakly with stellar mass and SFR, alpha is a stronger function of both stellar mass and especially star formation rate. The HISMF is a crucial tool that can be used to constrain cosmological galaxy simulations, test observational predictions of the HI content of populations of galaxies, and identify galaxies whose properties deviate from average trends.
تحميل البحث