Constraints on the Faint End of the Galaxy Stellar Mass Function at z ≃ 4–8 from Deep JWST Data

We analyze a sample of 3300 galaxies between redshifts z ≃ 3.5 and z ≃ 8.5 selected from James Webb Space Telescope (JWST) images in the Hubble Ultra Deep Field and UKIDSS Ultra Deep Survey field, including objects with stellar masses as low as ≃108 M ⊙ up to z ≃ 8. The depth and wavelength coverage of the JWST data allows us, for the first time, to derive robust stellar masses for such high-z, low stellar mass galaxies on an individual basis. We compute the galaxy stellar mass function, after complementing our sample with ancillary data from CANDELS to constrain the GMSF at high stellar masses ( M>M* ). Our results show a steepening of the low stellar mass end slope (α) with redshift, with α = −1.61 ± 0.05 at z ≃ 4 and α = −1.98 ± 0.14 at z ≃ 7. We also observe an evolution of the normalization ϕ * from z ≃ 7 to z ≃ 4, with ϕz≃4*/ϕz≃7*=130−50+210 . Our study incorporates a novel method for the estimation of the Eddington bias, which takes into account its possible dependence both on stellar mass and redshift, while allowing for skewness in the error distribution. We finally compute the resulting cosmic stellar mass density and find a flatter evolution with redshift than previous studies.


INTRODUCTION
Over the past decades, substantial progress has been made in our understanding of when and how galaxies have assembled their stellar masses.The redshift and stellar-mass ranges that could be studied have become steadily larger, thanks to the progressive availability of powerful galaxy surveys that probe increasingly wider areas and/or provide deeper photometric coverage (e.g., Ilbert et al. 2013;Caputi et al. 2015;Grazian et al. 2015;Davidzon et al. 2017;Deshmukh et al. 2018).
However, until now, the study of low stellar-mass galaxies (M < ∼ 10 8 −10 9 M ⊙ ) at high redshifts has been very difficult, as the existing galaxy surveys were not deep enough to enable the investigation of this region in parameter space.A notable exception were massive galaxy cluster fields, where the lensing magnification helped to detect intrinsically faint, distant objects Corresponding author: Rafael Navarro-Carrera navarro@astro.rug.nl in the background, which could not be seen in blankfield surveys of comparable depths (e.g., Karman et al. 2017;Bhatawdekar et al. 2019;Kikuchihara et al. 2020;Vanzella et al. 2021;Santini et al. 2022).Nonetheless, sufficiently deep imaging to detect these objects was, in most cases, only available with the Hubble Space Telescope (HST) up to ≃ 1.6 µm.This implied that, for galaxies at z > ∼ 3, only the rest-UV light could be traced and any derived stellar-mass estimate relied on long extrapolations of the galaxy spectral models.
The advent of JWST is now radically changing this situation.By providing high-sensitivity imaging up to mid-infrared wavelengths, JWST is allowing us for the first time to find the building blocks for galaxy formation at high redshifts and directly measure the light emitted by their evolved stellar populations, which is necessary to obtain robust estimates of their stellar masses.Therefore, we can now investigate early stellar mass assembly down to the typical stellar masses of local satellite galaxies.This is the main goal of this work.
The galaxy stellar mass function (GSMF) is one of the most important statistical tools to understand galaxy evolution across cosmic time.The evolution of this function with redshift reveals how galaxies of different stellar masses have assembled until producing the final distribution that we see in the Universe today.A Schechter function (Schechter 1976) is commonly used as a parametrization for the GSMF.This function is characterized by a power-law shape with slope α (lowmass end slope) up to a characteristic stellar mass (M), and an exponential decline at the high-mass end.The GSMF constitutes an important statistical tool to understand galaxy formation and reproducing the GSMF at different redshifts is one of the main challenges of galaxy formation models.Indeed a number of studies have provided strong constraints on the GSMF at low redshifts (e.g., Pozzetti et al. 2007;Ilbert et al. 2010;Wright et al. 2017;Bates et al. 2019).Some of these works have also investigated the GSMF's dependence on galaxy types, such as star-forming/quiescent galaxies (e.g., Davidzon et al. 2017), central/satellite galaxies (e.g., Yang et al. 2009), and galaxy morphology (Ilbert et al. 2010).Moreover, the evolution of the GSMF with cosmic time provides information on the different physical processes that shape galaxy assembly (e.g., supernova and super massive black-hole feedback and mergers) and clues to understand star-formation and quenching mechanisms in galaxies of different stellar masses (e.g., Peng et al. 2010;Weaver et al. 2022).
In particular, the value of the GSMF faint-end slope is a crucial constraint for galaxy formation models.Generally, these models predict a steepening of this low-mass end with redshift, meaning that low stellar-mass galaxies were relatively more abundant at earlier cosmic times.However, for the reasons explained above, determining observationally the GSMF low-mass end at high redshifts has been very challenging.A number of studies have attempted to constrain the GSMF faint-end slope, both from deep blank and lensing fields, but no consensus has been reached yet on its value and redshift evolution (e.g., Grazian et al. 2015;Song et al. 2016;Bhatawdekar et al. 2019;Kikuchihara et al. 2020).With ultra-deep JWST imaging in hand, we can now study the GSMF faint-end slope at high redshift in a much more robust way than ever before.
In its first year of operations, JWST has been able to study unknown or poorly constrained galaxy populations due to its unprecedented near and mid infrared sensitivity.Several studies are making use of these JWST unique capabilities to probe different aspects of galaxy evolution at high redshifts, both from the analysis of normal galaxies (e.g., Bagley et al. 2022;Iani et al. 2022;Finkelstein et al. 2023;Jin et al. 2023;Kokorev 2023;Rinaldi et al. 2023;Pérez-González et al. 2023), as well as active galactic nuclei (e.g., Yang et al. 2023).
In this work, we make use of very deep JWST Near Infrared Camera (NIRCam) observations to study the evolution of the GSMF (focusing on the low-mass end) at high redshifts up to the Epoch of Reionization, between z ≃ 3.5 and z ≃ 8.5.We consider the publicly available NIRCam images of the UKIDSS Ultra Deep Survey (UDS; Lawrence et al. 2007) and Hubble Ultra Deep Field (HUDF; Beckwith et al. 2006;Illingworth et al. 2013) fields for our analysis.We study a sample of ≃ 3300 galaxies between z ≃ 3.5 and z ≃ 8.5, using for first time JWST data to compute the GSMF, focusing our analysis on the evolution of the low-mass end slope with cosmic time and the derived evolution of the cosmic stellar mass density.
This paper is structured as follows: in Sect. 2 we describe the dataset and in Sect. 3 the SED fitting procedure and estimation of photometric redshifts and stellar masses.The methodology used for computing the GSMF is presented in Sect.4, followed by the results on the GSMF evolution in Sect. 5. Finally, we investigate the evolution of the low-mass end slope and cosmic stellar mass density in Sect.6 and Sect.7, respectively.We adopt a cosmology with Ω m,0 = 0.3, Ω Λ,0 = 0.7, and H 0 = 70 km s −1 Mpc −1 .Magnitudes are given in the AB system (Oke & Gunn 1983).Densities of the galaxies are measured on comoving scales.We adopt the Chabrier (2003) initial mass function (IMF) in a stellar mass range of M : 0.1 ≃ 100 M ⊙ to estimate stellar masses.In this paper, all of the stellar masses taken from the previous studies are converted to those estimated with the Chabrier (2003) IMF.

