COSMOS2020: UV-selected galaxies at z ≥ 7 . 5

This paper presents a new search for z ≥ 7 . 5 galaxies using the COSMOS2020 photometric catalogues. Finding galaxies at the reionisation epoch through deep imaging surveys remains observationally challenging. The larger area covered by ground-based surveys such as COSMOS enables the discovery of the brightest galaxies at these high redshifts. Covering 1 . 4deg 2 , our COSMOS catalogues were constructed from the latest UltraVISTA data release (DR4) combined with the ﬁnal Spitzer / IRAC COSMOS images and the Hyper-Suprime-Cam Subaru Strategic Program DR2 release. We identiﬁed 17 new 7 . 5 < z < 10 candidate sources, and conﬁrm 15 previously published candidates. Using deblended photometry extracted by ﬁtting surface brightness models on multi-band images, we selected four candidates which would be rejected using ﬁxed aperture photometry. We tested the robustness of all our candidates by comparing six di ﬀ erent photometric redshift estimates. Finally, we computed the galaxy UV luminosity function in three redshift bins centred at z = 8 , 9 , 10. We ﬁnd no clear evolution of the number density of the brightest galaxies M UV < − 21 . 5, in agreement with previous works. Rapid changes in the quenching e ﬃ ciency or attenuation by dust could explain such a lack of evolution between z ∼ 8 and z ∼ 9. A spectroscopic conﬁrmation of the redshifts, already planned with JWST and the Keck telescopes, will be essential to conﬁrm our results.


