# The Extreme Faint End of the UV Luminosity Function at $z \sim 6$ Through Gravitational Telescopes: a comprehensive assessment of strong lensing uncertainties

Hakim Atek, Johan Richard, Jean-Paul Kneib, Daniel Schaerer

► To cite this version:

Hakim Atek, Johan Richard, Jean-Paul Kneib, Daniel Schaerer. The Extreme Faint End of the UV Luminosity Function at $z \sim 6$ Through Gravitational Telescopes: a comprehensive assessment of strong lensing uncertainties. Monthly Notices of the Royal Astronomical Society, 2018, 479 (4), pp.5184-5195. 10.1093/mnras/sty1820. hal-01757202

HAL Id: hal-01757202

https://hal.science/hal-01757202v1

Submitted on 28 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The extreme faint end of the UV luminosity function at z ~ 6 through gravitational telescopes: a comprehensive assessment of strong lensing uncertainties

Hakim Atek,<fnref n="1" />★ Johan Richard,<fnref n="2" /> Jean-Paul Kneib<fnref n="3" />,<fnref n="4" /> and Daniel Schaerer<fnref n="5" />

Accepted 2018 July 5. Received 2018 July 4; in original form 2018 March 21

## ABSTRACT

With the Hubble Frontier Fields program, gravitational lensing has provided a powerful way to extend the study of the ultraviolet luminosity function (LF) of galaxies at $z \sim 6$ down to unprecedented magnitude limits. At the same time, significant discrepancies between different studies were found at the very faint end of the LF. In an attempt to understand such disagreements, we present a comprehensive assessment of the uncertainties associated with the lensing models and the size distribution of galaxies. We use end-to-end simulations from the source plane to the final LF that account for all lensing effects and systematic uncertainties by comparing several mass models. In addition to the size distribution, the choice of lens model leads to large differences at magnitudes fainter than $M_{UV} = -15$ AB mag, where the magnification factor becomes highly uncertain. We perform Markov Chain Monte Carlo (MCMC) simulations that include all these uncertainties at the individual galaxy level to compute the final LF, allowing, in particular, a crossover between magnitude bins. The best LF fit, using a modified Schechter function that allows for a turnover at faint magnitudes, gives a faint-end slope of $\alpha = -2.01_{-0.14}^{+0.12}$, a curvature parameter of $\beta = 0.48_{-0.25}^{+0.49}$, and a turnover magnitude of $M_T = -14.93_{-0.52}^{+0.61}$. Most importantly, our procedure shows that robust constraints on the LF at magnitudes fainter than $M_{UV} = -15$ AB remain unrealistic, as the 95 per cent confidence interval accommodates both a turnover and a steep faint-end slope. More accurate lens modeling and future observations of lensing clusters with the James Webb Space Telescope can reliably extend the ultraviolet (UV) LF to fainter magnitudes.

**Key words:** galaxies: evolution – galaxies: high-redshift – galaxies: luminosity function – gravitational lensing: strong.

## 1 INTRODUCTION

The epoch of cosmic reionization around $z \sim 6–10$ (e.g. Fan et al. 2006; Planck Collaboration XLVII 2016) has seen the neutral hydrogen content of the Universe become ionized, ending the period of dark ages. The identification of sources responsible for such dramatic phase transition remains a major question in extragalactic astronomy. Early star-forming galaxies could well be the best candidates to drive reionization (Robertson et al. 2013; Bouwens et al. 2015; Finkelstein et al. 2015) as the study of their luminosity function at $z > 6$ reveals a steep faint-end slope that indicates

a significant contribution from faint galaxies (Bunker et al. 2010; Oesch et al. 2010; Bouwens et al. 2011; McLure et al. 2013). Much of the progress in the field was achieved in deep blank fields observed with the Hubble Space Telescope (HST), which now reaches an observational limit around an absolute magnitude of $M_{UV} \sim -17$ AB. The total ultraviolet (UV) luminosity density emitted by galaxies down to such limit is insufficient to drive reionization and an extrapolation of the LF is used to account for the contribution of faint galaxies. However, the abundance of this population of faint galaxies strongly depends on the faint-end slope and the magnitude cut-off of the UV LF, which are uncertain.

Massive galaxy clusters act like cosmic telescopes, magnifying background galaxies in the strong gravitational lensing regime, and as such offer a viable route to reach beyond the current observational

<footnote n="1">¹Institut d'astrophysique de Paris, CNRS UMR7095, Sorbonne Université, 98bis Boulevard Arago, F-75014, Paris, France</footnote>

<footnote n="2">²Univ Lyon, Univ Lyon1, Ens de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69230, Saint-Genis-Laval, France</footnote>

<footnote n="3">³Laboratoire d'Astrophysique, Ecole Polytechnique Fédérale de Lausanne, Observatoire de Sauverny, CH-1290 Versoix, Switzerland</footnote>

<footnote n="4">⁴Aix Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille) UMR 7326, F-13388, Marseille, France</footnote>

<footnote n="5">⁵Observatoire de Genève, Université de Genève, 51 Ch. des Maillettes, CH-1290, Versoix, Switzerland</footnote>

<footnote>* E-mail: hakim.atek@iap.fr</footnote>
limit (Maizy et al. 2010; Kneib & Natarajan 2011; Postman et al. 2012; Sharon et al. 2012; Richard et al. 2014b; Coe, Bradley & Zitrin 2015). The Hubble Frontier Fields (HFF) program has delivered the deepest observations of lensing clusters to date, with 840 HST orbits reaching ~29 AB mag limit in seven optical and near-IR bands (Lotz et al. 2017), in addition to supporting space Spitzer, Chandra, and ground-based (VLT, ALMA) observations. As part of the HFF project, several groups produced mass models for each cluster, which were made public in order for the community to interpret high-redshift observations (e.g. Jauzac et al. 2014; Johnson et al. 2014; Grillo et al. 2015; Ishigaki et al. 2015; Hoag et al. 2016; Mahler et al. 2018).

The first studies of $z > 6$ UV luminosity functions relying on HFF data and these lensing models already discovered the faintest high-redshift galaxies (Atek et al. 2014; Yue et al. 2014; Zheng et al. 2014; Ishigaki et al. 2015), extending the LF by 2 mag. Combining more lensing clusters, later studies pushed the detection limits down to $M_{UV} \sim -13$ AB (Castellano et al. 2016; Kawamata et al. 2016; Laporte et al. 2016; Bouwens et al. 2017b and Ishigaki et al. 2018, hereafter, B17 and I18, respectively;) and even to $M_{UV} \sim -12$ AB (Livermore, Finkelstein & Lotz 2017, hereafter, L17). While luminosity functions from different groups tend to agree well on the global shape of the LF at $M_{UV} < -17$ mag, significant discrepancies appear at the very faint end, a territory where only sources with high magnification can be selected. For example, the most recent work by L17 shows significantly higher values of the LF at $M_{UV} > -17$ compared to results of B17, and I18 find a steeper faint-end slope compared to L17 and B17. While B17 find signs of a turnover in the LF at $M_{UV} > -15$ and a faint-end slope of $\alpha = -1.91 \pm 0.04$, L17 report a steep faint-end slope of $\alpha - 2.10 \pm 0.03$ and strong evidence against a possible turnover at $M_{UV} < -12.5$. This example highlights the challenges encountered by this type of studies to constrain the extreme faint-end of the UV LF at $z \sim 6$. The observed discrepancies in those studies might result from the differences in dropout catalogue selection, the size distribution adopted in the completeness simulations to compute the survey volumes (B17), and also from the different approaches in the treatment of lensing uncertainties.

