A Spectroscopic Survey of Lyα Emitters and Lyα Luminosity Function at Redshifts 3.7 and 4.8

We present a spectroscopic survey of Lyα emitters (LAEs) at z ∼ 3.7 and z ∼ 4.8. The LAEs are selected using the narrowband technique based on the combination of deep narrowband and broadband imaging data in two deep fields, and then spectroscopically confirmed with the MMT multifiber spectrograph Hectospec. The sample consists of 71 LAEs at z ∼ 3.7 and 69 LAEs at z ∼ 4.8 over ∼1.5 deg2, making it one of the largest spectroscopically confirmed samples of LAEs at the two redshifts. Their Lyα luminosities are measured using the secure redshifts and deep photometric data, and span a range of ∼1042.5–1043.6 erg s−1, so these LAEs represent the most luminous galaxies at the redshifts in terms of Lyα luminosity. We estimate and correct sample incompletenesses and derive reliable Lyα luminosity function (LFs) at z ∼ 3.7 and 4.8 based on the two spectroscopic samples. We find that our Lyα LFs are roughly consistent (within a factor of 2−3) with previous measurements at similar redshifts that were derived from either photometric samples or spectroscopic samples. By comparing with previous studies in different redshifts, we find that the Lyα LFs decrease mildly from z ∼ 3.1 to z ∼ 5.7, supporting the previous claim of the slow LF evolution between z ∼ 2 and z ∼ 6. At z > 5.7, the LF declines rapidly toward higher redshift, partly due to the effect of cosmic reionization.


INTRODUCTION
Spectroscopically confirmed galaxies at high redshift are important for us to understand galaxy properties and evolution in the distant Universe.In the past two decades, galaxies at z > 2 have been routinely found using the dropout (or Lyman break) technique (Steidel et al. 1996;Stark 2016;Harikane et al. 2023).This technique uses strong Lyman break in the spectra of star-forming galaxies to select Lyman-break galaxy (LBG) candidates.Follow-up spectroscopic confirmation of these objects are often observationally expensive because of their faint continuum emission.In addition to the Lyman break technique, the narrowband technique that uses strong Lyα lines has also played an important role in finding high-redshift galaxies.It was predicted 55 years ago by Partridge & Peebles (1967) that primeval galaxies undergoing their initial burst of intense star formation would appear very bright in the redshifted Lyα line.This is used by the narrowband technique to select Lyα emitting galaxy (Lyα emitter, or LAE) candidates.These candidates are typically much easier to be spectroscopically identified due to their strong Lyα emission.The two techniques are highly complementary.
High-redshift LAEs are usually young, compact, metal-poor, and low-mass (stellar mass ∼ 10 8−9 M ⊙ ) star-forming galaxies with star formation rates around ∼ 1 − 10 M ⊙ /yr (Ouchi et al. 2020).Since the discovery of the first high-redshift LAEs using the narrowband technique (Petitjean et al. 1996;Hu & McMahon 1996;Pascarelle et al. 1996), this technique has successfully found a large number of LAEs at z ≥ 2 (e.g., Ouchi et al. 2008;Ciardullo et al. 2012;Sobral et al. 2018;Hu et al. 2019;Guo et al. 2020;Ning et al. 2020Ning et al. , 2022)).Large samples of LAEs from low to high redshift allow us to study their Lyα luminosity function (LF) and evolution (e.g., Cassata et al. 2011;Drake et al. 2017;Konno et al. 2018;Herenz et al. 2019;Guo et al. 2020;Zhang et al. 2021), and physical properties (e.g., Jiang et al. 2013Jiang et al. , 2016;;Shibuya et al. 2019;Steidel et al. 2011).They also allow us to find protoclusters of galaxies (e.g., Jiang et al. 2018;Hu et al. 2021) and characterize the end of cosmic reionization (e.g., Kashikawa et al. 2006Kashikawa et al. , 2011;;Ning et al. 2022).The recently launched James Webb Space Telescope (JWST) allows us to study their rest-frame optical spectral properties and their correlations with the Lyα line (e.g., Roy et al. 2023).The Lyα LF describes the number density of LAEs as a function of Lyα luminosity, and is thus a basic statistical property of LAEs.Current studies show that the Lyα LF increases rapidly from z ∼ 0.3 to z ∼ 3, appears constant from z ∼ 3 to z ∼ 6, and then declines rapidly from z ∼ 6 towards higher redshift (Ouchi et al. 2020).However, the Lyα LF at z ∼ 4 − 5 has not been well explored, and previous studies for this redshift range were mostly based on photometric samples of LAEs.For example, Ouchi et al. (2003) calculated z ∼ 4.86 Lyα LF using a photometric sample consisting of 87 LAE candiadtes.Ouchi et al. (2008) derived z ∼ 3.7 Lyα LF based on 101 photometrically selected LAE candidates, and they spectroscopically confirmed 26 LAEs.Shioya et al. (2009) obtained z ∼ 4.86 Lyα LF with 79 photometrically selected LAE candidates.There are some studies based on spectroscopic samples, but the number is relatively small.For example, Zheng et al. (2013) combined results from the Large Area Lyman Alpha (LALA) survey and derived z ∼ 4.5 Lyα LF based on a large sample of 207 spectroscopically confirmed LAEs.Although the typical confirmation rate of LAE candidates selected by the narrowband technique is high (∼ 60% to ∼ 80%, e.g., Kashikawa et al. 2011;Zheng et al. 2013), a spectroscopically confirmed LAE sample is still important in studying Lyα LF by excluding contaminants and deriving more robust Lyα flux.
In this paper, we present spectroscopic surveys of LAEs at z ∼ 3.7 in the Subaru XMM-Newton Deep Survey (SXDS) field and LAEs at z ∼ 4.8 in the Subaru Deep Field (SDF).We select LAE candidates from deep broadband and narrowband images taken by the Subaru Suprime-Cam and spectroscopically observe them with the MMT Hectospec spectrograph.From these observations, we obtain large samples of LAEs at z ∼ 3.7 and z ∼ 4.8, and we further derive Lyα LFs at the two redshifts.The layout of this paper is as follows.
In Section 2, we introduce target selection and spectroscopic observations.In Section 3, we describe LAE samples and properties.Lyα LFs are derived in Section 4. Section 5 and Section 6 are discussion and summary.Throughout this paper, all magnitudes are in the AB system.We adopt a Λ-dominated flat cosmology with H 0 = 70 km s −1 Mpc −1 , Ω m = 0.3, and Ω Λ = 0.7.

