A simultaneous search for High-$z$ LAEs and LBGs in the SHARDS survey

SHARDS is an ESO/GTC deep optical spectro-photometric survey of the GOODS-N field acquired with 200 hours of observing time with the OSIRIS instrument on the 10.4 m Gran Telescopio Canarias (GTC). The surveyed area is $\sim 130$ arcmin², split into two pointings. This survey provides very deep photometry ($m\lesssim 26.5-27.0$ AB mag, at the $3\sigma$ level) in 25 medium-band filters, from 500 nm to 941 nm. The first 22 filters have $\sim 17$ nm FWHM, while the last three have a FWHM of 35, 25 and 33 nm, respectively.

The main purpose of SHARDS was to provide SEDs with good enough spectral resolution to study red and passive galaxies up to redshift $z\sim 2$. Nonetheless, they can be used to detect line emitting sources that otherwise would not be detected using broader filters (Cava et al., 2013; Rodríguez Espinosa et al., 2014; Hernán-Caballero et al., 2017; Lumbreras-Calle et al., 2018), as well as drop-out sources that can be well characterized by the many consecutive filters. Indeed, the relative narrowness of the SHARDS filters and large number of them result in good spectrally resolved SEDs ($R\sim 50$). This is what makes SHARDS an excellent survey to search and study LAEs and LBGs at redshifts between 3.35 and 6.8, whose Ly$\alpha$ emission line or Lyman break fall within the wavelength range of SHARDS. The SHARDS data reduction and calibration is explained in great detail in Pérez-González et al. (2013).

We have also made use of additional public data from HST/ACS (Giavalisco et al., 2004b; Riess et al., 2007), HST/WFC3 (Grogin et al., 2011; Koekemoer et al., 2011), Subaru/Suprime-Cam (as a part of the Hawaii Hubble Deep Field North project, Capak et al., 2004), Subaru/MOIRCS (Kajisawa et al., 2011), CFHT/WIRCam (Lin et al., 2012), Spitzer/IRAC (Fazio et al., 2004; Pérez-González et al., 2005, 2008; Ashby et al., 2015) and GALEX (Bianchi et al., 2014) as compiled in the Rainbow Cosmological Surveys Database¹¹1Operated by the Universidad Complutense de Madrid (UCM), partnered with the University of California Observatories at Santa Cruz (UCO/Lick, UCSC).
http://rainbowx.fis.ucm.es/Rainbow_navigator_public/ (Barro et al., 2011, 2011b, 2018, in prep.). We have therefore added valuable data in wavelengths beyond the SHARDS range. This is important to ensure that our candidates do not present significant emission bluewards of the Lyman break, and in general to discard interlopers.

We start by limiting our search field to the common area of all the SHARDS images, for each pointing, as each SHARDS image samples a slightly displaced area. Thus the detection code only takes into account sources that are inside a resulting effective area of 128.4 arcmin² (otherwise we get many spurious detections near the borders of the frames).

We also need to take into consideration the effective central wavelength (CWL) within each filter, which depends on the particular position of each source in the field of view. This is a pure geometrical effect, as a result of the angle of incidence of the GTC/OSIRIS light beam on the filter. Objects in the same SHARDS image are detected with similar effective wavelengths but not exactly the same, depending on their precise position in the image. Fortunately, this is an effect already calibrated (see Pérez-González et al., 2013), so we can get the exact CWL for each object in each filter directly from the SHARDS catalogue. The appropriate FWHM of the band sampling that CWL was also considered, as this varies slightly for each filter.

Here we explain in detail the process used to search for the candidate sources in the SHARDS catalogue. SHARDS contains 44,752 sources in the GOODS-N field, and for each object there is information on its flux in the SHARDS filters, plus the ancillary data available in Rainbow (see Section 2.1).

A first selection was made taking advantage of the photometric redshifts available in the Rainbow Database for 35,445 of the SHARDS objects in the CANDELS/GOODS-N field. These redshifts were obtained by Barro et al. (2018, in prep.), hereinafter Ba18, through SED fitting with EAZY (Brammer et al., 2008), using SHARDS combined with the WFC3 102 and 141 grisms and additional data from HST, Subaru, Spitzer and others. This results in a very well detailed SED characterization and redshift calculation, reaching an accuracy better than $\Delta z\sim 0.01$ (see also Ferreras et al., 2014; Domínguez Sánchez et al., 2016). From this sample, 1078 candidates were selected using their redshift estimation.

