© ESO 2020

1. Introduction

An increasing number of recent works focus on the study of high-redshift Lyman-α emitters (LAEs), objects showing prominent rest-frame Lyα emission within a spectrum (usually) devoid of other line features (e.g. Cassata et al. 2011; Nakajima et al. 2018). The spectral properties of LAEs are usually interpreted as coming from young (≲50 Myr) and low-mass (M_*​ <  ​10¹⁰ M_⊙) galaxies (e.g. Wilkins et al. 2011; Amorín et al. 2017; Hao et al. 2018; Santos et al. 2020) with small rest-frame UV half-light radii (R​ ≲ ​1 − 2 Kpc, as in e.g. Møller & Warren 1998; Lai et al. 2008; Bond et al. 2012; Guaita et al. 2015; Kobayashi et al. 2016; Ribeiro et al. 2016; Bouwens et al. 2017a; Paulino-Afonso et al. 2018) which are actively star forming (SFR​ ∼ ​1 − 100 M_⊙ yr⁻¹) and dust poor (dust attenuation A_V​ <  ​0.2, see e.g. Gawiser et al. 2006, 2007; Guaita et al. 2011; Nilsson et al. 2011; Bouwens et al. 2017b; Arrabal Haro et al. 2020). When observed at high redshift, isolated and grouped LAEs appear to represent the progenitors of present-day galaxies and clusters, respectively, providing extremely valuable insights into structure formation (e.g. Matsuda et al. 2004, 2005; Venemans et al. 2005; Gawiser et al. 2007; Overzier et al. 2008; Guaita et al. 2010; Mei et al. 2015; Bouwens et al. 2017b; Khostovan et al. 2019). A basic statistical tool to study the population of high-z LAEs is the description of their number density at a given redshift as a function of line luminosity (L_Lyα), namely the Lyα luminosity function (LF; see e.g. Gronke et al. 2015, for a theoretical approach). Several recent works have focused on the construction of the Lyα LF at z​ ≥ ​2 (Gronwall et al. 2007; Ouchi et al. 2008; Blanc et al. 2011; Clément et al. 2012; Konno et al. 2016; Sobral et al. 2017, 2018a) by making use of deep observations of narrow sky regions (up to few squared degrees, as in e.g. Matthee et al. 2014, 2017b; Cassata et al. 2015; Ono et al. 2018). Their findings describe a Lyα LF that follows a Schechter function (Schechter 1976) at relatively faint line luminosity (i.e. L_Lyα ≲ 10^42.5, see e.g. Ouchi et al. 2008; Konno et al. 2016; Sobral et al. 2016; Matthee et al. 2017a), a regime mostly occupied by low-mass star-forming galaxies (e.g. Hu et al. 1998; Kudritzki et al. 2000; Stiavelli et al. 2001; Santos et al. 2004; van Breukelen et al. 2005; Gawiser et al. 2007; Rauch et al. 2008; Guaita et al. 2011).

The bright end of the Lyα LF is instead populated by rare, bright, and actively star-forming Lyα-emitting systems (e.g. Matsuda et al. 2011; Bridge et al. 2013; Cai et al. 2017a, 2018) as well as active galactic nuclei and quasi-stellar object (AGN/QSOs; e.g. Sobral et al. 2018a; Calhau et al. 2020). Sources belonging to this latter class usually show multiple emission lines, including Lyα. The production and escape of Lyα radiation from AGN/QSOs results from the reprocessing of ionising fluxes produced by the accretion of matter onto their central massive black holes (see e.g. Haiman & Rees 2001; Bolton et al. 2011; Dijkstra 2017). The absorption of Lyα photons by the intra-galactic medium (IGM) surrounding AGN/QSOs and along their line of sight is used to infer the properties of the gaseous medium surrounding these sources (e.g. Wyithe et al. 2005; Bolton et al. 2011; North et al. 2012) and the progress of cosmic reionisation (e.g. Fan et al. 2006a,b; Becker et al. 2015; Eilers et al. 2018). Consequently, the interest in the Lyα LF bright end is motivated by both the search for bright and extreme star-formation-powered Lyman-α emitters (LAEs; e.g. Sobral et al. 2015, 2016; Hartwig et al. 2016; Shibuya et al. 2018) and the characterisation of high-z AGN/QSO populations (e.g. Konno et al. 2016; Sobral et al. 2018a; Matthee et al. 2017b; Calhau et al. 2020).

Current constraints on the Lyα LF at high luminosity are somewhat weak, given the relatively small cosmological volumes probed by past works (e.g. Fujita et al. 2003; Blanc et al. 2011; Herenz et al. 2019). In particular, recent measurements reveal a possible number-density excess with respect to an exponential (Schechter) decay at L_Lyα​ ≳ ​10⁴³ (e.g. Konno et al. 2016), which is possibly due to the contribution of a faint-AGN population (see e.g. Matthee et al. 2017b; Sobral et al. 2018a). This is supported by the tomographic analysis performed by Sobral et al. (2018a) in the COSMOS field using a combination of optical, infrared, and X-Ray data. The work of these latter authors suggests that the major contribution to the LF at L_Lyα​ ≳ ​10⁴³ is provided by X-ray-loud sources (i.e. likely AGNs; see also Matthee et al. 2017b; Calhau et al. 2020). This contribution completely vanishes at z​ ≳ ​3.5 in the data of Sobral et al. (2018a), thus paralleling the peak of AGN activity usually observed at z​ ∼ ​2​ − ​3 (e.g. Hasinger et al. 2005; Miyaji et al. 2015).

The current study exploits the first data release (DR1) of the Javalambre Photometric Local Universe Survey (J-PLUS, Cenarro et al. 2019), which provides imaging of the northern hemisphere in both narrow- and broad-bands (NBs and BBs; see Fig. 1 and Table 1). The DR1 covers an area of > 1000 deg², which is unprecedented for NB surveys of z​ >  ​2 luminous line emitters. Our goal is to exploit these characteristics in order to obtain large samples of photometrically selected Lyα-emitting sources, and probe the Lyα LF bright end at four different redshifts (see Table 2). Indeed, the combination of large survey area and multi-NB data allows us to sample the distribution of bright Lyα-emitting sources over a luminosity regime that remains poorly constrained (Gronwall et al. 2007; Guaita et al. 2010; Blanc et al. 2011; Konno et al. 2016; Sobral et al. 2018a). We complement our study with the results of two follow-up spectroscopic programs that assess the performance of our method.

Fig. 1. Measured transmission curves for the J-PLUS filter set after accounting for sky absorption, CCD quantum efficiency, and the total effect of the JAST/T80 telescope optical system. The four NB we exploit to look for bright Lyα emitters at z​ >  ​2 (namely, the J0395, J0410, J0430 and J0515 NBs) share their wavelength coverage with the g band and are shown here as filled-area curves.

This paper is organised as follows: Sect. 2 details the main features of the J-PLUS survey and the classes of the targeted sources. Our method for detecting NB excesses, our selection function, and our sample of LAE candidates are described in Sect. 3, along with our spectroscopic follow-up programs. Section 4 focuses on the computation of the four 2​ <  ​z​ <  ​3.3 Lyα LFs. Finally, we discuss our results in Sect. 5 and present our conclusions in Sect. 6. Throughout this paper, magnitudes are given in the AB system (Oke 1974; Oke & Gunn 1983), and we assume a flat ΛCDM cosmology described by PLANCK15 parameters (Planck Collaboration I 2016; Planck Collaboration XIII 2016), namely: H₀ = 67.3 Km s⁻¹ Mpc⁻¹, Ω_m; 0 = 0.315, Ω_Λ; 0 = 0.685.

2. Lyα emitters in the J-PLUS photometric survey

J-PLUS is an ongoing wide-area photometric survey performed at the Observatorio Astrofísico de Javalambre (OAJ, Cenarro et al. 2014) in Arcos de las Salinas (Teruel, Spain). Here we summarise its technical features (detailed in Cenarro et al. 2019) and we define the class of Lyα-emitting sources that we target.

2.1. Survey description and source catalogues

J-PLUS observations are being carried out by the T80Cam instrument on the JAST/T80 83 cm telescope (Marin-Franch et al. 2015). The JAST/T80 optical system provides a wide field of view (FoV​ ∼ ​1.96 deg²) while ensuring a high spatial resolution (0.55 arcsec pixel⁻¹, see Cenarro et al. 2019, for technical details). J-PLUS nominal depth is shallower than that of comparably wide optical surveys, namely r = 22 at a signal-to-noise ratio (S/N) = 3 (as compared to e.g. r′ = 23.1 at S/N = 5 for SDSS; see York et al. 2000). Nevertheless, it offers NB measurements over an unprecedentedly wide sky area, making it suitable for extensive searches for bright emission-line galaxies (ELGs). The J-PLUS filter set is composed of 12 pass-band filters (see Fig. 1) divided into five broad-bands and seven narrow-bands of width ∼ 800 − 2000 Å and ∼ 150 − 500 Å, respectively (Table 1). Their transmission curves are shown in Fig. 1, after accounting for optical elements, CCD quantum efficiency, and sky transparency.

J-PLUS images are automatically reduced to obtain public catalogues of sources¹. This work is based on the recent DR1, which was obtained on 511 pointings with stable pipelines for data reduction and source extraction specifically calibrated and tested on J-PLUS data (see e.g. Cenarro et al. 2019; López-Sanjuan et al. 2019a). We use the J-PLUS dual-mode catalogue, constructed with r as the band for source detection, sky positions, and photometric apertures definition for all the remaining filters. In particular, auto-aperture² photometry ensures recovery of the total Lyα line flux of point-like sources (see Sect. 4.1) and allows us to exploit the tabulated detection-completeness of each pointing (see Sect. 4.4). We note that relying on dual-mode catalogues has non-trivial implications on the completeness of our final samples. We address this issue in Sect. 4.4.3.

2.2. Targeting Lyα emission with J-PLUS

The design of J-PLUS filters potentially allows us to detect Lyα emission within seven redshift windows, one per NB, respectively centred at z​ ∼ ​2.11,  2.24,  2.38,  2.54,  3.23,  4.43 and 6.09. In particular, we employ the J0395, J0410, J0430 and J0515 filters (see Fig. 1) for targeting z​ ∼ ​2.24,  2.38,  2.54 and 3.23 (Table 2). Our selection is based on measuring NB excesses with respect to the continuum traced by BB photometry (see Sect. 3.1). Consequently, it is prone to contamination by prominent emission lines. In particular, we expect our samples to be contaminated by both nebular emission from star formation (e.g. H _β, [OIII]​_{4959 + 5007} and [OII]​₃₇₂₇ lines) and AGN/QSO ionising radiation (e.g. CIV​₁₅₄₉, CIII]​₁₉₀₈, MgII​₂₇₉₉ and SiIV​₁₃₉₇ lines, see also Stroe et al. 2017a,b). The latter lines and their associated redshift intervals in J-PLUS are listed in Table 2. We note that SiIV and MgII are minor sources of contamination because: (i) they are significantly fainter than Lyα (e.g. Telfer et al. 2002; Selsing et al. 2016), (ii) J-PLUS probes relatively small cosmological volumes at 0.35​ <  ​z​ <  ​0.85, and (iii) the number density of AGNs/QSOs at z​ <  ​1 is lower than at z​ >  ​2 (e.g. Palanque-Delabrouille et al. 2016; Pâris et al. 2018). We exclude the J0378 NB after verifying that our method does not reliably detect photometric excess in this NB (see Sect. 4.1). We also exclude the J0660 and J0861 NBs because they provide very scarce samples of candidates (≲100 sources) whose contamination cannot be reliably estimated, which is due to the absence of SDSS spectroscopic counterparts (see Sect. 3.3). We note that this is in agreement with the work of Sobral et al. (2018a), which shows no significant detection of bright (L_Lyα ≳ 10⁴³ erg s⁻¹) Lyα-emitting sources at z​ ≳ ​3.5, that is, at the redshift probed by the J0660 and J0861 NBs.

2.2.1. L_Lyα and EW_Lyα detection limits

In order to estimate the minimum flux of an emission line measurable with a given NB filter, it is mandatory to know the relative contribution of this line and of the source continuum to the measured NB flux: in other words, the line equivalent width (EW; see Appendix A). In addition, the knowledge of the exact wavelength position of the line peak within the filter passband (i.e. the line redshift z) is also necessary to estimate the line flux. Indeed, the latter might be severely damped if the line is redshifted close to the filter edges, where the transmission curve is low (see e.g. Fig. 1). Finally, the line redshift is also necessary to convert the minimum measurable line flux to a luminosity value. Unfortunately, neither the EW nor the z can be directly extracted from a single NB measurement; their estimates must rely on some working assumptions. For instance, to compute the minimum Lyα luminosity () for each J-PLUS NB, we first consider that faint sources are detected with higher probability at the wavelength of the transmission curve peak. Consequently, we assume that the faintest Lyα line measurable would be redshifted at λ_obs​ = ​λ_T^max. The choice of EW, on the other hand, is more arbitrary. Although EWs as low as 5 Å have been explored in the past (e.g. Sobral et al. 2017), high-z Lyα-emitting sources typically exhibit EW > 15 − 20 Å (as in e.g. Gronwall et al. 2007; Guaita et al. 2010; Santos et al. 2020). We therefore assume that the faintest Lyα line measurable by J-PLUS would have EW_min = 25 Å (in line with e.g. Ouchi et al. 2008; Santos et al. 2016; Konno et al. 2018) and combine this value to the J-PLUS detection limits (Table 1) to compute the minimum line-flux measurable with each NB (, see Sect. 3.1 and Appendix A for details). We then convert the latter to using .

The characteristics of J-PLUS filters and its observing strategy make the resulting data sensitive to very bright Lyα emission (L_Lyα​ >  ​10^43.3 erg s⁻¹; see Table 2). We note that few studies have explored this range of L_Lyα, mostly due to the limited sky areas of their associated deep photometric surveys (see e.g. Blanc et al. 2011; Konno et al. 2016; Matthee et al. 2017b; Sobral et al. 2018a). On the contrary, J-PLUS DR1 provides multi-band imaging over ∼1000 deg², which is unprecedented for studies targeting high-z Lyα-emitting sources. The effective survey area, after masking artifacts and bright stars, sums up to ∼900 deg², which corresponds to ≳​1 Gpc³ (comoving) in each z window we sample (see Table 2). This allows us to measure the Lyα luminosity function at 2.2​ ≲ ​z​ ≲ ​3.3 and L_Lyα ≳ 2 × 10⁴³ erg s⁻¹ with high precision.

2.2.2. AGNs/QSOs or star-forming galaxies

Recent compelling findings suggest that the majority of high-z Lyα-emitting sources at L_Lyα​ >  ​2 × 10⁴³ erg s⁻¹ are AGNs/QSOs (see e.g. Nilsson et al. 2011; Konno et al. 2016; Matthee et al. 2017b; Sobral et al. 2018b,a; Calhau et al. 2020). The work of Sobral et al. (2018b), in particular, pointed out the co-existence of two different classes of luminous z​ ∼ ​2 − 3 LAEs at roughly 3 L^*, namely dust-free, highly star-forming galaxies and AGNs. In addition, a significant fraction (at least ≳20%) of bright LAEs selected by Matthee et al. (2017b) and Sobral et al. (2018a) on the Boötes and COSMOS fields (with areas of ∼0.7 deg² and ∼2 deg²), respectively, show X-Ray counterparts, which strongly points towards confirming them as AGNs/QSOs. Finally, Calhau et al. (2020) show how the fraction of AGNs/QSOs within a sample of z​ >  ​2 Lyα-emitting candidates approaches ∼100% at L_Lyα ≳ 10^43.5 erg s⁻¹.