TARGET SELECTION AND SPECTROSCOPIC OBSERVATIONS
In this section, we will describe the imaging data in the SXDS and SDF fields, the selection of LAE candidates, the spectroscopic observations of these candidates, and our data reduction.

Subaru Suprime-Cam images in SXDS and SDF
The Suprime-Cam is a wide-field prime-focus imager with a field-of-view of 34 ′ × 27 ′ and a pixel scale of 0.202 ′′ per pixel for the 8.2m Subaru telescope.The Subaru XMM-Newton Deep Survey (SXDS, Furusawa et al. 2008, Figure 1) is centered on (02 h 18 m 00 s , −05 • 00 ′ 00 ′′ ).It consists of five contiguous subfields SXDS-C, N, S, E, and W (hereafter SXDS1, 2, 3, 4, and 5) corresponding to five Suprime-Cam pointings.The total area coverage is ∼ 1.2 deg 2 .The Subaru deep field (SDF, Kashikawa et al. 2004, Figure 1) is centered on (13 h 24 m 38.s 9, +27 arcmin 2 .Suprime-Cam has taken deep images of the SXDS and SDF fields in a series of broad and narrow bands.These images have been widely used to search for high-redshift LBGs and LAEs (e.g., Yoshida et al. 2006;Ouchi et al. 2008;Kashikawa et al. 2006Kashikawa et al. , 2011;;Konno et al. 2014Konno et al. , 2016)).We retrieved the raw images from the archival server SMOKA (Baba et al. 2002).The images were reduced, re-sampled, and co-added using a combination of the Suprime-Cam Deep Field REDuction package (Yagi et al. 2002) and the IDL routines by Jiang et al. (2013).Source extraction and astrometric and photometric solutions were then applied.The details are given in Jiang et al. (2013).We performed aperture photometry with SExtractor (Bertin & Arnouts 1996) in dual-image mode using the narrowband images as the detection images.Aperture photometry was measured in a 2 ′′ diameter aperture and an aperture correction was then determined from a large number of bright but unsaturated point sources and applied to correct for light loss.The depths (5σ in a 2 ′′ diameter aperture) of the imaging data in five broad bands BV Ri ′ z ′ reach 27.9, 27.6,27.4,27.4,and 26.2 AB mag in SXDS,and 28.0,27.2,28.0,27.8,and 26.8 AB mag in SDF.The typical PSF FWHM of the images is 0.6 ′′ in the i ′ band.Galactic extinction was corrected using the E(B − V ) values at the center of SXDS and SDF, respectively.
Narrowband imaging data were reduced in the same method.Figure 2 shows the transmission curves of the narrow bands used in this paper.No strong OH sky line lies within the passbands of these narrowband fil-ters.NB570 (λ c = 5703 Å; FWHM = 68 Å) is used to select z ∼ 3.7 LAE candidates in SXDS.The depth of the NB570 image is about 24.8 mag and it slightly varies across the five SXDS subfields (±0.2 mag).NB704 (λ c = 7042 Å; FWHM = 99 Å) and NB711 (λ c = 7120 Å; FWHM = 73 Å) are used to select z ∼ 4.8 LAE candidates in SDF.The depths of the NB704 and NB711 images in SDF are 26.2 mag and 25.5 mag, respectively.A smaller region in the north of SDF (hereafter SDFn) also has R, i ′ , NB704, and NB711-band images with depths of 26.9, 26.5, 26.1, and 25.3 mag, respectively, so z ∼ 4.8 LAE candidates are also selected in this field.

Target selection in SXDS and SDF
Based on the filter curves shown in Figure 2 and the positions of the redshifted Lyα emission lines, we apply the following criteria to select LAE candidates at z ∼ 3.7 and z ∼ 4.8.These criteria roughly select LAEs with Lyα rest-frame equivalent width EW 0 ≳ 20 Å.
For z ∼ 3.7 LAE candidates in SXDS, the Lyα line locates around 5700 Å, between the V and R bands.The selection criteria are as follows.
In the above selection procedure, we used slightly different narrowband detection limits (7σ or 8σ) for different fields.This is to optimize the target surface densities for follow-up spectroscopy, because different narrowband images have different depths.The selected LAE candidates were visually inspected to exclude spurious detections, such as those near the edges of the narrowband images and those contaminated by nearby bright stars.Finally, we selected 112 z ∼ 3.7 LAE candidates in the NB570 image of SXDS, 123 z ∼ 4.8 LAE candidates in the NB704 images of SDF (87) and SDFn (36), and 28 z ∼ 4.8 LAE candidates in the NB711 images of SDF (16) and SDFn (12).

