Evolution of the galaxy stellar mass function: evidence for an increasing $M^*$ from $z=2$ to the present day

Utilising optical and near-infrared broadband photometry covering $>5\,{\rm deg}^2$ in two of the most well-studied extragalactic legacy fields (COSMOS and XMM-LSS), we measure the galaxy stellar mass function (GSMF) between $0.1<z<2.0$. We explore in detail the effect of two source extraction methods (SExtractor and ProFound) in addition to the inclusion/exclusion of Spitzer IRAC 3.6 and 4.5$\mu$m photometry when measuring the GSMF. We find that including IRAC data reduces the number of massive ($\log_{10}(M/M_\odot)>11.25$) galaxies found due to improved photometric redshift accuracy, but has little effect on the more numerous lower-mass galaxies. We fit the resultant GSMFs with double Schechter functions down to $\log_{10}(M/M_\odot)$ = 7.75 (9.75) at z = 0.1 (2.0) and find that the choice of source extraction software has no significant effect on the derived best-fit parameters. However, the choice of methodology used to correct for the Eddington bias has a larger impact on the high-mass end of the GSMF, which can partly explain the spread in derived $M^*$ values from previous studies. Using an empirical correction to model the intrinsic GSMF, we find evidence for an evolving characteristic stellar mass with $\delta \log_{10}(M^*/M_\odot)/\delta z$ = $-0.16\pm0.05 \, (-0.11\pm0.05)$, when using SExtractor (ProFound). We argue that with widely quenched star formation rates in massive galaxies at low redshift ($z<0.5$), additional growth via mergers is required in order to sustain such an evolution to a higher characteristic mass.