We broadly expect the above findings to hold true over the much wider area of DR1 (bigger by a factor of ∼500), hence selecting a mixture of extremely Lyα-bright, rare star-forming galaxies (e.g. Sobral et al. 2016; Hartwig et al. 2016; Cai et al. 2017b, 2018; Shibuya et al. 2018; Marques-Chaves et al. 2019) and luminous AGNs/QSOs, numerically dominated by the latter. Indeed, our work selects objects showing strong and reliable NB excess, without employing any further criterion to disentangle its nature. Figure 2 shows typical spectra of high-z SF galaxies and QSOs, pointing out their significant diversity (see e.g. Hainline et al. 2011, for a comparison with narrow-line AGN spectra). Ideally, this difference should be mirrored by bi-modalities in the photometric properties of our selected samples, assuming that (i) both the Lyα-emitting source classes are significantly present in our selection and (ii) J-PLUS filters can effectively capture their spectral difference. For generality, we consider all the sources in our selected samples as Lyα-emitting candidates (LAE candidates). We then look for bi-modalities in their photometric properties in order to identify two distinct classes of Lyα-emitting objects. Where needed, we explicitly refer to either QSOs or SF LAEs to clearly state the distinction.

Fig. 2. Representation of our NB excess detection method. Grey lines in both panels show the observed spectra of typical z​ ∼ ​2 Lyα-emitting sources (from the publicly available VUDS DR1 spectroscopic dataset, see e.g. Le Fèvre et al. 2015; Tasca et al. 2017). Upper panel: SF LAE spectrum showing a single prominent Lyα line (here redshifted at λobs ∼ 3900 Å) and no other significant features. Bottom panel: QSO spectrum with evident CIV and CIII] lines in addition to Lyα (at λobs ∼ 4000 Å). We show the transmission curves and associated synthetic photometry of four J-PLUS bands as coloured lines and squares. From left to right: u (purple), J0395 NB (violet), g (green), and r (red). The additional lines of this caption were moved to Sect. 3.1

2.2.3. Morphology of Lyα-emitting sources in J-PLUS data

Due to resonant scatter of Lyα photons by neutral hydrogen, SF LAEs can be surrounded by faint Lyα-emitting halos and then appear more extended at Lyα wavelengths than in their continuum (e.g. Møller & Warren 1998; Fynbo et al. 2001, 2003; Nilsson et al. 2009a; Finkelstein et al. 2011; Guaita et al. 2015; Wisotzki et al. 2016; Shibuya et al. 2019, but see also Bond et al. 2010, 2012 and Feldmeier et al. 2013). As shown in Sect. 2.1, the DR1 dual-mode catalogue is based on r-band detections, which probe the rest-frame UV continuum of z​ ≳ ​2.2 sources. Observations in the UV show typical rest-frame half-light radii of about r₅₀​ ≲ ​2 kpc for z​ ≳ ​2 SF LAEs (see e.g. Venemans et al. 2005; Taniguchi et al. 2009; Bond et al. 2009, 2012; Kobayashi et al. 2016; Ribeiro et al. 2016; Paulino-Afonso et al. 2017, 2018). This translates into apparent sizes comparable to the spatial resolution of T80cam (R = 0.5″ pixel⁻¹) and to the typical J-PLUS seeing (i.e. s ≲ 1″, Cenarro et al. 2019). As QSOs are point-like by definition, we then expect both SF LAEs and QSOs at 2.2​ ≲ ​z​ ≲ ​3.3 to show compact r-band morphology. We exploit this assumption to look for low-z interlopers (Sect. 3.3.3).

Furthermore, the extended Lyα halos of SF LAEs are usually characterised by low surface brightness and can therefore be observed by means of very deep NB imaging (e.g. m_NB ≳ 26 − 27, see Leclercq et al. 2017; Bădescu et al. 2017; Erb et al. 2018) or IFU surveys (e.g. Bacon et al. 2015; Drake et al. 2017). This also applies to the peculiar class of high-z Lyα-emitting systems showing rest-frame very extended (d ≳ 20 − 30 kpc) and bright (L_Lyα​ >  ​10⁴³ erg s⁻¹) Lyα emission, namely Lyα-nebulae or blobs (i.e. LABs, see e.g. Matsuda et al. 2004; Bridge et al. 2013; Ao et al. 2015; Cai et al. 2017a; Cantalupo et al. 2019; Lusso et al. 2019). Despite extended Lyα emission usually being too faint for J-PLUS detection limits, extremely rare but sufficiently bright Lyα-emitting extended sources might still be observed within the very large area of J-PLUS DR1. These should be targeted by not relying on dual-mode catalogues but instead on analysing the 511 continuum-subtracted NB images of J-PLUS DR1 and applying specific source-extraction criteria (as in e.g. Sobral et al. 2018a). Nevertheless, we did not focus on these tasks since they merit a separate and detailed analysis which is beyond the scope of this work.

3. Lyα-emitting candidate selection

In order to select our candidates from the J-PLUS DR1, we first look for secure NB emitters (i.e. objects showing a reliable NB excess) for each of the four NBs we use. We then exploit cross-matches with external databases and the remaining J-PLUS NBs to remove low-z interlopers. Our selection rules are detailed in Sects. 3.2 and 3.3, while the following section explains how we target Lyα emission with J-PLUS NBs.

3.1. Detection of NB excess with a set of three filters

Our method to estimate the eventual NB excess for all DR1 sources and assess its significance is based on the works of Vilella-Rojo et al. (2015) and Logroño-García et al. (2019) which parallel well-established methodologies (see e.g. Venemans et al. 2005; Pascual et al. 2007; Gronwall et al. 2007; Guaita et al. 2010). We employ sets of three filters composed as: [NB; g; r], where NB stands for either J0395, J0410, J0430 or J0515. The g and r measurements are used for estimating a linear continuum (yellow dashed line in Fig. 2), which is then evaluated at the NB pivot wavelength (yellow square in Fig. 2; see Tokunaga & Vacca 2005, for the definition pivot wavelength). Finally, the ratio between the latter and the NB measurement (violet square) is used as a proxy for the Lyα line flux (see Eq. (1)). By using spectroscopically identified z​ >  ​2 QSOs, we checked that filter sets defined as [NB; u; g] provide less accurate Lyα flux measurements with respect to [NB; g; r] filter sets, due to the poor handling of the non-linear continuum in the region affected by the Lyα line profile.

As detailed in Vilella-Rojo et al. (2015), our method assumes that: (i) the emission line profile can be approximated by a Dirac-delta centred at a given wavelength λ_EL, and (ii) the source continuum is well traced by a linear function over the wavelength range covered by the three filters. The second hypothesis implies that NB measurements affected by an emission line should exhibit a photometric excess with respect to the straight line traced by g and r photometry (see Fig. 2). The goal of our method is to measure this excess and relate it to the line flux that produces it.

All the NBs we use share their probed wavelength ranges with the g filter, and therefore the eventual emission-line flux would also affect the g measurement and must be removed in order to estimate the source continuum. As detailed in Appendix A, we combine the NB, g, and r fluxes (respectively , and )³ to estimate the line-removed continuum-flux in the g and NB filters (respectively and ). In this way, we can estimate the eventual NB excess due only to an emission line as:

(1)

where the last equality is a consequence of the AB magnitude definition (Oke 1974; Oke & Gunn 1983). (m^NB) is the total NB flux (magnitude) including continuum and line contributions, while () is the continuum-only NB flux (magnitude), which is shown as a yellow square in Fig. 2. Here, Δm^NB is an indirect probe of F_Lyα, that is, the continuum-subtracted integrated line flux emitted by a given source. As fully detailed in Appendix A, by introducing the coefficients

(2)

which only depend on the transmission curve of a given filter “x” (i.e. ) and on λ_EL (i.e. the wavelength of the line peak in the NB), our method can directly estimate F_Lyα via the quantity:

(3)

We use Δm^NB for selecting reliable NB excesses (Sect. 3.2), and for computing the L_Lyα of our candidates (Sect. 4). In Eq. (3), the superscript 3FM (as in the three-filters method) points out that our method provides a photometric estimate of F_Lyα. The biases affecting are addressed in Sect. 4.1.

Figure 2 graphically explains our method when applied to both a SF LAE and a QSO spectrum⁴. In general, SF LAEs show narrow Lyα-line profiles as opposed to QSOs, whose emission can easily cover (observed) intervals of approximately a few hundred angstroms. This implies that part of the Lyα flux of a QSO can lie outside the NB wavelength coverage, and therefore might be undetected by J-PLUS NBs. The importance of this bias on depends for example on line profile details and the position of its peak in the NB. In turn, these are determined by a number of complex aspects, such as the accretion status of QSOs (e.g. Calhau et al. 2020), the transfer of Lyα photons in the hydrogen-rich ISM and IGM (e.g. Dijkstra 2017; Gurung-Lopez et al. 2019), or the metal and dust contents of the sources (e.g. Christensen et al. 2012). These details can be extracted from high-resolution spectroscopic data but not from J-PLUS photometry. For this reason, we apply Eq. (3) to all our selected candidates and then statistically correct to account for the line-flux loss, as detailed in Sect. 4.1.2.

3.2. Selection function

3.2.1. Parent sample definition

We extract our 2.2​ ≲ z ≲ ​3.3 LAE candidates from a parent sample of N​ ∼ ​1.1 × 10⁷ sources, obtained from the J-PLUS DR1 r-band selected, dual-mode catalogue (see Sect. 2.1 and Cenarro et al. 2019). Our selection targets strong NB excesses with respect to the BB-estimated continuum and removes secure contaminants (see Sect. 3.3). Its overall performance was significantly improved thanks to the spectroscopic follow-up programs described in Sect. 3.4. The selection results are presented in Sect. 3.5, while the implications of using dual-mode catalogues are addressed in Sects. 4.1 and 4.4.

Magnitude cut in g and r bands. The photometry of overly bright or faint objects is likely to be either saturated or severely affected by noise. Hence we apply a very broad cut on g and r magnitudes and their associated errors (σ_g and σ_r), namely:

We verified that these conditions do not significantly affect the final number of our candidates. Nevertheless, we account for eventual losses of continuum-faint sources with relatively bright Lyα emission (see Sect. 4.4). Spurious detections eventually included in these g and r intervals are removed by adequate S/N cuts (see below).

Detection confirmation in the three-filter set. We additionally require single-mode detection in each of the NB, g, and r bands, since all are necessary for our excess-detection method. For this, we exploit the detection flags provided by the DR1 database⁵. This condition implies that we are only sensitive to low EW at faint Lyα flux; we account for this in our completeness estimates (Sect. 4.4).

Effective exposure time cut. The normalised effective exposure time (provided in the DR1) can be used as a proxy for the number of exposures contributing to the photometry of each source. The limit excludes objects whose detection is affected by the dithering pattern of J-PLUS pointings, which might compromise the removal of cosmic rays or their extraction process.

MANGLE mask. Source photometry can be affected by optical artifacts or bright stars. J-PLUS makes use of the MANGLE software (Swanson et al. 2008) in order to mask out areas affected by these defects. For each of our selections, we apply the cumulative MANGLE mask associated to the three-filters [NB; g; r]. This reduces the total sky-coverage of our data to an effective area of A_eff ∼ 900 deg² (see Table 2 for details).

3.2.2. Pointing-by-pointing selection

The combined action of the previous cuts produces four different lists (one per NB) of N ≳ 2 × 10⁶ sources each (see Table 3). To proceed, we take into account that J-PLUS DR1 is composed of 511 different pointings (or tiles) which exhibit for example varying depths, source counts, and colours. Consequently, we apply the following conditions to each tile separately to build a selection function that is as uniform as possible.

NB excess significance. In order to select line-emitter candidates, we look for outliers in the Δm^NB versus m^NB distribution of each tile, after considering photometric uncertainties (as in e.g. Bunker et al. 1995; Fujita et al. 2003; Sobral et al. 2009; Bayliss et al. 2011; Matthee et al. 2017b, and Fig. 3). In particular, using Eq. (1) we compute the error:

(4)

Fig. 3. Example of a colour-magnitude diagram obtained for the NB filter J0410 on a DR1 pointing (out of 511). Our photometric cuts are summarised as follows: the blue dashed-dotted line shows the ΔmNB-significance threshold, while the vertical red line marks the NB S/N limit. We exclude sources below the blue dashed-dotted line and inside the grey shaded area. The orange horizontal dotted line shows ΔmNB associated to EW = 50 Å (see Eq. (7)). Grey-blue dots mark all the J-PLUS detections in the pointing, while red and purple crosses show z​ ∼ ​2.4 QSOs and low-z galaxies from SDSS DR14. Yellow triangles show J-PLUS mock data of z​ ∼ ​2.4 SF LAEs (Izquierdo-Villalba et al. 2019). Finally, our Lyα-emitting candidates are shown as green dots.

and identify reliable NB-emitters as the objects satisfying:

(5)

with Σ = 3. We account for pointing variations by anchoring our cut to the average colour ⟨Δm^NB⟩ of each tile, which acts as a rigid offset. Figure 3 shows the results of this procedure on a J-PLUS tile with ⟨Δm^NB⟩= − 0.27. As expected, only ≲10 − 15% of our parent-sample pass this cut (see Table 3).

NB signal-to-noise. We explicitly exclude objects with low-S/N NB measurements by imposing , where is the NB magnitude at which the average NB S/N of each pointing is equal to 5. This threshold has the same relative impact on the whole DR1 because only ∼35% of the sources are able to pass it. We verified that imposing ⟨S/N⟩ = 3 would lead to significantly higher contamination of our selected samples.

BB signal-to-noise. Clean BB photometry is a key requirement to estimate the NB excess of the sources. We exclude objects with g​ >  ​g_cut and r​ >  ​r_cut, where g_cut and r_cut are defined as the magnitudes at which ⟨S/N⟩ = 5 in each BB and pointing. Despite its effect on the parent samples being small (see Table 3), this cut might exclude genuine continuum-faint candidates with bright Lyα. We account for this as described in Sect. 4.4.

Minimum NB-colour cut. In principle, Lyα can be distinguished from for example CIV and CIII] of AGN/QSO spectra (e.g. Stroe et al. 2017a,b) or nebular H _β, [OIII]​_{4959 + 5007}, and [OII]​₃₇₂₇, by exploiting its generally higher intrinsic strength and EW (e.g. and Fig. 3; Vanden Berk et al. 2001; Hainline et al. 2011; Selsing et al. 2016; Nakajima et al. 2018). Therefore, we impose a NB-colour cut defined by assuming a minimum rest-frame EW for our candidates (as in e.g. Fujita et al. 2003; Gawiser et al. 2006; Gronwall et al. 2007; Hayes et al. 2010; Adams et al. 2011; Clément et al. 2012; Santos et al. 2016). Observed- and rest-frame EWs (respectively EW_obs and EW₀) are related via:

(6)

We set and obtain the corresponding from Eq. (6). We then link EW_obs and Δm^NB (defined in Eq. (1)) with the analytic expression:

(7)