In this paper we present for the first time a comprehensive assessment of the lensing uncertainties affecting the UV LF using end-to-end simulations from the source plane to the final UV LF. The focus of this study is the faint-end (beyond $M_{UV} = -15$) of the LF and how the choice of a given size distribution and a lensing model can impact the result. The paper is organized as follows. In Section 2, we describe the full HFF imaging data set used in the study and in Section 3 we present the selected sample of $z \sim 6$ galaxies. The lensing models are described in Section 4. Our end-to-end simulation procedure, used to both compute the effective survey volume (Section 6) and quantify the uncertainties, is detailed in Section 5. The final UV LF and the associated uncertainties are discussed in Section 7, while the UV luminosity density is computed in Section 8 before the conclusion given in Section 9. Throughout the paper we adopt standard cosmological parameters: $H_0 = 71$ km s$^{-1}$ Mpc$^{-1}$, $\Omega_{\Lambda} = 0.73$, and $\Omega_m = 0.27$. All magnitudes are expressed in the AB system (Oke & Gunn 1983).

## 2 HFF OBSERVATIONS

The Frontier Field program obtained deep HST imaging data in the optical and near-infrared for a total of six lensing clusters and flanking fields. The Advanced Camera for survey (ACS) was used for the optical coverage with three broad-band filters F435W, F606W,

and F814W. The Wide Field Camera Three (WFC3) was used for the NIR observations with four filters F105W, F125W, F140W, and F160W. HST obtains WFC3/IR and ACS optical simultaneously for each pair cluster-parallel field in the first epoch before switching instrument positions in the second epoch. Observations took place between 2013 and 2016 over HST cycles 21 to 23 devoting 140 orbits to each cluster/parallel pair. In Table 1 we give the limiting magnitude reached in each filter and measured in a circular aperture with 0.4"diameter (cf. Table 1).

In this analysis we use the high-level science products made available by the HFF data reduction team at the Space Telescope Science Institute<fnref n="1" /> (STScI). Data products include drizzled images, and weight maps from the FF program but also ancillary data from other observing programs. We used the ACS mosaics generated using the ‘self-calibration’ method and the WFC3/IR mosaics that were corrected for time-variable sky background. The detailed reduction procedure and science products can be found in Lotz et al. (2017).

All the images were first matched to the same Point Spread Function (PSF) with the F160W frame as reference. We used a PSF model computed for the F160W image with the TINYTIM software (Krist, Hook & Stoehr 2011). We created deep image stacks using a weighted combination of the four filters in the IR and the two bluest filters in the optical for the detection and non-detection, respectively, criteria for $z \sim 6$ dropout selection. For $z \sim 8$ galaxies the combined images include the three reddest IR filters for the detection and the three ACS filters for the non-detection criteria. The cluster fields are affected by intra-cluster light (ICL) and the bright cluster galaxies (BCGs) that significantly reduce our ability to detect galaxies around the cluster core. In order to mitigate these effects and improve the effective depth of the cluster fields, we performed a median filtering procedure to subtract the foreground light. Multiple filter sizes were used until we adopted the best trade-off between source recovery and the appearance of artefacts due to an over-subtraction of bright galaxies. We adopted a median filter of 2 arcsec $\times$ 2 arcsec for all individual and stacked images. On average, this procedure allows us to detect an additional five galaxies per cluster in the contaminated region. As described in Atek et al. (2015b), we ran simulations with mock galaxies exploring different sizes and profiles in order to compare the input and the computed magnitudes in the corrected images. We find that the flux errors become only significant for very faint galaxies ($H_{140} > 28AB$) where they reach 0.5 mag. Such uncertainties are comparable to the results obtained with more sophisticated techniques such as wavelet decomposition (L17). Other techniques have also been explored using GALFIT (Peng et al. 2002) combined with median filtering (B17), and other iterative modeling techniques of the BCGs (Shipley et al. 2018).

## 3 THE $z \sim 6$ GALAXY SAMPLE

We perform the source extraction using the SExTRACTOR software (Bertin & Arnouts 1996). In the cluster fields, the ICL-corrected images are used as a detection image while the photometry is performed in the original ones. The SExTRACTOR parameters were chosen to improve the detection of the faintest sources in the field while preserving the input flux used in the simulations. The following parameters were adopted: detect minarea=2 and detect thresh=1.5 for the detection and deblend nthresh=16

<footnote n="1">¹http://www.stsci.edu/hst/campaigns/frontier-fields/</footnote>
<table><caption>**Table 1.** Magnitude limits of HST observations in the six HFF clusters. The depth of the images are 3σ magnitude limits measured in a 0.4 arcsec circular aperture.</caption><thead><tr><th>Field</th><th>RA<br>J2000</th><th>Dec.<br>J2000</th><th>F435W</th><th>ACS<br>F606W</th><th>F814W</th><th>F105W</th><th>WFC3<br>F125W</th><th>F140W</th><th>F160W</th></tr></thead><tbody><tr><td>A2744</td><td>00:14:21.2</td><td>-30:23:50.1</td><td>28.8</td><td>29.4</td><td>29.4</td><td>28.6</td><td>28.6</td><td>29.1</td><td>28.3</td></tr><tr><td>MACS0416</td><td>04:16:08.9</td><td>-24:04:28.7</td><td>30.1</td><td>29.1</td><td>29.2</td><td>29.2</td><td>28.8</td><td>28.8</td><td>29.1</td></tr><tr><td>MACS0717</td><td>07:17:34.0</td><td>+37:44:49.0</td><td>29.5</td><td>28.6</td><td>29.3</td><td>28.9</td><td>28.6</td><td>28.5</td><td>28.8</td></tr><tr><td>MACS1149</td><td>11:49:36.3</td><td>+22:23:58.1</td><td>28.6</td><td>28.6</td><td>28.6</td><td>28.9</td><td>29.3</td><td>29.2</td><td>30.1</td></tr><tr><td>AS1063</td><td>22:48:44.4</td><td>-44:31:48.5</td><td>30.1</td><td>29.1</td><td>29.3</td><td>29.0</td><td>28.7</td><td>28.5</td><td>28.8</td></tr><tr><td>A370</td><td>02:39:52.9</td><td>-01:34:36.5</td><td>30.1</td><td>29.1</td><td>29.3</td><td>29.0</td><td>28.7</td><td>28.5</td><td>28.4</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>**Figure 1.** The median filtering procedure applied to MACS1149. The left frame is the F140W original image while the right frame shows the same image after subtracting the median filtered image described in Section 2 to correct for light contamination from bright cluster members. The same scaling is applied to both images.</figcaption></figure>

for source deblending. The isophotal magnitudes MAG\_ISO are used for the colour-colour selection whereas MAG\_AUTO values are adopted in calculating the total magnitudes. We adopt Lyman break selection criteria (e.g. Steidel et al. 1996) similar to those in previous studies of dropout galaxies at these redshifts (Bouwens et al. 2014; Finkelstein et al. 2015), and which were explained in detail in Atek et al. (2015a):

$$ 
\begin{align} 
(I_{814}-Y_{105}) &> 1.0 \nonumber \\ 
(I_{814}-Y_{105}) &> 0.6 + 2.0(Y_{105}-J_{125}) \tag{1} \\ 
(Y_{105}-J_{125}) &< 0.8 \nonumber 
\end{align} 
$$