The evolution of the galaxy stellar mass function
The galaxy stellar mass function (GSMF) is a measurement of the cumulative effects of physical processes that enhance or hinder star formation within galaxies.These processes include merger events, internal feedback mechanisms (both supernova and active galactic nuclei (AGN) driven) and environmental effects.Understanding the balance between these processes is key to understanding how galaxies have been assembled over cosmic time.Measurements of the local GSMF reveal a steep cut-off in the number of high-mass galaxies above a characteristic mass log 10 ( * / ) ∼ 10.7 (e.g.Baldry et al. 2012).In addition to this, the population of very high-mass galaxies becomes increasingly quenched with time (e.g.Davidzon et al. 2017;McLeod et al. 2020).Many theories have been proposed to explain why there is significant suppression in the star formation ★ E-mail: nathan.adams@physics.ox.ac.uk (often called quenching) of galaxies above this stellar mass, examples include but are not limited to: starvation/strangulation (Larson et al. 1980;Kawata & Mulchaey 2008;McCarthy et al. 2008;Feldmann et al. 2011;Bahé et al. 2013;Feldmann & Mayer 2015), virial shock heating (Birnboim & Dekel 2003;Dekel & Birnboim 2006;Cattaneo et al. 2006) and AGN feedback (Binney & Tabor 1995;Di Matteo et al. 2005;Springel et al. 2005;Bower et al. 2006;Croton et al. 2006;Cattaneo et al. 2009;Fabian 2012;Bongiorno et al. 2016;Beckmann et al. 2017).To adequately test these theories and increase our understanding of how these processes influence galaxy evolution, simulations are required.Key milestones in testing theories such as these include being able to accurately reproduce the observed galaxy population through Luminosity/Mass Functions (see Somerville & Davé 2015;Vogelsberger et al. 2020, for recent reviews).The evolving shape of the GSMF is dependent on all forms of stellar mass growth, including growth via merging events in addition to internal star formation (e.g.Rodriguez-Gomez et al. 2015;Qu et al. 2017;O'Leary et al. 2020).Consequently, to increase our understanding of both quenching and merger rates, stellar mass functions have become a key benchmark for simulations in the past few years (Henriques et al. 2013;Schaye et al. 2015;Pillepich et al. 2018;Lagos et al. 2018) and hence, require observations to be reliable in order to tune and test these simulations.
The advance in deep extra-galactic surveys over the past decade has led to a better measurement of the GSMF at high-redshift.Using predominantly photometric redshifts, several studies have now mapped out the evolution of the GSMF from  = 0-3 (e.g.Fontana et al. 2004;Pérez-González et al. 2008;Marchesini et al. 2009;Pozzetti et al. 2010;Ilbert et al. 2013;Muzzin et al. 2013;Davidzon et al. 2017;Wright et al. 2018) and even up to  7 (e.g.Verma et al. 2007;McLure et al. 2009;Stark et al. 2009;Grazian et al. 2015;Song et al. 2016;Thorne et al. 2020).Despite the increased cosmological volume probed by these studies, there is no consensus on the exact form of the GSMF over this epoch.For example Wright et al. (2018), McLeod et al. (2020) and Thorne et al. (2020) find an approximately constant characteristic mass with redshift, but offset from each other by up to  * ∼ 0.5 dex.Alongside this, other studies of high mass systems suggest merger events are required to explain the observed change in the size-mass relation of these objects (e.g.McLure et al. 2013).The rate of mergers required would consequently influence the evolution of the shape of the stellar mass function at the high mass end.To solve these issues, the precise number densities of galaxies with masses greater than  * are required.Uncertainties in both the stellar mass and number density of these objects can have a much larger impact on the measured shape and characteristic mass of the derived GSMF than the more numerous low-mass sources.This is due to the need to quantify a steep exponential fall off, with relatively few galaxies, that suffer from higher levels of cosmic variance (e.g.Moster et al. 2011).In this work we exploit deep, wide-area optical and IR imaging to measure the GSMF over > 5 deg 2 , leveraging the large area to understand, in particular, the evolution in the number density of the most massive galaxies between 0.1 <  < 2.0.

Measuring accurate photometric redshifts and stellar masses
In the coming decade, the Vera Rubin Observatory (Ivezić et al. 2019) and Euclid (Racca et al. 2016) survey programmes will be conducting broadband observations to further improve both the area and depths achieved within popular multi-wavelength fields such as COSMOS and XMM-LSS, among others.In parallel to these photometric surveys, a number of spectroscopic campaigns are also set to begin operations in the next couple of years.MOONS (Cirasuolo et al. 2012) and WAVES (Driver et al. 2019) are examples of multiobject spectrograph (MOS) surveys and both have an immediate need for high quality photometric redshift and physical parameter estimates in order to plan survey operations and develop final target catalogues ahead of commissioning.Consequently, there is demand for a maintained compilation of broadband photometry and photometric redshifts based on the improved optical and NIR data that has been obtained in recent years in order to best prepare for these upcoming projects (past examples including: Laigle et al. 2016;Alarcon et al. 2020).Moreover, with next-generation telescopes and survey programmes comes next-generation software and analysis pipelines.A workhorse in photometric source extraction for over 20 years has been SExtractor (Bertin & Arnouts 1996), but in recent years there has been a renewed effort in tackling the problem of obtaining accurate flux measurements from images with robust uncertainties and improved handling of source blending.One product of this ef-fort has been ProFound (Robotham et al. 2018), which potentially provides improved photometry, particularly for extended sources, over a variety of wavelengths due to the non-parametric apertures used (e.g.Davies et al. 2018;Hale et al. 2019;Bellstedt et al. 2020).In this study we produce multiple catalogues of broadband photometry and using these different source extraction methods.The use of two source extraction tools, as well as varying the use of Spitzer/IRAC data, allows for any potential bias' in the GSMF due to these different effects to be explored.
This paper is presented as follows: In Section 2 we describe the data used in producing our photometric catalogues.In Section 3 we describe the methods used in source extraction, fitting for photometric redshifts and derivation of basic physical properties of galaxies.In Section 4 we present the measured GSMF in the redshift range 0.1 <  < 2.0 and in Section 5 we present the results of modelling the GSMF, and compare the models with the results from previous studies.In section 6 we explore the time evolution of the GSMF and how our observations compare to results from simulations.We finally present our conclusions in Section 7. Throughout this paper we use the AB magnitude system (Oke 1974;Oke & Gunn 1983) and assume ΛCDM cosmology with  0 = 70 km s −1 Mpc −1 , Ω M = 0.3 and Ω Λ = 0.7.

DATA
The following sub-sections describe the data in the two fields used in the construction of our photometric catalogues and subsequent measurement of the GSMF.

COSMOS
The COSMOS field (Scoville et al. 2007) is one of the most widely studied multi-wavelength fields in extra-galactic astrophysics, with data spanning the X-ray through to the radio domains over ∼ 2 deg 2 of the sky centred on the J2000 coordinates of RA = 150.12deg (10:00:28.6)DEC = +2.21deg (+02:12:21.0).
The imaging data over this field that are used in the construction of our catalogues is derived from four telescopes.The bluest coverage is from the Canada-France-Hawaii Telescope Legacy Survey (CFHTLS; Cuillandre et al. 2012) which has an ultra-deep pointing in the central square degree of the field, and we restrict our analysis to this area for this reason.From this survey we make use of the  * -band.For the optical coverage we take data from the ultra-deep component of the HyperSuprimeCam Strategic Survey Programme DR1 (HSC; Aihara et al. 2018b,a).Near-infrared imaging is sourced from the UltraVISTA survey (McCracken et al. 2012).We make use of the fourth data release of UltraVISTA (DR4) which has a tiered observing strategy, leading to a striped pattern of near-infrared coverage across the field.The 'deep' component is approximately 1 magnitude shallower than the 'ultra-deep' region.3.6 and 4.5 micron photometry is obtained from the Spitzer Extended Deep Survey (SEDS; Ashby et al. 2013), Spitzer Matching Survey of the Ultra-VISTA ultra-deep Stripes survey (SMUVS; Ashby et al. 2018) and Spitzer Large Area Survey with Hyper-Suprime-Cam (SPLASH; Capak et al. 2012).The 5 detection limits for the 2 (2.8) arcsecond photometry in each optical/NIR (Spitzer) bands are outlined in Table 1.This table also highlights the impact of the tiered structure of the UltraVISTA survey.

XMM-LSS
With a much greater area (at 4.5 deg2 ), the XMM-Large Scale Structure field is one of 3 deep fields that make up the Vista Infrared Deep Extra-galactic Observations (VIDEO) survey (Jarvis et al. 2013).Located at RA = 35.5 deg (02:22:00.0)DEC = -4.8deg (-04:48:00.0),our study focuses on regions of the field with high quality HSC observations.We mask parts of the edge of the field due to the lack of overlapping coverage between HSC and VISTA (see Bowler et al. 2020, for more information on overlapping coverage from the different surveys and variable depths).The total area of the field that is used after considering the overlap of telescope footprints is 4.23 deg 2 (giving a total of 5.23 deg 2 when combined with COSMOS).We use the same photometric bands as those used in the COS-MOS field.Unlike in the COSMOS field, the optical coverage in XMM-LSS is not uniform while the near-infrared is more uniform.-band imaging over the full XMM-LSS field was obtained from the CFHTLS Wide survey in addition to the CFHTLS-D1 field which covers a 1 deg 2 patch where  * observations are 1 magnitude deeper than the rest of the field (Cuillandre et al. 2012).HSC SSP covers a different 1.5 deg 2 region of the field (centred on UKIDSS UDS Lawrence et al. 2007) which has 'ultra-deep' coverage.Nearinfrared photometry is derived from the final data release of VIDEO (Jarvis et al. 2013) and Spitzer data is sourced from the Spitzer Extragalactic Representative Volume Survey (SERVS; Mauduit et al. 2012).The depths of the images in each broadband filter are also outlined in Table 1.

PHOTOMETRIC CATALOGUES AND DERIVED PRODUCTS
For our study we define and produce four separate stellar mass estimates for all galaxies identified in the near-infrared.The various stellar mass estimates are determined by using two source extraction algorithms and including/excluding the Spitzer/IRAC 1 and 2 bands.The goal here is to examine the impact of these two variables in the methodology, such as how different measurements of total flux translates to differences in stellar-mass estimation and how the use of Spitzer/IRAC bands affects photometric redshift performance and the final distributions in stellar mass.We refer to these catalogues as 'SExtractor/SE' and 'Pro-Found/PF' when using them in a lengthened/shortened format.When the two Spitzer/IRAC bands are included in any analysis, the suffix '+IRAC' is added to the catalogue name.

SExtractor photometry
Source finding is performed in SExtractor (Bertin & Arnouts 1996) with the   -band image used for object selection.The fiducial photometry was derived using 2 arcsecond diameter circular apertures placed at the location of sources found by SExtractor.As Spitzer has a larger point spread function, for the IRAC bands we use 2.8 arcsecond diameter circular apertures (Bowler et al. 2020).The optical and NIR photometry is aperture corrected using a point-spreadfunction (PSF) generated in each band by PSFE (Bertin 2011) based on cut-outs of point-sources.This is calculated separately for the COSMOS field and over each of the three VISTA/VIDEO tiles in XMM-LSS.For Spitzer photometry, we use an aperture correction derived in the Spitzer handbook1 .Alongside this we also measure the MAG_AUTO parameter from SExtractor to estimate the total flux for each object.The flexible aperture used in the measurement of MAG_AUTO is generated from the   -band detection image.Due to the significantly larger PSF of the two Spitzer bands, carrying over the same aperture from the   band would lead to an underestimation of the total flux.To solve this problem, we derive a correction to translate aperture flux (  ) to total flux ().This correction follows the average value of  ,   =  ,   ×  ,  , for objects in bins of   magnitude and redshift used in the later calculation of the GSMF.The use of this correction is sanity checked against the methodology we use with the P F catalogues and is found to agree well for luminous objects.Errors on the SExtractor photometry are calculated based on local depth maps generated by inserting apertures in empty locations of the field (the same method as applied in Bowler et al. 2020).

Profound photometry
In addition to the SExtractor photometry, we produce P F (Robotham et al. 2018) catalogues selected on a weighted stack of the VISTA- , ,  and  s bands 2 .P F generates two photometric measurements for each object, a 'total' flux and a 'colour' flux.P F operates by iteratively dilating the aperture encompassing each galaxy in each image until it meets the local background.This results in a morphologically derived aperture around each identified object for each photometric band.For our purposes we only make use of the total flux measurements and the associated errors from the P F output.P F flux errors are calculated by combining the errors resulting from the sky root mean square (RMS), errors from the sky subtraction process and object shot noise.Because this dilation process is performed for each photometric band, the larger Spitzer/IRAC PSF is taken into account and thus no alterations are required to obtain total fluxes in these bands.To validate both the ProFound and SExtractor photometry.We compare our measured total photometry against measurements made in other catalogues which have targeted these fields.These comparison catalogues are the COSMOS2015 catalogue for the COSMOS field (Laigle et al. 2016) and the Subaru/XMM-Newton Deep Field (SXDF) catalogue (Mehta et al. 2018).We find that the vast majority of our photometric bands are consistent ( 0.1 mag median offsets) with the measurements from these reference catalogues.The methodology we use to determine total SExtractor Spitzer fluxes for bright, resolved sources is also found to be consistent with these catalogues.The greatest difference in photometry is found when running ProFound on Spitzer images in XMM-LSS.For bright sources ( < 20), the photometry is consistent with our SExtractor catalogue and the reference catalogue.However, for fainter objects ( > 20) the measured flux in ProFound is found to be fainter than measured with SExtractor (up to 0.2 mags offset at  ∼ 24).Similarly, colour distributions match those of the reference catalogues, with the exception of the aforementioned ProFound IRAC offsets in XMM-LSS.

Photometric Redshifts
The procedure we follow for obtaining photometric redshift estimates is the same as that used in Adams et al. (2020), with the only exception being the use of observer-frame NIR selection instead of optical selection.To summarise here, we make use of the minimising  2 code L P (Arnouts et al. 1999;Ilbert et al. 2006) to fit templates of galaxies, active galactic nuclei (AGN) and stars to our SExtractor derived aperture fluxes.These templates are modified for dust attenuation according to the Calzetti et al. (2000) extinction law with E(B-V) = 0, 0.05, 0.1, 0.15, 0.2, 0.3, 0.6, 1.0, 1.5.The result is a Probability Density Function (PDF) for the redshift and a simple classification as a likely galaxy, star or QSO-like object.The template sets used are those from Ilbert et al. (2009) and are sourced from Polletta et al. (2007) and from Bruzual & Charlot (2003).To conduct object classification, template spectra for AGN from Salvato et al. (2009) and stars from Hamuy et al. (1992Hamuy et al. ( , 1994)); Bohlin et al. (1995); Pickles (1998); Chabrier et al. (2000) were also fit.Photometric errors are set to a minimum of 5 per cent in flux during the template fitting process, this is to alleviate the consequences of using templates that probe the colour space discretely while the real galaxy population is continuous.

Zero-point calibration with spectroscopic samples
The two fields have been targeted by a large number of spectroscopic campaigns which can be used to calibrate our methods and examine photometric redshift accuracy.We make use of the spectroscopic catalogue compiled by the HSC team3 .Included are spectra from the VIMOS VLT Deep Survey (VVDS; LeFèvre et al. 2013), z-COSMOS (Lilly et al. 2009), Sloan Digital Sky Survey (SDSS-DR12; Alam et al. 2015), 3D-HST (Skelton et al. 2014;Momcheva et al. 2016), Primus (Coil et al. 2011;Cool et al. 2013) and the Fiber-Multi Object Spectrograph (FMOS; Silverman et al. 2015).From these we select only those with high quality flags (>95 per cent confidence) to ensure secure redshifts are being used.Together these provide a spectroscopic sample of 22,409 in COSMOS and 35,125 in XMM-LSS.
We use these spectroscopic redshifts to examine the accuracy of our photometric redshift estimates.In addition to this, we also make use of L P functionality to make iterative adjustments to the zero-points of each photometric filter in order to optimise the results against the spectroscopic sample.Minor shifts in the zero-points can occur as a result of inaccurate filter transmission functions, through biases within the choice of SED templates and from the calibration of the images.The inclusion of a very large sample of spectroscopically confirmed objects from a variety of different surveys minimises the risk of introducing an additional bias through calibration on a non-representative sample.For each catalogue, we run L P once on the spectroscopic sample to obtain the zero-point corrections, these offsets are applied and the entire field is then run.We show the results of this process in Table 2 and the majority are small compared to the errors ( 0.1 mags).

Photometric redshift accuracy
The quality of the photometric redshift catalogues can be described with two numerical values.The outlier rate, defined as percentage of objects satisfying abs(   −  ℎ )/(   + 1) > 0.15, and the Normalised Mean Absolute Deviation Hoaglin et al. (NMAD;1983), defined as 1.48 × median[|   − ℎ |/(1+   )].These two values quantify a) the rate at which the photometric redshift method produces a redshift value that is in significant tension to the measured spectroscopic redshift and b) the spread around    −  ℎ in a manner that is resistant to influence from the relatively small number of extreme outliers.Comparing the two fields, COSMOS has greater depth and uniformity while XMM-LSS is shallower and wider, and has around 1 magnitude variability in its optical coverage .It is therefore expected that XMM-LSS would produce photometric redshifts of lower quality.The results of this comparison are displayed in Table 3.For each field, we produce two sets of photometric redshift estimates, one with and one without the inclusion of the Spitzer IRAC 3.6 and 4.5 m bands.The addition of the Spitzer IRAC bands to COSMOS makes minimal difference in the quality of the photometric redshifts.However, redshifts are found to decrease in quality for the faintest objects in the XMM-LSS catalogue.
Table 3. Summary of the photometric redshift statistics.Displayed are the outlier rates and Normalised Median Absolute Deviation (NMAD) of the catalogues when compared to a large spectroscopic sample.We cut show these numbers in cuts of i-band magnitudes to enable comparison with results from Laigle et al. (2016).The cut of   < 22 corresponds to the brightest 50 per cent of the spectroscopic sample we compare to in COSMOS and the brightest 60 per cent in XMM-LSS.