(see Guaita et al. 2010, and Appendix A), where β_NB is defined in Eq. (2) and ⟨m^NB⟩ is the average colour in the tile. By requiring (orange horizontal dotted line in Fig. 3) we exclude ≳96% of our DR1 parent sample, because most sources do not show line emission. We note that the choice of has a certain degree of arbitrariness; indeed past works have explored a wide range of limiting values (see e.g. Gronwall et al. 2007; Ouchi et al. 2008; Bond et al. 2009; Nilsson et al. 2009a; Guaita et al. 2010; Konno et al. 2016; Matthee et al. 2016; Bădescu et al. 2017; Sobral et al. 2017). We fix EW₀ = 50 Å after checking our EW estimates on publicly available spectroscopic catalogues of z​ ≳ ​2 SF LAEs and QSOs (namely DR14, VUDS, and VVDS, Cassata et al. 2011; Le Fèvre et al. 2015; Pâris et al. 2018) and on the confirmed z​ ∼ ​2 QSOs in our follow-up data (see Table B.1 in Sect. 3.4). In particular, 50 Å provides a good compromise between the retrieval of z​ ≳ ​2 sources and the exclusion of z​ <  ​2 interlopers. We note that this relatively high is still close to the lower limits of EW distributions usually measured for high-z Lyα-emitting sources (e.g. Nilsson et al. 2009b; Bond et al. 2012; Amorín et al. 2017; Hashimoto et al. 2017; Santos et al. 2020). Besides, low EWs can be accessed with very narrow bands (FWHM ≲ 50 Å) and deep observations (r​ >  ​22, e.g. Sobral et al. 2017), which both act as limiting factors in our case. Finally, we stress that this condition is not directly applied on EW₀, and therefore it does not pose a strict limit on the measured EW of our candidates (see Ouchi et al. 2008, for a similar discussion).

The combined action of these cuts selects < ​1% of the DR1 parent sample for each NB (see Table 3 for the exact number). These samples are still likely to be contaminated by interlopers, such as lower-z QSOs, ELGs, and faint blue stars, which are usually targeted with BB-based colour cuts (e.g. Ross et al. 2012; Ivezić et al. 2014; Peters et al. 2015; Richards et al. 2015). We find that, in our case, these methods also significantly affect the number of selected z​ ≳ ​2 QSOs from SDSS DR14. We therefore decided to drop any colour cut because of its non-trivial effect on our selection.

3.3. Removal of residual contaminants

Despite efficiently identifying NB-emitters, the conditions in Sect. 3.2 might also select line-emitting interlopers (see Sect. 2.2). Previous works based on similar methods have usually explored limited sky regions already surveyed by deep multi-wavelength data, which supported the identification of contaminants (e.g. COSMOS, UDS, SXDS, SA22 and Boötes fields, see Warren et al. 2007; Scoville et al. 2007; Furusawa et al. 2008; Geach et al. 2008; Kim et al. 2011; Bian et al. 2012; Stroe & Sobral 2015). Unfortunately, few previous surveys uniformly cover the very wide area of J-PLUS DR1, which limits our ability to identify contaminants.

3.3.1. Cross-matches with public external databases

Interlopers with a secure identification (either spectroscopic, astrometric, or photometric) can be removed via cross-matches with public catalogues. We employ a radius of after finding that this provides a high matching completeness while keeping the number of multiple matches low for all the matched databases. Specifically, we recover 80% (95%) of all QSOs from SDSS DR14 (within the DR1 footprint) at r​ ∼ ​21.25 (r​ ∼ ​20.80) and Log (L_Lyα)​ ∼ ​44.25 (∼44.70).

SDSS DR14. We exploit the lists of spectroscopically identified galaxies (Bundy et al. 2015; Hutchinson et al. 2016), stars (Majewski et al. 2017), and QSOs (Pâris et al. 2018) provided by the recent SDSS-IV DR14 (DR14 hereafter, Blanton et al. 2017; Abolfathi et al. 2018). Given the wide overlap with J-PLUS DR1 and the higher depth of DR14 (Cenarro et al. 2019), this cross-match ensures the removal of secure contaminants from our selection. As discussed in Sect. 2.2, QSOs can act as both interlopers and genuine candidates depending on their z, and therefore we need to rely on a list of securely identified QSOs. The Pâris et al. (2018) catalogue includes N ≳ 5.3 × 10⁵ sources observed by BOSS and eBOSS surveys (Dawson et al. 2013, 2016) and confirmed as QSOs by careful inspection. We keep genuine Lyα-emitting sources at the z sampled by each NB, while the rest are identified as contaminants and removed. The cross-match with DR14 shows a generally low contamination (Table 4), with low-z galaxies accounting for 5.1%, 4.3%, 5.3%, and 3.1% of our J0395, J0410, J0430, and J0515 NB samples, respectively. On the other hand, the z​ ≲ ​2 QSO fraction drops from 11.1% to 0.3%, mirroring the drop of DR14z​ ≳ ​2.2 QSOs. Finally, SDSS stars account for ≲2% of our samples. These fractions are likely to be underestimated given the different depth of the two surveys and eventual mismatches between DR14 and DR1 catalogues. Nevertheless, being measured on spectroscopically confirmed sources, these are secure contamination estimates.

Gaia DR2. Our spectroscopic follow-up program 2018A (see Sect. 3.4) showed a non-negligible contamination from stars in our samples. To limit this issue, we built a specific criterion for excluding stars based on the very accurate measurements offered by Gaia DR2 data (Gaia Collaboration 2018). We identify secure stars using the significance of proper-motion assessments. Specifically, we exclude the J-PLUS sources with a counterpart in Gaia DR2, showing significant measurements (σ​ >  ​3) in each proper motion component, i.e.

(8)

where σ_pmra, σ_pmdec, and σ_μ are the errors on proper motion (RA and Dec) and parallax, respectively. With this cut, we explicitly remove objects showing significant apparent motion from our list of LAE candidates. The good performance of this criterion was confirmed by the results of our second follow-up program, whose targets were selected from the results of our updated pipeline (see Sect. 3.4 for details). The contamination from Gaia DR2 is presented in Table 4.

GALEX-UV. Lyα-emitting sources at z​ >  ​2 are generally expected to appear faint at (observed) UV wavelengths due to the dimming action of the Lyα-break and Lyman-break (e.g. Steidel & Hamilton 1992; Steidel et al. 1996, 1999; Shapley et al. 2003). On the contrary, z​ <  ​2 AGN/QSOs, blue stars, and low-z star-forming galaxies can show significant UV emission. We exploit this property in order to remove z​ <  ​2 interlopers by cross-matching our catalogues with GALEX all-sky UV observations (Gil de Paz et al. 2009). In particular, we remove sources with a S/N​ >  ​3 detection in either of the two FUV and NUV GALEX bands (see e.g. Ciardullo et al. 2012, and Table 4). In order to verify the validity of our assumption according to which only z​ <  ​2 sources are expected to be significantly observed in UV, we additionally matched the J-PLUS sources with counterparts in GALEX to the spectroscopic sample of DR14. This analysis confirmed that > 99.5% of sources with UV-bright GALEX detection show a spectroscopic z​ <  ​2, and therefore act as contaminants in our selection.

LQAC-3. The third release of the Large Quasar Astrometric Catalog (Souchay et al. 2015a,b) is a complete archive of spectroscopically identified QSOs. By combining data from available catalogues, it provides the largest complement to the DR14 list (Pâris et al. 2018). We exclude sources included in LQAC-3 with spectroscopic z lying outside the range probed by each NB. As expected, this step identifies only a few additional interlopers (see Table 4).

3.3.2. Multiple NB excesses

We target additional interlopers by exploiting the whole set of J-PLUS NBs. Indeed, we expect SF LAEs to not show any additional NB feature (e.g. Shapley et al. 2003; Nakajima et al. 2018), while QSOs at the targeted z can exhibit only particular combinations of NB excesses. Consequently, we remove candidates with multiple NB excesses incompatible with z​ >  ​2 spectral features (e.g. Matthee et al. 2017b). On the other hand, excesses consistent with z​ >  ​2 can also be produced by low-z interlopers. As an example, Fig. 4 shows the photo-spectra of a z​ ∼ ​0.05 galaxy (upper panel) and a z​ ∼ ​2.25 QSO (bottom panel) from the DR14 spectroscopic samples. Both sources show simultaneous excesses in J0395 and J0515 filters (purple and yellow empty squares, respectively) with respect to the linear continuum traced by g and r BBs (yellow dashed line). Furthermore, both photo-spectra exhibit comparable BB colours and might therefore be confused by our selection. Since J-PLUS data do not allow us to distinguish between these objects and measure this source of contamination, we estimate a statistical correction as explained in Sect. 4.3.

Fig. 4. Examples of multiple NB excess in J-PLUS photo spectra. Empty and filled squares mark NB and BB photometry, respectively, while the dashed yellow line shows the linear continuum we estimate through g and r BBs (in green and red, respectively). This comparison shows that both a z​ ∼ ​0.05 galaxy (upper panel) and a z​ ∼ ​2.25 QSO (lower panel) exhibit significant excesses in J0395 and J0515 NBs (second and fifth empty squares from the left, respectively). The J0515 excess is produced by H β at z​ ∼ ​0.05 and the CIV line at z​ ∼ ​2.25, but its nature is hardly distinguishable using J-PLUS photometry.

3.3.3. Morphological cut

We expect our z​ ≳ ​2.2 candidates to appear compact in J-PLUS data (see Sect. 2.2.3), and therefore sources showing extended morphology are likely to be low-z interlopers. The DR1 catalogue provides a parameter 𝒞 which allows us to discriminate between compact (𝒞​ ∼ ​1) and extended objects (𝒞​ ∼ ​0, see López-Sanjuan et al. 2019b, for details). By cross-matching the whole DR1 sample to SDSS spectroscopic catalogues of galaxies and QSOs, we checked that more than ≳90% of galaxies in SDSS (z​ ≲ ​1, see Hutchinson et al. 2016) and only ≲5% of DR14 QSOs (at any z) are found at 𝒞 ≤ 0.1. We then remove objects with 𝒞 ≤ 0.1 from our selection (see penultimate column to the right in Table 4).

3.4. Spectroscopic follow-up at the GTC telescope

This section presents two spectroscopic follow-up programs executed at the Gran Telescopio Canarias (GTC) telescope⁶ in the semesters 2018A and 2019A. The spectroscopic confirmation of a subsample of our candidates allowed us to assess the performance of our selection, to refine our methodology, and to estimate its residual contamination. Overall, these programs confirmed 45 sources selected among our J0395 NB emitters.

3.4.1. Program descriptions

To ensure uniform observations and comparable results, we performed the same target selection and required identical observing conditions for both programs (namely GTC2018A and GTC2019A). In particular, we randomly selected a sample of 24 (21) Lyα-bright candidates (L_Lyα​ >  ​10^43.5 erg s⁻¹) for program GTC2018A (GTC2019A), spanning the entire luminosity range covered by our candidates. We stress that targets for GTC2019A were selected after refining our selection with the help of GTC2018A results. We requested the use of the OSIRIS spectrograph and the R500B grism in order to exploit its good spectral resolution (R​ ∼ ​500, which translates to for the 0.8″ slit width we requested). The exposure times for our targets were computed assuming the observing conditions summarised in the header of Table B.1. These were calibrated to achieve S/N​ ≥ ​3 (in each λ bin) over the whole OSIRIS spectral range in order to identify eventual emission lines and measure their integrated flux.

We limited our programs length to < 20 h to ensure their completion. Due to the long observing times required by our targets, we only followed up candidates selected by J0395 NB. The target selection balanced the total observing time and the uniform sampling of our candidate L_Lyα distributions. Finally, we excluded objects with previous spectroscopic identifications (at any z). Our proposals were respectively awarded with 11.56 and 18.95 hours of observations and were both fully executed.

3.4.2. Spectroscopic results

The results of both programs are shown in Table B.1. Overall, we identified 29 of our 45 targets (64.4%) as genuine z​ ∼ ​2.2 Lyα-emitting sources, 8 of them (17.7%) as z​ ∼ ​1.5 QSOs emitting CIV at λ_obs ∼ 4000 Å, one (2.2%) Lyβ-emitting QSO at z​ ∼ ​2.76, five (11.3%) blue stars, and two of the 45 targets (4.4%) as low-z galaxies selected because of their narrow emission lines. As an example, Fig. 5 shows a spectrum for each different source class together with its associated J-PLUS photometry. Both z​ ∼ ​2.2 and z​ ∼ ​1.5 QSOs show prominent line emission at λ_obs ∼ 3950 and are consequently selected as genuine J0395 NB emitters. The same applies to the z​ ∼ ​2.8 QSO emitting Lyβ at λ_obs​ ∼ ​3950. On the contrary, the remaining sources do not show significant spectral features, indeed their selection is due to strong blue colours combined to a barely-significant NB excess (see e.g. third panel from above). In particular, the star and galaxy interlopers (i.e. third and fourth panels from the top) were picked as targets before we refined our selection rules and the J-PLUS DR1 was re-calibrated (López-Sanjuan et al. 2019a). With the current J-PLUS photometry and our updated selection these objects are not re-selected (right column of Table 5). Given the absence of emission lines at λ_obs ∼ 3900 Å for these objects, their low-significance NB excess is likely due to imperfections in their photometry. In the case of the z​​ ∼ ​​0.5 galaxy (fourth panel from the top in Fig. 5), we additionally observe a discrepancy between the spectrum and J-PLUS data. A number of possible explanations can account for this, such as errors in the spectrum extraction and calibration, overly low spectroscopic S/N at λ_obs ≲ 4500 Å, or artifacts biasing only the J0395 photometry. On the contrary, the excess of the z​ ∼ ​2.8 QSO (bottom panel in Fig. 5) is due to the Lyβ line redshifted at λ_obs ∼ 3950 in the observed spectrum, although this is in tension with J-PLUS photometry. In this case, QSO variability might play a role (e.g. Hook et al. 1994; Kozłowski 2016), as might photometric imperfections.

Fig. 5. Summary of our spectroscopic results showing one spectrum (grey line in each panel) for each source class identified in our target lists. From top to bottom: z​ ∼ ​2.2 QSO, z​ ∼ ​1.5 QSO, star, z​ ∼ ​0.5 galaxy, and Lyβ-emitting QSO. The corresponding J-PLUS photometry is shown as coloured squares. The star and galaxy targets show low-significance excesses in J0395 NB (third square from the left). Indeed, these interlopers were selected as targets by the first version of our methodology, i.e. before applying the improvements due to CTG2018A results and the re-calibration of J-PLUS data (López-Sanjuan et al. 2019a).

Overall, 40 of our 45 targets (88.9%) are genuine line emitters, confirming the efficiency of our selection. Moreover, the stars contamination is reduced from ∼17% to ≲5% between the two programs (see Table 5). Indeed, guided by the GTC2018A results, we (i) excluded sources with significant apparent motion according to Gaia DR2 and (ii) selected EW₀ = 50 Å as our limiting value for defining the Δm^NB cut (see Sect. 3.2). Our improved methodology retrieves 15 of the 24 original GTC2018A targets, with 11 of those 15 (∼74%) being z​ ∼ ​2.2 QSOs and 4 (∼26%) being z​ ∼ ​1.5 QSOs (i.e. no star is re-selected). Nevertheless, four out of 15 z​ ∼ ​2.2 QSOs are not re-identified as line emitters. This is partly due to the new calibration of the entire J-PLUS data (occurred after; GTC2018ALópez-Sanjuan et al. 2019a); indeed two out of the four z​ ∼ ​2.2 QSOs which are not re-selected fail to pass the NB S/N criterion due to their recomputed NB photometry. Finally, the fraction of genuine NB emitters significantly increased from ∼​74% for GTC2018A to over 95% for GTC2019A thanks to our methodology improvements.