Galaxies in the sample must also satisfy a detection significance above 5σ in at least two IR filters and 6.5σ in the IR stacked image. We also require a non-detection in both F435W and F606W filters and their stack. This colour selection includes both z ~ 6 and ~7 galaxies and is chosen to exclude low-redshift sources that have similar colours to those of high-z galaxies. The most important source of contamination consists of dust-obscured and evolved galaxies at lower redshift with extremely red colours and low-mass stars. In order to mitigate this contamination, we investigated the evolution of low-redshift elliptical galaxy templates from Coleman, Wu & Weedman (1980) and starburst galaxy templates from Kinney et al. (1996) in the colour-colour diagram as a function of redshift and extinction (see Fig. 2). We also simulated the colour track of stars from Chabrier et al. (2000) templates. An additional visual inspection is performed to rule out spurious detections and point-like sources. A more detailed review of all the potential interlopers is given in Atek et al. (2015a). As a natural consequence of strong lensing, we also expect that some of the background sources will have multiple images that need to be accounted for in the final galaxy number counts. The details of this procedure are given in Section 5.1. The final sample contains a total of 300 galaxies in all the cluster fields, with magnification factors ranging from μ ~ 1

<figure><img src="image_4.png" /><figcaption>**Figure 2.** The colour-colour selection for z ~ 6 dropout candidates. Sources that satisfy the selection criteria of equation (1) are shown with green circles. The expected colour tracks of low-redshift elliptical galaxies generated from Coleman et al. (1980) and starburst galaxies from Kinney et al. (1996) templates are shown in dotted and solid lines, respectively. The three colours represent the impact of the dust extinction in steps of A<sub>V</sub> = 1. The magenta points show the colour track of stars generated from Chabrier et al. (2000) templates.</figcaption></figure>

to few hundreds, including all models used here. The distribution of the observed apparent magnitudes (uncorrected for lensing) is given in Fig. 3.

## 4 LENSING MODELS

Computing an accurate luminosity function at z = 6–7 to the faintest luminosities requires a good knowledge of the lensing power of the HFF clusters. Lensing models are essential not only to estimate the magnification of the sources but also the effective survey volume and the completeness function (cf. Section 6). For the six clusters we adopt the publicly available lensing models (version 4.0) constructed by the CATS team (clusters as telescopes) using version 7.0 of the LENSSTOOL<fnref n="2" /> software, which follows a parametric approach to map the mass distribution in the cluster (Kneib 1993; Jullo et al. 2007; Jullo & Kneib 2009). The strong lensing analysis of the six HFF clusters using the full depth of HST observations is now complete. The mass reconstructions of the first three clusters were published in Jauzac et al. (2014, 2015); Limousin et al. (2016); Lagattuta et al. (2017). All the lens models are made publicly available

<footnote n="2">²http://projets.lam.fr/projects/lenstool/wiki</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Observed magnitudes distribution in the rest-frame UV at 1500 Å for the $z \sim 6$ dropout sample.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Cumulative surface area as a function of magnification. The area for each cluster is derived using the CATS model and each coloured area is indicative of the model uncertainties following the colour code in the legend. The magnification factor $\mu$ is given in units of magnitude.</figcaption></figure>

through the public Frontier Fields repository,<fnref n="3" /> which includes mass models submitted by several teams: the GLAFIC team (Kawamata et al. 2016, 2018), Bradac et al. (Bradač et al. 2005; Hoag et al. 2016), Merten & Zitrin (Merten et al. 2011; Zitrin et al. 2009, 2013), Sharon & Johnson (Johnson et al. 2014), and Diego et al. (2015). We note that at the time of conducting the present analysis, several teams have not updated their mass reconstruction models and do not include the full HFF observations. The impact of model uncertainties and systematic differences between the models on the computed UV LF is assessed in Section 7.3.

# 5 END-TO-END SIMULATION PROCEDURE

Alongside the galaxy number counts and their intrinsic luminosities, we need to estimate the effective survey volume in each cluster. Following Atek et al. (2014), we performed all the completeness estimates directly in the source plane. This approach has the advantage of naturally accounting for all the lensing effects that affect the galaxy properties (shape, luminosity) and the survey area. In Sections 7.1 and 7.2, we use the same technique in our end-to-end simulations to assess the reliability of the lensing models and size distribution in recovering the true UV luminosity function. The flow chart of Fig. 5 illustrates the different steps throughout the UV LF determination procedure.

## 5.1 Observed galaxies

The first step is the selection of dropout galaxies in the actual HST images (cf. Section 3) before de-lensing them to the source plane and determine their amplification factor using LENSTOOL. The magnification factor is computed accounting for the contribution of all the mass potentials of the cluster. Then the sources are lensed back to the image plane where we predict the position of the counter-images of each source in the multiple-image area. Then around this position we search for dropout galaxies that have a similar redshift probability, colour, and morphological symmetries in the case of resolved sources. We derive the intrinsic absolute magnitude based on the observed magnitude and the amplification factor. The galaxy redshifts are determined from the peak probability of the photometric redshifts using a modified version of the Hyperz code (Bolzonella, Miralles & Pelló 2000; Schaerer & de Barros 2009), which includes nebular continuum and emission lines contribution. Finally, the galaxy number counts are computed in magnitudes bins with a bin size of $\Delta\text{mag} = 0.5$ and corrected for the multiple images. As we will see in the assessment of lensing uncertainties of Section 7.2, we can use any lensing model to compute the intrinsic magnitudes and identify the multiple images.

## 5.2 Simulated galaxies

The second step consists of end-to-end simulations to determine the incompleteness function and the associated uncertainties. For each field, we simulated a set of 10 000 galaxies with randomly distributed redshifts and intrinsic absolute magnitudes and two different light profiles: an exponential disc and de Vaucouleurs profile (Ferguson et al. 2004; Hathi et al. 2008). The input galaxy sizes follow a log-normal distribution with a mean half light radius of 0.15 arcsec and a dispersion of 0.05 arcsec (Bouwens et al. 2004; Hathi et al. 2008; Grazian et al. 2012; Huang et al. 2013; Ono et al. 2013). However, the sizes are also luminosity-dependent following recent results on size-luminosity relation of high-redshift galaxies (Mosleh et al. 2012; Oesch et al. 2014; Kawamata et al. 2015). We applied a size-luminosity relation in the form of $r_{hl} \propto L^{0.5}$, derived for lensed galaxies (Kawamata et al. 2015; Bouwens et al. 2016). We discuss in further detail the influence of the size distribution in Section 7.1.

The simulated galaxies are randomly distributed directly in the source plane of each cluster, which was reconstructed from the WFC3 field of view using the lens model. The total survey area in the source plane and its evolution as a function of magnification is

<footnote n="3">³https://archive.stsci.edu/prepds/frontier/lensmodels/</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Illustration of the source plane method in computing the UV luminosity function and assessing the impact of lensing model uncertainties. Magnification maps provide the amplification of the dropout galaxies identified in the HST images, which are de-lensed into the source plane using the deflection maps and then projected back into the image plane to identify multiple-image systems. In parallel, we create in the source plane a uniform distribution of simulated sources, which are lensed to the image plane using the same mass model as the actual observations. The process includes the amplification and shape distortion of the galaxies and it automatically accounts for multiple images. Here we show directly the real WFC3 image (of MACS0717 in this example) where the galaxies are implemented. Note the rarity of sources in the region of highest magnification around the critical line at $z \sim 6$ shown by the magenta curve. The completeness function computed from the mock galaxies is combined with the reduced surface area to determine the survey volume, and ultimately the luminosity function.</figcaption></figure>

shown in Fig. 4. We clearly see that at very high magnifications ($\mu > 5$ mag), the survey area becomes so small that the probability to detect any galaxy is close to zero. A direct consequence of such effect is that the surface density of simulated galaxies in the image plane decreases rapidly moving closer to the critical line, as the amplification increases rapidly as well (cf. Fig. 5). For this procedure, other studies have adopted various methods. I18 adopted a similar approach to ours, where sources are simulated in the source plane using the GLAFIC software. B17 also use a source distribution in the source plane relying only on the magnification map (i.e. without using lensing deflection). An image plane method was adopted in L17, choosing to oversample the region of high magnification in the source plane. Our source plane approach with a non-uniform distribution of sources in the image plane is supposed to reflect reality since the distribution is uniform in the ‘physical’ plane. This method also ensures that galaxy shapes are distorted according to the shear potential, which directly affects the detection limit, hence the completeness function. Most importantly, the completeness function depends on several parameters, some of which are interdependent,