General overview
In this work we make use of photometric datasets from two different fields partly covered with JWST : HUDF (≃ 12 arcmin 2 ) and UDS (≃ 7 arcmin 2 ).The study of two different fields helps to mitigate the effects of cosmic variance, while providing a statistically larger galaxy sample populating the low and intermediate stellar-mass regimes.The region of the HUDF with JWST coverage includes the entire Hubble eXtreme Deep Field (XDF, Illingworth et al. 2013), which benefits from the deepest HST observations (typical depth of 30 ABMag at 5σ) and some of the deepest public JWST imaging.For both fields, we also make use of ancillary HST data to obtain our galaxy catalogs, as explained below.
Albeit very useful for the low stellar-mass end, the JWST fields considered here are too small to properly constrain the GSMF around the turnover stellar In our study we include both the XDF and nearby region of the HUDF field that lies within the JWST pointings.Two regions with different average depth have been defined: deep (XDF) and shallow (the rest of the field).The depths and areas are summarized in Fig. 1 and presented in Table 1.
Our procedure for the NIRCam data reduction is the same as the one described in Rinaldi et al. (2023).The version of the JWST pipeline used for is 1.8.2 and the Calibration Reference Data System (CRDS) pipeline mapping (pmap) version is 1018.We have added extra steps to the original pipeline provided by the Space Telescope Science Institute (STScI; Bushouse 2020), following a similar approach to the one of Bagley et al. (2022).These new steps include correction for stripping (1/f noise), stray-light effects (wisps) and residual cosmic rays (snow balls).As mentioned in Rinaldi et al. (2023), photometry has been checked for consistency against the default pipeline without any additional steps.
All the images have been resampled to pixel scale of 0. ′′ 03/pix and drizzled to a mosaic using the Hubble Legacy Fields catalog as reference to align the images.
Our procedure for the NIRCam data reduction is similar to that for the HUDF.For the UDS, however, due to presence of stray light in some of the filters we need to apply a mask to exclude regions of severely contaminated photometry.The reduction in area after the masking is significant, but less than 50%, accounting for the loss of all the detector area affected by the stripes.The effective area of the UDS-PRIMER dataset is then ≃ 7 arcmin 2 .
As a result of overlapping exposures within the field, the depth is non-uniform.We have identified two regions: deep and shallow.The area and depth in each band JWST can be found in Table 1 and the map of the field is presented in Fig. 1.
The HST ancillary imaging is obtained from CANDELS-UDS data (Galametz et al. 2013), including newer observations with HST/WFC3 that where not present in the original Galametz et al. (2013) analysis.The coverage is 0.4µm − 1.6µm, with 7 photometric bands.In the optical: ACS/WFC F435W, F606W, F814W.In the near-infrared: WFC3/IR, F105W, F125W, F140W and F160W.We have detected sources and measured photometry in both the HUDF and UDS field regions with JWST coverage, making use of the software Source Extractor (SExtractor, Bertin & Arnouts 1996).We considered a stacked image that combines all the available JWST bands as detection map, and then ran SExtractor in dual-image mode to measure photometry in every band.
As parameter values for source dection, we used those corresponding to the hot-mode described in Galametz et al. (2013), optimized for the detection of faint sources, which is ideal for low stellar-mass galaxies at high redshift.
We measured source photometry following Rinaldi et al. (2023), i.e., in fixed 0. ′′ 5 diameter circular apertures and Kron apertures (i.e., flux auto; Kron 1980).We corrected the aperture photometry to total considering the curve of growth of the corresponding point spread function in each filter, using WebbPSF (Perrin et al. 2014).For each source with apparent magnitude < 27 (as described in Rinaldi et al. (2023)), we chose the brightest flux (either the corrected aperture flux or flux auto).For all fainter sources, which can be con-sidered point-like even in JWST images, we adopted the corrected aperture fluxes for our photometric catalog.
We corrected all source fluxes for Galactic extinction with the help of dustmaps (Green 2018), which is in good agreement with the prescriptions from Schlafly & Finkbeiner (2011) for HST filters.To take into account possible flux calibration uncertainties and because SExtractor has been found to underestimate photometric errors (e.g., Sonnett et al. 2013), we have imposed a minimum error of 0.05 mag to all the photometry.
Flux upper limits have been estimated for all sources with image coverage, but not detected, in any given band, by measuring the local background r.m.s from multiple, circular apertures of 0. ′′ 5 diameter, randomly positioned around each source.In the case that no photometric coverage is available in a given band (because the source is not contained in the field of view), we ignored that band for the SED fitting explained below.
Our source catalogs have been cleaned for Galactic stars by cross-matching our catalog with Gaia Data Release 3 (Gaia Collaboration et al. 2022) and by excluding all sources with high stellarity parameter (i.e., CLASS STAR > 0.8) that lie on the stellar locus of the (F435W − F125W) vs (F125W − F444W) diagram, following the same criterion as in Caputi et al. (2011).