3.5. Selected samples of Lyα-emitting candidates

The final results of our selection procedure are four samples of z​ >  ​2 Lyα-emitting candidates which meet all the following requirements: (i) reliable excess in the NB used for their selection, (ii) secure detection and photometry in the filter triplet [NB; g; r], (iii) no spectroscopic counterparts in DR14 with redshift outside the ranges probed by each NB, (iv) no apparent motion according to Gaia DR2 data, (v) no significant observed-frame UV detection in GALEX, (vi) compact morphology and, eventually, (vii) multiple NB excesses compatible with being z​ ≳ ​2 sources. These lists account for 2547, 5556, 4994, and 1467 sources for J0395, J0410, J0430, and J0515 NBs, respectively (see Table 4), which translates into approximately 2.8, 6.2, 5.6, and 1.5 objects per square degree, respectively. We underline that these samples are the largest collections of Lyα-emitting candidates to date within the narrow redshift bins that we can access (see e.g. Guaita et al. 2010; Cassata et al. 2015; Konno et al. 2016; Matthee et al. 2017b; Sobral et al. 2018a).

The drop in the number of counts for J0515 NB can be ascribed to the combination of J-PLUS data depth and the cosmological decrease of bright SF LAEs and AGN/QSO number densities at z​ ≳ ​2.5 (e.g. Nilsson et al. 2009a; Ciardullo et al. 2012; Sobral et al. 2018a). Indeed, the right panels of Fig. 6 show that J0515 NB can only access ranges of Log (L_Lyα) and Lyα Log (EW₀) which are significantly higher than the other NBs. In general, filters sampling smaller wavelengths can access fainter Lyα luminosity and smaller EW₀, as a result of the combination between J-PLUS depth and the probed z interval.

Fig. 6. Left figure: EW0 distribution of our selected candidates as a function of r and g​ − ​r colour (left and right panels, respectively). Squared points and error bars show respectively the distribution median and 16th–84th percentiles in each magnitude and colour bin. Points have been artificially shifted for better visualisation. The values for J0395 filter at r​ <  ​20 (g − r​ <  ​0.75) are systematically below the theoretical cut we apply (Sect. 3.2). This is due to the small overlap between this NB and the g BB, which is reflected in a poor extrapolation of the linear continuum at the NB filter pivot wavelength (see e.g. Ouchi et al. 2008, and the discussion in Sect. 3.2). Right figure: normalised distributions of our candidates in EW and Lyα luminosity for each filter. This result clearly shows that filters sampling higher redshifts also sample brighter Lyα luminosity. This is a direct effect of J-PLUS detection limits which only allow to observe brighter and rare objects at higher redshifts. We address this issue by applying the completeness corrections described in Sect. 4.4.

3.5.1. EW₀ and L_Lyα distributions

The left panels in Fig. 6 shows the distribution of EW₀ measured in our samples as a function of both r magnitude and g​ − ​r colour. As mentioned in Sect. 3.2, our cut on Δm^NB derives from a theoretical expected limit of EW₀ = 50 Å. Nevertheless, not all selected sources display EW₀ > 50 Å (see also Ouchi et al. 2008, for a similar discussion). This is evident for the J0395 NB candidates, whose EW₀ distribution is systematically below 50 Å for r ≲ 20 and g − r ≲ 0.75. Indeed, the small overlap between J0395 and g transmission curves ultimately provides a relatively poor extrapolation of the linear continuum up to the pivot wavelength of J0395, which translates into an under-estimation of EW₀. This induces a bias on our Lyα luminosity measurement, which we account for as described in Sect. 4.1.2. On the other hand, no significant or systematic bias affects EW₀ with respect to colour, as shown by the flat g​ − ​r distribution in Fig. 6. We confirmed this by using the spectra of z​ ∼ ​2 DR14 QSOs, but we do not show the results for the sake of brevity.

Overall, our distributions are broadly consistent with previous determinations of the rest-frame EW of z​ ∼ ​2 − 3 LAEs (either SF LAEs or AGN/QSOs, see e.g. Gronwall et al. 2007; Guaita et al. 2010; Hainline et al. 2011; Ciardullo et al. 2012; Shibuya et al. 2014; Hashimoto et al. 2017; Santos et al. 2020). Interestingly, our samples include a moderate fraction of sources (≲7%, on average) showing EW₀ > 240 Å (e.g. Ouchi et al. 2008; Santos et al. 2020). High-EW LAEs have been studied with particular interest (e.g. Cantalupo et al. 2012; Kashikawa et al. 2012; Shibuya et al. 2014) since nebular emission of Pop-II stellar populations can only account for (e.g. Charlot & Fall 1993; Hernán-Caballero et al. 2017). At the same time, high Lyα EWs can easily be produced by AGNs/QSOs, which are likely to dominate our selected samples. Since analysing high-EW LAEs would require a careful separate analysis, we refrain from commenting further on this topic. Nevertheless, we underline that our lists of selected candidates can provide catalogues of high-EW LAE targets for upcoming studies.

3.5.2. Relative abundance of QSOs and SF LAEs

The design of J-PLUS filters potentially allows us to simultaneously capture a peculiar combination of high-z lines with different NBs (see Sect. 3.3.2). For instance, QSOs emitting Lyα at z​ ∼ ​2.3 could show simultaneous NB excesses in J0410 and J0515 NBs (the latter being due to CIV emission). This offers the possibility to investigate the relative fraction of AGNs/QSOs and SF LAEs in our samples, because the latter should not exhibit such double-NB emission. We therefore separate the DR14 QSOs selected with J0410 from the rest of J0410 candidates and plot the colour distribution of these two source classes.

Figure 7 shows the colour space defined by J0410 and J0515 NBs with respect to g BB. In this plane, both SF LAEs and QSOs should exhibit g​ − ​J0410​ >  ​0 due to Lyα emission, but the CIV line (usually much stronger in QSOs than in SF LAEs, see e.g. Hainline et al. 2011; Stroe et al. 2017a,b; Nakajima et al. 2018) should displace QSOs at g​ − ​J0515​ >  ​0. Nevertheless, our colour distribution (blue solid histogram) does not show any evident bimodality and no significant overdensity at g​ − ​J0515​ ∼ ​0, as expected for SF LAEs. This suggests that either (i) our J0410 candidates are mostly dominated by z​ ∼ ​2 QSOs or that (ii) the J-PLUS filter set does not reliably disentangle the different spectral features of high-z SF LAEs and AGNs/QSOs.

Fig. 7. Comparison between the colour–colour distributions of DR14 QSOs (green squares) within our J0410 sample and of the remaining J0410 genuine candidates (blue dots, after removing known interlopers; Sect. 3.3). The two source classes occupy comparable colour regions, suggesting that our selection results might be effectively dominated by z​ ∼ ​2.3 AGN/QSOs. This scenario is also supported by the results of our spectroscopic program (Sect. 3.4).

We further analyse the nature of our candidates by exploiting the cross-match with the all-sky WISE data (Wright et al. 2010) provided by the J-PLUS DR1 database. In particular, we compare the r − W1 colour of our candidates with WISE counterparts to the synthetic-photometry colour-tracks of galaxy and QSOs templates (from Polletta et al. 2007, and Hernán-Caballero et al. 2016, respectively). Figure 8 shows how the colours of our candidates are clearly compatible with those of QSOs, while being significantly different from the galaxy ones. We also highlight the SDSS QSOs (black-contoured dots) and the confirmed QSOs of our GTC programs (yellow-contoured dots, only for J0395 NB) to underline the comparability of the properties of our our candidates with those of spectroscopically confirmed QSOs.

Fig. 8. Coloured points: r − W1 colour vs. redshift for all our candidates with a counterpart in WISE. Our sources are uniformly spread within the z interval sampled by each NB. Grey lines: tracks for different galaxy templates (from the SWIRE library, Polletta et al. 2007) and the QSO template of Hernán-Caballero et al. (2016). Black-contoured and yellow-contoured points mark the SDSS QSOs in our selection and the GTC targets (only for J0395 NB), respectively. Our candidates are all compatible with the high-z QSO template, suggesting that the fraction of SF LAEs in our sample is very low. In addition, this suggests that we identify a large number of sources without previous spectroscopic identification as high-z QSOs.

Interestingly, by joining these various pieces of evidence with the results of our spectroscopic follow-up programs, we expect our samples of candidates to be dominated by z​ >  ​2 QSOs which still lack spectroscopic identification. Indeed, by considering the number of our genuine candidates without SDSS identification (namely 2057, 4959, 4494, and 1377 for J0395, J0410, J0430, and J0515 NBs, respectively) and conservatively applying a residual contamination of ≳35% (as suggested by our GTC follow-up) our method identifies ≳1300, 3200, 2900, and ∼900 z​ >  ​2 sources as z​ >  ​2 QSOs for these respective NBs, for the first time. This is shown in both Figs. 7 and 8 by the wide difference between the number counts of DR14 QSOs within our selection and our remaining genuine candidates. We interpret this as an effect of the NB-based selection we perform, which efficiently targets the line-emission features of these objects, eventually missed by previous target-selections based on BB colours (e.g. Richards et al. 2009; Ross et al. 2012; Ivezić et al. 2014). Nevertheless, a systematic and uniform spectroscopic confirmation of our samples is needed to validate these findings.

4. Construction of the Lyα luminosity function

The luminosity function Φ(L) of a given class of sources is usually defined as their comoving number density per unit luminosity (see e.g. Schmidt 1968). Following a common convention in the literature, we express our LFs in logarithmic units of luminosity and hence use the following definition:

(9)

where the sum at the numerator is extended to all the objects in a given bin of (logarithmic) luminosity Δ Log(L_Lyα), while the coefficients P_i and C_i are statistical weights that account respectively for the sample purity and completeness (as detailed below). We exploit our lists of candidates selected with J0395, J0410, J0430, and J0515 NBs to build four determinations of the Lyα LF at the redshifts given by Table 2. The following sections detail the steps we perform to assess the reliability of our Lyα flux measurements (Sect. 4.1), measuring the Lyα luminosity of our candidates and the cosmological volume probed by J-PLUS NBs (Sect. 4.2) and estimating the purity (Sect. 4.3) and completeness (Sect. 4.4) of our selection.

4.1. Retrieval of the total Lyα flux

