The Magellan M2FS Spectroscopic Survey of High-redshift Galaxies: The Brightest Lyman-break Galaxies at z ∼ 6

We present a study of a sample of 45 spectroscopically confirmed, UV luminous galaxies at z ∼ 6. They were selected as bright Lyman-break galaxies (LBGs) using deep multiband optical images in more than 2 deg2 of the sky, and subsequently identified via their strong Lyα emission. The majority of these LBGs span an absolute UV magnitude range from −22.0 to −20.5 mag with Lyα equivalent width (EW) between ∼10 and ∼200 Å, representing the most luminous galaxies at z ∼ 6 in terms of both UV continuum emission and Lyα line emission. We model the spectral energy distributions of 10 LBGs that have deep infrared observations from Hubble Space Telescope, JWST, and/or Spitzer, and find that they have a wide range of stellar masses and ages. They also have high star formation rates ranging from a few tens to a few hundreds of solar mass per year. Five of the LBGs have JWST or HST images, and four of them show compact morphology in these images, including one that is roughly consistent with a point source, suggesting that UV luminous galaxies at this redshift are generally compact. The fraction of our photometrically selected LBGs with strong Lyα emission (EW > 25 Å) is about 0.2, which is consistent with previous results and supports a moderate evolution of the intergalactic medium opacity at the end of cosmic reionization. Using deep X-ray images, we do not find evidence of strong active galactic nucleus (AGN) activity in these galaxies, but our constraint is loose, and we are not able to rule out the possibility of any weak AGN activity.


INTRODUCTION
High-redshift (z ≥ 6) galaxies are a key probe to understand the early Universe, including the early galaxy formation and evolution, the development of large-scale structures, and the history of cosmic reionization.In the past twenty years, a number of works have been conducted to search for high-redshift galaxies, and the number of known galaxies has increased dramatically thanks to the advances of instrumentation on the Hubble Space Telescope (HST) and large ground-based telescopes (e.g., Kashikawa et al. 2006;Vanzella et al. 2009;Hu et al. 2010;Wilkins et al. 2010;Stark et al. 2011;McLeod et al. 2016;Bouwens et al. 2021;Zheng et al. 2017;Ono et al. 2018;Ning et al. 2020;Wold et al. 2022).The majority of the earlier known galaxies at z ≥ 6 are photometrically selected Lyman-break galaxies (LBGs) or LBG candidates selected by the dropout technique.The narrowband technique (or the Lyα technique) is a complementary method to find high-redshift galax-ies.Ground-based narrowband imaging surveys are efficient in finding Lyα emitting galaxies (Lyα emitters, or LAEs) at certain redshift slices such as z ≈ 5.7, 6.6, and 7.0 that correspond to wavelength windows with weak OH sky lines.
Before the launch of James Webb Space Telescope (JWST), the majority of the photometrically selected high-redshift galaxies have not been spectroscopically observed.While narrowband selected LAE candidates typically have a high success rate of spectroscopic confirmation, only a small fraction of photometrically selected LBGs are spectroscopically confirmed (e.g., Toshikawa et al. 2012;Finkelstein et al. 2013;Watson et al. 2015;Roberts-Borsani et al. 2016;Song et al. 2016;Larson et al. 2022).Spectroscopic observations are important to understand LBGs, because they can not only tell us whether an object is a real LBG, but also provide a key parameter, redshift.Large samples of highredshift galaxies have been frequently used to measure galaxy properties, such as UV slopes (e.g., Dunlop et al. 2012;Finkelstein et al. 2012a;Bouwens et al. 2014), galaxy morphology (e.g., Guaita et al. 2015;Kawamata et al. 2015;Shibuya et al. 2015Shibuya et al. , 2016;;Curtis-Lake et al. 2016;Kobayashi et al. 2016;Liu et al. 2017;Naidu et al. 2022), stellar populations, and star-formation rates (e.g., Egami et al. 2005;Stark et al. 2013;González et al. 2014;Faisst et al. 2016;Castellano et al. 2017;Karman et al. 2017).These studies are mostly based on photometric samples.
High-redshift galaxies are also thought to be the major source of ionizing photons at the epoch of cosmic reionization.Recent studies have suggested that starforming galaxies, rather than quasars or AGN, provide most of the ionizing photons needed for reionization (e.g., Finkelstein et al. 2012b;Robertson et al. 2015;Onoue et al. 2017;Jiang et al. 2022).Spectroscopic samples of high-redshift LBGs can be used to constrain the properties of the IGM during this epoch.For example, Lyα photons are easily absorbed by neutral hydrogen, making the Lyα line a sensitive probe of the ionization state of the IGM.Therefore, the fraction of LBGs that exhibit strong Lyα emission lines can trace the fraction of neutral hydrogen in the IGM at z ≥ 6.This so-called Lyα visibility test from previous works has shown that the fraction of LBGs with strong Lyα emission increases steadily from z ∼ 2 to 6 and then declines sharply towards z ∼ 7 (e.g., Stark et al. 2010Stark et al. , 2011;;Ono et al. 2012;Schenker et al. 2012;Finkelstein et al. 2013;Tilvi et al. 2014;Vanzella et al. 2014;Cassata et al. 2015;Schmidt et al. 2016;Jung et al. 2018Jung et al. , 2022)).This is broadly consistent with the evolving neutral fraction in the IGM at the end of the reionization epoch.
In this paper, we present a sample of 45 LBGs with spectroscopically confirmed Lyα emission lines at z ∼ 6 from a survey of high-redshift galaxies using the Magellan telescope (Jiang et al. 2017).In Section 2, we introduce the survey program, target selection, spectroscopic observations, and data reduction.In Section 3, we construct our spectroscopic sample of LBGs.In Section 4, we present the properties of the galaxies using spectra and images from multi-wavelength observations.We discuss our results in Section 5, and summarize our paper in Section 6.Throughout the paper, we use a standard flat cosmology with H 0 = 70 km s −1 Mpc −1 , Ω m = 0.3 and Ω Λ = 0.7.All magnitudes refer to the AB system.