Spectroscopic observations and data reduction
We carried out follow-up spectroscopic observations of these LAE candidates from 2017 to 2021, using the optical fiber-fed spectrograph Hectospec on the 6.5 m telescope MMT (Fabricant et al. 2005).Hectospec has a large field-of-view of 1 • in diameter with 300 fibers.The diameter of each fiber is 1. ′′ 5 and adjacent fibers can be spaced as closely as 20 ′′ .The observations are summarized in Table 1.
We observed the z ∼ 3.7 LAE candidates in SXDS using four Hectospec pointings with two different configurations.This is because many fibers were shared by other projects.In 2017, we used the 600 lines mm −1 grating blazed at ∼ 6000 Å, providing a spectral resolu-tion of ∼ 2.1 Å and a wavelength coverage from 4050 Å to 6550 Å.In 2020 and 2021, we used the 270 lines mm −1 grating blazed at ∼ 5000 Å, providing a spectral resolution of ∼ 4.8 Å and a wavelength coverage from 3650 Å to 9200 Å.For the z ∼ 4.8 LAE candidates in SDF and SDFn, one Hectospec pointing was use to cover most of the sky area.They were observed in 2018 and 2019 with the 600 lines mm −1 grating that covered a wavelength coverage from 6050 Å to 8550 Å.The seeing of these observations varied from ∼ 0.7 ′′ to ∼ 1.5 ′′ .For each pointing, about 80 fibers were assigned to blank sky regions for background subtraction.
All z ∼ 4.8 LAE candidates and all but one z ∼ 3.7 LAE candidates were covered by the above Hectospec pointings, and the only one z ∼ 3.7 LAE candidate out of the pointings is excluded in the following analyses.The effective areas of the SXDS sub-fields SXDS1, 2, 3, 4, 5 covered by the pointings are 0.232, 0.235, 0.235, 0.167, and 0.181 deg 2 , respectively.The total effective area in SXDS is 1.050 deg 2 .The effective areas in SDF and SDFn covered by the pointings are 0.279 and 0.125 deg 2 , respectively, with a total area of 0.404 deg 2 .Due to the fiber collision, not all candidates were spectroscopically observed.In summary, 99 out of 112 z ∼ 3.7 LAE candidates, 96 out of 123 z ∼ 4.8 LAE candidates in NB704, and 12 out of 28 z ∼ 4.8 LAE candidates in NB711 were spectroscopically observed.The total exposure time for each target varies from 1 h to 7 h.This complexity is partly due to the fact that the fibers were shared by different programs, as mentioned earlier.However, our observing strategy ensures that fainter LAE candidates were assigned with longer exposure time so that their Lyα emission lines can be identified if they are real LAEs.
The Hectospec data were processed and reduced using the HSRED 1 reduction pipeline.For each exposure, science and lamp images were bias subtracted, flat-fielded, and cosmic ray rejected.Each fiber was traced, and each one-dimensional (1D) spectrum was extracted.Wavelength solution was then derived.Sky spectrum from each sky fiber was checked and poor sky spectra were rejected.An average "supersky" spectrum was obtained from good sky spectra, which was then scaled according to the strength of skylines in individual object's spectrum and subtracted from it.The resultant 1D, skysubtracted, and wavelength-calibrated spectra for different exposures of the same object were then re-sampled to the same wavelength grid and co-added with inverse variance weighting to achieve the final 1D spectrum.

LAE SAMPLES AND THEIR PROPERTIES AT
Z ∼ 3.7 AND Z ∼ 4.8 In this section, we construct the spectroscopic samples of LAEs at z ∼ 3.7 and z ∼ 4.8.We then derive their spectral properties, including redshift, UV continuum flux, and Lyα line flux and EW.