Our Lyα flux measurements ( hereafter) can be affected by systematic uncertainties associated with both the J-PLUS aperture photometry and our measuring method (see Eq. (3) and (A.13). In order to build our Lyα LFs, we first study the differences between and a corresponding spectroscopic measurement (i.e. F_Lyα), assuming that the latter provides a reliable estimate of the total emitted Lyα flux of the sources. We then compute statistical corrections which account for the bias between and its spectroscopic analogue using spectroscopically identified QSOs (Pâris et al. 2018) at the redshift sampled by each NB and their counterparts in the DR1 catalogue. We obtain the Lyα flux from QSOs spectra with the methodology described in Appendix B.

4.1.1. Aperture correction

To make sure that auto-aperture photometry (see Sect. 2.1) does not introduce any bias on the photometry of our candidates, we compare the synthetic flux of SDSS QSOs to the analogous measurements in J-PLUS DR1. For the sake of brevity, here we summarise our findings and present the details of this check in Appendix C. In general, no significant bias (≲0.2σ_r) affects the auto-aperture flux of point-like sources for each NB. This is expected because SDSS fiber-based photometry is normalised to PSF apertures⁷ and J-PLUS auto-apertures generally encompass the whole PSF of point-like objects, making the two measurements effectively comparable. Consequently, we do not apply aperture corrections to . On the other hand, the flux comparison points out the need for an additional statistic uncertainty on top of the J-PLUS photometric errors (see Appendix C for details). We then re-scale the r uncertainties and propagate them on . Finally, we account for the latter on our LF determinations as discussed in Sect. 4.4.4.

4.1.2. Filter width correction

A fraction of the flux of broad lines (i.e. broader than the FWHM of the measuring NB) can be systematically lost by photometric measurements, especially if the line-peak is displaced at the edge of the NB transmission curve. Star-forming LAEs usually show a narrow Lyα as opposed to the customary broad-line profile of QSOs (e.g. Vanden Berk et al. 2001; Telfer et al. 2002; Selsing et al. 2016, and Fig. 2). For these reasons, we expect this bias to significantly affect the measurements of QSOs, while not influencing those of SF LAEs. At the same time, no SF LAEs were observed among our followed-up targets (Sect. 3.4.2), in line with previous results suggesting that AGNs/QSOs dominate the samples of photometrically selected LAEs at L_Lyα ≳ 2 × 10⁴³ erg s⁻¹ (see e.g. Santos et al. 2004; Konno et al. 2016; Matthee et al. 2017b; Sobral et al. 2018a; Calhau et al. 2020). Furthermore, the (expected) low fraction of SF LAEs in our final selection cannot be reliably disentangled from QSOs by J-PLUS photometry (Sect. 3.5.2). This hinders the possibility of applying a flux correction exclusively to a subclass of our candidates. Consequently, we consider our method for measuring as valid and then apply a statistical correction to all our candidates. In particular, we obtain the corrected Lyα flux as follows (see Appendix C):

(10)

The quantity ΔF is a rigid offset obtained from the normalised distribution of flux difference: , where is our Lyα flux estimate and is its spectroscopic analogue measured on SDSS QSOs (see Appendix C for details). We obtain a ΔF for each NB and then use the corrected values for our luminosity function computation. With this analysis we also obtain a correction for the error on , which we propagate on our Lyα LF determination (see Sect. 4.4.4).

Finally, the obtained with J0378 NB are affected by a significant bias (ΔF = 1.75 ± 0.35). This can be ascribed to the wavelength separation between this NB and g, which leads to a poor extrapolation of the linear continuum approximation. Consequently, we exclude J0378 from the list of NBs we use.

4.2. Computation of L_Lyα and cosmological volume

We compute the Lyα luminosity as:

(11)

where d_L(z)​ = ​d_c(z)​ ⋅ ​(1 + z) and d_c(z) are the luminosity and comoving distances of our sources which are computed by assuming PLANCK2015 cosmology (Planck Collaboration I 2016; Planck Collaboration XIII 2016).

In order to compute the d_L(z) of our candidates without spectroscopic determination it is necessary to assume a value for z. Having no knowledge of their nature, we use the z_p obtained by shifting the Lyα rest-frame wavelength to the pivot wavelength (Tokunaga & Vacca 2005) of the NB used for selection (see Table 2). Consequently, the uncertainty σ_z is obtained from the half width of each NB (see Table 2). This does not apply to the candidates with a spectroscopic counterpart, as in these cases we use the DR14z and σ_z. Finally, we propagate the redshift errors on the total L_Lyα uncertainty and on our LF determinations (see Sect. 4.4.4).

The cosmological volume sampled by our data depends on the z windows associated to Lyα detection and the DR1 area not affected by masking, for each NB. In our case, the redshift intervals are given by the FWHM of each NB (see Table 2) and are converted to cosmological volumes by assuming the PLANCK2015 cosmology (Planck Collaboration I 2016; Planck Collaboration XIII 2016). On the other hand, the effective area observed in a given band can be obtained with the MANGLE software (Hamilton & Tegmark 2004; Swanson et al. 2008). As we require single detection in [NB; g and r], we computed the intersection between the three associated MANGLE masks (see Table 2). We assume negligible errors on volume estimates for our LF computation.

4.3. Estimate of the sample contamination

The steps detailed in Sect. 3.3 do not ensure the identification of all contaminants, as confirmed by our follow-up results (Sect. 3.4.2). For this reason, we estimate the residual contamination of our samples by computing a statistical purity weight for our candidates as a function of their r-magnitude:

(12)

Here, N_interlopers(r) and N_total(r) are the number of secure interlopers (see Sect. 3.3) and the total number of candidates at a given magnitude, respectively. We then fit P(r) with an error function and use the latter to obtain the statistical weight of each genuine candidate to the final LF. Figure 9 shows the error-function fits for each NB (solid coloured lines) and the computed P(r) values for the J0430 filter as an example (dotted grey line). The purple empty square shows the average purity measured on the complete sample of both our spectroscopic follow-up programs. This is in good agreement with the statistical weights of each NB (i.e. P(r)​ ≳ ​60% at r ≳ 18.5). The high values reached by J0515 (∼80%) are driven by the drop of interlopers with spectroscopic identification at z​ ≳ ​3.

Fig. 9. Statistical purity weight for each NB (coloured solid lines) as estimated by fitting an error-function to the computed purity. The grey-dotted line shows the computed purity of J0430 NB as an example. All filters show similar purity weights, rising to ≳60% at r ≳ 18.5. This is in agreement with the average purity of our 45 spectroscopic targets (Sect. 3.4.2) which is shown as a purple empty square.

4.4. Estimate of the sample completeness

Genuine line-emitting candidates might be lost by our selection because of the J-PLUS detection limits and source extraction, the effect of photometric errors, and the r-band pre-selection of our parent samples (see e.g. Geller et al. 2012; Loveday et al. 2012; Gunawardhana et al. 2013). In order to correct for these known issues, we estimate the completeness⁸ of our samples by considering three different components. Specifically, we account for (i) the DR1 source-extraction process (i.e. detection weightC^d), (ii) our selection methodology (i.e. selection weightC^s), and (iii) r band pre-selection of dual-mode catalogues (i.e. dual-mode weightC^dm). We obtain the total completeness weight of each candidate to the final LF as: .

4.4.1. Detection completeness

The detection completeness of each J-PLUS pointing (for each filter) is automatically computed by the standard source-extraction pipeline as:

(13)

where k_s and r_s are computed for each pointing; these are the decay-rate of and the magnitude at which reaches 50%, respectively. All details of this computation are provided by the J-PLUS DR1 datababse. We obtain from the [k_s; r_s] parameters and r corresponding to each DR1 source.

4.4.2. Selection completeness

Starting from r-detected catalogues, our selection makes use of NB-excess significance and a linear estimate of the continuum slopes of the sources, which are related to their g​ − ​r colour (see Sects. 3.1, 3.2 and Figs. 2 and 3). In order to capture its multiple dependencies, we test the retrieval efficiency of our selection as a function of r magnitude, Lyα flux, and g​ − ​r colour. In particular, we compute the recovery rate of simulated candidates over wide ranges of these three quantities by re-applying each of our selection rules. This accounts for source loss at different Lyα flux, continuum, and EW. We organise the measured recovery rates in a 3D grid which we interpolate at the measured position of each genuine candidate to compute its selection weight . The details of this computation are given in Appendix D, while Fig. 10 shows an example of the estimated selection completeness in a bin of r magnitude and g​ − ​r colour for the J0430 NB.

Fig. 10. Example of the recovery fraction of our selection as a function of Lyα flux, computed for J0430 NB in a bin of r magnitude and g​ − ​r colour (namely at r = 20.1 and g​ − ​r = 0.2). The full 3D grid is shown in Appendix D for the same NB. All filters show comparable values of recovery fractions, and therefore we only report the case of J0430, for brevity.

4.4.3. r−−L_Lyα bivariate completeness

The use of r-band detected catalogues makes our selection prone to the loss of continuum-faint z​ ≳ ​2 Lyα-emitting sources, with non-trivial effects on the EW distribution of our selected samples. At low Lyα flux, for example, the r-detection requirement might favour the selection of high-EW Lyα-emitting sources. This issue has already been pointed out by authors in the literature whose selection function was built on the convolution of r-band detection and NB-excess significance. In particular, Gunawardhana et al. (2015) showed that accounting for this effect requires a multi-variate approach. In other words, the fraction of undetected continuum-faint line-emitters can be estimated by modelling the full 2D luminosity function of candidates in the r vs. line-luminosity plane.

We closely follow the methods of Gunawardhana et al. (2015) applying their computations to the r versus L_Lyα space. The details of this procedure and its main equations are presented in Appendix D. Briefly, we assume that the 2D LF can be modelled by the product of two functions, which respectively describe the r and Log(L_Lyα) distributions (see also Corbelli et al. 1991). We combine a Schechter (in logaritmic form) and a Gaussian (in Log L_Lyα) function (as in Gunawardhana et al. 2015, see Appendix D). By fitting this 2D model to our measured 2D LF, we can model the number density of sources in regions of the r versus L_Lyα plane affected by our incompleteness. Finally, the ratio of our data to the 2D model (in the 2D space r vs. L_Lyα) allows us to compute the statistical weight C^dm(r, Lyα) for each source, which accounts for the loss of r-faint Lyα-emitting sources. Figure 11 shows the results of our 2D modelling for the J0430 filter. In particular, the top and right panels show the projection of both our 2D LF (green solid lines) and 2D model (red dashed line), respectively, along the Log(L_Lyα) and r axes. It is clear how the model extrapolates our measurements at r​ >  ​19.5 and Log(Lyα/erg s⁻¹)​ <  ​44.2.

Fig. 11. Central panel: full 2D luminosity function of our J0430 candidates as a function of r and Log(LLyα). Green solid lines in the top and right panels show the projections of the 2D LF along the Log(LLyα) and r axes, respectively. The red dashed lines show the projection of the 2D model along the same axes. This model was fitted to the 2D distribution shown in the central panel (see Appendix D for computational details) and allows us to extrapolate our data distribution at faint r and Lyα luminosity. We use the data-to-model ratio (in the r vs. Log(LLyα) 2D plane) to compute the weight of each candidate.

4.4.4. Errors on the Lyα luminosity function

The uncertainties on the redshift and Lyα flux of our sources, the binning in Lyα luminosity, and the internal variance of the samples (due to differences between each J-PLUS pointing) jointly contribute to the errors on our final LFs (e.g. Sobral et al. 2018b). We measure each source of uncertainty separately and finally sum in quadrature their different contributions. To account for L_Lyα uncertainties, we repeat the determination of our LF 1000 times by perturbing the sources flux each time according to its uncertainty. During this procedure, we keep the redshift of the sources fixed to z_p (see Sect. 4.2) in order to evaluate only the contribution of flux perturbations to the final errors on our LFs. We then compute the asymmetric errors from the percentiles of the LFs distribution as σ⁻ = 50th − 16th and σ⁺ = 84th − 50th, where 84th, 50th, and 16th are the corresponding distribution percentiles. The contribution of redshift errors is accounted for in an analogous way by fixing the flux measurements. To account for the internal variance of our LAE candidate sample due to field variations in J-PLUS DR1, we perform random realisations of the luminosity function by splitting our samples into ten independent subsamples and computing a LF for each subsample. We repeat this process 1000 times and ultimately extract the errors from the 16th and 84th percentiles of the LFs distribution (see above). Finally, we also add the Poissonian errors () associated to the number counts of the sources in each bin to the total LF uncertainties.

5. Results

In this section we present the four Lyα LFs that we compute. In particular, we compare our measurements to previous results in the literature in Sect. 5.2, we describe the computation of the Schechter parameters of the Lyα LFs in Sect. 5.3, and finally we estimate the fraction of AGNs/QSOs as a function of luminosity in Sect. 5.4.

5.1. The Lyα luminosity functions at 2 < z< 3.3

Figure 12 shows the four determinations of the Lyα LF we compute at z​ ∼ ​2.25,  2.37,  2.54, and z​ ∼ ​3.24 (coloured empty squares). For each NB, we only consider the candidates with a total completeness weight of C = C^d ⋅ C^s ⋅ C^dm​ >  ​0.85 (see Sect. 4.4 and Appendix D). This excludes sources whose contribution is severely affected by the completeness correction, especially at Log(L_Lyα/ erg s⁻¹)​​ ≲ ​​44. Overall, our results probe a luminosity interval of ∼​​1.5 dex, from Log(L_Lyα/ erg s⁻¹)​ ∼ ​44 to Log(L_Lyα/ erg s⁻¹)​ ∼ ​45.5. These regimes are expected to be significantly populated by Lyα-emitting AGNs/QSOs (e.g. Borisova et al. 2016; Matthee et al. 2017b; Sobral et al. 2018b; Calhau et al. 2020). Interestingly, our results extend by ∼​​1 dex into a previously unconstrained Lyα luminosity range, allowing us to probe it with high precision. In addition, our data extend down to ∼​​10⁻⁸ Mpc⁻³, a limit which is barely touched upon by previous studies (see e.g. Sobral et al. 2018a). These remarkable features are ultimately attained because of the very wide area covered by J-PLUS NB imaging (unprecedented for Lyα LF determinations), which balances the J-PLUS depth (r​ <  ​22).

Fig. 12. Lyα luminosity functions for each of the NB filters we used in our study (coloured squares). The grey shaded areas show the Lyα luminosity limit (vertical limit) and the limiting number density measurable by J-PLUS (horizontal limit). The wide area explored by the NBs of the J-PLUS survey allow us to remarkably extend the range of luminosity sampled by previous studies (coloured circles, triangles, hexagons, and diamonds in each plot) and to explore previously unconstrained LLyα intervals. Dashed lines mark the best-fit determinations by Sobral et al. (2018a), split respectively into a Schechter (light green) and a power-law contribution (dark green). Our results provide tight constraints at 44.5 ≲ Log (LLyα/erg s−1)≲45.5, a regime currently unexplored by previous Lyα LFs determinations. Our errors are dominated by the completeness correction at low luminosity, while poor statistics due to low number counts (i.e. Poissonian errors) dominate the bright tail of our distributions.

The shaded grey areas in each panel of Fig. 12 mark the regions which are not accessible by our data due to the limiting L_Lyα (vertical limit, see Table 2) and the survey area (horizontal limit). In particular, the latter marks the comoving number density (per Δ Log L_Lyα) obtained if only a single object were detected in the whole survey footprint. Errors on Φ_LLyα are computed as described in Sect. 4.4.4, and show a clear prevalence of the completeness correction at the lowest luminosity bins. On the other hand, the bright ends of our LFs are dominated by the internal variance of our samples as the number density of our candidates approaches the survey limit. To stress the impact of low sample sizes on the bright tail, we marked with faded colours the data points at ≤​1 dex above the density limit.

5.2. Comparison with previous determinations

We compare our Lyα LFs to a collection of previous determinations at similar z, after uniforming their underlying cosmology to the PLANCK2015 one. This task is complicated by the significant differences between the technical features of J-PLUS and previous high-z Lyα surveys (as noted in e.g. Blanc et al. 2011). Indeed, these can reach up to approximately five magnitudes in depth and a factor of 10³ on the surveyed area (see e.g. Ouchi et al. 2008; Konno et al. 2016). Nevertheless, the comparisons at z​ ∼ ​2.25 and z​ ∼ ​2.37 (respectively, J0395 and J0410 NBs) are remarkable, showing an overlap of our faint end with the works of Konno et al. (2016) and Matthee et al. (2017b).

We report the best fits from Sobral et al. (2018a) at each redshift, as these highlight both the Schechter and power-law components of the LFs (respectively, light-green and dark-green dashed lines in Fig. 12). These are obtained from a mixed Schechter/power-law model adapted to Log(L_Lyα/ erg s⁻¹)​ ≲ ​44.5 data, showing a transition between the two regimes at Log(L_Lyα/ erg s⁻¹)​ ∼ ​43.5. Despite the small overlap of luminosity regimes, our z​ ∼ ​2.25 LF shows remarkably good agreement with the power-law of Sobral et al. (2018a), as shown in the upper-left panel of Fig. 12. Interestingly, this component efficiently accounts for the population of X-ray bright objects in their samples, suggesting that these sources might belong to a separate class described by a different luminosity distribution from SF LAEs at Log(L_Lyα/ erg s⁻¹)​ ≲ ​43.3. On the contrary, a significant discrepancy between our data and the power-law components is evident at higher z. We ascribe this to the wider separation between the L_Lyα ranges probed by our data and those on which the fits of Sobral et al. (2018a) are obtained at these z.

We note that our z​ ∼ ​2.25 data also nicely complement the bright-end determination of Konno et al. (2016) (orange dots in the upper-left panel of Fig. 12). This latter work clearly shows an excess with respect to the exponential decay of a Schechter function at Log(L_Lyα/ erg s⁻¹)​ ≳ ​43. Their explanation relies on the contribution of a population of Lyα-emitting AGNs/QSOs, as in e.g. Matthee et al. (2017b) and Sobral et al. (2018a). By joining these hints to the results of our spectroscopic follow-up and our sample analysis (Sects. 3.4 and 3.5), our work further supports the picture according to which Lyα-emitting AGNs/QSOs are responsible for the bright-end excess observed on the 2​ ≲ ​z​ ≲ ​3 Lyα luminosity function at 43.3​ ≲ ​Log (L_Lyα/ erg s⁻¹)​ ≲ ​44.5.

Figure 12 additionally shows the Lyα LF of all the DR14 QSOs in the J-PLUS footprint (from Pâris et al. 2018), with spectroscopic redshift in the intervals sampled by each NB (red pentagons). We obtain this determination by performing synthetic photometry of SDSS QSOs with J-PLUS filters and applying the same flux corrections as those computed for our data (see Sect. 4.1). For simplicity, we only associate Poissonian errors to the SDSS LF.

Despite the comparison being only qualitative, the agreement between the SDSS QSO distribution and our data is good, especially at low z. Interestingly, the fraction of our genuine candidates showing SDSS QSO counterparts at the redshift probed by the NBs is ≲30% for each NB. Assuming that the Pâris et al. (2018) catalogue represents a ∼​100% complete sample of QSOs and considering the low fraction of SDSS QSOs in our data, the agreement between the two LFs could be explained in terms of a significant residual contamination of our samples (∼​70%). Nevertheless, this is in contrast with both our purity estimates and our spectroscopic follow-up (Sects. 4.3 and 3.4.2). A more interesting explanation is that our NB-based selection might actually be sensitive to high-z QSOs, which lack spectroscopic determination in SDSS (due e.g. to their BB colours; see Ross et al. 2012; Richards et al. 2015), such as those confirmed by our follow-up programs. Indeed, their previous classification based on SDSS photometry and morphology would identify most of them as compact objects (namely stars, see Table B.1). We suggest that this misclassification might originate from the SDSS target selection based on BB-colours, which might miss the presence of emission lines. On the contrary, our selection targets photometric excesses with respect to a continuum estimate, and therefore can efficiently select high-z line emitters.

5.3. Lyα LF parameters

5.3.1. The faint-end slope: power-law or double-Schechter?

As suggested by e.g. Konno et al. (2016), Matthee et al. (2017b), Sobral et al. (2018a,b) and Calhau et al. (2020), the population of bright Lyα-emitting sources at Log(L_Lyα/erg s⁻¹)​ >  ​43 is likely to be composed of a mixture of SF LAEs and AGNs/QSOs. In particular, Matthee et al. (2017b) and Sobral et al. (2018a) suggest that the two source classes might be described by substantially different distributions in terms of typical number density and Lyα luminosity. Interestingly, the power-law component of their studies can be explained as the faint end of a Schechter function (Schechter 1976, see also Eq. (D.2)) describing the luminosity distribution of QSOs. Our data can effectively support this hypothesis by providing the bright-end complement to the AGN/QSO Schechter distribution. At the same time, our analysis is limited by the J-PLUS depth which prevents us from constraining the faind-end slope at Log(L_Lyα/erg s⁻¹)​ ≲ ​44. This might significantly influence the determination of our Schechter parameters given their mutual correlation. Instead of fixing the faint-end slope to a fiducial value (as in e.g. Gunawardhana et al. 2015; Sobral et al. 2018a), we compute it by jointly exploiting our data and previous Lyα LF determinations over the whole interval 41.5​ ≲ ​Log(L_Lyα/erg s⁻¹)​ ≲ ​44. More specifically, we make use of the Schechter component from Sobral et al. (2018a) at each redshift to describe the Lyα LF at Log(L_Lyα/erg s⁻¹)​ ≲ ​43.3 and combine it with a second Schechter function to account for Log(L_Lyα/erg s⁻¹)​ ≳ ​44.. We then vary the faint-end slope of the latter, and for each α we jointly fit the complete double-Schechter model to both our data and all the literature determinations (see Fig. 13). Finally, for each NB we obtain α and its errors from the reduced χ² distribution of the double-Schechter fits. The α values are shown in Table 6.

Fig. 13. Joint fit of our Lyα luminosity functions and literature data with a double-Schechter model (grey solid lines in each panel). This is obtained by joining the best Schechter fit from Sobral et al. (2018a) at each redshift (green dashed lines in each panel) and a second Schechter function (coloured dashed-dotted lines). We jointly fit this double-Schechter model to both our data and those from literature leaving free the parameters of the second Schechter in order to constrain its faint-end slope α at each redshift.

We further assume no evolution of α with respect to redshift because neither our data nor those of previous works allow for this to be properly constrained. Under this assumption, we obtain our final α as the weighted average of the above values: α = −1.35 ± 0.84. This high uncertainty is expected given the limited amount of data populating the transition regime between the two Schechter functions at Log(L_Lyα/erg s⁻¹)​ ∼ ​43.5 (see Fig. 13). Nevertheless, our procedure consistently accounts for available data over ∼​3 dex in luminosity, providing one of the first estimates of α for the Schechter LF of Lyα-emitting sources at Log(L_Lyα/erg s⁻¹)​ ≳ ​44. A few works have estimated the LF shape at these very bright regimes, usually by performing a power-law fit (e.g. Matthee et al. 2017b; Sobral et al. 2018a). Interestingly, these works determined values of and at z​ ∼ ​2.2, respectively, which are both consistent with our faint-end slope determinations at z​ <  ​2.5 within 1σ. This suggests that the power-law component observed at the bright end by previous works might be explained as the faint end of a Schechter function describing the distribution of extremely luminous Lyα-emitting sources (i.e. AGNs/QSOs). In other words, the full Lyα luminosity function at 41.5​ ≲ ​Log(L_Lyα/erg s⁻¹)​ ≲ ​44 could be effectively described by a double-Schechter model.

5.3.2. Constraints on Φ^* and L^*

We employ the fixed α computed with the above procedure to fit our data with a single-Schechter model and constrain Φ^* and L^* at Log(L_Lyα/erg s⁻¹)​ >  ​44. We stress that for this step we explicitly use only our data points. The results of this procedure are compared to literature data in Fig. 14, while the left panel of Fig. 15 directly compares our four redshift bins. We account for correlations between α and the remaining parameters by sampling the error of α (assumed to be Gaussian) with 50 000 Monte-Carlo realisations of the single-Schechter fits, from which we extract our final values and errors for Φ^* and L^*. Our results are listed in Table 6 and shown in the right panel of Fig. 15.

Fig. 14. Final Schechter fits of our Lyα LFs (coloured solid lines in each panel) performed by keeping the faint-end slope fixed to α = −1.35 ± 0.84. The coloured shaded regions in each panel mark the 1σ confidence regions for the Φ* and L* parameters obtained by sampling their associated errors, which are obtained via Monte-Carlo simulations (see Sect. 5.3.2). The literature data shown in each panel are the same as in Fig. 12.

Under the hypothesis that our samples are greatly dominated by AGNs/QSOs, our results show that their LF is described by a clearly distinct distribution with respect to SF LAEs (see also Matthee et al. 2017b). In particular, by comparing our Φ^* and L^* to previous determinations at Log(L_Lyα/erg s⁻¹)​ <  ​43 (Gronwall et al. 2007; Ouchi et al. 2008; Konno et al. 2016), we measure a typical density and luminosity of AGNs/QSOs that is ∼​​3 dex lower and ∼​​2 dex higher, respectively, as already suggested by for example Matthee et al. (2017b) and Sobral et al. (2018a). In turn, this would suggest that the transition between the regime dominated by SF LAEs and AGN/QSOs, respectively, would fall at Log(L_Lyα/erg s⁻¹)​ ∼ ​43.5, as also highlighted by Sobral et al. (2018b) and Calhau et al. (2020).

Finally, our data do not allow us to constrain the evolution of our Lyα LFs determinations. Indeed the Φ^* and L^* we obtain are statistically consistent (at ∼2σ) among the four filters, with average values of Φ^*​ = ​(3.33 ± 0.19)×10⁻⁶ Mpc⁻³ and L^*​ = ​44.65 ± 0.65 erg s⁻¹. This is shown in the right panel of Fig. 15, where the faint and dark contours for each filter respectively mark the 2-σ and 1-σ levels (i.e. the ​86% and ​39% iso-contours) of the parameter distributions obtained from Monte-Carlo realisations. The wide overlap between the four filters shows the low constraining power of our data towards the evolution of the LF parameters with redshift. This was anticipated by the significant variation among the distributions of L_Lyα and EW at each z shown in Fig. 6, which ultimately hinders the possibility of disentangling the intrinsic variations of our sample properties from systematic effects.

Fig. 15. Left panel: single-Schechter fits to our data computed with the fixed faint-end slope α​ = ​ − 1.35 ± 0.84 obtained as in Sect. 5.3.1. We note that the differences among the four determinations (factor of about two in both luminosity and normalisation) are absorbed by the errors on the Schechter parameters (right panel). Right panel: distribution of Φ* and L* obtained from the Monte-Carlo sampling of α errors. The contours mark the levels including 86% and 39% of the Monte-Carlo realisations (faint and dark contours, respectively). This analysis shows that the parameters of the four determinations are statistically consistent; we do not observe hints for an evolution of the 2​ ≲ ​z​ ≲ ​3.3 Lyα LF at Log(LLyα/erg s−1) ≳ 43.5.

5.4. The AGN fraction of z ≳ 2 LAEs

By assuming that our Lyα LF describes the distribution of only AGNs/QSOs, we can build a simple toy model to estimate the relative fraction of AGNs/QSOs and SF LAEs as a function of Lyα luminosity. We define the latter as:

(14)

where LF^AGN/QSOs is one of the four determinations of the Schechter function computed from our data, while LF^{SF LAEs} is the best fit of Sobral et al. (2018a) at the corresponding redshift. We use the latter since it is obtained by excluding LAE candidates with X-ray counterparts from the determination of the Schechter fit. Consequently, we assume it provides a fair estimate for the luminosity distribution of only SF LAEs. We underline the fact that our estimate of q_AGN is an illustrative application of our results rather than a rigorous measurement, given the strong assumptions on which it is based.

The AGN/QSO fractions for all the redshifts we probe are shown in Fig. 16. Despite our simplifying assumptions, we find good agreement (within 1σ) with the measurements of Sobral et al. (2018b), which are obtained from a spectroscopic follow-up of Lyα-selected targets. On the contrary, the works of Matthee et al. (2017b) and Calhau et al. (2020) (also shown in Fig. 16 for comparison) are based on photometric selections that identify AGN/QSO candidates on the basis of their X-ray and/or radio-loudness. The latter are likely to be significant only for a subsample of AGNs/QSOs (as suggested by e.g. Sobral et al. 2018a, and Calhau et al. 2020), and therefore the discrepancy with our estimates might also be explained in terms of this incompleteness effect.

Fig. 16. AGN/QSO fraction as a function of luminosity for each NB. We estimated this quantity by assuming that our results are entirely dominated by AGNs/QSOs and that the best Schechter fit of Sobral et al. (2018a) describes the distribution of SF LAEs (see Eq. (14)). Our results are in agreement with the spectroscopic determination of Sobral et al. (2018b), which only employs Lyα emission pre-selection for their targets. On the other hand, the estimates of Matthee et al. (2017b) and Calhau et al. (2020) are based on the detection of either X-ray or radio counterparts for their Lyα-emitting candidates.

In conclusion, the good agreement between our AGN/QSO fraction estimates and the data of Sobral et al. (2018b) supports the scenario in which our samples are strongly dominated by Lyα-emitting AGNs/QSOs. Furthermore, the discrepancy with respect to X-ray/radio-selected AGN candidates suggests that the latter are likely a subsample of the whole high-z AGN/QSO population. Our selection, on the contrary, is only based on Lyα emission, and therefore it is likely to detect previously unidentified high-z AGN/QSOs. This is also in line with the results of our spectroscopic follow-up program (Sect. 3.4.2).

6. Conclusions

This work presents the determination of the bright end of the Lyα luminosity function at four redshifts in the interval 2​ ≲ z ≲ ​3.3, specifically z​ = ​2.25,  2.37,  2.54, and z​ = ​3.24 (see also Table 2). We obtain our LFs by employing four lists of Lyα-emitting candidates selected from the DR1 catalogue of the J-PLUS survey according to the significance of their photometric excess in the J0395, J0410, J0430, and J0515 NBs.

We select 2547, 5556, 4994, and 1467 bright candidates (L_Lyα >  2 × 10⁴³ erg s⁻¹), which jointly represent the largest sample of photometric Lyα-emitting candidates at 2​ ≲ z ≲ ​3.3 to date. We expect these lists to include both bright star-forming LAEs (SF LAEs) and Lyα-emitting AGNs/QSOs. To identify either of these source classes in our samples, we carried out a spectroscopic follow-up of a random subsample of our candidates (Sect. 3.4). The spectroscopic data confirmed 40 out of 45 targets as genuine high-z line emitters (with 29 out of 45 being z​ >  ​2 Lyα-emitting QSOs) and found no star-forming LAEs. In addition, we look for bi-modalities in the photometric properties of our candidates, such as Lyα luminosity and EW (Sect. 3.5.1) or colours (Sect. 3.5.2). Overall, the properties of our candidates are consistent with those of spectroscopically confirmed QSOs (Fig. 7) and high-z QSO templates (Fig. 8), suggesting that the fraction of SF LAEs in our samples is negligible.

We use our candidate samples to compute the Lyα LF at extremely bright luminosity regimes for the first time, namely at 44​ ≲ ​Log(L_Lyα/erg s⁻¹)​ ≲ ​45.5, and extend by ≳​1.5 dex the intervals covered by previous determinations. The extensive area observed by J-PLUS DR1 allows us to access wide cosmological volumes (≳​1 Gpc³), and therefore to probe number densities as low as ∼​​10⁻⁸ Mpc⁻³. This parameter-space region is unprecedented for surveys focused on bright photometrically selected Lyα-emitting sources. Interestingly, our Lyα LFs are in line with previous results at Log(L_Lyα/erg s⁻¹)​ ≳ ​43.5, prolonging their power-law end into a fully developed Schechter function. We derive the redshift-averaged parameters Φ^*​ = ​(3.33 ± 0.19)×10⁻⁶ Mpc⁻³, L^*​ = ​44.65 ± 0.65 erg s⁻¹, and α​ = ​ − 1.35 ± 0.84 for our Schechter best fits. This shows that the whole Lyα LF, from Log(L_Lyα/erg s⁻¹)​ <  ​42 up to Log(L_Lyα/erg s⁻¹)​ >  ​45, can be effectively described by a composite model of two Schechter functions, respectively accounting for the distribution of SF LAEs and bright AGNs/QSOs. These distributions appear structurally different, with , , and a transition-regime centred at Log(L_Lyα/erg s⁻¹)​ ∼ ​43.5 (in line with e.g. Konno et al. 2016; Matthee et al. 2017b; Sobral et al. 2018b; Calhau et al. 2020). Overall, our results support the scenario suggested by e.g. Konno et al. (2016), Matthee et al. (2017b) and Sobral et al. (2018a), according to which the excess of bright LAEs measured at Log  L_Lyα ≳ 43 with respect to a Schechter distribution is due to a population of AGNs/QSOs (see also Calhau et al. 2020). Our findings characterise for the first time this population as being approximately ​100 times more luminous and about ​1000 times less dense than the SF LAEs at comparable redshifts.

In addition, ∼​70% of our Lyα-emitting candidates lack any spectroscopic confirmation by current surveys. Based on our spectroscopic follow-up results, we suggest that our samples are dominated by high-z QSOs that are not yet identified as such but are rather misclassified as stars by current archival data based on their photometric colours. Indeed, even accounting for a conservative residual contamination of ∼35% in our final samples, the number of genuine z​ ≳ ​2 QSOs identified for the first time by our methodology would be approximately 1300, 3200, 2900, and 900 for J0395, J0410, J043,0 and J0515 J-PLUS NBs, respectively. We ascribe this possibility to the NB excess detection of our methodology, which can be particularly effective in targeting and selecting the strong line-emission features of z​ >  ​2 AGNs/QSOs. Indeed, these might be missed by spectroscopic target selection based only on BB colours (e.g. Richards et al. 2009, 2015; Ivezić et al. 2014). We stress that the confirmation of this speculative hypothesis must rely on a systematic and extensive confirmation of our candidates. The latter might be obtained via either spectroscopic analysis or by exploiting the very efficient source identification yielded by multi-NB imaging. Indeed, the upcoming J-PAS survey can provide a natural setting to extend our work.

Finally, our data do not show significant evolution of the LF over the probed redshifts. Despite X-ray studies suggesting little evolution of the 2​ <  ​z​ <  ​3.3 AGN/QSO population (e.g. Hasinger et al. 2007), our findings might also be affected by J-PLUS detection limits. This factor could be mitigated by deeper photometric imaging. Unfortunately, this is unlikely to be attainable by future J-PLUS data releases, because the technical features of the T80 (80cm) telescope hinder the possibility of reaching higher image depths over very wide sky areas. On the contrary, future multi-NB wide-area photometric surveys can provide valid tools to test the LF evolution at Log  L_Lyα ≳ 43.5.

Acknowledgments

We thank the anonymous Referee for the useful comments which significantly helped to improve the manuscipt quality. This work is based on observations made with the JAST/T80 telescope for J-PLUS project at the Observatorio Astrofísico de Javalambre in Teruel, a Spanish Infraestructura Cientifico-Técnica Singular (ICTS) owned, managed and operated by the Centro de Estudios de Física del Cosmos de Aragón (CEFCA). Data has been processed and provided by CEFCA’s Unit of Processing and Archiving Data (UPAD). Funding for the J-PLUS Project has been provided by the Governments of Spain and Aragón through the Fondo de Inversiones de Teruel; the Aragón Government through the Research Groups E96, E103, and E16_17R; the Spanish Ministry of Science, Innovation and Universities (MCIU/AEI/FEDER, UE) with grants PGC2018-097585-B-C21 and PGC2018-097585-B-C22; the Spanish Ministry of Economy and Competitiveness (MINECO) under AYA2015-66211-C2-1-P, AYA2015-66211-C2-2, AYA2012-30789, and ICTS-2009-14; and European FEDER funding (FCDD10-4E-867, FCDD13-4E-2685). The Brazilian agencies FINEP, FAPESP and the National Observatory of Brazil have also contributed to this project. The spectroscopic programs in this work are based on observations made with the Gran Telescopio Canarias (GTC), installed in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias, in the island of La Palma. R.A.D. acknowledges support from the Conselho Nacional de Desenvolvimento Científico e Tecnológico - CNPq through BP grant 308105/2018-4, and the Financiadora de Estudos e Projetos - FINEP grants REF. 1217/13 - 01.13.0279.00 and REF 0859/10 - 01.10.0663.00 for hardware funding support for the J-PLUS project through the National Observatory of Brazil.

References

Appendix A: Equations of the three-filter method

Here we derive the main equations we use to extract the integrated Lyα flux from J-PLUS photometry. These were originally detailed in Vilella-Rojo et al. (2015) and similar methods are described in e.g. Pascual et al. (2007) and Guaita et al. (2010). We start by defining the monochromatic flux density of an astrophysical source (or, simply, its intrinsic spectrum) as the emitted flux per unit frequency or wavelength: f_ν or f_λ, respectively in ν-units (erg cm⁻² s⁻¹ Hz⁻¹) and λ-units (erg cm⁻² s⁻¹ Å⁻¹). The two are connected by:

(A.1)

where c is the speed of light. In this work, we only use the λ-units formalism, although magnitudes are usually defined in terms of f_ν. Photometric measurements are usually performed through filters who probe f_λ over specific wavelength intervals or pass-bands. Consequently, photometric filters are defined by their transmission curves , which describe their response⁹ as a function of wavelength. All photons received within a given pass-band during the measuring process get integrated, and therefore the details of f_λ are lost. For this reason, the flux of a source measured in a given filter “x” is effectively defined as the average flux in the pass-band weighted by the filter , i.e. . For photon-counting devices (CCD), the latter quantity is given by (e.g. Tokunaga & Vacca 2005):

(A.2)

where we assume that f_λ can be written as the combination of line and continuum emission (respectively and ). In order to extract the line flux from the measurement, we need to disentangle from . The equivalent width of a line measures the relative contribution of line and continuum to :

(A.3)

where λ_min and λ_max encompass the whole line profile. By assuming that is constant between λ_min and λ_max, and denoting the wavelength of the line-profile peak as λ_EL, we have:

(A.4)

where . This also shows the definition of the continuum-subtracted, integrated line fluxF^EL (erg cm⁻² s⁻¹).

The above definitions allow us to derive the basic equations of our methodology. We stress that this is designed to extract F^EL using three photometric measurements (two BBs and one NB) and is based on two main hypotheses: (i) the emission-line profile can be approximated by a Dirac-delta, and (ii) the source continuum is well-traced by a linear function of wavelength over the whole interval covered by the three filters (see also Sect. 3.1):

(A.5)

(A.6)

where δ(λ − λ_EL) is centred at λ_EL, while A and B are two scalar coefficients. Equation (A.5) implicitly assumes that F^EL is entirely included within the NB pass-band. This might be false when part of the emission-line profile lies outside the NB pass-band, e.g. when the line-profile is wider than the NB pass-band (as for broad QSOs lines) or when its peak lies close to the NB pass-band edge. The implications of this bias on our results are discussed in Sect. 4.1. Using (A.6) in (A.2) we get:

(A.7)

where , while the last step makes use of Eq. (A.5) in the last term at the numerator and the properties of the Dirac-delta. To simplify the notation we introduce:

(A.8)

which depend only on and λ_EL. The latter is determined by each source redshift and cannot be measured without a spectroscopic observation, and therefore we must assume λ_EL a priori. For each NB, we choose the value which maximises the product T_λ ⋅ [(dn/dz)⋅(dz/dλ)], where dn/dz is the redshift distribution of SDSS QSOs (Pâris et al. 2018). This reflects that we expect most of our candidates to be z​ ≳ ​2 QSOs, as discussed in Sects. 2.2.2 and 3.5. We can now re-write (A.7) using (A.8):

(A.9)

which is valid for a generic filter. We note that if the targeted emission-line lies outside the pass-band, then implying β_x = 0. To determine A, B and F^EL, we apply (A.9) to a set of three filters: a NB, a line-contaminated BB (denoted LC here) and a line-uncontaminated BB (denoted LU):

(A.10)

(A.11)

(A.12)

By solving this linear system we finally obtain F^EL, A and B:

(A.13)

(A.14)

(A.15)

The coefficients A and B can be used to evaluate (A.6) at the NB λ-pivot (see Tokunaga & Vacca 2005) and get an estimate of the line-uncontaminated linear continuum in the NB . Equation (1) in Sect. 3.1 details how we use this continuum estimate to compute the NB excess of each J-PLUS source. Finally, we use Eq. (A.13) to estimate the total line flux of our LAE candidates and construct their LFs as explained in Sect. 4.

Appendix B: Measurement of GTC spectra and follow-up results

We measure the redshift of all the 37 sources identified as QSOs (either at z​ ∼ ​1.5 or at z​ ∼ ​2.2) in both spectroscopic programs. We do not aim to reach a higher precision than σ_z = 10⁻² because the main goal of our follow-up programs is the spectroscopic confirmation of our targets. We additionally extract the Lyα EW and integrated line flux F_Lyα for the 29 z​ ∼ ​2.2 QSOs, from which we compute the sources Lyα luminosity.

Following well-established procedures (see e.g. Pâris et al. 2011), we first identify the main spectral lines in our QSOs spectra, such as CIV and CIII]. We then use their profile-peaks to compute our redshift estimate. We discard the Lyα profile for this analysis because it provides a systematically biased z measure because of the complex radiative transfer of Lyα photons in the source rest frame and IGM (see e.g. Gronke et al. 2016; Dijkstra 2017; Gurung-Lopez et al. 2019). We fit a double Gaussian profile to both CIV and CIII] profiles in order to simultaneously trace its broad and narrow components. We use the λ position of the narrow-component peaks to obtain two z estimates, the average of which provides the final spectroscopic z of our sources. The Lyα line flux can only be obtained after estimating the sources continua. We then fit a power law to the wavelength regions of each spectrum that is not affected by any line feature, as shown by the yellow sections of the spectrum displayed in Fig. B.1. We use the following simple functional form:

(B.1)

Fig. B.1. Calibrated spectrum (grey solid line) of the GTC2018A_09 target, confirmed as z​ ∼ ​2.2 QSO. The four green regions highlight the intervals used for computing the power-law fit to the continuum (dashed yellow line). Finally, the Lyα integrated flux is highlighted in red.

where is the spectrum monochromatic flux density (in units of erg cm⁻² s⁻¹ Å⁻¹)) while k and α are fit parameters. Finally, we measure the total Lyα line flux by integrating the excess above the estimated continuum in the wavelength range affected by the Lyα line, which is shown in Fig. B.1 as the spectral region highlighted in dark-red. As a last step, we estimate the observed Lyα EW as:

(B.2)

in which is the value of the power-law fit to the continuum at the wavelength of the Lyα line-profile peak. Figure B.1 shows the spectrum of target GTC2018A_09 as a visual example of our measuring procedure. The results of these measurements are shown in Table B.1 together with a summary of the spectroscopic follow-up technical requirements and additional properties of the observed targets.

Appendix C: Retrieval of the total line flux

Here we describe how we characterise and correct the bias on our measurements of Lyα flux (namely ). We compute a correction for each NB by exploiting the comparison between J-PLUS sources and their counterparts in the SDSS QSOs catalogue (Pâris et al. 2018), within the redshift windows sampled by the NB. For the sake of brevity, we only show the results for the J0410 and J0430 filters. As a first step, we introduce the Lyα flux quantities obtained from SDSS spectra in order to perform our comparison (see Fig. C.1):