such as the magnification and the source position, and need to be simultaneously accounted for.

Using *LENSTOOL*, we lensed all the sources towards the image plane. Then we randomly assigned a starburst spectral energy distribution to galaxies among a library of stellar population models from Bruzual & Charlot (2003). The spectral templates span a range of ages and metallicities following normal distributions peaked at $0.2Z_{\odot}$ and 10 Myr, which are representative of the observed properties of typical high-$z$ galaxies (e.g. de Barros, Schaerer & Stark 2014). The spectral templates are redshifted and normalized to the observed magnitude, including magnification, in the F125W filter corresponding to the rest-frame UV, while the magnitudes in the rest of the filters are computed using their throughput curves. The mock galaxies are added directly to the real HFF images (10 galaxies per image at a time) of each filter and convolved with the PSF of the F160W image. Finally, we follow the same procedure used for the observations to extract the sources and select dropout galaxies at $z \sim 6$. The redshift selection function is shown in Fig. 6 together with the photometric distribution for our dropout sample.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Redshift selection function. The histogram represents the photo-metric redshift distribution of our $z \sim 6$ sample, while the blue solid curve shows the recovery efficiency in our completeness simulations marginalized over the redshift.</figcaption></figure>

# 6 EFFECTIVE SURVEY VOLUME

The completeness function is computed by comparing the output catalogue with the original input one as a function of the intrinsic magnitude. The effective survey volume is then given by

$$V_{\text{eff}} = \int_{0}^{\infty} \int_{\mu > \mu_{\text{min}}} \frac{dV_{\text{com}}}{dz} f(z, m, \mu) d\Omega(\mu, z) dz \quad (2)$$

where $V_{\text{com}}$ is the comoving volume, $f$ the completeness function, which depends on the redshift, the apparent magnitude $m$, and the magnification factor $\mu$. $d\Omega(\mu, z)$ is the surface element corresponding to a magnification $\mu$ at a given redshift $z$ and $\mu_{\text{min}}$ is the minimal magnification value at which a galaxy with a magnitude $m$ can be detected.

The resulting effective survey volume for each cluster field marginalized over the intrinsic absolute magnitude is shown in Fig. 7. The maximum volume depends mainly on the total surface area with no magnification (cf. Fig. 4) in each cluster, the maximum completeness being around 80 per cent. The curves drop rather quickly at the bright end before flattening out at the faint-end, with the help of lensing magnification. Some clusters are better lenses than others in detecting the faintest galaxies, which is a direct consequence of the balance between the magnification power and the reduction of survey area, but also the lens geometry.

# 7 THE UV LUMINOSITY FUNCTION AT $z \sim 6$

With the effective survey volume in hand, the number counts can now be used to compute the intrinsic UV luminosity function in each individual cluster as follows:

$$\phi(M_i) dM_i = \frac{N_i}{V_{\text{eff}}(M_i)}, \quad (3)$$

where $N_i$ is the number of galaxies in each magnitude bin and $V_{\text{eff}}(M_i)$ the effective survey volume in the $i$th bin of absolute magnitude $M_i$.

Before computing the final UV LF at $z \sim 6$ resulting from the combination of all the clusters, we will use our simulation procedure in the source plane to perform a realistic assessment of the different uncertainties that affect the UV LF.

<figure><img src="image_9.png" /><figcaption>**Figure 7.** The effective survey volume as a function of the intrinsic absolute magnitude in the UV. Each curve shows the result for an individual cluster. The volume is computed from the completeness function and the surface area of each cluster following equation 2 and using the CATS lens models. The coloured areas around the curves represent $1\sigma$ uncertainties.</figcaption></figure>

In order to fit the observed LF points, we adopt a modified Schechter function (cf. Bouwens et al. 2017b) to allow for a potential turnover of the LF at the faint end, which is in line with most of the theoretical models that predict a drop in the UV LF due to star formation inefficiency in small dark matter haloes (eg. Jaacks, Thompson & Nagamine 2013; Gnedin 2016; Yue, Ferrara & Xu 2016). The general Schechter form is given by

$$\phi(M) = \frac{\ln(10)}{2.5} \phi^* 10^{0.4(\alpha+1)(M^*-M)} \exp(-10^{0.4(M^*-M)}), \quad (4)$$

while for magnitudes fainter than $M_{AB} = -16$ we multiply the Schechter function by the turnover term

$$10^{-0.4\beta(M+16)^2}, \quad (5)$$

where $\beta$ is the curvature parameter. The UV LF will have a downward turnover for $\beta > 0$, but an upward turnover is also permitted for $\beta < 0$. In the following, in the case of $\beta > 0$, we define the turnover magnitude $M_T$, which corresponds to the magnitude for which $(d\phi/dM)_{M=M_T} = 0$.

## 7.1 Impact of source size distribution on the UV LF

Thanks to the sensitivity and the high spatial resolution of the ACS and WFC3 instruments onboard HST, important progress has been made during the past years in the measurements of sizes of high-redshift galaxies (e.g. Oesch et al. 2010; Huang et al. 2013;HWERDA et al. 2015). Besides the importance of these constraints on the early assembly and the structure evolution of galaxies, it directly affects the completeness calculations, hence the derived UV luminosity function. In our simulations presented in Section 5.2, we adopted a log-normal distribution with a mean half-light radius (hlr) of 150 mas and applied a size-luminosity relation as described in Section 5.2. Recent results from lensing fields suggest that the faintest galaxies in the high-redshift samples might have an even smaller size than what has been found in previous studies. In particular, Bouwens et al. (2017a) find that galaxies down to $M_{AB} \sim -15$ mag have near-point source profiles with typical half-light radii below 30 mas.
Considering the spatial distortion due to lensing shear, the ef-
fects of the size distribution might be more significant compared
to blank fields. Stacking very faint high-redshift galaxies in HFF
observations, Bouwens et al. (2016) find that they do not show
the apparent size expected for the calculated shear at their loca-
tion if one assumes the common sizes used in the literature (e.g.
the size-luminosity relation of Shibuya, Ouchi & Harikane 2015).
Rather, the observed sizes are compatible with intrinsic half-light
radii around 5 mas. To assess the impact of such size differences
on the shape of the UV LF, we adopted three different size distribu-
tions for the simulated galaxies in the source plane that are close to
what is commonly adopted in the literature. In combination with the
size-luminosity relation, galaxies fainter than M<sub>AB</sub> ~ -16 will have
hlr = 3 mas (Bouwens et al. 2017a), hlr = 20 mas (Atek et al. 2014;
Castellano et al. 2016; Ishigaki et al. 2018), and hlr = 50 mas (L17
adopted hlr = 90 mas). A comparison between the size-luminosity
relations used in different studies is given in fig. 12 of Kawamata
et al. (2018). The left-hand panel of Fig. 8 shows the effect of such
variation in galaxy sizes on the completeness function. Then such
completeness functions are applied to the observed number counts
in MACS1149 to derive the UV luminosity functions for each size
distribution (Fig. 8, right-hand panel). We use Markov Chain Monte
Carlo (MCMC) simulations to explore the full error space of each
galaxy and each effective volume in a given magnitude bin. At
the level of an individual galaxy, the MCMC simulations take into
account for the photometric scatter and magnification factor uncer-
tainties. Then, we compute the effective survey volume V<sub>eff</sub> for each
galaxy by randomly sampling the completeness functions shown in
the left-hand side of Fig. 8.