Photometric redshifts and stellar masses
We have estimated the photometric redshifts and stellar masses of our galaxies by performing SED fitting using the software LePHARE (Arnouts & Ilbert 2011).We considered a set of stellar population synthesis models from Bruzual & Charlot (2003, hereafter BC03), including a single stellar population and a series of exponentially declining star formation histories with τ = 0.01, 0.1, 0.3, 1.0, 3.0, 5.0, 10.0, and 15 Gyr.All models have been constructed for two possible metallicities, namely solar (Z ⊙ = 0.02) and sub-solar (Z = 0.2 × Z ⊙ = 0.004).
We also considered stellar population synthesis models from Leitherer et al. (1999, hereafter SB99) in our analysis.These additional models consist of a set of five young galaxy templates with ages spanning 10 6 -10 8 yr and constant star formation rates between 0.01 M ⊙ yr −1 and 10 M ⊙ yr −1 .All of them correspond to sub-solar metallicities of Z = 0.008 and Z = 0.001.All the relevant paramters used for the SED fitting both for BC03 and SB99 can be found in Table 2.
These templates are characterized by incorporating the contribution of nebular continuum emission, especially relevant for young star-forming galaxies.Some   From the plot, it can be seen that the difference between the limiting mass of HUDF and UDS is 1 dex or larger, consistent with the deeper data that we use in HUDF.
previous work explores the effect produced by nebular Note-1 In the case of SB99 models, an identical configuration with respect to BC03 for the extinction law, E(B-V), IMF, redshift interval, and cosmology has been used.However we opted for only 6 ages for the SB99 run because nebular emission is only significant for very low age stellar populations. 2The age grid is designed to provide a finer sampling for young ages.
emission in SED fitting using broadband photometry (e.g., Stark et al. 2013, de Barros et al. 2014).The effect of nebular emission is to enhance the measured flux, resulting in a overestimation of stellar masses (Bisigello et al. 2019), especially for young and low-mass galaxies.
There have been previous attempts to study this effect (e.g., Stark et al. 2013;Duncan et al. 2014;Salmon et al. 2015;Grazian et al. 2015. Stark et al. (2013) finds that the stellar masses obtained using SED fitting of Spitzer Space Telescope (Werner et al. 2004) data where overestimated by a factor of 2 to 4 when accounting for nebular line and continuum emission in redshifts z ≃ 6 − 7. Duncan et al. (2014) finds that stellar masses estimated including nebular emission are significantly lower.The models used for our work (SB99) only account for nebular continuum, so the effect on stellar masses is expected to be milder than the one including both contributions.
For all BC03 models we considered a Chabrier (Chabrier 2003) initial mass function (IMF), and the SB99 models we used were re-scaled to a Chabrier IMF.We convolved the spectral models with a reddening law following the prescription of Calzetti et al. (2000) and Leitherer et al. (2002), with color excess values of 0 ≤ E(B − V ) ≤ 1.5 equally distributed in steps of 0.1 for both BC03 and SB99.
We used our photometric catalog as the input catalog for SED fitting with LePHARE, considering 3σ flux upper limits in all cases, but source non detection in the particular band.All templates that produce fluxes greater than the 3σ upper limits are rejected by LePhare.The adopted stellar mass is obtained from the set of models (BC03 or SB99) that has the lowest value of χ 2 .The diagnostic of photometric versus spectroscopic redshift for the whole sample (Fig. 2) shows that the fraction of catastrophic outliers (defined in our case as (z phot − z spec )/(1 + z spec ) > 0.15) is 12.3%, and the normalized median absolute deviation of the sample is 0.0419.
Within our sample there are a total of 1804 galaxies that lie above the completeness limit (described in Sect.4.2) in the redshift range z ≃ 3.5 to z ≃ 8.5.Table 3 shows the number of galaxies used for computing the GSMF in each redshift interval.The resulting stellar mass versus redshift plot is shown in Fig. 3.

1/Vmax method
We apply the 1/Vmax method (Schmidt 1968, ,described in more detail in Sect.4.1) to compute the GSMF.The 1/Vmax method is non-parametric and does not assume any particular functional form for the GSMF, but requires binning in stellar mass.As shown by Davidzon et al. (2017), 1/Vmax provides good constraints to the GSMF, although other methods such as STY (Sandage et al. 1979) are also commonly used in the literature (e.g., Caputi et al. 2015).
To compute the GSMF with the 1/Vmax method, we binned our galaxy samples into redshift and stellar-mass bins.The corresponding comoving number density of galaxies for each stellar mass bin (M 0 , M 1 ) can be de-termined using , (1) where V max is the comoving volume of each galaxy corrected in the following way: Here V bin min and V bin max are the comoving volumes corresponding to the redshift bin (z bin min , z bin max ).V z lim is the comoving volume at z lim , the maximum (limiting) redshift at which each galaxy would lie within the flux completeness of the sample.The value of z lim for each galaxy can be obtained from the following relation between flux and redshift: The 1/Vmax method provides a set of data points to which a functional form can be fitted.We adopt a Schechter function: In order to determine the uncertainties which arise from the 1/Vmax fitting, we have implemented a χ 2 method, by calculating: between the observed Φ 1/V max (M) and fitted Φ(M|α, M * , ϕ * ) data for each stellar-mass bin.Where σ is the error (standard deviation) derived for each of the 1/Vmax points.For this purpose, we have developed a parameter space sampler, which provides the best fitting set of parameters (α, M * , ϕ * ), as well as the 1σ confidence region for each parameter pair.
After obtaining the GSMF with the 1/Vmax method considering the union of both fields for the calculation (taking into account the effective volume and completeness correction on a galaxy-by-galaxy basis), we fitted the resulting data points with a Schechter function by using a χ 2 minimization method developed for this purpose.
In our calculations, we considered the effects of flux completeness and derived stellar-mass completeness and applied the necessary corrections (Sects.4.2 and 4.3).
We also study the influence of Eddington bias (Eddington 1913), which can potentially affect the shape of the derived GSMF, as detailed in Sect.A.6.

Photometric completeness
For assessing completeness, we introduced simulated point sources in the detection image.After calculating the fraction of sources recovered for each bin in flux, we have derived a completeness curve (left panel of Fig. 4) which represents the fraction of sources detected versus the total number of sources in the field for each (stacked) magnitude.HUDF (shown in blue) and PRIMER-UDS (shown in red) do not have uniform coverage.Nonetheless, they can be divided into a deep region and shallow region with almost uniform depths.The limiting magnitudes for each band of JWST are shown in Table 1.The depth of HUDF field is almost uniform for JWST data, however this is not the case for HST ancillary observations, where the XDF region is significantly deeper with respect to the HUDF (e.g., Illingworth et al. 2013).
For the CANDELS complementary data (used to compute the high-mass end of the GSMF), we have obtained the completeness curves for CANDELS-UDS and GOODS-S field from Mortlock et al. (2015), which follows a similar methodology to ours for computing the photometric completeness.
In this study we have considered only galaxies that lie above 65% completeness.The number of galaxies with magnitude above the limiting magnitude is ≃1800, conforming the final sample that has been used for the computation of the GSMF.A correction derived from the completeness curves has been introduced when calculating the GSMF as a weight to each galaxy.The limiting stacked image magnitude at 65% completeness is shown for every region next to the each curve, which can be helpful to asses the difference in depth between the fields used here.
Although faint (low-mass) high-redshift galaxies are expected to be behave like point like sources in the detector, some studies explored the effects of galaxy sizes on the photometric completeness (e.g., Finkelstein et al. 2015).For this reason, we have re-obtained the completeness curve considering extended sources described by a Sérsic profile (Sérsic 1963) with index n = 1 and half-light radius following a uniform distribution between the range of usual values found in (Finkelstein et al. 2015, and references therein) (0.4 1.6 kpc).The completeness curve shows a similar limiting magnitude but slightly different shape in the faint and bright ends (at most reaching a 12% difference in completeness).When re-obtaining the GSMF using the completeness curve derived for extended sources, all 1/Vmax points  but the lowest stellar mass point are within the errorbars (the lowest stellar mass 1/Vmax point is the least weighted one in the fitting due to larger uncertainty).In summary, considering that our galaxies are extended instead than point-like does not lead to any significant change in the derived GSMF Schechter parameter values and does not affect any of our conclusions.

Stellar mass completeness
Estimating the minimum stellar mass at which the sample is complete is crucial for studying the low-mass end of the GSMFs.In the past decades, different methods to estimate the mass completeness have been developed.In particular, for our analysis we follow the method proposed by Pozzetti et al. (2010), which is commonly used among the literature (e.g., Davidzon et al. 2017;Weaver et al. 2022).
First, galaxies with similar mass-to-light ratios (M/L) to the faintest ones are selected by excluding the 40% brightest galaxies of the total sample.Their stellar mass M lim at the limiting flux of the survey F lim stack , which in our case is the stacked flux used for the detection (Sect.3.1), is determined following the relation (Weigel et al. 2016): The value of limiting stellar mass at the desired completeness and redshift is defined as the corresponding percentile of the M lim (z) distribution.The GSMF fit can be considered as secure down to this stellar mass, but below M lim (z) the mass estimations have nonnegligible incompleteness.
Fig. 4 shows the M lim (z) at 90% completeness curve for PRIMER-UDS (red), HUDF (blue) and for the combination of both (black).The limiting flux of HUDF field is considerably smaller compared to PRIMER-UDS, which results in a lower limiting stellar mass at 90% completeness.This effect is higher at low redshift, where the difference is ∼ 1 dex and milder at high redshift with ∼ 0.5 dex.
In this study we probe values of stellar mass as low as 10 8 M ⊙ up to z ≃ 5.In other words, we are able to probe down to ≃ 0.5 dex less massive galaxies than the deepest data of CANDELS (Grazian et al. 2015) and also down to ≃ 1.75 dex compared to COS-MOS2015 (Davidzon et al. 2017).We probe down to similar mass ranges as Song et al. (2016) (that uses very deep HST and SPITZER data in CANDELS/GOODS-S and HUDF) and (Stefanon et al. 2021) (that uses very deep SPITZER data in CANDELS/GOODS-N, CANDELS/GOODS-S, COSMOS/UltraVISTA).In literature, there have been studies (e.g., Kikuchihara et al. 2020) that use gravitational lensing to detect extremely low-mass galaxies, their limiting mass being M ≃ 10 7 M ⊙ at z ≃ 6 − 7 and M ≃ 10 8 M ⊙ at z ≃ 9.

RESULTS
By following the method described in Sect.4, we parameterize the evolution of the GSMF over cosmic time.First, we have derived the 1/Vmax points with the method described in Sect.4.1.Then, we fitted a Schechter Function parametrized as in Eq. 4 to the data points by applying a reduced χ 2 fitting method, as described in Sect.4.1.However, to account for the biases produced by the propagation of uncertainties to the galaxy number counts (Eddington bias), the Schechter function is first convolved with the stellar mass error distribution (described in Sect.A.6), to obtain the bias-free Schechter parameters.There is no clear consensus on the correct Eddington bias estimation.We refer the reader to Davidzon et al. (2017) for a thorough discussion on the effects of employing different methods for assessing the Eddington bias.
Due to the nature of our dataset, the parameter that is best constrained is the low-mass end slope (α).Our JWST observations are extremely deep but their area is not large, ≃ 20 arcmin 2 in total.Characterizing the high-mass end is possible by using ancillary data from larger area surveys (e.g., CANDELS).We will focus our discussion on the evolution of the low stellar-mass end slope and cosmic stellar mass density (CSMD), but we also test for consistency and evolution of the other parameters, M * and ϕ * .However we refer the reader to Davidzon et al. (2017) and Weaver et al. (2022) for a more comprehensive analysis of these other parameters.We stress that a direct comparison with the literature might result in different absolute values for each Schechter parameter, the low-mass end slope value is coupled to the value of M * .This is why a comparison can be made using other quantities such as the (evolution of the) CSMD or the direct 1/Vmax realizations, as shown in Fig. 5. .Parameter space for the Schechter fits using the χ 2 method described in Section 4.1.The point of minimum χ 2 is marked with a color dot, and 1σ contours (obtained imposing ∆χ 2 red = 1) are shown with solid color lines.The M * parameter value has been fixed for the redshift bin z: 7.5 -8.5 (black square), therefore the corresponding parameter space values are not shown in this figure, except in the middle panel.

Evolution of the GSMF with cosmic time
The results for the GSMF are shown in Table 3 and Fig. 5. Hexagons and crosses are 1/Vmax points calculated from JWST +HST and CANDELS ancillary data, respectively.The continuous line is the χ 2 fit to the 1/Vmax points assuming a Schechter function parametrization, with parameters given in Table 3.The shaded region is the maximum span of 30 Monte Carlo realizations.For this we recompute GSMF after scrambling the photometry within the error bars and reobtaining stellar masses and photometric redshifts (see Sect.A.2).This can be understood as an upper limit for the dispersion of the points arising from photometric uncertainties and SED fitting.As shown in Caputi et al. (2015), it can be also used to asses the effect of the propagation of photometric uncertainties to the GSMF (Eddington bias).
We find an evolution of > 1σ significance for both α and ϕ * between each redshift bin.Our results show that the low-mass end slope steepens with cosmic time, following a consistent trend between z ≃ 7 and z ≃ 4. The value of α changes from close to α ≃ −2.0 for z ≃ 7 to α ≃ 1.6 for z ≃ 4.This trend holds for the intermediate redshift bins and is discussed in detail in Sect.6.We find 1σ agreement in the 1/Vmax points with respect previous studies (e.g., Duncan et al. 2014;Grazian et al. 2015;Song et al. 2016;Stefanon et al. 2021).
A number of studies, from Pérez- González et al. (2008) to more recent ones like Stefanon et al. (2021), find that the normalization factor ϕ * evolves with cosmic time following the cosmic stellar mass buildup.Grazian et al. (2015)  from z ≃ 7 to z ≃ 4, compatible within the previous measurements within the uncertainties.
For M * , values are compatible with each other within 1σ.We emphasize that probing M * is not the main goal of this study, and that the high-mass is taken from ancillary CANDELS data.Recent studies that probe the high-mass end of the GSMF conclude that there is not a clear evidence on the evolution of M * (e.g., Davidzon et al. 2017), however evidences for such a evolution have been found for lower redshift samples (e.g., Adams et al. 2021).We find that M * values are, in general, consistent with the ones derived by Davidzon et al. (2017) and Grazian et al. (2015).However, some differences are expected to arise from the different methodology used for Eddington bias correction.
Compared to more recent studies, the resulting values of the GSMF best fit parameters are compatible to the ones derived by Song et al. (2016), that also made use of CANDELS and HUDF fields for their study.Stefanon et al. (2021) sees a non-steepening α and consistently lower values of M * , although we stress that our constraint in M * is not strong because of the small area covered by our dataset and comes mainly from existing CANDELS data.We thoroughly study the evolution of α in Sect.6.
Table 4.The GSMF computed with the 1/Vmax method.The highest redshift bin z ≃ 8 is twice as large as the other redshift bins, in order to increase the statistics and, thus, robustness of our result.The points that lie under the stellar mass completeness limits or that have ≤ 1 galaxy counts are not quoted in this Table.As seen in Fig. 5, we probe the GSMF until z ≃ 8.5.Due to the few number of galaxies in each stellar mass bin, for both JWST and CANDELS data, the value of M * has been constrained between the range of values found for the other redshift bins in the last redshift bin z ≃ 7.5 − 8.5.This is why the value of α derived for this redshift bin (z ≃ 8) is only given as a tentative value, and should be investigated further in future observations.
In Fig. 5 we compare our results to the GSMF computed from IllustrisTNG (Pillepich et al. 2018, TNG100) as a grey dashed line.The agreement is good for z ≃ 4 and z ≃ 5. Nonetheless, the low stellar-mass end has a steeper slope for the TNG100 determination at z ≃ 4. The low-mass end slopes of TNG100 are similar to the ones from our analysis for all redshift bins, but for z ≃ 6 and for z ≃ 7, the normalization of the GSMF is significantly different.TNG100 predicts a lower density of galaxies at all stellar masses.The low-mass end only behaves differently at z ≃ 8, where our observations also suffer from larger uncertainties.
We also include the GSMF obtained from Delphi (Dark Matter and the emergence of galaxies in the epoch of reionization; Dayal et al. 2014Dayal et al. , 2022;;Mauerhofer & Dayal 2023).The overall agreement of Delphi with our realization of the GSMF is better than the one of TNG100 for z ≃ 6, z ≃ 7.This is not true for the lower and higher redshift bins, where we find disagreement with Delphi for both high and low-mass end, specially for z ≃ 8. From the lack of data in the high-mass end, TNG100 volume is not large enough for capturing the galaxies that conform the massive end of the GSMF, specially above M > 10 9 M ⊙ and z ≳ 6.
In Fig. 6, we show the parameter space for Fig. 5 fits of the 1/Vmax points.It depicts the best-fit value as solid points and the 1σ confidence level as continuous contours.The 1σ confidence level has been determined by imposing for ∆χ 2 red = 1.This shows that the evolution of α and ϕ * is significant within contiguous redshift bins, and that there exists a steepening of α between z ≃ 4 and z ≃ 7. It can also be seen that there is statistically significant evolution of ϕ * to higher values toward lower redshifts.
We performed a careful analysis of the sources of uncertainty that may affect the GSMF.A detailed explanation is given in Appendix A. We conclude that these errors have a minor impact on our results and do not siginficantly change any of our conclusions.We include the best fit Schechter parameters not accounting for the Eddington bias in Table 5.
Table 5. Schechter function best-fit parameters for the GSMF obtained using the χ 2 method described in Sect.4.1 without Eddington bias correction.M * for z ≃ 8 has been fixed to the value of the previous redshift bin z ≃ 7, thus errors are not quoted for this parameter.

EVOLUTION OF THE LOW-MASS END SLOPE WITH COSMIC TIME
A key goal of this work is to constrain the GSMF low-mass end slope at high z, so here we analyse its evolution with redshift.While a direct comparison of α can be affected by numerous systematic effects such as Eddington Bias estimation, it can be helpful to visualize the evolution of this parameter with cosmic time.And although α is coupled to the other Schechter parameters (mainly M * ), we find a reasonable agreement with the literature for the M * values.Fig. 7 shows the comparison between the value of α obtained in our analysis and the values fromliterature.Errors have been estimated from the χ 2 fit by taking into account all contributions described in Sect.9. Points from previous studies have been included as black empty or gray shaded symbols: González et al. (2011), Duncan et al. (2014), Caputi et al. (2015), Grazian et al. (2015), Mortlock et al. (2015), Song et al. (2016), Davidzon et al. (2017), Stefanon et al. (2021), Santini et al. (2022).Predictions from Illustris TNG100 and Delphi are shown as dashed and dotted grey lines, respectively.
It is clear that α is becoming steeper towards high redshifts.This implies that the evolution of the GSMF is not driven by a pure (number) density evolution.
The agreement with Song et al. (2016) is remarkably good for all redshift bins, except for the last one (z ≃ 8), where α = −2.25 ± 0.5 for Song et al. (2016) and α = −1.93±0.22 for our analysis.Our results also agree with the ones from Stefanon et al. (2021), where the biggest difference between our and their results occurs at z ≃ 7 where α = −1.73 ± 0.15 for Stefanon et al. (2021) and α = −1.98±0.14 for our analysis.However, these values are still compatible within 1σ in both cases.
We find a good agreement with Grazian et al. ( 2015) for all redshift bins except for z ≃ 6, where our value of α is significantly higher: α = −1.55 ± 0.19 for Grazian et al. (2015) and α = −1.88 ± 0.09 for our analysis.Note that the Grazian et al. (2015) limiting mass is M lim ≃ 10 9.2 M ⊙ while ours is M lim ≃ 10 8.4 M ⊙ , i.e. our constraint on the low-mass end slope α is more robust.Davidzon et al. (2017) find an overall higher value of α with respect to our analysis but also with respect to Grazian et al. (2015) and the rest of literature.This can be due to the large area of the COSMOS field (Scoville et al. 2007) but shallower data.The limiting stellar mass for z ≃ 5 is M lim ≃ 10 10 M ⊙ .The effect of the Eddington Bias estimation can also play a role in their estimations, as discussed in Davidzon et al. (2017).The α values obtained by Davidzon et al. (2017) are not included in the comparison in Fig. 7.
From our analysis it follows that the low-mass end slope steepens with cosmic time within z ≃ 4 − 7, from α = −1.61± 0.05 for z ≃ 4 to α = −1.98 ± 0.14 for z ≃ 7. We find a good agreement with recent studies that probe the low-mass end of the GSMF (e.g., Song et al. 2016;Bhatawdekar et al. 2019;Stefanon et al. 2021) in the steepening of α with redshift.We cannot confirm or rule out the presence of a turnover above z ≃ 8. Altough the 1/Vmax points for z ≃ 8 agree with other studies, we find different Schechter parameters.We consider the value of α at z ≃ 8 obtained from our analysis as tentative (we include it as empty symbols in Figs.7 and 8).
However, in Kikuchihara et al. ( 2020) there is no significant evolution of α between z ≃ 6 − 7 and z ≃ 9 for their study of the GSMF in the Hubble Frontier fields (HFF, Lotz et al. 2017).Bhatawdekar et al. (2019) do find a steepening of α (within the HFF cluster MACSJ0416.1-2403)between z ≃ 6 and z ≃ 9, this is, identical redshift range studied by Kikuchihara et al. (2020).Both of them probe a similar range of stellar masses, reaching M ≃ 10 8 M ⊙ for z ≃ 8, although arguably Kikuchihara et al. ( 2020) uses a larger dataset.
Finally, both TNG100 and Delphi predict a steepening low-mass end slope with cosmic time within z ≃ 4 (z ≃ 5 for Delphi ) and z ≃ 8.The values from Delphi are located within the error bars of our measurements for all redshifts.In the case of TNG100, we observe a rigid shift towards steeper values of α that holds roughly constant for all redshift bins analyzed herein.

EVOLUTION OF THE COSMIC STELLAR MASS DENSITY
The cosmic stellar mass density (CSMD, ρ ⋆ ) accounts for the total mass in form of stars that has been formed until a certain redshift or cosmic time.The CSMD provides a more global overview of the stellar mass buildup compared to the GSMF.We can derive the CSMD from the GSMF computed for a given redshift bin (with mean redshift z i ) by integrating the GSMF (Φ zi ) for reach redshift bin (z i ) over a range of stellar masses (M inf , M sup ) (Fontana et al. 2006): In the literature, the values of stellar masses to compute the integral are usually chosen as M inf = 10 8 M ⊙ and M sup = 10 13 M ⊙ (e.g., Elsner et al. 2008;Stark et al. 2013;Kikuchihara et al. 2020).This is the con-  2022)).We include the CSMD resulting from integrating the cosmic star formation history of Madau & Dickinson (2014) with a recycled fraction of 41% as a black solid line.We also compare our results with the predictions of theoretical models, namely Delphi (dotted black line), Illustris TNG100 (dashed black line) and Di Cesare et al. ( 2023) (dash-dotted black line).
vention followed in this work as well.Values for the CSMD are shown in Table 6.We used the GSMF bestfit parametrizations shown in Table 3 after correction for Eddington bias for calculating the CSMD.Fig. 8 shows the evolution of the CSMD from z ≃ 3.5 to z ≃ 8.5.We find an overall good agreement with previous studies.In the common redshifts, z ≃ 4 and z ≃ 5, our values are very similar to the ones of Davidzon et al. (2017).Song et al. (2016) cover the same redshift range as us, finding compatible results for all of them within the uncertainties, although we find an overall flatter evolution of the CSMD.
Our results on the CSMD are very similar to those by Bhatawdekar et al. (2019) for redshifts z ≃ 6, z ≃ 7 and z ≃ 8. Nevertheless, the values of the CSMD are more uncertain at higher redshifts.This is also the case for Stefanon et al. (2021), except for z ≃ 7 were their CSMD is more similar to the one by Song et al. (2016), but compatible with our estimation within 1σ.
All the CSMD values show the same behavior at z ≃ 8, but in our case incompleteness can be affecting the results.To compute the CSMD down to M inf = 10 8 M ⊙ we had to extrapolate below the mass completeness limit.Also, the value of M * was fixed for the best-fit Schechter function.For all the reasons above, our estimation of the GSMF and CSMD at z ≃ 8 has to be taken as a tentative one.
We have computed the CSMD integrating the cosmic star formation history by Madau & Dickinson (2014), considering a recycled fraction of R = 0.41 (41%), compatible with the Chabrier IMF Chabrier (2003).We can reproduce the flattening of the CSMD, but our point at z ≃ 8 falls below the prediction of Madau & Dickinson (2014).However, only very few measurements were available to constrain the CSMD at z ≃ 8 at the time of the Madau & Dickinson (2014) analysis.
Fig. 8 also shows the CSMD calculation for TNG100 simulation as a black dashed line, integrating over the GSMF presented in Fig. 5.The agreement is good for redshifts z ≃ 4 and z ≃ 5. From z ≃ 6, however, the CSMD of TNG100 decreases rapidly meanwhile the our results and the ones from literature suggest that the shape of the CSMD flattens with higher redshifts.
We have also included the recent results by Di Cesare et al. (2023) and Delphi in the comparison, as dotted and dashed-dotted black lines.The CSMD of Di Cesare et al. ( 2023) is significantly higher than the one of TNG100 at all redshifts, and the agreement with our results is better at high redshifts (z ≃ 6, z ≃ 7 and z ≃ 8).Di Cesare et al. ( 2023) make use of the hydrodynamical code DustyGadget (Graziani et al. 2020).It accounts for the dust produced by stellar populations and results in a better agreement with the observed CSMD obtained in this work.The CSMD predicted by Delphi is in good agreement with the predictions of (Di Cesare et al. 2023) for redshifts z ≳ 7.For lower redshifts, Delphi behaves more similarly to the best-fit line produced by Madau & Dickinson 2014 and our own data points.

