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

We carried out a large spectroscopic survey of galaxies at 5.5 < z < 6.8, using the large FoV, fiber-fed, multiobject 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^{\prime}$-diameter FoV, M2FS is very efficient in identifying relatively bright high-redshift galaxies (e.g., Oyarzún et al. 2016, 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². These fields have deep Subaru Suprime-Cam imaging data in a series of broad bands ${BVR}(r^{\prime} )I(i^{\prime} )z^{\prime}$ 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.

Table 1. Survey Fields

Fields	Coordinates	Area	Filters	$R/r^{\prime}$	$I/i^{\prime}$	$z^{\prime}$	Candidates	Good	Possible
	(J2000)	(deg2)		(mag)	(mag)	(mag)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)
SXDS	02:18:00–05:00:00	1.12	$R,i^{\prime} ,z^{\prime}$	27.4	27.4	26.2	368	25	7
A370a	02:39:55–01:35:24	0.16	$R,I,z^{\prime}$	27.0	26.2	26.3	86	3	2
ECDFS	03:32:25–27:48:18	0.22	$r^{\prime} ,i^{\prime} ,z^{\prime}$	27.4	27.5	26.7	105	3	2
COSMOS	10:00:29 + 02:12:21	1.26	$r^{\prime} ,i^{\prime} ,z^{\prime}$	26.7	26.3	25.5	145	3	0

Note. Column (1): field names. Column (2): center coordinates. Column (3): area coverage. Column (4): broadband filters. Columns (5)–(7): magnitude limits (5σ detections in a 2'' diameter aperture). Column (8): numbers of LBG candidates. Columns (9)–(10): numbers of spectroscopically confirmed LBGs in the "good" sample and the "possible" sample.

Download table as:  ASCIITypeset image

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 are presented in Section 2.3.

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.36.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 wave front 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 minutes, 45 minutes, or 1 hr, depending on airmass and weather conditions.