Stellar mass determination
With photometric redshifts and object classification determined for each source, we proceed to measure the stellar mass.This is performed by fixing the redshift to the best-fit value (template and redshift with minimum  2 ) and rerunning LePhare using the total flux measurements, rather than aperture fluxes.For SExtractor we use fluxes from the MAG_AUTO parameter and for ProFound this is magt.Compared to the aperture fluxes, the total fluxes are essential to making accurate measurements of the normalisation of the SED for resolved objects and hence the total luminosity of the galaxy and stellar mass.In the case that an object has a spectroscopic redshift from one of the surveys described in Section 3.2.1, this value is used instead of the photometric redshift.

MEASURING THE GSMF
We select objects for use in measuring the GSMF based on a number of criteria to maintain purity and completeness.
(i) The source exists in both SExtractor and ProFound derived catalogues.The majority of sources that fail this are either artefacts on the edge of manual masking or are a consequence of the initial ProFound selection on a VISTA stack verses just the   -band.Such sources are also removed by the   magnitude cuts detailed below.
(ii) The source has a 2 arcsecond aperture magnitude following the condition   < 24.5 in COSMOS and   < 23.2 in XMM-LSS.This corresponds to a SNR cut of 8 and is employed to minimise the potential for contamination while enabling  * to be well constrained up to  ∼ 2.
(iii) The source has a best fit SED that is a galaxy or AGN with a redshift between 0.1 and 2.0 ( 2 Gal/QSO <  2 Star ).In the case a source has a spectroscopic redshift, that value is used in place of the photometric redshift.
(iv) We apply an upper limit on the quality of the photometric redshift of  2 < 250 (removing the worst 0.5 per cent of objects).
This results in a sample of ∼ 320, 000 galaxies in the combined COSMOS and XMM-LSS fields used in measuring the GSMF.We present the galaxy number counts in each redshift bin in Tables A1 and A2.The inclusion of Spitzer bands leads to a significant reduction in the number of highly massive galaxies at redshifts 0.1 <  < 2. This is the result of Spitzer bands breaking degeneracies between stellar templates and certain combinations of galaxy template, leading to a number of these massive objects being reclassified as stars.The use of ProFound photometry over SExtractor leads to minimal difference to the general population of objects, cases will be discussed in Section 5.1.1.

The 1/𝑉 max method
We first compute the GSMFs using the 1/Vmax methodology (Schmidt 1968;Rowan-Robinson 1968).We determined the  max for each galaxy by redshifting the best-fitting template (from the 2 arcsecond aperture photometry) until the object no longer meets our  -band magnitude limit.
The GSMF is then determined using: where  is stellar mass, Δ is the bin width which we set to 0.25 dex and  max, is the maximum volume for which galaxy  could have been successfully detected.

Stellar mass completeness
Towards lower stellar masses, galaxies become intrinsically less luminous.This ultimately leads to a regime where the detection limits of the data are reached and galaxy number counts begin to fall as they are lost to noise.As the science goals of this study focus on the massive end of the GSMF, we adopt a conservative approach while still probing a significant mass range.In our model fitting procedures we elect to only use bins of redshift and mass where over 95 per cent of the galaxy sample are brighter than the 8 2 arcsecond aperture detection limit of the respective field.Due to the shallower coverage in XMM-LSS, the criteria for completeness is approximately 0.5-0.75dex higher in stellar mass than in COSMOS.We apply a simple correction to the survey areas based on the fraction of the field that is occupied by other sources or masked by foreground stars.For COSMOS this is 15 per cent and XMM it is 7 per cent.This corrects the GSMF for the probability that sources are highly blended (> 50 per cent overlap) with other sources or significantly effected by bright foreground stars.

Cosmic Variance
Our measurements of the GSMF are based on data that only probe a limited volume of the Universe.As a result, they are susceptible to biases that are a consequence of large-scale fluctuations in Here we display the cosmic variance for each field if treated independently and the result of combining the fields using the cosmic variance calculator from Moster et al. (2011) and our measured number counts.COSMOS is shown in gold, XMM-LSS in blue and the combination of the two fields in black.Where XMM-LSS becomes incomplete the cosmic variance value for the combined case is just the COSMOS cosmic variance.density in the galaxy distribution.This is commonly referred to as cosmic variance ( 2  ).As we are measuring the GSMF across a wide range of mass and redshift, there is no single quantitative value that can be used to describe this effect.In order to model the effects of cosmic variance on our measurements, we use the treatment from Moster et al. (2011), which provides an estimate of the cosmic variance as a function of both stellar mass and redshift.Our dataset consists of two fields with differing area and dimensions, thus allowing us to mitigate some of the effects of cosmic variance.Where both the XMM-LSS and COSMOS fields are used in measuring the GSMF we calculate the cosmic variance for each field independently ( 2  , ) and combine the values together with a co-moving volume weighting (Equation 7 in Moster et al. 2011).The output is a percentage error on the GSMF, and so to combine the fields this value is converted back into variance ( 2  ) using our galaxy number counts.When data from the XMM-LSS field drops out of consideration due to its shallower depth, cosmic variance is determined from the area of the COSMOS field alone.Fig. 1 shows the results of our cosmic variance calculations for two redshift bins (0.2 <  < 0.3 and 0.75 <  < 1.0), highlighting how the increased area from including the XMM-LSS field, coupled with combining two widely separated fields, results in a factor of ∼ 2 decrease in the cosmic variance uncertainty.