For the 9,307 objects in the SHARDS catalogue, without photometric redshift in Ba18, a selection criterion based on colour excesses between filters was designed. Since we are interested in detecting both Ly$\alpha$ in emission and/or the Lyman break feature, we first constructed colour-magnitude diagrams using consecutive filters. Thus we were able to detect all relevant flux changes from one filter to the next (regardless of whether they were emission lines or drop-outs). In order to only select objects with reliable flux changes, we studied the dependence of magnitude errors with magnitude in the SHARDS catalogue for each one of the filters. We then estimated the mean error for a given magnitude and determined whether a flux change between filters was large enough. If, on the contrary, this error was within the expected errors for that magnitude, we discarded the source. The procedure used first order logarithmic fits of the magnitude error for each filter. We used $1\sigma$ error limits to identify relevant emission changes (see Fig. 1).

We then built SHARDS SEDs for the objects that satisfied the colour-magnitude criterion, and started a more selective filtering. For this purpose we implemented a code that analyses the SEDs of the objects and selects candidates according to the following criteria:

The last step was to visually inspect all the candidate objects, not only looking at their SEDs but also directly viewing the images. In this process we used images mainly from GTC/SHARDS, HST/ACS, and HST/WFC3, though we also used other images available in the Rainbow Database from Subaru/Suprime-Cam and Spitzer/IRAC. The SEDs were further completed with useful data at other wavelengths available in Rainbow (see Section 2.1). Finally, we visually discarded those objects which were deemed either spurious or presented strong emission blueward of their supposed Lyman break wavelength. This check is necessary since lower $z$ objects with a high slope SED could match our SHARDS preselection criteria when their emission in one of the redder filters is substantially higher than the average emission in the filters blueward of that one.

The whole process was iterated many times, identifying the wavelength where the Lyman break appears for each object, paying special attention to galaxies with photometric redshift estimations from previous authors, making sure we did not skip possible candidates while not including unwanted objects. Using this last method we added 492 objects. These together with the 1078 selected objects from Ba18 mentioned before constitute our final selection sample, consisting of 1570 well characterized candidates, pending a final interlopers screening.

Once the selection was finished, we carried out some additional theoretical models fits to our SEDs to make a further test on those objects that had met our selection criteria but could still be lower $z$ interlopers. In the search for high-$z$ LBGs the main interlopers are lower $z$ galaxies where the D4000 break is wrongly identified as the Lyman break. Additionally, M, L, and T dwarf stars with a very steep slope can be mistakenly taken as drop outs. In these, it is possible to measure an abrupt flux change, specially with broad band filters. However, these last cases are in general easily recognizable when using many filters. So we did not expect an important number of them among our interlopers. The fitting process was made with the code Le Phare (Arnouts et al., 1999; Ilbert et al., 2006). We used different galaxy models (e.g., Coleman et al., 1980; Kinney et al., 1996; Bruzual & Charlot, 2003), assuming a Salpeter (1955) Initial Mass Function (IMF) and a Calzetti extinction law (Calzetti et al., 1994). For the stars rejection, we used SEDs from Bohlin et al. (1995), Kinney et al. (1996) and Chabrier et al. (2000), with a special attention to cold M, L and T stars from Pickles (1998), Burgasser et al. (2004) and others.

As in every other survey, we did not get good SED fitting for all the 1570 candidate galaxies in our sample, but only for 644 (with $\chi^{2}\lesssim 40$). Therefore, we expected some unwanted objects within the worst fitted sources of the sample. Only 12 out of the well fitted part of the sample were found to be lower $z$ interlopers. None of them were found to be steep red-slope stars, but lower $z$ galaxies. To estimate the final interloper fraction in our high-$z$ galaxy selection we extrapolated the interloper fraction found for the well fitted sample to the whole sample, obtaining $\sim 2$% fraction of interlopers in our final selection. Thus, after eliminating the identified interlopers, our final sample consists of 1558 high-$z$ galaxies, classified as 528 LAEs and 1030 LBGs with no emission line. 124 of the emitters are pure LAEs with no continuum detected in SHARDS (although it can be measured in the HST/ACS broad band images), while the rest (404) are LBGs with Ly$\alpha$ line emission. Figures 2 and 3 show examples of a good candidate SED and its view on each SHARDS filter, respectively.