SUMMARY AND CONCLUSIONS
In this paper we study the GSMF and its evolution with redshift between z ≃ 4 and z ≃ 8.By making use of new observations from JWST up to ∼ 5 µm in the PRIMER-UDS and HUDF fields, we are able to probe individual galaxies down to stellar masses M ∼ 10 8 M ⊙ up to z ≃ 8 for the first time.Until now, this low stellarmass regime at high redshifts was only reached using lensing clusters (e.g., Bhatawdekar et al. 2019;Rinaldi et al. 2022), or statistically with stacking analysis (e.g., Song et al. 2016), and in most cases the stellar-mass determinations were based on long spectral extrapolations.In contrast, the stellar masses derived here are based on the directly measured rest-frame optical light of galaxies up to the highest redshifts.To compute the GSMF, we adopt a technique based on the 1/Vmax method, using a χ 2 minimization method that samples the parameter space and provides us with both best-fit values and uncertainties for the GSMF Schechter-function parameters.Before doing so, we assess all the systematics and sources of uncertainty that can affect the GSMF, such as the propagation of photometric uncertainties to the photometric redshifts and stellar masses (the so-called Eddington bias) and cosmic variance.
To account for the Eddington Bias, we first derive the error distributions on the stellar mass for each redshift and stellar mass bin.We choose a combination of Gaussian and Student-T kernels, allowing for skewness in the final distribution.We then convolve the pure Schechter function with the error distributions to correct for the effects of Eddington bias.This affects the shape of the GSMF mainly in the high and low-mass ends.In particular, the value of the low-mass end slope becomes less steep after accounting for the Eddington bias, and M * moves towards lower masses.
The GSMF obtained here is in good agreement with the literature based on pre-JWST datasets (e.g., Grazian et al. 2015;Song et al. 2016;Bhatawdekar et al. 2019).This overall good agreement could be considered to some extent fortuitous, as our study is one of the first one to directly probe down to the typical stellar masses of (local) satellite galaxies at high redshifts with significant level of statistics.
Also, in agreement with some recent results (albeit not all), we find a significant (> 1σ) evolution of the faint-end slope α with cosmic time, corresponding to a steepening of the low-mass end of the GSMF toward earlier times.Instead, we do not see any significant evolution in M * .
When comparing to TNG100, we find that the agreement is good for the lowest redshifts analysed here (z ≃ 4 − 5).The faint-end slope is in broad agreement up to higher redshifts, but the normalization of TNG100 falls significantly below ours for at z > ∼ 5, which suggests that these models systematically underestimate the overall number density of galaxies at high redshifts.With respect to Delphi we find a overall better agreement specially in the higher redshifts z ≳ 5, in both the GSMF and the CSMD.
The number density of galaxies with masses between log 10 (M/ M ⊙ ) = 8.25 − 8.75 (the lowest mass regime within our completeness in estellar mass) grows by a factor of 1.5 between redshifts z ≃ 4 and z ≃ 7 according to our realization of the GSMF, meanwhile Delphi predicts a increase factor of 2.8.However, the number density of these galaxies in Delphi is systematically higher at the low redshifts z ≃ 4.5.
Finally, we compute the cosmic stellar mass density (CSMD) and study its evolution.We find a broad agreement with the recent literature, but at the same time predict a tentative shallower evolution from z ≃ 4 to z ≃ 7. Unfortunately, the large errors for the CSMD data point at z ≃ 8 do not allow us to confirm whether the flattening trend continues up to such high redshift.TNG100 shows a similar behaviour to our results for at z ≃ 4 − 5, however their CSMD falls rapidly below ours (and most literature) for z ≳ 6.New models such as those by Di Cesare et al. ( 2023) seem to reproduce better the evolution of the CSMD for these higher redshifts.Poisson statistics describe random processes with small number counts, which is the case of the 1/Vmax method used here to compute the GSMF for our galaxy samples.The error associated with having N counts in a stellar mass and redshift bin is √ N .For bins with no detections at all, upper and lower limits can be obtained from the Poisson distribution, as described in Gehrels (1986) and in Ebeling (2003) (where a more robust estimation of the upper limits is provided).However, we have opted not to include them in our calculation, as very few bins have no detections in both our data and the ancillary CANDELS dataset.
Figure 9. Left: the bottom panel shows the effect of Eddington Bias correction on the Schechter fit for the GSMF, note the effect on both the high and low-mass end by making the GSMF less steep and reducing the value of M * .The middle panel shows the 30 realizations of the GSMF after randomizing the photometry within the errors, where the scatter is higher for the lower and highest mass points.The upper panel shows the effect of accounting for nebular emission, all the 1/Vmax points are in good agreement within 1σ errors.The lowest mass point showing the biggest difference.Right: the effect of cosmic variance for redshift bins z ≃ 4, 5, 6, 7, by plotting the 1/Vmax for PRIMER-UDS (red) and HUDF (blue) separately.