Introduction
Star-forming galaxies are expected to be major contributors to the reionisation process in the primordial Universe (e.g. Dayal & Ferrara 2018;Finkelstein et al. 2019). During this time, when the Universe was less than 1 Gyr old, Lyman continuum photons emitted by newborn massive stars ionised the intergalactic medium around them. The efficiency of this process is expected to depend on the density of galaxies, their star formation rate (SFR), and the fraction of ionising photons escaping from their interstellar medium (e.g. Bouwens et al. 2015a). The latest results from the Planck survey (Planck Collaboration et al. 2020) give the reionisation redshift mid-point z re = 7.67 ± 0.73, while the reionisation probably ended around redshift z = 6 (e.g. Fan et al. 2006). However, an accurate timeline showing the interplay between the formation of the first galaxies and reionisation has not been established yet. At high redshifts, the galaxy UV luminosity function (UVLF) is the main observable to establish the cosmic star formation rate density (SFRD). The num-ber of ionising photons generated by galaxies can be derived from the SFRD, assuming a number of Lyman continuum photons per second produced per unit SFR, and a fraction of them escaping from galaxies (e.g. Robertson et al. 2015;Vanzella et al. 2018). Furthermore, because of the fluctuation of the intergalactic medium (IGM) opacity on large scales (e.g. Kulkarni et al. 2019), the observed redshift of the end of the reionisation may depend on sight lines. Therefore, a census of star-forming galaxies at z > 6 in deep fields is essential to understand the origin and physics of reionisation, as well as the formation of the first galaxies.
Neutral hydrogen in the IGM, the circumgalactic medium, and the interstellar medium absorbs the flux blueward of the Lyman limit at 912 Å, with a flux expected to be consistent with zero below this limit (e.g. Bouwens et al. 2015b). Moreover, the IGM absorbs the light bluer than the Lyman alpha emission line at 1215.67 Å. These features have been used to select highredshift galaxies in deep imaging surveys (Steidel et al. 1996).
The search for these Lyman-break galaxies (LBGs) is carried out by using fluxes straddling the Lyman emission wavelength at any given redshift and detecting the amplitude of this feature (the height of the Lyman break). The colour redward of the break constrains the rest-frame UV slope, which is sensitive to star formation, dust, and metallicity. Photometric redshifts are also commonly used to select high-redshift candidates, through spectral energy distribution (SED) fitting (e.g. McLure et al. 2013;Finkelstein et al. 2015;Bowler et al. 2015). This method presents the advantage of using more than three bands if available, as well as providing a redshift probability distribution function (hereafter zPDF) for each candidate.
Recent deep and wide-area surveys have provided the necessary multi-wavelength data to select high-redshift galaxies, in particular at near-infrared wavelengths (1 − 2 µm). The Wide Field Camera 3 (WFC3) onboard the Hubble space telescope (HST), with its unmatched sensitivity at near-infrared wavebands, has revolutionised our knowledge of the distant Universe. One of the major surveys in this context is the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS; Grogin et al. 2011;Koekemoer et al. 2011), covering a total of 750 arcmin 2 with deep HST imaging. The CANDELS fields include the Great Observatories Origins Deep Survey (GOODS; Giavalisco et al. 2004) with two fields (GOODS-South and GOODS-North), COSMOS , the Extended Groth Strip (EGS; Davis et al. 2007), and the UKIDSS Ultradeep Survey field (UDS; Cirasuolo et al. 2007). Furthermore, the Hubble Ultra-Deep Field (HUDF) covering 4.7 arcmin 2  has allowed for the deepest ever astronomical images to be taken, reaching depths of 30 mag in the optical and near-infrared. Alternatively, the Hubble Frontier Fields (HFFs; Lotz et al. 2017) include deep imaging of six massive galaxy clusters, which are used to study extremely faint background galaxies magnified through gravitational lensing. Using these multi-wavelength surveys, over a thousand galaxies with z > 6 have been identified, reaching z ∼ 10 (e.g. McLeod et al. 2016;Oesch et al. 2014Oesch et al. , 2018Bouwens et al. 2015bBouwens et al. , 2019. Based on galaxy candidates detected in these HST deep fields, there is evidence that the Schechter function (Schechter 1976) still provides a reasonable fit to the UVLF at high redshift (e.g. Schmidt et al. 2014;Bouwens et al. 2015b). The shape of the UVLF appears to evolve with redshift, with a steepening of the faint-end slope α from −1.6 at z = 4 to −2.3 at z = 10 (e.g. Bouwens et al. 2015b;Finkelstein et al. 2015;Oesch et al. 2018). The observed steep slope for the UVLF at high redshifts predicts a high abundance of faint galaxies in the Universe during the reionisation epoch. However, there is no consensus whether the observed evolution of the UVLF at z > 6 is driven by changes in the characteristic absolute magnitude or in the normalisation (e.g. McLure et al. 2013;Finkelstein et al. 2015). This may be due to an evolving shape of the LF towards a non-Schechter form. When converted in SFRD, most of the studies are consistent with a steady decrease in the SFRD from its peak at z = 2−3 to z = 10 (e.g. Madau & Dickinson 2014;Finkelstein 2016). Oesch et al. (2018) find that this decrease is even more rapid between z = 8 and z = 10 by more than a factor ten, possibly explained by the fast buildup of the underlying dark matter mass function in the primordial Universe.
While an extremely powerful technique to study the faint galaxy populations, the current deep pencil-beam surveys do not cover sufficient cosmological volumes to capture the rare and bright sources which would constrain the bright end of the UVLF. For instance, Oesch et al. (2018) analysed 800 arcmin 2 of archival data, but did not constrain the LF brighter than M UV < −21.2. Bridge et al. (2019) found eight galaxies brighter than M UV < −21.5 at 7 < z < 8 in the Brightest of Reionising Galaxies (BoRG) survey, a parallel HST survey specially designed to observe the brightest galaxies. In this regime, deep ground-based imaging surveys are invaluable. The combination of deep optical (27 − 28 mag) and near-infrared (25 − 26 mag) imaging over a degree-scale area has made it possible to isolate the brightest galaxies in the z > 7 − 8 Universe (e.g. Bowler et al. 2014). Sufficiently deep Spitzer coverage with the InfraRed Array Camera (IRAC) at 3.6 and 4.5 µm complements these observations, allowing for the efficient removal of the intermediateredshift contaminants (e.g. Oesch et al. 2012), as well as the direct detection of high-redshift sources.
Using the available data from the UltraVISTA survey (Mc-Cracken et al. 2012) in the COSMOS field, multiple samples of galaxy candidates have already been identified at z > 6 (Bowler et al. 2014(Bowler et al. , 2015 and z ≥ 7.5 (Stefanon et al. 2017 hereafter S19), imposing strong constraints on the bright end of the galaxy UVLF. Bowler et al. (2020, hereafter B20) found 27 LBGs over a 6 deg 2 area combining data from the COSMOS and XMM-LSS fields. The bright end (−23 < M UV < −21) of the resulting galaxy number density is in excess compared to the exponential decline predicted from a Schechter parametrisation, suggesting a double power law may be more appropriate. Early models have linked the UVLF evolution with the one of the dark matter halo mass functions (Cooray & Ouchi 2006;Bouwens et al. 2008;Tacchella et al. 2013). The change of shape at the bright end of the UVLF compared to the halo mass function is interpreted by the combination of several physical processes, including quenching and dust extinction (e.g. Harikane et al. 2022). Indeed, for a given massive galaxy, its host dark matter halo or its own intrinsic stellar content may not be massive enough to trigger internal star-formation quenching (Peng et al. 2010). Moreover, dust extinction is known to bend the bright end of the UVLF at intermediate redshift (Reddy et al. 2010). The formation of dust in the interstellar medium (ISM) of galaxies at z > 6 may still be immature, thus high SFR galaxies at this epoch may experience less attenuation (Finkelstein et al. 2015). This would result in a UVLF that is well populated at the bright end, describing a double power law (Bowler et al. 2014(Bowler et al. , 2015. With the COSMOS2020 photometric catalogues (Weaver et al. 2022), we have the opportunity to improve the search and identify the brightest star-forming galaxies at the epoch of reionisation. Candidates are identified using deep near-infrared imaging from UltraVISTA DR4, complemented with IRAC images from the Cosmic Dawn Survey (Euclid Collaboration et al. 2022) to detect the galaxy rest-frame optical emission. While B20 used the latest UltraVISTA release (DR4), we used, in addition, the latest deep optical images from the public DR2 of the Hyper Suprime-Cam (HSC) Subaru Strategic Program (HSC-SSP; Aihara et al. 2018) and the final Spitzer images in the mid-infrared. We also validated our results using the recently released DR3 of the HSC-SSP (Aihara et al. 2022). These data are essential to improve the purity of the high-redshift sample and to extend the area under investigation. Furthermore, COSMOS2020 photometric extractions have been made with two different techniques: a traditional approach using aperture photometry and a surface brightness profile-fitting technique using multi-band images. In this paper, we combine both photometric catalogues, and multiple photometric redshift codes, leading to a more robust final sample of candidates.
For this study, we searched for galaxies at z ≥ 7.5 in the COSMOS field using the COSMOS2020 photometric catalogue, and we estimated the galaxy rest-frame UVLF at 8 ≤ z ≤ 10. The paper is structured as follows. Section 2 presents the imaging data used in this work, as well as the estimated source photometry. Section 3 describes the high-redshift galaxy selection. The sample of candidates is presented in Sect. 4 and compared with previous studies. The updated constraints on the highredshift UVLF resulting from the selected galaxies are presented in Sect. 5. We discuss these results in Sect. 6 and summarise our conclusions in Sect. 7. We adopted the standard Λ Cold Dark Matter cosmology with Ω m = 0.3, Ω Λ = 0.7, and H 0 = 70 km s −1 Mpc −1 . We used the initial mass function (IMF) from Chabrier (2003). Magnitudes are given in the AB system (Oke 1974).

Data
This work is based on the COSMOS2020 photometric catalogues (Weaver et al. 2022), which includes the currently deepest optical to mid-infrared imaging in the COSMOS field. The imaging data, as well as the extracted photometry are briefly described here. Full details can be found in Weaver et al. (2022).

Imaging
The COSMOS field is covered by several multi-wavelength deep imaging surveys. While the photometric catalogue consists of observations in more than 35 bands, here we only discuss the deepest and reddest broad-band imaging, which are the most relevant datasets to search for high-redshift galaxies. The photometric depths of these broad-band images are given in Table 1. We stress that the depths are not homogeneous over the full field and differ by wavelength. We discuss the impact on the selection of the high-redshift candidates in Sect. 4.
The UltraVISTA survey (McCracken et al. 2012) provides deep near-infrared imaging over 1.5 deg 2 of the COSMOS field, as shown in Fig. 1. The exposure time is not homogeneous over the entire field. Four ultra-deep stripes across the field covering 0.62 deg 2 have deeper exposure times (vertical dark grey stripes in Fig. 1). We use the UltraVISTA DR4 data, with a near-infrared coverage in four broad bands, Y JHK s . The additional narrowband NB118 covers the ultra-deep stripes.
The mid-infrared data comes from the Cosmic Dawn Survey (Euclid Collaboration et al. 2022), in which all the Spitzer observations in the IRAC/[3.6], [4.5] bands are processed (the full list of included programs is given in their  Ashby et al. 2018) and Completing the Legacy of Spitzer/IRAC over COSMOS (COM-PLETE, COMPLETE2;Labbe et al. 2016;Stefanon et al. 2018). The IRAC depths reported in Table 1 are separately computed inside the deep and ultra-deep UltraVISTA stripes. The 3σ depths of S-COSMOS data outside the SPLASH field (shown in Fig. 1) reach 25.3 mag in [3.6] and 25 mag in [4.5].
The optical data provided by the Hyper Suprime-Cam (HSC) Subaru Strategic Program (HSC-SSP; Aihara et al. 2018), include deep imaging in the g, r, i, z, y broad bands over 2.2 deg 2 . We use the HSC-SSP DR2 images (Aihara et al. 2019), stacked for the COSMOS2020 catalogue. The central region of the COS-MOS field includes ultra-deep HSC imaging within a 0.75 deg radius circle (blue circle in Fig. 1), and deep imaging in the extended COSMOS survey, about 0.5 mag shallower than in the centre. Since the COSMOS2020 catalogue was created, a third release of the HSC-SSP became available (Aihara et al. 2022). While these data are not included in our analysis, we used these images to insure the robustness of our candidates. In addition, we include the reprocessed Subaru Suprime-Cam images with 12 medium and two narrow bands in optical (Taniguchi et al. 2007(Taniguchi et al. , 2015, the u and u * bands from the CFHT Large Area Uband Deep Survey (CLAUDS; Sawicki et al. 2019), and the UV photometry from GALEX (Zamojski et al. 2007).

Photometry
Source detection is performed in the combined CHI_MEAN image constructed with SWarp (Bertin et al. 2002) from the HSC/i, z and UltraVISTA/Y, J, H, K s bands. This stacked detection image provides an advantage for faint, high-redshift sources compared to single images. Photometry is extracted following two independent approaches, leading to two separate photometric catalogues, called the Classic catalogue and The Farmer catalogue. In the Classic approach, the photometry of the high-resolution images is performed with SExtractor (Bertin & Arnouts 1996) in dual-image mode. The photometry is extracted in fixed 2 diameter apertures. To ensure that the apertures include the same features at all wavelengths, the point-spread functions (PSF) of the science images are homogenised. Multiple corrections are applied to the measured magnitudes, including the magnitude error scaling and the aperture-to-total magnitude corrections. The photometry of the low-resolution IRAC images is performed with the software IRACLEAN (Hsieh et al. 2012), using the highresolution images as prior to extract the photometry of the confused sources. In this case, the IRAC PSF is iteratively subtracted from the IRAC images centred within the source boundary, as defined by the high-resolution prior. In the second approach, all bands are extracted with The Farmer (Weaver et al., in prep.), which uses The Tractor software (Lang et al. 2016) to obtain more accurate photometry by fitting galaxy profiles with parametric models. The morphology of the sources is determined through a decision tree, separating point and extended sources. In contrast with the Classic approach, this method directly provides total magnitudes, performs an improved deblending in the high-resolution images, and extracts all the images consistently.

Galaxy selection criteria
Here we describe different steps taken to select a complete sample of galaxies at high redshifts. The resulting candidates are discussed in Sect. 4.
We exclusively search for candidates over the UltraVISTA area, by rejecting sources located in the masked regions near bright stars in the HSC images. This corresponds to an area of 1.404 deg 2 . We require the candidates to be detected in three bands among the H, K s , [3.6] and [4.5] bands. We require at minimum a 5σ detection in one of these bands, one at 3σ in a second band, and at 1σ in a third band. This ensures that at least two colours are reliable to estimate photo-z. Sources only detected in IRAC images are not included. We apply these selection criteria separately on both the Classic and The Farmer catalogues.
We reject candidates detected in any broad band blueward of the Lyman alpha break based on visual inspection, rather than explicit magnitude cuts. The measured flux may be inconsistent with zero in the case of a noise local maximum at the source coordinates. Moreover, nearby sources likely contaminate aperture fluxes for blended candidates, and generate a significant flux in the photometric apertures. In these cases, we estimate whether the observed flux originates from the detected source or from nearby sources. In the latter case, we flag these sources as blended. Artefacts are also rejected through the visual inspection of the science images, as discussed in Appendix A.

Photometric redshift selection
Our initial candidate lists are selected taking into account both the photometric redshifts and the posterior probability redshift distributions (zPDF) computed for all sources in COSMOS2020 using LePhare (Arnouts et al. 2002;Ilbert et al. 2006) with both the Classic and The Farmer catalogues. We first use the default COSMOS2020 configuration of LePhare, as done in Ilbert et al. (2013), Laigle et al. (2016) and Weaver et al. (2022), with the templates, dust attenuation curves and the recipe to add emission lines from Ilbert et al. (2009). The galaxy library includes 33 templates, covering various SED types, from quenched to starbursting galaxies. The most relevant ones for this work are 12 templates generated with the Bruzual & Charlot (2003) (hereafter BC03) stellar population synthesis models, with ages ranging from 30 Myr to 3 Gyr, and including sub-solar metallicities (Z = 0.004, 0.008, 0.02). We add the dust attenuation as a free parameter, allowing E(B − V) to vary between 0 and 0.5 for two different attenuation curves: Calzetti et al. (2000) and Prevot et al. (1984). Several emission lines are added using an empirical relation between the UV luminosity corrected for dust attenuation and H α emission line flux. Physically constrained ratios are considered between the intrinsic emission lines. We added [OII], H β , [OIII], and H α emission lines, as well as Lyman α , despite the large uncertainties that potentially affect the modelling of this line. Dust attenuation is applied to the emission line fluxes using the same dust model as for the stellar continuum. The normalisation of the emission line fluxes are allowed to vary by a factor of two (using the same ratio for all lines) during the fitting procedure. IGM absorption is implemented following the analytical correction of Madau (1995).
As in Laigle et al. (2016), we perform the fit using fluxes (and not magnitudes), even when a source is extremely faint or non-detected in one band. This type of approach is suitable as long as uncertainties are measured consistently. Therefore, we do not have to include upper-limits in the fitting procedure.
LePhare provides the redshift likelihood distribution for each source, after a marginalisation over the galaxy templates and the dust attenuation. We use it as the posterior redshift probability density function (zPDF), assuming a flat prior. The photo-z point estimate, z phot , is defined as the median of the zPDF. We also consider the photo-z which minimises the χ 2 over the full template library as an alternate photo-z point estimate, to strengthen our selection of high-redshift sources. We select any source with a solution at z ≥ 7.5 with LePhare, in either the Classic or The Farmer catalogue. This is defined by imposing that both the median of the zPDF and the redshift which minimises the χ 2 are at z ≥ 7.5.

Complementary template-fitting procedures
We produce photometric redshifts with other template-fitting procedures for the selected candidates. These estimates are not directly used to select the candidates, but are used in Sect. 4.1 to assess their robustness.
The photometric redshifts computed for all sources using EAZY (Brammer et al. 2008) are available in COSMOS2020. The adopted strategy is equivalent to LePhare and corresponds to the technique described in Weaver et al. (2022). The fitted galaxy library consists of 17 templates derived from the Flexible Stellar Population Synthesis models (Conroy et al. 2009;Conroy & Gunn 2010). Moreover, EAZY allows the combination of the templates in the fitting procedure. Furthermore, we compute new photometric redshifts for the selected candidates using a different LePhare configuration, hereafter noted LePhare BC03. Such a configuration is optimised to model sources with extreme colours. Since we select a really specific population with extreme properties, a larger coverage of the parameter space (in terms of dust, ages, star formation histories) should be explored to reject potential intermediate redshift contaminants. In this case, we include a set of BC03 templates assuming different star formation histories (exponentially declining and delayed), as described in Ilbert et al. (2015). Each of these 12 templates is generated at 43 ages from 50 Myr to 13.5 Gyr. During the fitting procedure, no template with an age older than the age of the Universe is considered. We assume two attenuation curves (Calzetti et al. 2000;Arnouts et al. 2013) with E(B − V) varying from 0 to 2. The dust attenuation reaches A V ∼ 8, enabling potentially extremely dusty sources at lower redshift to be rejected ). We add the emission line fluxes with a recipe described in Saito et al. (2020). For each template, we derive the number of ionising photons by integrating the SED blueward of the Lyman break, and convert it into H β luminosity following Schaerer & de Barros (2009). Then, the ratios between H β and other emission lines are given by Osterbrock & Ferland (2006); Anders & Fritze-v. Alvensleben (2003). Emission line fluxes are allowed to vary by a factor of two during the fitting procedure, to reproduce potential variations around the expected value. We do not use this configuration to compute photometric redshifts for the full COSMOS2020 catalogue, since the large parameter space covered by the templates increases the risk of degeneracy in the colour-redshift space and creates a larger fraction of catastrophic failures for the general population 1 .

Brown dwarf contamination
The selected high-redshift candidates may be contaminated by cool Milky Way brown dwarfs, because of their similar nearinfrared colours (e.g. Wilkins et al. 2014). Given the predictions by Ryan & Reid (2016), we could expect 277/deg 2 brown dwarfs at J < 25.
To isolate these contaminants, we fit brown dwarf templates to COSMOS2020 photometry. We include the simulated highresolution brown dwarf templates from Baraffe et al. (2015, BT-Settl/CIFIST2011_2015) and Morley et al. (2012Morley et al. ( , 2014 in LePhare. The modelled emission extends to at least 10 µm in the mid-infrared. In addition, the nonphysical templates are rejected following constraints from Saumon & Marley (2008) in the effective temperature versus surface gravity space, based on the predicted evolution of cool brown dwarfs. Then, we compare the minimum χ 2 obtained with the galaxy templates (we use the standard LePhare configuration described in Sect. 3.1) and the one obtained with the stellar library including brown dwarf templates. We compute ∆χ 2 = χ 2 star − χ 2 gal with χ 2 gal being the minimum χ 2 obtained with the galaxy templates, and χ 2 star the one obtained with the stellar templates. We select candidates with ∆χ 2 > 0.

The case of ID720309
The source ID720309 has the highest photo-z among the candidates extracted from the COSMOS2020 catalogue. However, we identified a bright star at the coordinates R.A.=09h59m10.81s, Dec.=+2d11m04.29s which could potentially generate a cross-1 Reducing the parameter space could be seen as a prior for the general population. Fig. 2: Images of the source ID720309 in the Y, J, H, K s UltraV-ISTA broad bands (from left to right), in the full stack (top) and in the stack including only the images not contaminated by the potential cross-talk artefact (bottom). The circles represent the 2 diameter aperture used in the Classic photometry.
talk signal (see Appendix A) at a position 09h59m10.81s +2d06m42.40s (assuming a native pixel scale of 0.341 expected for the detector 7 of the VIRCAM camera). Each ultradeep stripe is the combination of VIRCAM observations taken at three different declinations, each called "paw" (see McCracken et al. 2012). This source presents two components in one paw (paw 1), and just one component in the other (paw 2). It suggests that ID720309 in COSMOS2020 is the superposition of a cross-talk and a real source separated by 0.7 . It contributes to the signal in the 2 aperture of the Classic catalogue. As The Farmer photometry is centred on the real source, the cross-talk has less impact on the profile-fitting photometry. We also notice that the [3.6] and [4.5] positions are well aligned with the K s position of the source in paw 2.
The star responsible for the artefact is not included in 80% of the images taken in the paw 2. Therefore, we created new stacks in Y, J, H, K s using only these observations. This impacts the sensitivity of our dataset, but it ensures that the source is not affected by the cross-talk. Figure 2 shows the difference between the images on the full stack (top panels) and the stack not affected by the cross-talk (bottom). The main difference appears in H with a change in the shape of the source. We recomputed the flux using The Farmer on this new stack. The resulting magnitudes are given in Table 2, showing fainter magnitude in H and K s , and a non-detection in the Y-band.

Resulting galaxy sample
Our final galaxy sample comprises 32 candidates at z ≥ 7.5. Based on the LePhare/The Farmer photometric redshifts, the sample includes 15 candidates in the range 7.5 < z < 8.5, 11 candidates at 8.5 < z < 9.5 and 1 candidates at z > 9.5 (some photo-z are below z < 7.5 in the LePhare/The Farmer configuration). Tables 2 and 3 present the coordinates, photometry, photometric redshifts and absolute magnitudes of each object. The identifiers are from the Classic catalogue, and are always indicated starting with the letters "ID" in the following discussion. We report the magnitudes from The Farmer, corrected for the Milky Way extinction and systematic zero-point offsets. Appendix B describes the photometry, the best-fitting templates and the zPDF estimated with LePhare together with a detailed discussion for each candidate. Figure 1 illustrates the coordinates of the identified candidates over the COSMOS field. We emphasise an important aspect of our selection: four candidates are located in the UltraVISTA "deep stripes". These areas are one magnitude shallower than the ultra-deep stripes in Y JH. This makes it more difficult to meet selection criteria requiring an unambiguous three-band detection.
We identify 17 new unpublished candidates at z ≥ 7.5, many of which are located within the westernmost UltraVISTA ultradeep stripe (at low right ascension), as shown in Fig. 1. This region is fully covered with the "deep" HSC-SSP DR2 survey, about 0.5 mag shallower on the outer part of the field compared to the centre. In contrast, the Suprime-Cam data did not cover this part of the field, which made impossible the search of z > 7 galaxies in this region in previous works, because of the lack of a sufficiently deep optical photometry.
We find some differences between the total apparent magnitudes between the Classic and The Farmer photometric COS-MOS2020 catalogues, as shown in Appendix B. The Farmer catalogue has the advantage of photometric measurements performed consistently over the full wavelength range. Moreover, The Farmer measures fluxes more accurately in crowded fields. Therefore, we adopt this photometry by default for our UVLF analysis.

Robustness of the selected candidates
We produce several tests to establish the robustness of our candidates. Based on this step, we define a gold sample of 22 sources satisfying all our criteria as z > 7.5 galaxies (see Table 3). But we cannot exclude a possible contaminant for the remaining sources, based on these additional criteria.
First, we compare several photo-z estimates to assess the sensitivity of our selection to the template-fitting approach. Table 3 presents the six photometric redshift estimates for each candidate, computed using the three different template-fitting procedures (LePhare, EAZY and LePhare BC03) applied to both the Classic and The Farmer catalogues, as described in Sect. 3.1 and 3.2. All template-fitting codes consolidate our selection for the first 16 candidates listed in Table 3. They represent the core of our gold sample. We find some discrepancies for the remaining candidates, separated by categories in Table 3.
Blended galaxies: for four of these candidates, the three photo-z estimates with the Classic catalogue are at z < 3. These sources are flagged as blended in the images (flag B in Table 3), as a bright foreground source is clearly identified and contaminates the 2 aperture. Based on the Classic aperture photometry, all four of these candidates show a 3σ detection in at least one HSC optical band. The flux captured in the fixed aperture present a spatial offset. The Classic redshift probability distributions peak at z ∼ 2 for all candidates, although one of them presents a secondary z > 7 solution. In contrast, the majority of the zPDF weights are located at z > 7 with The Farmer catalogue. It is precisely in these cases that The Farmer photometry is expected to be more reliable than the Classic aperture photometry, since these candidates are well deblended thanks to the profile-fitting photometry. Therefore, we consider these estimates with The Farmer as robust and we include these candidates in our gold sample. We note that ID859061 falls in the CANDELS region and Stefanon et al. (2019) confirmed the lack of emission in the visible from HST data (see their Figure 5).
Low-redshift galaxies based on EAZY photo-z: for two sources, we find one solution at z < 3 using EAZY. Since EAZY produces a high-redshift solution in either the Classic or The Farmer catalogue, we decide to consider these candidates as robust and keep them in our gold sample.
Low-redshift galaxies based on LePhare photo-z: EAZY points towards a robust z > 9 solution with a single peak for ID1356755. However, all the solutions obtained with LePhare present a significant peak at z < 3, including more than 30% of the zPDF. Hence, we do not include this candidate in our gold sample. The configuration The Farmer/LePhare is our reference to derive the UVLF, but corresponds to a photo-z of 2.55 +0.19 −0.18 . Therefore, we do not indicate the corresponding absolute magnitude in Table 3 and this source is not used to derive the UVLF. For the source ID720309, all photo-z indicate a z > 9 solution with the original COSMOS2020 photometry. However, our photometry could be contaminated by a cross-talk artefact (see Sect. 3.4). We recompute the photo-z with LePhare using the new stack not contaminated by the cross-talk. The redshifts which minimise the χ 2 are 9.46 +0.05 −6.59 and 9.47 +1.53 −0.18 for the COS-MOS and BC03 templates, respectively. We obtain 2.89 +6.52 −0.43 and 9.46 +0.75 −1.6 using the median of the zPDF. This galaxy remains a high-redshift candidate in our sample, but with a significant peak in the zPDF at z ∼ 2.8.
Dusty star-forming galaxies at intermediate redshift: we checked that none of the selected candidates has H − [3.6] > 2 which would indicate a low-redshift galaxy with strong Balmerbreak or high dust obscuration . None of the candidates has coincident Spitzer/MIPS 24 µm emission nor at longer wavelengths in the COSMOS Super-deblended catalogue (Jin et al. 2018), and are thus unlikely to be contaminated by low-redshift dusty galaxies. We also checked that none of these sources were detected in SCUBA2 (Simpson et al. 2019) at less than 6 . Three of the candidates present a lower redshift z < 6 solution when fitted using the LePhare BC03 configuration. In this case, we allow for a possible attenuation reaching A V ∼ 8, which could reveal a possible solution associated with a dusty star-forming galaxies. These three candidates have also secondary peaks in their zPDF using other template-fitting configurations. They are located outside the ultra-deep HSC region, and outside the SPLASH coverage. Their mid-infrared photometry relies on S-COSMOS, which is about 1 mag shallower than in the centre of the field. The use of shallower data may explain the lower constraint on the zPDF resulting in multiple peaks. In addition, the width of the zPDF is systematically larger when using the Classic catalogue compared to The Farmer, because of larger flux errors for the same objects (as explained in Weaver et al. 2022). While our fits favour a high-redshift solution for these sources, we flag them as possible intermediate redshift contaminants and we do not include them in our gold sample.
Secondly, we also assess the robustness of the star-galaxy classification. We already rejected all candidates with a better χ 2 using the brown dwarf templates than with the galaxy templates. However, such a sharp cut does not quantify the risk of degeneracy in the classification. Therefore, we generate 500 realisations of the candidates, and for each realisation we add noise to the observed flux corresponding to the photometric uncertainty in each band 2 . We compute ∆χ 2 for each realisation. The analysis of these 500 realisations allows us to quantify the robustness of the star-galaxy classification for a given source. For three sources ID485056, ID545752, ID1346929, we find that 31%, 28%, 23% (27%, 33%, 43%) of the ∆χ 2 distribution falls at ∆χ 2 < 0 with The Farmer (Classic). Therefore, these sources present a significant probability to be brown dwarf Table 2: Coordinates and observed photometry of the selected z ≥ 7.5 candidates. The first columns indicate the ID and coordinates from the Classic catalogue. The other columns give the photometry from The Farmer catalogue, corrected for Milky Way extinction and systematic zeropoint offsets. The upper-limits correspond to the 3σ depth of the image at the source coordinates (10σ for IRAC), for sources with S/N < 1. contaminants, and we do not include them in our gold sample.
Thirdly, we inspect the HSC-SSP DR3 images in g, r, i and z for the 32 candidates. The increase in depth for DR3 compared DR2 is minor in the COSMOS field (Aihara et al. 2022). DR3 correspond to the final images taken on COSMOS and the sky subtraction has been improved in the SSP pipeline.
In the z-band image at the position of each of the 32 candidates, we measure the photometric flux in a 2 diameter aperture and the corresponding sky flux in a 2 − 4 diameter annulus using the python package photutils. To minimise the impact of objects in the sky annulus, we compute a sigma-clipped mean value. The noise per pixel is computed by aggressively detecting and removing all objects in the postage-stamp around the object and computing the standard deviation of the remaining pixels. We note that the precise signal-to-noise (hereafter S/N) values reported here are dependent on the exact details of object thresholding and background computation. However, we have verified that the conclusions presented below are largely robust to the exact parameter choice. In addition to the S/N values presented below, we examine both the pixels present in the measurement aperture and in the sky annulus around the object. Six candidates have S/N > 1. One of them is already identified as potential z < 6 dusty galaxy (ID241443). Two Table 3: Photometric redshifts estimated with LePhare and EAZY as in Weaver et al. (2022), and with LePhare using BC03 templates, for the z ≥ 7.5 candidates. The UV absolute magnitudes are derived from a Monte Carlo sampling of the galaxy zPDF, using The Farmer and LePhare. The three last columns list identifiers from S19 and B20, flags a and spec-z from Bouwens et al. (2022 .. a The flags are the following. 1: outside the HSC ultra-deep region, 2: outside the IRAC SPLASH region, B: blended, D: inside the UltraVISTA deep stripes. b LePhare photo-z computed using the UltraVISTA photometry uncontaminated by the cross-talk (see Sect. 3.4). c Absolute magnitudes computed using the ALMA spec-z. We report the error associated with the H-band apparent magnitude. galaxies (ID1297232 and ID441697) have a z-band S/N of 1.4 and 2.7, respectively, and are listed separately in Table 3. The Article number, page 8 of 25 O. B. Kauffmann et al.: COSMOS2020: UV-selected galaxies at z ≥ 7.5 three remaining S/N > 1 candidates are not considered since the bright pixels are shifted by 1 with the expected source position, or the signal comes from only one pixel.
To summarise, the first 16 candidates satisfy all the criteria to be selected as z > 7.5 galaxies. Six additional candidates present discrepancies between the different photo-z procedures, explained by blending in the Classic catalogue and one discrepant photo-z in the EAZY run. We still consider them as robust and include them in our gold sample. For the last 10 candidates in Table 3, there is a significant probability that these sources are either intermediate redshift galaxies or brown dwarf contaminants.  Table 3. One of them was independently confirmed using a detection of the Lyman-α emission line (Valentino et al. 2022). Only one ALMA source (ID564423/REBELS-24) does not present a [CII] 158µm line detection, either because of too weak SFR, or observations which are still missing for this source (Bouwens et al. 2022). These five spec-z confirm the high-redshift nature z > 7 of our candidates. We find a systematic over-estimate of the photo-z. For the Classic catalogue, these spec-z are consistent with the zPDF (with a PIT 4 of 0.3, 0.05, 0.10, 0.13 for ID336101, ID978062, ID1103149, ID1151531, respectively). For The Farmer catalogue, their PIT is below 0.02 and the uncertainties for these five sources are underestimated. While this comparison could suggest a systematic over-estimate of the photo-z, we note that these sources could be really specific, which could explain their detection of [CII] 158µm with ALMA. Therefore, we do not extrapolate this trend to the full population at this stage. In Sect. 5, we replace the photo-z of these five sources by their spec-z.

Comparison with
We compare our candidates to those previously identified by S19 and B20 in the COSMOS field based on ground-based imaging. S19 (see also Stefanon et al. 2017) used the near-infrared broad and narrow bands from UltraVISTA DR3, all the available CFHT/MegaCam, Subaru/Suprime-Cam optical bands and Spitzer/IRAC channels 1 to 4 in the mid-infrared. This study also benefited from HST/WFC3 coverage from the Drift And SHift mosaic (DASH; Momcheva et al. 2016;Mowla et al. 2019), with an improved spatial resolution. The HSC imaging was not available at the time of their initial sample selection but the authors used the HSC-SSP DR1 data to validate their candidate list, checking that there was no significant detection in the optical. S19 identified 16 galaxy candidates in the COSMOS field, including 10 at z ∼ 8 and 6 at z ∼ 9. B20 included the UltraVISTA DR4 data and the Spitzer/IRAC [3.6] and [4.5] images from SPLASH, SEDS, and SMUVS in the infrared. Optical data consisted of the CFHT/MegaCam u * , g, r, i, z broad bands from CFHTLS, and the HSC-SSP DR1 g, r, i, z, y broad bands. The Suprime-Cam z band, deeper than the HSC/z band in DR1 release, was also 3 https://sites.google.com/view/a3cosmos/data 4 Probability Integral Transform (PIT, Dawid 1984) defined as PIT = zs 0 P(z) dz. used. The search was performed in the HSC ultra-deep area. B20 recovered seven candidates from S19 and selected nine new candidates for a total of 16 LBG at 7.5 < z < 9.1 in the COSMOS field, including 14 at z ∼ 8 and 2 at z ∼ 9.
We select 15 out of the 25 high-redshift candidates from S19 and B20. The majority of these candidates present single-peaked redshift probability distributions located at z > 7, and are classified as galaxies in both COSMOS2020 catalogues. The identifiers from S19 and B20 are shown in Table 3. The 10 rejected candidates from S19 and B20 are described in details in Appendix B.3. Four of those candidates include strong low-redshift solutions, and six are not detected in the combined izY JHK s image. Since these galaxies are not expected to be visible in the i and z bands, the signal may be diluted in the combined izY JHK s detection image. To test the impact of such approach, we matched these six sources with the official release catalogues from UltraVISTA DR4 (McCracken et al. 2012) with a detection performed in each VIRCAM band. These sources do not have any counterpart in the K s catalogue. One source from S19 is detected only in the H selected catalogue (but not identified in COSMOS2020 or B20 because of the proximity of a bright source J = 20.7). Another source from S19 is only detected in the J-band selected catalogue with a magnitude of J = 26.3 ± 0.1. At best, we could have retrieved one of these six candidates with a detection performed in each VIRCAM band, rather than our combined izY JHK s detection image.
Finally, we checked the z ≥ 7 galaxy candidates identified in CANDELS (Bouwens et al. 2015b). None of these sources are in our sample considering our selection criteria at z ≥ 7.5. Although deeper, CANDELS covers a smaller area of 151.9 arcmin 2 , which is less efficient when selecting the brightest sources at z > 7.5.

Lensing magnification
Any high-redshift galaxy sample may be subject to gravitational lensing, and in particular lensing magnification, from massive, low-redshift galaxies (e.g. Mason et al. 2015;Barone-Nugent et al. 2015;Roberts-Borsani et al. 2016). While gravitational lensing preserves surface brightness, the apparent solid angle of the background source may increase, leading to an increased apparent flux.
We investigate the possible impact of lensing magnification on our selected galaxy candidates. We search the full COS-MOS2020 catalogue for massive low-redshift galaxies within a 20 radius for each candidate galaxy 5 . Lens galaxies are modelled as singular isothermal spheres (SIS). Stellar velocity dispersions are estimated from photometry through the Faber-Jackson relation (FJR; Faber & Jackson 1976) of Barone-Nugent et al. (2015), based on the rest-frame B-band absolute magnitude and calibrated using early-type galaxies with redshifts spanning 0 < z < 1.6. Velocity dispersion uncertainties are dominated by the intrinsic scatter in the FJR, estimated at 46 km s −1 . We checked that these velocity dispersion estimates are consistent with the FJR of Bernardi et al. (2003) using i-band absolute magnitudes. Photometric redshifts and absolute magnitudes in the rest-frame Suprime-Cam/B are taken from The Farmer catalogue. We restrict the lens selection to galaxies with a velocity dispersion of at least σ v = 200 km s −1 , because the spectroscopic samples used to calibrate the FJR become incomplete at lower A&A proofs: manuscript no. output values (Barone-Nugent et al. 2015). In addition, we only include lenses with a magnification of µ ≥ 1.1.
We find that eight candidates out of 32 are probably magnified with 1.1 < µ < 1.2, and 5 galaxies with µ ≥ 1.2. This includes five already identified candidates from S19 and B20. However, we find no evidence of strongly lensed galaxies with multiple images. The most magnified candidate is ID441697, with a cumulative magnification of µ = 2.36 ± 0.80 from five lenses within 14 , representing a boost of 0.9 mag. The main contribution comes from a z = 0.50 +0.01 −0.01 galaxy located within 4.6 and with a σ v = 231 km s −1 velocity dispersion, leading to a µ = 1.37 ± 0.20 magnification. The second most magnified candidate is ID442053 because of two sources, one z = 1.50 +0.03 −0.04 galaxy at 6.0 with σ v = 218 km s −1 and µ = 1.14 ± 0.07, and one z = 1.74 +0.04 −0.04 galaxy at 4.0 with σ v = 207 km s −1 and µ = 1.17 ± 0.09. The candidate Y8 (ID1209618) has a z = 1.26 +0.02 −0.02 galaxy at a 4.4 distance with a velocity dispersion of σ v = 215 km s −1 , which gives a µ = 1.21±0.11 magnification. This candidate was already identified as magnified by S19, who found a 40 km s −1 higher velocity dispersion and a 0.3 smaller angular separation 6 , resulting in a higher magnification µ = 1.39. Three other candidates (Y1, Y10, Y12) from S19 are moderately magnified, with 1.1 < µ < 1.2.
Consequently, every object in our z ≥ 7.5 sample is potentially affected by lensing magnification. However, as we shall see in Sect. 6.2, the impact of this magnification remains limited on our luminosity function measurements.
The velocity dispersion values quoted above are extracted from the FJR and their uncertainties reach 46 km s −1 due to the scatter of the relation ( Barone-Nugent et al. 2015). Moreover, the velocity dispersion may be overestimated, because of the calibration based on early-type galaxies, leading to overestimated magnifications. Inversely, we note that the velocity dispersion lower limit of σ v = 200 km s −1 significantly restricts the number of lens galaxies, so that the reported total magnifications may be underestimated. This will need to be further investigated with more precise velocity dispersion estimates at fainter luminosities.

The UV luminosity function
In this section, we use our z ≥ 7.5 candidates to make a new measurement of the bright end of the UVLF. We select objects brighter than H = 25.6, equivalent to the 5σ depth in the ultradeep stripes. We use the 0.812 deg 2 of the HSC-masked Ultra-VISTA ultra-deep area, so we do not include the 4 candidates selected in the deep area. This ensures a homogeneous selection function across the field. We take the physical parameters estimated with the LePhare/The Farmer configuration. This is necessary because the blended candidates are not included in the Classic selection.
We split our galaxy sample into three redshift bins, centred at z = 8, 9, 10 with a ∆z = 1 width. Since the photometric redshift estimates mainly rely on infrared broad-band imaging, the redshift probability distributions are relatively broad in the interval 7 < z < 10, and so are the photometric redshift uncertainties. As a result, the candidates with photometric redshifts at the limit between two adjacent bins may be scattered in one bin 6 We find an astrometric shift of 0 . 34 in declination compared to S19 for the candidate Y8 (ID1209618). This is more than the shift of about 0 . 15 induced by the change of astrometric calibration from UltraVISTA DR3 to DR4. We note that this candidate remains faint with H > 25.6, which may have led to uncertain coordinates in S19 who used DR3. or the other. Therefore, we develop a Monte Carlo simulation to propagate the photometric uncertainties into the absolute magnitude estimates. We generate 500 realisations of the high-redshift candidates. For each realisation, we add noise to the observed flux, according to the photometric uncertainty measured in the considered band. We base our estimate of the UVLF on the 500 catalogues for which the photo-z and absolute magnitudes have been derived.
We adopt the spec-z derived with REBELS, when possible. As a result, two sources among the five are included in our first redshift range, for which we recompute the FUV absolute magnitudes.

Completeness correction
To estimate the completeness of our sample of high-redshift galaxies, we simulate point-like sources with a Moffat profile (Moffat 1969) with parameter β = 3 and a FWHM of 0.9 . We generate these objects with a uniform distribution in magnitude over 22 < H < 27 and in position over the ultra-deep stripes area and add them to the J, H, K s bands, avoiding placing objects on masked areas or where object are present in the segmentation map derived from the CHI_MEAN image (see Sect. 2.2). We then generate a new CHI_MEAN detection image from these simulated J, H, K s images combined with the original i, z, Y images. Object detection is repeated using the same parameters as the main catalogue and using this new CHI_MEAN image. The resulting catalogue is cross-correlated with the input simulated catalogue to the find the fraction of recovered sources as a function of magnitude.
We must assume colours for the simulated galaxies. We use LePhare and the BC03 templates to predict the expected colours. For E(B − V) = 0.1, we find average colours of J − H = 0.05, 0.49, 4.94 at z = 8, 9, 10, respectively, and K s − H = 0.00, −0.10, −0.09 at z = 8, 9, 10, respectively. The predicted K s − H colours do not vary by more than 0.25 mag at z > 7, depending on the assumed attenuation. We simplify the simulation by assuming J = H = K s for galaxies at z ∼ 8, and H = K s for galaxies at z ∼ 9 and 10, with no flux contribution in the other bands.
We estimate the completeness as a function of the H-band magnitude, defined as the fraction of recovered sources in the CHI_MEAN image. At z ∼ 8, we find a drop in completeness from 84% to 60% between H = 25 to 25.6, the latter being our selection limit. At z ∼ 9 and 10, the completeness drops from 72% to 37% in the same magnitude range.
We do not attempt to correct for contamination or incompleteness in our selection method, considering that possible biases in the redshift estimate (as shown by the comparison with ALMA) dominate the uncertainty budget and cannot be captured by a simulation. Moreover, the use of the zPDF to generate multiple realisations of the UVLF would make such correction difficult, since a single source is split in several redshift bins.

Binned luminosity function
Absolute UV magnitudes are computed with LePhare in the GALEX far-UV filter 7 , as follows: with m f being the observed magnitude in the filter f , DM the distance modulus. KC is the sum of the k-correction k f and a rest-frame colour term, all derived from the best-fit SED. The filter f is chosen among J, H, K s , [3.6] to minimise the SED dependency of the KC term (Ilbert et al. 2005). For each source, we list in Table 3 the median of the M UV distribution, as well as the associated 0.16 and 0.84 quantiles. The binned luminosity function is calculated using the V max estimator (Schmidt 1968). This estimator is non-parametric, although the number of bins and the bin widths are set, and implicitly assumes a uniform spatial distribution of galaxies. The number density in a given magnitude bin depends on the maximum volume V max in which each galaxy could have been selected. This volume, for a given galaxy i, is computed as: where z min,i and z max,i are the lower and upper redshift limits in which a galaxy i can be included in the sample, and dV is the differential co-moving volume. The co-moving volume is defined as the shell between the limits of the redshift bins. The maximum redshift z max,i is the redshift at which a galaxy with the same intrinsic properties would not be selected in our sample, with a flux limit set to H = 25.6. Thus, the luminosity function φ can be expressed as: where N(M) is the number of galaxies in the magnitude bin centred at M and of width ∆M, considering all the N real = 500 realisations of the catalogues. w i is the inverse of the completeness estimated at the H-band magnitude of the galaxy i. The associated Poisson uncertainties σ φ are computed as Marshall (1985): The total uncertainties are computed as the quadratic sum of the Poisson errors and cosmic variance errors, estimated following Trenti & Stiavelli (2008). Since the galaxy samples become incomplete at faint magnitudes, the LFs are computed brightward of M UV = −21.5. The galaxy UV luminosity functions at z = 8, 9, 10 estimated from the selected galaxy candidates are represented in Fig. 3, 4, and 5, and tabulated in Table 4. In magnitude bins with no galaxies, we put an upper limit with the number density computed from one galaxy. Cosmic variance represents about 14 % of the Poisson uncertainties at M UV < −23, 20 % at M UV = −22.75, and 30 % at M UV = −22.25.
We show in Fig. 6 the three UVLFs. Given the magnitude limit at H < 25.6, we cover the brightest end of the UVLF at −23 < M UV < −21.5. There is no clear evolution of the number densities from z = 8 to z = 9, with an equivalent number of candidates at −22.5 < M UV < −21.5 in both redshift bins. At z = 10, the density in the single absolute magnitude bin is lower than at z = 9, but the median M UV is shifted to brighter magnitudes and could follow a simple extrapolation of the LF at lower redshift. However, the constraint on the UVLF at z ∼ 10 remains weak.
We provide a constraint on the brightest part of the UVLF and over 1.5 mag range, thus the shape of the UVLF cannot be Table 4: Galaxy luminosity functions derived from all the z ≥ 7.5 candidates. The columns indicate redshifts, central absolute magnitudes, magnitude bin widths, median absolute magnitude in the bin, number of galaxies averaged over the 500 realisations (see Sect. 5.2) and co-moving number densities.  constrained from our data alone. Hence, we do not attempt a fit with a power law or a Schechter function. Nevertheless, we discuss the consistency with the various fitting forms when comparing with the literature in Sect. 5.3.

Comparison with UVLFs from the literature
The calculated UVLFs are in good agreement with results from the literature in particular at M UV < −22. This suggests that our sample is complete for the brighter magnitude bins. The most straightforward comparison is with B20 who used the same nearinfrared images as this work (see Sect. 4.2). In the COSMOS field, we have 12 galaxy candidates in common. However, the authors analysed about 6 deg 2 of imaging data from both the COSMOS and the XMM-LSS fields, sampling larger co-moving volumes than in this work. All candidates from B20 at z < 8.5 and M UV < −22.5 are outside the COSMOS field, and this is the reason for their better constraint on the bright end of the UVLF. Nonetheless, our number densities at M UV < −22.5 are in excellent agreement with B20 and S19. At redshift z ∼ 10, the number densities computed from the three candidates selected in this work are in agreement with the double power-law evolution from B20. The candidate XMM3-3085 from B20, identified in the XMM-LSS field with a photometric redshift of z phot = 10.8 ± 1.0, is extremely bright with an absolute magnitude of M UV = −23.7 and H = 23.9. It is the brightest z > 7 galaxy candidate ever found in the literature, although spectroscopic confirmation is still required. The z ∼ 10 candidate 2140+0241-303 from Morishita et al. (2018) has an HST/F160W flux of 24.4 mag and an absolute magnitude of M UV = −22.6. The authors used the Brightest of Reionising Galaxies (BoRG[z9]) survey, including HST optical and nearinfrared imaging in five broad bands over 370 arcmin 2 , in addition to IRAC/[3.6] imaging. The resulting number density from that paper is an order of magnitude higher than our results at M UV < −22.5. We compare our results with Bouwens et al. (2021), who used the most comprehensive compilation of HST data taken on deep fields. Our results at z ∼ 8 and 9 are consistent within the uncertainties. The faintest of our points at z ∼ 8 falls below their fit by a factor two in density. At z ∼ 10, they considered the results from Oesch et al. (2018) which indicate a sharp decline of the galaxy density at M UV < −22. Our data probe brighter absolute magnitudes (our only point is brighter than M UV < −22 while their brightest point is at M UV = −21.25). However, the lower density of our fit is unconstrained and consistent with zero. Therefore, we cannot draw clear conclusions regarding differences with Oesch et al. (2018).

Sources of contamination
The most obvious explanation for the high density of bright sources at z ≥ 7.5 is the contamination by low-redshift galaxies or brown dwarfs wrongly classified as z > 9 candidates. Even if we combine several photometric redshift codes, and two different methods to extract the photometry, contamination of our sample remains a possibility. In Sect. 4.1, we found 10 candidates with a significant probability of being either dusty starforming galaxies at intermediate redshift or brown dwarf contaminants. Seven over ten of these candidates fall on the westernmost ultra-deep stripe with a lower HSC and IRAC coverage (see Fig. 1). The larger surface density of sources in this stripe compared to the others points out towards a significant population of contaminants among these ten galaxies.
In Fig. 3, Fig. 4 and Fig. 5, we present the UVLF after having removed these sources (gold squares). We find that the measure-ment at 7.5 < z < 8.5 is severely affected by this removal, as expected since 7 sources have a photo-z estimated in this redshift range. Therefore, the density could be affected by a factor 2 if all these sources are in fact contaminants.
A systematic overestimate of the photo-z could also impact the UVLF. The 5 sources with a spec-z have a systematically higher photo-z (see Sect. 4.1). If this trend affects the whole sample, galaxies could move in a lower redshift bin. Given the shape of the redshift distribution, it would lead to an overestimation of the UVLF density. A better modelling of the galaxy properties (e.g. dust attenuation law, emission lines) could alleviate these biases. We are particularly sensitive to the modelling because the major feature used to derive the photo-z at this epoch, namely the continuum break around 1216 Å, cannot be accurately located. Indeed, this break is redshifted within the gap between the ground-based Y and J-bands for galaxies between z ∼ 7.6 and z ∼ 8.7. Spectroscopic confirmation of the candidate redshifts is the solution to alleviate these uncertainties (Muñoz & Loeb 2008), as well as future space mission without such gap between near-infrared filters.
In addition, the uncertainties in the absolute UV magnitudes may affect the bright end through the Eddington bias (Eddington 1913). Because of the steep slope of the luminosity function, there are statistically more faint galaxies scattered in the brighter bins than the reverse, resulting in a flattened slope. The Eddington bias is also stronger for steep luminosity functions. In the selected galaxy sample, photometric redshift uncertainties can be relatively large, leading to large uncertainties once propagated to absolute magnitudes. This bias may be limited by using large magnitude bins.
The presence of Active Galaxy Nuclei (AGN) in the selected galaxy sample may also affect the estimated galaxy UVLF at the bright end. High-redshift AGN have Lyman alpha break features similar to star-forming galaxies without an AGN component. At intermediate redshifts 4 < z < 6, the contribution from faint AGN dominates the number densities at M UV < −23, whereas it becomes negligible at M UV > −22 (Ono et al. 2018;Harikane et al. 2022). At z > 6, the number density of faint AGN is still uncertain. The evolution of the quasar spatial density is often parametrised as ρ(z) ∝ 10 k(z−6) , with k −0.47 from z = 3.5 to z = 5 (Fan et al. 2001) and k −0.72 from z = 5 to z = 6 (Jiang et al. 2016). Recently, Wang et al. (2019) measured a consistent value k −0.78 from z = 6 to z = 6.7. With this accelerated redshift evolution, high-redshift AGN are sufficiently rare such that they have a negligible impact on the galaxy number density at M UV = −23 in a survey of this size. The faint-end slope of the quasar UV luminosity function is nonetheless poorly constrained at high redshift (Matsuoka et al. 2019).

Impact of magnification
As a consequence of gravitational lensing at very high redshifts, and in particular for steep luminosity functions (Mason et al. 2015), the bright end of the luminosity function is expected to be artificially higher. As discussed in Sect. 4.3, we find no evidence of strongly magnified galaxies in the selected sample; although, candidates are still affected by multiple lenses, leading to cumulative magnifications up to µ = 2.4. We estimate the UVLF after having corrected the UV absolute magnitudes from the magnification (not attempting to correct the density). The binned UVLF at z = 8 and z = 9 remain unchanged. Thus, magnification bias does not explain the lack of evolution of the UVLF bright end at the probed magnitudes. In the z = 10 galaxy sample, the can-didate ID441697 with a cumulative magnification of µ = 2.36 becomes fainter than M UV = −22 after removing the effect of lensing. Therefore, it would even lower down the density found at z ∼ 10.

Shape of the UVLF
It is well established that the UVLF at z < 6 is well-fit with a Schechter function (e.g. Moutard et al. 2020;Bouwens et al. 2021). At the same time, the very different shape of the bright and massive ends of the dark matter (DM) halo mass function from cosmological models and the observed galaxy luminosity function is currently explained with "quenching", or any process which can halt star formation in these haloes and galaxies. This cessation of the star-forming activity can be due to numerous processes like AGN feedback (e.g. Croton et al. 2006;Hopkins et al. 2006) or halo quenching due to shock-heated gas in halos more massive than 10 12 M (e.g. Somerville et al. 2008;Gabor & Davé 2015). Peng et al. (2010) showed that if the quenching rate is proportional to the SFR, it would naturally lead to a galaxy luminosity function with a Schechter shape.
While the Schechter form is well established at z < 6, whether this holds at higher redshifts is unclear. Bowler et al. (2015) find a high density of sources at M UV < −22 and conclude that the UVLF at z > 8 is better fit with a power law than a Schechter function, without displaying an exponential cutoff at the bright end. We confirm the density measured by B20 and our points are superposed on their extrapolation of the UVLF with a power law. Our measurements are inconsistent with the Schechter function obtained by McLeod et al. (2016), and the density of galaxies we measure is ten times higher than expected by their extrapolation to bright magnitudes (or brighter by 0.4 mag). The brightest galaxies observed by McLeod et al. (2016) is at M UV = −20.7, not allowing them to constrain the UVLF in the bright regime probed by our data.
The latest estimates by Bouwens et al. (2021) compile all HST measurements and find that the UVLF are consistent with a Schechter function at z ∼ 8 − 9. The new Schechter parameters obtained by Bouwens et al. (2021) shift the exponential cut-off at brighter magnitude than previous publications (e.g. McLeod et al. 2015). This is more consistent with the density of galaxies we find at M UV < −22. At z ∼ 8, our faintest point falls below their fit; however, such an effect could be explained by some residual incompleteness that is not corrected in our measurement. We find that our points are still consistent with their fit at z ∼ 9. So, we cannot reject the Schechter fit by Bouwens et al. (2021).

Lack of evolution at the bright end
The UVLF evolves rapidly at magnitudes fainter than M UV > −20.5. At a given density of φ ∼ 10 −4 mag −1 Mpc −3 , we expect a brightening of about 0.7 mag between z ∼ 9 and z ∼ 8 following the latest compilation from Bouwens et al. (2021). Between z ∼ 9 and z ∼ 10, this trend is even more extreme with a brightening of 0.9 mag. This evolution is interpreted as the galaxy population following the growth of the dark matter halos with a constant star-formation efficiency (Oesch et al. 2018).
In contrast with the faint end, the density of galaxies at M UV ∼ −22 is consistent with no evolution from z ∼ 8 to z ∼ 9 as shown in Fig. 6. This is consistent with the findings of B20.
One interpretation is that quenching efficiency increases between z = 9 and z = 8 at high mass and as a consequence we observe the building of the exponential cut-off of the UVLF. Peng et al. (2010) introduced a quenching rate proportional to the SFR to preserve a constant characteristic stellar mass around 10 10.6 M . If we extrapolate the relation between M UV and stellar mass found by Stefanon et al. (2021), we do not reach this regime at z = 8. Also, halo quenching is efficient above 10 12 M (Cattaneo et al. 2006). According to Stefanon et al. (2021), galaxies in our sample are hosted by lower mass halos. So, mass quenching is not expected to be efficient yet. However, the physical conditions at z > 8 may be different from those considered in lower redshift studies and the quenching mechanisms may be already effective for lower mass galaxies or DM halos.
Another possibility is that dust is generated in significant quantities between z = 9 and z = 8 (Finkelstein et al. 2015), sufficiently to decrease the light emerging from the brightest UV galaxies at z = 8, leading to the wrong interpretation that the density of bright galaxies remains the same between z = 9 and z = 8. Indeed, if dust is increasing with decreasing M UV , it would bend the bright end of the luminosity function, with a more pronounced dimming of the UV absolute magnitudes at bright magnitudes 8 . The mean dust content is generally expected to decrease with increasing redshift in particular for bright galaxies (e.g. Bouwens et al. 2016), because of the low metallicity of young stars in galaxies. Hence, the mean dust attenuation at z > 9 is often assumed to be zero (e.g. Bouwens et al. 2015b), so that the bright end of the UVLF is expected to reflect the recent star formation in the galaxy population.

Summary and conclusions
This paper presents one of the first science results from the COS-MOS2020 catalogues: a search for candidate z ≥ 7.5 galaxies at the epoch of the reionisation of the Universe. The deep optical, near-infrared, and mid-infrared imaging over the 1.4 squaredegree field of COSMOS enables the detection of rare and bright galaxies, a complementary population to the much fainter galaxies found in deep pencil beam surveys such as CANDELS or HUDF.
COSMOS2020 uses the latest and deepest optical, midinfrared, and near-infrared imaging in COSMOS. It includes the most recent UltraVISTA DR4 data together with new optical data from the HSC-SSP PDR2 release. We also used new images comprising all of the Spitzer/IRAC [3.6], [4.5] bands' data ever taken on COSMOS from the Cosmic Dawn Survey. From these data, we extracted two photometric catalogues using independent approaches: the Classic catalogue measures colours using circular aperture photometry, whilst The Farmer catalogue colours are measured based on a fit of surface brightness models. Sources were detected in a combined izY JHK image.
The galaxy selection was primarily based on photometric redshifts estimated from SED fitting with LePhare and their associated zPDF. We selected sources with a robust solution at z ≥ 7.5 with LePhare using either the Classic or The Farmer catalogues. The final sample consists of 32 candidates at z ≥ 7.5, including 17 unpublished candidates. To assess the robustness of the COSMOS2020 photo-z estimates, we computed new sets of photometric redshifts with three different template-fitting procedures using LePhare and EAZY, with both the Classic and The Farmer catalogues. We isolated a gold sample of 22 out of 32 candidates, which is more robust. Among them, four blended candidates were identified thanks to the profile-fitting photometry from The Farmer, where the contamination of light from nearby sources can be identified. These candidates have low photo-z solutions of z < 3 with the Classic catalogue because of the optical flux of the nearby sources. This illustrates the effectiveness of profile-fitting techniques in deblending confused sources typically found in the search of distant galaxies. The final sample of high-redshift galaxies is therefore more complete.
From this unique list of z ≥ 7.5 star-forming galaxies, we made a new determination of the bright end of the UV luminosity functions in three redshift bins, centred at z = 8, 9, 10. There is no clear evolution of the number densities from z = 8 to z = 9 between −23 < M UV < −21.5, which is in excellent agreement with B20. One interpretation is that quenching efficiency increases between z = 9 and z = 8. Another possibility is that dust is generated in a significant quantity, which is sufficient enough to decrease the light emerging from the brightest UV galaxies at z = 8. Another explanation for the high density of bright sources at z ≥ 7.5 is the contamination by lowredshift galaxies or brown dwarfs, or a systematic overestimation of photometric redshifts. Spectroscopic confirmation is thus essential. Follow-up observations with JWST-NIRSpec (Weaver et al. 2021) are planned for five of our brightest candidates (see Table 3). We have also begun a systematic follow-up of these candidates with the WERLS Key Strategic Mission Support Program on Keck (Casey et al., in prep.).
This work has demonstrated the great potential of COS-MOS2020 for finding rare, luminous objects in the distant Universe and the benefit from using multiple photometric extraction techniques and photometric redshift codes to assess the robustness of the results. In the longer term, deep Euclid near-infrared observations of COSMOS and other Cosmic dawn survey fields, together with our approved 218+80 h, 0.6 deg 2 COSMOS-Webb proposal, will provide the ultimate survey of bright, luminous objects in the early Universe.

Appendix A: Artefacts
The cross-talk effect is an electronic artefact in which bright sources reappear at different locations in the detector focal plane (e.g. Bowler et al. 2017). Such artefacts are particularly problematic for high-redshift galaxy searches as they appear in the near-infrared channels only. We identified several highredshift candidates which are most probably these inter-channel cross-talk artefacts in the UltraVISTA images. They are caused by bright sources located at the same RA with a Dec differing by ±k × 128 × 0 . 339 with integers k = 1, ..., 15 in the VIRCAM images. Here, 128 is the number of detector pixels in each of the 16 readout channels of the detector, and 0 . 339 is approximately the size of the original pixels. These sources are typically 11.9 − 13.8 mag brighter than the artefacts, and were identified comparing single pawprint stacked images with similar total exposure times, in the UltraVISTA bands. Moreover, we find no evidence of an optical counterpart with HSC for these sources, or in the mid-infrared with IRAC. The coordinates of these candidates are reported in Table A.1. Figure A.1 shows the UltraVISTA stamps at the corresponding coordinates.

Appendix B: Detailed description of each candidate
From Fig. B.1 to B.6, we show the observed photometry of each candidate, together with the best-fitted galaxy and stellar templates, for both the Classic and The Farmer catalogues. In the case where the flux is smaller than the flux uncertainty, for clar-ity the photometric measurement is replaced by a 3σ upper limit. The Classic photometry and posterior redshift distributions are also shown. Figures B.2, B.3, B.5, B.7 and B.8 display stamp images centred on the candidate coordinates. We note that there may be some offsets between the total apparent magnitudes between the Classic and The Farmer photometric COSMOS2020 catalogues, even though the observed colours match each others. One of the main reasons for this is that IRACLEAN and The Farmer provide total fluxes, whereas fixed aperture fluxes are used in the Classic catalogue. Hence, these aperture fluxes are rescaled using aperture-to-total corrections applied to all the aperture-extracted bands, and computed from the weighted mean difference between fixed aperture and pseudo-total fluxes (using Kron apertures) from SExtractor. While the corrections remain low (-0.29, -0.11, -0.01 mag for the 25%, 50%, 75% percentiles), this procedure introduces some additional noise for faint sources.
In IRAC, the flux is already total in both catalogues. Still, we find inconsistencies between the IRACLEAN and The Farmer photometry (e.g. for ID720309 and ID1103149), partially explained by the two different approaches in the algorithms. The Farmer assumes a parametric light-profile for a galaxy, convolved with the IRAC PSF before the fit. IRACLEAN repeatedly removes point-like source contributions from the residual map, until reaching a threshold with no pixel having a flux above a given signal-to-noise within the detection area, defined by the segmentation map. We identify that the presence of closeby objects induces large differences between the two sets of IRAC photometry. The different deblending procedures in The Farmer and Classic introduce small differences in the segmentation map, resulting in this case in large variations of the IRAC flux. We are not able to conclude the superiority of one catalogue over the other, highlighting the importance of using several methods to assess the robustness of the results.
in this area allows us to select this candidate in a region which has already been studied by S19 and B20.
The observed photometry and images of best-fitted galaxy templates of two blended candidates ID442053 and ID1313521 are shown in Fig. B.4 and Fig. B.5, respectively.These candidates have magnitudes 25.0 < H < 25.5. For these objects we find a photo-z solution at low redshift z < 3 with the Classic catalogue. From aperture photometry alone, all of these candidates have a 3σ detection in at least one HSC optical band. The Classic redshift probability distributions peak at z ∼ 2 for all of them, although one shows a secondary z > 7 solution. In contrast, the majority of the zPDF weights are located at z > 7 with The Farmer catalogue. .1 for the z ≥ 7.5 candidates from S19 and B20 which are recovered in the COSMOS2020 catalogue. The identifiers from those papers are indicated in parentheses.