We present the Lyman Ly and ultraviolet UV luminosity function LF in bins of redshift for quasars selected in the Physics of the Accelerating Universe Survey PAUS. A sample of 915 objects was selected at 2.7 lt z lt 5.3 within an effective area of 36 degSUP2SUP observed in 40 narrowband NB filters FWHM  120 . We cover the intermediatebright luminosity regime of the LF inlineformula idFI1texmath idtexeq1 1043.5 lt LmathrmLyalphamathrmergs1 lt 1045.5 texmathinlineformula 29 lt MSUBUVSUB lt 24. The continuous wavelength coverage of the PAUS NB set allows very efficient target identification and precise redshift measurements. We show that our method is able to retrieve a relatively complete C  85 and pure P  90 sample of Lyemitting quasars for LSUBLySUB gt 10SUP44SUP erg sSUP1SUP. In order to obtain corrections for the LF estimation and assess the accuracy of our selection method we produced mock catalogs of 0 lt z lt 4.3 quasars and galaxies that mimic our target population and their main contaminants. Our results show a clear evolution of the Ly and UV LFs with a declining tendency in the number density of quasars toward increasing redshifts. In addition the faintend powerlaw slope of the Ly LF becomes steeper with redshift suggesting that the number density of Lybright quasars declines faster than that of fainter emitters. By integrating the Ly LF we find that the total Ly emitted by bright quasars per unit volume rapidly declines with increasing redshift being subdominant to that of starforming galaxies by several orders of magnitude by z  4. Finally we stack the NB pseudospectra of a visually selected golden sample of 591 quasars to obtain photometric composite SEDs in bins of redshift enabling us to measure the mean intergalactic medium absorption using the Lyman forest as a function of redshift yielding results consistent with previous spectroscopic determinations.
['line: identification', 'methods: observational', 'galaxies: high-redshift', 'galaxies: luminosity function', 'mass function', 'quasars: emission lines', 'Astrophysics - Astrophysics of Galaxies']
The PAU Survey The quasar Ly and UV luminosity functions at 2.7 lt z lt 5.3
{'The PAU Survey: The quasar Ly α and UV luminosity functions at 2 . 7 < z < 5 . 3': "Alberto Torralba-Torregrosa 1 , 2 ⋆ , Pablo Renard 3 , Daniele Spinoso 3 , Pablo Arnalte-Mur 1 , 2 , Siddhartha Gurung-López 1 , 2 , Alberto Fernández-Soto 4 , 5 , Enrique Gaztañaga 6 , 7 , 8 , David Navarro-Gironés 7 , 8 , Zheng Cai 3 , Jorge Carretero 9 , 10 , Francisco J. Castander 7 , 8 , Martin Eriksen 11 , 10 , Juan Garcia-Bellido 12 , Hendrik Hildebrandt 13 , Henk Hoekstra 14 , Ramon Miquel 11 , 15 , Eusebio Sanchez 9 , Pau Tallada-Crespí 9 , 10 , Juan De Vicente 9 , and Enrique Fernandez 11  \n- 1 Observatori Astronòmic de la Universitat de València, Ed. Instituts d'Investigació, Parc Científic. C / Catedrático José Beltrán, n2, 46980 Paterna, Valencia, Spain\n- 2 Departament d'Astronomia i Astrofísica, Universitat de València, 46100 Burjassot, Spain\n- 3 Department of Astronomy, Tsinghua University, Beijing 100084, China\n- 4 Instituto de Física de Cantabria (CSIC-UC), Avda. Los Castros s / n, 39005 Santander, Spain\n- 5 Unidad Asociada 'Grupo de Astrofísica Extragaláctica y Cosmología', IFCA-CSIC / Universitat de València, València, Spain\n- 6 Institute of Cosmology and Gravitation (University of Portsmouth), Portsmouth, UK\n- 7 Institute of Space Sciences (ICE, CSIC), E-08193 Barcelona, Spain\n- 8 Institut d'Estudis Espacials de Catalunya (IEEC), Edifici RDIT, Campus UPC, 08860 Castelldefels (Barcelona), Spain\n- 9 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Avenida Complutense 40, E-28040 (Madrid), Spain\n- 10 Port d'Informació Científica (PIC), Campus UAB, C. Albareda s / n, 08193 Bellaterra (Barcelona), Spain\n- 11 Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, 08193 Bellaterra (Barcelona), Spain\n- 12 Instituto de Fisica Teorica (UAM / CSIC), Nicolas Cabrera 13, E-28049 Madrid, Spain\n- 13 Ruhr University Bochum, Faculty of Physics and Astronomy, Astronomical Institute (AIRUB), German Centre for Cosmological Lensing, 44780 Bochum, Germany\n- 14 Leiden Observatory, Leiden University, Einsteinweg 55, 2333 CC, Leiden, The Netherlands\n- 15 Institució Catalana de Recerca i Estudis Avançats (ICREA), 0810 Barcelona, Spain  \nAccepted XXX. Received YYY; in original form ZZZ", 'ABSTRACT': "We present the Lymanα (Ly α ) and ultraviolet (UV) luminosity function (LF), in bins of redshift, of quasars selected in the Physics of the Accelerating Universe Survey (PAUS). A sample of 915 objects was selected at 2 . 7 < z < 5 . 3 within an e ff ective area of ∼ 36 deg 2 observed in 40 narrow-band filters (NB; FWHM ∼ 120 Å). We cover the intermediate-bright luminosity regime of the LF (10 43 . 5 < ( L Ly α/ erg s -1 ) < 10 45 . 5 ; -29 < M UV < -24). The continuous wavelength coverage of the PAUS NB set allows a very e ffi cient target identification and precise redshift measurements. We show that our method is able to retrieve a fairly complete ( C ∼ 85%) and pure ( P ∼ 90%) sample of Ly α emitting quasars for L Ly α > 10 44 erg s -1 . In order to obtain corrections for the LF estimation, and assess the accuracy of our selection method, we produced mock catalogs of 0 < z < 4 . 3 quasars and galaxies that mimic our target population and their main contaminants. Our results show a clear evolution of the Ly α and UV LFs, with a declining tendency in the number density of quasars towards increasing redshifts. In addition, the faint-end power-law slope of the Ly α LF becomes steeper with redshift, suggesting that the number density of Ly α -bright quasars declines faster than that of fainter emitters. By integrating the Ly α LF we find that the total Ly α emitted by bright quasars per unit volume rapidly declines with increasing redshift, being sub-dominant to that of star-forming galaxies by several orders of magnitude by z ∼ 4. Finally, we stack the NB pseudo-spectra of a visually selected 'golden sample' of 591 quasars to obtain photometric composite SEDs in bins of redshift, enabling to measure the mean IGM absorption by the Lymanα forest as a function of redshift, yielding results consistent with previous spectroscopic determinations.  \nKey words. Methods: observational - Quasars: emission lines - Galaxies: luminosity function - Galaxies: high-redshift - Line: identification", '1. Introduction': "Active galactic nuclei (AGN) are extremely energetic events powered by e ffi cient accretion into super massive black holes (SMBH) in the cores of galaxies. These objects are among the most luminous in the Universe across the full electromagnetic spectrum (e.g., Shen et al. 2020), and can be easily detected at  \nvery high redshift from ground-based observatories ( z > 6; e.g., White et al. 2003; Jiang et al. 2016; Wang et al. 2021). Recent studies confirmed that most massive galaxies host a SMBH in their center (e.g., Askar et al. 2022; Tremmel et al. 2024), and multiple lines of evidence have been found that these SMBH coevolve with the bulge of the host galaxy (e.g., Yang et al. 2019; Mountrichas 2023). Additionally, most galaxy formation models require feedback from AGNs to prevent the interstellar gas to  \novercool, quenching the star formation (e.g., Dubois et al. 2016; Davé et al. 2019; Ward et al. 2022) - otherwise, the models predict an excessively large number of luminous galaxies, in contradiction with the observations (Fabian 2012). Hence, to investigate the properties of the AGN population, and how it evolves with cosmic time is a key piece to fully understand the big picture of galaxy evolution.  \nThe luminosity function (LF) of AGNs, defined as the number density of these objects per comoving volume as a function of luminosity, evolves with cosmic time (e.g., Kulkarni et al. 2019; Shen et al. 2020), both in shape and normalization: thefaint end of the bolometric AGN LF peaks at z ≈ 2-2 . 5 - coinciding with the peak of star formation history (e.g., Madau & Dickinson 2014) whilst the bright-end peaks at slightly earlier cosmic time z ≈ 3, and evolves slower toward the lowest redshifts, in what is often called 'cosmic downsizing' (e.g., Hasinger et al. 2005; Hopkins et al. 2007; Croom et al. 2009; Fanidakis et al. 2012).  \nThe Lymanα (Ly α ) emission-line is typically the most luminous in the rest-frame ultraviolet (UV) spectra of quasars (QSO; e.g., Vanden Berk et al. 2001; Ivashchenko et al. 2014; Harris et al. 2016), and star forming galaxies (SFG; e.g., Partridge & Peebles 1967; Pritchet 1994; Nakajima et al. 2018). The observed wavelength of Ly α is redshifted into the optical range at z ∼ 2-6, allowing to detect this line in ground-based surveys. Sometimes, secure identification of Ly α is possible even without an explicit detection of the spectral continuum (e.g., Bacon et al. 2015). On the one hand, the faint-end of the Ly α LF is dominated by a population of typically small, young, metal and dust-poor SFGs ( L Ly α ≲ 10 43 . 5 erg s -1 ; e.g., Guaita et al. 2010; Drake et al. 2017; Arrabal Haro et al. 2020; Santos et al. 2020). The Ly α LF deviates from a typical Schechter shape for L Ly α ≳ 10 43 . 5 erg s -1 , and the detection of X-ray counterparts of the sources at this range of Ly α luminosities suggests that most of them are powered, at least partially, by AGNs (Sobral et al. 2018; Matthee et al. 2017). On the other hand, the brightest end of the Ly α LF ( L Ly α > 10 43 . 5 erg s -1 ) is dominated by QSOs (e.g., Calhau et al. 2020; Spinoso et al. 2020). The Ly α LF in the whole range is well described by a double Schechter function (Schechter 1976). The intermediate Ly α luminosity regime ( L Ly α ∼ 10 43 . 5 erg s -1 ) is particularly interesting, because it might contain objects where the Ly α emission from the AGN and star-formation activity is similarly relevant (e.g., Sobral et al. 2018). While the Ly α LF of SFGs has been widely explored in the literature, (e.g., Cowie & Hu 1998; Hu et al. 1998; Martin & Sawicki 2004; Malhotra & Rhoads 2004; Dawson et al. 2007; Ouchi et al. 2008; Blanc et al. 2011; Adams et al. 2011; Bacon et al. 2015; Karman et al. 2015; Drake et al. 2017; Santos et al. 2016; Konno et al. 2018; de La Vieuville et al. 2019; Hu et al. 2019; Khusanova et al. 2020; Ono et al. 2021; Santos et al. 2021; Morales et al. 2021; Wold et al. 2022; Xu et al. 2023; Thai et al. 2023), the AGN-dominated, bright-end of the Ly α LF has only been constrained in the recent years (Spinoso et al. 2020; Zhang et al. 2021; Liu et al. 2022; Torralba-Torregrosa et al. 2023).  \nIn the same manner as the Ly α LF, the UV LF is dominated by SFG at the faint-end ( M UV ≳ -23) and AGN in the brightend ( M UV ≲ -23; e.g., Adams et al. 2023). Contrary to the Ly α LF, the bright-end of the UV AGN LF has been broadly explored (e.g., Glikman et al. 2011; Ross et al. 2013; McGreer et al. 2013; Palanque-Delabrouille et al. 2016; Jiang et al. 2016; Yang et al. 2016; Kulkarni et al. 2019; Akiyama et al. 2018; Matsuoka et al. 2018; Schindler et al. 2019; Niida et al. 2020; Pan et al. 2022). However, the faintest end of the AGN UV LF is still uncertain. In particular, a population of faint dust-obscured AGNs has been  \nunveiled by the most recent JWST observations (e.g., Labbé et al. 2023; Matthee et al. 2024; Kokorev et al. 2024; Kocevski et al. 2024; Greene et al. 2024; Akins et al. 2024), whose number density appears to be larger by more than one order of magnitude than the expected when extrapolating the bright-end of the AGN UV LF. The connection between this population of faint dusty AGNs and bright QSOs is still unclear, and the determination of the AGN LF in the intermediate luminosity regime ( M UV between -24 and -20) is crucial in order to solve this problem. The global determination of the AGN LF is key to determine the role of this population in the cosmic reionization. Hydrodynamical and radiative transfer simulations suggest that the contribution of AGNs to hydrogen reionization is sub-dominant (e.g., Dayal et al. 2020; Trebitsch et al. 2021). This is supported by the rapid decline of the number of bright quasars at high redshifts (e.g., Kulkarni et al. 2019). However, the role of AGNs can be determinant in helium reionization, which happened at later cosmic time (e.g., Compostella et al. 2014; Worseck et al. 2016).  \nIn this context, in this work we seek to identify and characterize the population of quasars with bright ( L Ly α > 10 43 . 5 erg s -1 ) Ly α emission using the narrow-band (NB) data from the Physics of the Accelerating Universe Survey (PAUS). Multi-NB surveys such as PAUS are extremely suitable for this task. On the one hand, a continuous coverage of great part of the optical range with NB photometry yields a low-resolution spectrum for every source in the field-of-view without any target pre-selection, allowing for very complete samples, if carefully selected. On the other hand, it becomes possible to probe a wide and continuous range of redshift, unlike surveys that only use a few NBs placed in strategic wavelengths. In addition, the availability of pseudo-spectra can be used to obtain precise photometric redshifts (Martí et al. 2014; Eriksen et al. 2020; Alarcon et al. 2021; Hernán-Caballero et al. 2021; Laur et al. 2022; Navarro-Gironés et al. 2023; Hernán-Caballero et al. 2024).  \nIn Torralba-Torregrosa et al. (2023) (hereafter TT23), we presented a method to select Ly α emitters in the miniJPAS (Bonoli et al. 2021) and J-NEP (Hernán-Caballero et al. 2023) multi-NB data. In this work we adapt and improve the selection method to the PAUS data, and use it to estimate the UV and Ly α LFs of QSOs selected through Ly α emission. We also use the obtained sample of Ly α -selected QSOs to estimate the mean absorption by the intergalactic medium (IGM) due to the Ly α forest as a function of redshift (e.g., Faucher-Giguère et al. 2008; Becker et al. 2013; Inoue et al. 2014).  \nThis paper is structured as follows. In Sect. 2 we describe the observational data of the PAU Survey, employed throughout this work. In Sect. 3, we first describe the mock catalogs, used to test and validate our selection method; we define our Ly α selection method, and the procedures to estimate the photometric redshifts of the sources, the Ly α luminosity and to obtain the Ly α and UV LFs. In Sect. 4 we present the Ly α -selected sample, and produce photometric composite QSO spectra in bins of redshift. We present the Ly α and UV LFs, and their integrated values as a function of redshift. In Sect. 5 we compare our measured Ly α and UV LFs with previous determinations in the literature, and discuss the evolution of these with redshift. We also comment on the e ffi ciency of our selection method, comparing with other AGN samples in the literature, and cross-matching our catalog with spectroscopic surveys. Finally, Sect. 6 summarizes the contents of this work.  \nThroughout this work we use a Λ CDM cosmology as described by Planck Collaboration et al. (2020), with ΩΛ = 0 . 69, Ω M = 0 . 31, H 0 = 67 . 7 kms -1 Mpc -1 ; unless specified other-  \nAll the magnitudes are given in the AB system (Oke & Gunn 1983).", '2. Observations: The PAU Survey': "The Physics of the Accelerating Universe Survey ( PAUS ) is a wide-field multi-narrowband survey carried out by the PAUCam instrument (Padilla et al. 2019), mounted at the 4.2 m William Herschel Telescope (WHT) at the Roque de los Muchachos Observatory in the Canary Islands. PAUCam is a 18-detector camera, with a FoV of diameter ∼ 1 · , designed to cover a wide area of the sky making use of a filter set composed of 40 NBs (FWHM ∼ 115 Å). The unique set of NB filters provides continuous wavelength coverage over observed wavelengths 45008500 Å, allowing to obtain a low resolution pseudo-spectrum for every source in the FoV. The detection of di ff use Ly α emission in PAUS has already been evaluated in Renard et al. (2021), and other prominent spectral features have already been studied with PAUS data, such as emission line ratios (Alarcon et al. 2021), or the D4000 spectral break (Renard et al. 2022).  \nPAUS has targeted the CFHTLS (Cuillandre et al. 2012) W1, W3andW4fields; the GAMA G09 field (Driver et al. 2022); and the COSMOS field (Scoville et al. 2007). In this work we use the data of W1, W3 and G09 (hereafter dubbed W2). In total, these three fields are covered with the full set of 40 filters in an e ff ective area of 35.68 deg 2 (after masking). For PAUS catalogs, all objects with magnitudes brighter than i = 23 are selected; at fainter magnitudes, the photometry S / N is too low to reliably measure galaxy redshifts (Alarcon et al. 2021; Navarro-Gironés et al. 2023).  \nThe standard PAUS photometry is computed using forced photometry (Serrano et al. 2023). For each object and narrowband, its fluxes are separately measured on each single-exposure image. Three exposures are taken per pointing, though the number of exposures on individual objects may vary if they are positioned near the edges of the tiling, due to overlaps and dithering. The flux apertures are ellipses that contain 62.5% of the flux, according to their Sérsic profiles convolved with the seeing of each exposure; the total flux for a given object and narrow band is the inverse-variance weighted average of all exposures. This method results in a custom aperture for each object, fitted to their size and morphology, which in turn yields a higher S / Nthan standard fixed apertures (Serrano et al. 2023). For this work, we computed custom photometric runs of W1, W2 and W3 utilizing the same data reduction outlined in Serrano et al. (2023), but using fixed circular apertures of diameter 3 '' instead. This fixed aperture was chosen in order to be consistent with the photometry used in TT23.  \nThe positions of all objects, their i magnitudes and morphological parameters are extracted from broad-band reference catalogs; no detections are performed on PAUS narrow-band images. The reference catalog used for the W1 and W3 fields is the CFHTLS data release from Heymans et al. (2012), and for the W2 (GAMA G09) field it is the fourth KiDS data release (Kuijken et al. 2019).", '3. Methods': 'In this Section we describe the methodology employed to select high redshift QSOs in the PAUS wide fields, and the estimation of the Ly α and UV luminosity functions. We first describe the mock catalogs developed to mimic our observations, which we use to test and assess the e ff ectiveness of our methods (Sect. 3.1). Next, we describe our target selection method,  \nbased on Ly α emission detection (Sect. 3.2), we characterize the purity and completeness of our method (Sect. 3.3) and describe the procedure to estimate the Ly α and UV luminosity functions (Sect. 3.4).', '3.1. Mock catalogs': 'We generate mock catalogs and photometric NB pseudo-spectra of Ly α emitting QSOs (our target sample) and their common contaminants (lowz QSOs and galaxies), with the aim of characterizing the e ff ectiveness of our selection method. The mocks are employed to train machine learning algorithms to improve the selection function (Sect. 3.2) and to measure accurate photometric redshifts (Sect. 3.2.6); to estimate the purity and completeness of the selected sample (Sect. 3.3) and to obtain corrections for the estimation of the Ly α and UV luminosity functions (Sect. 3.4).', '3.1.1. QSO mock catalog': 'The mock catalog of the QSO population at 0 < z < 4 . 5 is generated following the same procedure as in TT23, which at the same time is based on the method presented in Queiroz et al. (2023) for a QSO mock population in the J-PAS survey. We select spectra in this redshift range from the SDSS DR16 SupersetQ (Lyke et al. 2020) with good median signalto-noise over all pixels ( SN\\_MEDIAN\\_ALL > 0), measured redshift with high level of confidence ( ZWARNING = 0); and classified as QSO ( IS\\_QSO\\_FINAL = 1). Next, we sample values over a grid of redshift and r -band magnitude from the luminosity functions of Palanque-Delabrouille et al. (2016) (PLE + LEDE model). We look for the closest objects in redshift within the SDSS DR16 SupersetQ sub-sample, and apply the corresponding multiplicative factors to the observed wavelength array and fluxes to match the sampled values of ( z , r ). In TT23 we show that this procedure yields a population of QSOs that, for z > 2 is in good agreement with the Ly α LF measured in the literature (see Fig. 2 in TT23). We produce a QSO mock with a size equivalent to 1000 deg 2 , and a × 5 larger area for objects with L Ly α > 10 44 erg s -1 to overcome the limitations due to the low number of objects at such high luminosities. The sizes of the simulated areas are chosen as arbitrarily large values, enough to obtain smooth distributions in the mock properties (for more details about the QSO mock see TT23). The number of objects in these mocks as a function of redshift is shown in Fig. 2.  \nAs a final step, we compute synthetic photometry of all the selected objects using the transmission curves of all the PAUS filters to obtain pseudo-spectra. Since all the SDSS spectra used to generate the mocks are selected to have good S / N, the error of our synthetic photometry is negligible compared to the error of real PAUS observations. We model the PAUS photometric errors separately. For this, we model the uncertainties as a function of photometric magnitude with an exponential (see TT23), fitting the exponential parameters for every PAUS filter in each one of the W1, W2 and W3 fields. We add random Gaussian perturbations to the synthetic fluxes of the mock according to this model to obtain the final mock QSO catalog.', '3.1.2. Low-redshift galaxy mock catalog': 'Among the main contaminants in the search for Ly α emission there are nebular emission-lines associated to star-formation such as H α ( λ 6565 Å), H β ( λ 4861 Å), [O iii ] ( λλ 4959 , 5007 Å)  \nFig. 1. Transmission curves of PAUS NBs and BBs of the reference catalogs. On top of the transmission curves, horizontal lines are displayed showing the FWHM of each filter. We show the Ly α redshift corresponding to the pivot wavelength of each NB. As shown in Table A.1, only NB455-NB755 are used in this work, since we obtained no reliable detections in NBs of higher wavelengths.  \n<!-- image -->  \nFig. 2. Redshift distribution of galaxies within the QSO and lowz galaxy mocks. The di ff erent mock populations were simulated in different areas: 9, 1000 and 5000 deg 2 for the lowz galaxy, QSO and high L Ly α QSO mocks, respectively. In all mocks we applied conservative magnitude cuts of r < 25. Despite the much larger area used in the QSO mocks, the lowz galaxies are significantly more numerous.  \n<!-- image -->  \nor [O ii ] ( λλ 3727 , 3729 Å) emitted by lowz galaxies (i.e., z < 2, in this context). On top of their eventual bright line-emission, lowz galaxies may contaminate our LAE sample due to random fluctuations in their NB fluxes. This source of contamination may be relevant because lowz galaxies account for the bulk of the objects in the parent photometric catalog. Despite this, TT23 found little contamination of nebular emission lines form z < 2 galaxies, with the minor exception of a few extreme [O ii ] emitters.  \nIn order to assess the impact of lowz galaxies as contaminants in our samples, we make use of a galaxy mock, built by using the semi-analytic model (SAM) L-Galaxies on top of the dark-matter merger trees of the Millennium N-body dark matter (DM)-only simulation (Springel et al. 2005). In particular, the version of L-Galaxies we use is derived from the public  \nrelease of Henriques et al. (2015), with the inclusion of detailed models for the formation, mass-growth and dynamical evolution of SMBHs recently presented in Izquierdo-Villalba et al. (2020, 2022) and Spinoso et al. (2023).  \nOn top of being one of the most detailed SAMs in the literature to describe the evolution of galaxies and SMBHs, L-Galaxies includes the possibility to produce photometric mock-data on-the-fly. Indeed, Izquierdo-Villalba et al. (2019) presented a methodology to simulate an observed light-cone for a generic multi-band photometric survey, showing the example of the J-PLUS survey (Cenarro et al. 2019). In particular, their method relies on simulating the continuum emission of galactic stellar populations by means of a galaxy-templates grid which depends on metallicity, stellar age and a given choice for the stellar initial mass function (Izquierdo-Villalba et al. 2019 opted for a standard Chabrier IMF; Chabrier 2003). On top of the stellar continua, nine di ff erent nebular emission lines 1 from starforming regions are simulated by using the method presented in Orsi et al. (2014), where the emission line flux is obtained as a function of stellar age, gas density, metallicity and ionization parameter (see Izquierdo-Villalba et al. 2019, for further details on the implementation of the emission-lines model in L-Galaxies ).  \nIn order to build our lowz galaxy mock we follow the methodology outlined above, adapting it to the PAUS photometric filter-set and survey volume. In particular, we simulate a 3 × 3 deg = 9 deg 2 sky patch; that is, roughly one-fourth of the e ff ective PAUS area we use (see Sect. 3.2.1). The choice of this area is mainly limited by computation time for the lightcone. Despite the much smaller area as compared to the QSO mock, the number of galaxies in the lightcone is similar to the number of objects in the QSO mock (see Fig. 2). We apply a conservative magnitude cut (i.e., i < 25) to the photometry of galaxies within the light-cone. Given the e ff ective depth of PAUS data, this magnitude cut allows us to significantly reduce the size of our lowz galaxy mock data without compromising their com-  \nleteness. The PAUS photometric errors are simulated as described in Sect. 3.1.1 for the QSO mock. Our final mock catalog contains ∼ 5 × 10 5 galaxies with i < 25 at redshifts 0 < z ≲ 6 (median z = 1 . 2) (see Fig. 2).', '3.2. Lymanα selection method': 'In this Section we describe the selection method used to obtain a sample Ly α selected sample of QSOs in the footprint of PAUS.', '3.2.1. Parent sample': "We use the diameter 3 '' aperture photometry catalogs of the PAUS wide fields W1, W2 and W3 as discussed in Sect. 2. For our parent sample, we use only the regions of these fields where there is coverage for all 40 PAUS NBs. The reference catalogs used by these fields are taken from surveys with deeper broad-band imaging and better PSF than the average of PAUS NBs. Due to this, some sources from the reference catalogs are blended in the PAUS NB images, leading to contamination in the fixed aperture photometry. Hence, we remove objects with neighbours with a separation closer than 3 '' . This cut removes the 8% of the objects across the W1, W2 and W3 catalogs, hence this will be taken into consideration in the computation of the LF.  \nAfter removing potentially blended sources, and applying the PAUS field masks for the area covered by all 40 NB filters, we are left with a parent sample containing 1 022 196 objects (239 312, 267 393 and 515 491 in W1, W2 and W3, respectively) in a total area of 35.7 deg 2 (e ff ective area after masking).", '3.2.2. NB flux excess selection': 'We search for Ly α emission in each of the PAUS narrow bands through a NB flux excess selection. Each narrow band probes Ly α lines at redshifts within δ z = ± 0 . 06 around the central redshifts listed in Table A.1, and displayed in Fig. 1. This redshift interval is defined by the FWHM of the NB filters. A detailed description of the LAE candidate selection method can be found in TT23. In the following Sections we summarize the method, which was adapted for the particularities of PAUS data.  \nIn the first place, we look for tentative Ly α detections in all the NBs. For each NB, an estimation of the continuum flux is obtained by averaging the photometric fluxes of the 5 closest NBs to the one of interest, at bluer and redder wavelengths, ignoring both directly adjacent NBs to avoid contamination from the tail of broad Ly α lines (see e.g., Greig et al. 2017). That makes a total of 10 NB filters to compute the continuum level. In addition, the flux of the NBs blueward from Ly α is a ff ected by the Ly α forest. The mean absorption of the Ly α forest is corrected in the continuum estimation using the curve for the mean IGM transmission in Faucher-Giguère et al. (2008). After the estimation of the continuum flux level of every source in our parent sample, we look for NBs with S / N > 6 and a 3 σ flux excess over the continuum:  \nf λ NB -f λ cont > 3 q σ 2 NB + σ 2 cont . (1)  \nWe also impose that this excess is compatible with a Ly α line with EW0 ≥ EW min 0 = 20 Å by applying  \nf λ NB f λ cont > 1 + (1 + z NB) · EW min 0 FWHMNB , (2)  \nwhere z NB is the Ly α redshift corresponding to the target NB. (see TT23, also e.g., Santos et al. 2016; Sobral et al. 2018;  \nFig. 3. Four examples of candidates selected by our pipeline. The colored squares represent the NB fluxes of PAUS. We compare with the SDSS spectrum in gray (resampled by a factor 10 to reduce visual noise). We mark the most important QSO emission lines (dotted gray vertical lines) and the predicted observed wavelength of Ly α by our selection pipeline (vertical red line). The estimation of the continuum level is marked with a black square. The first three panels are examples of correctly selected Ly α emitters. The last panel is an example of a confusion between Ly α and C iv ( z C iv = 2 . 1), the most common source of contamination.  \n<!-- image -->  \nSpinoso et al. 2020). If two or more contiguous NBs meet these criteria, the Ly α detection is assigned to that of higher flux.  \nAfter the preliminary Ly α detection, we similarly search for other lines in the selected sources. For this secondary line search, following the optimal search criteria in TT23, we look for NBs with > 3 σ flux over continuum excess and EWobs > 100 Å. These emission lines are searched using a continuum estimate without the Ly α forest transmission correction. The algorithm checks the secondary line detections and discards candidates with lines not compatible with the expected position of the most prominent QSO lines other than Ly α ; in particular Ly β + O vi , C iv , C iii ] and Mg ii . Figure 3 shows four hand-picked examples of selected candidates using the procedure described above. The first three panels show correctly selected QSOs with coinciding spectroscopic Ly α redshifts of their SDSS counterparts. In these three examples, we also detect the secondary C iv line, which en-  \nFig. 4. Classification of the candidates with spectroscopic counterpart by our ANN algorithm. The solid light teal histogram represents the full spectroscopic sample, and the solid dark blue the sub-sample classified as LAE by the NN. The empty brown and green histograms represent the candidates classified as contaminant QSOs and lowz galaxies, respectively. For our target redshift range ( z ≳ 2 . 75), the algorithm correctly classifies most objects ( ∼ 90% completeness). The red dotted vertical line marks the minimum redshift targeted by the NB with the lowest wavelength ( z = 2 . 7, NB455).  \n<!-- image -->  \ncorrect determination of the redshift. The bottom panel shows an example of the most common contaminant: a confusion between Ly α and a lower redshift ( z = 2 . 3) C iv emission line.', '3.2.3. Artificial neural network classifier': "In TT23 a visual inspection was performed through all the candidates in order to remove obvious contaminants such as lowz galaxies with extended morphology, or NBs a ff ected by cosmic rays and other artifacts. In this work, we trained a multi-layer perceptron artificial neural network (ANN) model in order to identify secure candidates and discard contaminants. This is motivated by the notable increase in the size of the candidate sample with respect to that in TT23, and the potential application of the method to even larger datasets in the future.  \nIn order to train the ANN, we generate training and test sets from the subsets of objects in our mocks selected by NB flux excess (see Sect. 3.2.2). Each set contains equal number of the three classes to consider: (1) QSO with the correct narrow-band redshift, (2) QSO with wrong narrow-band redshift and (3) lowz galaxy. The first of these three classes corresponds to the true positive LAE detections, while the other two represent the main contaminants.  \nThe input for the ANN classifier is composed of 46 features:  \n- -All 40 NB fluxes. Normalized to the maximum NB flux of each object.\n- -g , r , i and z BB fluxes. Normalized to the maximum BB flux. The band u is omitted for being partially out of the SDSS range, used to create the mocks.\n- -r synth. Synthetic broad-band magnitude computed from the NB whose wavelength coverage reproduces the classical r band ( λ ≈ 5750-7000 Å). Normalized with a factor 1 / 30 (so every value is < 1).\n- -NBLy α : the NB of the Ly α detection, numbered from 0 to 39. Normalized with a factor 1 / 39.  \nThe NB and BB arrays of fluxes are normalized separately since the BB data is obtained from reference catalogs external to PAUS (see Sect. 2), and their fluxes might not be directly comparable to PAUS 3 '' aperture NB fluxes. The r synth feature gives information about the apparent magnitude of the object, and NBLy α establishes an initial guess of the redshift of the source. A synthetic composite r synth broad-band is used instead of r from the reference catalogs for consistency.  \nWeperform a search for optimal parameters of the ANN over a wide parameter space. The optimal parameters obtained in the grid-search, and more details about the ANN architecture can be found in Appendix B. As shown in Fig. B.1, the classifier yields a ∼ 90% pure sample of QSO LAEs over the test set ( ∼ 95% for L Ly α > 10 44 erg s -1 ), with 7% (2%) contamination from QSOs (galaxies) at lower redshifts. The completeness of the sample classified as QSO LAE is 96% and 98% in the test set, for the general sample and for L Ly α > 10 44 erg s -1 , respectively.  \nWe also determine the accuracy of the ANN classifier with our spectroscopic sub-sample (presented in Sect. 5.4 below). In Fig. 4 we show the classification over this sub-sample. For z ≳ 2 . 75, our target redshift range (see Table A.1), the classifier correctly identifies as LAEs ∼ 90% of the objects, in line with the completeness estimated in the test set. For lower redshift, z ≲ 2 . 75 (i.e., the range of the contaminants), the algorithm successfully identifies most of the objects as contaminants. The confusion between contaminant QSO and galaxies is irrelevant in this work, since we are only interested in the QSO LAE subsample. The purity derived from this result is ∼ 86% ( > 90% for L Ly α > 10 44 erg s -1 ), again in line with that of the test set.", '3.2.4. Visual inspection': "The mock data used to train the ML model described in Sect. 3.2.3 only contains objects with clean synthetic photometry, performed over either simulated SEDs (in the case of the lowz galaxy mock) or spectra with no warnings raised by the SDSSpipelines (in the case of the QSO mocks). Hence, the training set does not account for objects a ff ected by instrumental artifacts, cosmic rays or photometry calibration problems, among other issues that can a ff ect the NB fluxes leading to a false detection. For this reason, despite the good performance shown by the MLmodel, a visual inspection is still necessary in order to catch objects misclassified by the algorithm due to unexpected features. Nonetheless, the ML classification significantly reduced the candidate sample for the visual analysis, from an initial sample of 4356 objects to 999 candidates.  \nThe NB photometry of some objects is clearly a ff ected by the extended emission of a nearby object, leading to unusually high NBfluxes that are mistaken by Ly α emission lines. In some other objects, the images were a ff ected by cosmic rays, instrumental artifacts or similar e ff ects that alter the NB fluxes. Finally, some images of the PAUS fields are a ff ected by calibration issues, resulting in unusually large zero-points. The zero-points are computed with a cross-match of bright SDSS stars on a per-CCD basis (Serrano et al. 2023, each CCD zero-point being the median of the reference stars). Hence, flux underestimation of reference stars (e.g., due to scatter-light or cross-talk between CCDs) could result on abnormally high zero-points for a few exposures. This issue artificially enhances the measured flux of some NBs in some regions, leading to wrong Ly α line detections. All these problems cannot be trivially reproduced in our mocks, thus the a ff ected objects must be removed from the selected sample via visual classification.  \nWe visually inspected the initial sample of 999 QSO LAE candidates selected by our method, and agreed on the final classification. The inspection yielded a 'golden sample' of 591 (59%) objects selected as true positive QSO LAEs at the correct redshift. On the other hand, a non-negligible amount of 102 (10%) objects were removed from the final sample due to issues with their pseudo-spectra or images. Some of these removed objects showed contamination by cosmic rays in their selected NB, evident instrumental artifacts, extended emission from nearby extended objects or clear issues with the photometric calibration that led to unusual pseudo-spectra. The remaining 316 (31%) objects are are left as 'uncertain', not being classified as true positives nor as contaminants.", '3.2.5. Highz sample': 'In addition to the search of candidates in the interval 2 . 7 < z < 4 . 2, we look for objects at higher redshifts. However, since our machine learning algorithm was trained using SDSS data, it is limited to the above mentioned redshift range, given that the SDSS sample becomes too scarce for higher redshifts. We perform a candidate search for z > 4 . 2, without applying the ANN classification, and inspect all selected objects to identify secure Ly α emitting QSOs. This search produced a sample of 22 visually selected QSOs ( ∼ 0 . 6 deg -1 ) at 4 . 2 ≲ z ≲ 5 . 3 (see Table A.1), from an initial sample of 549 highz candidates.', '3.2.6. NB redshift correction': 'Each NB can probe an observed Ly α wavelength in an interval defined by the width of the transmission curve. The mean redshift probed by the PAUS NBs used in this work are summarized in Table A.1. In addition, the full 40-NBs pseudo-spectra provide information that can be used to better constrain the redshifts of the detected sources. This is especially relevant for the correct assessment of the redshift of sources selected via Ly α detection. The Ly α emission-lines of these sources typically present broad profiles (width ∼ 4500 km / s; Greig et al. 2016), and the absorption due to the Ly α forest (for λ obs < λ Ly α obs ) typically induces a bias in the Ly α redshift towards slightly higher values.  \nWe improve the redshift accuracy of our selected objects using machine learning algorithms. In order to obtain train / test sets, we run our selection method on our mocks. Then we keep only those candidates whose di ff erence between the NB redshift ( z NB; i.e., the mean redshift of the selected NB) and the true redshift of the mock object ( z spec) is smaller than the width of a NB in redshift space ( ∆ z ≲ 0 . 12; see Table A.1). We train ANN and random forest (RF) regressors with the aim of retrieving the spectroscopic redshift for the objects in our test set. The input for the ML redshift regressors is the same as for the ANN classifier described in Sect. 3.2.3. We perform a search of optimal parameters for the ANN and RF algorithms over an wide grid. We find that the RF regressor greatly outperforms the ANN and is less a ff ected by over-fitting for the optimal parameters of both.  \nBottom panel of Figure 5 shows the di ff erence between the redshift obtained by the RF and the spectroscopic redshift of all the selected objects in the PAUS catalogs with spectroscopic counterparts in SDSS, DESI and HETDEX. The initial NB redshifts are biased with respect to those of the spectroscopic counterparts, and present larger variance than the redshifts obtained by the RF: the bias and standard deviation are ( µ, σ )( ∆ z ) = (0 . 04 , 0 . 03) for the initial NB redshifts, and ( µ, σ )( ∆ z ) = (0 . 003 , 0 . 013) and (0 . 010 , 0 . 012) for the redshifts  \nz  \n<!-- image -->  \nFig. 5. Top: Comparison between photometric and spectroscopic redshifts of the candidates in our sample with counterparts in SDSS DR18 (green circles) HETDEX (orange triangles) and DESI DR1 (blue crosses). The straight lines mark the confusion between Ly α and other bright QSO lines. Thus, objects laying on the red line are true positive detections of LAEs. Bottom: Di ff erence between the photometric and spectroscopic redshifts of the objects shown in the top panel. The black dots represent the same objects using the photometric redshift of the selection NB, before applying the machine-learning algorithm presented in Sect. 3.2.6. The colored dots represent the same objects after the random forest redshift correction, for each spectroscopic subsample. The bottom right panel shows normalized histograms, result of collapsing the left panel in the horizontal axis. The HETDEX histogram is not shown due to poor statistics.  \n<!-- image -->  \nproduced by the RF, for the SDSS and DESI spectroscopic subsamples, respectively. The redshift bias of the DESI spectroscopic sub-sample is slightly larger than that of the SDSS sample. This can be explained by the fact that the mocks employed to train the RF algorithm are based on SDSS spectra. However, the RFestimation significantly improves the precision of the redshift measurements within the global spectroscopic sample.', '3.2.7. Ly α luminosity correction': 'The luminosity of the selected Ly α lines is estimated as  \nL Ly α ≈ GLYPH<16> f λ NB -f λ cont GLYPH<17> · FWHMNB · 4 π d 2 L ( z ) , (3)  \nwhere f λ NB and f cont are the NB and continuum fluxes; FWHMNB the full width at half maximum of the NB filter, and dL the luminosity distance according to our cosmology. Several e ff ects can alter the measurement of the Ly α flux. In the first place, the estimation of f cont is performed assuming a mean IGM transmission for the filters at shorter wavelengths than Ly α (see Sect. 3.2), yet there can be significant variability in the Ly α forest for individual  \nsources, especially at the lowest redshift bins (e.g., Rollinde et al. 2013). Emission lines from QSOs typically show broad profiles, so the tails of the Ly α profile may be partially out of the coverage of a given NB. This e ff ect is also enhanced if the intrinsic Ly α peak is not well aligned with the NB e ff ective wavelength. In addition, in Eq. 3 a top-hat filter profile is assumed 2 . Finally, while most spectroscopic surveys are capable of disentangling the N v QSO line ( λ 1240 . 81 Å) from Ly α , this task is not possible with NB photometry, since the Ly α and N v are cannot be separated within the single flux measurements obtained with PAUS NBs. The flux of N v is typically fainter than Ly α , but it can significantly a ff ect the Ly α measurement (see e.g., Vanden Berk et al. 2001; Selsing et al. 2016).  \nWe empirically correct for these biases based on the comparison of our measured Ly α fluxes with the spectroscopic determinations of the SDSS DR16 counterpart sub-sample. For this, we determine the systematic o ff set between our flux measurements and those of SDSS as a function of L Ly α , and subtract it from the flux measurements of our general sample (up to ∆ log 10 L Ly α ∼ 0 . 2; see Fig. 4 in TT23).', '3.3. Purity and completeness': "Following TT23, we compute 2D maps of purity and number count corrections in order to estimate the Ly α LF. We apply our selection method (Sect. 3.2) to the mocks, and distribute the selected mock sample over a 2D grid of measured Ly α luminosity and synthetic r band magnitude (i.e., L obs Ly α , r synth). We then compute the 2D purity as  \nP 2D = N sel ., true N sel . (4)  \nand the number count correction as  \nw 2D = N tot . N sel ., true (5)  \nwhere N sel . is the number of selected objects, Nsel ., true the number of correctly selected objects (i.e., Ly α emitting quasars at the correct redshift), and N tot . the total number of objects in the mock, which would all be selected if the method's completeness was 100%. This quantities are evaluated in each cell of the 2D grid. These estimates of P 2D and w 2D implicitly include the purity and completeness components of the NB flux excess selection and the ANN classifier. The visual inspection of our selected sample was performed only to remove artifacts unforeseen by the mocks, and thus we assume it does not a ff ect the completeness. We also neglect the incompleteness arising from the intrinsic limitations of the reference catalogs, as these catalogs are based on significantly deeper data and are complete up to i ∼ 24 . 5 (see Sect. 2).  \nBy marginalizing the 2D purity and completeness maps over the r synth axis we estimate the total purity and completeness of our sample in terms of L Ly α . The 1D purity and completeness estimations from our mocks is shown in Fig. 6. Figure C.1 is the 2D version of this figure in bins of r and L Ly α . The purity of the sample reaches ∼ 90% at log 10 ( L Ly α/ erg s -1 ) ≈ 44 for all the considered redshift bins. Conversely, the completeness reaches a maximum value of ∼ 85%. Interestingly, the maximum completeness values correspond to the higher redshift bins. Fig. 6  \nFig. 6. 1D purity and completeness of our sample as a function of measured L Ly α . We show this values in the bins of redshift used to estimate the Ly α and UV LFs (see Sect. 3.4). These curves were estimated by applying the selection method to the mocks, in bins of measured Ly α luminosity. The purity of the sample increases with L Ly α , reaching values of ∼ 1 for the high-end at all redshifts. The completeness, on the other hand, increases monotonically with L Ly α , reaching maximum values of ≈ 85% at L Ly α ≈ 10 44 erg s -1 at all redshifts.  \n<!-- image -->  \nalso shows that the purity increases more rapidly with luminosity at the highest redshift bins. The reason for this is that the 7 NBs with lowest e ff ective wavelengths are more likely to su ff er from higher contamination from QSO C iv lines; the rest of the NBs probe redshifts in which the Ly α line is observed in combination with C iv , allowing our method to identify the redshift correctly. This e ff ect reflects in Fig. 5 (top panel). This figure also illustrates that the purity of the sub-sample with spectroscopic counterparts is very close to unity for a large fraction of our z phot range, especially for z phot > 3 . 25, due to the decrease of line confusions.", '3.4. Ly α luminosity function computation': "In order to compute the Ly α LF we follow the same scheme detailed in TT23, adapting those methods to the specific case of this work. In order to obtain the Ly α LF, we need to characterize the purity and completeness of our LAE candidate sample, and apply the adequate corrections. Each candidate k is assigned a weight w = w 2D , being w 2D the number count correction from Eq. 5. On the other hand, we introduce the parameter pk whose value is zero if the k -th candidate is rejected, and unity if it is accepted to be included in the LF computation. Each candidate is accepted for the LF computation with a probability equal to the value of the purity P 2D k (Eq. 4).  \nThe Ly α LF for each considered NB set is computed an arbitrarily large number of times ( n = 1000) to achieve significant  \nstatistics. The i -th determination is computed as  \nΦ i [log 10 L Ly α ] = P k pk · wk V · ∆ GLYPH<16> log 10 L Ly α GLYPH<17> , (6)  \nwhere the index k runs through the LAE candidates within the redshift interval defined by the given NB set, in bins of Ly α luminosity; V is the maximum volume our survey can probe (Schmidt 1968), determined by the covered interval of wavelengths by the NBs and the survey area 3 . On each iteration, the Ly α luminosity and r flux is randomly perturbed using Gaussian uncertainties, and the parameters pk and wk are recomputed accordingly. In addition, a multiplicative factor (1 -0 . 08) -1 was applied to the LF in order to account the incompleteness due to the removal of sources with neighbours closer than 3 '' from the parent sample (see Sect. 3.2.1).  \nIn order to include the e ff ect of cosmic variance and fieldto-field variation e ff ects, the LF is computed for bootstrapped sub-regions of the wide PAUS fields. We define tiles of equal area ( ∼ 0 . 12 deg 2 ) for each of the three PAUS wide fields. These tiles are obtained by generating random points, from an uniform distribution in spherical coordinates, within the PAUS masks and defining regions with equal number of random points inside. Finally, the LF is obtained as the median of all the Φ i determinations obtained from the bootstrapped regions. The final errors on our LF are computed as the 16th and 84th percentiles of the same Φ i distribution. Finally, the UV LF can be computed in a similar way, using Eq. 6 substituting log 10 L Ly α with the absolute UV magnitude computed at rest-frame 1450 Å, M UV.", '4. Results': 'In this Section, we first present the composite QSO spectrum, obtained from the stack of NB fluxes of the objects in our sample (Sect. 4.1). Next, we present the Ly α and UV luminosity functions estimated from our data (Sects. 4.2 and 4.3).', '4.1. QSO composite spectrum': 'We stack the NB SEDs of our golden sample (see Sect. 3.2.4) in order to obtain a photometric composite QSO spectrum. For this task we use the code stonp 4 , developed to stack multi-band photospectra in a common rest-frame wavelength grid. Before stacking, each spectrum is normalized dividing by the r synth band (see Sect. 3.2.3). The composite spectrum is shown in Fig. 7, and it reproduces the expected features of a QSO spectrum: power-law continuum with a break caused by the Ly α forest ( λ 0 < 121 . 56 nm), and total absorption past the Lyman limit ( λ 0 < 91 . 2 nm); and the typical Ly β + O vi , Ly α + N v , Si iv + O iv , C iv and C iii broad emission lines.  \nWe fit a power-law to the QSO UV continuum, using the rest-frame wavelength windows 134 < ( λ 0 / nm) < 135 & 146 < ( λ 0 / nm) < 147 (Fig. 7). The choice of these windows is motivated by the goal of capturing the power-law shape of the spectrum in rest-frame wavelength intervals that do not contain emission lines nor any prominent non-continuum emission like the Fe ii emission complex (see, e.g., Vanden Berk et al. 2001).  \nRest-frame wavelength [nm]  \n<!-- image -->  \n×  \nFig. 7. Stacked composite spectrum of all QSOs in our golden sample. We show the mean composite spectrum (teal solid line) and bootstrap 1 σ interval (shaded region). The shaded grey regions mark the rest-frame wavelength intervals used to fit a power-law continuum. We represent the fitted power-law as a dashed red line, whose slope is β UV = -1 . 55 ± 0 . 02. We show the labels of the most important QSO emission lines and the LyC limit (91.2 nm).  \nTable 1. Intervals of NBs used for the Ly α and UV LF computation.Notes. The NB indices shown in the first two columns refer to the NB IDs defined in Table A.1.  \nThe spectral index (i.e., 1 + β UV, being β UV the continuum powerlaw slope) of QSOs is known to vary with the choice of the restframe wavelength window, as its measurement is susceptible to by a ff ected broad-emission features (e.g., thermal emission from the accretion disk, Fe ii emission series etc.; see e.g., Shull et al. 2012). For this fit, which will serve as reference, we use the stack of every QSO in the full redshift range in order to maximize the S / N. We obtain a best-fitting UV slope of β UV = -1 . 55 ± 0 . 02 for the global stack ( f λ ∝ λ β UV ). This result is in line with the slope measured by Vanden Berk et al. (2001) for a SDSS composite spectrum using a fitting window of 135 < ( λ 0 / nm) < 136 . 5 ( β UV = -1 . 56). More recently, Lusso et al. (2015) measured a UVslope of β UV = -1 . 61 ± 0 . 01, using a di ff erent choice for the wavelength windows (109.5-111.0, 113.5-115.0, 145.0-147.0, 197.5-201.0, and 215.0-220.0 nm).  \nFig. 8. Ly α luminosity function in 9 intervals of redshift defined in Table 1. We compare with Ly α LFs in the literature at similar redshifts. The green solid line and blue shaded area show the best Schechter fit with free parameters and the 1 σ confidence interval, respectively. The dashed blue line shows the Schechter best-fit with for fixed L ∗ = 45 erg s -1 . The Ly α LF of the highest redshift bin shown in this figure was computed assuming purity and completeness equal to 1, hence it should be regarded as a lower limit (See Sect. 4.2). The horizontal gray dotted line marks the 1 object per bin level. Data-points below this limit signal LF bins which contain less than 1 object after we average-out the 1000 LF determinations used to compute our final LF.  \n<!-- image -->', '4.2. The Ly α luminosity function at 2 . 7 ≲ z ≲ 5 . 3': 'Figure 8 shows the Ly α LF in the 9 overlapping intervals of redshift as defined in Table 1. We discuss the evolution of our LF in Sect. 5.1 below. The highest redshift bin in Fig. 8 is computed from the highz sample (see Sect. 3.2.5), and assuming purity and completeness equal to 1. The reason for this is that our 2D purity and number-count correction estimations cannot be accurately computed for z ≳ 4 . 5, where the SDSS DR16 sample begins to be scarce. As discussed in Appendix C of TT23, purity  \ninaccuracies do not significantly a ff ect the estimation of the LF, and completeness smaller than 1 induces a positive correcting factor in the LF. Hence, the results of this highestz bin should be taken as lower limits to the Ly α LF. We also note that, due to the limited number of candidates for such high redshifts, this last bin is considerably wider than the rest, covering z ≈ 4 . 26-5 . 26.', '4.2.1. Schechter fit to our data': 'We fit a Schechter function (Schechter 1976) to our data, of the form  \nΦ [ L ]d(log 10 L ) = log 10 · Φ ∗ · GLYPH<18> L L ∗ GLYPH<19> α + 1 · e -L / L ∗ d(log 10 L ) . (7)  \nWe perform the fitting using the Monte-Carlo Markov Chain (MCMC) technique, aided by the python package autoemcee 5 . Fit our Ly α LF in the full measured range for consistency. In all panels of Fig. 8, a deviation from the characteristic power-law shape of the Schechter function is evident at the faint end of the LFs. This e ff ect is likely due to inaccuracies in the completeness and purity models at such low luminosity intervals. However, our method e ff ectively accounts for these inaccuracies by generating large uncertainties, which reduce the influence of these bins on the fitted curve.  \nThe Schechter fit is performed using a normal loglikelihood and wide flat priors for all three free parameters: log 10 ( Φ ∗ / Mpc -3 ) ∈ [ -9 , -3], log 10 ( L ∗ / erg s -1 ) ∈ [43 , 46] and α ∈ [ -3 , -1]. The best-fit parameters are summarized in Table D.1. The best-fitting value for the normalization Φ ∗ marginally shows anti-correlation with redshift. The e ff ect of this evolution can be easily seen as the overall decline towards higher redshift of the Ly α LF in Fig. 8. The faint-end power-law slope α also shows steepening values with redshift, while the characteristic luminosity L ∗ has no clear evolutionary trend in the probed redshift range.  \nWe perform further Schechter fits to our Ly α LFs fixing log 10 L ∗ = 45 erg s -1 . Fixing this parameter removes the possible correlations with Φ ∗ ( Φ ∗ defines the value of the LF at L = L ∗ , see Eq. 4.2.1). The evolution of the parameter Φ ∗ becomes more evident, and the faint-end slope shows a very similar steepening trend with redshift. We discuss the fitted parameters more in detail in Sect. 5.1.', '4.2.2. Integrated Ly α luminosity function': 'In Fig. 10, we show the results of integrating the best-fit of the Ly α luminosity function in the range log 10 ( L Ly α/ erg s -1 ) > 43 . 5, for the redshift intervals shown in Fig. 8. The resulting total volume density of Ly α emission ρ Ly α shows a clear declining trend towards increasing redshifts. We fit this trend to log 10 ( ρ Ly α/ erg s -1 Mpc -3 ) = b + a · z , obtaining the parameters a = -0 . 53 ± 0 . 05 and b = 39 . 90 ± 0 . 16. For this fit, we excluded the highest redshift bin (empty circle in Fig. 10, however this lower-limit is in agreement with the extrapolation of the fitted evolution of ρ Ly α .  \nOur results are compatible with the measurements obtained by integrating the best-fits of the data presented in TT23, for most of the redshift bins. The surveyed area employed in TT23 was significantly smaller (by a factor ∼ 30), and the discrepancies can be due to cosmic variance. The extrapolation of our ρ Ly α best-fit is considerably lower than the result of integrating the double power-law fits in Liu et al. (2022), in their z = 2 . 25 , 2 . 37 and 2 . 54 bins, in the same luminosity range. These di ff erences are also reflected in the discrepancies with their faint-end of the Ly α LF that can be seen in Fig. 8. In their work, they used blind spectroscopy to select a sample of much fainter AGNs, which also include narrow-line AGNs and sources with significant contribution from star formation in the host galaxy, while our sample is mainly composed of luminous broad-line quasars.  \nWe also compare the presented measurements of ρ Ly α for the AGN Ly α LF with the measurements of Sobral et al. (2018), for the faint-end of the Ly α LF, populated by star-forming galaxies. The value of ρ Ly α for this population increases with redshift, showing the opposite trend to the bright Ly α emitting AGN population. Our results show that for z ≈ 4, the total emissivity of AGN Ly α is subdominant by a factor of ≈ 200, with respect to the star-forming LAEs. Moreover, we also display the estimated ρ Ly α from unresolved emitters and di ff use IGM emission from Lin et al. (2022); this value is expected to be significantly larger than than Ly α luminosity from resolved emitters (Byrohl & Nelson 2023), but can be explained by integration of faint LAE LFs (e.g., Drake et al. 2017; Sobral et al. 2018) to an arbitrarily faint end (Renard et al. 2024).', '4.3. UV luminosity function': 'One of the advantages of our multi-NB data is the possibility to measure SEDs for all the objects in the field without preselection, in an extensive range of wavelengths. We measure the UV absolute magnitude ( M UV) for every object in our selected sample, defined as the luminosity at rest-frame wavelength 145 nm (e.g., Glikman et al. 2010). The UV magnitude is thus computed as  \nM UV = m NB -DM( z ) -K , (8)  \nwhere m NB is the apparent magnitude of the NB corresponding to the rest-frame wavelength 145 nm of each object, DM( z ) the redshift dependent distance modulus according to our cosmology, and K the K -correction defined as  \nK = -2 . 5 log 10 (1 + z ) α + 1 , (9)  \nwhere α is the spectral index (i.e., α = -2 -β UV; see e.g., Wisotzki 2000). For the computation of the K -correction we adopt the UV slope measured in our NB composite spectrum ( β UV = -1 . 55, α = -0 . 45; see Sect. 4.1). We note that the continuous coverage of the PAUS filters allow us to photometrically measure M UV in the adequate window without having to apply bandpass corrections.  \nWe obtain the UV LFs for our sample following a very similar procedure to that used for obtaining the Ly α LF, described in Sect. 3.4. The median purity and number count corrections based on r synthL Ly α are applied, because the selection of the sources which make part of the UV LF is still based on the same selection function as the Ly α LF. Hence, our results presented in Fig. 11, correspond to the UV LFs of our Ly α -selected AGN sample. We compare our results with the best-fitting curves in the literature for similar redshift intervals (Kulkarni et al. 2019; Niida et al. 2020; Schindler et al. 2019; Pan et al. 2022).  \nThe last redshift bin of our UV LF (4 . 26 < z < 5 . 26) is computed in the same way as for the Ly α LF (see Sect. 4.2), using visually selected objects from a pre-selection by our method (see Sect. 3.2.5). As for the case of the Ly α LF (see introduction of Sect. 4.2), the purity and completeness of this redshift bin is not well determined, thus this result may be taken as a lower limit for the UV LF. Nevertheless, the UV LF in this last redshift bin is in good agreement with z ∼ 5 determinations in the literature (Kulkarni et al. 2019; Niida et al. 2020; Pan et al. 2022).', '4.3.1. Double power law fit': 'The faint-end of our UV LF ( M UV ≳ -24) is a ff ected by the poor purity and completeness (see Fig. 6). For M UV ≲ -24, the UV  \nFig. 9. Best-fit Schechter parameters for the LFs in the 9 intervals of redshift. The horizontal dashed lines represent the average best-fit Schechter parameters across all redshifts. When leaving all three parameters free (green squares), the normalization Φ ∗ and the characteristic luminosity L ∗ show a slight decline within the intervals of confidence, whilst the faint-end slope α shows no clear trend. We also show the fitted parameters when fixing α and L ∗ to the average values marked with the green dashed lines (blue triangles). The values of the latter are slightly shifted in the horizontal axis for visual clarity.  \n<!-- image -->  \nFig. 10. Ly α luminosity density per unit comoving volume. These points were obtained by integrating the Ly α LF for L Ly α > 10 43 . 5 erg s -1 . Our results show a declining trend with increasing redshifts, compatible with the results presented in TT23 and with the integrated Ly α LF from Liu et al. (2022). We also compare with the integrated LF of the SFG population in Sobral et al. (2018), which covers L Ly α ≲ 10 43 . 5 erg s -1 , showing that the integrated AGN Ly α luminosity is increasingly smaller than the galaxy contribution at higher redshift. We also compare with the di ff use Ly α emission measurement of Lin et al. (2022).  \n<!-- image -->  \nLF is typically described by a double-power law (DPL):  \nΦ [ M ] = Φ ∗ 10 -0 . 4( M ∗ -M )(1 -β ) + 10 -0 . 4( M ∗ -M )(1 -γ ) . (10)  \nWe fit our UV LF using Eq. 10, using flat priors: log 10 ( Φ ∗ / Mpc -1 ) ∈ [ -10 , -5], M ∗ ∈ [ -28 , -22], β ∈ [ -3 , -1], γ ∈ [ -6 , -1]. We obtain the best-fit parameters shown in Fig. 12. The fitted parameters are discussed in Sect. 5.2 below.', '4.3.2. Integrated UV emissivity': 'The integral of the UV LF directly yields the total UV emissivity of the AGN population per unit volume. We integrate our best fits to obtain ϵ 1450, the emissivity at 145 nm (i.e., the wavelength where M UV is measured), by extrapolating from M UV = -21 to an arbitrarily luminous magnitude. This allows for a direct comparison with Fig. 8 in Kulkarni et al. (2019) 6 , as we show in our Fig. 13. Our Results are compatible with the fitted trend in Kulkarni et al. (2019) as well as with the integrated UV LFs of several past works in the literature. We computed ϵ 1450 in the same range for our best fits with fixed faint end slope ( β = -1 . 25), obtaining similar values. Therefore, a small change in the faint-end slope has little impact on the total emissivity of intermediate-luminosity quasars.', '4.4. Mean Lymanα forest transmission': 'In Sect. 4.1 we stacked our sample in order to obtain a composite QSO spectrum in the full redshift range ( z = 2 . 7-5.3). Here we repeat the same procedure in bins of redshift. In Fig. 14 we show the stacked spectra in 8 bins of redshift. At wavelengths redward Ly α , our QSO spectra are self-similar at all redshifts the average UV spectrum of quasars is known to not significantly vary across cosmic time (e.g., Onorato et al. 2024). The EW of the most luminous emission lines (Ly α and C iv ) is notably scattered, however this is probably due to fainter sources having higher S / N at lower redshifts, dominating the sample. The uncertainties in the systemic redshift can induce biases in the emission-line EWs of the stacked spectra. In addition, the spectral resolution of the NB pseudo-spectra increases with redshift, given that the rest-frame wavelength interval that is covered by a single NB is proportional to ∝ (1 + z ). For this reason, the shape of the Ly α line might also vary in the composite spectrum, once the latter is brought to the rest-frame (e.g., the Ly α lines appear to be narrower in the two highest redshift bins).  \nFig. 11. UV luminosity function of our sample selected via Ly α emission. We show the LF in the same redshift bins of the Ly α LF (Fig. 8). The filled green squares are used to fit a double power law, while the empty squares are discarded in the fit, due to low completeness (average completeness < 50%). The solid green line and shaded area show the best-fit and 1 σ uncertainties to a double power-law (teal solid line shows the fit with fixed β ). We compare with double power-law fits in the literature, for the redshifts that approximately fall in our bins. The UV LF of the highest redshift bin shown in this figure was computed assuming completeness equal to 1, hence it should be regarded as a lower limit (See also Fig. 8 and Sect. 4.2). We mark the 1 object per bin limit as a horizontal dotted line. Data points below this limit signal LF bins which contain less than 1 object after we take the average of the 1000 LF determinations used to compute our final LF.  \n<!-- image -->  \nOn the other hand, the continuum in the rest-frame wavelength range within 91 . 2 < ( λ 0 / nm) < 121 . 57 shows an evolution with redshift consistent with increasing mean IGM absorption by the Ly α forest ( + Ly β forest for λ 0 < 102 . 5 nm). We average the Ly α forest transmitted flux of the composite spectra in each redshift bin, in the range 113 < ( λ 0 / nm) < 115 (wavelength interval in which no prominent emission lines are expected; see Fig. 14), and divide it by the extrapolated power-law flux fit-  \nove, in order to obtain the mean IGM absorption due to the Ly α forest. We show the mean IGM optical depth ( τ IGM) in Fig. 15, which relates to the IGM transmission as T IGM = e -τ IGM . For this analysis, we only consider redshifts where the observed window used to estimate the transmitted flux is covered at least by the second NB in our set (NB465), to avoid border e ff ects of the stack. Our measurements are in line with previous determinations in the literature (Faucher-Giguère et al. 2008; Becker et al.  \nFig. 12. Best-fit double power-law parameters for the UV LF in bins of redshift. We show the results obtained by leaving all four parameters free (green squares) and fixing the faint-end slope ( β = -1 . 25, shown as a horizontal dashed blue line; blue circles). We compare our results with those of Kulkarni et al. (2019). Our results show a clear decrease in the value of the normalization parameter Φ ∗ with redshift, in line with the fitted Schechter parameter for the Ly α LF. The faint-end slope β also shows hints of a declining evolution, while the bright-end slope γ is poorly constrained by our data.  \n<!-- image -->  \nFig. 13. UV emissivity at rest-frame wavelength 1450 Å. We show the emissivity obtained from the integral of the double power-law fit with four free parameters (green circles) and fixed faint end slope -1.25 (blue circles) - the latter are slightly shifted in the horizontal axis for clarity. The total emissivity at 145 nm rapidly declines with redshift. This quantity is directly related to the ionizing photon emissivity (i.e., λ 0 < 91 . 2 nm), given that the AGN continuum is well approximated by a powerlaw. We compare our results with the integral of other UV LFs in the literature for z > 2 (Matsuoka et al. 2018; Akiyama et al. 2018; Schindler et al. 2019; Kulkarni et al. 2019; Niida et al. 2020; Pan et al. 2022; Liu et al. 2022).  \n<!-- image -->  \n2013; Inoue et al. 2014; Ding et al. 2024) in the range z = 3-5. We recall that the best-fit in Faucher-Giguère et al. (2008) was used in our selection method in order to correct the EWLy α measurement, thus Fig. 15 shows that our selected sample retrieves the same result.', '5. Discussion': 'In this section, we compare our estimated redshift-dependent Ly α and UV LFs with past measurements in the literature, and discuss the evolution of the QSO number densities as a function of cosmic time (Sects. 5.1, 5.2). We comment on the e ffi ciency of our selection method to retrieve a complete and pure sample of Ly α -selected QSOs (Sect. 5.3), and cross-match our sample with spectroscopic surveys (Sect. 5.4).', '5.1. Ly α LF: Evolution and comparison with literature': 'As described in Sect. 4.2.1, we used a Schechter function to fit our Ly α LFs (e.g., Spinoso et al. 2020, TT23). Contrary to previous works, our data probing the faint end ( L Ly α < L ∗ ) enables to sample three Schechter parameters with good significance. Nevertheless, the correlations between these parameters are strong, thus we choose to fix the turnover luminosity parameter L ∗ , in order to have a better understanding of the evolution of α and Φ ∗ . As shown in Fig. 9, the normalization Φ ∗ decreases with redshift - this being more evident for a fixed L ∗ . This result particularly aligns with the evolution of the total SMBH mass function predicted by semi-analytical models (e.g., Fanidakis et al. 2012; Izquierdo-Villalba et al. 2020). In addition, in our results the faint-end slope α evolves toward steeper values with redshift. This trend suggests that, at lower redshift, faint AGNs tend to be proportionally less abundant than the brighter analogs ( L Ly α > L ∗ ). Interestingly, our results about the evolution of Φ ∗ and α appear to be robust against the choice of fixing L ∗ , suggesting that the trends are likely not driven by a shift in the overall luminosity of the AGN population, but instead to be due an evolving balance between faint and bright AGN.  \nIn general, our results are consistent with previous estimations of the AGN Ly α LF in the literature. In particular, the results of this work are in agreement with TT23 within the large intervals of confidence, demonstrating the consistency of our method for di ff erent datasets. The sky area in TT23 is considerably smaller (1.14 deg 2 ) than the combined ∼ 35 deg 2 of the PAUS wide fields, and the depth of both surveys is similar. For this reason, in this work we are able to measure the LF in a range of higher Ly α luminosities; pushing toward the bright end and higher redshifts. The main limiting factor at the LF bright-end is the rapid drop of bright AGNs at increasing redshifts (e.g., Kulkarni et al. 2019; Shen et al. 2020). The discrepancies between this work and some redshift bins in TT23 can be explained by the cosmic variance due to the small area, and higher contamination in the TT23 sample product of a less refined selection method.  \nOur estimate of the binned LF di ff ers from that of Sobral et al. (2018) and Liu et al. (2022) in the intermediate luminosity regime (log 10 ( LLy α ) ∼ 43 . 5-44 erg s -1 ). Sobral et al. (2018) uses a sample selected in deep medium and narrow bands (magnitude ∼ 23 . 5 and ∼ 25 . 5-26 at 5 σ , respectively) within a smaller area of ∼ 2 deg 2 . Their sample contains mainly SFGs, with a small fraction of sources that present radio or X-rays counterparts ( ∼ 1 . 2-4 . 4%), suggestive of AGN activity. Our best-fits for faint end slope (ranging from α = -1 . 45 to -1 . 91; mean  \nFig. 14. Composite QSO spectra in di ff erent bins of redshift. The composite spectra were obtained stacking the objects in our visually selected sample using the stonp code. The UV continuum slope show negligible evolution with redshift. The visible di ff erences in the typical EW of the most prominent lines is likely due to intrinsically fainter sources dominating the lowest redshift bins, and the resolution increase in the rest-frame with redshift. The dotted vertical lines mark the rest-frame wavelengths of the most prominent QSO lines, and the Lyman limit (91 . 2 nm; LyC). The dashed pink line shows the best power-law fit for the UVcontinuum, with slope β UV = -1 . 55 ± 0 . 02. The shaded orange region is used to compute the average IGM optical depth (see Sect. 4.4).  \n<!-- image -->  \nFig. 15. Mean IGM transmission due to the Ly α forest. These estimations were obtained from the stacked composite spectra shown in Fig. 14. As comparison, we show the estimations in Becker et al. (2013), Inoue et al. (2014), Ding et al. (2024) and Faucher-Giguère et al. (2008), the latter employed as continuum correction in our selection method (See Sect. 3.2).  \n<!-- image -->  \n-1 . 59) are in line with the slopes yielded by the results in Sobral et al. (2018) which employed a combination of a Schechter ( α = -1 . 7 + 0 . 3 -0 . 2 ) and a power-law ( A -1 = -1 . 75 ± 0 . 17), where the latter only fitted the Ly α LF of the objects showing AGN radio or X-ray signatures. These fits were obtained from a sample that mainly contains objects in the Ly α luminosity range of 42 . 5 ≲ log 10 ( L Ly α/ erg s -1 ) ≲ 43 . 5 (for comparison, our luminosity range starts at ∼ 43 . 5 erg s -1 ).  \nIn Liu et al. (2022) the Ly α LF is estimated from a sample of objects selected through blind spectroscopy in HETDEX (Gebhardt et al. 2021). Our lya LF at z ∼ 3 . 24 is significantly lower than the one of Liu et al. (2022). Their sample contains fainter (up to r ∼ 26) objects with typically narrower lines, that might contain significant contribution from star-formation activity and Type-II AGN (e.g., Hainline et al. 2011). However, the selection method employed by Liu et al. (2022) relies on the detection of two emission lines inside a limited spectral range, which limits the completeness in specific intervals of redshift where two or more QSO lines are not visible. In addition, their selection method is less e ffi cient in detecting very  \nbroad lines, and thus their sample su ff ers from incompleteness in the brightest Ly α luminosity bins. Our results capture better the exponential Schechter decay of the Ly α LF at the bright end (log 10 ( L Ly α/ erg s -1 ) > 44 . 5).  \nFinally, our bright-end of Ly α LF strongly disagrees with that of Spinoso et al. (2020) in the z = 3 . 24 bin (their highest probed redshift). In Spinoso et al. (2020) the LF is estimated from a sample of NB-selected QSOs in the ∼ 1000 deg 2 of the J-PLUS survey (Cenarro et al. 2019) data-release one. The discrepancies with our estimation of the Ly α LF at z = 3 . 24 can be explained by the fact that, as stated in Spinoso et al. (2020), their z ∼ 3 . 24 QSO sample is highly a ff ected by mis-classification of C iv or C iii QSO lines (at z < 3) as Ly α at z ∼ 3 . 24. In addition, the J-PLUS data used in their work had a limiting depth of r ≲ 21 . 5, which made di ffi cult a reliable identification of z ∼ 3 . 24 objects. Finally, as reported in their work, their purity computation was based on a cross-match with SDSS QSO catalogs, which become increasingly incomplete at z > 3, likely leading to an overestimated purity for the highest-redshift bin they analyzed.', '5.2. UV LF: Evolution and comparison with literature': "The UV LF of quasars has been intensively studied in the literature, mostly focusing on spectroscopically selected samples (e.g., Matsuoka et al. 2018; Akiyama et al. 2018; Kulkarni et al. 2019; Schindler et al. 2019; Niida et al. 2020; Pan et al. 2022). The UV LF clearly evolves across cosmic time, and several models have been proposed to explain this evolution: from 'pure luminosity evolution' (PLE) models, where the change in the LF is mainly driven by a shift in the luminosity distribution; to 'luminosity evolution and density evolution' models where both the luminosity distribution and global number density of AGNs evolve (see e.g., Boyle et al. 2000; Ross et al. 2013).  \nOur best-fitting DPL parameters reflect an evolution of the UVLF(see Fig. 12). The normalization parameter Φ ∗ in general decreases with redshift, in parallel with what is found for our Ly α LF Schechter fits (Sect. 5.1). This decrease is in line with the evolution found by works in the literature that studied the UV LF at di ff erent redshifts, as is also made evident in Fig. 13. This trend also becomes more evident when fixing the faint-end slope. This result is consistent with a global picture in which cosmic AGN activity decreases from z ∼ 2 toward z ∼ 5 (e.g., Shen et al. 2020). When left free, the faint-end slope β shows a correlation with redshift, evolving toward flatter slopes. However, we note that the faint-end of the UV LF might not be well constrained by our data. In general, this parameter does not show  \nconsistent evolution in the literature (see Fig. 12). The other two DPL parameters ( γ, M ∗ ) do not present an evident evolution. In particular the bright-end slope γ shows rather shallow values, as compared to the ones obtained by previous works. The excess in our UV LF is not significant given the uncertainties, and can be explained by cosmic variance. This excess cannot be explained by magnification bias, as this e ff ect is predicted to become relevant at ∼ 8 (Wyithe et al. 2011). In general all the DPL parameters present significant correlations and large variance.  \nIn Akiyama et al. (2018); Niida et al. (2020) and Matsuoka et al. (2018), the UV LF of quasars is computed using a combination of deep imaging data for faint AGNs ( M UV ≳ -26) and SDSS samples of quasars for the bright-end ( M UV ≲ -26). These three works fitted significantly flatter faint-end slopes (mean β = -1 . 25; see Fig. 12, with no significant evolution with redshift. On the other hand, Kulkarni et al. (2019) obtained steeper faint-end slopes, by only employing QSOs from SDSS and BOSS, which in general are more luminous than the sample of HSC deep imaging. In order to try to reconcile our results with the faint-end slope constrained by HSC, we fit a DPL model to our UV LFs using a fixed β = -1 . 25. The fitted curve is displayed in Fig. 11 as a teal solid line, and it is generally compatible with our binned LF, showing that our data does not tightly constrain the faint-end slope. The fits yield slightly fainter values for the turnover magnitude M ∗ , which are more in agreement with the fainter values for M ∗ derived from the HSC + SDSS data. By fixing β = -1 . 25 we also obtain better constrained values of the bright-end slope, and more in line with the HSC + SDSS determinations. Moreover, the anticorrelation of Φ ∗ with redshift increases, reflecting an evolution consistent with that of the Ly α LF (see Sect. 4.2.1). The fact that fixing β = -1 . 25 yields fainter turnover magnitudes might mean that our UV LF only covers a specific M UV range (i.e., M UV ∼ M ∗ ), hence hindering our possibility to well constrain the faint-end of the power-law LF, which is instead well-sampled in the HSC data.", '5.3. Efficiency of the selection method': "The selection method employed in this work (see Sect. 3.2) e ff ectively selects quasars based on strong line emission, selected mainly through NB flux excess. We do not apply any direct color cut in our selection, which removes biases that affect most spectroscopic samples that rely on photometric preselection on color-color diagrams (e.g., Schneider et al. 2010; Ross et al. 2013; Lyke et al. 2020). Selection methods based on QSO variability have been proved to be capable to obtain more pure and complete samples than color-based selections, however loose color cuts are often also applied in order to minimize contamination from stars (e.g., Palanque-Delabrouille et al. 2011). The problems introduced by color cuts especially a ff ect the redshift range 2 . 5 ≲ z ≲ 3 . 5, where most QSOs occupy the sample region as the stellar locus in optical color-color planes (Fan 1999), and the selection functions become ine ffi cient (see e.g., Richards et al. 2006; Ross et al. 2013). Moreover, variability in quasars anti-correlates with continuum luminosity and redshift, and the variability of emission lines is found to be only ∼ 10% of the variability of the continuum (Meusinger et al. 2011). In Schindler et al. (2017) it is shown that the classical QSO selection methods employed by SDSS are prone to miss a significant part of z > 3 quasars in the brightest end, mainly due to the ine ffi ciency of the color-based selection.  \nIn Fig. 16 we show the number counts per unit area in different g -band cuts for 2 . 75 < z < 4 . 3. Our number counts are in line with those of Palanque-Delabrouille et al. (2016) for  \nFig. 16. Di ff erential number counts of QSOs per unit area of sky. The dashed histograms show our complete selected sample, the solid histograms show the distribution of our visually selected 'golden sample'. We compare with the spectroscopic sample of Palanque-Delabrouille et al. (2016) (PD16), selected through QSO variability. Our number counts are significantly higher than those of PD16, especially at the highest bins of redshift.  \n<!-- image -->  \nRedshift  \nz < 3 . 5, from a variability-selected QSO sample. For higher redshifts ( z > 3 . 5), our method appears to be more e ffi cient selecting fainter QSOs ( g > 22 . 5), due to the intrinsic incompleteness of the SDSS data at such magnitudes (the completeness of SDSS imaging data rapidly decays for r > 22; Aihara et al. 2011).", '5.4. Cross-match with spectroscopic surveys': 'We cross-match our catalog of selected LAE candidates with the catalogs of SDSS DR18 (Almeida et al. 2023), DESI EDR (DESI Collaboration et al. 2024) and HETDEX Public Source Catalog 1 (Mentuch Cooper et al. 2023). The match with these catalogs yield a sub-sample of 370 objects with spectroscopic counterpart: 229, 132 and 9 in SDSS, DESI and HETDEX, respectively. The comparison between our measured redshifts and the spectroscopic redshifts measured by the respective surveys is shown in Fig. 5, and discussed in Sect. 3.2.6 above.  \nIn Fig. 17 we show the di ff erential number counts per unit area in our selected sample as compared to the sub-sample having spectroscopic counterparts. While the overlap between the coverage of SDSS DR16 and all three PAUS wide fields is total, there is only partial overlap between those and HETDEX Public Catalog 1 and DESI EDR. Our sample is only slightly more complete than that of SDSS DR16 for log 10 ( L Ly α/ erg s -1 ) ≳ 44 . 25, and contains a significantly larger number of objects for fainter Ly α luminosities ( ∼ 1 dex larger at log 10 ( L Ly α/ erg s -1 ) ≈ 43 . 75). On the other hand, the ratio between the number density of our sample and the QSO sample in the DESI EDR appears to be constant down to log 10 ( L Ly α/ erg s -1 ) ≈ 43 . 5.', '6. Summary': 'In this work we have presented a method to select QSOs through Ly α emission lines in the PAU Survey filter system. This method is a refinement of that developed in Torralba-Torregrosa et al. (2023). We built mock catalogs of the target sample, QSOs with  \n<!-- image -->  \nFig. 17. Di ff erential number counts per square degree as a function of Ly α luminosity (top) and M UV (bottom). We show the counts for the full sample (dashed red line), the visual "golden sample" (solid red line) and the sub-sample with confirmed spectroscopic redshift in either SDSS, HETDEX or DESI (teal solid histograms). Most of the objects in the spectroscopic sub-sample have matching photometric redshifts ( ∼ 90%).  \n<!-- image -->  \nstrong Ly α emission, and the expected contaminants: that is QSOs at redshifts where Ly α cannot be observed with our NB filters, and lowz galaxies.  \nWe probed a total e ff ective area of 35 . 68 deg 2 in three PAUS wide fields (PAUS W1, W2 and W3), and obtained a sample of 915 Ly α selected QSOs, in the redshift range 2 . 7 < z < 5 . 3, with Ly α luminosities of 43 . 5 < log 10 ( L Ly α/ erg s -1 ) < 45 . 5. Our method obtains, in the first place, a sample of candidates selected through prominent flux excess in one or more NBs, that can be compatible with QSO lines at such redshifts. Then we used a neural-network classifier to obtain a purer sample, trained with mock catalogs. This classifier is able to correctly pick up Ly α emitting QSOs from our pre-selected sample, with an accuracy of ≈ 95%, while maintaining a low false-relative rate ( ≈ 7%).  \nA random forest regressor was used, trained with the same mock catalogs in order to obtain precise redshifts. The determination of the redshift using the Ly α detection NB only is slightly inaccurate, since the intrinsic width of the observed wavelength coverage of the NBs correspond to a redshfit uncertainty of ∆ z ≈ 0 . 12. Moreover, since QSO lines can be fairly wide in velocity space (up to ∼ 5000 kms -1 ), the selected NB could not correspond to the peak of the Ly α line, in addition to the expected bias of the QSO Ly α redshift with respect to the systemic redshift. Our algorithm e ff ectively obtains more precise  \nredshifts, when compared to the spectroscopically measured redshifts.  \nThe purity of our sample reaches ≈ 90% at log 10 ( L Ly α/ erg s -1 ) ≈ 44 . 5 for the lowest redshift intervals (2 . 75 < z < 3 . 2) and at log 10 ( L Ly α/ erg s -1 ) ≈ 44 for 3 . 2 < z < 4 . 2, as estimated using our mocks. The purity of the sub-sample with spectroscopic counterparts (cross-matching with SDSS, DESI and HETDEX) is ∼ 88% for the full selected sample, and ≈ 91% ( ≈ 100%) for Ly α luminosity log 10 ( L Ly α/ erg s -1 ) > 44 (log 10 ( L Ly α/ erg s -1 ) > 44 . 5). The main source of contamination, as predicted from our mocks and confirmed by the spectroscopic cross-matched sub-sample, are lowerz QSOs ( z < 2 . 7). In particular, the strong C iv QSO line is the main source of misidentification with Ly α .  \nWe estimated the Ly α and UV ( M UV, at rest-frame 145 nm) LFs. For this task, we obtained estimations of purity and number count corrections over a grid of r synth (a synthetic broad-band from the combination of NBs, mimicking the classical band r ) and L Ly α , in di ff erent intervals of redshift. Both the Ly α and UV LFs show evident evolution, with declining normalization parameter Φ ∗ , fitting a Schechter function, and a double power-law, respectively. The faint-end of the Ly α LF shows a clear correlation with redshift in the probed interval, suggesting that the number density of bright quasars declines faster than the faint counterparts, in agreement with the SMBH growth models at z ≳ 2.  \nWe obtained composite QSO spectra from NB data, for different intervals of redshift. We measured a UV continuum slope of β UV = -1 . 55 ± 0 . 02. We employ the composite spectra to measure the mean IGM transmision of the Ly α forest, obtaining results compatible with previous spectroscopic determinations.  \nThe total Ly α luminosity density, obtained from the integral of the Ly α LF, shows a rapid decline with increasing redshift, being sub-dominant by several orders of magnitude already by z ≈ 4, with respect to that of star-formation Ly α emission.', 'Data availability': "The catalog of selected objects elaborated in this work will be available at CosmoHub (Tallada et al. 2020; Carretero et al. 2017), together with the whole PAUS public data release, through this link: https: // cosmohub.pic.es / catalogs / 319. More information on the data release, as well as access to the PAUS database and raw and reduced images, can be found at https: // pausurvey.org / public-data-release / . The stacking code stonp is publicly available at https: // github.com / PAUsurvey / stonp.  \nAcknowledgements. The authors acknowledge the insightful feedback from the scientific referee, which has significantly contributed to enhancing the quality of this paper. This work has been funded by project PID2019109592GBI00 / AEI / 10.13039 / 501100011033 from the Spanish Ministerio de Ciencia e Innovación (MCIN)-Agencia Estatal de Investigación, by the Project of Excellence Prometeo / 2020 / 085 from the Conselleria d'Innovació Universitats, Ciència i Societat Digital de la Generalitat Valenciana. This work is part of the research Project PID2023-149420NB-I00 funded by MICIU / AEI / 10.13039 / 501100011033 and by ERDF / EU. This work is also supported by the project of excellence PROMETEO CIPROM / 2023 / 21 of the Conselleria de Educación, Universidades y Empleo (Generalitat Valenciana). The authors acknowledge the financial support from the MCIN with funding from the European Union NextGenerationEU and Generalitat Valenciana in the call Programa de Planes Complementarios de I + D + i (PRTR 2022). Project (VALJPAS), reference ASFAE / 2022 / 025. PR and DS acknowledge the support by the Tsinghua Shui Mu Scholarship. PR, DS and ZC acknowledge the funding of the National Key R&D Program of China (grant no. 2023YFA1605600), the science research grants from the China Manned Space Project with No. CMS-CSST2021-A05, and the Tsinghua University Initiative Scientific Research Program (No. 20223080023). PR acknowledges additional funding from  \nthe National Science Foundation of China (grant no. 12350410365). EG acknowledges grants from Spain Plan Nacional (PGC2018-102021-B-100) and Maria de Maeztu (CEX2020-001058-M). The PAU Survey is partially supported by MINECO under grants CSD2007-00060, AYA2015-71825, ESP201789838, PGC2018-094773, PGC2018-102021, PID2019-111317GB, SEV-20160588, SEV-2016-0597, MDM-2015-0509 and Juan de la Cierva fellowship and LACEGALand EWC Marie Sklodowska-Curie grant No 734374 and no.776247 with ERDF funds from the EU Horizon 2020 Programme, some of which include ERDF funds from the European Union. IEEC and IFAE are partially funded by the CERCA and Beatriu de Pinos program of the Generalitat de Catalunya. Funding for PAUS has also been provided by Durham University (via the ERC StG DEGAS-259586), ETH Zurich, Leiden University (via ERC StG ADULT-279396 and Netherlands Organisation for Scientific Research (NWO) Vici grant 639.043.512), University College London and from the European Union's Horizon 2020 research and innovation programme under the grant agreement No 776247 EWC. The PAU data center is hosted by the Port d'Informació Científica (PIC), maintained through a collaboration of CIEMAT and IFAE, with additional support from Universitat Autònoma de Barcelona and ERDF. We acknowledge the PIC services department team for their support and fruitful discussions. EG acknowledges grants from Spain Plan Nacional (PGC2018-102021-B-100) and Maria de Maeztu (CEX2020-001058-M). JC acknowledges support from the grant PID2021-123012NA-C44 funded by MCIN / AEI / 10.13039 / 501100011033 and by 'ERDF A way of making Europe'. FJC is supported by grants PID2022-141079NB and CEX2020-001058-M funded by MCIN / AEI / 10.13039 / 501100011033. ME acknowledges funding by MCIN with funding from European Union NextGenerationEU (PRTR-C17.I1) and by Generalitat de Catalunya. IFT participation is supported by the grant PID2021-123012NB-C43P funded by MCIN / AEI / 10.13039 / 501100011033. H. Hildebrandt is supported by a DFG Heisenberg grant (Hi 1495 / 5-1), the DFG Collaborative Research Center SFB1491, as well as an ERC Consolidator Grant (No. 770935). IFAE participation is supported by grant PID2021123012NB funded by MCIN / AEI / 10.13039 / 501100011033. CIEMAT participation is supported by the grant PID2021-123012NB-C42P funded by MCIN / AEI / 10.13039 / 501100011033.", 'References': "Adams, J. J., Blanc, G. A., Hill, G. J., et al. 2011, ApJS, 192, 5 Adams, N. J., Bowler, R. A. A., Jarvis, M. J., Varadaraj, R. G., & Häußler, B. 2023, MNRAS, 523, 327 Aihara, H., Allende Prieto, C., An, D., et al. 2011, ApJS, 193, 29 Akins, H. B., Casey, C. M., Lambrides, E., et al. 2024, Submitted to ApJ, arXiv:2406.10341 Akiyama, M., He, W., Ikeda, H., et al. 2018, PASJ, 70, S34 Alarcon, A., Gaztanaga, E., Eriksen, M., et al. 2021, MNRAS, 501, 6103 Almeida, A., Anderson, S. F., Argudo-Fernández, M., et al. 2023, ApJS, 267, 44 Arrabal Haro, P., Rodríguez Espinosa, J. M., Muñoz-Tuñón, C., et al. 2020, MNRAS, 495, 1807 Askar, A., Davies, M. B., & Church, R. P. 2022, MNRAS, 511, 2631 Bacon, R., Brinchmann, J., Richard, J., et al. 2015, A&A, 575, A75 Becker, G. D., Hewett, P. C., Worseck, G., & Prochaska, J. X. 2013, MNRAS, 430, 2067 Blanc, G. A., Adams, J. J., Gebhardt, K., et al. 2011, ApJ, 736, 31 Bonoli, S., Calderone, G., Abramo, R., et al. 2021, in Nuclear Activity in Galaxies Across Cosmic Time, ed. M. Povi'c, P. Marziani, J. Masegosa, H. Netzer, S. H. Negu, & S. B. Tessema, Vol. 356, 12-16 Boyle, B. J., Shanks, T., Croom, S. M., et al. 2000, MNRAS, 317, 1014 Byrohl, C. & Nelson, D. 2023, Mon. Not. R. Astron. Soc., 523, 5248 Calhau, J., Sobral, D., Santos, S., et al. 2020, MNRAS, 493, 3341 Carretero, J., Tallada, P., Casals, J., et al. 2017, in Proceedings of the European Physical Society Conference on High Energy Physics. 5-12 July, 488 Cenarro, A. J., Moles, M., Cristóbal-Hornillos, D., et al. 2019, A&A, 622, A176 Chabrier, G. 2003, PASP, 115, 763 Compostella, M., Cantalupo, S., & Porciani, C. 2014, MNRAS, 445, 4186 Cowie, L. L. & Hu, E. M. 1998, AJ, 115, 1319 Croom, S. M., Richards, G. T., Shanks, T., et al. 2009, MNRAS, 399, 1755 Cuillandre, J.-C. J., Withington, K., Hudelot, P., et al. 2012, in Society of PhotoOptical Instrumentation Engineers (SPIE) Conference Series, Vol. 8448, Observatory Operations: Strategies, Processes, and Systems IV, ed. A. B. Peck, R. L. Seaman, & F. Comeron, 84480M Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827 Dawson, S., Rhoads, J. E., Malhotra, S., et al. 2007, ApJ, 671, 1227 Dayal, P., Volonteri, M., Choudhury, T. R., et al. 2020, MNRAS, 495, 3065 de La Vieuville, G., Bina, D., Pello, R., et al. 2019, A&A, 628, A3 DESI Collaboration, Adame, A. G., Aguilar, J., et al. 2024, AJ, 168, 58 Ding, J., Madau, P., & Prochaska, J. X. 2024, MNRAS, 532, 2082 Drake, A. B., Garel, T., Wisotzki, L., et al. 2017, A&A, 608, A6  \nArticle number, page 18 of 22  \nDriver, S. P., Bellstedt, S., Robotham, A. S. G., et al. 2022, MNRAS, 513, 439 Dubois, Y., Peirani, S., Pichon, C., et al. 2016, MNRAS, 463, 3948 Eriksen, M., Alarcon, A., Cabayol, L., et al. 2020, MNRAS, 497, 4565 Fabian, A. C. 2012, ARA&A, 50, 455 Fan, X. 1999, AJ, 117, 2528 Fanidakis, N., Baugh, C. M., Benson, A. J., et al. 2012, MNRAS, 419, 2797 Faucher-Giguère, C.-A., Prochaska, J. X., Lidz, A., Hernquist, L., & Zaldarriaga,  \nM. 2008, ApJ, 681, 831  \nGebhardt, K., Mentuch Cooper, E., Ciardullo, R., et al. 2021, ApJ, 923, 217  \nGlikman, E., Bogosavljevi'c, M., Djorgovski, S. G., et al. 2010, ApJ, 710, 1498  \nGlikman, E., Djorgovski, S. G., Stern, D., et al. 2011, ApJ, 728, L26  \n- Greene, J. E., Labbe, I., Goulding, A. D., et al. 2024, ApJ, 964, 39  \nGreig, B., Mesinger, A., McGreer, I. D., Gallerani, S., & Haiman, Z. 2016, MN-  \nRAS, 466, 1814  \n- Greig, B., Mesinger, A., McGreer, I. D., Gallerani, S., & Haiman, Z. 2017, MNRAS, 466, 1814  \nGuaita, L., Gawiser, E., Padilla, N., et al. 2010, ApJ, 714, 255  \n- Hainline, K. N., Shapley, A. E., Greene, J. E., & Steidel, C. C. 2011, ApJ, 733, 31\n- Harris, D. W., Jensen, T. W., Suzuki, N., et al. 2016, AJ, 151, 155  \nHasinger, G., Miyaji, T., & Schmidt, M. 2005, A&A, 441, 417  \n- Henriques, B. M. B., White, S. D. M., Thomas, P. A., et al. 2015, MNRAS, 451, 2663\n- Hernán-Caballero, A., Akhlaghi, M., López-Sanjuan, C., et al. 2024, A&A, 684, A61\n- Hernán-Caballero, A., Varela, J., López-Sanjuan, C., et al. 2021, A&A, 654, A101\n- Hernán-Caballero, A., Willmer, C. N. A., Varela, J., et al. 2023, A&A, 671, A71 Heymans, C., Van Waerbeke, L., Miller, L., et al. 2012, Mon. Not. R. Astron. Soc., 427, 146  \nHopkins, P. F., Richards, G. T., & Hernquist, L. 2007, ApJ, 654, 731 Hu, E. M., Cowie, L. L., & McMahon, R. G. 1998, ApJ, 502, L99 Hu, W., Wang, J., Zheng, Z.-Y., et al. 2019, ApJ, 886, 90 Inoue, A. K., Shimizu, I., Iwata, I., & Tanaka, M. 2014, MNRAS, 442, 1805 Ivashchenko, G., Sergijenko, O., & Torbaniuk, O. 2014, MNRAS, 437, 3343 Izquierdo-Villalba, D., Angulo, R. E., Orsi, A., et al. 2019, A&A, 631, A82 Izquierdo-Villalba, D., Bonoli, S., Dotti, M., et al. 2020, MNRAS, 495, 4681 Izquierdo-Villalba, D., Sesana, A., Bonoli, S., & Colpi, M. 2022, MNRAS, 509, 3488  \nJiang, L., McGreer, I. D., Fan, X., et al. 2016, ApJ, 833, 222 Karman, W., Caputi, K. I., Grillo, C., et al. 2015, A&A, 574, A11 Khusanova, Y., Le Fèvre, O., Cassata, P., et al. 2020, A&A, 634, A97 Kocevski, D. D., Finkelstein, S. L., Barro, G., et al. 2024, Submitted to ApJ, arXiv:2404.03576  \nKokorev, V., Caputi, K. I., Greene, J. E., et al. 2024, ApJ, 968, 38  \nKonno, A., Ouchi, M., Shibuya, T., et al. 2018, PASJ, 70, S16  \nKuijken, K., Heymans, C., Dvornik, A., et al. 2019, A&A, 625, A2  \nKulkarni, G., Worseck, G., & Hennawi, J. F. 2019, MNRAS, 488, 1035  \nLabbé, I., van Dokkum, P., Nelson, E., et al. 2023, Nature, 616, 266  \nLaur, J., Tempel, E., Tamm, A., et al. 2022, A&A, 668, A8  \nLin, X., Zheng, Z., & Cai, Z. 2022, ApJS, 262, 38  \nLiu, C., Gebhardt, K., Cooper, E. M., et al. 2022, ApJ, 935, 132  \nLusso, E., Worseck, G., Hennawi, J. F., et al. 2015, MNRAS, 449, 4204  \nLyke, B. W., Higley, A. N., McLane, J. N., et al. 2020, ApJS, 250, 8  \nMadau, P. & Dickinson, M. 2014, ARA&A, 52, 415  \nMalhotra, S. & Rhoads, J. E. 2004, ApJ, 617, L5  \nMartí, P., Miquel, R., Castander, F. J., et al. 2014, MNRAS, 442, 92  \nMartin, C. L. & Sawicki, M. 2004, ApJ, 603, 414  \nMatsuoka, Y., Strauss, M. A., Kashikawa, N., et al. 2018, ApJ, 869, 150  \nMatthee, J., Naidu, R. P., Brammer, G., et al. 2024, ApJ, 963, 129  \nMatthee, J., Sobral, D., Best, P., et al. 2017, MNRAS, 471, 629  \nMcGreer, I. D., Jiang, L., Fan, X., et al. 2013, ApJ, 768, 105  \nMentuch Cooper, E., Gebhardt, K., Davis, D., et al. 2023, ApJ, 943, 177  \nMeusinger, H., Hinze, A., & de Hoon, A. 2011, A&A, 525, A37  \nMorales, A. M., Mason, C. A., Bruton, S., et al. 2021, ApJ, 919, 120  \nMountrichas, G. 2023, A&A, 672, A98  \n- Nakajima, K., Fletcher, T., Ellis, R. S., Robertson, B. E., & Iwata, I. 2018, MNRAS, 477, 2098  \nNavarro-Gironés, D., Gaztañaga, E., Crocce, M., et al. 2023, Submitted to MN-  \nRAS, arXiv:2312.07581  \nNiida, M., Nagao, T., Ikeda, H., et al. 2020, ApJ, 904, 89  \nOke, J. B. & Gunn, J. E. 1983, ApJ, 266, 713  \nOno, Y., Itoh, R., Shibuya, T., et al. 2021, ApJ, 911, 78  \n- Onorato, S., Hennawi, J. F., Schindler, J.-T., et al. 2024, Submitted to MNRAS, arXiv:2406.07612  \nOrsi, Á., Padilla, N., Groves, B., et al. 2014, MNRAS, 443, 799  \nOuchi, M., Shimasaku, K., Akiyama, M., et al. 2008, ApJS, 176, 301  \nPadilla, C., Castander, F. J., Alarcón, A., et al. 2019, AJ, 157, 246  \nPalanque-Delabrouille, N., Magneville, C., Yèche, C., et al. 2016, A&A, 587,  \nA41  \nPalanque-Delabrouille, N., Yeche, C., Myers, A. D., et al. 2011, A&A, 530, A122 Pan, Z., Jiang, L., Fan, X., Wu, J., & Yang, J. 2022, ApJ, 928, 172 Partridge, R. B. & Peebles, P. J. E. 1967, ApJ, 147, 868 Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825  \nPlanck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6 Pritchet, C. J. 1994, PASP, 106, 1052  \n- Queiroz, C., Abramo, L. R., Rodrigues, N. V. N., et al. 2023, MNRAS, 520, 3476\n- Renard, P., Gaztanaga, E., Croft, R., et al. 2021, MNRAS, 501, 3883\n- Renard, P., Siudek, M., Eriksen, M. B., et al. 2022, MNRAS, 515, 146\n- Renard, P., Spinoso, D., Sun, Z., et al. 2024, arXiv e-prints, arXiv:2406.18775\n- Richards, G. T., Strauss, M. A., Fan, X., et al. 2006, AJ, 131, 2766\n- Rollinde, E., Theuns, T., Schaye, J., Pâris, I., & Petitjean, P. 2013, MNRAS, 428, 540\n- Ross, N. P., McGreer, I. D., White, M., et al. 2013, ApJ, 773, 14  \nSantos, S., Sobral, D., Butterworth, J., et al. 2021, MNRAS, 505, 1117  \n- Santos, S., Sobral, D., & Matthee, J. 2016, MNRAS, 463, 1678\n- Santos, S., Sobral, D., Matthee, J., et al. 2020, MNRAS, 493, 141\n- Schechter, P. 1976, ApJ, 203, 297\n- Schindler, J.-T., Fan, X., McGreer, I. D., et al. 2019, ApJ, 871, 258  \nSchindler, J.-T., Fan, X., McGreer, I. D., et al. 2017, ApJ, 851, 13  \nSchmidt, M. 1968, ApJ, 151, 393  \nSchneider, D. P., Richards, G. T., Hall, P. B., et al. 2010, AJ, 139, 2360  \nScoville, N., Aussel, H., Brusa, M., et al. 2007, ApJS, 172, 1  \n- Selsing, J., Fynbo, J. P. U., Christensen, L., & Krogager, J. K. 2016, A&A, 585, A87\n- Serrano, S., Gaztañaga, E., Castander, F. J., et al. 2023, MNRAS, 523, 3287 Shen, X., Hopkins, P. F., Faucher-Giguère, C.-A., et al. 2020, MNRAS, 495, 3252\n- Shull, J. M., Stevans, M., & Danforth, C. W. 2012, ApJ, 752, 162 Sobral, D., Santos, S., Matthee, J., et al. 2018, MNRAS, 476, 4725\n- Spinoso, D., Bonoli, S., Valiante, R., Schneider, R., & Izquierdo-Villalba, D. 2023, MNRAS, 518, 4672\n- Spinoso, D., Orsi, A., López-Sanjuan, C., et al. 2020, A&A, 643, A149 Springel, V., White, S. D. M., Jenkins, A., et al. 2005, Nature, 435, 629 Tallada, P., Carretero, J., Casals, J., et al. 2020, Astronomy and Computing, 32,\n- 100391\n- Thai, T. T., Tuan-Anh, P., Pello, R., et al. 2023, A&A, 678, A139\n- Torralba-Torregrosa, A., Gurung-López, S., Arnalte-Mur, P., et al. 2023, A&A, 680, A14\n- Trebitsch, M., Dubois, Y., Volonteri, M., et al. 2021, A&A, 653, A154 Tremmel, M., Ricarte, A., Natarajan, P., et al. 2024, The Open Journal of Astrophysics, 7, 26\n- Vanden Berk, D. E., Richards, G. T., Bauer, A., et al. 2001, AJ, 122, 549 Wang, F., Yang, J., Fan, X., et al. 2021, ApJ, 907, L1\n- Ward, S. R., Harrison, C. M., Costa, T., & Mainieri, V. 2022, MNRAS, 514, 2936 White, R. L., Becker, R. H., Fan, X., & Strauss, M. A. 2003, AJ, 126, 1 Wisotzki, L. 2000, A&A, 353, 861\n- Wold, I. G. B., Malhotra, S., Rhoads, J., et al. 2022, ApJ, 927, 36\n- Worseck, G., Prochaska, J. X., Hennawi, J. F., & McQuinn, M. 2016, ApJ, 825, 144  \nWyithe, J. S. B., Yan, H., Windhorst, R. A., & Mao, S. 2011, Nature, 469, 181 Xu, C., Smith, A., Borrow, J., et al. 2023, MNRAS, 521, 4356 Yang, G., Brandt, W. N., Alexander, D. M., et al. 2019, MNRAS, 485, 3721 Yang, J., Wang, F., Wu, X.-B., et al. 2016, ApJ, 829, 33 Zhang, Y., Ouchi, M., Gebhardt, K., et al. 2021, ApJ, 922, 167  \nTable D.1. Best-fit Schechter parameters for the Ly α LF.  \nIn Fig. B.1 we show the confusion matrices of the three classes, in diferent bins of measured Ly α luminosity.  \nFor the RF redshift regressor we use the method RandomForestRegressor , and perform a grid-search using the same method employed for the NN. We obtain the following optimal parameters:  \nmin\\_samples\\_split = 3,  \nmin\\_samples\\_leaf = 2,  \nmax\\_depth = 20,  \nbootstrap = True", 'Appendix C: 2D purity and completeness': "In Fig. C.1 we display the purity and completeness of the selection method in bins of r and measured L Ly α . These 2D maps are obtained by applying the selection method on the mocks (Sect. 3.1) and comparing the selected sample with the parent sample using Eqs. 4 and 5.  \nTable A.1. NB filters used to select Ly α emitting quasars.  \nNotes. a The name of the filters is composed of 'NB' + the approximate e ff ective wavelength in nm. b Average limit magnitude at the 5 σ level. c Interval of Ly α redshifts probed. The amplitude of this interval is defined by the FWHM of the filters. d Number of selected Ly α emitting quasars in the final sample by each NB.", 'Appendix A: Ly α redshift of the PAUS narrow bands': 'In Table A.1 we display the list of NBs used to select Ly α emitting candidates, the associated Ly α redshift of each NB and the number of candidates selected.', 'Appendix B: ANN classifier details': "For the machine learning models used in this work, we use the python library Sci-kit learn (Pedregosa et al. 2011). The ANN classifier is built using the MLPClassifier method. We performed a grid-search of optimal parameters with GridSeachCV , over an extense grid of parameters. The ANN source classifier optimal parameters are:  \nhidden\\_layer\\_sizes = (60, 60, 60),  \n```\nbatch\\_size = 64, alpha = 1e-4, learning\\_rate = 'adaptive', max\\_iter = 1e4, solver = 'adam' .\n```", 'Appendix D: Schechter and double power-law fit parameters': "In Tables D.1 and D.2 we list the best-fit parameters to a Schechter function for the Ly α LF, and double powerlaw for the UV LF, respectively. For the Ly α LF Schechter fit we list the results of fitting with all three free Schechter parameters, and fixing log 10 ( L ∗ / erg s -1 ) = 45, as discussed in Sect. 4.2.1. Similarly, for the UV LF fit, we list the results of leaving all four DPL parameters free, and fixing the faint-end slope β = -1 . 25, as discussed in Sect. 4.3.1  \nFig. B.1. Confusion matrix for the ANN classifier used to select LAEs. The three classes considered are: QSOs selected by a feature other than Ly α ('QSO cont.'), QSOs correctly selected by their Ly α line, in the correct redshift ('QSO LAE'), and lowz ( z ≈ 0-2) galaxy contaminants ('Low -z galaxy'). The annotations in the cells are the normalized counts so the vertical columns add up to 1. We show the confusion matrices for the full test set (top, L Ly α > 10 43 . 5 erg s -1 ) and splitting in low (middle, L Ly α < 10 44 erg s -1 ) and high luminosity bins (bottom, L Ly α > 10 44 erg s -1 ). This figure shows that the ANN e ff ectively selects QSO LAEs with an accuracy of 93% for the general sample, while maintaining a low false negative rate (6% and 1% in the case of objects labeled as QSO cont. and lowz galaxies, respectively.)  \n<!-- image -->  \nTable D.2. Best-fit double power-law parameters for the UV LF.  \nFig. C.1. 2D maps of purity and completeness in bins of L Ly α and magnitude r . The values displayed in Fig. 6 are obtained by collapsing these 2D maps over the horizontal axis.  \n<!-- image -->"}