A.2. Photometric uncertainty
Photometric uncertainty is the one of the main random uncertainties contributors to the error budget.Errors in photometry propagate into the stellar mass and photometric redshift obtained from the SED fitting.This errors are especially relevant for faint galaxies, and their correct assessment is crucial for an appropriate description of the GSMF low-mass end.
In order to estimate the effect of this uncertainties on the GSMF, we have performed 30 Monte Carlo realizations.First, we have scrambled the flux in each band assuming a normal distribution given by the errors in photometry.Then, photometric redshifts and stellar masses have been extracted for each source from LePHARE.The methodology is the same followed for estimating the Eddington Bias, as explained in Sect.A.6.
Having done this, we then calculated a set of 1/Vmax points following the methodology described in Sect.4.1.In panel 2 of Fig. 9 we show the scatter of the 1/Vmax points, with the 30 Monte Carlo 1/Vmax points shown in gray and the real value superimposed in green.The effect is greater toward stellar masses M < ∼ 10 8.5 M ⊙ , with a scatter of 0.4 − 0.5 dex.It becomes milder for galaxies with stellar masses M ≳ 10 8.5 M ⊙ , with typical scatters of 0.1 − 0.2 dex.

A.3. SED fitting templates
Systematic uncertainties can additionally be introduced by the templates used to fit the observed photometry and star formation histories.In Marchesini et al. (2009) a thorough analysis of the effect of uncertainties on the GSMF is made.It is found that testing a set of different set of stellar templates, IMFs and metallicities can have a non negligible effect on the realization of the GSMF.
In the present work we will only take into account random uncertainties.The high quality of the photometry we use together with the wide spectral coverage of HST and JWST minimize the offsets in photometric redshift and stellar mass even when using different templates for SED fitting.