Eddington bias
The steep drop in the GSMF beyond the knee can lead to a bias in the derived number densities of the most massive galaxies due to Eddington bias (Eddington 1913).All galaxies in the sample have an uncertainty in the derived stellar mass, however as low-mass galaxies are significantly more numerous this leads to more galaxies scattering to higher stellar mass than the number that scatter to lower masses.To account for this effect and hence, determine an estimate of the intrinsic GSMF, we require an estimate of the uncertainty in the stellar masses derived for our sample.With this as a part of our study.We show the probability of a certain mass scatter (   =  original −  perturbed ) as a result of photometric errors in the SED fitting process.The lower plot shows identical data to the upper plot, however we have logged the -axis to reveal the low-amplitude wings of the distribution.The non-parametric distribution derived directly from our data is shown as the thick red line.The black dashed line is the best fit Gaussian to this data and the orange line is for a Gaussian multiplied by a Lorentzian.distribution we can then deconvolve (or in reality, fit a convolved double Schechter form) to our observations to determine the intrinsic GSMF.To measure the uncertainty in the stellar masses derived in this study, we repeat both the photometric redshift measurement and the SED fitting process after perturbing the photometry of our sources according to the photometric errors in each band.This process is repeated multiple times to produce the distributions shown in Fig. 2. Based on this, we examine three possible methodologies to uncover the intrinsic stellar mass function from this observed distribution in our analysis described in Section 4.5.

The measured galaxy stellar mass functions
The GSMFs, as measured from the samples produced from our four catalogues, are presented in Fig. 3.They probe stellar masses from 7.75 < log 10 ( ) < 11.75 and are split into nine redshift bins between 0.1 <  < 2.0.
To each GSMF we fit a double Schechter function (Schechter 1976, Eqn. 3) using a Markov-Chain Monte Carlo (MCMC) implemented in emcee (Foreman- Mackey et al. 2013).In this redshift range, past studies have found the double Schechter functional form better describes the GSMF due to the underlying bimodality in the galaxy population (e.g.Strateva et al. 2001;Driver et al. 2006;Baldry et al. 2012;Ilbert et al. 2013;Davidzon et al. 2017) and we clearly see an upturn at the low mass end of our GSMFs.
A series of priors are applied to prevent parameters from flipping due to the symmetry of the double Schechter functional form shown in Equation 3. The normalisation (Φ 1 & Φ 2 ) follows the condition Φ 1 > Φ 2 , the low-mass slopes ( 1 &  2 ) are limited to being between [-1.8,1.5] and [-3.0,-0.9]respectively for the two components.To compare against past studies, we make the same assumption that each Schechter component has a single, shared  A1 and A2).The black line shows the median of the MCMC results for the SExtractor+IRAC data when fit with a double Schechter function convolved with our Eddington correction.The grey shaded region shows the region contained within 1 of the model fit and is based on 10,000 random samples of the final posterior.The gold line and shaded region follows the same process with the ProFound+IRAC data.The high redshift bins in the right column have the x-axis truncated to higher masses to focus on the complete regime.In the final panel the Eddington corrected ProFound models are shown simultaneously and cut where data is incomplete.
value of  * in the range 10 < log 10 ( * /M ) < 12.We only fit to data points where the bins in stellar mass are greater than 95 percent complete.The MCMC is set up with 500 walkers that burn in for 100,000 steps before conducting a further 20,000 steps for use in mapping the posterior distributions.Each walker is distributed uniformly in the parameter space and limited by the above priors.
We perform the fitting procedure four times, once on the observed GSMF and three times with the double Schechter function modified with the convolution with one of three Eddington bias methods that we describe below.First, we modify the fit to be a convolution of the double Schechter function with a Gaussian distribution, where the standard deviation (  ) of the distribution is calculated by fitting a Gaussian to the measured scatter in masses shown in Fig. 2.This is a commonly used method in studies of the GSMF (e.g.Wright et al. 2018).In the second method we extend this model by multiplying the Gaussian with a Lorentzian in the same manner as described in Ilbert et al. (2013); Davidzon et al. (2017).This adds extended wings to the function, which more adequately reproduces the distribution, however it continues with the assumption that the mass scatter is symmetric.Any asymmetry is important to account for, as it means that there is a greater probability of scattering to lower masses than towards higher masses (due to the photometric redshift uncertainty), and this would impact on the measured GSMF.Therefore, in the third case, we do not fit any parametric distribution to the mass scatter, instead we use the smoothed histogram of the scatter convolved directly with the double Schechter function.This should improve upon the use of the analytic forms because it captures the observed asymmetry in the mass scatter ().
Visual inspection of the distribution in mass scatter derived in Section 4.4 reveals there to be minimal dependence on redshift.The two lowest redshift bins are slightly broader as a consequence of photometric redshift uncertainty, but these bins contain a much higher fraction of galaxies with spectroscopic redshifts (30 per cent).As a result, the real scatter within these bins is likely much smaller.Using the method which convolves a Gaussian distribution with the measured GSMF, we find values with 0.08 <   < 0.10 across all redshifts.For the second method, which uses a Gaussian distribution multiplied with a Lorentzian, we find 0.43 <   < 0.58.We note that the   values for the two methodologies are not directly comparable due to the different functional forms.Since minimal evolution was found, we remove the redshift dependence on the Lorentzian component that was introduced in Ilbert et al. (2013) to minimise total fit parameters.The resultant formula for the extended wings is thus , which is equivalent to fixing the redshift to 1.0 in the original formula from Ilbert et al. (2013).The range of 0.43 <   < 0.58 for the Gaussian component is in agreement with the findings of Ilbert et al. (2013) and slightly larger than the value found by Davidzon et al. (2017).In addition, Grazian et al. (2015) and Davidzon et al. (2017) both find evidence for redshift and stellar mass dependence on the measured mass scatter when approaching the completeness limited regime.The probable explanation for the lack of such a dependence in our data is the conservative SNR cuts that have been implemented.
Following previous studies, we fix the   values in our final fits to 0.09 in the Gaussian case and 0.5 for the Gaussian × Lorentzian case.The shapes of these distributions are presented in Fig. 2. We discuss the impacts of this correction on the measured GSMF in Section 5.2.
Our preferred method for correcting for Eddington bias is to use the histogram presented in Fig. 2 directly in the convolution.This method directly uses the results of the perturbed catalogues and captures the subsequent asymmetry found in the distribution.Such an asymmetry has previously been described in recent studies exploring the Eddington bias (Grazian et al. 2015).The best-fitting double Schechter function fit parameters using this Eddington bias correction are presented in Table 4. Corner plots showing the posterior probability distribution for the parameters in this model will be provided in an online resource.For completeness, we also present our results without the application of the Eddington bias correction in Appendix A alongside the results obtained using the various parametric fits to the mass scatter.
Our double Schechter model is only fit up to stellar masses of 10 11.75  .While there are a small number of objects with stellar masses greater than this limit in our sample, these are increasingly likely to be subject to forms of contamination such as AGN activity, source blending, misclassification of stars or artefacts.In the high redshift regime we are also unable to constrain the lowmass Schechter component.If we instead fit with a single Schechter function for  > 1.25, we find the shift in the fit parameters to be minimal (< 1) compared to results obtained with a double Schechter function.For consistency we proceed with the double Schechter functional form for all redshift bins.
Visual inspection of the measured GSMFs reveal a change in the shape of the massive-end between the redshift bins of 0.75 <  < 1.0 and 1.0 <  < 1.25 (see the last panel of Fig. 3).This evolution is present in the results obtained using both source extraction methods (SExtractor/ProFound) and both sets of photometric redshifts (including/excluding Spitzer data).Inspection of the redshift distributions reveal no significant features, such as peaks or troughs, within these redshift bins.The total shift in the high-mass component amounts to around 2 sigma in  * and the normalisation Φ 1 between these two redshift bins.So, a combination of statistical errors/cosmic variance or an unknown systematic could be the driver of such a change.