It is clear that adopting smaller sizes for very faint galaxies leads
to a larger recovery fraction, hence a larger effective survey volume,
which in turn yields a shallower faint-end slope and even a turnover
for the smallest size distributions. On the other hand, as it can be
seen in Fig. 8, larger simulated galaxies with 50 mas lead to an
upward turnover in the LF. This effect might explain the steepening
of the LF at the faint-end found by L17, who assume a normal
size distribution with a peak half-light radius of 90 mas at z ~
6. Assuming sizes below 10 mas, Bouwens et al. (2017b) find a
shallower slope as indicated in their fig. 11 because of a higher
completeness at the faint end. Similarly, Kawamata et al. (2018)
also investigated the impact of the size distribution on the faint-
end slope of the UV LF. They find a clear correlation between the
slope of the size-luminosity relation and the faint-end slope. For a
steeper size-luminosity relation, faint galaxies are smaller, hence
easier to detect. Therefore, the completeness correction becomes
smaller, which makes the faint-end slope shallower. For instance,
their derived slope for the size-luminosity relation is 0.46<sup>+0.08</sup><sub>-0.09</sub>,
leads to a faint-end slope of α = -1.86<sup>+0.17</sup><sub>-0.18</sub>. Most importantly, our
Fig. 8 shows that uncertainties on the size distribution and on the
extrapolation of the size-luminosity relation beyond M<sub>AB</sub> ~ -16
lead to large uncertainties on the faint end of UV LF, in particular
beyond M<sub>AB</sub> ~ -15, even in the case where extremely small half-
light radii (2–3 mas) are adopted. For the final LF determination
in Section 7.3, we will use the results obtained for MACS1149 to
account for the uncertainties related to the size distribution.

## 7.2 Impact of lensing uncertainties on the UV LF

The present end-to-end simulation procedure allows us to incor-
porate any lensing model in several steps of the process depicted
in Fig. 5. The choice of the lensing model will affect not only
the amplification factor but also the survey area, the completeness

function, and the multiple-images identification and positions. For
the first time, the end-to-end nature of our forward modeling han-
dles all these aspects, hence providing realistic estimates of the
systematic uncertainties by comparing the results of different lens-
ing models.

In addition to the CATS models, we selected among the public
models three different teams who provided all the necessary in-
formation to our simulations: amplification maps, deflection maps,
and shear values projected along both directions (γ₁ and γ₂). These
high-resolution maps have been ingested into LENSTOOL and interpo-
lated to treat them with the same procedures as the CATS models.
Based on the published maps from the CATS team, we have tested
that this interpolation does not affect our results compared to the
LENSTOOL parametric model.

The first comparison is the reconstruction of the source plane of
the same cluster (the example shown here is MACS1149) using each
mass model. Fig. 9 highlights large differences in the size of the
source plane, therefore the total survey area, where for example the
Diego et al. model yields a source area nearly twice as large as that of
Sharon et al. Furthermore, the differences in the shape of the source
plane also lead to different magnification and different positions for
the multiple images when lensed to the image plane. Both actual
and simulated images are impacted by the choice of lensing model.
First, the adopted model determines the amplification factor of the
observed galaxies and the position of the counter-images. Secondly,
it determines the amplification, the position, and the distortion of
the simulated galaxies (cf. Section 5), which are all linked to each
other. Overall, in combination with the source plane area, the mass
model will significantly affect the resulting completeness function.

In Fig. 10, we show the results of the different lensing models on
both the effective survey volume as a function of intrinsic absolute
magnitude and on the final UV luminosity function of MACS1149.
We can see that for the same cluster, the total survey volume can
vary by a factor of 2. Differences in intrinsic model uncertainties are
also reflected by the dispersion in each curve. The overall shape as a
function of intrinsic magnitude also varies from a model to another.
The right-hand panel of the same figure shows that the combination
of systematic and intrinsic uncertainties leads to large differences
in the final UV luminosity function beyond an intrinsic magnitude
of M<sub>UV</sub> = -15 mag. This is important since it now becomes clear
that, with the current depth and lensing models, we will not be
able to put robust constraints on the UV LF shape beyond such
magnitude limit. It is noteworthy that despite significant differences
in the survey volume between the models, the luminosity functions
agree well for magnitudes brighter than M<sub>UV</sub> = -16 mag. This
is essentially due to the fact that for relatively small amplification
factors, the variations of the amplification and the effective volume
are anticorrelated and tend to maintain the same faint-end slope.
Furthermore, relative errors on the effective volume are smaller
than what is observed at M<sub>UV</sub> > -16 mag.

In the next section, we use the same procedure in each cluster to
assess the lensing uncertainties in the final LF that combines the six
clusters.

## 7.3 Final constraints on the UV LF

We now compute the final UV luminosity function at z~6 result-
ing from the combination of all the individual clusters. In order
to account for all the uncertainties described in the previous sec-
tions, we use MCMC simulations to explore the full error space
of each galaxy and each effective volume in a given magnitude
bin. At the level of an individual galaxy, the MCMC simulations
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Impact of galaxy size distribution on the effective survey volume and the UV LF of MACS1149. In all cases the simulated sizes are based on a log-normal distribution in the source plane that is lensed to the image plane. The legend indicates the mean half light radius of galaxies fainter than $M_{AB} = -16$ mag in each simulation. The left-hand panel shows the effective survey volume for each size distribution while the right-hand panel shows its impact on the final UV LF at $z \sim 6$. The error regions for each size distribution are represented in grey in the LF plot and overlap with each other.</figcaption></figure>

<figure><img src="image_11.png" /></figure>

MACS1149

<caption>**Figure 9.** Source plane reconstruction for MACS1149 cluster using four different lensing models (Ammons et al. 2014; Johnson et al. 2014; Richard et al. 2014a; Diego et al. 2016; Jauzac et al. 2016) for which deflection maps were submitted to the HFF lensing project. The maps represent the inverse of the magnification ($1/\mu$), which is how the surface area is reduced in the source plane. The mass models affect not only the magnification maps but also the total surface area in the source plane and the spatial distribution of the sources when lensed back to the image plane.</caption>

take into account the magnification factor uncertainties by using multiple lensing models. Similarly to the previous section, we randomly sample the completeness functions and lensing models shown in the left-hand side of Fig. 10 for the corresponding cluster. For each galaxy magnitude, the MCMC procedure randomly samples the four effective volumes, which are based on four different lensing models, and their associated uncertainties. For each

iteration, the total sample is distributed in magnitude bins $M_i$, with a size of 0.5 mag, for which $\phi(M_i)$ is computed with the associated cosmic variance errors estimated from the $z=7$ lensing curve in fig. 2 of Robertson et al. (2013). These simulations have the advantage of accounting for uncertainties in $M_i$ when estimating uncertainties on $\phi(M_i)$ because individual galaxies can change bins at each iteration, and at the same time, rectifying the limitations of
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Impact of cluster mass models on the UV luminosity function in MACS1149: Left: the effective survey volume (and associated $1\sigma$ uncertainties) as a function of absolute intrinsic magnitude for each lensing model. The effective volume depends mainly on the geometrical configuration, the magnification map, and the surface area (cf. Fig 9) of each model. Right: Derived UV luminosity function using each cluster model and associated uncertainties. All models support a turnover in the faint-end of the LF but at different magnitudes and with different slopes. Beyond an absolute magnitude of $M_{UV} = -15$ mag, the differences are too large to robustly constrain the LF shape.</figcaption></figure>

the binning procedure. Simultaneously, we also fit the resulting LF results at each iteration to determine the likelihood of the modified Schechter parameters. We here combine our results with the bright-end constraints derived by Bouwens et al. (2015) from wide-area blank field surveys. To account for the difference in the redshift selection function between their $z = 7$ sample and our combined $z = 6-7$ one, we rescale their UV LF by 0.15 dex.