LAE samples
We identify LAEs as follows.For each 1D spectrum, we first search for a possible emission line in the wavelength range covered by the narrowband filters.If a line with signal-to-noise ratio SNR> 5 is detected, this object is regarded as a possible LAE.Low-redshift interlopers are mostly likely [O II] λλ3727, 3729, [O III] λλ4959, 5007, or Hβ emitters.Hα emitters are also possible for z ∼ 4.8 candidates.These contaminants are identified and excluded using the following steps.If an emission line is [O III], Hβ, or Hα, the spectrum would cover more than one of these lines, and such an interloper can be easily identified.If an emission line is [O II] and if the spectrum was taken by the 270 lines mm −1 grating, the spectrum would cover [O III], Hβ, and Hα, and thus the line is easy to identify.If the spectrum was taken by the 600 lines mm −1 grating, the wavelength No  coverage is short, but the spectral resolution is higher enough to identify the [O II] doublet.In rare cases, a Lyα line profile can also exhibit a double-peak feature (e.g., Verhamme et al. 2006), which mimics the [O II] doublet if its SNR is low.In this case, Lyα and [O II] can be distinguished using the broadband photometry, given the fact the two redshifts would be very different.If the SNR is high, they can be immediately distinguished by the line profile.
Active galactic nuclei (AGNs) are also identified and excluded.For the z ∼ 3.7 LAE candidates in SXDS, deep X-ray images taken by Chandra (Kocevski et al. 2018) and XMM-Newton (Chen et al. 2018) are available.If an object is detected in X-ray, it is considered as an AGN.In addition, if an object has broad emission lines (line FWHM greater than 1000 km s −1 ) in the spectra, it is also considered as an AGN.One Xray AGN and one broad-line AGN are identified in our sample.For the z ∼ 4.8 LAE candidates in SDF and SDFn, there are no X-ray data available, and AGNs are identified based on broad emission lines in the spectra.We find 2 AGNs in our z ∼ 4.8 sample.Thus, the AGN fractions in narrowband selected, LAE photometric samples at both z ∼ 3.7 (2/99) and z ∼ 4.8 (2/(96 + 12) = 2/108) are around 2%, which is consistent with the value reported in previous studies at similar redshifts (e.g., Ouchi et al. 2008).Interestingly, we find that one z ∼ 4.8 LAE candidate (13 h 25 m 32.s 73, +27 • 41 ′ 34.′′ 4) is a supernova (SN) happened in 2001 in a host galaxy at z = 0.515 (see Figure 7 in Morokuma et al. (2010)).This object was selected as a LAE candidate because the narrowband image (NB711) was taken in 2001 about one month after the maximum light of the SN, while the broadband images were all taken after 2002 so that the broadband magnitudes were fainter.
From the above procedure, 71 z ∼ 3.7 LAEs and 69 z ∼ 4.8 LAEs (including 63 NB704-selected LAEs and 6 NB711-selected LAEs) are spectroscopically confirmed.In the remaining 28 (39) objects at z ∼ 3.7 (4.8), we find 1 (3) low-redshift [O II] emitters, 10 (2) low-redshift [O III] emitters, and 2 (2) AGNs.In addition, 15 (32) objects do not have SNR> 5 emission lines in the spec-tra, and they are likely red stars/galaxies, transients, or spurious detections.In the following analysis, only the confirmed LAEs are considered.Their spatial distributions are shown as the blue and red points in Figure 1.Figures 3 presents the stamp images and spectra of the confirmed z ∼ 3.7 LAEs.Figures 4 presents the stamp images and spectra of the confirmed z ∼ 4.8 LAEs.Properties of these LAEs are listed in Table 2  and Table 3.

Lyα redshifts
We calculate the redshifts of the LAEs using the Lyα emission lines.For each LAE, its redshift is determined by fitting a composite Lyα line profile to the Lyα line in its spectrum.The composite Lyα line profile is obtained as follows.We first assume that the peak of the Lyα line is at 1215.67 Å in the rest frame, and transform all spectra into the rest frame.We then take the average of all spectra to construct a composite Lyα line profile.When we fit the composite line profile to the individual Lyα lines, we vary its redshift and amplitude.After obtaining new redshifts from the fitting results, we use them to transform the spectra into the rest frame again.We repeat the above procedure a few times.The final products include the average Lyα line profiles at z ∼ 3.7 and z ∼ 4.8 and redshifts for all our LAEs.Because the Lyα line is typically redshifted compared with other strong emission lines (e.g., Shapley et al. 2003;Matthee et al. 2021), possibly due to the back scattering of Lyα photons in outflowing gas (e.g., Verhamme et al. 2006), Lyα redshifts derived from the Lyα line are higher than systemic redshifts.
Figure 5 shows the redshift distributions of the LAEs at z ∼ 3.7 and z ∼ 4.8.We can see from the figure that, although the narrowband filter response curves are symmetric, the redshift distributions are not.There are more LAEs at lower redshifts.The reason is that the UV continuum at the blue side of Lyα line is much weaker than that at the red side due to the absorption of λ rest < 1215.67Å photons by ISM and IGM.Therefore, a relatively lower-redshift LAE tends to be brighter in the narrowband and has a higher possibility to be selected as an LAE candidate.

Lyα line and UV continuum flux
Because the rest-frame UV continuum of these LAEs is too weak to be detected in our spectra, we calculate the Lyα line flux and UV continuum flux based on the redshifts and deep broadband and narrowband photometry (Jiang et al. 2013).The UV continuum of star-forming galaxies can typically be represented by a power law (e.g., Meurer et al. 1995;Bouwens et al. 2009Bouwens et al. , 2014)).In addition, there is no strong emission line other than Lyα in the considered wavelength range (e.g., see a composite spectrum in Shapley et al. 2003).Therefore, we can model the UV spectrum of an LAE using where f λ (in unit of erg s −1 cm −2 Å−1 ) is the sum of the Lyα emission line flux (A × S Lyα ) and the powerlaw continuum (B × λ β ).Here S Lyα is the average Lyα line profile from Section 3.2, β is the UV continuum slope, and A and B are two scaling factors in unit of erg s −1 cm −2 Å−1 .The Lyα profile has negligible impact on the flux estimation because the broad and narrow bands are much wider than typical Lyα line widths.
For the wavelength range blueward of Lyα, we apply an average IGM absorption to the UV continuum according to Madau (1995) (the composite Lyα line has already taken into account the IGM absorption).We then fit the model spectra to the observed magnitudes and obtain A, B, and β.
For the z ∼ 3.7 LAEs, the R, i ′ , and z ′ bands do not cover their Lyα emission lines.Because f λ ∝ λ β means AB magnitude m AB ∝ (β + 2) × log(λ), we first use the R, i ′ , and z ′ magnitudes to do a linear fit and obtain B and β.In rare cases that an object is not detected in a band, we use a 2σ upper limit for this band in the fitting procedure.We then apply the IGM absorption to the model spectrum and convolve the spectrum with the narrowband transmission curve.Finally, we calculate A by matching the narrowband magnitude from the model spectrum and the observed magnitude.In Figure 6 we present a model fit to a z ∼ 3.7 LAE and illustrate the procedure.
For the z ∼ 4.8 LAEs in SDF and SDFn, the z ′ -band is the only available band redward of Lyα that does not cover the Lyα line, so we are not able to directly calculate continuum slopes.Therefore, we adopt an average UV slope β = −2.17according to Hashimoto et al. (2017) for the z ∼ 4.8 LAEs.We calculate A and B using the i ′ and narrowband magnitudes by comparing   the observed magnitudes and the magnitudes measured from the model spectrum.A 2σ upper limit is used for non-detections, as we did above.
Based on the fitting results, the UV magnitude at 1500 Å, Lyα line flux, Lyα luminosity, and the restframe Lyα EW (EW 0 ) are measured.Tables 2 and 3 present the results for the z ∼ 3.7 and z ∼ 4.8 LAEs, respectively.The uncertainties are estimated using the Monte Carlo method.For each LAE, we simulate 10000 mock sources whose broadband and narrowband magnitudes follow the observed magnitudes and error distributions.The 1σ uncertainties of the UV continuum and Lyα line quantities are determined from the 16 th and 84 th percentiles of the cumulative distributions.Photometry of a few LAEs is severely affected by nearby objects, so their properties are not measured and they are excluded in our analyses.Figure 7 shows the distribution of the UV slopes β for the z ∼ 3.7 LAEs.Our selection criteria are mainly based on the Lyα EW but not sensitive to β, so the β distribution suffers little selection effect and can be regarded as an intrinsic distribution for LAEs at this redshift.By fitting a Gaus- sian function to the β distribution, we obtain an average β = −1.72 with a scatter of 0.80.