Changes in the GSMF with varying methodology
Within each redshift bin, we measure a total of four stellar mass functions.First we compare the results with and without the inclusion of Spitzer/IRAC [3.6] and [4.5] data, and second we compare the GSMF derived from SExtractor based photometry in comparison to that derived with Profound.While the low mass component of the GSMF is consistent between these different methods, we find some differences in the results at stellar masses greater than  * .

The impact of including Spitzer/IRAC photometry on the measured GSMF
The most immediately apparent feature is the offset between the GSMFs that include or exclude Spitzer/IRAC data in the redshift bins of 0.5 <  < 0.75 and 1.0 <  < 1.25.This is due to two effects.Firstly, many high-mass objects that lie in the 0.5 <  < 0.75 bin have a different redshift solution ( ∼ 0.1) when Spitzer/IRAC is included.They consequently have smaller masses in this lower redshift bin and their influence is negated by the high number densities within this mass-redshift range.The cause of the different redshift solutions is the uncertainty of the SED slope redder than the   band, with a red slope favouring the higher redshift solutions and a blue slope favouring low redshift solutions.Secondly, there is a likely contamination from stars in the 0.5 <  < 0.75 and 1.0 <  < 1.25 bins.Many high mass objects have smooth red optical slopes that turn over around the  or   bands.With a limited wavelength range, some black body-like spectra and red galaxies are indistinguishable.Introducing the Spitzer/IRAC bands vastly increases the  2 of the galaxy models and reduces the  2 of the stellar models for a significant number of these high mass objects.Consequently these objects no longer meet our selection criteria (either the  2 increases above 250 or the classification switches to a star) when using the values associated with the inclusion of the mid-infrared bands.Both of these cases are examples of degeneracies between template sets as a result of the use of a limited number of broadband filters.Thus we find that the inclusion of the [3.6] and [4.5] bands makes a significant difference to the derived number density of the most massive galaxies in our data.Therefore, we focus on the GSMFs measured with the inclusion of the Spitzer/IRAC when discussing redshift evolution and comparisons to simulations.

The impact of varying the choice of source extraction software
While the global mass distributions between SExtractor and Pro-Found based catalogues are broadly very similar, there are individual cases where mass estimates can vary widely between objects.Issues can typically be put down to artefacts affecting the photometry, proximity to bright sources and disagreement between the two source extraction measurements in the Spitzer bands.Instances of significant differences in mass estimations tend to occur in high redshift and/or low mass objects that fall within our incomplete regime and so do not affect our final results.
At low redshifts, galaxies become increasingly resolved and so any systematic variation between the SExtractor and ProFound photometry would be expected to be more apparent.In our measurements we find any such variation between the two to be very small, with the low-mass components ( <  * ) between 0.1 <  < 1.5 being highly consistent between the two source extraction methods.We find that at the lowest redshift (0.1 <  < 0.2) the ProFound based GSMF produces more galaxies of very high mass compared to the SExtractor derived measurements.This is then reversed at higher redshifts  > 0.5, where the ProFound GSMF produces fewer galaxies of very high mass.However, these differences are of relatively low significance (of order 1) and demonstrate that over the redshift and mass ranges probed in this study, the choice of source extraction software makes no significant difference to GSMF results.

The intrinsic GSMF corrected for Eddington bias
To recover the intrinsic GSMF from our observations, which is affected by the uncertainties in the estimate of stellar mass, we trial three different methods described in Section 4.4.In the first method, we assume that the scatter induced by uncertainties on stellar mass are described by a Gaussian distribution of   = 0.09.We find that this imposes minimal changes to the shape of the high-mass end of the MF, corresponding to a shift in  * of order 0.05 dex lower when compared to the observed GSMF.With the second method, where we convolve a Gaussian × Lorentzian combination with   = 0.5, we find the shift between the observed and intrinsic parameters to be more significant with  * shifting to lower masses by ∼ 0.1 dex compared to the observed GSMF.For completeness, we display the The impact that differing models of the Eddington Bias have on the measurement of the intrinsic GSMF at high masses.Shown is the intrinsic GSMF at 0.5 <  < 0.75 as measured with SExtractor photometry when recovered using the three methods applied in this study.In blue we show the results of using a simple Gaussian to model mass errors, in red we expand the model to include Lorentzian wings and in black we show the results when using a non-parametric approach.Shaded regions indicate the 1 sigma uncertainty and are derived from 10,000 random samples of the posterior probability distribution.The edges of the shaded regions are made more bold to assist in readability.
MCMC results for these two methods in Tables A3 and A4 respectively.Each of these methods make a fundamental assumption in that the scatter in the derivation of the stellar mass is symmetric in logarithmic stellar mass space.We find this assumption to work best when redshifts are confident e.g. if we conduct the photometry scatter procedure on just the mass calculations and assume the photometric redshift is correct or the objects have spectroscopic redshifts.However, when the uncertainty on the photometric redshifts is included, the measured distribution is found to be broader and more asymmetric.This results in the Gaussian × Lorentizan method underestimating the amount of scattering in certain regimes (small up-scatter and most of the down-scattering), even with the extended wings of the function.It is for this reason that we elect to instead use the measured distribution of mass scatter to convolve with the intrinsic double Schechter function (see Fig. 4 for an example of the impact).
The GSMF, when corrected with this distribution in mass scatter, undergoes a similar shift in the Schechter parameters to that of the Gaussian × Lorentzian.The shift in  * is 0.12 dex compared to the observed GSMF and all the fit parameters lie within 1 of the results found with the Gaussian × Lorentzian method.While the broad wings of the distribution of mass scatter are very small in probability beyond shifts of 0.3 dex (see Fig. 2), the nature of the GSMF spanning many orders of magnitude in number density around the knee requires these wings to be modelled in order to capture the impact of a small number of objects scattering into the less populated, high-mass bins.This highlights that the intrinsic GSMF is very sensitive to the strength of the Eddington bias correction and could be a key source of inconsistency between results of observational studies.The best-fit parameters for the double Schechter function when using our preferred non-parametric Eddington correction are shown in Table 4.
The use of high completeness spectroscopic surveys would reduce the uncertainty on stellar mass due to photometric redshift uncertainties (Pozzetti et al. 2010;Moustakas et al. 2013;Leauthaud et al. 2016).While such studies have been limited to the brighter and more massive objects (e.g.log 10 (/ ) > 10.5), it is these objects that are most at risk at having their number counts inflated by Eddington Bias.With new surveys and instruments using multi-object spectrographs coming online in the coming years (e.g.DEVILS and WAVES), such studies will soon be capable of measuring the GSMF to lower masses and higher redshifts (Davies et al. 2018;Driver et al. 2019)