The final LF results are shown in Fig. 11, together with the LF values listed in Table 3. We show on the left-hand panel the best-fitting constraints using a classical Schechter function with the associated confidence regions. It is noteworthy that, unlike LF measurements that extend to fainter magnitudes with a nearly constant slope (I17, L17), there is clear drop in our LF at faint magnitudes. Therefore, it is sensible to truncate the luminosity range to bins brighter than $M_{UV} = -15$ when using a Schechter function with a constant faint-end slope. We find a faint-end slope of $\alpha = -1.98_{-0.09}^{+0.11}$, which is steeper than the B17 estimate of $\alpha = -1.91 \pm 0.02$, but shallower than estimates of L17 and I18 with $\alpha = -2.10_{-0.03}^{+0.03}$ and $\alpha = -2.15_{-0.06}^{+0.08}$, respectively. Importantly, the uncertainties of the best-fitting parameters derived in this work are significantly larger than those studies.

A better comparison is achieved with the best-fitting function that allows for a curvature at very faint magnitudes following equation (5). Now we obtain a best-fitting faint-end slope of $\alpha = -2.01_{-0.14}^{+0.12}$ and a curvature parameter of $\beta = 0.48_{-0.25}^{+0.49}$. Our results show a stronger turnover than what is found by B17, who extend the LF to $\sim -13$ mag with a slight turnover of $\beta = 0.17 \pm 0.2$. From equation 5, we determine the turnover magnitude of the LF at $M_T = -14.93_{-0.52}^{+0.61}$, which is close to the result of B17 at $M_T = -14.9$. Using empirical results from four HFF clusters and a luminosity function model, Yue et al. (2017) define $M_T$ as the magnitude at which the classical Schechter LF drops by a factor of 2 and report $M_T \ge -14.3$. Comparing to our definition and their fig. 5, they find a similar turnover magnitude but with a shallower curvature. The results of L17 show a strong evidence against a turnover at magnitudes brighter than $M_{UV} = -12.5$ at $z = 6$.

The discrepancies observed between our results and other studies could partially be explained by the larger number of galaxies reported in the sample of L17 due to a less drastic selection

criteria and a more sophisticated correction for ICL contamination, or by a larger size distribution in their completeness simulations (cf. Section 7.1). B17 find a much shallower curvature in the LF at a similar magnitude compared to our study. They use very small sizes (around 3 mas) for faint galaxies in their completeness simulations based on the latest results of high-z studies including the HFF clusters. Such distribution naturally leads to a shallower slope at the faint end. Most importantly, the uncertainties derived from our end-to-end simulations prevent us from putting robust constraints on the very faint-end shape of the LF. Most of the lensing and sizes uncertainties investigated by this procedure impact the faint end of the LF for which the $2\sigma$ confidence region allows a wide range of curvature and slope parameters.

## 8 THE $z \sim 6$ UV LUMINOSITY DENSITY

Using our UV LF determination we now compute the total galaxy UV luminosity density $\rho_{UV}$ at $z = 6$. The MCMC simulations performed for the constraints on the LF parameters are propagated to compute the confidence intervals of the luminosity density. We investigate different truncation magnitudes (the faint integration limit of the LF) to determine to what extent galaxies can produce enough UV radiation to sustain reionization. The UV luminosity density as a function of the limiting magnitude is shown in Fig. 13. Beyond a magnitude of $M_{UV} = -15$, which corresponds to the turnover of the luminosity function, the luminosity density becomes nearly flat (on a log scale), reaching $\text{Log}(\rho_{UV}/\text{erg s}^{-1} \text{Mpc}^{-3}) = 26.15 \pm 0.09$ at $M_{\text{lim}} = -15$ and $\text{Log}(\rho_{UV}/\text{erg s}^{-1} \text{Mpc}^{-3}) = 26.21 \pm 0.13$ at $M_{\text{lim}} = -10$. In comparison, B17 report a luminosity density of $\text{Log}(\rho_{UV}/\text{erg s}^{-1} \text{Mpc}^{-3}) = 26.33$ at $M_{\text{lim}} = -15$, slightly larger than our determination due to their larger value of $\phi^* = 0.58 \times 10^{-3} \text{Mpc}^{-3}$. In the case where a Schechter function with a constant faint-end slope is used, I18 and L17 find smaller values resulting from different constraints used in the bright end of the LF, which yield smaller values for $\phi^*$. When we use a Schechter fit, we find $\text{Log}(\rho_{UV}/\text{erg s}^{-1} \text{Mpc}^{-3}) = 26.41 \pm 0.14$ at $M_{\text{lim}} = -10$. Again, the confidence region can accommodate a wide range of values for the total luminosity density and the shape of its evolution with the truncation magnitude and remains compatible with most of the literature results.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** The final UV luminosity function at $z \sim 6$. The left-hand panel shows the best fit of the UV luminosity function using a Schechter form (gold curve) with the associated $1\sigma$ (2$\sigma$) uncertainties in dark (light) blue areas. Black points represent the bright-end results of Bouwens et al. (2015). Measurements (gold data points with $1\sigma$ errors) at the very faint-end clearly depart from a simple Schechter form. Therefore, the fit is restricted to magnitudes brighter than $M_{UV} = -15$ (marked by the vertical line). For comparison, $z \sim 6$ UV LF results from B17, L17, and I18 are also shown with the red, green, and brown curves, respectively. For consistency, following B17, we shift the B17 and L17 LFs down by 0.15 dex, corresponding to the expected redshift evolution of the LF between their $z = 6$ sample and our combined $z = 6-7$ sample. Right: same as left-hand panel but measurements are now fitted with a modified Schechter function (cf. equation 5) that allows a curvature fainter than $M_{UV} = -16$. Note that among the literature results only Bouwens et al. (2017b) allow a turnover at faint magnitudes.</figcaption></figure>

<table><caption>**Table 3.** Best-fitting constraints on the $z \sim 6$ UV luminosity function.</caption><thead><tr><th>Reference</th><th><i>M</i><sup>*</sup><sub>UV</sub><br>(AB mag)</th><th>α</th><th>log<sub>10</sub>φ<sup>*</sup><br>(Mpc<sup>-3</sup>)</th><th>β<sup>c</sup></th><th><i>M</i><sub>T</sub><sup>d</sup><br>(AB mag)</th></tr></thead><tbody><tr><td>This work<sup>a</sup></td><td>-20.74<sup>+0.21</sup><sub>-0.20</sub></td><td>-1.98<sup>+0.11</sup><sub>-0.09</sub></td><td>-3.43<sup>+0.21</sup><sub>-0.21</sub></td><td>-</td><td>-</td></tr><tr><td>Atek et al. (2015b)</td><td>-20.90<sup>+0.90</sup><sub>-0.73</sub></td><td>-2.01<sup>+0.20</sup><sub>-0.28</sub></td><td>-3.55<sup>+0.37</sup><sub>-0.57</sub></td><td>-</td><td>-</td></tr><tr><td>Bouwens et al. (2017b)<sup>e</sup></td><td>-20.94</td><td>-1.91 ± 0.02</td><td>-3.18 ± 0.03</td><td>-</td><td>-</td></tr><tr><td>Ishigaki et al. (2018)</td><td>-20.89<sup>+0.17</sup><sub>-0.13</sub></td><td>-2.15<sup>+0.08</sup><sub>-0.06</sub></td><td>-3.78<sup>+0.15</sup><sub>-0.15</sub></td><td>-</td><td>-</td></tr><tr><td>Livermore et al. (2017)</td><td>-20.82<sup>+0.04</sup><sub>-0.03</sub></td><td>-2.10<sup>+0.03</sup><sub>-0.03</sub></td><td>-3.64<sup>+0.04</sup><sub>-0.03</sub></td><td>-</td><td>-</td></tr><tr><td>This work<sup>b</sup></td><td>-20.84<sup>+0.27</sup><sub>-0.30</sub></td><td>-2.01<sup>+0.12</sup><sub>-0.14</sub></td><td>-3.54<sup>+0.06</sup><sub>-0.07</sub></td><td>0.48<sup>+0.49</sup><sub>-0.25</sub></td><td>-14.93<sup>+0.61</sup><sub>-0.52</sub></td></tr><tr><td>Bouwens et al. (2017b)<sup>e</sup></td><td>-20.94</td><td>-1.91 ± 0.04</td><td>-3.24 ± 0.04</td><td>0.17 ± 0.2</td><td>-14.9</td></tr></tbody></table>