LYα LUMINOSITY FUNCTION
In this section, we calculate Lyα LFs from the z ∼ 3.7 and z ∼ 4.8 LAEs in our samples.We first estimate sample completeness and correct for selection effects.We then derive the Lyα LFs using the 1/V a method and the maximum likelihood method.

Sample completeness
Incompleteness was introduced into our samples from several sources.The first one is the object detection in the images.Brighter LAEs in the narrow band have higher probabilities to be detected, while some faint LAEs can be missed in this step.We conduct Monte Carlo simulations to estimate object detection probabilities as a function of the narrowband magnitude.We first co-add narrowband images of the spectroscopically confirmed LAEs to obtain typical LAE morphology and sizes in different subfields.We then generate 2500 mock LAEs with different magnitudes that cover our samples, and randomly distribute them into the narrowband images.Finally, we detect these mock LAEs using SExtractor with the same configuration parameters as we did previously, and estimate the detection rate as a function of magnitude.Figures 8 shows the results for all subfields.These results will be used to correct detection incompleteness later.
The second source of the sample incompleteness comes from our LAE candidate selection criteria.For example, the color criterion V R − NB570 > 0.9 for z ∼ 3.7 LAE candidates roughly corresponds to the rest-frame Lyα EW 0 ≳ 20 Å, so LAEs with lower EW 0 have lower probabilities be selected.In addition, LAEs with relatively low Lyα luminosities and weak UV continuum can have high EW 0 , but they can be faint in the narrow band and will not satisfy our magnitude cut NB > 7σ.
In narrowband surveys so far, very few LAEs with Lyα luminosities lower than ∼ 10 41.5 erg/s are found even with deep images (Cassata et al. 2011).We carry out Monte Carlo simulations to estimate the selection completeness.We generate a grid of 10000 mock LAE spectra with different Lyα luminosities and redshifts (log(L Lyα ), z) according to equation ( 1).The UV slopes β and EW 0 are randomly chosen from their intrinsic distributions.For z ∼ 3.7 LAEs, we use the β distribution from our sample (Gaussian distribution with β = −1.72 ± 0.80) to represent its intrinsic distribution.Note that β is not sensitive to the selection criteria.For z ∼ 4.8 LAEs, we are not able to calculate β from our data, so we use β = −2.17± 1.57 from Hashimoto et al. (2017).For EW 0 , its distribution can often be fitted by an exponential distribution, N = N 0 × exp(−w/w 0 ) (e.g., Gronwall et al. 2007) with a scale length w 0 .It is usually difficult to obtain the intrinsic EW 0 distribution.We use the results from Hashimoto et al. ( 2017) Liu et al. (w 0 = 113 Å and w 0 = 68 Å for z ∼ 3.7 and z ∼ 4.8 LAEs, respectively), whose LAE sample is not selected by the narrowband technique but by a blind Integral Field Unit (IFU) spectroscopic survey.Their sample reaches a great depth for Lyα luminosity without strong selection effects on EW 0 , so it can be used to represent the intrinsic EW 0 distribution.After applying the IGM absorption to the mock spectra, we convolve them with broad-and narrowband filter transmission curves to obtain their magnitudes as if they were actually observed.Errors for magnitudes are estimated based on the magnitude-error relation of each band in each subfield.Then, we select these mock LAEs with our selection criteria and calculate the selection rate as a function of Lyα luminosity and redshift.
Third, our spectroscopic observations and target identification also brought incompleteness to our samples.As previously mentioned, due to the fiber collision and other constraints, 99 out of 111 LAE candidates at z ∼ 3.7 and 108 out of 151 LAE candidates at z ∼ 4.8 were assigned with fibers.This incompleteness is independently of luminosity.As we also mentioned earlier, our strategy ensures that each real LAE gets enough exposure time for line identification, so we assume that target identification did not introduce any incompleteness.Figure 9 shows the total completeness (combination of object detection, target selection, and spectroscopic observations) in different subfields for the z ∼ 3.7 and z ∼ 4.8 LAE samples.We do not apply any completeness cut when we calculate the Lyα LF in the next two subsections.

