© ESO 2021

1 Introduction

Luminous Lyman alpha emitters (LAEs) at high redshift offer valuable probes of galaxy evolution and cosmology in the early Universe. The existence of young galaxies in the early phase of their evolution with strong Lyman alpha (Lyα) emission (1216) was first hypothesised by Partridge & Peebles (1967) and first observed by Hu & McMahon (1996) and Pascarelle et al. (1996). Dedicated narrowband, spectroscopic, and integral field unit surveys have since discovered thousands of LAEs at redshifts 2 < z < 7, opening up a new way of studying the high-redshift Universe (e.g. Cowie & Hu 1998; Rhoads et al. 2000; Ouchi et al. 2003; Ciardullo et al. 2011; Shibuya et al. 2012; Konno et al. 2014; Sobral et al. 2018a). Previous studies have identified two populations of LAEs: blue and faint LAEs with low masses and metallicities (e.g. Bacon et al. 2015; Sobral et al. 2015a, 2019; Nakajima et al. 2016; Ono et al. 2010) and red, massive, and luminous LAEs (e.g. Chapman et al. 2005; Sandberg et al. 2015; Matthee et al. 2016) often observed to host an active galactic nucleus (AGN; e.g. Ouchi et al. 2008; Konno et al. 2016; Sobral et al. 2017).

Lyman alpha emitters have been widely used to study luminosity functions (LFs) and clustering properties of galaxies in the early Universe (e.g. Ouchi et al. 2003, 2010; Shimasaku et al. 2004; Kusakabe et al. 2018; Khostovan et al. 2019). Recent studies, for example Khostovan et al. (2019), suggest that LAEs are progenitors of a wide range of galaxy types, where the brightest LAEs are located in the most massive halos and are highly clustered. The LAEs that have been detected at radio wavelengths have been found to have steep radio spectral indices (Calhau et al. 2020).

Furthermore, studies have shown that the AGN fraction rises with Lyα luminosity, leading to the AGN dominating (over the star-forming population) at the bright end of the Lyα LF (e.g. Sobral et al. 2018b,a; Matthee et al. 2017; Wold et al. 2014, 2017; Calhau et al. 2020). As suggested by Sobral et al. (2018b) and demonstrated by Calhau et al. (2020), the AGN fraction of LAEs declines towards higher redshifts at a fixed Lyα luminosity. However, these conclusions are limited by the number of detected LAEs at high redshift, and larger samples are required to confirm this claim. Constraining the AGN fraction in the Lyα LF is vital for our understanding of early super massive black hole (SMBH) formation (e.g. Calhau et al. 2020), but it is also critical for understanding the source of ionising photons of LAEs if they are to be used as probes of reionization (e.g. Santos et al. 2016; Matthee et al. 2015).

The Low Frequency Array (LOFAR Two-metre Sky Survey (LoTSS; Shimwell et al. 2017, 2019) is entering a new regime of deep, low-frequency surveys by pushing noise levels to below 100 μJy beam⁻¹ at 150 MHz across the entire northern sky and complementing this with targeted, deeper observations in the degree-scale northern deep fields. The LoTSS Deep Fields first data release (Tasse et al. 2021; Sabater et al. 2021; Kondapally et al. 2021; Duncan et al. 2021) covers a total area of 25 deg² in the Lockman Hole, European Large Area Infrared Space Observatory Survey-North 1 (hereafter ELAIS-N1), and Boötes fields, at 6″ resolution, reaching a root mean square (rms) depth of S_150MHz ~ 20 μJy beam⁻¹ in the deepest field, ELAIS-N1.

One of the LoTSS target fields, ELAIS-N1 (RA = 242.75, Dec = 54.95 degrees), has also been the target field for a dedicated LAE survey named the Systematic Identification of LAEs for Visible Exploration and Reionization Research Using Subaru HSC (SILVERRUSH; Ouchi et al. 2018). The SILVERRUSH program makes use of the narrowband observations of the Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP; Aihara et al. 2018a). Three major scientific goals of the SILVERRUSH program are studying LAE properties at high redshift, using LAEs to probe the low-mass young galaxy population, and studying the Lyα LF towards the epoch of reionization (EoR; Ouchi et al. 2018). The SILVERRUSH program identified ~ 2000 LAE candidates at z = 5.7 and z = 6.6 across four different fields, ~600 of which are located in ELAIS-N1 (Shibuya et al. 2018a). The overlapping LoTSS and SILVERRUSH observations of ELAIS-N1 open up a new opportunity to study the as yet unexplored low-frequency (150 MHz) radio properties of LAEs. Furthermore, the identification of LAE radio AGN would allow us to study the AGN fraction as a function of Lyα luminosity and follow up on the claim that the LAE AGN fraction declines towards higher redshifts. Finally, according to model predictions by Saxena et al. (2017), more than ten radio-loud AGN at z > 6 are expected to be detected in all the LOFAR Deep Fields combined with the current sensitivity. Currently, the most distant radio galaxy has a redshift of z = 5.72 (Saxena et al. 2018a). The SILVERRUSH catalogue is a potentially excellent source for finding radio galaxies at z > 6, even thoughonly a small fraction of the volume is probed due to the detection range of the narrowband filters. The detection of such a high-redshift radio galaxy would not only be a substantial step forwards in the study of the formation and evolution of massive galaxies into the EoR, but it would also be the most distant radio galaxy discovered to date.

This paper is structured as follows. In Sect. 2, we describe the details of the SILVERRUSH program, the LoTSS Deep Fields data used, and the LAE sample selection. In Sect. 3, we present the characteristics of the radio-detected SILVERRUSH population. Subsequently, in Sect. 4, we analyse the optical and infrared (IR) properties and multi-wavelength stacks of the wider LAE population. In Sect. 5, we discuss the LOFAR detection rate and LAE selection. Finally, in Sect. 6, we summarise our findings.

In this work, a flat lambda cold dark matter (Λ-CDM) cosmology isassumed using H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_M = 0.3, and Ω_Λ = 0.7. Furthermore, all magnitudes presented are given in the AB system (Oke & Gunn 1983).

2 Data and sample selection

2.1 SILVERRUSH program

To identify LAEs, the SILVERRUSH survey uses four narrowband (hereafter NB) filters, which are mounted on the Hyper Suprime-Cam (HSC). These filters are NB387, NB816, NB921, and NB101, corresponding to central wavelengths (λ_cen) of 3858, 8169, 9204, and 10 092 ¹, which allow for the identification of LAEs at redshifts of z = 2.17 ± 0.02, 5.72 ± 0.05, 6.57 ± 0.05, and 7.30 ± 0.04, respectively, in addition to rest frame optical emission line galaxies at lower redshifts (see Aihara et al. 2018a). Using NB816 and NB921 imaging data from the HSC-SSP survey, the SILVERRUSH program identified ~2000 LAE candidates at z = 5.7 and z = 6.6 in five fields, with a total area of 14 and 21 deg², respectively (Shibuya et al. 2018a). These resulting LAE catalogues are publicly available² and will be referred to as the NB816 and NB921 LAE catalogues in this work. The public data release does not include the results of LAE studies based on NB387 and NB101 observations. The detailed SILVERRUSH program strategy is described by Ouchi et al. (2018), and the first SILVERRUSH catalogues and properties are presented by Shibuya et al. (2018a). Other main scientific results following from SILVERRUSH, such as clustering properties and Lyα LFs, have been published in Shibuya et al. (2018b), Konno et al. (2018), Harikane et al. (2018, 2019), Inoue et al. (2018), Higuchi et al. (2019), and Kakuma et al. (2019). In the following sections, we focus on the data and sample selection in the ELAIS-N1 field. The data obtained by the HSC-SSP survey are outlined in Sect. 2.1.1, and the SILVERRUSH procedure for LAE selection from this survey is summarised in Sect. 2.1.2.

2.1.1 HSC-SSP survey