Comparisons to previous studies
To put this study into greater context we compare against a number of past studies.The studies selected for these comparisons are Davidzon et al. (2017) which utilised only the COSMOS field in their study with the Laigle et al. ( 2016) catalogue, Wright et al. (2018) which uses a combination of the GAMA (Driver et al. 2011), COSMOS (Davies et al. 2015;Andrews et al. 2017) and 3D-HST surveys (Brammer et al. 2012;Skelton et al. 2014;Momcheva et al. 2016) and the recent study conducted by McLeod et al. (2020) that uses components of the COSMOS & XMM-LSS fields.The results of the the double Schechter fits conducted in these studies are presented alongside our own in Fig. 5 and all were calculated with the same cosmological model as used in this study.
Examining the Schechter parameters individually, the strongest agreement is that of the high mass normalisation (Φ 1 ), where there is very close agreement between our results and those from McLeod et al. (2020).In contrast, the slope of the high mass component  1 is found to agree more with the results from Wright et al. (2018), however this is subject to degeneracy based effects from our poorer capability to constrain the low-mass slope  2 .Inspection of the corner plots show that the more negative values of  1 are driven by the steeper slopes found in our broad uncertainties of the low-mass slope (see accompanying online resources).
The Eddington bias corrections implemented in previous work all use different functional forms.Wright et al. (2018) andMcLeod et al. (2020) use corrections that are independent of mass and redshift as in our study, whereas Davidzon et al. (2017) 2020), who fold their mass uncertainties into their model fitting procedures.They find a constant  * = 10.78 and mildly evolving values between 10.8 and 10.9.The work conducted in Leja et al. (2020) is an example of a study that included redshift evolution terms within a model fitting procedure that did not bin by redshift.Such a methodology can result in a smoothing of the evolution of the model parameters, which has its pros and cons depending on the timescales examined and assumptions made (e.g.fixing the slope of  1 ).For the propose of this study, we have elected to assume no particular evolutionary form for the GSMF.
These past studies, combined with our work trialling different functional forms to account for Eddington bias, suggest that up to 0.15 dex in variation of the characteristic mass  * can be attributed to implementing different methodologies.While significant enough to cause changes of a few , it is not enough to explain the ∼ 0.5 dex discrepancy among  * values found between these studies, 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00  2020) is a modified value from Baldry et al. (2012) indicating that there are further systematics at play between the data sets.