In order to properly correct the Luminosity Functions (LF) and to compare our sample with previous ones, we estimated our completeness. To achieve this, we took the rest-frame UV magnitude at 1500 Å as reference magnitude and made a plot of the logarithm of the cumulative number of LBGs with $m_{1500}$ below some apparent magnitude. As we do not reach faint enough magnitudes to sample the potential part of the LF, our selection will be mostly dominated by the exponential term of the Schechter function. It is thus reasonable to assume that we can calculate our completeness looking at the $m_{1500}$ at which the slope of our logarithmic cumulative distribution starts decreasing. In this way we obtain the value beyond which we start missing candidate objects. This will be our completeness apparent magnitude (see Fig. 4). We obtain a 90% completeness at magnitude $\sim 25.87$ AB mag, although our faintest sources reach up to $m_{1500}\sim 28$. This completeness value corresponds to a different absolute magnitude for each different redshift.

Section 2.3 showed how the SED of each object was used to identify the filter wavelength at which that object presents either an emission line or a continuum drop. The main assumptions, for obtaining the redshift of the objects, is that the line is Ly$\alpha$ and the emission drop is the Lyman break. These premises are supported by the selection criteria described in Sections 2.3 and 2.4. With these assumptions, we get their break wavelength, taking into account the precise position in the image to correct the CWL shift effect described in Section 2.2. Our preliminary calculations pointed out that this effect is relevant when the object is located far from the central region. We also take care of using the appropriate filter width corresponding to the pre-break filter. It is important to clarify that the Lyman alpha break wavelength used is not 912 Å but just the blue side of Ly$\alpha$, since the Ly$\alpha$ forest due to intergalactic absorption erases almost any continuum emission blueward of 1215.7 Å at high redshifts (Rauch, 1998; Fan et al., 2006b).

To prove the goodness of our photometric redshifts we check every case where a spectroscopic redshift was available either in the Rainbow Database or in the NASA/IPAC Extragalactic Database (NED)²²2https://ned.ipac.caltech.edu/ (from Steidel et al., 2003; Kakazu et al., 2007; Barger et al., 2008; Kajino et al., 2009; Adams et al., 2011; Conselice et al., 2011; Shim et al., 2011; Stark et al., 2011, 2013; Schenker et al., 2013, among others) and compare them with our photometric redshifts (see Fig. 5). The excellent correlation obtained indicates that we can trust the calculated redshifts for our candidates with an average error of around $\Delta z=0.07$ (0.14 for objects above $z\sim 6.2$, where the filters measuring the Lyman break are wider).

The redshift distribution of our LAE and LBG candidates is shown in Fig. 6. The number of LBGs decreases with $z$, which may be partly due to the lower number of SHARDS filters providing relevant information in the red part of the spectrum, apart from the Malmquist bias, as farther away objects are fainter. We also notice that the relative number of LAEs is larger beyond $z\sim 5.0$-5.5. The fraction of pure LAEs and LBGs and their behaviour are discussed in depth in Arrabal Haro et al. (2018, in prep.). The complete list of redshifts is in the on-line version.

Two previously reported GOODS-N overdensities at $z\sim 4.05$ Daddi et al. (2009) and $z\sim 5.2$ Walter et al. (2012) have also been recovered in this work. The number of objects within our sample which may belong to those overdensities according to their redshifts are 87 and 55, respectively. For the $z\sim 5.2$ proto-cluster, those authors find 13 objects at $z=5.198\pm 0.015$ in the GOODS-N field, 11 of which are also detected by us. Of the two not detected, one is too faint and the other one, HDF850.1, was detected via cold molecular gas lines (Walter et al., 2012), not showing any emission in the optical/NIR. We have further detected 44 additional objects whose redshifts are compatible with $z=5.198\pm 0.015$, within the errors. They are all detected with the Ly$\alpha$ line/break in the SHARDS filter f755w17 and have an estimated redshift error of 0.07. Although we have to wait for spectroscopy of these sources to further constrain their redshifts, a total of 55 possible cluster members for a single proto-cluster is something unseen beyond z=5 (cf., e.g., Capak et al., 2011; Toshikawa et al., 2012, 2014; Higuchi et al., 2018), making this a very interesting overdensity.

Traditionally, two different methods have been used to obtain SFRs for the LAEs and LBGs respectively. The emitters SFR can be obtained from the Ly$\alpha$ line emission following Kennicutt (1998), who assumed case B recombination (Brocklehurst, 1971) and a Salpeter IMF (Salpeter, 1955), through the relation:

where $L(\mathrm{H}\alpha)=L(\mathrm{Ly}\alpha)/8.7$.

The other method is based on the UV continuum luminosity density at 1500 Å ($L_{1500}$) with the assumptions of solar metallicity and a Salpeter IMF, as in Madau et al. (1998):

To correct for galactic dust extinction, $A_{\lambda}$ values from the Schlafly & Finkbeiner (2011) dust maps are used. The internal dust extinction is calculated following Calzetti et al. (1994), Calzetti et al. (2000) and Castellano et al. (2012), assuming $F_{\mathrm{\lambda}}\propto\lambda^{\beta}$, where $\beta$ is the UV slope. This $\beta$ slope is estimated via linear fit to the magnitudes measured in the filters sampling the rest-frame UV wavelength from 1300 Å through to 2600 Å:

A problem with the Kennicutt method is that the Ly$\alpha$ line is very sensitive to both HI resonant scattering and dust absorption. Therefore the actual Ly$\alpha$ emission of the galaxy can be greatly affected. Indeed, the photon escape fraction is very uncertain (see, e.g., Hayes et al., 2010). On the other hand, the continuum is much less affected by absorption and scattering due to neutral gas so we assume that SFRs calculated using the rest-frame UV continuum luminosity at 1500 Å are more reliable. In those cases where we have LBGs with emission line we can compare their SFRs calculated via one method and the other. In every case, the SFRs calculated using eq. (3) are much larger, reaching three orders of magnitude difference for 11 candidates with bright $M_{1500}$ but a small Ly$\alpha$ line (2.7% of the sample presenting both line emission and well measured UV continuum). This large difference can be due to: 1) the uncertainties in measuring the Ly$\alpha$ line extinction, and 2) the fact that the continuum measurement integrates not only the current starburst but also the previous history of star formation of that source. Moreover, for those LBGs with line emission, the ratio between the two SFRs seems to be stochastic, since, as we explain in detail in Arrabal Haro et al. (2018, in prep.), the intensity of the line emission relative to the continuum heavily depends on the strength of the current starburst, which quickly decays within a few Myr. That said, the ratio between SFRs, measured with the Kennicutt and the Madau methods, will only depend on the current stage of the young starburst and the number of former recent starbursts undergone by a given source.

The SFRs, corrected for extinction following the above prescription, are given in Table 1, where we show the SFR based on $L(Ly\alpha)$ and/or $L_{1500}$ depending on whether the galaxy presents emission line or continuum emission. As expected for these early galaxies, the SFRs are fairly high (e.g., Bouwens et al., 2007, 2015; Mashian et al., 2016; Casey et al., 2017), with those calculated with $L_{1500}$ being $\sim 58\ M_{\odot}\ \mathrm{yr^{-1}}$ on average. The complete LBGs SFR distribution can be seen in Fig. 7, where there are 7 objects forming stars at a rate above 10³ $M_{\odot}\ \mathrm{yr^{-1}}$. These galaxies are very dusty star-forming galaxies (see Casey et al., 2017) whose SEDs show a very red UV slope. An example is given in Fig. 8.

The EW of the Ly$\alpha$ line tells us the relative strength of the young star formation burst. It relates the line emission due to the current starburst to the continuum produced by the more evolved population of the galaxy. For every LBG with line emission we obtain the EW measuring the continuum either in the adjacent SHARDS filter or, if this filter has a large uncertainty, calculating a weighted mean of the fluxes in nearby filters. If it is not possible to get a continuum from the SHARDS filters, the continuum is obtained from HST broad band data. To establish a minimum reliably detected line flux we impose that the filter measuring the line emission should be at least $1\sigma$ above the continuum flux level. Applying this, we reach a minimum measured rest-frame $\mathrm{EW}=5.14\pm 4.70$ Å.

Measuring a continuum for pure LAEs is difficult as these objects are characterized by having very low or no detectable continuum emission. In some of the cases we can use broad band data but there are a few ones (28) where we cannot reliably measure their continua and therefore we cannot calculate the EW. As expected, pure LAEs have larger EWs than the rest of LAEs, although their uncertainties are also larger, since their continuum is much fainter or even absent. In our sample, the pure LAEs EWs are in the range 35-353 Å with a median value of 109.13 Å, while the LBGs/LAEs EWs are in the range 5.14-217 Å with a median of 27.83 Å. Nonetheless, the number of sources decreases as we move to higher EWs (see right panel in Fig. 9). Notice that for the LBGs we plot only those 404 that are also LAEs, which means that there are 1030 additional LBGs with no emission line detected, i.e., with EW$<5.14$ Å.