The 1/V a method
We calculate the non-parametric binned Lyα LF of our z ∼ 3.7 and z ∼ 4.8 LAEs using the 1/V a method (Avni & Bahcall 1980), a modified version of the 1/V max method.V a is the available comoving volume that a LAE with Lyα luminosity L can be detected in our sample.If we use p(log L, z) to denote the sample completeness as a function of Lyα luminosity and redshift, then V a can be expressed as where ω is the area of the field, z min and z max are determined by the narrowband filters, and V c is comoving volume.The binned LF is given by where V a,i is V a for the i th LAE in the sample and ∆ log L is the bin size.The 1σ error for each bin is given by the Poisson statistics, If several subfields are considered, equation 2 can be written as where ω j is the area of the j th subfield, and p j (log L, z) is the completeness for the j th subfield.The effective areas for all subfields are given in Section 2.

Maximum likelihood method
To parameterize the Lyα LF, we use the Schechter function (Schechter 1976), where ϕ * is the characteristic volume density, L * is the characteristic luminosity, and α is the faint-end slope.We do not fit the Schechter function to the binned Lyα LF because the chosen size of the luminosity bin may potentially affect the fitting results.Instead, we adopt a maximum likelihood method (e.g., Ciardullo et al. 2013) to get the best-fit Schechter function parameters.
Given (ϕ * Lyα , L * Lyα , α), the expected number of LAEs with L < L(Lyα) < L + ∆L that should be discovered in our sample is denoted by λ = f (L)∆L, where The actual number n of LAEs found with L < L(Lyα) < L + ∆L follows Poisson statistics with P (n) = e −λ λ n n! .If we reduce ∆L enough so that in each luminosity bin either 0 LAE or 1 LAE is found, since P (1) = λe −λ and , the likelihood function based on our sample can be written as where L i is the Lyα luminosity for each confirmed LAE and L k refers to luminosity bin without any discovered LAE.After taking small enough ∆L, the logarithm of the likelihood function can be expressed by where i represents the i th LAE and j represents the j th subfield.In this formula, all three Schechter function parameters ϕ * Lyα , L * Lyα , and α are free parameters.By varying (ϕ * Lyα , L * Lyα , α), we find the best-fit parameters that maximize ln P , or minimize S = −2 × ln P .
Because our samples are not deep enough to constrain the faint-end slope α of the Lyα LF, we fix α from −1.5 to −2.0 in a step of 0.1 and then fit the other two parameters.Table 4 shows the best-fit results.The 1σ error for each parameter is determined by marginalizing over the other parameter and finding the 16 th and 84 th percentiles of the cumulative distribution.In Figures 10   and 11, we show the Lyα LFs with α = −1.5 that was commonly used in previous studies (e.g., Gronwall et al. 2007;Ouchi et al. 2008;Ning et al. 2022).The best-fit Schechter functions are consistent with the binned LF, which suggests that the Lyα LF of LAEs is approximated well by the Schechter function.
In the above procedure, we did not apply any luminosity cut when we derived the LFs.In addition, the Eddington bias is small compared to the uncertainties of the LFs.We estimate the effect of the Eddington bias using Monte Carlo simulations below.For each LAE, we randomly choose a Lyα luminosity that follows the observed Gaussian distribution of its Lyα luminosity.We then derive the LFs at the two redshifts with α fixed to −1.5, as we did earlier.We repeat the simulation 100 times and find that the LF results are all consistent with the previous results within one sigma uncertainties.Therefore, we did not consider the Eddington bias above.

Comparison with previous studies
In Figures 10 and 11, Lyα LFs at similar redshifts from the literature are overploted.These studies use different methods (including the narrowband technique, IFU spectroscopy, and slit spectroscopy) to obtain LAE samples, so their sample sizes, survey volumes, redshift ranges, and Lyα luminosity limits are different.Nevertheless, our Lyα LFs at the two redshifts are roughly consistent with most of the previous measurements.
We first compare Lyα LFs at z ∼ 3.7.Ouchi et al. (2008) used the same narrowband filter NB570 and selected LAEs in the same field SXDS as we did.They measured their Lyα LF based on the photometric sample (they fixed α = −1.5).Their LF at the faint end agrees well with our faint-end Lyα LF, but their brightend LF is a factor of 2 − 3 higher than ours.The reason is unclear, since bright LAEs are relatively easy to detect and identify, with little selection bias or incompleteness.On the other hand, our spectroscopically confirmed sample ruled out contaminants that may exist in previous photometric samples.The Lyα LF from Sobral et al. ( 2018) is systematically lower that ours from the faint end to the bright end, possibly because they used a medium band (FWHM ∼ 280 Å) to select LAE candidates with large Lyα EWs (EW 0 > 50 Å).The Cassata et al. (2011) LAE sample was obtained using slit spectroscopy and the Drake et al. (2017) sample was obtained using MUSE IFU spectroscopy.Their Lyα LFs are not consistent with each other, and both are much lower than ours at the bright end.It is likely because their survey volumes are relatively small (≲ 10 5 Mpc 3 ) and thus there are not enough bright LAEs.At the faint end, Lyα LFs from Drake et al. (2017) and Herenz et al. (2019) are higher than our result.As Herenz et al. (2019) pointed out, LAEs exhibit diffuse extended low surface-brightness halos.If taking this into consideration in the calculation of the Lyα luminosity, Herenz et al. (2019) found that the Lyα LFs at the faint end are underestimated in the narrowband studies.We, like previous narrowband studies, did not take this into account.The Zhang et al. (2021) LAE sample is from a large HETDEX IFU spectroscopic survey of bright LAEs, and their Lyα LF at z ∼ 2 − 3.5 is consistent with (or slightly below) our LF.The consistency may reflect the fact that there is a moderate increase of Lyα LF from z ∼ 2 to 3 and almost no evolution from z ∼ 3 to 6 (e.g., Konno et al. 2016).
We then compare z ∼ 4.8 Lyα LFs.We first compare our LF with Shioya et al. (2009) who used the same narrowband filter NB711 as we did to obtain a photometric LAE sample in the COSMOS field.Their Lyα LF is similar to ours at the bright end, but significantly smaller at the faint end.The reason is likely that Shioya et al. (2009) derived their Lyα LF without considering completeness correction, because a sample suffers larger incompletenesses at fainter luminosities.Our z ∼ 4.8 LF agrees well with the z ∼ 4.5 Lyα LF in Zheng et al. (2013), which used a spectroscopically confirmed LAE sample and incorporated the LAE sample from Dawson et al. (2007).The Lyα LF by Drake et al. (2017) based on the MUSE data is smaller at the bright end due to their smaller survey volume.The Lyα LFs by Cassata et al. (2011) based on the slit spectroscopy and by Herenz et al. (2019) based on the MUSE data are similar to our LF at the bright end (because at z ∼ 4.8 our survey volume is comparable to theirs) and higher than our Lyα LF at the faint end.This is likely because their slit spectroscopy and IFU spectroscopy data probed fainter Lyα luminosities.In addition, Herenz et al. (2019) took diffuse low surface brightness halos into consideration, as we explained earlier.