A.4. Nebular emission
As described in Sect.3, we assesed the effect of nebular emission on the GSMF by including SB99 stellar templates (Leitherer et al. 1999) which account for young (age< 10 Myr) and star forming (with constant star forming history) galaxies.The final masses and photometric redshifts have been obtained from the LePHARE run between BC03 and SB99 that has a smaller reduced χ 2 .
In our case, about 15% and 18% of the total number of galaxies have lower reduced χ 2 when using SB99 models for HUDF and PRIMER-UDS, respectively.To test the significance of this result to 1σ, we impose that the difference of reduced χ 2 is greater than 1 (∆χ 2 red > 1).The amount of galaxies that are better fit by SB99 templates with 1σ significance is of 4% and 1% respectively.The effect on including nebular emission is shown in panel 3 of Fig. 9.The scatter produced by using different templates is within the error bars for all masses and redshifts, so we do not expect the inclusion of nebular lines (following this methodology) to significantly alter the shape of the GSMF.
As studied in Rinaldi et al. (2023), strong line emitters can show a flux excess in their photometry, particularly in medium bands.For assessing the effect of line contamination in the estimation of stellar masses, we re-run the SED fitting removing the possible medium bands affected by typically strong lines or lines complexes (i.e., Hα and Hβ + OIII) on a galaxy-by-galaxy basis.After doing so, we find that the fraction of galaxies that show a difference on stellar mass greater than 0.3dex is around 12%, with a mean offset of around 0.1dex.To all effects, within the where σ and df are the dispersion parameters of the Gaussian and Student-T error distributions, respectively, and ∆ is the mean of the Gaussian component, allowing for asymmetry in the error distributions.
The main reasons that have motivated our reparametrization of the kernel are that we allow for non-symmetrical distributions (which can be relevant for the lowest stellar mass bins), and that the Student-T distribution provides a better fit to the actual errors of our data.In addition, we carefully study the dependance of the errors with both stellar mass and redshift.
The method that we developed can be summarized in the following steps: 1.We divided the sample into stellar mass and redshift bins.2. We obtained 30 realizations of stellar mass and photometric redshift after randomizing the photometry within the error bars using LePHARE with the same configuration as the original fitting.3. We derived the stellar mass error distribution for each photometric redshift and stellar mass bin by studying the differences between the original masses and the 30 LePHARE realizations.4. We fitted the Gaussian × Student-T kernel K(σ, df, ∆) to the error distribution.Here df are the degrees of freedom of the Student-T distribution and ∆ is the skewness of Gaussian component of the error distribution.Panel 1 of Fig. 9 shows the Eddington bias correction effects on the shape of the GSMF for the redshift range z : 3.5 ≤ z ≤ 4.5.It can be seen that our approach produces a significant effect in both low and high-mass end of the GSMF, thus altering the values of α and M * .Fig. 11 shows the error distributions for three stellar-mass ranges obtained by considering 30 realizations of photometric redshift and stellar mass from LePHARE using as input scrambled photometry within the error bars.In the plot we only show the distributions for 3.5 ≤ z ≤ 4.5, but for both fields HUDF (red) and PRIMER-UDS (blue).The Gaussian times Student-T kernel fittings are shown in dotted and dashed black lines, respectively.
Our error distributions are skewed towards lower masses.Such an effect has been also found in Grazian et al. (2015).We also show that the deeper data in HUDF benefits from smaller photometric uncertainties, which results in an overall smaller scatter in stellar mass and less skewness in the error distribution.Fig. 11 shows the error distributions for the combined HUDF + PRIMER-UDS sample.We conclude that a clear evolution of the error distribution with both stellar mass and redshift is present for the total sample.However, combining both fields results in a less skewed error distribution, especially compared to the one from PRIMER-UDS (this is a shallower dataset with larger photometric uncertainties at any given stellar mass).