Zoom In Zoom Out Reset image size

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.

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^{\prime}$ color (here, i means either $i^{\prime}$ 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^{\prime}$ band. We applied the following color cuts:

$$\begin{eqnarray*}&&i-z^{\prime} \gt 1.3\ \mathrm{or}\ i\lt 3\sigma \ \mathrm{detection},\end{eqnarray*}$$

$$\begin{eqnarray*}&&r-z^{\prime} \gt 2.0\ \mathrm{or}\ r\lt 3\sigma \ \mathrm{detection}.\end{eqnarray*}$$

Zoom In Zoom Out Reset image size

To illustrate the evolution of the $i-z^{\prime}$ 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^{\prime}$ 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.

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" sample 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 for 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.

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 interstellar medium (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 high-redshift 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. (2020, 2022), COSMOS1 and COSMOS3 have some alignment problems during the observations, so they only have a few LBGs confirmed.

Zoom In Zoom Out Reset image size

Table 2. Confirmed LBGs with Lyα Emission Lines

No.	Field	R.A.	Decl.	Redshift	i	$z^{\prime}$	log10 L(Lyα)	MUV	Lyα EW	βUV
		(J2000)	(J2000)		(mag)	(mag)	(erg s−1)	(mag)	(Å)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)
Good Sample
1	A370a	02:40:21.24	−01:46:37.52	6.043	>27.2	25.43	42.62	−21.31	14.77	⋯
2	A370a	02:40:18.51	−01:43:51.66	6.042	>27.2	24.84	42.73	−21.93	10.66	⋯
3	A370a	02:40:21.60	−01:33:42.16	6.319	>27.2	25.24	43.10	−21.83	27.14	⋯
4	ECDFS	03:31:56.08	−27:59:06.50	5.766	>28.5	25.61	42.58	−20.86	39.40	1.38
5	ECDFS	03:32:18.92	−27:53:02.81	5.561	>28.5	24.48	43.37	−21.96	44.83	−1.84
6	ECDFS	03:32:30.85	−27:38:22.20	5.808	27.84	24.94	43.26	−22.25	44.38	0.55
7	COSMOS	09:59:41.18	+01:55:51.75	6.083	>27.3	24.78	42.65	−21.22	14.87	−2.44
8	COSMOS	09:59:56.54	+02:12:27.15	5.992	27.12	24.76	43.27	−21.62	42.13	−2.52
9	COSMOS	10:00:12.59	+02:28:52.08	6.146	>27.3	25.03	43.42	−21.36	87.75	⋯
10	SXDS1	02:18:41.01	−05:12:47.18	6.474	27.9	25.72	42.93	−21.76	19.85	⋯
11	SXDS1	02:18:00.91	−05:11:37.77	6.029	27.67	24.75	43.35	−21.48	52.53	−2.98
12	SXDS1	02:18:38.92	−05:09:44.02	6.192	>28.4	25.06	43.17	−21.55	41.04	−1.93
13	SXDS1	02:17:43.25	−05:06:47.60	6.087	>28.4	25.90	43.25	−20.08	189.33	−1.98
14	SXDS1	02:17:48.40	−05:04:10.79	6.108	>28.4	25.67	43.25	−19.80	229.40	−2.22
15	SXDS1	02:18:23.34	−04:58:58.18	6.071	>28.4	25.70	43.20	−20.95	53.24	−3.64
16	SXDS1	02:18:07.14	−04:58:41.50	6.050	>28.4	24.84	43.10	−21.61	22.74	−3.64
17	SXDS1	02:18:30.93	−04:50:58.73	5.850	27.04	24.65	42.67	−21.27	13.86	−2.77
18	SXDS2	02:17:28.71	−04:45:19.61	6.104	>28.4	25.67	43.45	−19.48	527.92	⋯
19	SXDS2	02:18:44.06	−04:35:54.24	6.078	>28.4	25.39	43.26	−20.98	85.69	⋯
20	SXDS2	02:18:34.48	−04:34:25.86	5.745	26.59	25.19	42.92	−22.21	17.81	−0.21
21	SXDS2	02:17:25.27	−04:34:12.85	6.003	>28.4	25.69	42.98	−20.87	49.57	⋯
22	SXDS3	02:18:24.72	−05:35:41.82	6.059	27.55	25.67	42.84	−20.97	33.18	⋯
23	SXDS3	02:18:23.90	−05:34:49.34	6.192	28.19	24.95	42.94	−22.06	24.16	0.32
24	SXDS3	02:18:05.68	−05:31:21.90	6.064	>28.4	25.55	43.29	−20.65	125.98	⋯
25	SXDS3	02:18:34.47	−05:20:20.96	6.048	>28.4	25.41	43.05	−21.16	44.85	⋯
26	SXDS3	02:18:06.96	−05:16:21.93	5.928	27.07	25.26	43.04	−20.95	55.84	−1.64
27	SXDS4	02:20:20.49	−05:01:45.06	6.046	>28.4	25.56	43.02	−20.99	49.65	⋯
28	SXDS4	02:19:26.97	−04:50:17.71	5.848	26.67	25.27	42.67	−21.40	15.23	⋯
29	SXDS5	02:15:24.62	−05:06:07.27	5.753	26.34	25.07	42.86	−21.57	19.90	⋯
30	SXDS5	02:16:16.53	−05:02:17.80	6.051	>28.4	25.51	43.15	−20.88	73.98	−1.82
31	SXDS5	02:15:50.69	−05:00:57.88	6.076	27.32	25.50	43.06	−21.62	48.69	0.42
32	SXDS5	02:17:02.70	−04:56:59.28	5.793	26.05	24.28	42.85	−21.96	13.90	−1.73
33	SXDS5	02:16:27.81	−04:55:34.25	6.016	27.36	24.35	43.37	−21.13	77.92	−2.89
34	SXDS5	02:16:02.36	−04:49:08.81	6.043	>28.4	25.60	43.03	−20.94	52.51	⋯
Possible Sample
1	A370a	02:39:26.57	−01:36:03.84	6.463	>27.2	25.59	42.43	−22.05	4.76	⋯
2	A370a	02:39:28.99	−01:33:31.00	6.383	26.73	25.84	42.68	−21.51	13.99	⋯
3	ECDFS	03:32:50.43	−27:50:40.23	5.536	28.45	25.63	42.89	−21.16	55.39	0.85
4	ECDFS	03:31:41.99	−27:46:30.91	5.795	>28.5	25.56	42.63	−21.09	18.11	⋯
5	SXDS2	02:18:14.62	−04:46:07.52	5.778	27.5	26.03	42.47	−20.62	19.60	⋯
6	SXDS2	02:17:49.61	−04:22:41.36	6.074	>28.4	25.74	43.10	−20.65	81.85	⋯
7	SXDS3	02:17:13.11	−05:37:53.43	5.768	26.32	24.76	42.64	−21.88	9.04	⋯
8	SXDS4	02:20:22.93	−05:04:31.83	5.849	25.81	24.23	43.08	−22.74	18.18	0.40
9	SXDS4	02:19:29.35	−04:53:38.27	5.862	27.83	25.38	42.47	−21.69	13.55	1.11
10	SXDS4	02:19:37.99	−04:52:41.57	6.246	28.14	25.65	43.10	−20.98	52.87	−2.41
11	SXDS5	02:15:48.66	−05:11:58.28	5.938	27.32	25.02	43.45	−21.53	79.60	⋯

Notes. Column (1): object numbers. Column (2): field names. Columns (3)–(4): source coordinates. Column (5): spectroscopic redshifts (their errors are smaller than 0.001). Columns (6)–(7): Subaru Suprime-Cam i-band and $z^{\prime}$-band magnitudes. The upper limits indicate 2σ detections. Column (8): Lyα line luminosities. Column (9): absolute UV magnitudes at rest-frame 1500 Å. Column (10): Lyα line EWs. Column (11): UV slopes calculated from multiband data.

Download table as:  ASCIITypeset image

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 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 χ² 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. The 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, 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 fields and were observed at the same time, the LBGs in our sample should follow the same relation. We apply this relation to measure the line flux of our LBGs, and the measured Lyα luminosities L(Lyα) are listed in column (8) of Table 2. The IGM absorption blueward of Lyα is not considered.

Zoom In Zoom Out Reset image size

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 an 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. (2020, 2022) found that more luminous galaxies tend to have higher Lyα FWHMs, because more massive host halos possess higher H i 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.

Zoom In Zoom Out Reset image size

In this section, we collect multiwavelength 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 (HSC) Subaru Strategic Program (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), the VISTA Deep Extragalactic Observations Survey (VIDEO; Jarvis et al. 2013), the Cosmic Assembly Near-IR Deep Extragalactic Legacy Survey (CANDELS; Grogin et al. 2011; Koekemoer et al. 2011), and the Cosmic Dawn Survey by Spitzer Space Telescope (only the 3.6, 4.5 μm data are used, Euclid Collaboration et al. 2022). The No. 13 LBG is covered by the JWST Cycle-1 program (GO 1837), Public Release IMaging for Extragalactic Research (PRIMER; Dunlop et al. 2021), and has NIRCam imaging in eight bands (F090W, F115W, F150W, F200W, F277W, F356W, F410M, and F444W). Ning et al. (2023) carried out a pilot study of seven spectroscopically confirmed z ∼ 6 LBGs, including No. 13 LBG, using the JWST NIRCam imaging data. The No. 7 LBG is covered by another JWST Cycle-1 program COSMOS-Web (GO 1727, Casey et al. 2023), and has NIRCam imaging in four bands (F115W, F150W, F277W, and F444W).

Most LBGs in our sample except those in the A370a field have multiwavelength data that can be used to determine UV continuum slope β_UV and other physical properties via spectral energy distribution (SED) modeling. The IRAC images of many LBGs are blended with nearby sources. In this case, GALFIT (Peng et al. 2002, 2010) is used to model and subtract contaminant sources. Photometry of the target is then done after the deblending.

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 ${\beta }_{\mathrm{UV}}=-{1.73}_{-0.20}^{+0.14}$ 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^{\prime}$-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^{\prime}$ band, its flux is subtracted from the $z^{\prime}$-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⁴² erg s⁻¹. They represent the most luminous galaxies at z ∼ 6 in terms of both UV continuum luminosity and Lyα luminosity.

Zoom In Zoom Out Reset image size

Ten LBGs in the "good" sample have sufficient multiwavelength imaging data (especially IRAC 3.6 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 the 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 cannot 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 subtracting the Lyα contribution from the broadband photometry.

We mainly constrain stellar mass M_*, age, 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 the 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⁻¹. While most of them have SFRs around 100 M_⊙ yr⁻¹, two LBGs, (No. 5 and No. 13) have SFRs close to 1000 M_⊙ yr⁻¹. 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 NIRCam and will be further discussed later.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 the increasing gas and dust content in more massive galaxies and thus the 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 have higher Lyα EWs

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 002 for the short wavelengths (SW; F090W, F115W, F150W, and F200W) and 004 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. ⁴ The cutout images of the five LBGs are shown in Figure 9.

Zoom In Zoom Out Reset image size

We measure the 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 point-spread function (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 0063 and 0148 (≈0.36 and 0.84 kpc at z = 6), respectively. The PSF FWHM of the HST images is about 0223 (≈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, ${\rm{S}}\acute{{\rm{e}}}\mathrm{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 ${\rm{S}}\acute{{\rm{e}}}\mathrm{rsic}$ profile. No. 7 LBG seems to have two separate components and will be further discussed later. The fitting results of the ${\rm{S}}\acute{{\rm{e}}}\mathrm{rsic}$ profiles are reported in Table 3.

Table 3. Galaxy Morphological Parameters

No.	Bands	n	Re (kpc)	b/a
13	F090W, F115W, F150W, F200W	1.76 ± 0.45	0.25 ± 0.03	0.67 ± 0.08
	F277W, F356W, F410M, F444W	1.10 ± 0.27	0.38 ± 0.03	0.99 ± 0.12
7 a	F115W, F150W	0.49 ± 0.13	0.59 ± 0.04	0.51 ± 0.04
	F277W, F444W	0.67 ± 0.11	1.02 ± 0.05	0.61 ± 0.04
5	F125W, F160W	1.18 ± 0.19	0.94 ± 0.09	0.92 ± 0.08
11	F125W, F160W	1.00 ± 0.29	1.16 ± 0.11	0.59 ± 0.09
26	F125W, F160W	0.44 ± 0.19	1.42 ± 0.20	0.70 ± 0.11

Note.

^a The parameters of ${\rm{S}}\acute{{\rm{e}}}\mathrm{rsic}$ profile of the main component are reported.

Download table as:  ASCIITypeset image

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. 2024). Its large Lyα EW, high SFR, and compact size are 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 circumgalactic medium (CGM), which helps us better understand their high-redshift 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 northwest of the major one. The major component is more than 5 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 (036). 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 ${\rm{S}}\acute{{\rm{e}}}\mathrm{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.

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 toward higher redshift (e.g., Stark et al. 2010, 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 definitions of LAEs and LBGs since the definitions are slightly different in the literature. LBGs mean the 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.

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 P(>EW) 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).

Zoom In Zoom Out Reset image size

In the lower panel of Figure 10, we show the LAE fractions ${\chi }_{\mathrm{Ly}\alpha }^{25}$ 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 a complex sample bias and incompleteness. The figure also shows that ${\chi }_{\mathrm{Ly}\alpha }^{25}$ increases from lower redshifts to z ∼ 6 and then decreases toward higher redshifts, consistent with the previous claim. The mild decline of ${\chi }_{\mathrm{Ly}\alpha }^{25}$ 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.

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 ${\chi }_{r}^{2}$ 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.

Zoom In Zoom Out Reset image size

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 broad-line 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 images 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 15 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¹⁹cm⁻² (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 one-third (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.