The HSC filter transmission curves of broadbands (BBs) g, r, i, z, and y and narrowbands NB816 and NB921 are shown in Fig. 1, including the central wavelengths λ_cen and full width half maxima (FWHM) for the i and z bands and NB filters. The NB816 and NB921 LAE catalogues were created using HSC-SSP S16A data taken using these five BB and NB filters (Ouchi et al. 2018). Deep observations of ELAIS-N1 were conducted as part of the HSC-SSP survey, reaching 5σ limiting magnitudes of 25.3 for the NBs and ~24–26 for the BBs (Shibuya et al. 2018a). The image reduction was performed using the HSC pipeline (Bosch et al. 2018). The source detection and photometric measurements are obtained using a “forced” and “unforced” method. In the unforced method, the coordinates, shape, and flux of each source are determined individually in each band, whereas in the forced method the coordinates and shape of the sources are fixed in a reference band and applied to all other bands to determine the flux. In this work, the photometry from the both forced and unforced methods is used to ensure the largest possible sample. Further details on the HSC-SSP survey is available in Aihara et al. (2018a) and Shibuya et al. (2018a).

Fig. 1 Filter transmission curves of HSC BB (grizy) and NB (NB816 & NB921) filters. The central wavelength (λcen) and full width at half maximum (in nm) are indicated for the i, z, NB816, and NB921 filters. The filter transmission curves and λcen values have been obtained from the SVO Filter Profile Service (http://svo2.cab.inta-csic.es/theory/fps/index.php?mode=browse&gname=Subaru&gname2=HSC).

2.1.2 LAE selection

The LAEs in the SILVERRUSH program were selected using colour selection criteria (presented in Shibuya et al. 2018a and based on Ouchi et al. 2008, 2010), which ensure a non-detection in bands blue-wards of the Lyman break at a certain redshift and significant detected flux excess in the NBs. These criteria are given by: (1)

and (2)

for the z = 5.7 and z = 6.6 LAEs respectively, where g_3σ, r_3σ, and z_3σ are the 3σ limiting magnitudes in the g, r, and zbands. In the forced catalogue, a stricter colour selection criterion of z−NB921 > 1.8 is used, and the limiting magnitudes are lowered to 2σ in the g and r bands. These stricter colour criteria are also used in Konno et al. (2018) for studying the Lyα LF. Besides the colour selection, several parameters and flags were used to denote imaging problems or artefacts, such as sources containing saturated pixels and images with short exposure times (see Shibuya et al. 2018a; Aihara et al. 2018b). Further details on the SILVERRUSH LAE selection algorithm are presented in Shibuya et al. (2018a).

The final SILVERRUSH samples in ELAIS-N1 contain 229 and 349 sources at z = 5.7 and 6.6, respectively, in a 6 deg² area (Konno et al. 2018). The reliability of their LAE selection in all fields was checked using spectroscopic observations of a sub-sample of the selected LAEs. In total, redshifts of 96 of the HSC LAEs were spectroscopically confirmed, and they yielded a contamination rate of 0–30%, depending on the magnitude ranges (Shibuya et al. 2018b).

2.2 LOFAR radio observations of ELAIS-N1

LoTSS aims to cover the whole Northern Hemisphere, reaching ≤100 μJy beam⁻¹ noise levels; a part of the data was released in the first data release (LoTSS-DR1; Shimwell et al. 2019). LoTSS is being complemented by a series of LoTSS Deep Fields, which aim to ultimately cover a sky area of ~50 deg² down to noise levels of ~10 μJy beam⁻¹, probing the fainter and higher redshift radio population. The LoTSS Deep Fields DR1 targets the Lockman Hole, Boötes, and ELAIS-N1 fields, which are by design at declinations optimal for LOFAR observing sensitivity. These fields have already been extensively observed across optical and IR wavelengths, enabling the determination of photometric redshifts and physical galaxy properties. New imaging and calibration algorithms have been developed by Tasse et al. (2021) to enable the construction of thermal noise-limited images, and the LoTSS Deep Fields radio images reach ~20 μJy beam⁻¹ rms sensitivities in ELAIS-N1 (Sabater et al. 2021). Details on the ELAIS-N1 radio imaging and catalogue are presented in Sabater et al. (2021).

2.3 Deep optical-IR photometry in ELAIS-N1

An extensive range of multi-wavelength observations from ultraviolet (UV) to far-infrared (FIR) are available in the ELAIS-N1 field. An overview of the observations used in this work are given in Table 1. Here, we provide a summary of these observations and the generation of the multi-wavelength catalogue. We refer to Kondapally et al. (2021) for a detailed description of the catalogues used.

All LOFAR-detected radio sources were cross-matched with the multi-wavelength catalogues, and photometric measurements extracted, by Kondapally et al. (2021); full details of the procedure are available in that paper. In summary, a multi-wavelength catalogue was created using forced, matched aperture photometry on pixel-matched images from all surveys. To achieve this, all individual images were resampled to a pixel scale of 0.2″ and sky background is subtracted before the individual images were added together using SWarp (Bertin et al. 2002). The flux was adjusted to a common scale using the zero-point magnitude, exposure time, and Vega-AB conversion factors, where needed. Optimal signal to noise (or χ²) detection images were created to be able to detect the faintest sources. They were created by stacking multiple bands usingSWarp for the optical to near-infrared (NIR) and Spitzer Infrared Array Camera (IRAC) observationsseparately due to the lower resolution of the Spitzer data. The weight assigned to each band in the χ² detection image varies according to the colour of the source. The sources were extracted from the χ² detection image using SEXTRACTOR (Bertin & Arnouts 1996), and the two source lists were combined to make a single catalogue. The fluxes in the different bands were obtained from all sources detected in either of the χ² detection images using aperture sizes with diameters of 1″–7″ (in steps of 1″) and 10″. The varying point spread functions (PSF) in each filter were corrected for by aperture corrections, determined using the curve of growth estimated from moderately bright sources (see Kondapally et al. 2021). In this work, we used the flux measurements from the 3″ apertures for the optical-NIR filters and 4″ apertures for the Spitzer-IRAC filters, which have both been aperture corrected. The 3″ aperture for the optical-NIR filters is less affected by PSF variations than a 2″ aperture and will therefore result in more robust colours (Kondapally et al. 2021). The final ELAIS-N1 multi-wavelength catalogue contains over 2.1 million sources; 1.5 million of these sources in the overlapping PanSTARRS, UKIDSS-DXS, and Spitzer-SWIRE surveys are used for radio-optical cross-matching in an area of 6.7 deg². The cross-matching was carried out using an adaptation of the technique developed in LoTSS DR1 and presented by Williams et al. (2019). In short, optical and IR counterparts of radio sources were identified by either the statistical likelihood ratio method or by visual classification schemes, as determined by a decision tree described in Williams et al. (2019) and Kondapally et al. (2021). Sources with extended and/or complex radio emission were associated and classified using a combination of the LOFAR Galaxy Zoo (Williams et al. 2019) and an expert-user work-flow designed for de-blending radio sources. A detailed description of the procedure followed to create the catalogues is available in Kondapally et al. (2021).

In addition, the Spectral and Photometric Imaging Receiver (SPIRE) and Photodetector Array Camera and Spectrometer (PACS) FIR fluxes were added to the catalogue obtained from the Herschel Extragalactic Legacy Project (HELP; Oliver et al., in prep.). The procedure for the determination of the FIR fluxes is described in McCheyne et al. (in prep.). In short, mid-infrared (mid-IR) and FIR fluxes were derived from the prior driven de-blending of the Spitzer Multiband Imaging Photometer (MIPS) and Herschel PACS/SPIRE imaging using the XID+ deblending code (Hurley et al. 2017).

Photometric redshifts for all optical sources in the Kondapally et al. (2021) catalogues were estimated using a hybrid template and machine-learning method. Full details of the method and characterisation of the redshift performance are presented in Duncan et al. (2021, hereafter D20).

3 Characterising the radio-detected SILVERRUSH population

The LOFAR-detected SILVERRUSH population in ELAIS-N1 was determined by cross-matching sources from the LoTSS catalogue and the SILVERRUSH NB921 and NB816 LAE catalogues within 3″ separation (half of the LOFAR 6″ beam size). The location of the optical counterparts of the radio sources were used for cross-matching. The resulting cross-matched sources are inspected by eye, leading to five LOFAR-detected SILVERRUSH sources out of the eight potential matches, which will from now on be referred to as the “LOFAR-detected sample”. All five of these sources matched with SILVERRUSH sources in the NB921 LAE catalogue (z = 6.6), whereas no sources have been found to match with the NB816 LAE catalogue. The multi-wavelength cutouts of the LOFAR-detected sample are shown in Fig. 2.

Table 2 summarises the photometric source properties, including the photometric redshift (photo-zs) derived by D20. Only one source (ILTJ160658.74+550607.0) has a photo-z z > 6, whereas estimates for other sources place them at 1.4 < z < 2.1, suggesting that they are likely to be low-redshift interlopers. It is furthermore noticeable that z−NB921 ≥ 1.0 within 1σ in the LoTSS photometry, but it is not as strong as the HSC values. This is likely due to the difference in underlying photometry used, as well as the large uncertainties in z magnitude, which could partly be due to the chosen 3″ aperture size. To further investigate the nature of these five sources and their probable redshifts, we examined their spectral energy distributions (SEDs). In SED fitting, we also included the FIR data, which is not included in the photo-z estimates.

3.1 SED fitting

We performed SED fitting on the LOFAR-detected sample using the code Bayesian Analysis of Galaxies for Physical Inference and Parameter EStimation³ (or BAGPIPES; Carnall et al. 2018). In summary, BAGPIPES is a Bayesian spectral fitting code developed for generating complex model galaxy spectra from photometric and spectroscopic observations using the MultiNest nested sampling algorithm (Feroz et al. 2019). The input observations can range from far-UV to microwave wavelengths. The model parameters used for fitting are the star-formation history (SFH), the nebular component, and the dust attenuation and emission components. An overview of the parameter values used in this work is given in Table A.1, and the available input photometric data are summarised in Table 1. The SFH is modelled using a double power law with parameters for the total stellar mass formed and the metallicity, and the dust attenuation model of Calzetti et al. (2000) is used. A more extensive description of BAGPIPES and its performance is available in Carnall et al. (2018).

Besides Lyα, other strong nebular emission lines that could account for a NB excess in the NB921 filter are Hα, [OIII], Hβ, and [OII] nebular emission lines at z~ 0.40, 0.84, 0.89, and 1.47, respectively. The photo-z accuracy is not high enough to be able todistinguish the Hβ emission line from the [OIII] emission line, so throughout this work we refer to it as the [OIII] emission line at z ~ 0.84 since this line is more prevalent in emission line galaxies (e.g. Hayashi et al. 2018; Sobral et al. 2015b; Khostovan et al. 2015). The combination of the implausibly high dust extinction required to produce such red SEDs at z = 0.4 and the proportionally low co-moving volume probed means that the contamination of this sample by Hα emitters at z = 0.4 is thought to be negligible. Therefore, specific models for [OIII], [OII], and Lyα emission lines are fitted to the LAE candidates.

It must be noted that BAGPIPES does not include AGN templates, which could lead to an underestimation of the mid-IR flux. An AGN contribution can be therefore be identified by the inability of BAGPIPES to fit the mid-IR flux. The possible effects of the lack of AGN templates are therefore discussed together with the results.

To investigate the probability of a source residing at the above-mentioned redshifts, the SED fits are performed with fixed z = 0.84, 1.47, 6.57 as well as at 80 redshifts in the range z = [0.1, 8], in steps of 0.1. For each redshift, the SED fit resulting in the minimum reduced chi-squared value () is determined, to be used as a comparison for the different models. Additionally, non-detections, with signal-to-noise ratio (S∕N) < 3.0, are used by BAGPIPES and in the calculation. However, negative fluxes are not taken into account. The SED fitting results of the high-z LAE candidate ILTJ160658.74+550607.0 are shown in Fig. 3. The value lies within a small range and is lowest at z ~ 4.9 when fitting with the original errors. However, in panchromatic SED fittings, there are many model uncertainties and the flux errors arising from flux calibration issues or correlated noise are usually underestimated (see e.g. Marquez et al. 2014). We therefore re-fitted the LOFAR-detected sample, incorporating an additional 10% flux uncertainty in quadrature to assess what effect this may have on the preferred redshift solution (see Fig. 3). The overall shape of the curve remains largely the same with the additional 10% flux uncertainty, including the decrease around z ~ 4−5, though it does improve the fit at z = 1.47, making it comparable within the errors to the z ~ 4−5 fit. However, around this redshift of z ~ 4−5, there is no emission line with a typical equivalent width (EW) high enough to be detected in the HSC NB921 filter. Also, because of the large number of data points, the relative contribution of the NB measurement to the overall χ² value is limited. The decrease in χ² can therefore be attributed to the BB SED shape. In the right-hand panel of Fig. 3, the best-fitting posterior median models at z = 0.84, 1.47, and 6.57 for [OIII], [OII], and Lyα, respectively, are shown in yellow. Here, the SED fit at z = 1.47 is the preferred solution. The SED fitting results of the four other sources in the LOFAR-detected sample are shown in Fig. B.1. The χ² values are minimised around z ~ 1.5, providing strong evidence that these are [OII] emitters at z = 1.47.

To compare the preference for each model, we used the Bayesian information criterion (BIC), which in its general form is given by (3)

where L is the maximised likelihood, k the number of parameters, and N the number of observations (see Kass & Raftery 1995). In the case of Gaussian distributed model errors, this becomes (4)

The model with the highest probability minimises this value, and, when comparing models, a ΔBIC > 2 gives positive evidence for one of the models being preferred over the other, whereas −2 < ΔBIC < 2 indicates that there is no evidence for either model being preferred. A ΔBIC > 6 gives strong positive evidence (as defined by Kass & Raftery 1995). When comparing our SED fits, the number of observations and parameters remains constant; therefore, BIC becomes χ² and we can use the Δχ² values of the fits at z = 0.84, 1.47, and 6.57 to determine the preference for each model. The and Δχ² values of the Lyα, [OII], and [OIII] emission line fits to the LOFAR-selected sample are summarised in Table 3. Here, the original errors are used to compare the models since the 10% added flux error causes over-fitting of the data, resulting in values below 1. For all five sources, Δχ²(Lyα −[OII]) and/or Δχ²(Lyα − [OIII]) > 2, meaning the models for [OII] or [OIII] provide a better fit to the data and Lyα is likely not to be the detected emission line. In addition, Δχ²([OIII] −[OII]) > 2 for all sources, suggesting the detected emission line is [OII], in turn implying that the sources are situated at z = 1.47. This conclusion is also supported by comparing the co-moving volume being probed at z = 0.84 and z = 1.47, which is ~2 times higher for z = 1.47 when assuming a flat Λ-CDM cosmology. The evidence for Lyα not being the source of emission is less strong in the case of ILTJ160658.74+550607.0 (shown in Fig. 3); however, the FIR SED shape, consisting of SPIRE and PACS observations, is best reproduced by the z = 1.47 model and therefore seems to be the actual redshift. Follow-up observations are needed to confirm the redshift of this source. If this source is a radio galaxy at z = 6.6, it would be the most distant radio galaxy known to date.

We present the derived stellar mass and star-formation rate for each fixed-redshift SED fit in Table 4. When redshift is fixed at z = 6.57, the derived stellar masses are 10^11.7 < M_* < 10^11.9 M_⊙ and SFRs are ~650−4400 M_⊙ yr⁻¹. Even given the large uncertainties on the high mass end of the stellar mass function at these redshifts, the likelihood of masses being this high is small (see e.g. Duncan et al. 2014; Grazian et al. 2015; Song et al. 2016). The median Lyα SFR obtained by, for example, Calhau et al. (2020) is 9.8 M_⊙ yr⁻¹ for their sample at z > 3.5, when excluding AGN. The SFRs derived from their FIR-detected LAEs (of which 76% are X-ray or radio AGN) using their IR fluxes yield higher values of 200 M_⊙ yr⁻¹, with a few sources >600 M_⊙ yr⁻¹.

As noted before, BAGPIPES does not include AGN templates. If there were a strong AGN contribution, there should be a warm dust contribution at mid-IR wavelengths. However, none of these five sources show excess in the MIPS 24 μm band. In addition, all these sources are classified as star-forming from both their optical and radio emissions via the hybrid SED fitting method from Best et al. (in prep.). Furthermore, in the photo-z fitting procedure of D20, AGN templates are used, and these results are consistent with the redshifts obtained from SED fitting. Unlike in D20, in the BAGPIPES modelling FIR fluxes have been included in the fit, which, for the source ILTJ160658.74+550607.0, helps break the degeneracy between the z = 1.5 and 6.6 solutions. The remaining small differences in the resulting redshifts between this work and the photo-zs of D20 can be attributed to the use of a different set of templates (see Duncan et al. 2018) and the machine-learning contribution in D20.

Fig. 2 Multi-wavelength cutouts of the LOFAR-detected sample including, from left to right: HSC r band, HSC z band, HSC NB921 band, UKIDSS K band, IRAC 4.5 μm, and LOFAR 150 MHz. The sources show negligible or no emission below the Lyman break in the r band; however, clear emission is seen in both the UKIDSS K band and IRAC.

Fig. 3 SED fitting results for ILTJ160658.74+550607.0, obtained using  BAGPIPES. Left: reduced χ2 for SED fits in the range 0 < z < 8, with locations of the Lyα, [OII], and [OIII]+Hβ emission lines indicated by dashed black lines. The dashed grey lines indicate the redshift range probed by the NB filter for the strong emission lines [OIII], [OII], as well as Lyα. Right: observed flux of ILTJ160658.74+550607.0 (blue) in mJy with the original errors and posterior median model (yellow) for fits with redshifts fixed at z = 0.84, 1.47, and 6.57. The source is best fitted by a template with z = 1.47 when considering the possible emission lines. No other emission line is strong enough for the source to be at z ~ 4.

3.2 Emission line and radio properties

From the photometric data, we determined the physical properties of the LOFAR-selected sample, such as radio luminosity, line luminosity, and EW (see Table 5). As each of these quantities depend on the redshift of the source, the values are given for fixed z = 0.84, 1.47, and 6.57. A spectral index α of −0.7 ± 0.7 (defined by S_ν ∝ ν^α with S_ν the flux density in Jy) was assumed when determining the radio luminosity, where the error is based on the spectral index distribution derived by Saxena et al. (2018b). No spectral index could be obtained for the sources in the LOFAR-selected sample since no sufficiently deep radio data (e.g. from the Very Large Array or the Giant Metrewave Radio Telescope) was publicly available. For the calculation of the EW and line luminosity, we used the equations presented in Sobral et al. (2012). Subsequently, the line luminosity and radio luminosity were calculated using the line flux and radio flux, respectively. The rest-frame EW (EW₀) was determined by dividing the observed EW by 1 + z to correct for the Hubble expansion. The values obtained for each assumed redshift are shown in Table 5.

A large range of Lyα rest EW values are quoted in the literature. Lyα EW values of, for example, ~0−400 (Ono et al. 2010) and 10−1000 (Leclercq et al. 2017) have been found. Furthermore, the Lyα LFs presented in Calhau et al. (2020) and Konno et al. (2018) show Lyα luminosities ranging ~10^42.5−10^44.0 erg s⁻¹. A global study of LAEs at 2.5 < z < 5.8 by Sobral et al. (2018a) in the COSMOS field resulted in a characteristic luminosity L_{Ly α} ~ 10^{43.8 ± 0.1} erg s⁻¹ of the Schechter function. The values for the EW and L_Lyα found in this work are consistent with these literature results.

The EW₀ derived when we assumed that the detected lines were [OII] and [OIII] are comparable to the EW distribution found in Khostovan et al. (2016), where the majority of the [OIII]+Hβ line emitters at z = 0.84 have 20 < EW₀ < 500, and the majority of [OII] line emitters at z = 1.47 have 20 < EW₀ < 300. The derived [OII] and [OIII] luminosities are within the ranges of 41.4 < log₁₀(L_[OIII]) < 42.3 and 40.5 < log₁₀(L_[OIII]) < 42.1. However, it must be noted that the sources in this study could be different from the general population of [OII] and [OIII] emitters since this sub-sample of sources satisfies the LAE selection criteria. These criteria (see Eqs. (1) and (2)) select for sources with high magnitude differences between the z band and NBs, as well as non-detections in filters below the supposed Lyman break, suggesting that these are potentially sources with strong Balmer breaks.

The obtained radio luminosities and SFRs from BAGPIPES can be compared to the L_150MHz-SFR relation derived by Gürkan et al. (2018), assuming the relation does not evolve with redshift. This assumption is supported by the study of Duncan et al. (2020), where no redshift evolution is found out to z ~2.6. If the sources are located at z = 1.5, all five sources are slightly offset to higher radio luminosities than would be expected from the L_150MHz-SFR relation, with a difference of 0.04–0.69 dex from the 1σ upper limit, which indicates that they could be low luminosity AGN. However, these deviations are non-significant when taking the 0.3 dex error on the radio luminosities and the ~0.3 dex scatter of observations in the L_150MHz-SFR relation into account. Moreover, the SED fitting did not suggest an AGN contribution. The Gürkan et al. (2018) relation is consistent with the new derived L_150MHz-SFR relation from the LOFAR Deep Fields data by Smith et al. (2021).

4 Radio and optical-IR properties of the wider LAE population

The available evidence suggests that the LAE candidates in the LOFAR-selected sample are not located at z ~ 6.6. Therefore, we cannot give an upper limit on the AGN fraction of LAEs at z > 6. These results, however, highlight different possible problems of the impurity of LAE samples as well as the low reliability of measurements on the LAE properties and the claimed AGN fractions at these redshifts. To investigate the scale of this problem, we further analysed the wider z ~ 5.7 and z ~ 6.6 LAE samples using our full panchromatic photometry.

The compiled optical and IR photometric catalogue for ELAIS-N1 also allows for the study of the multi-wavelength properties of SILVERRUSH sources that do not have radio detections. The sources from the LoTSS multi-wavelength catalogue and the SILVERRUSH catalogues are cross-matched within 1″ and inspected by eye. Of the 349 and 229 LAEs in ELAIS-N1 in the SILVERRUSH NB921 and NB816 LAE catalogues, 58 and 53 sources, respectively, are detected in the LoTSS multi-wavelength catalogue. This sample will from now on be referred to as the “optically selected sample”. The LoTSS catalogue, selected from the χ² detection image, may not be as optimised for Lyα candidates as the SILVERRUSH extraction; furthermore, SILVERRUSH used the deeper internal intermediate data release HSC-SSP 16A, so faint sources (i mag > 26) are often missed in the LoTSS catalogue. We performed forced photometry on all SILVERRUSH sources to compare the two samples. The K and NB magnitudes from the optically selected and non-detected SILVERRUSH sample are shown in Fig. 4. The non-detected SILVERRUSH sources with a flux measurement at > 1σ have fainter K band magnitudes (K = 23.0/23.2 for NB816/NB921) compared to the optically selected sources (K = 22.1/22.3 for NB816/NB921). The median NB magnitudes are only slightly fainter for the non-detected SILVERRUSH sources (NB = 24.2/24.6 for NB816/NB921) compared to the optically selected sources (NB = 23.9/24.4 for NB816/NB921). This indicates that the optically selected sample contains sources that are brighter and redder than the non-detected SILVERRUSH sample. We note that these magnitude differences are estimates because of the faintness of the sources and the choice of a 1σ detection limit. The bias of our sample is further discussed in Sect. 5.2.

Properties of the optically selected sample obtained from photometry and SED model fitting are discussed in Sect. 4.1. Furthermore, we stack the radio and the IR and optical images of the full sample to determine the limiting magnitudes, which are presented in Sect. 4.2.

Fig. 4 K and NB magnitudes (NB816 or NB921) from the optically selected sample and the non-detected SILVERRUSH sample. Limiting magnitudes (1σ) are given for sources below 1σ detection (filled and unfilled triangles for sources in the optically selected and non-detected SILVERRUSH samples, respectively). The histograms only include sources with a >1σ detection. The median K and NB magnitudes are indicated by the solid and dashed lines for the optically selected and non-detected samples, respectively. In the optically selected sample, 23% of sources are not detected in K above 1σ, while for the non-detected SILVERRUSH sample, this fraction is 82%. The optically selected sample is brighter and redder than the original SILVERRUSH sample.

4.1 SED fitting optically selected sample

The stacked full photometric redshift probabilities from D20 for the optically selected sample are shown in the left-hand panel of Fig. 5. This is a sum of the photo-zs of the individual sources. Of the two probability distributions, 44% of the posterior is at z > 5 for the NB816 sources, and less than 10% of the posterior is at z > 6 for the NB921 sources. This already suggests a high contamination fraction in both of these samples.

These photo-zs could potentially be biased towards lower redshifts due to the applied magnitude prior, which favours low-z solutions (z < 2) even for faint sourceswhile not disallowing high-z solutions (Duncan et al. 2018). We therefore extend this with more detailed SED fitting by repeating the procedure in Sect. 3.1 for the full optically selected sample of LAEs by fixing the redshifts at z = 0.63, 1.19, 5.72 for the NB816 catalogue and z = 0.84, 1.47, 6.57 for the NB921 catalogue. Of the 58 candidate LAEs at z = 5.7, nine have FIR detections. For the z = 6.6 candidates, six out of 53 have FIR detections. The other sources have only UV-NIR data. We fitted all available data for each source. For each source, we identified the strongest emission line with the best fitting SED (lowest ) and assigned thesources the corresponding redshift. We present the numbers of sources at each of the redshifts in Fig. 5(right-handpanel). For five sources in NB921 and one source in the NB816 LAE catalogue, the SED fits are poorly constrained due to detections in fewer than nine photometric bands, hence these results are excluded from further analysis. The best SED fits (right-hand panel in Fig. 5) suggest contamination rates of 90% and 93% for the NB816 and NB921 catalogue matches, respectively. To quantify this further, we repeated the model comparison method for this larger non-radio-detected sample, again with fixed redshifts. The resulting Δχ² values for Δχ² (Lyα − [OII]), Δχ² (Lyα−[OIII]), and Δχ² ([OIII]−[OII]) are presented in Fig. 6, where the dotted lines indicate Δχ² values of 2 and −2. Sources likely to be Lyα emitters at z = 6.57/5.72 satisfy Δχ² < −2 for either Δχ²(Lyα−[OII]) or Δχ²(Lyα− [OIII]), under the condition that the other Δχ²  2. This is the case for four sources in the NB816 and five sources in the NB921 catalogue. On the other hand, 43 and 47 of the sources have Δχ² > 2 for either of the two Δχ² values, for both NB816 and NB921, so these are more likely to be z = 0.63/0.84 or z = 1.19/1.47 contaminants, respectively, rather than z = 5.72/6.57 LAEs. The model comparison for sources with Δχ² between −2 and 2 is inconclusive and can be regarded as an estimate of the uncertainty. The contamination rate determined from SED fittings using the three emission lines [OII], [OIII], and Lyα is therefore 81–92% (43–49 out of the 53 sources) and 81–91% (47–53 out of the 58 sources) for NB816 and NB921, respectively. Similarly, comparing the χ² of [OII] and [OIII], [OII] is likely (Δχ² > 2) to be the emission line for 63 ± 31% and 55 ± 28% of the possible contaminating sources for NB816 and NB921, respectively, whereas [OIII] is likely to be the corresponding emission line for 6 ± 31% and 17 ± 28% of the sources for NB816 and NB921, respectively (with the remainder inconclusive). Comparing these results to the primary and secondary photo-zs of D20, we obtain that, within the error, 21/53 (40%) of the photo-zs in the NB816 sample and 41 of the 58 (71%) photo-zs in the NB921 sample are consistent. The remaining discrepancy between the two photo-z estimates can be attributed to the faint magnitudes of the sources and the high redshifts, which significantly increase the photo-z uncertainties and result in a broad photo-z posterior (see D20). Comparing our contamination rates to the stacked photo-z probability distribution in the left-hand panel of Fig. 5 shows that the photo-z estimates are consistent for the NB921 sample, but contamination is estimated to be lower for the NB816 sample. The photometric redshifts of D20 and the SED fitting results together therefore suggest the contamination of our subset is 56–92% and 81–91% for the NB816 and NB921 samples, respectively.

Physical properties of the samples, derived from SED fitting (SFR and stellar mass) and from the photometric data (line luminosity and EW) are shown in Fig. 7 for each of the three proposed redshifts. Assuming these sources are LAEs, the median Lyα luminosity is 10^{43.1 ± 0.2} and 10 erg s⁻¹ for NB816 and NB921, in line with other works from, for example, Calhau et al. (2020) and Konno et al. (2018). Fixing redshifts at z = 1.47/1.19 and z = 0.84/0.63 gives [OII] and [OIII] line luminosities that are also in good agreement with Khostovan et al. (2016) and Hayashi et al. (2018), as described in Sect. 3.2.

The EW₀ and SFR values derived for all three redshifts lie within the range of reported values for other samples in the literature. However, these physical parameters vary substantially both in this work and the literature, which is likely due to a wide range of galaxy types and survey selection functions (see e.g. Khostovan et al. 2019; Calhau et al. 2020). If we assume the candidates to be low-z interlopers, we are again studying a subset of sources of the [OII] and [OIII] emitter population, which is likely biased, and thus these sources do not necessarily have the same physical properties. We therefore refrain from making any further comparisons between the obtained EW₀ and SFRs in an attempt to support the evidence for the low-z nature of a high fraction of these sources. Remarkable, however, are the stellar masses for fixed z = 6.57, which are extremely high (M_* > 10¹¹ M_⊙), especially for the sources in the NB921 catalogue. The majority of these sources are therefore unlikely to be LAEs.

The derived source properties are also subject to systematic effects from the SED fitting models used and the possible biases introduced by the sample selection. BAGPIPES is constrained by the set input parameters, prior distributions, and model dependencies. Despite the limitations, the models can often successfully recreate the observed SEDs. However, spectroscopic observations are necessary to confirm the redshifts and derived physical properties.

Fig. 5 Left: stacked photometric redshift probability distribution (sum of individual sources) of the optically cross-matched sample of the NB921 and NB816 LAE catalogues from D20. Right: histograms of the best-fitting solutions for SED fits performed with fixed-redshift solutions corresponding to the three candidate emission lines. The known possible redshifts from the NB emission line detections and more detailed SED fitting allows us to exclude other high-z solutions that are sometimes favoured by the photometric redshift code.

Fig. 6 Comparison of χ2 values obtained from SED fitting models of Lyα, [OII], and [OIII] for the optical cross-matched LAE sample of NB921 and NB816. The dashed lines indicate Δχ2 values of 2 and − 2. For values higher than 2 and lower than −2, there is a preferred model, but in between the dashed lines there is no preference for either model. The bin width is set to 2.

4.2 Radio and IR-optical stacking

To investigate the radio properties of the LAE candidates not detected by LOFAR, first we stacked the LOFAR images to set observational constraints on the average radio properties of the samples. SILVERRUSH sources outside the area covered by the multi-wavelength LoTSS data (see Sect. 2.3) are removed from the catalogue (22 and three sources for the NB921 and NB816 catalogues, respectively). The five sources in the LOFAR-detected sample in the NB921 catalogue are also removed from the stacking data set (see Sect. 3). For both NB816 and NB921, first the full sample was stacked (containing 226 and 322 sources, respectively) and then the optically selected samples (containing 53 and 49 sources) were stacked separately; these results are shown in the top and middle panels of Fig. 8. We assumed the sources to be unresolved in the 6″ resolution radio images, and could therefore obtain the median flux from the peak pixels. The average rms noise level for individual radio images is ~25 μJy beam⁻¹. The stackedradio flux density, rms, and resulting S/N obtained in each sample are summarised in Table 6. None of the stacked radio flux densities are significant (for all S∕N < 2). Therefore, we can place a 3σ upper limit on the radio flux density of ~4.8 μJy for the sources in the full sample and ~12.0 μJy for sources in the optically selected sample.

Furthermore, we stacked the four and five optically selected sources from NB816 and NB921, respectively, that seemed to be likely LAEs from SED fitting together with the 173 and 274 sources in the non-matched SILVERRUSH samples, under the assumption that they are LAEs. The non-matched sources are fainter and bluer and therefore less likely to be low-redshift contaminants (see Fig. 4 and Sect. 5.2). The stacks again yield no detection with a limiting radio flux density of 5.7 and 4.8 μJy (3σ) for NB816 and NB921, respectively (see the bottom panels of Fig. 8). When assuming their redshifts to be z = 5.7 and 6.6, the radio flux density can be converted to 2σ radio luminosity upper limits of 8.2 × 10²³ and 8.7 × 10²³ W Hz⁻¹, respectively. Using the low-z L_150MHz-SFR relation from Gürkan et al. (2018), an estimated 2σ upper limit on the SFR can be placed on these LAEs. Our derived upper limits are ~53 and 56 M_⊙ yr⁻¹ on the z = 5.7 and z = 6.6 LAEs, respectively. The currently known relations between L_150MHz and the SFR are derived from low-z galaxies; therefore, this is our current best SFR estimate. However, the work of Smith et al. (2021) suggests that there is no strong redshift evolution. Similar SFRs are obtained with the L_150MHz-SFR relations from Wang et al. (2019), with differences of a few solar masses per year. This upper limit is higher than the expected SFR of LAEs (see e.g. Calhau et al. 2020), so, to be able to place stronger constraints on the SFR, a larger sample of LAEs, or deeper observations, are necessary.

Next, we stacked the optical bands, ranging from u to y, and find no detection in the stacked full sample, with SpARCS u, HSC g, and HSC r band 3σ limiting magnitudes of ~27.6, 25.1, and 24.5, respectively (using a 3″ aperture), indicating that this population of likely z ~ 6 SILVERRUSH galaxies without radio detection is consistent with z = 6.6 down to this magnitude limit, as expected from the initial selection method. These non-detections, however, donot rule out the possibility of them being low-z interlopers. We find similar results when we stack optical images for the optically selected sample, with measurements of u, g, and r down to 3σ limiting magnitudes of 26.8, 24.2, and 23.7, respectively. In contrast, SED fitting on the full panchromatic photometry in the optically selected sample suggested a large contamination rate (see Sect. 4.1). These sources in the optically selected sample are typically faint in the optical and are primarily detected in the χ² image because of their high IR flux, indicating that IR observations are critical for identifying low-redshift interlopers.

Finally, we investigated the NIR and mid-IR data for the full LAE candidate sample (322 and 226 sources in NB921 and NB816, respectively) by stacking the J and K images (1.25 and 2.21 μm) obtained by UKIDSS, as well as the mid-IR (3.6, 4.5, 5.8, and 8.0 μm) images obtained by Spitzer-SWIRE. We find a detection (S∕N > 3) in all, except for the 5.8 μm channel, with only a weak detection (S∕N ~ 3) in the 8 μm channel; this, again, indicates an absence of strong AGN contribution.

Fig. 7 Derived line luminosity, EW, SFR, and stellar mass of the optical cross-matched LAE sample with the NB921 and NB816 LAE catalogues at different fixed redshifts of z = 0.63/0.84, 1.19/1.47,and 5.72/6.57. The dashed blue and orange lines indicate the median for NB921 and NB816, respectively; the values given in the legend include errors derived from the central 68% range of the distribution. The lowest bin in the SFR histograms show the sources with SFR < 10−1 M⊙ yr−1.

5 Discussion

In this section, we discuss the LOFAR detection rate (Sect. 5.1), the LAE sample selection, and possible biases of and improvements to the selection process (Sect. 5.2).

5.1 LOFAR detection rate of SILVERRUSH LAEs

No z > 6 radio galaxy candidates have been found in the cross-matched LoTSS and SILVERRUSH survey. In this section, we discuss this in light of predictions made for the number densities of such sources. The expected number of radio galaxies in LoTSS in the redshift slices 5.72 ± 0.05 and 6.57 ± 0.05, from models presented by Saxena et al. (2017), are 5.3 ± 0.4 and 1.5 ± 0.1 sources per deg², respectively.From these models, we would expect about nine radio galaxies with 6.52 < z < 6.62 to be detected in the ~6 deg² overlap region between SILVERRUSH and LoTSS. However, the LAE fraction in brighter galaxies at redshifts 6 < z < 8 is shown to be only 10–20% in the study by Schenker et al. (2012), and the distribution of Lyα luminosities in high-z radio galaxies (z > 3) is poorly constrained due to small sample sizes (see e.g. Saxena et al. 2019). Assuming this LAE fraction, only ~1 out of the nine predicted sources is expected to be selected by SILVERRUSH and detected by LOFAR. Given this expectation of only one detection, it is not surprising that we do not find any radio emitting LAEs at z ~ 6. Furthermore, there are still many uncertainties in the expected number of radio detections, especially at high redshift. The LoTSS Deep Field observations probe the faint end of the radio luminosity function (RLF) due to the small volume and high sensitivity, which is poorly constrained by the models (see Saxena et al. 2017). Therefore, the uncertainties in the RLF are critical for the number density predictions in this paper. These points together provide a plausible explanation for the absence of z > 6 radio galaxies in the studied sample. Based on these numbers, a larger area of overlapping coverage would be necessary to find many z > 6 radio galaxies.

Fig. 8 Stacked LOFAR cutouts of LAE candidate samples. Left: NB816 full sample (226 sources; upper left), optically selected sample (53 sources; middle left), and likely Lyα emitters (177 sources; bottom left). Right: NB921 full sample (322 sources; upper right), optically selected sample (49 sources; middle right), and likely Lyα emitters (279 sources; bottom right). The LOFAR-detected sources and sources not overlapping with LOFAR coverage have been removed.

5.2 LAE selection methods

It is possible that, with the optical and IR selected samples, we are probing galaxies that are more likely to be low-redshift interlopers than the SILVERRUSH samples that do not appear in our optical and IR catalogue. Very dusty low-redshift galaxies could, for instance, have a Balmer break at 4000, mimicking the Lyman break (see e.g. Matthee et al. 2014, 2017). These galaxies often have red J− K colours, so constraining J−K to be flat or blue (J−K ≤ 0) prevents contamination from these dusty low-redshift interlopers. The J − K colours of the optically selected samples and the LOFAR-detected sample are shown in Fig. 9. For both the NB816 and NB921 optically selected samples, eight sources (out of the 53 and 58 sources, respectively) do not have J or K magnitudes. Only three sources in NB921 satisfy J− K ≤ 0, and these three all have high photo-z values (see Fig. 9). No sources within the LOFAR-detected sample satisfy this criterion. The effect is even more apparent for NB816, where the diagonal line divides the low photo-z candidates from the high photo-z ones, except for one extremely low-z source. Therefore, NIR observations are of great value in identifying the low-z interlopers. However, we note thata strict J−K ≤ 0 cut could cause red LAEs to be rejected from the sample.

Shibuya et al. (2018b) have shown that the contamination rate depends on the NB magnitude, where populations with lower magnitudes have higher contamination rates. The sources matched with our optical catalogue are brighter and likely redder on average than the full LAE sample due to the shallower HSC observations and the source extraction from the deep χ² image in the optical catalogue. The source extraction is based on significant detection in the χ² stack of either optical and NIR or IRAC Channel 1 and 2 images. The SILVERRUSH sources that are missing from our optical sample will therefore be the ones that are fainter in the optical bands and not red enough to be detected in the NIR or mid-IR. The χ² detection image did not include any y-band data; this could significantly contribute to whether or not the real LAEs are detected since, at 6 < z < 7, sources are not detected blue-wards of the y band where the Lyman break is situated. Therefore, the non-matched sources are more likely LAEs. If we assume that all the non-matched sources are true LAEs, then we obtain contamination rate lower limits of 23 and 20% for the full NB816 and NB921 catalogues, respectively, which would agree with the quoted 30% contamination rate from the spectroscopic follow-up by SILVERRUSH (see Shibuya et al. 2018b).

To highlight the difference between the likely LAEs and the contaminants in the NB921 catalogue, we show the stacked SEDs in Fig. 10. The SED of the stacked assumed LAE sample is consistent with z = 6.6, whereas the SED of the contaminant sample is consistent with z = 1.47. For LAEs at z = 6.6, the IRAC 3.6 μm filter corresponds to the [OIII]+Hβ emission line (see Smit et al. 2015). Since the high EWs of the Lyα and [OIII]+Hβ lines might not be represented sufficiently well by the nebular emission line models in BAGPIPES, the NB921 filter and IRAC 3.6 μm filter have not been included in the fits. The J−K colours are −0.02 ± 1.05 and 1.27 ± 0.11 for the stacked assumed LAEs and contaminants, respectively. This again highlights the advantage of selecting galaxies for which the NIR slope satisfies J− K ≤ 0 when searching for high-redshift LAEs. From the BAGPIPES models, we obtained stellar masses of 10^{9.6 ± 0.2} and 10^{10.5 ± 0.1} M_⊙ for the SED fit of stacked LAEs and stacked contaminants respectively, and SFRs of 8.5 and <0.4 M_⊙ yr⁻¹, where the upper limit is the 84th percentile of the SFR posterior. This is consistent with the upper limit on the LAE SFR obtained from the corresponding radio stacks (see Sect. 4.2).

As we concluded in Sect. 4.1, the majority of the LAE candidates that have a match in our optical catalogues are likely to be low-z interlopers of [OII] and [OIII] emitters, or possibly Hα at z = 0.40. If the sample is indeed contaminated by a large fraction of lower-redshift line emitters, other NB studies might have also been affected by large fractions of contaminants, possibly leading to a misrepresentation of the physical properties of LAEs.

To minimise the effects of contamination, one option for future studies is to consider the optical colour selections used to select LAE candidates. In the SILVERRUSH studies, a forced-aperture catalogue is used for statistical studies (see Shibuya et al. 2018a). This catalogue is generated with apertures at fixed positions, with a typical 3σ detection requirement that is lowered to 2σ for the g and r bands. Colour selections are then imposed to select LAEs. The selections made are z−NB921 > 1.8 for z = 6.6 LAEs, which corresponds to EW_0,Lyα > 14 for the z ~6.6 sample. We now compare the numbers that would be derived from the forced SILVERRUSH catalogue, as opposed to the combined forced and unforced catalogue used in the rest of this study. The forced SILVERRUSH catalogue contains 48 sources in NB921 and 130 in NB816. The number of sources in the forced catalogues compared to the unforced catalogue is substantially lower; therefore only two LOFAR-detected sources are found in NB921, and no sources are found in NB816. These two matched sources are ILTJ160227.18+544759.8 and ILTJ160930.42+544435.9, which are also likely z = 1.5 [OII] emitters (see the top two panels of Fig. B.1). From the optical sample, there are three and 11 in the forced catalogues of NB921 and NB816, respectively. All three sources in NB921 have photo-zs of ~1.7 ± 0.7, suggesting that these are low-redshift interlopers. However, of the 11 sources in the NB816 forced catalogue, ten have a photo-z in the range 4.5-5.8 ± 1.5 (the one remaining source has a photo-z of 0.5 ± 0.9), suggesting that the LAE fraction is likely high. This implies that, at least in the case of z = 5.7 LAEs, the contamination rate could be significantly lowered by making stricter cuts.

Many LAE NB surveys apply similar cuts for selection, but differ in their strictness on the cuts and the availability of multiple bands (e.g. Ouchi et al. 2008; Santos et al. 2016; Matthee et al. 2015). Other studies have also shown that the majority of LAEs have EW₀ > 50 (e.g. Ouchi et al. 2008). However, Sobral et al. (2018a) also show, by using a large compilation of spectroscopic redshifts, that, even though their sample of high EW candidate line-emitters is dominated by likely LAEs, there is still a significant population of Hα, [OIII]+Hβ, and [OII] emitters. Since the observed EW₀ distribution in these lines is similar to Lyα (see e.g. Khostovan et al. 2016; Hayashi et al. 2018), it can be challenging to clean the samples further using just NB and BB optical surveys. The contamination rates we obtain are comparable to those found by Matthee et al. (2014), who find a 90% contamination in LAE samples selected on the Lyman break criteria. However, their survey targeted higher redshifts of z > 7, so contamination is expected to be higher.

This work further highlights the fact that ancillary IR observations are critical for the LAE sample selection. Although the J to 4.5 μm data in ELAIS-N1 are significantly shallower than those available in smaller fields, such as COSMOS, they can still play a important role in minimising interlopers in wide field LAE samples. Decreasing the contamination rate is critical for improving the constraints on the Lyα LF and clustering measurements, as well as for reducing the time and expense required for spectroscopic follow-up observations.

Fig. 9 UKIDSS J and K band magnitude of the optically selected sample and LOFAR-detected sample, with colour-coded photo-zs from D20. The diagonal (dashed) line indicates the J−K ≤ 0 criterion. The sources above this line satisfy this criterion and also have notably higher photo-z values.

6 Summary

In this paper, we have studied the LOFAR and optical and IR properties of SILVERRUSH LAEs in ELAIS-N1 using the latest LoTSS Deep Fields observations and the cross-matched multi-wavelength catalogue that ranges from UV to FIR. The publicly available SILVERRUSH catalogues of ELAIS-N1 consist of 229 and 349 LAE candidates at z = 5.7 and z = 6.6, respectively, identified with the HSC NB816 and NB921 narrowbands. We found five LOFAR-detected LAEs in the NB921 SILVERRUSH catalogue. These five sources have LOFAR radio fluxes in the range of 0.09–0.35 mJy. We performed SED fitting on the radio-detected sources using BAGPIPES and find that all five sources are likely z = 1.5 [OII] emitters, based on χ² analysis and derived physical properties. We therefore do not find any promising z > 6 radio galaxycandidates and cannot put any constraints on the AGN fraction of LAE at high redshift.

In light of the high contamination rate of radio-selected SILVERRUSH sources, we performed the same SED fitting procedure on the wider LAE population by cross-matching the SILVERRUSH NB816 and NB921 catalogues with the LoTSS Deep Field multi-wavelength catalogue. We find 53 and 58 sources with matches in the LoTSS multi-wavelength catalogue, out of the 229 and 349 SILVERRUSH LAE candidates, respectively. The resulting Δχ² values suggest veryhigh contamination rates of 81% (92%) and 81% (91%) for the NB816 and NB921 catalogue sample, respectively, where the values in the parentheses include the sources that are also likely contaminants but not statistically significant. This is consistent with the photo-zs from D20 given the large photo-z uncertainties at these high redshifts and faint magnitudes. Of all the contaminated sources, 63 ± 31% are likely to be [OII] and 6 ± 31% are likely to be [OIII] emitters for NB816, and similarly 55 ± 28% are likely to be [OII] and 17 ± 28% [OIII] emitters for NB921 (with the remainder inconclusive). However, the non-matched sources are fainter and bluer and therefore more likely to be true LAEs, which would be in line with the 30% contamination rate obtained from spectroscopic follow up by SILVERRUSH (Shibuya et al. 2018b).

Most importantly, in agreement with Matthee et al. (2014, 2017), we find that a J− K ≤ 0 cut can identify a large fraction of red and dusty low-z interlopers, which are not removed by the standard LAE criteria due to the low optical S/Ns and the Balmer break mimicking the Lyman break. These results highlight the need for ancillary IR observations in LAE sample selection to minimise the amount of interlopers.

We stacked radio image cutouts in the positions of all LAE candidates and find no significant detection with LOFAR when removing the radio-detected sources. Removing the likely contaminating sources and stacking the LOFAR data for the four and five remaining matched sources in NB816 and NB921 together with the 173 and 274 non-matched sources yields a 2σ upper limit on the SFR of ~53 and 56 M_⊙ yr⁻¹, respectively. In the optical stacked u, g, and r measurements, we do not find a detection down to ~27.6, 25.1, and 24.5 magnitudes, respectively, which is consistent with the Lyman break at z = 6.6. The stacked IRAC SWIRE 3.6 and 4.5 μm and UKIDSS J and K images do yield a detection (S∕N > 3), but there is no significant detection (S∕N ~ 2) in the IRAC SWIRE 5.8 and 8.0 μm measurements.

The current 180 h of observations of ELAIS-N1 are expected to reach 500 h within the next two years. This will result in factor of ~2 deeper observations down to 11 μJy beam⁻¹, enabling the study of the fainter and higher-redshift radio population. Furthermore, the WEAVE-LOFAR spectroscopic survey (Smith et al. 2016; commissioned for early 2021) will obtain around 10⁶ optical spectra of LOFAR-detected radio sources, allowing for an accurate redshift determination and source classification. This survey will provide the opportunity to reveal high-redshift (z > 6) radio galaxies, which are necessary for advancing our current understanding of galaxy formation and evolution and can be used as probes of the EoR. Follow-up observations with existing 10 m class optical telescopes and upcoming facilities, such as the James Webb Space Telescope and the European Extremely Large Telescope, will enable detailed studies of these high-redshift galaxies for the first time.

Fig. 10 Stacked SEDs (blue) with posterior median SED models (yellow). Left: stacked likely z = 6.6 LAE candidates in the NB921 catalogue, consisting of the non-matched sample and probable LAEs in the matched sample. Right: stacked likely contaminants in the NB921 matched optical sample, fitted with fixed z = 1.47. The stack of likely contaminants shows brighter NIR and mid-IR fluxes, highlighting the importance of using NIR data in LAE selection.

Acknowledgements

K.J.D. and H.J.A.R. acknowledge support from the ERC Advanced Investigator programme NewClusters 321271. J.S. and P.N.B. are grateful for support from the UK STFC via grant ST/R000972/1. I.P. and M.B. acknowledges support from INAF under PRIN SKA/CTA FORECaST and MAIN STREAM SAUROS projects. K.M. has been supported by the National Science Centre (UMO-2018/30/E/ST9/00082). Ma.B. acknowledges support from the Ministero degli Affari Esteri della Cooperazione Internazionale – Direzione Generale per la Promozione del Sistema Paese Progetto di Grande Rilevanza ZA18GR02. R.K. acknowledges support from the Science and Technology Facilities Council (STFC) through an STFC studentship via grant ST/R504737/1. This paper is based (in part) on data obtained with the International LOFAR Telescope (ILT) under project codes LC0 015, LC2 024, LC2 038, LC3 008, LC4 008, LC4 034 and LT10 01. LOFAR (van Haarlem et al. 2013) is the Low Frequency Array designed and constructed by ASTRON. It has observing, data processing, and data storage facilities in several countries, which are owned by various parties (each with their own funding sources), and which are collectively operated by the ILT foundation under a joint scientific policy. The ILT resources have benefitted from the following recent major funding sources: CNRS-INSU, Observatoire de Paris and Université d’Orléans, France; BMBF, MIWF-NRW, MPG, Germany; Science Foundation Ireland (SFI), Department of Business, Enterprise and Innovation (DBEI), Ireland; NWO, The Netherlands; The Science and Technology Facilities Council, UK; Ministry of Science and Higher Education, Poland. The Hyper Suprime-Cam (HSC) collaboration includes the astronomical communities of Japan and Taiwan, and Princeton University. The HSC instrumentation and software were developed by the National Astronomical Observatory of Japan (NAOJ), the Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research Organisation (KEK), the Academia Sinica Institute forAstronomy and Astrophysics in Taiwan (ASIAA), and Princeton University. Funding was contributed by the FIRST program from Japanese Cabinet Office, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the Japan Society for the Promotion of Science (JSPS), Japan Science and Technology Agency (JST), the Toray Science Foundation, NAOJ, Kavli IPMU, KEK, ASIAA, and Princeton University.

Appendix A BAGPIPES model parameters

Table A.1 shows the parameters we used to fit to our data in BAGPIPES (Carnall et al. 2018). The SFH is modelled using a double-power-law given by (A.1)

where α is the falling power slope, β is the rising power slope, and τ is the time of peak star formation. A prior distribution uniform in log₁₀ is used for α and β. Further information on the model components can be obtained online⁴.

Appendix B SED fits of LOFAR-detected sources

The SED fits and reduced χ² as a function of redshift are given in Fig. B.1 for four sources in the LOFAR-detected sample.

Fig. B.1 SED fitting results obtained using  BAGPIPES for ILTJ160227.18+544759.8, ILTJ160930.42+544435.9, ILTJ16209.13+544818.7, and ILTJ161447.04+550401.8 (in order from top to bottom). Left: reduced χ2 as a function of redshift. The redshifts corresponding to the [OIII], [OII], and Lyα emission lines are indicated, and both resulting reduced χ2 values using normal error (pink) and 10% extra flux error (grey) are shown. Right: posterior median model obtained from BAGPIPES (yellow) with z = 1.47 and over-plotted photometric fluxes as a function of wavelength (blue).

References

All Tables

All Figures

	Fig. 1 Filter transmission curves of HSC BB (grizy) and NB (NB816 & NB921) filters. The central wavelength (λcen) and full width at half maximum (in nm) are indicated for the i, z, NB816, and NB921 filters. The filter transmission curves and λcen values have been obtained from the SVO Filter Profile Service (http://svo2.cab.inta-csic.es/theory/fps/index.php?mode=browse&gname=Subaru&gname2=HSC).
In the text

	Fig. 2 Multi-wavelength cutouts of the LOFAR-detected sample including, from left to right: HSC r band, HSC z band, HSC NB921 band, UKIDSS K band, IRAC 4.5 μm, and LOFAR 150 MHz. The sources show negligible or no emission below the Lyman break in the r band; however, clear emission is seen in both the UKIDSS K band and IRAC.
In the text

Fig. 3 SED fitting results for ILTJ160658.74+550607.0, obtained using  BAGPIPES. Left: reduced χ2 for SED fits in the range 0 < z < 8, with locations of the Lyα, [OII], and [OIII]+Hβ emission lines indicated by dashed black lines. The dashed grey lines indicate the redshift range probed by the NB filter for the strong emission lines [OIII], [OII], as well as Lyα. Right: observed flux of ILTJ160658.74+550607.0 (blue) in mJy with the original errors and posterior median model (yellow) for fits with redshifts fixed at z = 0.84, 1.47, and 6.57. The source is best fitted by a template with z = 1.47 when considering the possible emission lines. No other emission line is strong enough for the source to be at z ~ 4.

In the text

Fig. 4 K and NB magnitudes (NB816 or NB921) from the optically selected sample and the non-detected SILVERRUSH sample. Limiting magnitudes (1σ) are given for sources below 1σ detection (filled and unfilled triangles for sources in the optically selected and non-detected SILVERRUSH samples, respectively). The histograms only include sources with a >1σ detection. The median K and NB magnitudes are indicated by the solid and dashed lines for the optically selected and non-detected samples, respectively. In the optically selected sample, 23% of sources are not detected in K above 1σ, while for the non-detected SILVERRUSH sample, this fraction is 82%. The optically selected sample is brighter and redder than the original SILVERRUSH sample.

In the text

Fig. 5 Left: stacked photometric redshift probability distribution (sum of individual sources) of the optically cross-matched sample of the NB921 and NB816 LAE catalogues from D20. Right: histograms of the best-fitting solutions for SED fits performed with fixed-redshift solutions corresponding to the three candidate emission lines. The known possible redshifts from the NB emission line detections and more detailed SED fitting allows us to exclude other high-z solutions that are sometimes favoured by the photometric redshift code.

In the text

	Fig. 6 Comparison of χ2 values obtained from SED fitting models of Lyα, [OII], and [OIII] for the optical cross-matched LAE sample of NB921 and NB816. The dashed lines indicate Δχ2 values of 2 and − 2. For values higher than 2 and lower than −2, there is a preferred model, but in between the dashed lines there is no preference for either model. The bin width is set to 2.
In the text

Fig. 7 Derived line luminosity, EW, SFR, and stellar mass of the optical cross-matched LAE sample with the NB921 and NB816 LAE catalogues at different fixed redshifts of z = 0.63/0.84, 1.19/1.47,and 5.72/6.57. The dashed blue and orange lines indicate the median for NB921 and NB816, respectively; the values given in the legend include errors derived from the central 68% range of the distribution. The lowest bin in the SFR histograms show the sources with SFR < 10−1 M⊙ yr−1.

In the text

Fig. 8 Stacked LOFAR cutouts of LAE candidate samples. Left: NB816 full sample (226 sources; upper left), optically selected sample (53 sources; middle left), and likely Lyα emitters (177 sources; bottom left). Right: NB921 full sample (322 sources; upper right), optically selected sample (49 sources; middle right), and likely Lyα emitters (279 sources; bottom right). The LOFAR-detected sources and sources not overlapping with LOFAR coverage have been removed.

In the text

	Fig. 9 UKIDSS J and K band magnitude of the optically selected sample and LOFAR-detected sample, with colour-coded photo-zs from D20. The diagonal (dashed) line indicates the J−K ≤ 0 criterion. The sources above this line satisfy this criterion and also have notably higher photo-z values.
In the text

Fig. 10 Stacked SEDs (blue) with posterior median SED models (yellow). Left: stacked likely z = 6.6 LAE candidates in the NB921 catalogue, consisting of the non-matched sample and probable LAEs in the matched sample. Right: stacked likely contaminants in the NB921 matched optical sample, fitted with fixed z = 1.47. The stack of likely contaminants shows brighter NIR and mid-IR fluxes, highlighting the importance of using NIR data in LAE selection.

In the text

Fig. B.1 SED fitting results obtained using  BAGPIPES for ILTJ160227.18+544759.8, ILTJ160930.42+544435.9, ILTJ16209.13+544818.7, and ILTJ161447.04+550401.8 (in order from top to bottom). Left: reduced χ2 as a function of redshift. The redshifts corresponding to the [OIII], [OII], and Lyα emission lines are indicated, and both resulting reduced χ2 values using normal error (pink) and 10% extra flux error (grey) are shown. Right: posterior median model obtained from BAGPIPES (yellow) with z = 1.47 and over-plotted photometric fluxes as a function of wavelength (blue).

In the text