The evolution of the galaxy stellar mass function over the last twelve billion years from a combination of ground-based and HST surveys

We present a new determination of the galaxy stellar mass function (GSMF) over the redshift interval $0.25 \leq z \leq 3.75$, derived from a combination of ground-based and Hubble Space Telescope (HST) imaging surveys. Based on a near-IR selected galaxy sample selected over an area of 3 deg$^{2}$ and spanning $\geq 4$ dex in stellar mass, we fit the GSMF with both single and double Schechter functions, carefully accounting for Eddington bias to derive both observed and intrinsic parameter values. We find that a double Schechter function is a better fit to the GSMF at all redshifts, although the single and double Schechter function fits are statistically indistinguishable by $z=3.25$. We find no evidence for significant evolution in $M^{\star}$, with the intrinsic value consistent with $\log_{10}(M^{\star} / M_{\odot})=10.55\pm{0.1}$ over the full redshift range. Overall, our determination of the GSMF is in good agreement with recent simulation results, although differences persist at the highest stellar masses. Splitting our sample according to location on the UVJ plane, we find that the star-forming GSMF can be adequately described by a single Schechter function over the full redshift range, and has not evolved significantly since $z\simeq 2.5$. In contrast, both the normalization and functional form of the passive GSMF evolves dramatically with redshift, switching from a single to a double Schechter function at $z \leq 1.5$. As a result, we find that while passive galaxies dominate the integrated stellar-mass density at $z \leq 0.75$, they only contribute $\lesssim 10$ per cent by $z\simeq 3$. Finally, we provide a simple parameterization that provides an accurate estimate of the GSMF, both observed and intrinsic, at any redshift within the range $0 \leq z \leq 4$.


INTRODUCTION
An accurate determination of the evolving galaxy stellarmass function (GSMF) is crucial for improving our understanding of galaxy evolution. In addition to tracing the history of stellar-mass assembly, the evolving shape of the GSMF encodes vital information about the impact of different feedback mechanisms and the physical processes through which star formation is quenched. As a consequence, together with the cosmic star-formation rate density, the evolving GSMF is arguably one of the most fundamental observational constraints that all theoretical models of galaxy evolution must be able to reproduce.
Over the last two decades, an enormous amount of ef-Email: mcleod@roe.ac.uk fort has been invested exploring the evolution of the GSMF. At low redshifts, numerous studies have exploited the large areas and spectroscopic redshifts provided by the 2dF-GRS (Colless et al. 2001), SDSS (York et al. 2000) and GAMA surveys (Driver et al. 2011) to study the form of the local GSMF (e.g. Cole et al. 2001;Bell et al. 2003;Blanton et al. 2003;Li & White 2009;Baldry et al. 2008Baldry et al. , 2012Weigel et al. 2016). At intermediate redshifts, studies have exploited a combination of photometric and spectroscopic data to study the evolution of the GSMF out to z 1 (e.g. Drory et al. 2009;Pozzetti et al. 2010;Moustakas et al. 2013), while others have used a combination of increasingly deep ground-based and HST near-IR imaging to push the study of the GSMF to z 4 − 5 and beyond (e.g. Fontana et al. 2006;Ilbert et al. 2009Ilbert et al. , 2013Muzzin et al. 2013;Tomczak et al. 2014;Mortlock et al. 2015;Davidzon et al. 2017;Wright et al. 2018;Leja et al. 2020). At higher redshifts still, attempts have been made to constrain the GSMF using the deepest available HST imaging over the redshift range 5 < z < 8 (e.g. Duncan et al. 2014;Grazian et al. 2015;Song et al. 2016;Bhatawdekar et al. 2019;Kikuchihara et al. 2020).
Based on the wealth of literature studies, several characteristics of the GSMF have been firmly established. Firstly, it is clear that the local GSMF is well described by a double Schechter function (Schechter 1976), with a characteristic mass of log 10 (M / M ) 10.6, a low-mass slope of α 2 −1.4 and a high-mass slope of α 1 α 2 + 1.0. Secondly, when the local GSMF is split into star-forming and passive galaxy subsamples, it is clear that the passive GSMF requires a double Schechter function, whereas the star-forming GSMF is usually found to be adequately described by a single Schechter function (e.g. Baldry et al. 2012). Moreover, the majority of previous studies have concluded that the evolution of the star-forming GSMF is remarkably modest, at least out to z 2 (e.g. Tomczak et al. 2014;Davidzon et al. 2017).
A useful insight into the physical information that can be extracted from the GSMF is provided by the analytic model proposed by Peng et al. (2010), which was motivated by the observed stability of the star-forming GSMF and evidence that the effects of mass and environmental quenching appear to be fully separable in the local Universe (e.g. Baldry et al. 2006). In the Peng et al. (2010) model, the exponential cut-off and M of the star-forming GSMF is established and maintained by a mass quenching rate proportional to star-formation rate (SFR). If the slope of the main sequence of star-formation is close to unity (i.e. SFR ∝ M ) then a natural consequence is the build-up of the high-mass component of the passive GSMF, with the same value of M and a low-mass slope of α 1 α + 1.0, where α is the lowmass slope of the star-forming GSMF. In this model, the quenching of high-mass (log 10 (M / M ) ≥ 10.5) galaxies is dominated by mass quenching, usually attributed to some form of active galactic nuclei (AGN) feedback, at all epochs and in all environments.
However, at lower stellar masses, environmental quenching, a combination of galaxy mergers and satellite quenching, becomes increasingly important and dominates at late times (i.e. z < 1). Crucially, because environmental quenching is independent of stellar mass, it naturally produces a second passive-galaxy Schechter function component whose shape, but not normalization, mirrors that of the star-forming GSMF. This apparently simple model can accurately reproduce the key characteristics of the low-redshift GSMF, and illustrates how accurately determining the evolution of the GSMF offers the prospect of constraining the relative timing and importance of different quenching mechanisms.
How well the Peng et al. (2010) model performs at higher redshifts is not entirely clear and the observational constraints are inevitably somewhat less stringent. At z ≤ 1 there is a general consensus that the total GSMF maintains a double Schechter functional form and that the star-forming GSMF remains approximately constant (e.g. Pozzetti et al. 2010;Ilbert et al. 2013;Mortlock et al. 2015). However, at higher redshifts, studies arrive at different conclusions regarding the shape and evolution of the total GSMF and, in particular, the detectability, or otherwise, of an environmentally induced upturn in the number densities of low-mass passive galaxies at z ≥ 1 (e.g. Tomczak et al. 2014;Davidzon et al. 2017;Wright et al. 2018).
Within this context, the primary motivation for this study is to use a combination of the best available ground and space-based photometry, covering a sufficiently large cosmological volume and dynamic range in stellar mass, to accurately determine both the high and low-mass shape of the GSMF out to z 4. To achieve this, we exploit the best available near-IR ground-based imaging over an area of 3 deg 2 and combine it with the publicly available data over the five separate HST CANDELS survey fields (Grogin et al. 2011;Koekemoer et al. 2011). Crucially, in addition to the deepest available optical and near-IR data, the survey fields used in this study also feature the deep mid-IR data from the Spitzer Space Telescope that is necessary to derive robust stellar masses at z ≥ 1.
The ground-based data alone allows us to accurately determine the high-mass end of the GSMF, by accessing a consistent co-moving cosmological volume of 10 7 Mpc 3 in six redshift bins, spanning the range 0.25 ≤ z ≤ 3.75. However, the addition of the HST imaging ensures that we have access to sufficient dynamic range in stellar mass, 2.5 − 3.0 dex below M at all redshifts, to also accurately determine the evolution of the low-mass end of the GSMF.
The structure of the paper is as follows. In Section 2 we discuss the suite of ground-based imaging data utilised in this study, along with the publicly available HST CAN-DELS catalogues. In Section 3 we describe the production of the photometric catalogues, the determination of the photometric redshifts and the construction of the final galaxy sample. In Section 4 we present our determination of the evolving GSMF over the redshift interval 0.25 ≤ z ≤ 3.75 and compare our new results to those of previous studies in the literature and the predictions of the latest theoretical models. Based on our results, we provide a simple evolving parameterization that can accurately reproduce the total GSMF at any redshift within the range 0 ≤ z ≤ 4. In Section 5 we explore the evolution of the star-forming and passive GSMFs and compare to the predictions of the Peng et al. (2010) model. In Section 6 we investigate the evolution of the integrated stellar-mass density and compare with previous literature results, theoretical models and the integral of the cosmic star-formation rate. Finally, we present a summary of our results and conclusions in Section 7. All magnitudes are expressed in the AB system (Oke 1974;Oke & Gunn 1983) and we assume the following cosmology: Ω 0 = 0.3, Ω Λ = 0.7 and H 0 = 70 kms −1 Mpc −1 .