**Notes:**

<sup>a</sup>Using a Schechter function fit

<sup>b</sup>Allowing for a turnover in the LF fainter than $M_{UV} = -16$ mag

<sup>c</sup>The turnover parameter presented in equation 5 that quantifies the curvature of the LF at magnitudes fainter than $M_{UV} = -16$ mag

<sup>d</sup>The turnover magnitude at which the LF departs from the simple Schechter function

<sup>e</sup>Using the CATS model

Therefore, it is clear that accounting for the lens model uncertainties no robust constraints can be inferred regarding the existence of a large population of faint ($M_{UV} > -15$) galaxies that could provide the required energy budget to reionize the Universe.

## 9 CONCLUSION

The lensing clusters of the HFF recently made it possible to explore the distant Universe to the faintest magnitude limits ever achieved. Based on the unprecedented quality of the HST images, several groups have constructed mass models in order to interpret high-redshift observations. For instance, the galaxy luminosity function and the ability of faint galaxies to reionize the Universe was the focus of multiple studies from different groups (e.g. Atek et al. 2015b; Bouwens et al. 2017b; Livermore et al. 2017; Ishigaki et al. 2018; Yue et al. 2017). With the completion of the HFF program,

significant differences were found between publications at the very faint end of the LF, where the lensing magnification, and their associated complex uncertainties become important.

In this paper, we used the six HFF clusters data to construct a robust UV LF at $z \sim 6$. We have used for the first time a comprehensive approach based on end-to-end simulations that includes all lensing uncertainties and their impact on the study of the UV luminosity function at $z \sim 6$. Our approach uses simulated galaxies directly in the source plane, which has the advantage of accounting for all lensing uncertainties present in both effective survey volume and the number counts.

We first used this procedure to assess the impact of the source size distribution on the survey volume, hence on the UV LF at $z \sim 6$. We find that uncertainties on the size–luminosity relation has a significant impact on the faint end of the UV LF. Adopting large sizes with an hlr = 50 mas leads to a small recovery fraction
<figure><img src="image_14.png" /><figcaption>**Figure 12.** The final UV luminosity function at $z \sim 6$ compared to recent determinations in the literature. Black points represent the bright-end results of Bouwens et al. (2015). The two best-fitting curves are shown in gold while LF measurements of Livermore et al. (2017) and Bouwens et al. (2017b) are shown with green squares and red diamonds, respectively. These two literature curves were shifted down by 0.15 dex to account for the difference in the mean redshift of the LF described in the text. Note that the data points of L17 were not corrected for Eddington bias, while their best-fitting LF determination was corrected.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 13.** The UV luminosity density at $z \sim 6$ as a function of the lower integration limit of the UV LF is illustrated by the gold curve. The shaded light (dark) blue regions are the $1\sigma$ ($1\sigma$) uncertainties. Other $\rho_{UV}$ determinations of B17, L17, and I18 are shown with the red, green, and brown curves, respectively. Note that the LF of B17 and L17 used to compute $\rho_{UV}$ was shifted down by 0.15 dex to account for a redshift evolution (cf. Section 7.3).</figcaption></figure>

therefore to a steeper faint-end slope (e.g. Livermore et al. 2017). Smaller sizes with half light radii below 10 mas lead to a shallower slope (Bouwens et al. 2017b). Overall, such uncertainties create a wide range of slopes beyond $M_{UV} = -15$ mag.

Our procedure has the second advantage of using any mass model and therefore provides means to assess the systematic lensing uncertainties. We achieve this goal by comparing the results, using four different models to compute the effective survey volume and

the UV LF. We show that the combination of systematic and intrinsic uncertainties leads to important differences in the final UV luminosity function beyond an intrinsic magnitude of $M_{UV} = -15$, where most galaxies will have a magnification factor greater than 10.

Finally, we computed the UV luminosity function while incorporating the different uncertainties discussed in this paper. Adopting a simple Schechter fit, we find a faint-end slope of $\alpha = -1.98_{-0.09}^{+0.11}$, whereas a modified Schechter function that permits curvature in the LF at $M_{UV} > -16$ mag yields a turnover in the LF with a faint-end slope of $\alpha = -2.01_{-0.14}^{+0.12}$ and a curvature parameter of $\beta = 0.48_{-0.25}^{+0.49}$. Most importantly, while galaxies were detected down to an intrinsic magnitude of $M_{UV} \sim -13$, we were unable to reliably extend the UV LF beyond $M_{UV} \sim -15$ because of the large confidence interval. Consequently, the existence of a large reservoir of faint galaxies that significantly contribute to the total UV luminosity density is still uncertain.

The existence of such a turnover at faint magnitudes is also suggested by near-field studies in order to match the slope of the local galaxy luminosity function. A possible explanation provided by such studies and numerical simulations is that star formation becomes inefficient in dark matter haloes smaller than $10^9 M_{\odot}$ at those early epochs, which would correspond to a rollover of the UV LF around $M_{UV} \sim -14$ (Boylan-Kolchin, Bullock & Garrison-Kimmel 2014; Boylan-Kolchin et al. 2015). The local group studies also suggest that the break in the UV LF should be shallower to allow the existence of ultra-faint dwarfs down to $M_{UV} \sim -3$ (Weisz & Boylan-Kolchin 2017).

Overall, we demonstrated that with the current depth of observations and current state of the art in mass modeling of lensing clusters, robust constraints on the UV luminosity function fainter than $M_{UV} = -15$ mag remain unrealistic. Future observations of lensing clusters with the upcoming *James Webb Space Telescope* will push observed flux limits by about 2 mag and at the same time provide hundreds of spectroscopic redshifts of multiple images to improve the accuracy of lensing models. Such observations will therefore bring a definitive answer to the potential turnover in the UV LF and the contribution of extremely faint galaxies to cosmic reionization.

## ACKNOWLEDGEMENTS

We thank the STScI and the HFF team for their efforts in obtaining and reducing the *HST* data. HA is supported by the Centre National d'Etudes Spatiales (CNES). JPK is supported by the European Research Council (ERC) advanced grant 'Light in the Dark' (LIDA). JR acknowledges support from the ERC starting grant CALENDAS.

This study is based on observations made with the NASA/ESA *Hubble Space Telescope*, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with programs 13495, 11386, 13389, and 11689. STScI is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555. This work utilizes gravitational lensing models produced by PIs Bradac, Natarajan & Kneib (CATS), Merten & Zitrin, Sharon, and Williams, and the GLAFIC and Diego groups. This lens modeling was partially funded by the *HST* Frontier Fields program conducted by STScI. The Hubble Frontier Fields data and the lens models were obtained from the Mikulski Archive for Space Telescopes (MAST).
REFERENCES

Ammons S. M., Wong K. C., Zabludoff A. I., Keeton C. R., 2014, ApJ,
781, 2