Evolution of the Lyα LF
We further study the evolution of the LF parameters.In Figure 12, we plot the likelihood contours in the ϕ * Lyα − L * Lyα space derived from the z ∼ 3.7 and 4.8 LFs with fixed α = −1.5.The 1σ (68%) and 2σ (95%) confidence regions are determined by respectively setting ∆S = 2.30 and ∆S = 6.17 when fitting two free parameters (ϕ * Lyα and L * Lyα ) in the maximum likelihood method.Although the 1σ contours barely overlap, the 2σ contours overlap with each other.The best-fit ϕ * at the two redshifts are almost the same and the best-fit L * shows a small decrease (∼ 0.1 dex) from z ∼ 3.7 to z ∼ 4.8.Therefore, the Lyα LF evolves little from z ∼ 3.7 to 4.8.
In Figure 13, we plot the Lyα LFs at several redshifts from z ∼ 3.1 to z ∼ 6.6, including results from our work, Guo et al. (2020), Ning et al. (2022), and Zheng et al. (in preparation, hereafter Z23).Guo et al. (2020) calculated a Lyα LF at z ∼ 3.1 using a large spectroscopic sample of 166 LAEs over ∼ 1.2 deg 2 .Ning et al. (2022) derived a Lyα LF at z ∼ 6.6 based on 36 spectroscopically confirmed LAEs over ∼ 2 deg 2 .Z23 calculated a Lyα LF at z ∼ 5.7 based on a large sample of 260 spectroscopically confirmed LAEs constructed by Ning et al. (2020).The imaging data used in these studies were reduced in the same manner as in our work (Jiang et al. 2013), the LAEs were selected using the narrowband technique in a consistent way, the areas covered are all large enough to avoid the effect of cosmic variance, and the Lyα LFs were all based on spectroscopically confirmed LAE samples.Therefore, the evolution of the Lyα LF can be measured reliably from these studies.All these studies are based on narrowband selected LAE samples, and are not deep enough to reach faint Lyα luminosities, so α was fixed to about -1.5 in these studies.Figure 13 shows that the Lyα LF evolves slowly from z ∼ 3.1 to 5.7, with a factor of ∼ 3 decrease in ϕ * Lyα and a factor of ∼ 1.5 decrease in L * Lyα .The similar conclusion that the Lyα LF has no significant evolution from z ∼ 3 to z ∼ 6 is also supported by Ouchi et al. (2008), Cassata et al. (2011), Drake et al. (2017), and Herenz et al. (2019).
In contrary to the slow evolution from z ∼ 3.1 to 5.7, the Lyα LF declines rapidly from z ∼ 5.7 to z ∼ 6.6 at the faint end.Such a decrease can be explained (at least partly) by cosmic reionization.Cosmic reionization is the phase transition of the IGM when the neutral hydrogen (HI) was gradually ionized by ionizing photons at z ≥ 6 (Fan et al. 2006).The universe became largely ionized after reionization.Therefore, Lyα photons emitted by LAEs at z ∼ 6.6 are more absorbed/scattered (compared with z ∼ 5.7 LAEs) in the IGM, and thus LAEs at z ∼ 6.6 have fainter Lyα luminosities, especially at the faint end.At the bright end, however, the Lyα LF at z ∼ 6.6 seems unaffected by the IGM absorption, probably due to the existence of large ionized bubbles around luminous LAEs that allow Lyα photons to escape (e.g., Ning et al. 2022).