Fig. C.1. Graphic definition of the quantities used for comparing Lyα flux. The spectrum of a z​ ∼ ​2.2 QSO from SDSS DR14 is used as example in all panels (grey lines). Each Lyα flux definition is outlined by a grey shaded area (see text for details). Yellow lines in each panel show the reconstruction of the source continuum (power-law and linear approximation, respectively, in the first two and last two panels from the top). Finally, coloured squares and crosses (third and last panel from the top, respectively) show J-PLUS measurements and synthetic photometry, respectively, performed on the SDSS spectrum with J-PLUS transmission curves (Eq. (A.2)).

• : synthetic flux-density, in erg cm⁻² s⁻¹ Å⁻¹, measured by convolving SDSS spectra with the transmission curve of a given J-PLUS filter “x”, as in Eq. (A.2). Coloured crosses in the bottom panel of Fig. C.1 mark the synthetic photometry of the SDSS QSO taken as example.

• : spectroscopic measurement of the wavelength-integrated Lyα flux, in erg cm⁻² s⁻¹. This is obtained by fitting the QSOs continuum and integrating the spectra above it, on the wavelength range affected by the whole Lyα line profile (see Appendix B for details). This measurement does not involve the convolution of SDSS spectra with the filters transmission curves.

• : a version of obtained by integrating the spectra above the continuum-fit exclusively over the wavelength range covered by a J-PLUS NBs. If the Lyα line of a given QSO is wider than the NB, this measurement is just a fraction of , because the Lyα flux lying outside the transmission curve would not be accounted for (see Fig. C.1).

• : this photometric quantity is analogous to , with the only difference being that it is computed on the synthetic photometry of SDSS QSOs (coloured crosses in the bottom panel of Fig. C.1).

The photometric quantity is directly comparable to the spectroscopic because the method of Vilella-Rojo et al. (2015) is designed to remove the effect of the filter transmission curve (see also Appendix A). Moreover, any photometric measurement is only sensitive to the flux received within a given band, hence to and not to , in our case.

C.1. r-band auto-aperture flux

Figure C.2 shows the comparison between r-band flux and for J-PLUS sources and their QSOs counterparts with Lyα in the J0410 filter. The left panel displays no systematic shift, hence no strong aperture bias (≲0.2 σ_r) affects the auto-aperture flux of point-like sources. Since z​ ≳ ​2 sources appear point-like in J-PLUS (see Sect. 2.2.3), we conclude that auto-aperture photometry collects the total light of our Lyα-emitting candidates. On the other hand, the spread of the distribution in the right panel is significantly greater than one, and therefore we need to account for this additional statistic uncertainty in addition to the J-PLUS photometric errors when computing . We re-scale the photometric errors of r band photometry and propagate the resulting σ_r on . The latter is accounted for in the errors of our LFs as discussed in Sect. 4.4.4.

Fig. C.2. Blue solid lines in both panels show the histograms of the difference between J-PLUS r flux (i.e. ) and obtained from SDSS QSOs spectra (see Fig. C.1). The distributions are normalised by respectively (left panel) and its photometric error (right panel). Both distributions are centred at zero (see plot legends), meaning that and values are statistically equivalent. On the other hand, the distribution spread in the right panel is significantly bigger than one, and therefore photometric errors do not fully account for the flux difference. Consequently, we re-scale to the value obtained by the Gaussian fit (right-panel legend).

C.2. Filter-width effect on the line flux

As the observed Lyα line profile of QSOs is generally wider than the FWHM of J-PLUS NBs, we need to account for line-flux losses. For this, we obtain two corrections for and its errors¹⁰ , respectively. These affect the Lyα luminosity of our candidates and its errors, and therefore we account for them on our final LFs (as discussed in Sect. 4.4.4). Figure C.3 shows the comparison between and for J0430 filter. Their flux difference is presented in the left and middle panel, respectively normalised by and by . The left panel shows clear evidence of a systematic offset ΔF between the two flux quantities. We measure ΔF with a Gaussian fit (see Fig. C.3) and use it to correct as follows:

(C.1)

Fig. C.3. Cyan solid lines in all panels show the histograms of differences between line flux measurements performed on J-PLUS photometry (Sect. 3.1) and on SDSS spectroscopy (Fig. C.1), for the case of J0430. The difference between and is normalised respectively by and by in the left and middle panels. We use the systematic shift of the distributions in the left panel to statistically correct (see Sect. 4.1.2). The spread of the distribution in the middle panel, on the other hand, allows the rescaling of the uncertainty associated with measurements in order to account for residual statistical errors not included in . Section 4.4.4 details how we make use of the corrected and in the computation of the errors associated to our final LFs. Finally, right panel shows that well compares to the spectroscopic measure .

We stress that by directly comparing to , our statistical correction accounts for any systematic bias of our measurements, such as the linear-continuum approximation or the line-peak position of Lyα within the NB. The spread of the distribution in the middle panel of Fig. C.3 is significantly bigger than unity. This shows that cannot fully account for the difference between and . Consequently, we re-scale according to the measured spread of the distribution in the middle panel of Fig. C.3. Finally, we account for these errors on our final LFs (see Sect. 4.4.4). In conclusion, the comparison between middle and right panels of Fig. C.3 clearly shows how better compares to the fraction of spectroscopic line-flux measured only on the wavelength range covered by the J-PLUS NB (i.e. , see Fig. C.1 and Appendix C). This is a direct effect of the filter-width bias, because our methodology is only sensitive to the flux captured by J-PLUS photometry within the NB wavelength range, as discussed in Appendix C.

Appendix D: Multi-variate completeness computation

Here we describe the computation of the corrections accounting for incompleteness due to our selection methodology (Appendix D.1) and the use of r-band detected catalogues (Appendix D.2).

D.1. Selection completeness

To simplify, our selection depends on (i) the linear approximation of the sources continuum, (ii) their Lyα flux, and (iii) their Lyα EW. To account for these dependencies, we measure the recovery rate of our methodology as a function of (i) g​ − ​r colour, (ii) Lyα flux, and (iii) r-band magnitude. Indeed, g​ − ​r colour can be thought of as a proxy of our linear continuum approximation (see e.g. Fig. 2), and therefore the EW dependence of our selection is accounted for by independently varying the Lyα flux with respect to g​ − ​r and r magnitude. Specifically,

1. we first subtract the measured from the sources photometry of our candidates in order to get a list of non-emitters, i.e. sources without a significant NB excess according to our measuring method.

2. We then artificially re-add increasing values of line flux to the sources photometry (i.e. to the NB and the g band, as both are affected by the emission line).

3. For each value of the re-added flux, we apply our complete set of selection rules (Sect. 3.2) and we store the number of re-selected sources as a function of their r magnitude and g​ − ​r colour.

4. We finally compute the sources-recovery rate of our selection as C = N_selected/N_total in each bin of r magnitude, g​ − ​r colour and artificially-injected Lyα flux.

With this method we obtain a 3D grid of recovery rates which can be interpolated in order to compute a selection weight for each source. The latter accounts for the loss of candidates due to our selection within the total completeness correction we apply to the final LFs. The computation of depends on Lyα-flux since this is the observable tackled by our measuring method (see e.g. Sects. 3.1 and Appendix A). Nevertheless, this dependence can be directly converted into a L_Lyα dependence by assuming a redshift for every source (see Sect. 4.2). Figure D.1 shows the recovery-rates grid for J0430 filter as an examples. The values are projected in the planes F_Lyα vs. r (left panel) and F_Lyα vs. g​ − ​r (right panel).