Imaging data
The imaging data utilized in this study primarily consist of ground-based UV+optical+near-IR imaging of the UKIDSS Ultra Deep Survey (UDS), COSMOS and CFHTLS-D1 survey fields. In addition to the ground-based imaging data, we have also made extensive use of the deep Spitzer Space Telescope mid-IR imaging available in all three fields. In Table  1 we list the data used in each field, along with our determinations of the median global 5σ-depths in each filter. The depths have been calculated within a circular aperture Table 1. The median global 5σ-depths for each of the filters used in this study. For all UV, optical and near-IR filters the 5σ-depths have been calculated using a circular aperture with a 2-arcsec diameter and corrected to total assuming a point-source aperture correction. The median depths in the two IRAC bands were calculated using the photometric uncertainties produced by the tphot deconfusion software package (Merlin et al. 2015). Note that the IRAC mosaics used in the COSMOS and UDS fields consist of data from a number of different observing programmes, leading to significant spatial variations in the depth.  with a 2-arcsec diameter and have been corrected to total assuming a point-source aperture correction.

The UKIDSS Ultra Deep Survey
In this study we utilized the JHK near-IR imaging from the latest data release (DR11) of the UKIDSS Ultra Deep Survey (Lawrence et al. 2007;Almaini et al. in prep). Additional Y −band near-IR imaging data was taken from the DR4 release of the VISTA VIDEO survey (Jarvis et al. 2013). The UV and optical coverage of the UDS field consists of CFHT MegaCam u * −band imaging and Subaru Suprime-Cam imaging in the BV Riz and NB921 filters (Furusawa et al. 2008;Koyama et al. 2011;Sobral et al. 2016). Additional z−band imaging (z new ), taken following the refurbishment of Suprime-Cam with CCDs with improved red sensitivity, was also employed (Furusawa et al. 2016). All of the UV, optical and near-IR imaging data in the UDS field was PSF-homogenized to a Moffat profile with a FWHM of 0.9 .
At mid-IR wavelengths, Spitzer IRAC mosaics of the UDS field at 3.6 and 4.5 µm were constructed by combining the data from the SPLASH (PI Capak; see e.g. Mehta et al. 2018), SEDS (Ashby et al. 2013) and S-CANDELS (Ashby et al. 2015) programmes using mopex (Makovoz & Marleau 2005). The overlap region covered by the full set of UV-tomid-IR data in UDS is 0.8 deg 2 .

UltraVISTA
The UltraVISTA survey (McCracken et al. 2012) provides near-IR Y JHK s imaging over an area of 1.5 deg 2 within the COSMOS (Scoville et al. 2007) survey field. The data utilized in the study is comprised of the 1 deg 2 overlap region between the latest UltraVISTA data release (DR4) and the optical u * griz imaging of the CFHTLS-D2 field provided by the T0007 data release of the CFHT Legacy Survey (Hudelot et al. 2012). In addition to the CFHTLS z− band imaging, we also employed deeper Subaru Suprime-Cam z new −band imaging (Furusawa et al. 2016).
While the CFHTLS-D2 and z new imaging is homogeneous, the UltraVISTA imaging is divided into "deep" and "ultra-deep" stripes, that account for approximately 45% and 55% of the total area, respectively. The Y JH imaging in the ultra-deep stripes is typically 0.7 mag deeper than in the deep stripes, whereas the K s −band imaging, thanks to an on-going homogenization programme, is only 0.2 magnitudes deeper (see Table 1). The UV, optical and near-IR imaging data in the UltraVISTA field was PSF-homogenized to a Moffat profile with a FWHM of 1.0 .
As in the UDS field, Spitzer IRAC mosaics at 3.6 and 4.5 µm were constructed by combining data from the SPLASH, SEDS and S-CANDELS surveys, in addition to data from the SMUVS (Ashby et al. 2018) and S-COSMOS (Sanders et al. 2007) programmes.

CFHTLS-D1
In the CFHTLS-D1 survey field we utilized the 1 deg 2 overlap region between the u * griz imaging from the CFHT Legacy Survey (Hudelot et al. 2012) and the Y JHK s near-IR imaging from the VISTA VIDEO survey (Jarvis et al. 2013). The UV, optical and near-IR imaging in the CFHTLS-D1 field was PSF-homogenized to a Moffat profile with a FWHM of 1.0 . The mid-IR Spitzer imaging at 3.6 and 4.5 µm was provided by the SERVS programme (Mauduit et al. 2012).

CANDELS catalogues
In order to increase the available dynamic range in stellar mass, we have used the publicly available photometric redshifts and stellar masses derived for each of the five HST CANDELS fields (Guo et al. 2013;Galametz et al. 2013;Santini et al. 2015;Stefanon et al. 2017;Nayyeri et al. 2017;Barro et al. 2019). Although the CANDELS catalogues only cover a total area of 0.27 deg 2 , the depth of the HST near-IR imaging plays a crucial role in our ability to properly constrain the low-mass end of the GSMF. A brief description of the relevant properties of each CANDELS catalogue is provided in Table 2.

CATALOGUE PRODUCTION AND SAMPLE SELECTION
In this section we provide an overview of how the photometry catalogues for the UKIDSS UDS, UltraVISTA and CFHTLS-D1 survey fields were produced. We also provide an overview of how the photometric redshifts and stellar  (Grogin et al. 2011;Koekemoer et al. 2011) and, in the case of GOODS-S, the ultra-deep near-IR data available in the Hubble Ultra Deep Field (HUDF). The final column lists the references for the catalogues containing the photometry, photometric redshifts and stellar mass information; which correspond to (1) Galametz et al. (2013), (2) Guo et al. (2013), (3) Santini et al. (2015), (4) Nayyeri et al. (2017), (5) Stefanon et al. (2017) and (6)  masses were calculated and the processes employed to select the final galaxy sample.

Photometry catalogues
In order to generate photometric catalogues, sextractor (Bertin & Arnouts 1996) was run in dual-image mode with the K−band image serving as the detection image across all three fields. For objects detected in the K−band, photometry was extracted from the PSF-homogenized images in all UV-to-near-IR filters using circular apertures with a 2-arcsec diameter. The photometry was extracted from the PSF-homogenized images to minimize aperture correction effects for the subsequent spectral energy distribution (SED) fitting. To further reduce any colour systematics, additional flux corrections were made (at the 1 − 2% level) based on the curves of growth of point sources in each filter. Accurate flux errors were calculated for each object, in each individual filter, by measuring the aperture-to-aperture r.m.s. of 150 − 200 nearby blank sky apertures (see e.g. McLeod et al. 2015), where the local value of σ is calculated using the robust median absolute deviation (MAD) estimator.
The photometry in the lower spatial resolution Spitzer IRAC imaging at 3.6 and 4.5 µm was measured using the tphot software package (Merlin et al. 2015). Given that the tphot algorithm uses the isophotal footprint of the objects in a higher spatial-resolution image as a prior (the original K−band detection image in this case), the fluxes generated by tphot can be regarded as isophotal. As a consequence, we aperture match the PSF-homogenized photometry at shorter wavelengths to match the tphot fluxes, multiplying by f = K iso /K 2 , where K iso is the isophotal flux extracted by sextractor from the high-resolution image, and K 2 is the 2-arcsec diameter flux extracted from the PSF-homogenized K−band image.