Figure 1 .
Figure 1.Top: PRIMER-UDS field.Bottom: JWST HUDF field.Sources are marked in black and background in red.Lighter shades indicate deeper areas.The XDF can be seen as the lightest patch on the left-hand side of the field.mass M * .Therefore, we complement our galaxy catalogs in the JWST -covered fields with the publicly available source catalogs from the Cosmic Assembly Near-Infrared Deep Extragalactic Legacy Survey (CAN-DELS; Grogin et al. 2011; Koekemoer et al. 2011) in the GOODS-S field (≃ 170 arcmin 2 ) (Galametz et al. 2013) and the UDS field (≃ 200 arcmin 2 ) (Guo et al. 2013).The larger area covered by CANDELS allows us to design a tiered strategy, with very deep JWST observations to probe the faint end (in small areas) and wider-area data to constrain the GSMF at the high-mass end.2.1.1.HUDFJWST data for HUDF (including XDF) have been obtained byWilliams et al. (2023), as part of a General Observers Cycle-1 program (Proposal-ID: 1963; PI: Christina C. Williams).The obtained datasets consist of five JWST -NIRCam medium bands: F182M, F210M, F430M, F460M, and F480M.We also make use of NIR-

Figure 2 .
Figure 2. Stellar mass and photometric redshift distribution of the galaxy sample used for our study.Catastrophic outliers are shown in red.The total number of galaxies in the sample is 988, out of which 12.3% are catastrophic outliers.The median absolute deviation is 0.0419, or 0.0356 after removing the catastrophic outliers from the sample.