Time evolution of the high-mass component of the GSMF
Including Spitzer/IRAC photometry is found to be important in obtaining accurate galaxy classification and photometric redshifts, we therefore discuss the redshift evolution in the context of the GSMFs that make use of these data.We show the time evolution for the best-fit Schechter parameters (Φ 1 ,  * ,  1 ,  2 ), corrected for Eddington bias via our non-parametric method, in Fig. 5. Examining these best-fit Schechter parameters reveals evolution to be initially driven by the normalisation (Φ 1 ) of the double Schechter components, which increases with time.This evolution is found to be much stronger at higher redshifts, with a change of over 0.5 dex between  = 2 through to  = 1.25.For  < 1.25, the normalisation Φ 1 stabilises and the evolution of the mass function becomes dominated by evolution in  * (there is also a possible flattening of  1 , but this is degenerate with  2 which is poorly constrained at higher redshift).
Under the assumption that  * is constant for 0 <  < 2, the best fitting value of log 10 ( * / ) = 10.73 ± 0.04 is found for both SExtractor and ProFound derived photometry.However, we find that a constant  * model gives a poor quality of fit, with  2 red = 3.4 (3.4) for SExtractor (ProFound).The observed rise in the value of  * at  < 1.25 is present in all versions of the intrinsic GSMF that we produce with the differing Eddington corrections.We therefore fit a simple linear model to the evolution of  * of the form log 10 ( * / ) =  + , where  and  are free parameters and fit using minimisation of  2 .Introducing this simple time dependence into  * produces significantly better fits ( 2 red = 1.4 (2.4) with a Pearson correlation of 0.7 (0.75) for SExtractor (ProFound)) .The evolutionary fits can be described by, Δ log 10 ( * / ) Δ = −0.164± 0.047 SExtractor −0.113 ± 0.054 ProFound (4) a 2-3 disagreement with the evolution described by a constant  * .However, the observed evolution in  * could also be attributed to a steeper evolution over a smaller redshift range, rather over the entire redshift range probed (see top-right panel of Fig. 5).Such an evolution could be indicative of a break in the balance of growth channels that maintains a constant  * above  ∼ 1 while Φ 1 increases.Two such channels include star formation and the rate of major mergers.Observations and simulations show that the contribution towards new stellar mass from both mechanisms falls at  < 1.5 in massive galaxies (Rodriguez-Gomez et al. 2015;Tomczak et al. 2016;Qu et al. 2017;O'Leary et al. 2020).The natural consequence of this would be to stabilise the high-mass end of the GSMF towards lower redshifts, as we find (e.g. the bottom right panel of Fig. 3 for  < 0.75).

The evolving quenched fraction of galaxies
The effects of the evolving in situ star formation rate and merger rates of galaxies results in the changing shape of the GSMF.In order to probe this further, we examine the population of massive galaxies over the redshift range 0.5 <  < 1.5 in more detail.We separate our sample of massive galaxies into two broad specific star formation rate (SSFR) bins, based on the same SED template used to derive the stellar mass.We term these bins 'Star-forming' and 'Quenched', which are defined as bins of −8.5 > log 10 ( /yr −1 ) > −10.0 and −10.0 > log 10 ( /yr −1 ) respectively.The star forming bin is selected to cover the majority (∼ 95 per cent) of the star forming main sequence at  ∼ 1.5 as found in both observational and simulation based studies (e.g.Sparre et al. 2015;Tomczak et al. 2016).The quenched galaxies are thus defined as falling below this star forming main sequence.We do not enforce our definitions to be redshift dependent (i.e. a certain level above or below the evolving star-forming main sequence) to capture the global fall in star formation post-cosmic noon and produce results similar to those that use rest-frame colour selection (e.g.Davidzon et al. 2017;McLeod et al. 2020).We restrict the redshift range to 0.5 <  < 1.5, as SSFR estimates for galaxies at redshifts  < 0.5 would become increasingly unreliable as the rest-frame UV exits the wavelength coverage of our observations.
We reproduce the GSMF for these two populations in Fig. 6.Our definitions for star-forming and passive galaxies exclude a small number of galaxies undergoing starbursts −8.5 < log 10 ( /yr −1 ), and consequently the GSMFs in Fig. 6 do not sum directly to the total GSMF presented in Fig. 3.We find the change in total number density of galaxies with  >  * to be driven by a strong evolution in the passive population.As a consequence of this trend, the most massive component of the GSMF becomes increasingly dominated by more passive galaxies while the number densities of star-forming galaxies are more constant, agreeing with Davidzon et al. (2017) and McLeod et al. (2020).Assuming a constant star formation rate, a borderline quenched galaxy with mass  ≥  * , log 10 ( /yr −1 ) = −10 at  = 1.5 would grow of order 0.15 dex in stellar mass during the ∼ 4 Gyr that passes between 0.5 <  < 1.5.The majority of the objects in our high mass sample will produce fewer stars than predicted from this overly simple estimation, as most start with SSFR values which are lower and there is a global trend for the SSFR to decrease with time.Consequently, the growth of galaxies through solely main sequence evolution is insufficient to produce the observed increase in high-mass galaxy number densities.
Since these quenched systems are not expected to form enough stars to move between bins of stellar mass ( 0.25 dex), the two primary ways of increasing the quenched number densities at high masses are consequently merger events and the addition of newly quenched systems taken from the star-forming population.As our star-forming number densities at high mass remain near constant, this implies that these galaxies must be replenished by slightly lower mass star-forming galaxies at approximately the same rate at which they undergo quenching.Major mergers (∼ 1 : 4 mass ratio) are capable of generating mass gains of order ∼ 0.1 dex for passive galaxies on top of internal star formation.Mergers themselves can trigger starbursts in the aftermath of the event.However, at later times ( < 1), such mergers are increasingly 'dryer' as massive galaxies tend to be gas poor (e.g.De Lucia & Blaizot 2007;Lin et al. 2008;Davidzon et al. 2016;Tomczak et al. 2017).Merger events come at the cost of reducing the total number of galaxies and its effects will be embedded in the evolving shape of the GSMF.In addition to providing pathways for stellar mass growth of passive galaxies, merger events can also be used to solve issues surrounding the radial size of passive galaxies.At higher redshifts ( ∼ 1.5), quenched galaxies are found to be compact relative to passive systems observed at  = 0 (e.g.Williams et al. 2010;Wuyts et al. 2011;McLure et al. 2013).The passive nature of these galaxies requires that external processes must be responsible for the observed shift in galaxy radius.The study by McLure et al. (2013) finds that combinations of major and minor mergers can simultaneously reproduce Figure 6.The GSMF when broken down into star forming and passive components by SSFR.The quenched population is shown with the red lines and the blue lines show the star-forming population.Darker shading indicates mass functions at higher redshifts.Shading widths indicates observational errors through both cosmic variance and   .The bin 0.5 <  < 0.75 is not shown to reduce clutter and is found to be near identical to the bin 0.75 <  < 1.0.The star-forming population is found to be near constant at higher masses while the passive population is found to evolve significantly between 1.5 <  < 0.75.both the growth in stellar mass and radius observed in these systems between  ∼ 1.5 and today.

Comparison to simulations
In this section we compare our measured GSMF with the semianalytic simulation (SAM) SHARK (Lagos et al. 2018) and the hydrodynamical simulations Simba (Davé et al. 2019) and EAGLE (Schaye et al. 2015).In Fig. 7 we show the GSMF results from these alongside our intrinsic (Eddington-bias corrected) double Schechter functions.The SHARK-SAM has the GSMF as one of physical measurements that it is tuned to reproduce at  = 0, 1 and 2. As a result it is unsurprising to find excellent agreement with our results with the low mass slope and location of the primary 'knee'.However, there is an excess in the number density of the most massive objects (log 10 (/ ) > 11.5) in the simulation.With the close matching to all other components of the GSMF, the likely source for this discrepancy could be in part the choice of GSMF which was used to tune the simulation.The Lagos et al. (2018) models elect to calibrate against studies such as Muzzin et al. (2013) and Wright et al. (2018), all studies with  * values on the higher end of the parameter space covered by observational studies (10.8 < log 10 ( * / ) > 11.0).For our comparisons we utilise the GSMF derived from total masses provided from the primary SHARK run.The study by Lagos et al. (2018) finds that implementing a 30kpc aperture to measure stellar mass leads to minimal impact for log 10 (/ ) < 12.0.
A similar situation is present when comparing against the results of the Simba hydrodynamical simulation (Davé et al. 2019).Simba implements a new, torque-limited accretion model (Anglés-Alcázar et al. 2017) for cold gas alongside more conventional Bondi accretion (Bondi 1952) for hot gas.The energy built up in accretion is used to fuel feedback and quench galaxies.While this model for black hole growth and AGN activity was considered to be more physically motivated over previous simulations, it is found to still over produce galaxies of very high stellar mass.Our observations  reinforce the findings of Davé et al. (2019), where the largest discrepancy in high mass galaxies is at  ∼ 2 (even when compared to the SHARK results) and that there is a slight underproduction of galaxies around  * at low redshifts ( < 1).This results in a much shallower 'knee' than would be described with an exponential cut off.The authors explain possible causes of over-production of massive galaxies include over-merging of galaxies that blend, due to the use of a friend-of-friends (FoF) algorithm to count star particles, and the over-production of large halo masses.The EAGLE simulation (Schaye et al. 2015) implements AGN feedback through inputting a fraction of the accreted gas as thermal energy into the local surroundings.It is found to produce mass functions with slightly higher low-mass normalisation but much sharper cut offs at the high mass end.This results in an over production of galaxies on the low mass slope by a factor of around 2, but a closer match at  ≥  * .The inconsistencies at the lowmass end are thought to originate from the lack of quenching of lower mass galaxies at higher redshifts, or from the need to implement more burst-like star formation histories (Furlong et al. 2015).These GSMF measurements implement three-dimensional apertures of radius 30kpc when calculating stellar masses, the size of which is found to reduce the mass assigned to the largest, most massive systems in the simulation and prevent some of the overmerging issues described in Davé et al. (2019).This was found to reduce the number density of galaxies at log 10 (/ ) ∼ 11.5 by around half a dex, imposing a steeper cut off in the number of high-mass systems.In the work conducted in Furlong et al. (2015), the EAGLE GSMF is fit with a double Schechter function and is found to exhibit a strong evolution in  * from 0.1 <  < 2.0.The value of log 10 ( * / ) = 10.44 ± 0.08 at  = 2 is notably small, even lower than the results from McLeod et al. (2020) which lie at the lower extrema of the range of observational constraints.The increase from log 10 ( * / ) = 10.74 ± 0.05 at  = 1 to log 10 ( * / ) = 10.93 ± 0.03 at  = 0.1 matches our ProFoundbased GSMF to within 1, while the value of  * is lower at  = 0.1 from our SExtractor-derived GSMF.However, it is worth noting that Furlong et al. (2015) suggest that such an evolution in the characteristic mass could be the result of overly strong AGN feedback limiting the production of massive galaxies between 2 <  < 4. It remains difficult to determine if faults lie within observations (mass errors, systematics etc.) or with the simulations (e.g.fine tuning AGN feedback) due to the sensitivity off all of these factors on the exponential decline in number density.
Within EAGLE, the origin of this stellar mass growth in the most massive galaxies arises from multiple sources.Qu et al. (2017) show that in the redshift range of  < 1.5, around 68 per cent of the massive galaxies (log 10 (/ ) > 11) in EAGLE undergo at least one major merger (at least 1:4 mass ratio, leading to > 0.1 dex growth) and have the average fraction of stellar mass originating from outside the primary galaxy increase by a factor of two (from a median of 10 per cent to 20 per cent), with most of this external contribution fuelled by these major merger events.At  = 1, over half of the simulated galaxies with log 10 (/ ) > 11 are defined as quenched (log 10 ( /yr −1 ) < −10), falling significantly below EAGLE star-forming main sequence (Furlong et al. 2015).This leads to a similar fraction of quiescent galaxies to our sample.These consequently experience mass gains of less than 0.15 dex over the 4Gyr between 0.5 <  < 1.5 under the simple assumption of a constant star formation rate.Creating this evolving value of  * while maintaining a near constant value of Φ 1 thus likely requires a balance in merger events and internal star formation in order to sustain the growing number density of  >  * systems.

CONCLUSIONS
Utilising new photometric catalogues generated from optical and near-infrared data in the COSMOS and XMM-LSS fields, we have measured the GSMF over the redshift range of 0.1 <  < 2.0, covering 60 per cent of the age of the Universe.The use of the two fields greatly reduces the impact of cosmic variance, due to both the wider area and the fact that they are widely separated on the sky, and allows for tight constraints on the GSMF at  >  * .Simultaneously, the depth of the photometry available allows for constraints on on the GSMF of log 10 (/ ) = 7.75 for galaxies at  0.1 and log 10 (/ ) = 9.75 for galaxies at  2.0.We have investigated the use of ProFound and SExtractor source extraction software and also the impact of Spitzer/IRAC [3.6] and [4.5] photometry on the GSMF.Our main conclusions are summarised as: (i) The inclusion of Spitzer/IRAC [3.6] and [4.5] photometry alleviates the degeneracies in select cases where the SEDs of red, massive galaxies in the redshift ranges 0.5 <  < 0.75 and 1.0 <  < 1.25 are confused with low redshift or stellar templates.Thus the inclusion of these data leads to fewer contaminants and a small decrease in the derived  * (of up to 0.1 dex) in these redshift ranges.
(ii) Both SExtractor and Profound derived photometry produce consistent faint-end components of the GSMF.Differences between the mass functions at higher masses are greater when examining the extreme parts of our results (e.g. at the lowest and highest redshifts) but the resultant double Schechter fits are found to agree to within ∼ 1.
(iii) The measured GSMF is found to disagree with the assumption that the characteristic mass  * is constant with time between 0.1 <  < 2.0 at the 3 (2) sigma level for the SExtractor (ProFound) derived results.Such an evolution in  * between 0.5 <  < 1.5 can be seen in some (but not all) previous work and also in the predictions from the EAGLE hydrodynamical simulation.However, significance is low and caveats in both the understanding of observational systematics and AGN feedback strength in simulations (Furlong et al. 2015) mean a claim of an evolving  * is presently very mild.
(iv) Eddington bias, and the methodology used to correct for it, is found to be highly influential on the shape of the high-mass end when attempting to retrieve the intrinsic GSMF from observations.There is presently still no consensus on the handling of observational and systematic errors that can impact stellar mass estimates (see also discussions within Grazian et al. 2015;Davidzon et al. 2017).For our data, we find the correction required to be asymmetric and poorly described by commonly used analytic forms.When comparing to the observed GSMF, applying a simple Gaussian treatment to the Eddington bias is found to reduce  * by around 0.05 dex.With the Lorentzian wings added into the description, the shift in  * doubles to around 0.1 dex.Utilising a non-parametric form, based on the measured error distribution, leads to  * shifts of 0.12 dex when compared to the observed GSMF.
(v) When splitting our galaxy sample by specific star formation rate (SSFR), our results confirm the findings of previous studies that show increasing number densities of quenched galaxies are responsible for the rise in the GSMF for  >  * (e.g.Davidzon et al. 2017;McLeod et al. 2020).Examining growth channels for stellar mass in the EAGLE simulation show this to be the result of a combination of internal star formation and merger events.This is because internal star formation alone would amount to less than 0.15 dex in stellar mass growth for the majority of the population of massive galaxies ( >  * ).Between 0.5 <  < 1.5, the constant number densities of star-forming galaxies at  >  * indicate that these galaxies quench at approximately the same rate that lower mass galaxies replace them through their own in situ star formation.
(vi) Comparisons to simulations reveal that the semi-analytic model SHARK and hydrodynamical simulation Simba both overproduce massive galaxies at low and intermediate redshifts.Likewise the EAGLE is found to over-predict the number of low-mass galaxies by a factor of around 2. This highlights that even in the era where such simulations can be fine tuned to better replicate observations, discrepancies still exist.It also further emphasises the conclusions drawn by these studies that additional work is required on all sides, from the observational data with which to calibrate/compare simulations against, to the development of the physical models used in simulations.However, we find the evolution of the high-mass component of the EAGLE simulations replicates our observations of an evolving  * .
Table A1.The median and standard deviations of the MCMC samples for the double Schechter function when using SExtractor based total photometry.Priors are set to ensure that  1 & log 10 (Φ 1 ) refer to the high mass component and  2 & log 10 (Φ 2 ) refer to the low mass component.Alongside are the number of galaxies from each of the two fields that contribute to the GSMF as shown in Fig. 3.These mass functions are based on the raw observations and are not corrected for Eddington bias.

Figure 1 .
Figure 1.The estimated cosmic variance for two example redshift bins.The solid lines show 0.2 <  < 0.3, the dashed lines show 0.75 <  < 1.0.Here we display the cosmic variance for each field if treated independently and the result of combining the fields using the cosmic variance calculator fromMoster et al. (2011) and our measured number counts.COSMOS is shown in gold, XMM-LSS in blue and the combination of the two fields in black.Where XMM-LSS becomes incomplete the cosmic variance value for the combined case is just the COSMOS cosmic variance.

Figure 2 .
Figure2.The shapes of the three Eddington bias corrections implemented as a part of our study.We show the probability of a certain mass scatter (   =  original −  perturbed ) as a result of photometric errors in the SED fitting process.The lower plot shows identical data to the upper plot, however we have logged the -axis to reveal the low-amplitude wings of the distribution.The non-parametric distribution derived directly from our data is shown as the thick red line.The black dashed line is the best fit Gaussian to this data and the orange line is for a Gaussian multiplied by a Lorentzian.

Figure 3 .
Figure3.The GSMF from our analysis of the COSMOS and XMM-LSS data.In blue and red are the results from SExtractor and ProFound, while navy and brown data points denote the points measured after including Spitzer/IRAC photometry.Unfilled data points indicate those with less than 95 per cent completeness and are not used in any fitting procedure.Displayed data points are based on raw observations and do not correct for Eddington bias (see TablesA1 and A2).The black line shows the median of the MCMC results for the SExtractor+IRAC data when fit with a double Schechter function convolved with our Eddington correction.The grey shaded region shows the region contained within 1 of the model fit and is based on 10,000 random samples of the final posterior.The gold line and shaded region follows the same process with the ProFound+IRAC data.The high redshift bins in the right column have the x-axis truncated to higher masses to focus on the complete regime.In the final panel the Eddington corrected ProFound models are shown simultaneously and cut where data is incomplete.
Figure4.The impact that differing models of the Eddington Bias have on the measurement of the intrinsic GSMF at high masses.Shown is the intrinsic GSMF at 0.5 <  < 0.75 as measured with SExtractor photometry when recovered using the three methods applied in this study.In blue we show the results of using a simple Gaussian to model mass errors, in red we expand the model to include Lorentzian wings and in black we show the results when using a non-parametric approach.Shaded regions indicate the 1 sigma uncertainty and are derived from 10,000 random samples of the posterior probability distribution.The edges of the shaded regions are made more bold to assist in readability.
implement a minor redshift dependence in the Lorentzian terms.Furthermore, Wright et al. (2018) use a Gaussian of width   = 0.1 dex while Ilbert et al. (2013); Davidzon et al. (2017) use the Gaussian × Lorentzian distribution of widths 0.5 dex and 0.35 dex respectively.McLeod et al. (2020) uses a log-normal distribution of width   = 0.15 dex to correct for the Eddington bias, leading to larger changes in their fit parameters, with  * changing by up to 0.2 dex compared to the raw observations.Other examples of recent studies include Thorne et al. (2020) and Leja et al. (

Figure 5 .
Figure 5.The time evolution of the double Schechter parameters as measured with photometry derived by SExtractor (blue) and ProFound (red) and including the use of Spitzer/IRAC photometry in both redshift and mass calculations.These data points are for the fits that utilise the non-parametric correction for Eddington bias.Since the parameter Φ 2 is largely unconstrained for a significant portion of the redshift range, we do not show those results in this figure.Alongside our results we display the results of Davidzon et al. (2017); Wright et al. (2018); McLeod et al. (2020) in black, gold and green respectively.The lowest redshift data point from McLeod et al. (2020) is a modified value fromBaldry et al. (2012)

Figure 7 .
Figure 7.The intrinsic Galaxy Stellar Mass Function (GSMF) from Table4with Eddington bias removed.SExtractor derived results are in blue while ProFound is in red with the shading indicating the 1 confidence interval.Presented alongside are examples of previously conducted simulations.For data from external sources, we present the closest redshift bin available for our comparisons.The studies used are the primary results from the SHARK run conducted inLagos et al. (2018) as the black dotted lines, EAGLE(Schaye et al. 2015) as yellow dashed lines and Simba as the dot/dashed purple lines(Davé et al. 2019).

Table 1 .
Summary of the 5 detection depths within the COSMOS and XMM-LSS fields.Depths are calculated in 2 diameter circular apertures (2.8 for IRAC bands due to Spitzer's larger point spread function) that were placed on empty regions of the image.Values are grouped into Deep (D) and Ultra-Deep (UD) sub-regions.The XMM-LSS field has a single ultra-deep pointing from HSC centred on the UDS field and so a third of XMM-LSS has this deeper optical coverage.Similarly the  * coverage has a single square degree of ultra-deep coverage in a separate part of XMM-LSS as part of the CFHT Legacy Survey observing programme.

Table 2 .
The list of zero-point offsets applied to the 2 arcsecond aperture photometry after optimising the template fitting process against a large spectroscopic sample.Offsets listed are in magnitudes and are calculated separately for each field.

Table 4 .
The best-fitting double Schechter function parameters derived from our observed GSMF when corrected for Eddington bias through the direct implementation of the mass scatter over fitting an analytic form.Priors are set to ensure that  1 & log 10 (Φ 1 ) refer to the high mass component and  2 & log 10 (Φ 2 ) refer to the low mass component.

Table A2 .
The median and standard deviations of the MCMC samples for the double Schechter function when using ProFound based total photometry.Priors are set to ensure that  1 & log 10 (Φ 1 ) refer to the high mass component and  2 & log 10 (Φ 2 ) refer to the low mass component.Alongside are the number of galaxies from each of the two fields that contribute to the GSMF as shown in Fig.3.These mass functions are based on the raw observations and are not corrected for Eddington bias.

Table A3 .
The median and standard deviations of the MCMC samples for the double Schechter function when corrected for Eddington bias by convolving with a Gaussian distribution using a width of   = 0.09.Priors are set to ensure that  1 & log 10 (Φ 1 ) refer to the high mass component and  2 & log 10 (Φ 2 ) refer to the low mass component.

Table A4 .
The median and standard deviations of the MCMC samples for the double Schechter function when corrected for Eddington bias through the scattering of photometry in the photo-z process, leading to a convolution with a Gaussian × Lorentzian functional form with a strength of   = 0.50.Priors are set to ensure that  1 & log 10 (Φ 1 ) refer to the high mass component and  2 & log 10 (Φ 2 ) refer to the low mass component.