Photometric redshifts
In order to derive robust photometric redshifts, we undertook six different photometric redshift runs, using three different codes. Three photometric redshift runs were performed using the LePhare (Arnouts & Ilbert 2011) SED fitting code, using the Bruzual & Charlot (2003), Pegase2 (Fioc & Rocca-Volmerange 1999) and COSMOS (Ilbert et al. 2009) template libraries. In each of these runs a Calzetti et al. (2000) dust attenuation curve was adopted, with colour excess in the range E(B − V) = 0 − 0.6. Emission lines were included in the fits and IGM absorption was accounted for using the Madau (1995) prescription.
Two further photometric redshift runs were performed using the eazy SED fitting code (Brammer et al. 2008), with the default PCA and Pegase2 libraries. One final run was performed with the BPZ code (Benítez 2000), using the default set-up and CWW (Coleman et al. 1980) templates.
Before running on the full photometry catalogues, the three different SED fitting codes were trained by fitting to the photometry of objects with robust spectroscopic redshifts. This process allowed us to apply the necessary zeropoint off-sets (e.g. Dahlen et al. 2013) and to quantify the performance of each code/template combination using σ z and the catastrophic outlier rate. The value of σ z is our preferred measurement of the photometric redshift accuracy, and is defined as 1.483 × MAD(dz), where MAD is the median absolute deviation and dz = (z phot −z spec )/(1+z spec ). Any object with |dz| > 0.15 is classified as a catastrophic outlier.
Our best-estimate z phot for each object was taken as the median of our six different photometric redshift estimates (hereafter z med ). These z med measurements were tested against spectroscopic redshift samples for each of the three fields to ensure their accuracy.
For the UltraVISTA/COSMOS field we compiled a catalogue of 11000 high-quality spectroscopic redshifts, the vast majority of which were drawn from the zCOSMOS (Lilly et al. 2007(Lilly et al. , 2009 For this sub-sample, the photometric redshift accuracy was σ z = 0.019, with a catastrophic outlier rate of 2.5%. Finally, to test the photometric redshifts in the CFHTLS-D1 field, we used a sample of 4200 robust spectroscopic redshifts from the VIMOS VLT Deep Survey (VVDS, Le Fèvre et al. 2013). The photometric redshift accuracy for this sub-sample was σ z = 0.019, with a catastrophic outlier rate of 2.3%.
In Fig. 1 we show a comparison between the spectroscopic and photometric redshifts for our training set of 18000 objects across all three survey fields. In summary, it can be seen that our z med photometric redshifts are both robust and consistent across all three fields, with a typical accuracy of σ z = 0.021 and a catastrophic outlier rate of 2.4%. We note that the publicly available photometric red-  Figure 1. A comparison between photometric and spectroscopic redshifts for 18000 galaxies with robust spectroscopic redshifts across the three ground-based data sets. The corresponding values of σ MAD and the fraction of catastrophic outliers are provided in the legend. A density map has been used to aid legibility.
shifts available for the five CANDELS fields (see Table 2) are of very comparable quality to those derived here.
Once the photometric redshift training process had been completed, the six different photometric redshift code/template combinations were run on the full photometric catalogues for the UDS, UltraVISTA and CFHTLS-D1 fields. As before, the final adopted photometric redshift for each object was taken to be z med , the median of our six different estimates.

Sample construction
In order to construct the final sample to be used in the GSMF determination, the ground-based catalogues were initially cut at their global 5σ limit in the K−band and then restricted to those objects in the photometric redshift range 0.25 ≤ z phot ≤ 3.75. The lower redshift cut is imposed to ensure that our survey encloses sufficient cosmological volume to constrain the GSMF, whereas the upper redshift cut is imposed to ensure that our K−band selection always corresponds to rest-frame wavelengths long-ward of the 4000Å break. Following the initial cuts based on the near-IR signalto-noise ratio and photometric redshift, the ground-based catalogues were then cleaned to remove stars, AGN, artefacts and objects with contaminated photometry.
In order to reduce the contamination by stars, objects were excluded from the final catalogue using the stellar locus in two colour-colour plots. For the UDS sub-sample, stars were rejected using the stellar locus on the B − z vs z − K colour-colour plot, following Baldry et al. (2010). For the UltraVISTA and CFHTLS -D1 sub-samples, stars were rejected using the stellar locus on the g − i vs J − K s colour plot, following Jarvis et al. (2013).
Potential AGN were removed from the final sample based on a combination of X-ray and Spitzer 24 µm information. The initial rejection of potential AGN was performed using the publicly available X-ray catalogues which cover all three ground-based fields. In the COSMOS field we utilized the Chandra X-ray catalogue published by Elvis et al. (2009), while in the UDS field we used the XMM-Newton Xray catalogue produced by the Subaru/XMM-Newton Deep Survey (Ueda et al. 2008). Finally, in the CFHTLS-D1 field we utilized the XMM-XXL north survey catalogue from Liu et al. (2016). Based on the information available in these catalogues, potential AGN were excluded using the soft X-ray to optical ratio (X/O), as described in Salvato et al. (2011).
Following the exclusion of potential AGN based on their X-ray characteristics, we performed a second round of AGN rejection based on the 24 µm imaging available in the UDS (spUDS, PID 40021, PI Dunlop), COSMOS (S-COSMOS, Sanders et al. 2007) and CFHTLS-D1 (SWIRE, Lonsdale et al. 2003) fields. Objects were removed based on the specific star-formation rate criteria: sSFR≥ 10.0 Gyr −1 , where the specific star-formation rate was calculated using the 24 µm prescription of Rieke et al. (2009), which was designed to remove those objects whose 24 µm flux is dominated by AGN heated dust. The fraction of objects removed from our sample as potential AGN is of order 1%.
Artefacts and objects with contaminated photometry were removed using a two-step process. Firstly, all objects within the haloes of bright/saturated stars were removed from the catalogue, with the effective survey area recalculated to compensate. Secondly, within each of the GSMF redshift bins, those objects whose SED fits produced the worst 5% of χ 2 values were also excluded. The SED fits for this population of objects were statistically unacceptable, and visual inspection confirmed that they were dominated by artefacts and objects with badly compromised photometry.
Finally, those objects from the five CANDELS catalogues detected at ≥ 5σ significance in the H 160 filter and within the redshift range 0.25 ≤ z ≤ 3.75 were included in our GSMF sample. Objects identified as stars or AGN were once again excluded, based on the flags provided.

Stellar masses
The stellar masses were calculated by re-fitting the photometry in our UKIDSS UDS, UltraVISTA and CFHTLS-D1 catalogues using a standard set of Bruzual & Charlot (2003) templates, based on a Chabrier (2003) initial mass function (IMF), τ−model star-formation histories and the same dust and IGM absorption prescriptions as described in Section 3.2. During this SED fitting process the redshift for each galaxy was held fixed at z med .
At this stage, the stellar masses returned by the SED fitting process were based on isophotal photometry (see Section 3.1). To convert the stellar masses to total, they were multiplied by K auto /0.9K iso , where K iso and K auto are the flux iso and flux auto fluxes measured by sextractor in the K−band, respectively. The factor of 0.9 is necessary to account for the fact that flux auto typically only captures 90% of the total flux 1 . The left-hand panel shows the 90% mass-completeness limits as a function of redshift for each of the three ground-based surveys fields, where the UltraVISTA/COSMOS field has been separated into the deep and ultra-deep components. The grey vertical dashed lines show the limits of the six redshift bins adopted for the determination of the GSMF. The right-hand panel shows the same information for the CANDELS fields, where the GOODS-South field has been separated into the wide, deep and HUDF components. Note that, for clarity, we do not show the mass-completeness limits for the wide and deep components of the GOODS-North field, which are assumed to be the same as for the equivalent components of GOODS-South.
The stellar masses for the objects within the five CAN-DELS fields were taken from the public catalogues (see Table 2). The stellar masses across all five CANDELS fields were also calculated using BC03 stellar population templates, based on a Chabrier IMF. Cross checks performed using objects in common with our ground-based photometry catalogues in the UDS and COSMOS fields confirmed that our stellar-mass measurements are in excellent agreement with those derived by the CANDELS team, following a tight 1:1 relation with a typical scatter of ±0.05 dex.

GSMF DETERMINATION
In this section we present our basic determination of the evolving GSMF, taking into account the effects of stellarmass completeness and an assessment of the impact of cosmic variance. We also provide a full description of how we fit the observed and intrinsic GSMF, after accounting for the effects of Eddington bias.

Number densities
After first splitting the data into six redshift bins, we employed the 1 V max estimator (Schmidt 1968) to determine the number densities: where is the number density of galaxies per dex, per unit comoving volume, ∆M is the logarithmic stellar mass bin and C i (M, z) is the completeness calculated for each galaxy. In Eqn. 1 we have defined M ≡ log 10 (M / M ) and will repeatedly adopt this shorthand throughout the rest of the paper. Incompleteness in each of the three ground-based surveys was accounted for separately via simulations, in which artificial galaxies with a wide range of physical properties (z, M K , M, r e ) were injected into the K−band imaging and recovered using the same sextractor set-up used to construct the original photometry catalogues. In performing the simulations we converted between M K and stellar mass by drawing randomly from the SED templates fitted to the real galaxies and adopted half-light radii predicted by the sizemass-redshift distributions derived by Shibuya et al. (2015). To calculate the effective stellar-mass limit of our groundbased survey data, we followed the procedure proposed by Pozzetti et al. (2010) and calculated the distribution of limiting masses for the galaxies at each redshift, where the limited mass is defined as: and (K −K lim ) is the difference between the apparent K−band magnitude of a galaxy and the 5σ magnitude limit. Using this method we define the 90% mass-completeness limit at each redshift as the stellar mass below which 90% of the limiting stellar masses lie. Unlike Pozzetti et al. (2010), at each redshift we calculate a more conservative limiting-mass based on the full galaxy sample (i.e. using the full range of mass-to-light ratios), rather than adopting the faintest 20% of galaxies. In Fig. 2 we plot the 90% mass-completeness limit versus redshift for both the ground-based and CAN-DELS data. All of the GSMF plots shown in this paper include only those galaxies that lie above the appropriate 90% mass-completeness limit.

Individual ground-based GSMF determinations
In Fig. 3 we show our determinations of the GSMF, based on the ground-based data in the UDS, COSMOS and CFHTLS-D1 and fields alone. The first five redshift bins all have a width of ∆z = 0.5 and are centred on z = 0.5, 1.0, 1.5, 2.0 and 2.5, whereas the final redshift bin spans the range 2.75 ≤ z < 3.75 in order to maintain the statistics at a similar level to the lower-redshift bins. It can be seen that the independent GSMF determinations for the three, degree-scale, fields are generally in good agreement, although significant differences are present, most noticeably at the high-mass end. For example, at z = 0.5−1.0 8.5 9.0 9.5 10.0 10.5 11.0 11.5 12.0 12.5 the COSMOS field can be seen to be somewhat overdense compared to the other two fields, a fact which was also noted by Moustakas et al. (2013). In terms of cosmic variance, the advantage of determining the high-mass end of the GSMF from three independent degree-scale fields is therefore clear. Indeed, a further advantage of having three non-contiguous fields is that we are able to empirically quantify the level of cosmic variance in the GSMF. The number density uncertainties plotted in Fig. 3 are simply the Poissonian counting errors. However, in all subsequent plots and tables, the quoted uncertainties are based on the quadrature addition of σ poisson , σ boot and σ cv , where σ boot is the error contribution calculated from a bootstrap analysis based on many thousands of GSMF realisations. For σ cv , we consider both an empirical estimate based on the field-to-field variance and that estimated using Moster et al. (2011). The cosmic variance contribution (σ cv ) to each bin is taken as the greater of the Moster et al. (2011) estimate and the measured field-to-field variance of those fields contributing to the bin. The number density uncertainties for the HST CANDELS data were calculated in an identical fashion.

The combined HST and ground-based GSMF
In Table 3 and Fig. 4 we present our determination of the observed GSMF over the redshift range 0.25 ≤ z < 3.75, based on a combination of the full ground-based and HST data set. The process adopted for producing the combined GSMF determination was as follows. Firstly, we produced a combined ground-based GSMF by merging the three ground-based catalogues into a single catalogue, calculating the numbers of objects and the cosmological volume contributing to each redshift-mass bin based on the 90% mass-completeness limits for each field. Secondly, we produced a combined HSTbased GSMF by applying an identical methodology to the data for the five CANDELS fields. The third step in the process was to match the ground-based and HST-based GSMFs in each redshift bin, but adjusting the normalization of the HST-based GSMF to match that of the ground-based GSMF in the overlap region between the two. The typical adjustment required was at the ±0.03 dex level. In each redshift bin, the final split between the ground-based (black data points) and HST-based (blue data points) shown in Fig. 4 is based on the 90% completeness limit of the deepest groundbased survey field (typically the UDS). For clarity, the blue data points shown in Fig. 4 are entirely based on HST data, whereas the black points are entirely based on ground-based data.

Fitting the observed GSMF
In each redshift range we derive maximum likelihood fits to the binned GSMF data shown in Fig. 4 using both a single and double Schechter function parameterization (Schechter 1976). The single Schechter function has the following functional form: where φ(M) is the number density of galaxies per Mpc 3 per dex stellar mass, M ≡ log(M / M ), where M is the characteristic stellar mass and α is the low-mass slope. The double Schechter function has the form: where both components have the same characteristic mass and we define α 1 to be the high-mass slope and α 2 to be the low-mass slope. The best-fitting parameters and their corresponding uncertainties are presented in Table 4, which also Table 3. The observed GSMF as a function of redshift, based on the combined ground-based and HST data set. The first column lists the adopted stellar mass bins, where M ≡ log 10 (M / M ), while columns 2 − 7 list the logarithm of the number densities (φ k ) within six redshift bins. The units of φ k are dex −1 Mpc −3 . The data presented in this table are plotted in Fig. 4.  Table 4. The best-fitting single (upper section) and double (lower section) Schechter function parameters to the observed GSMF, where M ≡ log 10 (M / M ) and the units of φ , φ 1 and φ 2 are dex −1 Mpc −3 . The final two columns list the χ 2 and χ 2 ν values of the fits, respectively. The best-fitting double Schechter function parameters for the local observed GSMF derived by Baldry et al. (2012) are provided for comparison.   . The observed GSMF as a function of redshift, based on the combined ground-based and HST data set. Also plotted are the best-fitting single (dashed grey) and double (solid black) Schechter function fits. The black data points are derived from the ground-based data set alone, while the blue data points are derived from the HST data set alone. The split between the ground-based and HST data is highlighted by the dashed grey vertical line. Over the redshift range 0.25 ≤ z ≤ 2.75, the double Schechter fit is seen to be a better representation of the observed GSMF than the single Schechter fit. However, in the final redshift bin centred on z = 3.25, the single and double Schechter function fits are basically indistinguishable.
includes, for comparison, the double Schechter function parameters for the local GSMF derived by Baldry et al. (2012). The best-fitting single and double Schechter functions are plotted as the dashed grey curves and solid black curves in Fig. 4, respectively. It can be seen from Fig. 4 that a double Schechter function appears to provide a better description of the data in all redshift bins out to z = 2.5, whereas in the final redshift bin at z = 3.25 the difference between the single and double Schechter function fit is negligible. This impression is confirmed by the information displayed in Table 4, which shows that the double Schechter function provides a better statistical description of the data in all six redshifts bins. Notably, in the first four redshift bins, covering the redshift range 0.25 ≤ z < 2.25, the single Schechter function does not provide a statistically acceptable fit to the data, whereas the double Schechter function fit is statistically acceptable at all redshifts. However, it is also worth noting that in the last two redshift bins, covering the range 2.25 ≤ z < 3.75, the single and double Schechter function fits are both statistically acceptable.
Overall, it can be seen from Table 4 that the best-fitting Schechter function parameters exhibit remarkably smooth and modest evolution over the redshift range studied here, a subject we will return to in Section 4.6. Focusing on the double Schechter function parameters, it can be seen that the characteristic stellar mass, in particular, remains remarkably constant, lying within the range M = 10.75 ± 0.1 at all redshifts. Moreover, it is noteworthy that the difference in the fitted Schechter function power-law indices (i.e. ∆α = α 1 − α 2 ) is consistent with unity at all redshifts, with a variance weighted mean difference of ∆α = 1.09 ± 0.21. This result is in good agreement with the phenomenological model of Peng et al. (2010), and is a subject we will return to when we investigate the individual star-forming and passive GSMFs in Section 5. The evolution of the observed GSMF can be seen more clearly in Fig. 5, which shows an overlay of both the data and the corresponding best-fitting double Schechter function fits over the full 0.25 ≤ z < 3.75 redshift range. To allow comparison with the local GSMF, we have also included the double Schechter function fit from Baldry et al. (2012). This plot very clearly illustrates that there is very little evolution in the observed GSMF at either the low-mass (8.0 < M < 9.0) or high-mass (M > 11.5) end. In contrast, substantial evolution is apparent in the number density of galaxies close to the characteristic stellar mass (i.e. M 10.75).

Eddington bias
The observed GSMF results presented in the previous section will inevitably be subject to Eddington bias, whereby the combination of stellar-mass uncertainties and the steep exponential fall-off of the GSMF leads to a net bias towards higher number densities at the high-mass end of the GSMF. As a consequence, the Schechter function parameters derived from a direct fit to the observed GSMF data points will be biased with respect to the intrinsic Schechter function parameters. Recovering the intrinsic form of the GSMF is of particular interest for direct comparison to galaxy simulation results (see Section 4.5) and for deriving accurate measurements of the integrated stellar-mass density (see Section 6). In order to recover the intrinsic form of the GSMF it is necessary to determine the effective uncertainties in the stellar-mass measurements, which we model as a log-normal distribution with σ M . We adopted two approaches to quantifying σ M . The first approach was to run a series of simulations in which the photometry for each galaxy was scattered according to its errors, before the photometric redshifts and stellar masses were re-calculated. The results of these simulations indicated that σ M 0.2 dex and was not a strong function of either redshift or stellar mass.
The second approach was to use the binned GSMF data itself to constrain the value of σ M . This process involved refitting the observed GSMF data with a double Schechter function as before, but convolving the intrinsic Schechter function with a log-normal distribution with σ M , where σ M is a free parameter in the fit. The results of this fitting process demonstrated that σ M had a mean value of 0.15 dex and lay within the range 0.15 ± 0.04 dex in all six redshift bins. As a result, we adopted a value of σ M = 0.15 dex and re-ran fits to the observed GSMF data including the convolution due to stellar-mass uncertainties. The best-fitting Schechter function parameters from this fitting process represent our best estimates of the intrinsic form of the GSMF. The typical impact of Eddington bias on the form of the best-fitting Schechter function is illustrated in Fig. 6.

The intrinsic GSMF
Our determination of the best-fitting intrinsic Schechter function parameters is presented in Table 5. As with the best-fitting observed Schechter function parameters shown in Table 4, it can be seen that the double Schechter function provides a better description of the data at all redshifts, and that the single Schechter function fits are statistically unacceptable at z < 2.25. Once again, the single and dou- Table 5. The best-fitting single (upper section) and double (lower section) Schechter function parameters to the intrinsic GSMF, where M ≡ log 10 (M / M ) and the units of φ , φ 1 and φ 2 are dex −1 Mpc −3 . The final two columns list the χ 2 and χ 2 ν values of the fits, respectively. As described in the text, the best-fitting intrinsic Schechter function parameters have been derived by incorporating a convolution of σ M = 0.15 dex when fitting the observed GSMF data. For comparison, we also include our estimate of the intrinsic local GSMF, derived by fitting the data from Baldry et al. (2012) assuming that σ M = 0.1 dex (Wright et al. 2018). The evolution of the best-fitting observed and intrinsic double Schechter function parameters is shown in Fig.  7. As is to be expected, following what is essentially a deconvolution process, the normalizations of the intrinsic double Schechter function shift to slightly higher values, and the slopes shift to slightly shallower values. Unsurprisingly, the largest difference between the observed and intrinsic Schechter function parameters is the best-fitting value of the characteristic stellar mass. In both cases the characteristic stellar mass displays remarkably little redshift evolution but, after accounting for the effect of Eddington bias, the bestfitting intrinsic value of M is shifted 0.2 dex lower, lying within the range M 10.55 ± 0.1 at all redshifts. In Fig. 8 we show a comparison between our intrinsic double Schechter function parameters and those derived by recent studies in the literature (Ilbert et al. 2013;Tomczak et al. 2014;Davidzon et al. 2017;Wright et al. 2018;Kawinwanichakij et al. 2020). We note that all of the studies we compare to in Fig. 8 quote intrinsic Schechter function parameters, with the exception of Tomczak et al. (2014) who quote observed parameters. It can be seen from Fig. 8 that our intrinsic Schechter function parameters are in reasonable agreement with previous determinations, although the parameter estimates in the literature span a significant range. It is also clear from Fig. 8 that the unique combination of the dynamic range and cosmological volume sampled by this work has led to significantly improved parameter constraints.
It is noteworthy that our determination of the evolving characteristic stellar mass is significantly lower than most previous studies. However, it can be seen from Fig. 8 that our determination of M is in good agreement with Davidzon et al. (2017), and is also likely to be in good agreement with Tomczak et al. (2014), assuming that the Eddington bias correction for their data set is similar to our estimate.

Comparison to simulations
In Fig. 9 we show a comparison between our observed and intrinsic GSMFs and the results of the EAGLE (Furlong et al. 2015), SIMBA (Davé et al. 2019) and Munich galaxy formation models (Henriques et al. 2015).
At z = 0.5 and z = 1.0 the EAGLE hydrodynamical simulation is in generally good agreement with our GSMF determinations, particularly at M ≥ M . However, it can also be seen that EAGLE systematically over-predicts the numbers of M ≤ M galaxies in the redshift range 0.5 ≤ z ≤ 2.0, by a factor of 1.5 − 2.0. By z = 3, it is notable that EAGLE predicts a significantly lower value of M than we observe.
The Munich semi-analytic model is in good agreement with our determination of the M ≤ M number densities at z = 1.0 and z = 3.0, but over-predicts the M ≤ M number densities at z = 2.0 by a factor of 1.5. It is notable that the Munich model systematically over-predicts the number of galaxies with the highest stellar masses (i.e. M ≥ 11.5).
The systematic over-prediction of the high-mass end of the GSMF is a problem that is shared by the SIMBA hydrodynamical simulation. That said, the SIMBA results at z = 1.0 and z = 2.0 are in generally good agreement with our observational results at M ≤ M .
In Section 6, we compare the redshift evolution of the integrated stellar-mass density predicted by the three theoretical models with our observational results. It can be seen that, following what is essentially a deconvolution process, the normalizations of the intrinsic double Schechter function shift to slightly higher values and the slopes shift to slightly shallower values. The largest difference between the observed and intrinsic parameters is M , which shifts 0.2 dex lower after Eddington bias has been accounted for.

An evolving fit to the galaxy stellar mass function
For the purposes of comparing to a variety of different theoretical and observational results, it is clearly desirable to be able to derive an accurate estimate of the total GSMF at any redshift. Guided by the smoothly evolving double Schechter function parameters produced by our maximum likelihood fitting (see Fig. 7), we adopt the following functional forms: log(φ 1 ) = a 7 + a 8 z + a 9 z 2 (8) log(φ 2 ) = a 10 + a 11 z.
It can be seen that all of the parameters follow a simple linear evolution with z, with the exception of log(φ 1 ), which includes an additional z 2 term to account for the steep decline at z > 2.5.
To fit this 11-parameter model to the data we used dynesty (Speagle 2020), a python implementation of the Bayesian Nested Sampling algorithm (Skilling 2006), that can be used to estimate posteriors on model parameters after specifying an appropriate likelihood function and prior. We fitted the model across all 96 (z, log 10 (φ), M) values in our data set (Table 3), assuming a standard Gaussian likelihood of the form: where φ i and σ φ i are the observed number densities and their corresponding errors at a given redshift and stellar mass, and φ i (θ) are the corresponding model number densities, based on the parameter set θ = (a 1 , . . . , a 11 ). We assumed flat priors on all parameters, within the ranges specified in Table 6. We ran fits to both the observed and the intrinsic GSMF. For the intrinsic GSMF we assumed a constant σ M = 0.15 dex based on our maximum likelihood analysis. The median posterior parameter values and their corresponding 68% confidence intervals for both fits are listed in Table 6, and a corner plot showing the 1-D and 2-D marginalized posteriors for the fit to the intrinsic GSMF is shown in Fig. A1.
The best-fitting evolving model is shown in Fig. 10 and provides an excellent fit to the data, with χ 2 ν = 1.00 and 0.85 for the observed and intrinsic fits, respectively. As a result, this simplified evolving prescription can be used to provide an accurate estimate of the GSMF (both observed and intrinsic) at any desired redshift within the range 0.0 ≤ z ≤ 3.75. In fact, comparison with GSMF constraints at higher redshift (e.g. Duncan et al. 2014;Grazian et al. 2015;Song et al. 2016) suggests that the parameterization presented here remains in reasonable agreement with observational constraints out to z 5.  Figure 10. The left-hand panel shows a comparison between the observed GSMF data and the fit produced using our simple evolving parameterization (Eqns. 5 − 9, Table 6) over the full redshift range of this study. The right-hand panel shows our simple parameterization of the evolving intrinsic GSMF, after accounting for Eddington bias (see text for details.) Table 6. Details of our simple parameterization of the evolving GSMF. The first column lists the parameters (see Eqns. 5 − 9) and the second column lists the range of the corresponding flat priors adopted during the fitting process. The final two columns list the best-fitting parameter values from the fit to the observed and intrinsic GSMF, respectively.
Parameter Prior range GSMF GSMF observed intrinsic

THE PASSIVE AND STAR-FORMING GALAXY STELLAR MASS FUNCTIONS
In this section we proceed to split our galaxy sample into its star-forming and passive components, in order to explore how the GSMF and integrated stellar-mass density of each component evolves with redshift. Fundamentally, the differential evolution of the star-forming and passive GSMFs provides crucial constraints on the impact of mass and environmental quenching as a function of redshift and stellar/halo mass.

UVJ selection
We separate our galaxy sample into its star-forming and passive components using the UVJ colour-colour criteria proposed by Williams et al. (2009). Specifically, we follow the results of Carnall et al. (2018Carnall et al. ( , 2020 and apply the following criteria: at all redshifts. Although not entirely model independent, applying this set of criteria is robust, and has the advantage of being easy to apply to a wide variety of observed and simulated data sets. The rest-frame colours for the groundbased component of our final galaxy sample were generated from the SED fitting described in Section 3.4. For the HST CANDELS component, we adopted the rest-frame colours from the relevant publicly available catalogue (see Table 2), and used the overlap with our ground-based data to adjust for off-sets due to different filter definitions (typically at the ±0.1 mag level). The number densities of the passive and star-forming populations were calculated according to the description provided in Section 4.1. However, when calculating the uncertainties associated with the number densities in each (M, z) bin, an additional contribution was included to account for objects scattering into and out of the UVJ selection box due to photometric uncertainties. This additional contribution was calculated as part of the bootstrap simulations described previously, by scattering the U − V and V − J colours according to their uncertainties. The number densities for the star-forming and passive galaxies are presented in Table  7 and Table 8, respectively.

Schechter function fits
We performed maximum likelihood fitting to the starforming and passive galaxy number densities, using both Table 7. The observed GSMF for UVJ-selected star-forming galaxies (see text for details). The first column lists the adopted stellar mass bins, where M ≡ log 10 (M / M ), while columns 2 − 7 list the logarithm of the number densities (φ k ) within six redshift bins. The units of φ k are dex −1 Mpc −3 . The data presented in this table are plotted in Fig. 11. 0.25 ≤ z < 0.75 0.75 ≤ z < 1.25 1.25 ≤ z < 1.75 1.75 ≤ z < 2.25 2.25 ≤ z < 2.75 2.75 ≤ z < 3.75 M log 10 (φ k ) log 10 (φ k ) log 10 (φ k ) log 10 (φ k ) log 10 (φ k ) log 10 (φ k )  M log 10 (φ k ) log 10 (φ k ) log 10 (φ k ) log 10 (φ k ) log 10 (φ k ) log 10 (φ k )  Schechter function blue Schechter function red Figure 11. The redshift evolution of the observed GSMF for star-forming (blue data points) and passive galaxies (red data points). The solid blue curves show the best-fitting single Schechter functions to the star-forming GSMF. The red curves show the best-fitting double Schechter functions to the passive GSMF in the first three redshift bins and the best-fitting single Schechter functions to the passive GSMF in the final three redshift bins. In each redshift bin, the dashed blue and red curves show the best fits to the star-forming and passive GSMFs using the 5-parameter Peng et al. (2010) model (see text for discussion).
single and double Schechter functional forms, as before. The best-fitting Schechter function parameters are presented in Tables 9 & 10 for the star-forming and passive galaxies, respectively. In Fig. 11 we plot the separate star-forming and passive GSMFs, along with the best-fitting single and double Schechter function fits.

The star-forming GSMF
We find that a single Schechter function provides a statistically acceptable description of both the observed and intrinsic star-forming galaxy GSMF, over the full 0.25 ≤ z < 3.75 redshift range covered by our data set. This is in contrast to some other studies (e.g. Drory et al. 2009;Ilbert et al. 2013;Tomczak et al. 2014;Davidzon et al. 2017), who concluded that a double Schechter function is a superior fit to the starforming GSMF, particularly at z ≤ 2. Based on our data, we find that although double Schechter function fits do return lower values of reduced χ 2 , the improvement in the quality of the fit is not generally sufficient to justify the inclusion of two additional degrees of freedom.
It can be seen from Table 9 and Fig. 11 that the starforming GSMF is remarkably stable. Over the redshift range 0.0 ≤ z < 1.25 the intrinsic Schechter function parameters are effectively constant (within the errors). The same statement can be made about the observed Schechter function parameters over this redshift range, with the 2σ shift in M between z 0 and z = 0.5 being largely attributable to increased Eddington bias. The star-forming GSMF only evolves gradually over the redshift interval 1.25 ≤ z < 2.75, with both the observed and intrinsic Schechter function fits displaying a 0.2 dex drop in φ , while the values of M and α remain essentially unchanged. It is only in the highest redshift bin at 2.75 ≤ z < 3.75 that we see a significant change in the shape of the star-forming GSMF, with both the observed and intrinsic Schechter functions showing a 0.2 dex increase in M , a further 0.7 dex drop in φ and a significant 0.3 steepening in the faint-end slope. The remarkably gradual evolution of the star-forming GSMF is illustrated by the left-hand panel of Fig. 12, which shows an overlay of the data and best-fitting Schechter function for all six redshift bins.

The passive GSMF
In contrast to the star-forming GSMF, it can be seen from Table 10 and Fig. 11 that the passive GSMF evolves dramatically over the redshift range studied here. As discussed in the introduction, it has long been established that a double Schechter function is required to match the shape of the passive galaxy GSMF in the local Universe, due to a distinct upturn in the number densities of passive galaxies at low stellar masses, that is usually interpreted as a clear signature of environmental quenching (e.g. Peng et al. 2010).
Our results indicate that the double Schechter functional form of the passive GSMF persists until at least z 1.0, and very likely until z 1.5. If confirmed, the upturn in the number densities of low-mass passive galaxies seen in the z = 1.5 redshift bin would argue that the impact of some form of environmental quenching, presumably galaxy-galaxy mergers rather than satellite quenching, is becoming apparent at a look-back time of ≥ 9 Gyr. This result is in agreement with the previous GSMF study of Tomczak et al. (2014), who also concluded that the passive GSMF required a double Schechter at z ≤ 1.5. Moreover, although Mortlock et al. (2015) only fitted a single Schechter function to the passive GSMF at z ≥ 1, there is an indication of an upturn at low stellar masses in their 1.0 < z < 1.5 redshift bin. In contrast, the recent study by Davidzon et al. (2017) only detects evidence of an upturn in the passive GSMF at z ≤ 0.8, although this is almost certainly explained by a lack of dynamic range in stellar mass.
At redshifts z ≥ 1.75, our determination of the passive GSMF is well described by a single Schechter function. Due to the increase of our stellar-mass completeness limit with redshift, based on the current data set, it is not possible to determine whether this change in shape is intrinsic, or simply due to insufficient dynamic range in stellar mass. Accurately determining the shape of the passive GSMF is clearly a task which can be addressed with the unique near-IR sensitivity offered by the James Webb Space Telescope (JWST).
The dramatic evolution of the passive GSMF is illustrated by the right-hand panel of Fig. 12, which shows an overlay of the data and the best-fitting Schechter function fits over the full redshift range. In addition to the change in shape, this figure also illustrates the dramatic ( 2 dex) decrease in the number density of M M passive galaxies from z 0 to z 3. The evolving contribution of passive galaxies to the integrated stellar-mass density is explored in Section 6.

A comparison with the Peng et al. model
In the introduction we discussed how the empirical model proposed by Peng et al. (2010) can provide useful insights into how different quenching mechanisms control the shape of the GSMF. Given that the Peng et al. model can accurately reproduce the shape of local star-forming and passive GSMFs (e.g. Baldry et al. 2012), it is interesting to explore how well the model continues to perform at higher redshifts.
To investigate this question, we performed maximum Table 9. The best-fitting observed (upper section) and intrinsic (lower section) Schechter function parameters for the star-forming GSMF, where M ≡ log 10 (M / M ) and the units of φ are dex −1 Mpc −3 . The final two columns list the χ 2 and χ 2 ν values of the fits, respectively. We have included the parameters from our own fits to the Baldry et al. (2012) data at z < 0.06 for comparison. To derive the intrinsic Schechter function parameters for the Baldry et al. (2012) data we assumed that σ M = 0.1 dex (Wright et al. 2018 likelihood fits to the star-forming and passive GSMF data in Tables 7 & 8  In the z = 0.5 redshift bin the Peng et al. (2010) model continues to produce an excellent qualitative, and statistically acceptable, match to the star-forming and passive GSMFs. In the next two redshift bins, the Peng et al. model continues to produce a good qualitative match to the observed data, although the fits become progressively poorer in a statistical sense. In the final three redshift bins, the Peng et al. model struggles to match the shape of the passive GSMF, although it arguably still produces a respectable qualitative match to the observed data.
It is worth remembering that the comparison between the model and the observed data in the higher redshift bins is complicated by the increasing difficulty in cleanly separating the star-forming and passive galaxy populations. Moreover, while in this study we have adopted UVJ criteria to separate the star-forming and passive populations, the Peng et al. model assumes the populations are separated based on an evolving rest-frame U − B colour cut. Once again, it is clear that the unique near-IR sensitivity provided by JWST will be crucial for confirming or refuting the steep fall-off in the number density of M ≤ M passive galaxies currently indicated by our z > 2 data.

THE INTEGRATED STELLAR-MASS DENSITY
In order to investigate the redshift evolution of the assembled stellar-mass density, we integrate our intrinsic double Schechter function fits 2 presented in Table 5 between the stellar-mass limits M = 8 and M = 13. To calculate the local stellar-mass density we integrated the double Schechter fit provided by Baldry et al. (2012) between the same limits. We show a comparison between our stellar-mass density results (purple data points) and those from comparable previous literature studies in the left-hand panel of Fig. 13. Where necessary, we have converted the literature results to a Chabrier IMF and the same integration limits. It is important to note that stellar-mass densities based on integrating the Schechter function parameters presented in Tables 4, 5, 9 & 10 are living main sequence stellar masses densities, and do not include stellar remnants. To illustrate the difference, the black data points in the left-hand panel of Fig. 13 show our stellar-mass density results including the contribution of stellar remnants.
It can be seen from Fig. 13 that our stellar-mass density results are consistent with the majority of previous studies but, thanks to the large volume and dynamic range in stellar mass provided by the current data set, carry significantly smaller uncertainties. The dashed purple line shown in the left-hand panel of Fig. 13 is a log-linear fit to our living main sequence stellar-mass densities, and has the form: log 10 (ρ / M Mpc −3 ) = −0.28(±0.01)z + 8.33(±0.01).
We note that this relationship is very similar to the pre- Table 10. The best-fitting observed (upper section) and intrinsic (lower section) double Schechter function parameters for the passive GSMF, where M ≡ log 10 (M / M ) and the units of φ , φ 1 , φ 2 are dex −1 Mpc −3 . The final two columns list the χ 2 and χ 2 ν values of the fits, respectively. The intrinsic double Schechter function parameters were derived assuming σ M = 0.15 dex at all redshifts. We have included the parameters from our own fits to the Baldry et al. (2012) data at z < 0.06 for comparison. To derive the intrinsic Schechter function parameters for the Baldry et al. (2012) data we assumed that σ M = 0.1 dex (Wright et al. 2018). 0.49 ± 0.59 3.00 0.75 Table 11. Integrated stellar-mass densities for the total, star-forming and passive GSMFs over the redshift interval 0.25 ≤ z < 3.75. In all cases the GSMFs have been integrated between the limits of M = 8 and M = 13. For comparison, we have also included our calculation of the integrated stellar-mass densities for the data at z < 0.06 from Baldry et al. (2012). The units of ρ are M Mpc −3 .
total total star-forming star-forming quiescent quiescent Redshift range log(ρ obs ) log(ρ int ) log(ρ obs ) log(ρ int ) log(ρ obs ) log(ρ int ) vious determination of Tomczak et al. (2014), but with a somewhat shallower slope. The solid black line in the lefthand panel of Fig. 13 is a log-linear fit to our stellar-mass densities including stellar remnants, and has the form: In the right-hand panel of Fig. 13 we show a comparison between our stellar-mass density results (including remnants) and the predictions of the SIMBA (Davé et al. 2019), Munich (Henriques et al. 2015) and EAGLE (Furlong et al. 2015) theoretical models. On the same plot, we also show the predicted stellar-mass densities from Madau & Dickinson (2014), based on integrating their fitting function to the evolving cosmic star-formation rate density and converting to a Chabrier IMF.
Overall, it can be seen that there is good agreement between our observational results and the predictions from the latest hydrodynamical and semi-analytic galaxy evolution models. Over the majority of the redshift range explored by this study, the observational and theoretical results are in agreement to within 0.15 dex although, as discussed in Section 4.5, some difference do exist with regard to the precise shape of the evolving GSMF.
Given that the stellar-mass densities predicted by integrating the Madau & Dickinson (2014) fit to the cosmic star-formation rate density include the contribution from stellar remnants, it is clearly of interest to compare them to our direct results. It can be seen from the right-hand panel of Fig. 13 that the Madau & Dickinson (2014) curve begins to overshoot our observed data at z 2.0, reaching a maximum discrepancy of 0.1 dex at z 1.0, before closing again to fall into excellent agreement at z = 0. This effect is well Table 12. Integrated stellar-mass densities for the total, star-forming and passive GSMFs over the redshift interval 0.25 ≤ z < 3.75, but this time including stellar remnants. In all cases the GSMFs have been integrated between the limits of M = 8 and M = 13. For comparison, we have also included our calculation of the integrated stellar-mass densities for the data at z < 0.06 from Baldry et al. (2012) (2018). Where necessary, the literature results have been converted to a Chabrier IMF and recalculated to match our adopted integration limits. The purple data point at z = 0 is based on the local GSMF derived by Baldry et al. (2012). We also show the evolution of log 10 (ρ ) when including stellar remnants, which follows the relation log 10 (ρ / M Mpc −3 ) = −0.31(±0.01)z + 8.44(±0.01) (black points, solid line).
The right-hand panel shows a comparison between the integrated stellar-mass density (including stellar remnants) derived here and the predictions of the SIMBA (Davé et al. 2019), Munich (Henriques et al. 2015) and EAGLE (Furlong et al. 2015) theoretical models. For comparison, we also plot the predicted stellar-mass density from Madau & Dickinson (2014), based on integrating their fitting function to the evolving cosmic star-formation rate density. The Madau & Dickinson (2014) curve has been converted to a Chabrier IMF.
known, and the potential reasons for the discrepancy are discussed at length in Madau & Dickinson (2014). However, it is noteworthy that based on our new observational results, using a Chabrier IMF and including the contribution of stellar remnants, the discrepancy is much smaller than has often been reported in the literature. Given the excellent agreement between our observational data, the Madau & Dickinson (2014) curve and all three theoretical models at z = 2.0 − 2.5, combined with the continued agreement between our results and the theoretical models at z ≤ 2.0, it is tempting to speculate that the SFR estimates used to study the evolution of the cosmic starformation rate density are systematically over-estimated in the redshift interval 0.5 ≤ z ≤ 2.5. Within this context, we note that an evolving off-set between observed and in-trinsic SFRs of this form is predicted by the universemachine model of Behroozi et al. (2019). Alternatively, it is clearly possible that our stellar masses could be systematically under-estimated, perhaps due to a failure to correctly account for the contribution of older stellar populations (e.g. Leja et al. 2020;Carnall et al. 2019). However, if this is the case, any systematic increase in the stellar masses must be limited to a relatively modest 0.1 dex.

The evolving stellar-mass density of passive galaxies
In the top-left panel of Fig. 14 we compare the redshift evolution of the integrated stellar-mass density with the evolu-  Tables 5 & 6). The bottom-left panel shows the same information as the top-right panel, but plotted against lookback time rather than redshift. The bottom-right panel again shows passive fraction versus lookback time, but here the red and blue data points have been shifted by the required amount to match the passive fraction of the purple data points at z = 0.5 (i.e. lookback time 5 Gyr).
tion of the separate star-forming and passive galaxy contributions. As with the total stellar-mass densities, the starforming and passive contributions have been calculated by integrating the best-fitting intrinsic Schechter function fits between the stellar-mass limits M = 8 and M = 13. This panel illustrates the dramatic rise in the stellar-mass density of passive galaxies, from providing an essentially negligible contribution at z ≥ 3, to reaching parity with the star-forming galaxy population between z = 1.0 and z = 0.5, to dominating the stellar-mass density in the local Universe. We note that the data shown in the top-left panel of Fig. 14 is in good agreement with the results derived by Tomczak et al. (2014).
In the top-right panel of Fig. 14 we explore this issue further by plotting the redshift evolution of the passive fraction, f pass = ρ pass /ρ tot , in three different stellar-mass ranges.
The purple data points show the evolution of f pass calculated within the full 8 ≤ M ≤ 13 stellar-mass range. The purple curve is a fit of the form: f pass = a exp (−bz) + c, with bestfitting parameters a = 1.095 ± 0.183, b = 0.271 ± 0.093 and c = −0.392±0.215. Based on the fitted curve, passive galaxies dominate the total stellar-mass density budget at z ≤ 0.75, but contribute 10% by z 3.
As discussed in Section 4.4, the characteristic stellar mass of the double Schechter function fit to the intrinsic GSMF is remarkably stable at M = 10.55 ± 0.1, over the full redshift range studied here. To capture the different behaviour of the passive fraction either side of M , the blue data points in the top-right panel of Fig. 14 show f pass integrating over the stellar-mass range 8 ≤ M < 10.55 (i.e. M < M ) and the red data points show f pass integrating over the stellar-mass range 10.55 ≤ M ≤ 13 (i.e. M ≥ M ). Comparison of the red and blue data points shows a clear downsizing signature, with a significantly higher f pass amongst the M ≥ M galaxies at all redshifts. The fit to the red data points is a linear relation of the form: f pass = az + b, with best-fitting parameters of a = − 0.193 ± 0.014 and b = 0.805 ± 0.031. This relation indicates that the M ≥ M galaxies already have f pass 0.15 by z = 3.5 and reach a passive fraction of f pass = 0.5 by z 1.6. The curve fitted to the blue data points has the form: f pass = a exp (−bz) + c, with best-fitting parameters a = 0.487 ± 0.017, b = 0.759 ± 0.064 and c = −0.026 ± 0.011. In contrast to the high-mass galaxies, this indicates that the M < M galaxies have a passive fraction of essentially zero at z = 3.5 and only reach approximate parity between the stellar-mass contributions of the star-forming and passive populations at z 0.
In the bottom-left panel of Fig. 14 we re-plot the f pass data as a function of lookback time (t lb ). Plotting the data in this fashion serves to highlight the rapid build-up of the passive fraction at high stellar-masses, where the dominant quenching mechanism is thought to be mass quenching (c.f. Peng et al. 2010). Moreover, plotting the data versus lookback time also highlights the fact that the increase in f pass at lookback times of t lb 8 Gyr appears to follow a very similar slope in all three stellar-mass ranges. This is confirmed by the bottom-right panel of Fig. 14, in which the red and blue data points have been shifted vertically in order to match the passive fraction calculated over the full stellar-mass range at a lookback time of t lb 5 Gyr (i.e. z = 0.5).
The results shown in the bottom-right panel of Fig. 14 indicate that the rate of increase in f pass appears to be largely independent of stellar mass at lookback times of t lb 8 Gyr. When combined with the observed stability of the starforming GSMF over this epoch (see Section 5.3), this suggests that the quenching rates (i.e. dρ dt ) at the low and high-mass end of the GSMF must be broadly comparable.

CONCLUSIONS
In this paper we have presented a new derivation of the GSMF over the redshift interval 0.25 ≤ z ≤ 3.75, based on a near-IR selected galaxy sample covering an area of 3 deg 2 and spanning ≥ 4 dex in stellar mass. The powerful combination of a large dynamic range in stellar mass and a large, non-contiguous, cosmological volume has allowed us to robustly constrain both the high and low-mass end of the GSMF. Moreover, by carefully accounting for Eddington bias, we have been able to derive best-fitting Schechter function parameters for both the observed and intrinsic GSMF. By splitting our galaxy sample into its constituent parts, we have investigated the differential evolution of the starforming and passive GSMFs and explored their evolving contribution to the integrated stellar-mass density. Where appropriate, we have compared our new results with previous observational constraints and the predictions of both phenomenological models and galaxy evolution simulations. Our main conclusions can be summarized as follows: (i) We find that a double Schechter function is a better fit to both the observed and intrinsic GSMF over the full redshift range explored in this study, although by z 3.25 the single and double Schechter function fits are indistinguishable.
(ii) The redshift evolution of the GSMF is remarkably smooth in general, and we find no evidence for significant evolution in M . Over the full redshift range explored, the best-fitting values of M are consistent with M = 10.75 ± 0.1 and M = 10.55 ± 0.1 for the observed and intrinsic GSMFs, respectively.
(iii) Motivated by the smooth evolution of the GSMF, we derive a simple evolving parameterization that can provide an accurate estimate of either the observed or intrinsic GSMF at any desired redshift within the range 0 ≤ z ≤ 4.
(iv) Our new determination of the GSMF is in generally good agreement with the predictions of the EAGLE, SIMBA and Munich galaxy evolution models although, in detail, differences still exist. In particular, all SIMBA and Munich models have a tendency to over-predict the number densities of high-mass (M ≥ 11.5) galaxies over the full redshift range.
(v) Splitting our galaxy sample into its constituent starforming and passive galaxy components, we find that the star-forming GSMF is adequately described by a single Schechter function at all redshifts. Moreover, we find that the star-forming GSMF has not evolved significantly since z 2.5.
(vi) In contrast, we find that the passive GSMF has evolved significantly over the redshift range explored by this study, both in normalization and functional form. We find that the passive GSMF is best described by a double Schechter function at z ≤ 1.5, but can be described by a single Schechter function at higher redshifts. Based on our current data set, it is not possible to determine if this change in functional form is intrinsic, or the result of insufficient dynamic range in stellar mass at z ≥ 1.5.
(vii) We find that the Peng et al. (2010) phenomenological model does a qualitatively good job of reproducing the functional form of our star-forming and passive GSMFs at z ≤ 1.5, but appears to perform less well at higher redshifts. That said, more dynamic range in stellar mass will be required to robustly confirm that the components of the highredshift GSMF deviate significantly from the predictions of this continuity-based model.
(viii) Based on our new determinations of the evolving GSMF, we find that the redshift evolution of the integrated stellar-mass density (including stellar remnants) is well described by a log-linear relation of the form: log 10 (ρ / M Mpc −3 ) = −0.31(±0.01)z + 8.44(±0.01) out to z 4. This functional form is in agreement with, although much better constrained than, previous literature results, and in excellent agreement with the predictions of recent theoretical galaxy evolution models.
(ix) We find that the passive galaxy contribution to the integrated stellar-mass budget ( f pass = ρ pass /ρ tot ) evolves by an order of magnitude over the redshift range explored in this study. Within the stellar-mass range 8 ≤ M ≤ 13, we find that passive galaxies dominate the total integrated stellar-mass budget at z ≤ 0.75, but only contribute 10% at z 3. (x) By exploring the evolution of f pass within low stellarmass (M < M ) and high stellar-mass (M > M ) subsamples, we find that the rate of increase in the passive fraction appears to be largely independent of stellar mass at lookback times of t lb 8 Gyr. This suggests that at this epoch the quenching rates (i.e. dρ dt ) at the low and highmass end of the GSMF are broadly comparable.