In this work we have detected LAEs up to quite low EWs. However, in order to compare with previous studies, we have analysed the fraction of them with EW larger than a standard commonly used value of Ly$\alpha$ $\mathrm{EW}>25$ Å. The comparison is done with various authors at different redshifts (Pentericci et al., 2011; Stark et al., 2011; Curtis-Lake et al., 2012; Ono et al., 2012; Treu et al., 2013; Schenker et al., 2014; Tilvi et al., 2014; Cassata et al., 2015; De Barros et al., 2017; Caruana et al., 2018). Results are shown in Fig. 10, where we divided the sample into two brightness groups. Our fraction of LAEs with $\mathrm{EW}>25$ Å ($X_{\mathrm{Ly}\alpha}$) matches well those previously obtained, although we register a slightly smaller fraction of the faintest sources at $z\sim 4$ and $z\sim 5$. The decrease of $X_{\mathrm{Ly}\alpha}$ beyond $z\sim 6$ found by other authors (e.g., Caruana et al., 2014; Schenker et al., 2014) suggests an increase of the H i abundance as we move to a non-fully reionized universe (see, e.g., Fan et al., 2006), which would strongly affect the scattering of the Ly$\alpha$ line.

To find objects with close neighbours, we look for galaxies at the same redshift (within redshift errors) that are inside a 30 kpc radius circular region, assuming that objects farther away would not be part of the group, as considered, e.g., in Mundy et al. (2017). The distance between objects in the group is calculated from their angular separation using the angular diameter distance at their redshift as it appears in Hogg (1999).

202 of the candidates are in 92 close groups of two or more galaxies, which represent $\sim 13$% of the total sample. In particular, we register, 79 pairs, 10 trios, 2 quartets and a sextet. The remaining sources are isolated. Previous works looking for pairs or very close groups at these redshifts find similar fractions (Conselice & Arnold, 2009; Mundy et al., 2017).

The distance between the objects in these groups vary from 59 kpc to 3.26 kpc. Three of the trios and one of the quartets are confined in a quite small area of around 20 kpc diameter. Many of the other trios are composed of a very close pair separated a larger distance from the third object of the group. Something similar happens with the sextet, which looks like four very close objects (one of which seems to be three different non-resolved sources) in a $\sim 14$ kpc region, with the fifth and sixth ones somewhat more separated, showing a large and conspicuous tail between them. Table 2 lists the groups we found with their assigned names, indicating the distance between pairs or, in the case of groups with more members, the largest distance between the central object and the rest. The images of the quartet (G6) in the SHARDS bands sampling the Lyman break are shown in Fig. 11. We will discuss more about these close groups in Section 4.3.

Keeping in mind our completeness magnitude, we compare our high-$z$ galaxy sample results with previous works. As reference, we took the big survey carried out by Bouwens et al. (2015), hereinafter B15. In that work the LBGs search is made using HST broad band colour criteria. They find a much higher number of objects than we do. Indeed, in the GOODS-N field, B15 find 3917 LBGs, out of which 3455 fall within our redshift range according to their classification. Instead we only selected 1558 objects. To understand this discrepancy we analyse the B15 GOODS-N sample.

First of all we remove the GOODS-N B15 objects that fall out of our field of view. Their images survey 133.7 arcmin² in the GOODS-N field, while our effective search field is only $\sim 128.4$ arcmin². Moreover, the centring and exact shape of the two fields differ slightly. We miss 405 galaxies (11.4%) just because of this effect, leaving still 3050 objects to explain. However, only 2757 of them have photometric detection in at least one of the SHARDS filters. This means we are missing 293 sources, or 9.6% of the B15 sample falling within our field, of which there are no detections in SHARDS.