Atek H. et al., 2014, ApJ, 786, 60

Atek H. et al., 2015a, ApJ, 800, 18

Atek H. et al., 2015b, ApJ, 814, 69

Bertin E., Arnouts S., 1996, A&AS, 117, 393

Bolzonella M., Miralles J.-M., Pelló R., 2000, A&A, 363, 476

Bouwens R. J. et al., 2011, ApJ, 737, 90

Bouwens R. J. et al., 2014, ApJ, 803, 49

Bouwens R. J. et al., 2015, ApJ, 803, 34

Bouwens R. J., Illingworth G. D., Blakeslee J. P., Broadhurst T. J., Franx
M., 2004, ApJ, 611, L1

Bouwens R. J., Illingworth G. D., Oesch P. A., Atek H., Lam D., Stefanon M., 2016, ApJ, 843, 41

Bouwens R. J., van Dokkum P. G., Illingworth G. D., Oesch P. A., Maseda M., Ribeiro B., Stefanon M., Lam D., 2017a, ApJL, preprint (arXiv: 1711.02090)

Bouwens R. J., Oesch P. A., Illingworth G. D., Ellis R. S., Stefanon M., 2017b, ApJ, 843, 129

Boylan-Kolchin M., Bullock J. S., Garrison-Kimmel S., 2014, MNRAS, 443, L44

Boylan-Kolchin M., Weisz D. R., Johnson B. D., Bullock J. S., Conroy C., Fitts A., 2015, MNRAS, 453, 1503

Bradač M., Schneider P., Lombardi M., Erben T., 2005, A&A, 437, 39

Bruzual G., Charlot S., 2003, MNRAS, 344, 1000

Bunker A. J. et al., 2010, MNRAS, 409, 855

Castellano M. et al., 2016, ApJ, 823, L40

Chabrier G., Baraffe I., Allard F., Hauschildt P., 2000, ApJ, 542, 464

Coe D., Bradley L., Zitrin A., 2015, ApJ, 800, 84

Coleman G. D., Wu C.-C., Weedman D. W., 1980, ApJS, 43, 393

de Barros S., Schaerer D., Stark D. P., 2014, A&A, 563, A81

Diego J. M. et al., 2016, MNRAS, 456, 356

Diego J. M., Broadhurst T., Zitrin A., Lam D., Lim J., Ford H. C., Zheng
W., 2015, MNRAS, 451, 3920

Fan X. et al., 2006, AJ, 132, 117

Ferguson H. C. et al., 2004, ApJ, 600, L107

Finkelstein S. L. et al., 2015, ApJ, 810, 71

Gnedin N. Y., 2016, ApJ, 825, L17

Grazian A. et al., 2012, A&A, 547, A51

Grillo C. et al., 2015, ApJ, 800, 38

Hathi N. P., Jansen R. A., Windhorst R. A., Cohen S. H., Keel W. C., Corbin M. R., Ryan R. E., Jr., 2008, AJ, 135, 156

Hoag A. et al., 2016, ApJ, 831, 182

Holwerda B. W., Bouwens R., Oesch P., Smit R., Illingworth G., Labbe I., 2015, ApJ, 808, 6

Huang K.-H., Ferguson H. C., Ravindranath S., Su J., 2013, ApJ, 765, 68

Ishigaki M., Kawamata R., Ouchi M., Oguri M., Shimasaku K., Ono Y., 2015, ApJ, 799, 12

Ishigaki M., Kawamata R., Ouchi M., Oguri M., Shimasaku K., Ono Y., 2018, ApJ, 854, 73

Jaacks J., Thompson R., Nagamine K., 2013, ApJ, 766, 94

Jauzac M. et al., 2014, MNRAS, 443, 1549

Jauzac M. et al., 2015, MNRAS, 452, 1437

Jauzac M. et al., 2016, MNRAS, 457, 2029

Johnson T. L., Sharon K., Bayliss M. B., Gladders M. D., Coe D., Ebeling
H., 2014, ApJ, 797, 48

Jullo E., Kneib J.-P., 2009, MNRAS, 395, 1319

Jullo E., Kneib J.-P., Limousin M., Elíasdóttir Á., Marshall P. J., Verdugo
T., 2007, New J. Phys., 9, 447

Kawamata R., Ishigaki M., Shimasaku K., Oguri M., Ouchi M., 2015, ApJ,
804, 103

Kawamata R., Oguri M., Ishigaki M., Shimasaku K., Ouchi M., 2016, ApJ,
819, 114

Kawamata R., Ishigaki M., Shimasaku K., Oguri M., Ouchi M., Tanigawa
S., 2018, ApJ, 855, 4

Kinney A. L., Calzetti D., Bohlin R. C., McQuade K., Storchi-Bergmann
T., Schmitt H. R., 1996, ApJ, 467, 38

Kneib J.-P., 1993, PhD thesis, Université Paul Sabatier, Toulouse

Kneib J.-P., Natarajan P., 2011, A&A Rev., 19, 47

Krist J. E., Hook R. N., Stoehr F., 2011, in Optical Modeling and Perfor-
mance Predictions V, Proceedings of SPIE - The International Society
for Optical Engineering, p. 81270J

Lagattuta D. J. et al., 2017, MNRAS, 469, 3946

Laporte N. et al., 2016, ApJ, 820, 98

Limousin M. et al., 2016, A&A, 588, A99

Livermore R. C., Finkelstein S. L., Lotz J. M., 2017, ApJ, 835, 113

Lotz J. M. et al., 2017, ApJ, 837, 97

Mahler G. et al., 2018, MNRAS, 473, 663

Maizy A., Richard J., de Leo M. A., Pelló R., Kneib J. P., 2010, A&A, 509, A105

McLure R. J. et al., 2013, MNRAS, 432, 2696

Merten J. et al., 2011, MNRAS, 417, 333

Mosleh M. et al., 2012, ApJ, 756, L12

Oesch P. A. et al., 2010, ApJ, 709, L16

Oesch P. A. et al., 2014, ApJ, 786, 108

Oke J. B., Gunn J. E., 1983, ApJ, 266, 713

Ono Y. et al., 2013, ApJ, 777, 155

Peng C. Y., Ho L. C., Impey C. D., Rix H.-W., 2002, AJ, 124, 266

Planck Collaboration XLVII, 2016, A&A, 596, A108

Postman M. et al., 2012, ApJS, 199, 25

Richard J. et al., 2014a, MNRAS, 444, 268

Richard J. et al., 2014b, MNRAS, 444, 268

Robertson B. E. et al., 2013, ApJ, 768, 71

Schaerer D., de Barros S., 2009, A&A, 502, 423

Sharon K., Gladders M. D., Rigby J. R., Wuyts E., Koester B. P., Bayliss M. B., Barrientos L. F., 2012, ApJ, 746, 161

Shibuya T., Ouchi M., Harikane Y., 2015, ApJS, 219, 15

Shipley H. et al., 2018, ApJS, 235, 14

Steidel C. C., Giavalisco M., Pettini M., Dickinson M., Adelberger K. L., 1996, ApJ, 462, L17

Weisz D. R., Boylan-Kolchin M., 2017, MNRAS, 469, L83

Yue B. et al., 2017, ApJ, preprint (arXiv:1711.05130)

Yue B., Ferrara A., Vanzella E., Salvaterra R., 2014, MNRAS, 443, L20

Yue B., Ferrara A., Xu Y., 2016, MNRAS, 463, 1968

Zheng W. et al., 2014, ApJ, 795, 11

Zitrin A. et al., 2013, ApJ, 762, L30

Zitrin A., Broadhurst T., Rephaeli Y., Sadeh S., 2009, ApJ, 707, L102

This paper has been typeset from a TeX/LA TeX file prepared by the author.