TARGET SELECTION AND SPECTROSCOPIC OBSERVATIONS
2.1.A Magellan M2FS Survey of z ∼ 6 Galaxies We carried out a large spectroscopic survey of galaxies at 5.5 < z < 6.8, using the large field-of-view (FoV), fiber-fed, multi-object spectrograph M2FS (Mateo et al. 2012) on the 6.5 m Magellan Clay telescope.The overview of the survey program was provided in Jiang et al. (2017).The survey was designed to build a large and homogeneous sample of high-redshift galaxies, including LAEs at z ≈ 5.7 and 6.6, and LBGs with strong Lyα emission at 5.5 < z < 6.8.Based on this sample, we can study the properties of these galaxies, the Lyα luminosity function (LF) and its evolution at high redshift, high-redshift protoclusters, cosmic reionization, etc. Taking advantage of a 30 ′ -diameter FoV, M2FS is very efficient in identifying relatively bright high-redshift galaxies (e.g., Oyarzún et al. 2016Oyarzún et al. , 2017)).The fields that we chose to observe are well-studied deep fields, including the Subaru XMM-Newton Deep Survey (SXDS), A370, the Extended Chandra Deep Field-South (ECDFS), COSMOS, etc.They cover a total of ∼ 2 deg 2 .These fields have deep Subaru Suprime-Cam imaging data in a series of broad bands BV R(r ′ )I(i ′ )z ′ and narrow bands (e.g., NB816 and NB921).The images have been used to select LAEs at z ≈ 5.7 and 6.6 and LBGs at z ∼ 6 in the survey program.The fields are summarized in Table 1.Some of the fields were later observed by Subaru Hyper Suprime-Cam (HSC) (Aihara et al. 2022), but we did not use the HSC data by the time of our observations.The M2FS observations of the program have been completed and the data have been reduced.From the survey, Ning et al. (2020) built a large sample of 260 LAEs at z ≈ 5.7 and Jiang et al. (2018) discovered a giant protocluster at z ≈ 5.7.Ning et al. (2022) then assembled a sample of 36 LAEs at z ≈ 6.6.By comparing with the z ≈ 5.7 Lyα LF, they confirmed that there is a rapid LF evolution at the faint end, but there is a lack of evolution at the bright end.Wu et al. (2020) detected diffuse Lyα halos around galaxies at z ≈ 5.7 by stacking 310 spectroscopically confirmed LAEs.In this paper, we focus on luminous LBGs at z ∼ 6.The details of the target selection is presented in Section 2.3.

Spectroscopic Observations and Data Reduction
We carried out the M2FS observations in 2015-2018.The resolving power of the spectra is R ∼ 2000, and the wavelength coverage is roughly from 7600 to 9600 Å, corresponding to a redshift range of z = 5.3-6.8 for the detection of the Lyα emission line.The selection of M2FS pointing centers was limited by the number and spatial distribution of bright stars that were used for guiding, alignment, and primary-mirror wavefront corrections.The layout of the pointings is shown in Figure 1.The effective integration time per pointing was about 5 hr on average, and the individual exposure time was 30 min, 45 min, or 1 hr, depending on airmass and weather conditions.
We used our own customized pipeline for data reduction, including bias (overscan) correction, dark subtraction, flat-fielding, cosmic ray identification, wavelength calibration, and sky subtraction.For detailed steps, see Jiang et al. (2017).The pipeline produces combined 1D and 2D spectra as well as the 2D spectra of individual exposures.This pipeline is performed for all targets in the same fields, including the candidates of LAEs and LBGs.The final 1D and 2D spectra were used to identify Lyα emission lines.