Another effect to consider is the loss of sources due to neighbour contamination. This is so because of the higher spatial resolution of the HST in comparison to our ground-based GTC/OSIRIS resolution. To quantify this effect, we mask the stars and other big objects in the SHARDS images and analyse how many of the B15 candidates would be affected by them. For the smaller objects, we assume that any B15 galaxy with a brighter (or at most 0.5 mag fainter) object at a distance of $\theta=0.9$ arcsec or less (typical SHARDS seeing) is at least partially contaminated in our images, and therefore we could be missing it. 192 objects (7%) are affected by this issue. The $m_{1500}$ distribution of the objects missed because of non-detection or neighbour contamination is shown in 12. At this point there are still 2565 galaxies of the B15 that should be detected.

We then limit the remaining B15 sample to those objects that should be in our selection. That is sources brighter than our completeness magnitude, within our field, not affected by proximity contamination and whose redshifts fall in the range we can measure with our 25 SHARDS filters. As we are going to exclude galaxies falling outside our completeness $m_{1500}$, we first calculate the wavelength corresponding to 1500 Å rest-frame for each B15 object, which depends on $z$, to identify the HST filter that better samples that wavelength in each case. Then we look for every object in its corresponding filter in the HST/ACS catalogue to check whether or not we should detect it with the SHARDS data. Up to $m_{1500}=25.87$ AB, there are 575 sources meeting all our limiting conditions, however, out of these, 115 (20.0%) are not actually high-$z$ galaxies neither in our study nor according to the photometric redshifts from Ba18.

That said, we are conscious that the magnitude limit and spatial resolution in our search cause a loss of sources. Nevertheless, even considering these effects, we still find a non negligible fraction of objects in the B15 sample that meet our limiting magnitude condition but are not selected applying our criteria. These missing sources are visually revised in the SHARDS images to confirm they are bright enough and isolated enough not to be contaminated by any neighbour. We find that these sources do not look like LBGs when seen with our higher spectral resolution. Instead, these objects seem to be either lower $z$ Balmer break galaxies, cool stars or just red galaxies. Indeed, when these are seen through a few broad band filters only, they could be easily misclassified as LBGs, as they may present a reasonable flux drop from one broad filter to the next. We therefore conclude that, thanks to the much better wavelength discrimination of the SHARDS survey, we identify these objects as interlopers.

Note that we have not calculated photometric redshifts for objects not in our sample. However, as they have not been selected with our criteria, they are considered lower $z$ sources. In order to study the redshift distribution of these not selected sources, we make use of the redshifts calculated in Ba18. In addition, we repeat the analysis previously described in this section using different values of magnitude limit and flux contamination distance (seeing). The number of B15 galaxies that meet our limiting conditions but whose SHARDS photometric $z$ is much lower than that given in B15 is shown in Table 3, for different values of $m_{1500}$. In this Table, we use both 0.9 arcsec, corresponding to the typical SHARDS seeing, and 1.2 arcsec, a bit above that seeing, as the minimum distance between objects at which they can be considered isolated. There is little change in the interloper fraction for limiting magnitudes up to $\sim 26.13$ AB, mostly constant at $\sim 20$%. Beyond that, this fraction slightly increases (up to 22.3% at $\sim 27.5$ AB), which makes sense since fainter objects are more difficult to characterize and identify as either good galaxy candidates or interlopers.

Therefore we conclude that $\sim 20$% of the B15 sample, up to $m_{1500}\sim 25.87$ AB, is actually misclassified as high-$z$ LBGs (quite above the $\sim 2$-6% found by Vulcani et al., 2017). It is interesting to see the redshift distribution of these misidentified objects. Figure 13 is a histogram representing them as a function of redshift for different $m_{\mathrm{lim}}$ and isolation distances. We find a clear concentration of sources at very precise redshifts. Moreover, the shape of the interlopers redshift distribution is practically the same independently of the depth and seeing considered. This implies that these interlopers are not due to seeing or brightness issues, but to the spectral resolution quality. In fact, the distribution peaks around $z\sim 0.5$, and a big fraction ($\sim$70%) of the misclassified sources are distributed between $z=0.3$ and $z=0.9$, corresponding to a range between $z\sim 3$ and $\sim 5$ if the Balmer break is mistakenly taken as the Lyman break. We would actually expect a similar fraction of interlopers in other broad band studies with similar colour selection criteria. Fig. 14 shows an example of one of these misidentified objects selected in B15.

