Exploring over 700 massive quiescent galaxies at z = 2-7: Demographics and stellar mass functions

High-redshift ($z>2$) massive quiescent galaxies are crucial tests of early galaxy formation and evolutionary mechanisms through their cosmic number densities and stellar mass functions (SMFs). We explore a sample of 743 massive ($\rm M_*> 10^{9.5}M_\odot$) quiescent galaxies from $z=2-7$ in over 800 arcmin$^2$ of NIRCam imaging from a compilation of public JWST fields (with a total area $>$ 5 $\times$ previous JWST studies). We compute and report their cosmic number densities, stellar mass functions, and cosmic stellar mass density. We confirm a significant overabundance of massive quiescent galaxies relative to a range of cosmological hydrodynamical simulations and semi-analytic models (SAMs). We find that no simulations or SAMs accurately reproduce the SMF for massive quiescent galaxies at any redshift within the interval $z=2-5$. This shows that none of these models' feedback prescriptions are fully capturing high-z galaxy quenching, challenging the standard formation scenarios. We find a greater abundance of lower-mass ($\rm M_*<10^{10}M_\odot$) quiescent galaxies than previously found, highlighting the importance of sSFR cuts rather than simple colour selection. We show the importance of this selection bias, alongside individual field-to-field variations caused by cosmic variance, in varying the observed quiescent galaxy SMF, especially at higher-z. We also find a steeper increase in the cosmic stellar mass density for massive quiescent galaxies than has been seen previously, with $\rho_*\propto (1+z)^{-7.2\pm0.3}$, indicating the dramatic increase in the importance of galaxy quenching within these epochs.