Figure 3 .
Figure 3. Stellar mass and photometric redshift distribution of the galaxy sample used for our study.UDS and HUDF are shown as red and blue dots, respectively.From the plot, it can be seen that the difference between the limiting mass of HUDF and UDS is 1 dex or larger, consistent with the deeper data that we use in HUDF.

Figure 4 .
Figure 4. Left.Completeness fraction versus stacked magnitude, as obtained from introducing simulated point sources to the detection images.PRIMER-UDS is shown in red and HUDF in blue.Empty symbols represent the shallow part and filled ones the deep part for both fields.The horizontal dashed line shows 65% completeness.Right.Limiting mass at 90% completeness over the two fields HUDF and UDS using the method described in Pozzetti et al. (2010).The value of M lim is shown as points for all galaxies excluding the brightest 40% of the sample.The M lim (z) curve is shown as a line.PRIMER-UDS is shown in red, HUDF in blue and the combination of them in black.The shaded area represents the redshift range z ≃ 4 − 8.
Figure6.Parameter space for the Schechter fits using the χ 2 method described in Section 4.1.The point of minimum χ 2 is marked with a color dot, and 1σ contours (obtained imposing ∆χ 2 red = 1) are shown with solid color lines.The M * parameter value has been fixed for the redshift bin z: 7.5 -8.5 (black square), therefore the corresponding parameter space values are not shown in this figure, except in the middle panel.

Figure 7 .
Figure 7. Evolution of the faint-end slope as a function of redshift.Points from our study are shown as red hexagons.Errors have been estimated from the χ 2 fit by taking into account all contributions described in Sect.9. Points from previous studies have been included as black empty or gray shaded symbols:González et al. (2011),Duncan et al. (2014),Caputi et al. (2015),Grazian et al. (2015),Mortlock et al. (2015),Song et al. (2016),Davidzon et al. (2017),Stefanon et al. (2021),Santini et al. (2022).Predictions from Illustris TNG100 and Delphi are shown as dashed and dotted grey lines, respectively.
like to thank P. Cáceres-Burgos and J. Méndez-Gallego for useful scientific discussions; and to Pratika Dayal and Valentin Mauerhofer for providing numerical predictions of the DELPHI galaxy formation model.RN, VK and KIC acknowledge funding from the Dutch Research Council (NWO) through the award of the Vici Grant VI.C.212.036.This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope.The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST.These observations are associated with JWST programs GO #1963, GO #1895, GO #1837.The authors acknowledge the team led by coPIs C. Williams, M. Maseda and S. Tacchella, and PI P. Oesch, for developing their respective observing programs with a zero-exclusive-access period.Also based on observations made with the NASA/ESA Hubble Space Telescope obtained from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 526555.Software: AstroPy (Collaboration et al. 2022), dustmaps (Green 2018), LePHARE (Arnouts & Ilbert 2011), Matplotlib (Hunter 2007), NumPy (Harris et al. 2020), Photutils (Bradley et al. 2022), Python (van Rossum 1995), SciPy (Virtanen et al. 2020), Source Extractor (Bertin & Arnouts 1996), TOPCAT (Taylor 2017), WebbPSF (Perrin et al. 2014).Facilities: HST, JWST APPENDIX A. UNCERTAINTIES IN THE DETERMINATION OF THE GSMF A.1.Poisson uncertainty

Figure 11 .
Figure11.Left: Stellar mass error distributions for 3.5 ≤ z ≤ 4.5 galaxies, obtained after 30 realizations of stellar mass and photometric redshifts using randomized photometry.The actual error distributions for HUDF (red) and PRIMER-UDS (blue) are shown.The fits of our kernel are shown as dotted and dashed black lines, respectively.Vertical grey bands show the central value of each stellar mass bin, namely log10(M/ M⊙) : 7.5, 8.5, 9.5.Right: Stellar mass error distributions for the combined PRIMER-UDS + HUDF fields.Each color shows a different redshift bin, z : 4, 5, 6.The grey bands show the central value of each mass interval (log10(M/ M⊙) = (7.5, 8.5, 9.5)).A clear evolution of the errors with stellar mass and redshift is present in our data.

Table 1 .
Area and depth (5σ limiting magnitudes measured in apertures of 0. ′′ 25/pix radius) for JWST fields used for our analysis.The depths in JWST bands is uniform in the case of the HUDF, while the HST data has different depths in the XDF and the rest of HUDF regions.

Table 2 .
SED parameters employed for the LePHARE run, both for BC03 and SB99 stellar population models used in this work.

Table 3 .
Schechter function best-fit parameters for the GSMF obtained by applying the χ 2 method described in Sect.4.1 and corrected for Eddington bias (following the method described in Sect.A.6). M * for z ≃ 8 has been fixed to the value of the previous redshift bin z ≃ 7, thus errors are not quoted for this parameter.

Table 6 .
Results for the CSMD obtained from our best fit Schechter functions, integrating from M = 10 8 M⊙ to M = 10 13 M⊙.