Selection of z ∼ 6 LBG Candidates
Our M2FS observations included a sample of luminous LBG candidates at z ∼ 6. Deep Subaru Suprime-Cam images were used for the selection of these candidates.The selection was mainly based on the i − z ′ color (here i means either i ′ or I).The upper panel of Figure 2 shows the filters that were used for our target selection.The selection limit is 7σ detections in the z ′ band.We applied the following color cuts, To illustrate the evolution of the i − z ′ color with redshift, we produced a simple galaxy spectrum consisting of a power-law continuum and a Lyα emission line.The power-law continuum has a fixed slope of -2.0, and the Lyα equivalent width (EW) has two values, 20 and 200 Å, respectively.We then apply the Inoue et al. (2014) IGM attenuation models to the spectra.The result is shown in the lower panel of Figure 2. The value of i − z ′ increases rapidly at z ∼ 6.With our selection criteria, galaxies with z ≳ 5.6 can be selected.The actual selection function is more complex and depends on the galaxy continuum emission and Lyα line emission.We also required no detections (< 2σ) in any bands bluer than r (or R), assuming that no flux can be detected at the wavelength bluer than the Lyman limit.Each candidate was visually inspected.
Compared to the previous selection criteria of LBGs (e.g., Ono et al. 2018;Bouwens et al. 2021), our selection criteria are relatively conservative.Previous studies usually applied redder i−z colors (e.g., i−z > 1.5 in Ono et al. 2018), and often considered infrared photometry.For example, many contaminants of high-redshift galaxies are late-type dwarf stars or low-redshift red (dusty)   galaxies, and thus tend to have redder colors in the near-IR.We did not use near-IR data, despite the fact that some of our fields are covered by deep near-IR images.Our conservative criteria allow us to include less promising candidates and thus achieve higher sample completeness.On the other hand, it means a relatively lower efficiency, i.e., a larger fraction of contaminants.It is not a concern in our program, since we have enough fibers to cover almost all these candidates.In Table 1, Column 8 shows the numbers of the candidates in each field that were observed by our program.

Lyα Line Identification and a Sample of z ∼ 6 LBGs
In this subsection, we will describe how to identify Lyα emission lines in the spectra, and present two samples of spectroscopically confirmed LBGs, including a 'good' sample with high confidence and a 'possible' sam-ple with lower confidence.We will then derive some basic properties of the Lyα lines.
We use both 1D and 2D spectra to identify Lyα emission lines.We first carry out an initial search.For each galaxy, we bin its 1D spectrum by 3 Å and search for an emission line with signal-to-noise ratio S/N > 5 in the whole wavelength range.For a line in the 'good' sample, it needs to cover at least 5 contiguous pixels with S/N > 1 in each pixel in the binned spectrum.This is determined by the typical Lyα emission line width in our spectra.For the 'possible' sample, we allow a line to cover only 3 or 4 contiguous pixels with S/N > 1.The S/N of a line is estimated by summing up the corresponding pixels of the line in the original spectrum.The reason of adopting a criterion of S/N > 1 is that most of the regions in the LBG spectra are contaminated by strong OH sky lines.If a Lyα line is close to OH lines, it could be severely affected by sky line residuals after the sky background subtraction.5.8 6.0 6.2 6.4 6.6 6.8 redshift After the initial screening above, we visually inspect each identified emission line in the 2D spectra, including the combined spectrum and the spectra of individual exposures.Obvious false detections are removed.The Lyα emission line of a z ∼ 6 galaxy is often much broader than other emission lines at much lower redshift in the observed frame due to the (1+z) broadening.It also shows an asymmetric profile due to strong IGM absorption and ISM kinematics.Possible contaminant lines from low-redshift galaxies usually appear narrow and symmetric, including Hβ, [O III] λ5007, or Hα emission lines.In addition, our resolving power of ∼ 2000 can nearly resolve the [O II] λλ3727, 3729 doublet.To quantify the asymmetry of the Lyα line, we calculate the weighted skewness S W introduced by Shimasaku et al. (2006).This is the skewness (or third moment) of the line multiplied by (λ 10,r − λ 10,b ), where λ 10,r and λ 10,b are the wavelengths where the flux drops to 10% of its peak value at the red and blue sides of the line, respectively.Since the Lyα emission of highredshift galaxies tends to be broader than other emission lines of lower-redshift galaxies, this line width factor enhances the difference between Lyα and other lines.The Lyα S W values of the LBGs in our 'good' sample and 'possible' sample range from 2.3 to 18.4 and from 2.1 to 12.2, respectively.In comparison, we also calculate S W of 10 emission lines that are identified as low-redshift interlopers based on their double-peak (e.g., an [O II] λλ3727, 3729 doublet at z ∼ 1.3), or based on the detection of other lines at the same redshift in the spectra (e.g., an [O III] λ5007 and corresponding [O III] λ4959 at z ∼ 0.72).Their S W values span a range from −12.5 to 3.3.Usually S W > 3 has been used as a threshold to distinguish Lyα from other lines, but the asymmetric shape of a high-redshift Lyα emission line may not be obvious when its S/N is relatively low (e.g., Ning et al. 2022).For three LBGs with 2.1 < S W < 3, they were not detected in two blue bands B and V , so we put them in the 'possible' sample.
We emphasize that in our LBG candidate selection procedure addressed in Section 2.3, we required nondetections for the candidates in the B-and V -band images.These images reach a great depth of ≥ 28.5 mag.They can efficiently exclude low-redshift interlopers.
Finally, the 'good' sample of 34 LBGs has convincing, luminous Lyα emission lines with S/N > 5 and asymmetric profiles.Figure 3 shows the 1D and 2D spectra in the Lyα region of the 36 LBGs in the 'good' sample.We can see that strong emission lines clearly show asymmetric line shapes due to the IGM absorption and ISM kinematics.The 'possible' sample contains 11 LBGs with less promising Lyα emission lines.We show all the confirmed LBGs in Table 2. Column 5 lists the spectroscopic redshifts measured from the Lyα lines, and their errors are smaller than 0.001.Figure 1 illustrates the positions of the targets in the five fields, including the observed candidates (all points) and the spectroscopically confirmed LBGs (color-coded points).The large circles represent the M2FS pointings.Despite the fact that the exposure time and depth of the individual pointings are similar, the numbers of the confirmed LBGs in these pointings are quite different, suggesting the existence of significant cosmic variance.As mentioned in Ning et al. (2020Ning et al. ( , 2022)), COSMOS1 and COS-MOS3 have some alignment problems during the observations, so they only have a few LBGs confirmed.