SUMMARY
We have presented one of the largest spectroscopically confirmed sample of LAEs at z ∼ 3.7 and z ∼ 4.8.Using the narrowband technique, we selected LAE candidates with deep broadband and narrowband images taken by the Subaru Suprime-Cam: 112 z ∼ 3.7 LAE candidates in SXDS were selected in the narrowband NB570, and 151 z ∼ 4.8 LAE candidates in SDF and SDFn were selected in the narrowbands NB704 and NB711.Based on the spectra taken by the MMT Hectospec spectrograph, we finally confirmed 71 z ∼ 3.7 LAEs and 69 z ∼ 4.8 LAEs.
We determined the Lyα redshifts of the LAEs by fitting a composite Lyα line profile to individual LAE spectra.From secure redshifts and deep broadband and narrowband photometry, we calculated Lyα luminosity, Lyα EW 0 , and M UV for each LAE.The z ∼ 3.7 LAEs and the z ∼ 4.8 LAEs span Lyα luminosity ranges of ∼ 10 42.6 − 10 43.7 erg s −1 and ∼ 10 42.4 − 10 43.4 erg s −1 , respectively.They represent the most Lyα-luminous galaxies at the two redshifts.
We have derived Lyα LFs at z ∼ 3.7 and 4.8 based on the two LAE samples, after considering sample incompletenesses.The binned Lyα LFs are well approximated by the Schechter function.We determined the best-fit ϕ * Lyα and L * Lyα by varying the slope α from -1.5 to -2.0 in the maximum likelihood method.The Lyα LFs show little evolution between the two redshifts.Our Lyα LFs are broadly consistent with previous measurements within a factor of 2 − 3.By comparing with Lyα LFs at other redshifts in the literature, we found that Lyα LFs have no significant evolution from z ∼ 3.1 to z ∼ 5.7, but decline rapidly from z ∼ 5.7 to z ∼ 6.6 at the faint end.2013,2018,2022), SExtractor (Bertin & Arnouts 1996), HSRED

Figure 1 .
Figure 1.LAEs in the SXDS, SDF, and SDFn fields.The left panel shows LAEs at z ∼ 3.7 in SXDS and the right panel shows LAEs at z ∼ 4.8 in SDF and SDFn.In each panel, the grey regions indicate the coverage area of the Subaru Suprime-Cam images.The large circles represent the pointings of our Hectospec observations.The color-coded points represent spectroscopically confirmed LAEs.

Figure 3 .
Figure 3. Example of 10 LAEs spectroscopically identified at z ∼ 3.7.The left panel shows the stamp images in five broad bands and the narrowband NB570.The image size is 6 ′′ × 6 ′′ .The right panel shows the optical spectra of the LAEs taken by MMT Hectospec.The red vertical lines indicate the positions of the redshifted Lyα lines.A figure of the full version is shown in the Appendix.

Figure 4 .
Figure 4. Same as Figure 3, but for our z ∼ 4.8 LAEs selected in NB704.A figure of the full version is shown in the Appendix.

Figure 5 .
Figure 5. Redshift distribution of the spectroscopically confirmed LAEs at z ∼ 3.7 (upper panel) and z ∼ 4.8 (lower panel).The dashed lines show the response curves of the narrowband filters.

Figure 6 .
Figure 6.Example to illustrate our procedure to model an LAE spectrum and measure its Lyα and UV continuum properties.The red circles are photometric data points.The vertical error bars indicate the photometry errors and the horizontal bars indicate the FWHM of the narrow-and broadband filters.This LAE is at z = 3.673.Its power-law continuum is fitted with the broadband R, i ′ , and z ′ photometry.After applying the IGM absorption, we calculate the scaling factor A in equation (1) by matching the Lyα line profile with the narrowband photometry.The inset shows the region around the Lyα line.

Figure 7 .
Figure 7. Distribution of the UV slope β for z ∼ 3.7 LAE sample.The dashed line represents the best-fit Gaussian distribution with an average β = −1.72 and σ β = 0.80.

Figure 8 .
Figure 8. Source detection rates in the 5 SXDS subfields (upper panel) and SDF and SDFn (lower panel) as a function of the narrowband magnitudes.Differences between the curves reflect different depths of the narrowband images.

Figure 9 .
Figure 9.The total completeness fraction (combining source detection, target selection, and spectroscopic observations) in different subfields as a function of Lyα luminosity and redshift.The three contours represents the fractions of 10%, 50%, and 70% in the upper panel, and 10%, 30%, and 50% in the lower panel.The red points with error bars represent the LAEs in our samples.
3. The maximum available comoving volume is ∼ 6 × 10 5 Mpc 3 for our z ∼ 3.7 LAE sample and ∼ 4 × 10 5 Mpc 3 for our z ∼ 4.8 LAE sample.The red points in Figures 10 and 11 show the binned Lyα LF.The vertical error bars represent Poisson errors and the horizontal bars indicate bin sizes.
Astropy (Astropy Collaboration et al. • 29 ′ 25. ′′ 9) and covers an area of ∼ 800 Response curves of the narrowband filters (black solid lines) and broadband filters (dashed lines) used in this work.Instrument response has been included for the broad bands.The maximum responses of the narrowband filters are set to 0.8 for the purpose of clarity.The filter NB570 is used to select z ∼ 3.7 LAEs, and the NB704 and NB711 filters are used to select z ∼ 4.8 LAEs.

Table 1 .
Summary of the MMT Hectospec observations.
Figure 11.Lyα LFs at z ∼ 4.8.The red points with error bars represent the binned LF of the z ∼ 4.8 LAE sample.The vertical error bars denote the 1σ Poisson error and the horizontal bars indicate the bin size of ∆ log L = 0.2 dex.The black solid line is the best-fit LF with the Schechter function and α fixed to -1.5 for the z ∼ 4.8 LAE sample.

Table 4 .
Best-fit parameters of the Schechter function for the z ∼ 3.7 and 4.8 Lyα LFs.