Fig. D.1. 3D grid of recovery rates for the J0430 filter shown as an example. Left and right panels respectively show the projections of recovery rates in the Lyα flux vs. r plane and Lyα flux vs. g​ − ​r plane. We note that the recovery rates show noisy values at r​ <  ​17 and g​ − ​r​ <  ​ − 1.25 due to the low number of sources in these magnitude and colour bins. Nevertheless, these regions of the 3D parameter space are excluded from the LF computation by the purity weight (Sect. 4.3).

D.2. Bivariate completeness model

Since we use r-band detected catalogues, we need to take into account the loss of undetected continuum-faint Lyα-emitting sources. In other words, we need to estimate the distribution of sources in regions of the L_Lya vs. r plane that lie outside the thresholds on r and NB-excess significance imposed by J-PLUS detection limits. To carry out this analysis we closely follow the methods of Gunawardhana et al. (2015), who tackled a very similar issue with a multi-variate approach. Specifically, we build the bivariate luminosity function of our candidates, which is defined as the number density of sources in each bin of Log (L_Lyα) and r, weighted by the purity and completeness corrections:

(D.1)

where j and k indexes identify the 2D bins of Log(L_Lyα) and r, while the index i runs over the total number of sources in each 2D bin. P_i is the purity weight of each candidate (see Sect. 4.3), while and are respectively its detection-completeness (Sect. 4.4.1) and selection-completeness (Sect. 4.4.2 and Appendix D.1) weights. Finally, V is the survey effective volume for a given NB filter (see Sect. 4.2 and Table 2). By following Gunawardhana et al. (2015), we assume that the 2D LF can be modelled by the product of two functions, describing respectively the r and Log(L_Lyα) distributions (see also Corbelli et al. 1991). We choose to employ the combination of a Schechter function (in logaritmic form) for r and a Gaussian in Log L_Lyα (as in Gunawardhana et al. 2015):

(D.2)

where , r^*, α_r are the ordinary Schechter parameters (see Schechter 1976), while , L^*, σ_L describe respectively the number-density normalization, the average luminosity and the spread of the L_Lyα distribution in each r bin. In order to obtain the bivariate model, we follow Gunawardhana et al. (2015) and join the two univariate distributions presented in Eq. (D.2) with an equation between their structural parameters:

(D.3)

where A and B are free parameters to be determined by fitting the 2D model to our data, while (as in Gunawardhana et al. 2015). By substituting Eq. (D.3) into (D.2) and multiplying the two univariate distributions, we obtain the full bivariate model:

(D.4)

We fit the seven free parameters of this function using our measured 2D luminosity function. In particular, we only use the candidates with a completeness weight in order to avoid biasing our fit with regions affected by the incompleteness of our selection. Finally, the ratio between the measured 2D LF of our candidates and the fitted model (over the whole r − L_Lyα plane) provides an estimate of the incompleteness of our selection in each [Log(L_Lyα),r] bin. We use this ratio to compute a bivariate weight for each of our candidates. These weights are then combined to and in order to obtain our total completeness correction (see Sect. 4.4).

All Tables

All Figures

	Fig. 1. Measured transmission curves for the J-PLUS filter set after accounting for sky absorption, CCD quantum efficiency, and the total effect of the JAST/T80 telescope optical system. The four NB we exploit to look for bright Lyα emitters at z​ >  ​2 (namely, the J0395, J0410, J0430 and J0515 NBs) share their wavelength coverage with the g band and are shown here as filled-area curves.
In the text

Fig. 2. Representation of our NB excess detection method. Grey lines in both panels show the observed spectra of typical z​ ∼ ​2 Lyα-emitting sources (from the publicly available VUDS DR1 spectroscopic dataset, see e.g. Le Fèvre et al. 2015; Tasca et al. 2017). Upper panel: SF LAE spectrum showing a single prominent Lyα line (here redshifted at λobs ∼ 3900 Å) and no other significant features. Bottom panel: QSO spectrum with evident CIV and CIII] lines in addition to Lyα (at λobs ∼ 4000 Å). We show the transmission curves and associated synthetic photometry of four J-PLUS bands as coloured lines and squares. From left to right: u (purple), J0395 NB (violet), g (green), and r (red). The additional lines of this caption were moved to Sect. 3.1

In the text

Fig. 3. Example of a colour-magnitude diagram obtained for the NB filter J0410 on a DR1 pointing (out of 511). Our photometric cuts are summarised as follows: the blue dashed-dotted line shows the ΔmNB-significance threshold, while the vertical red line marks the NB S/N limit. We exclude sources below the blue dashed-dotted line and inside the grey shaded area. The orange horizontal dotted line shows ΔmNB associated to EW = 50 Å (see Eq. (7)). Grey-blue dots mark all the J-PLUS detections in the pointing, while red and purple crosses show z​ ∼ ​2.4 QSOs and low-z galaxies from SDSS DR14. Yellow triangles show J-PLUS mock data of z​ ∼ ​2.4 SF LAEs (Izquierdo-Villalba et al. 2019). Finally, our Lyα-emitting candidates are shown as green dots.

In the text

Fig. 4. Examples of multiple NB excess in J-PLUS photo spectra. Empty and filled squares mark NB and BB photometry, respectively, while the dashed yellow line shows the linear continuum we estimate through g and r BBs (in green and red, respectively). This comparison shows that both a z​ ∼ ​0.05 galaxy (upper panel) and a z​ ∼ ​2.25 QSO (lower panel) exhibit significant excesses in J0395 and J0515 NBs (second and fifth empty squares from the left, respectively). The J0515 excess is produced by H β at z​ ∼ ​0.05 and the CIV line at z​ ∼ ​2.25, but its nature is hardly distinguishable using J-PLUS photometry.

In the text

Fig. 5. Summary of our spectroscopic results showing one spectrum (grey line in each panel) for each source class identified in our target lists. From top to bottom: z​ ∼ ​2.2 QSO, z​ ∼ ​1.5 QSO, star, z​ ∼ ​0.5 galaxy, and Lyβ-emitting QSO. The corresponding J-PLUS photometry is shown as coloured squares. The star and galaxy targets show low-significance excesses in J0395 NB (third square from the left). Indeed, these interlopers were selected as targets by the first version of our methodology, i.e. before applying the improvements due to CTG2018A results and the re-calibration of J-PLUS data (López-Sanjuan et al. 2019a).

In the text

Fig. 6. Left figure: EW0 distribution of our selected candidates as a function of r and g​ − ​r colour (left and right panels, respectively). Squared points and error bars show respectively the distribution median and 16th–84th percentiles in each magnitude and colour bin. Points have been artificially shifted for better visualisation. The values for J0395 filter at r​ <  ​20 (g − r​ <  ​0.75) are systematically below the theoretical cut we apply (Sect. 3.2). This is due to the small overlap between this NB and the g BB, which is reflected in a poor extrapolation of the linear continuum at the NB filter pivot wavelength (see e.g. Ouchi et al. 2008, and the discussion in Sect. 3.2). Right figure: normalised distributions of our candidates in EW and Lyα luminosity for each filter. This result clearly shows that filters sampling higher redshifts also sample brighter Lyα luminosity. This is a direct effect of J-PLUS detection limits which only allow to observe brighter and rare objects at higher redshifts. We address this issue by applying the completeness corrections described in Sect. 4.4.

In the text

Fig. 7. Comparison between the colour–colour distributions of DR14 QSOs (green squares) within our J0410 sample and of the remaining J0410 genuine candidates (blue dots, after removing known interlopers; Sect. 3.3). The two source classes occupy comparable colour regions, suggesting that our selection results might be effectively dominated by z​ ∼ ​2.3 AGN/QSOs. This scenario is also supported by the results of our spectroscopic program (Sect. 3.4).

In the text

Fig. 8. Coloured points: r − W1 colour vs. redshift for all our candidates with a counterpart in WISE. Our sources are uniformly spread within the z interval sampled by each NB. Grey lines: tracks for different galaxy templates (from the SWIRE library, Polletta et al. 2007) and the QSO template of Hernán-Caballero et al. (2016). Black-contoured and yellow-contoured points mark the SDSS QSOs in our selection and the GTC targets (only for J0395 NB), respectively. Our candidates are all compatible with the high-z QSO template, suggesting that the fraction of SF LAEs in our sample is very low. In addition, this suggests that we identify a large number of sources without previous spectroscopic identification as high-z QSOs.

In the text

Fig. 9. Statistical purity weight for each NB (coloured solid lines) as estimated by fitting an error-function to the computed purity. The grey-dotted line shows the computed purity of J0430 NB as an example. All filters show similar purity weights, rising to ≳60% at r ≳ 18.5. This is in agreement with the average purity of our 45 spectroscopic targets (Sect. 3.4.2) which is shown as a purple empty square.

In the text

	Fig. 10. Example of the recovery fraction of our selection as a function of Lyα flux, computed for J0430 NB in a bin of r magnitude and g​ − ​r colour (namely at r = 20.1 and g​ − ​r = 0.2). The full 3D grid is shown in Appendix D for the same NB. All filters show comparable values of recovery fractions, and therefore we only report the case of J0430, for brevity.
In the text

Fig. 11. Central panel: full 2D luminosity function of our J0430 candidates as a function of r and Log(LLyα). Green solid lines in the top and right panels show the projections of the 2D LF along the Log(LLyα) and r axes, respectively. The red dashed lines show the projection of the 2D model along the same axes. This model was fitted to the 2D distribution shown in the central panel (see Appendix D for computational details) and allows us to extrapolate our data distribution at faint r and Lyα luminosity. We use the data-to-model ratio (in the r vs. Log(LLyα) 2D plane) to compute the weight of each candidate.

In the text

Fig. 12. Lyα luminosity functions for each of the NB filters we used in our study (coloured squares). The grey shaded areas show the Lyα luminosity limit (vertical limit) and the limiting number density measurable by J-PLUS (horizontal limit). The wide area explored by the NBs of the J-PLUS survey allow us to remarkably extend the range of luminosity sampled by previous studies (coloured circles, triangles, hexagons, and diamonds in each plot) and to explore previously unconstrained LLyα intervals. Dashed lines mark the best-fit determinations by Sobral et al. (2018a), split respectively into a Schechter (light green) and a power-law contribution (dark green). Our results provide tight constraints at 44.5 ≲ Log (LLyα/erg s−1)≲45.5, a regime currently unexplored by previous Lyα LFs determinations. Our errors are dominated by the completeness correction at low luminosity, while poor statistics due to low number counts (i.e. Poissonian errors) dominate the bright tail of our distributions.

In the text

Fig. 13. Joint fit of our Lyα luminosity functions and literature data with a double-Schechter model (grey solid lines in each panel). This is obtained by joining the best Schechter fit from Sobral et al. (2018a) at each redshift (green dashed lines in each panel) and a second Schechter function (coloured dashed-dotted lines). We jointly fit this double-Schechter model to both our data and those from literature leaving free the parameters of the second Schechter in order to constrain its faint-end slope α at each redshift.

In the text

Fig. 14. Final Schechter fits of our Lyα LFs (coloured solid lines in each panel) performed by keeping the faint-end slope fixed to α = −1.35 ± 0.84. The coloured shaded regions in each panel mark the 1σ confidence regions for the Φ* and L* parameters obtained by sampling their associated errors, which are obtained via Monte-Carlo simulations (see Sect. 5.3.2). The literature data shown in each panel are the same as in Fig. 12.

In the text

Fig. 15. Left panel: single-Schechter fits to our data computed with the fixed faint-end slope α​ = ​ − 1.35 ± 0.84 obtained as in Sect. 5.3.1. We note that the differences among the four determinations (factor of about two in both luminosity and normalisation) are absorbed by the errors on the Schechter parameters (right panel). Right panel: distribution of Φ* and L* obtained from the Monte-Carlo sampling of α errors. The contours mark the levels including 86% and 39% of the Monte-Carlo realisations (faint and dark contours, respectively). This analysis shows that the parameters of the four determinations are statistically consistent; we do not observe hints for an evolution of the 2​ ≲ ​z​ ≲ ​3.3 Lyα LF at Log(LLyα/erg s−1) ≳ 43.5.

In the text

Fig. 16. AGN/QSO fraction as a function of luminosity for each NB. We estimated this quantity by assuming that our results are entirely dominated by AGNs/QSOs and that the best Schechter fit of Sobral et al. (2018a) describes the distribution of SF LAEs (see Eq. (14)). Our results are in agreement with the spectroscopic determination of Sobral et al. (2018b), which only employs Lyα emission pre-selection for their targets. On the other hand, the estimates of Matthee et al. (2017b) and Calhau et al. (2020) are based on the detection of either X-ray or radio counterparts for their Lyα-emitting candidates.

In the text

	Fig. B.1. Calibrated spectrum (grey solid line) of the GTC2018A_09 target, confirmed as z​ ∼ ​2.2 QSO. The four green regions highlight the intervals used for computing the power-law fit to the continuum (dashed yellow line). Finally, the Lyα integrated flux is highlighted in red.
In the text

Fig. C.1. Graphic definition of the quantities used for comparing Lyα flux. The spectrum of a z​ ∼ ​2.2 QSO from SDSS DR14 is used as example in all panels (grey lines). Each Lyα flux definition is outlined by a grey shaded area (see text for details). Yellow lines in each panel show the reconstruction of the source continuum (power-law and linear approximation, respectively, in the first two and last two panels from the top). Finally, coloured squares and crosses (third and last panel from the top, respectively) show J-PLUS measurements and synthetic photometry, respectively, performed on the SDSS spectrum with J-PLUS transmission curves (Eq. (A.2)).

In the text

Fig. C.2. Blue solid lines in both panels show the histograms of the difference between J-PLUS r flux (i.e. ) and obtained from SDSS QSOs spectra (see Fig. C.1). The distributions are normalised by respectively (left panel) and its photometric error (right panel). Both distributions are centred at zero (see plot legends), meaning that and values are statistically equivalent. On the other hand, the distribution spread in the right panel is significantly bigger than one, and therefore photometric errors do not fully account for the flux difference. Consequently, we re-scale to the value obtained by the Gaussian fit (right-panel legend).

In the text

Fig. C.3. Cyan solid lines in all panels show the histograms of differences between line flux measurements performed on J-PLUS photometry (Sect. 3.1) and on SDSS spectroscopy (Fig. C.1), for the case of J0430. The difference between and is normalised respectively by and by in the left and middle panels. We use the systematic shift of the distributions in the left panel to statistically correct (see Sect. 4.1.2). The spread of the distribution in the middle panel, on the other hand, allows the rescaling of the uncertainty associated with measurements in order to account for residual statistical errors not included in . Section 4.4.4 details how we make use of the corrected and in the computation of the errors associated to our final LFs. Finally, right panel shows that well compares to the spectroscopic measure .

In the text

Fig. D.1. 3D grid of recovery rates for the J0430 filter shown as an example. Left and right panels respectively show the projections of recovery rates in the Lyα flux vs. r plane and Lyα flux vs. g​ − ​r plane. We note that the recovery rates show noisy values at r​ <  ​17 and g​ − ​r​ <  ​ − 1.25 due to the low number of sources in these magnitude and colour bins. Nevertheless, these regions of the 3D parameter space are excluded from the LF computation by the purity weight (Sect. 4.3).

In the text