Properties of the Lyα Emission Lines
We measure redshifts and Lyα line flux for the LBGs using a Lyα profile template.The template is a composite Lyα spectrum of a large number of LAEs at z ≈ 5.7 (Ning et al. 2020).For each LBG, we first estimate an initial redshift using the wavelength of the Lyα line peak (1215.67Å).We then refine this redshift and model the Lyα profile as follows.From the Lyα template, we generate a set of model spectra for a grid of peak values, line widths, and redshifts.The peak values vary from  2. 0.9 to 1.1 (times of the observed peak value) with a step size of 0.01.The line widths are from 0.5 to 2.0 times the original width with a step size of 0.1.The redshift values vary within the initial redshift ±0.002 with a step size of 0.0001.We find the best-fit model for each LBG using the χ 2 method.The best-fit redshifts are listed in Column 5 of Table 2.
We measure the Lyα line flux of each LBG by integrating the best-fit Lyα model profile and applying a correction.Line flux directly calculated from an observed line in faint targets (like our LBGs) can be significantly underestimated due to several reasons, such as fiber loss, alignment problems, skyline contamination, etc.For narrowband-selected LAEs, their Lyα line flux can be well determined using the combination of spectroscopic redshifts, narrowband photometry, and broadband photometry.This has been applied to the LAEs in our fields (Ning et al. 2020(Ning et al. , 2022)).For 260 LAEs at z ≈ 5.7 in our fields, we compare their Lyα flux derived from the photometric data and from the best-fit Lyα line profiles, and the result is shown in Figure 4. We use a hierarchical Bayesian method to do a linear regression on the Lyα luminosities calculated by these two methods, considering the errors on both axes.The best-fit function y = 0.92x+2.32 is depicted by the black dashed line in Figure 4, and the 1σ and 3σ confidence intervals are shown in the shaded regions.Given that these LAEs and the LBGs in this work are in the same   2.The IGM absorption blueward of Lyα is not considered.We use LBGs in the 'good' sample to study the dependence of the Lyα line width (full width at half maximum, or FWHM) on redshift and Lyα luminosity.We exclude three galaxies whose Lyα lines are likely largely affected by strong residuals of the sky subtraction.We then divide the remaining galaxies into two Lyα luminosity bins, and produce a stacked spectrum for each bin.We adjust the number of spectra in each bin so that the two stacked spectra have similar S/N ∼ 50.The result is shown in the upper panel of Figure 5.We also produce two stacked spectra for two redshift bins using the same method, and the result is shown in the lower panel of Figure 5.The figure does not show apparent dependence of the line FWHM on redshift or luminosity.Using a large sample of LAEs at z ≈ 5.7 and 6.6, Ning et al. (2020Ning et al. ( , 2022) ) found that more luminous galaxies tend to have higher Lyα FWHMs, because more massive host halos possess higher HI column densities and higher gas velocities.We do not see this trend in our LBG sample, likely because the line shape has been severely affected by OH skylines.Unlike LAEs at z ≈ 5.7 and 6.6, the Lyα lines of the LBGs in our sample are often contaminated by skylines.