Finally, for some galaxies in the B15 sample that also meet our selection criteria, the redshift assigned in B15 does not match ours within a $\Delta z=0.3$ error (approximately twice our highest photo-$z$ error). An example of this is shown in Fig. 15. To better characterize this discrepancy, we study the subsample of objects with photo-$z$ calculated both in B15, Ba18 and this work, consisting on 933 galaxies. The B15 redshift of 138 of these sources ($\sim 15$%) are found to differ from both Ba18 and our calculated redshifts. We consider our redshift determination much more accurate since we have a much better spectral resolution in our data, and our redshifts are in an excellent agreement with spectroscopic values as shown in Fig. 5. Indeed, the photo-$z$ obtained with broad band filters show a considerable uncertainty.

Interestingly, there are other 312 objects in our sample that are not selected by B15. Of them, 121 are placed in two little gaps of the B15 GOODS-N field, as shown in Fig. 16. Of the remaining 191, 17 have spectroscopic $z$ available in Rainbow and/or NED that matches the photo-$z$ calculated in this work. A closer look to their SEDs does not give us any conclusive hint about why they are not selected using broad band filters. An example of these objects is shown in Fig. 17. Some of these galaxies are pure LAEs with strong line emission and a very faint continuum that could have been missed with broad band filters. This explanation is however valid for just 16 objects out of the 191. We do not know why the rest are not selected with the broad band criteria.

Summarizing, from the 3455 objects detected by B15 in GOODS-N within our $z$ range, after considering limiting brightness difference, small field of view differences and possible neighbouring light contamination, we would expect to detect 575 of them up to our 90% completeness magnitude. Of them, 115 are not selected and are classified as low $z$ objects, implying that $\sim 20$% of the B15 broad band selected sample up to $m_{1500}=25.87$ AB mag seem to be interlopers. On the other hand, we have further found 191 galaxies that we consider good candidates and were not selected in B15.

We proceed now to build the LBGs LFs, which are shown in Fig. 18. Our 90% completeness absolute magnitude limit for each redshift is marked with a dashed vertical line. A correction for those objects that we miss, either because of non-detection in the SHARDS images or because of neighbouring light contamination of their SEDs, is considered for each $z$ range, based on their $m_{1500}$ distribution (see Fig. 12). The magnitude values are measured using HST photometry. To correct our data in the faint region of the LF, we use the $V/V_{\mathrm{max}}$ correction to adequately consider the effective volume sampled by each detected source. We divide our sample in three main redshift bins to facilitate comparing with the literature, namely $z\sim 4$ ($3.5\leq z<4.5$), $z\sim 5$ ($4.5\leq z<5.5$) and $z\sim 6$ ($5.5\leq z<6.5$). A Schechter function is used to fit the data:

Since we are missing a fair number of the faintest objects, we cannot get information about the power-law part of the Schechter function. Therefore, the $\alpha$ parameter cannot be directly calculated, so we fix its value to a value consistent with the literature. In particular, we set it at each redshift to the mean value obtained from the following authors: Bouwens et al. (2007), van der Burg et al. (2010) and B15 at $z\sim 4$, Bouwens et al. (2007), Iwata et al. (2007), McLure et al. (2009), van der Burg et al. (2010) and B15 at $z\sim 5$ and Bouwens et al. (2007), McLure et al. (2009), Bouwens et al. (2012), Bowler et al. (2015) and B15 at $z\sim 6$. It is important to note that the $\alpha$ values derived from the literature may be overestimated according to the large interloper fraction we find in typical surveys made through broad band observations. The bright region can be fitted with the exponential term of the Schechter function to estimate $\phi^{*}$ and $M^{*}$. The results obtained are listed in Table 4. Furthermore, we notice that the brightest points of our LFs tend to be slightly over their fitted Schechter function. This is in agreement with recent big field studies that propose a power-law behaviour for the bright region of the LFs, instead of the classical Schechter function (Sobral et al., 2017, 2017b). Nonetheless, our field is not large enough to claim conclusive results in this regard.

A cosmic variance calculator, developed by Trenti & Stiavelli (2008) using halo models, is also considered to evaluate the effect of cosmic variance and Poisson counts errors, as done e.g., in B15, Conselice et al. (2016) and Vulcani et al. (2017). We thus calculate that the number counts at $z\sim 4$, $z\sim 5$ and $z\sim 6$ has a fractional error due to cosmic variance of 16%, 20% and 28%, respectively. These uncertainties are relatively large in our study as we are limited to a relatively small field.