UV slopes and Lyα EWs
In this section, we collect multi-wavelength observations to study the physical properties and stellar populations of the LBGs in our sample.We make use of the following data, including the Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP, Aihara et al. 2022), the UKIRT Infrared Deep Sky Survey (UKIDSS, Lawrence et al. 2007), the Great Observatories Origins Deep Survey (GOODS, Giavalisco et al. 2004) Most LBGs in our sample except those in the A370a field have multi-wavelength data that can be used to determine UV continuum slope β UV and other physical properties via SED modeling.The IRAC images of many LBGs are blended with nearby sources.In this case, GALFIT (Peng et al. 2002(Peng et al. , 2010) is used to model and subtract contaminant sources.Photometry of the target is then done after the de-blending.
We then measure β UV and the absolute UV magnitude M UV (at rest-frame 1500 Å) for the LBGs using the near-IR data mentioned above.The median β UV value is −1.86, which is consistent with the median β UV = −1.73+0.14  −0.20 for z ∼ 6 LBGs within the absolute magnitude bin M UV = −20.75± 0.5 mag from Bouwens et al. (2014).We also measure the UV slopes of bright LAEs at z ≈ 5.7 in Ning et al. (2020) and find a median value of β UV ∼ −2, which is consistent with the value for the LBGs in this sample.For the LBGs in A370a that do not have near-IR data, we use the z ′ band photometry to estimate the UV continuum flux.Since we are not able to determine β UV , we adopt the median β UV = −1.86 from the other fields.For simplicity, we assume that the flux blueward of Lyα is completely absorbed.If a Lyα line is in the z ′ band, its flux is subtracted from the z ′ -band photometry.The calculated M UV and Lyα rest-frame EW values are listed in Columns 9 and 10 of Table 2.
Figure 6 shows the distribution of the LBG redshifts and Lyα EWs with respect to M UV in our sample.These LBGs span a UV magnitude range of −22.5 to −19.0 mag, corresponding to 0.3 − 3.4 times the characteristic luminosity of galaxies at z ∼ 6 (Bouwens et al. 2021), and a Lyα luminosity range of (0.4 − 3) × 10 42 erg s −1 .They represent the most luminous galaxies at z ∼ 6 in terms of both UV continuum luminosity and Lyα luminosity.Ten LBGs in the 'good' sample have sufficient multiwavelength imaging data (especially IRAC 3.6 µm and 4.5 µm photometry) that allow us to perform SED modeling.We model their SEDs and derive their stellar populations using the Code Investigating GALaxy Emission (CIGALE; Noll et al. 2009).Our spectroscopic redshifts remove one critical free parameter in the fitting process.We use the stellar population synthesis models of Bruzual & Charlot (2003) and adopt the Salpeter initial mass function.Given the limited number of available photometric data points, we use as few free parameters as possible.We use constant SFRs and fix metallicity to be 0.2 Z ⊙ .A nebular emission template based on Inoue ( 2011) is included, with ionization parameter U , metallicity, and line width fixed.We do not include AGN components during the model fitting because these LBGs do not show any AGN features in the observed bands.Furthermore, Lyα emission and radiative transfer in galaxies are very complex (e.g., Charlot & Fall 1993).Our LBGs are among galaxies with the strongest Lyα emission, and the nebular emission model built in CIGALE can not reproduce the Lyα emission in the SEDs of many LBGs unless we assume an extremely young stellar population (see Section 5.2).Therefore, we mask out the Lyα line in the nebular template of CIGALE when fitting and subtract- ing the Lyα contribution from the broadband photometry.
We mainly constrain stellar mass M * , age, star formation rate (SFR), and dust reddening E(B − V ).We use the dust extinction law of Calzetti et al. (2000) with E(B − V ) in a range of 0 − 0.6.We allow age to vary between 1 and 800 Myr (the age of the Universe at z = 6 is 914 Myr).The distributions of the derived properties are shown in Figure 7 and the SED modeling results are shown in Figure 8.As seen from the figures, the derived values of the properties all span wide ranges, indicating that these LBGs have a variety of stellar populations.The LBGs are strongly star-forming galaxies with SFRs roughly between 30 and 1000 M ⊙ yr −1 .While most of them have SFRs around 100 M ⊙ yr −1 , two LBGs, (No. 5 and No. 13) have SFRs close to 1000 M ⊙ yr −1 .This is not surprising, because No. 5 is very luminous with M UV ≈ −22 mag and No. 13 is also luminous (M UV ≈ −20 mag) with a very high Lyα EW (∼ 190 Å).The two galaxies also have high stellar masses, relatively large ages, and relatively high dust extinction.No. 13 has been observed by JWST NIR-Cam and will be further discussed later.
Observational studies of stellar populations in Lyα galaxies have suggested that young and low-mass galaxies tend to hold Lyα lines with large EWs (e.g., Ono et al. 2010;Hayes et al. 2023;Roy et al. 2023).This can be explained by increasing gas and dust content in more massive galaxies and thus lower Lyα escape fractions, or intrinsically low Lyα production by more evolved populations.The middle panel of Figure 8 shows that our galaxies generally follow this trend that less massive galaxies having higher Lyα EWs

Galaxy Morphology
Two LBGs in our sample (No. 7 and No. 13) have high-resolution JWST NIRCam images, and three LBGs (No. 5,No. 11,and No. 26) have HST WFC3 near-IR images.The JWST images are retrieved from the Mikulski Archive for Space Telescopes (MAST) and then reduced by the JWST Calibration Pipeline (version 1.11.4).We basically follow the procedure of the CEERS team for the image reduction (Finkelstein et al. 2023).The final resampled and combined images have a pixel scale of 0. ′′ 02 for the short wavelengths (SW; F090W, F115W, F150W, and F200W) and 0. ′′ 04 for the long wavelengths (LW; F277W, F356W, F410M, and F444W).The HST data of the three LBGs are from the CANDELS survey and the images are downloaded from the website 1 .The cutout images of the five LBGs are shown in Figure 9.
We measure morphological parameters for the five LBGs using GALFIT.For either of the two LBGs with NIRCam images, we stack its SW images and its LW images to produce two combined images.This is to increase S/N, assuming that the object morphology and image quality are similar in individual input bands.Likewise, we stack the F125W and F160W images for each of the three LBGs with HST data.A PSF is obtained for each combined image by stacking several unsaturated stars in the same image, and the resultant FWHM of PSFs are shown as the red circles in Figure 9.The FWHM values of the PSFs in the NIRCam SW and LW bands are about 0. ′′ 063 and 0. ′′ 148 (≈ 0.36 and 0.84 kpc at z = 6) , respectively.The PSF FWHM of the HST images is about 0. ′′ 223 (≈ 1.3 kpc at z = 6).We then use GALFIT to model the galaxies on the stacked images, and constrain a few basic parameters, including the effective radius R e along the major axis, Sérsic index n (allowed to vary between 0.3 and 4), and projected minor-to-major axis ratio (b/a).We find that all galaxies except No. 7 can be well fitted by one Sérsic profile.No. 7 LBG seems to have two separate components and will be further discussed later.The fitting results of the Sérsic profiles are reported in Table 3.
Figure 9 and Table 3 show that most galaxies in our sample are compact.For example, No. 13 is barely resolved.Its effective radius is R e = 0.25 kpc in the SW bands and R e = 0.38 kpc in the LW bands, making it one of the most compact galaxies in this redshift and luminosity range (see, e.g., Sun et al. 2023).Its large Lyα EW, high SFR, and compact size is consistent with Green Pea (GP) galaxies found at low redshift (Cardamone et al. 2009) that likely have compact star-forming regions embedded in diffuse, older stellar components (Amorín et al. 2012;Clarke et al. 2021).The GP galaxies are thought to have numerous counterparts at the reionization epoch and some of them may have been discovered recently by JWST (Hall 2023;Rhoads et al. 2023).Studies of local GP galaxies have provided valuable information on how Lyα photons escape from ISM and CGM, which help us better understand their highredshift counterparts (e.g., Henry et al. 2015;Yang et al. 2017).Another possible explanation for the compactness of the No. 13 galaxy is that this galaxy holds an AGN that partly powers its strong Lyα emission.
No. 7 LBG has a quite extended structure and appears to show two components, including a major component and a much weaker component to the north-west of the major one.The major component is more than five times brighter than the minor component in all four bands, so the major component dominates the radiation in these bands.It is unclear whether the weak component is part of the galaxy or a foreground object.In addition, mergers and interacting systems occur more frequently at high redshift (Rodriguez-Gomez et al. 2015), so it is possible that No. 7 LBG is a merger.To verify that the two components belong to the same source, we estimate the probability that the distance between any two random sources is less than the distance between the two components of No. 7 LBG (0.36 ′′ ).Based on the NIRCam images used above, we find that the probability is only 0.02, suggesting that the minor component is unlikely a foreground object.Therefore, we use an additional single Sérsic profile to model the minor component of No. 7, and we obtain a much better fitting result.Table 3 shows the result for the major component.

Fraction of Lyα emitters in LBGs
The fraction of LAEs among LBGs is an important tracer of the IGM state during the epoch of reionization, and many works have been done to constrain the fraction of the neutral hydrogen and its redshift evolution.This method is also called the Lyα visibility test.Previous studies have shown that the fraction of LAEs in LBG samples increases steadily from low redshift to z ∼ 6, and drops dramatically towards higher redshift (e.g., Stark et al. 2010Stark et al. , 2011;;Ono et al. 2012;Bian et al. 2015).It is difficult to explain this using a sudden change in the physical properties of galaxies, but it can be easily explained by the increase of the fraction of the neutral hydrogen in the IGM that attenuates Lyα emission.Before we calculate the LAE fraction in our LBG sample, we would like to clarify the We measure the LAE fraction among LBGs at z ∼ 6 using our spectroscopically confirmed LBG sample.As mentioned earlier, our LBG selection criteria are slightly looser than those used in previous studies, which may result in a slightly higher contamination rate.Therefore, our result can be regarded as a lower limit.Because of the presence of strong OH skylines in the wavelength range that we probe, the identification of the Lyα line is incomplete.We run simulations to estimate the completeness for the Lyα line identification in the spectra.We generate a grid of model Lyα lines with a range of redshifts z = 5.3-6.8 and line strengths L(Lyα) = 10 42.5 -10 43.5 erg/s using the Lyα template from Ning et al. (2020).For each spectrum without a Lyα detection, we insert an artificial Lyα line.We then search for these lines using the same method that was used for real spectra in Section 3.1.The derived completeness function is the recovery rate of the lines in each redshift and Lyα luminosity bin.
We then apply the completeness correction to the LAE fraction and the result is shown in Figure 10.The upper panel of the figure shows the cumulative distribution of the Lyα rest-frame EWs in the 'good' LBG sample.The fraction in a fainter subsample (−21.75 < M UV < 20.25 mag) is marginally higher than that in the whole LBG sample (−21.75 < M UV < −20.25 mag).The results are generally consistent with previous studies (Pentericci et al. 2011;Stark et al. 2011;Curtis-Lake et al. 2012;De Barros et al. 2017).
In the lower panel of Figure 10, we show the LAE fractions χ 25 Lyα in both 'good' and the whole samples, and compare them with previous studies (Stark et al. 2011;Curtis-Lake et al. 2012;Ono et al. 2012;Schenker et al. 2014;Tilvi et al. 2014;Cassata et al. 2015;De Barros et al. 2017) in a redshift range of 4 − 8.The UV magnitude range of the LBGs is −21.75 < M UV < −20.25 mag, which corresponds to the fainter part of our sample and the brighter part of the previous samples.As seen from the figure, our result is broadly consistent with previous results, albeit with a large scatter.The large scatter and error bars reflect complex sample bias and incompleteness.The figure also shows that χ 25 Lyα increases from lower redshifts to z ∼ 6 and then decreases towards higher redshifts, consistent with the previous claim.The mild decline of χ 25 Lyα may suggest that the evolution of the IGM opacity from z ∼ 7 to z ∼ 6 is not very dramatic, if the Lyα visibility mostly depends on the IGM state.Nevertheless, larger and more complete samples are needed.

SED modeling with strong Lyα emission
In Section 4.2 when we modeled the SEDs of the LBGs, we did not consider the Lyα emission due to its complexity.Here we discuss how the results change if we consider the Lyα emission in our SED modeling.As has been pointed out, the SED modeling for z ∼ 6 galaxies has a strong degeneracy between young galaxies with prominent nebular emission and old galaxies with strong Balmer breaks (e.g., Jiang et al. 2015).To include the Lyα emission, we make a pseudo narrowband filter with a spectral resolving power of R = 100 at the wavelength of Lyα for each LBG.The narrowband photometry is then calculated from the Lyα flux and implemented for SED modeling.In the modeling procedure, the ionization parameter U , gas metallicity in the nebular emission model, and the ratio of continuum to line reddening are considered as free parameters.
The derived properties of the stellar populations and the comparison with the previous results from Section 4.2 are shown in Figure 11.The best-fit χ 2 r values of the two methods are about the same.Figure 11 clearly shows that when the Lyα emission is considered in the SED modeling procedure, these galaxies are much younger with significantly lower stellar masses.The reason is straightforward.In order to produce strong Lyα emission seen in the galaxies, the stellar populations must be very young with high SFRs, which further results in lower stellar masses.It should be noted, however, that it is difficult to decide how to include the Lyα emission for SED modeling due to the complexity of the Lyα emission.For example, the measurement of the Lyα flux is almost always a lower limit, this is mainly because the resonant scattering of the Lyα photons often forms a large diffuse halo, but the actual measurement is usually done on the central part.It is also because the Lyα escape fraction varies a lot among galaxies.In addition, there are other sources of uncertainties that are particularly important for z ≥ 6 galaxies, including the variation of the IGM absorption along the line-of-sights.

Search for AGN activity
The LBGs in our sample have strong Lyα and UV continuum emission.In addition, some of them are morphologically very compact.It is possible that they harbor a central AGN that partly contributes to the total emission.We search for possible evidence of AGN activity in these galaxies.First of all, they do not have typical Type I AGN components, as seen from their narrow Lyα emission lines.Figure 5 shows that the stacked Lyα line displays a narrow profile with a velocity dispersion of ∆v ∼ 250 km/s (IGM absorption is neglected here), which is much smaller than those in typical broadline AGN (≥ 2000 km/s).
The wavelength coverage of our spectra is from 7600 Å to 9600 Å, corresponding to a rest-frame range of 1090-1370 Å for z ∼ 6 objects.This range covers the N V 1240 Å emission line that is often regarded as evidence of the presence of AGN.We do not detect any N V emission in the combined spectrum of the LBGs.For a typical AGN, however, the N V flux is only a few percent of the Lyα flux (e.g., Vanden Berk et al. 2001), so even if any of our LBGs contain an AGN component, its N V flux is below our detection limit.
We also check these galaxies in deep Chandra X-ray images, including the Chandra Deep Field-South field (CDF-S) (Luo et al. 2017), the UKIDSS Ultra Deep Survey Field (Kocevski et al. 2018), and the Chandra Cosmos Legacy Survey field (Civano et al. 2016).None of them is individually detected.X-ray detections were actually excluded in the first steps of the target selection.The No.5 and No.6 LBGs are in the 7 Ms CDF-S field, and their stacked image have an effective exposure time of over 12 Ms.For the stacked 0.5-7 keV image, we measure the X-ray photon counts in a 1. ′′ 5 radius aperture, which encloses 90% of the total energy at 6.4 keV.The background counts are computed in an annular region with the inner and outer radii of 4 ′′ and 8 ′′ .The 1σ error of the net counts is calculated based on the Poisson errors on the extracted source and background counts.The calculated 1σ upper limit is 6.3 counts.To estimate the upper limit of a possible X-ray detection, we assume a typical type 1 quasar with a similar UV magnitude at the same redshift as the No.5 and No.6 LBGs (z ∼ 5.7, M UV = −22 mag).We also assume the Shen et al. (2020) SED template for the quasar.We adopt the Galactic neutral hydrogen column density along the line-of-sight to the CDF-S n H = 8.8 × 10 19 cm −2 (Stark et al. 1992).No other absorption is applied.In this case, the 12 Ms exposure would yield a 2σ detection (19.2 counts in 0.5-7 keV band) for this quasar.This simple calculation means that AGN contributes less than onethird (1σ upper limit) of the total luminosity in the two LBGs.Therefore, we are not able to rule out the possibility of weaker AGN activities in these galaxies.The stacked images in another two fields are much shallower and thus cannot provide useful constraints.

SUMMARY
We have presented a sample of 45 spectroscopically confirmed LBGs at z ∼ 6 in four well-studied fields, including SXDS, A370, ECDFS, and COSMOS.This sample is one of the largest samples of spectroscopically confirmed LBGs at this redshift.The LBG candidates were selected as i-band dropout objects from deep broadband images.The spectroscopic observations were carried out using M2FS on the Magellan Clay telescope.We identified the Lyα lines of the galaxies based on their 1D and 2D M2FS spectra.Their redshifts and Lyα luminosities were measured by fitting a composite Lyα line template to the individual 1D spectra.The absolute UV magnitude and UV slope were calculated by fitting a powerlaw to the available multi-wavelength photometric data.The UV absolute magnitudes of these LBGs are between -22.5 and -19.0 mag, corresponding to 0.3 − 3.4 times the characteristic luminosity of galaxies at z ∼ 6, with a Lyα EW range from ∼10 Å to ∼400 Å.They represent the most UV luminous galaxies known at z ∼ 6.
The SED modeling of ten LBGs in the sample reveals that they have a wide range of stellar masses and ages with high SFRs.We note that whether to consider the Lyα emission in the SED modeling procedure has a significant impact on the measurement of the stellar populations in these galaxies, indicating the complexity of the Lyα emission.Five LBGs have high-resolution JWST NIRCam images or HST WFC3 near-IR images.Four of them show compact morphology, and one appears to have two components that might be a merger.We calculated the LAE fraction among the photometrically selected LBGs after correcting for the Lyα-identification completeness.The fraction of ∼ 0.2 is consistent with previous works, and it supports a moderate evolution of the IGM opacity at the end of cosmic reionization when previous results at lower and higher redshifts were combined.Finally, we did not find evidence of strong AGN activities in these galaxies.

Figure 1 .
Figure 1.Deep fields used in this work.The large circles represent the M2FS pointings.All points represent the LBG targets observed by our M2FS survey.The color-coded points represent the spectroscopically confirmed LBGs with Lyα emission: the large points represent LBGs in the 'good' sample and the small points represent LBGs in the 'possible' sample.The colors indicate their spectroscopic redshifts measured from the Lyα lines.The histogram on the right shows the the redshift distribution of the 'good' sample (blue) and the 'possible' sample (orange).

Figure 2 .
Figure 2. Upper panel: transmission curves of the Suprime-Cam broadband filters that were used for our LBG candidate selection.Lower panel: expected values of i ′ − z ′ or I − z ′ as a function of redshift.The simulated LBG spectra have a UV continuum slope of -2.0 and a Lyα emission line with two EW values, and the IGM absorption model from Inoue et al. (2014) is applied.

Figure 3 .
Figure 3. M2FS 1D and 2D spectra of the 34 LBGs in our z ≈ 6.6 sample.In each panel, the 1D spectrum has been binned to 3 Å per pixel in arbitrary units.The light grey line indicates the 1σ uncertainty.The best-fit Lyα line template is over-plotted in blue.The red dashed line shows the peak position of the Lyα emission line.The source numbers correspond to those shown in Column 1 in Table2.

Figure 4 .
Figure 4. Relation between the Lyα luminosities L(Lyα) calculated from the spectral fitting method and from the photometric data (see details in Section 3.2).The blue points represent the z ≈ 5.7 LAE sample in Ning et al. (2020) and the red points represent our LBG sample.The bestfit relation from the z ≈ 5.7 LAE sample is depicted by the black dashed line, and the shaded regions indicate its 1σ and 3σ confidence regions, respectively.We apply this relation to estimate the Lyα luminosities of the LBGs in this work.

Figure 5 .
Figure 5. Stacked Lyα line spectra for the LBGs in two Lyα luminosity bins (upper panel) and two redshift bins (lower panel).The spectra have been normalized so that the peak values are 1.The dashed lines represent the instrument resolution assuming a Gaussian profile.

Figure 6 .
Figure 6.Distributions of the LBG redshifts (upper panel) and Lyα EWs (lower panel) with respect to MUV.The 'good' LBG sample and 'possible' LBG sample are both presented.The median errors are shown in the lower right corner, and the error bar of redshift is too small to be seen.

Figure 7 .
Figure 7. Properties of the ten LBGs measured from the SED modeling.They all cover wide ranges.
Figure 9. JWST NIRCam and HST ACS and WFC3 cutout images of five LBGs.The image size is 2 ′′ × 2 ′′ .The red circles in each stamp image indicate the FWHM of the PSF.

Figure 10 .
Figure 10.Upper panel: cumulative distribution of the Lyα rest-frame EWs P(> EW) for the z ∼ 6 LBGs.The red symbols represent LBGs with −21.75 < MUV < −20.25 and the blue circles represent all LBGs within −21.75 < MUV < −20.25.The results are consistent with previous measurements.The lines are the best fits to the distributions.Lower panel: fraction of LAEs with EW > 25 Å(χ 25 Lyα ) for −21.75 < MUV < −20.25 galaxies at 4 ≤ z ≤ 8.The red filled circles are calculated from the 'good' LBG sample and the red open circles are calculated from both 'good' and 'possible' samples.The results are compared with previous studies.
Figure 11.Comparison of the stellar populations obtained by two SED modeling methods, with and without Lyα constraints (see details in Section 5.2).

Table 1 .
Survey Fields Numbers of spectroscopically confirmed LBGs in the 'good' sample and the 'possible' sample.

Table 3 .
Galaxy morphological parametersThe parameters of Sérsic profile of the main component is reported.definitions of LAEs and LBGs since the definitions are slightly different in the literature.LBGs mean galaxies selected by the Lyman-break technique, regardless of the existence of Lyα emission.For LBGs at z ≥ 6, they are usually very faint, and spectroscopic identifications primarily rely on their Lyα emission.This means that spectroscopically confirmed LBGs typically have strong Lyα emission, and also means that the vast majority of the known LBGs at z ≥ 6 are not spectroscopically confirmed (or just candidates).LAEs are galaxies with strong Lyα emission (typically EW ≥ 25 or 50 Å).Narrowband-selected galaxies are mostly LAEs by definition.If an LBG has strong Lyα emission, it is also an LAE, meaning that most spectroscopically confirmed z ≥ 6 LBGs to date are also LAEs.In the following analyses, we use EW ≥ 25 Å to define LAEs.