At $z\sim 6$ we have large uncertainties, due to the scarce statistics, cosmic variance effects and poissonian errors. For the other two redshift ranges, the results are comparable to previous studies (see Fig. 19). Notice that both at $z=4$ and 5, our LF estimations fit the bright part fairly well, decaying quickly as we move to fainter magnitudes.

We do match previous studies in the bright region. However, we find a discrepancy in the intermediate magnitudes region. We are indeed missing candidates because of seeing contamination in this region, but this effect is already corrected, though not enough to explain the difference we find, since the seeing contamination only affects a $\sim 7$% of the cases, as described in Section 4.1. On the other hand, the reliability of our candidates is very high, suggesting that the reason for the discrepancy is not only due to a lack of faint objects in our sample, but to an overestimation in the detection of LBGs when the selection is made with broad band filters. As we have seen before, the number of interlopers in broad band surveys is large and increases with magnitude. This is a possible reason for the increasing deviation of our LFs as we move to fainter magnitudes. We cannot estimate the precise effect that the interloper fraction produces at magnitudes beyond our completeness limit, but we expect it to increase, implying an uncertain calculation of the Schechter $\alpha$ parameter in surveys done with broad band photometry. We would like to stress the importance of having good spectral resolution for building a high-$z$ galaxy sample free from interlopers. As an estimation, we would expect the “real” LF to be somewhere between ours and that given in previous HST works, to account for the seeing and limiting magnitude difference. If this study is correct, previous broad band LFs should be downscaled by a factor $\sim 20$%, which could be even more important at fainter magnitudes. This would presumably lead to a decrease in the faint end slope of the LFs.

As commented in Section 3.4, 92 groups are found within our sample. Almost half of them (48) are pairs separated by longer distances than the rest, though they could still be interacting sources according to the typical distances given in Duc (2014). In the other 44 groups, the galaxies are separated by distances of around 30 kpc (3.90-4.72 arcsec, depending on redshift) or less. For these very close sources it is difficult to distinguish whether they are independent galaxies or star-forming knots in a single clumpy galaxy (Elmegreen et al., 2013). Fitting their SEDs with simple models results in a mass of the order of 10⁸ $M_{\odot}$, in which case they should be high-$z$ galaxies rather than clumps. The typical size of these objects is $\sim 0.4$ arcsec, which is consistent with typical sizes of high-$z$ galaxies as found by (Bouwens et al., 2004; Oesch et al., 2010; Ono et al., 2013; Shibuya et al., 2015; Holwerda et al., 2015; Liu et al., 2017). Moreover, luminous clumps in clumpy galaxies are usually surrounded by a common isophote, embedding the entire galaxy, below the bright level of the clumps but above the background sky (see Hinojosa-Goñi et al., 2016). This effect is not seen in our objects. Because of these reasons we claim that these objects should actually be individual galaxies in compact groups.

In particular, the group G5 is an already spectroscopically confirmed pair of galaxies at $z\sim 5.1$ in which the northern source presents a nebulosity around it that could be either a tidal tail caused by the interaction with the other galaxy or the remnant of a recent merger with a third object (Rodríguez Espinosa et al., 2014). We find similar tail-like morphologies in most of our groups (e.g., the westernmost object of G6 shown in Fig. 11). These structures support the hypothesis that our groups are gravitationally bounded systems. In addition to the tidal tails we also notice some other little spots near the trios and one more near G6 in the HST/ACS images that could be additional members of their respective groups. However, these extra objects do not appear either in the SHARDS or the ACS catalogues nor in the Rainbow database, and they cannot be resolved in the GTC images, so we do not have conclusive photometric information for them. In fact, some of the objects forming these groups seem to be more than one unresolved source.

Note that the number of groups found, though significant, is not large enough to warrant studies of their dependence on redshift. They seem to follow a similar $z$ distribution as the whole sample. Concerning their SFR, we check whether they present particularly high ones. We find that they have similar SFR distribution as the rest of the sources in the sample.

In addition to these very close groups, we have also detected 87 and 55 objects possibly belonging to two already known overdensities at $z=4.05$ (Daddi et al., 2009) and $z=5.198$ (Walter et al., 2012), respectively, in the GOODS-N field. Given our small typical redshift errors ($\Delta z\sim 0.07$) we do suggest they belong to those overdensities. This increases by a factor of $\sim 4$ the number of sources previously found for the $z=5.198$ overdensity. If confirmed, this proto-cluster at $z\sim 5.2$ will be one of the richest proto-clusters beyond $z=5$.