Luminous Lyman-alpha Emitters with Very Blue UV-continuum Slopes at Redshift 5.7<= z<= 6.6

We study six luminous Lyman-alpha emitters (LAEs) with very blue rest-frame UV continua at $5.7\le z \le 6.6$. These LAEs have previous HST and Spitzer IRAC observations. Combining our newly acquired HST images, we find that their UV-continuum slopes $\beta$ are in a range of $-3.4\le \beta \le -2.6$. Unlike previous, tentative detections of $\beta \simeq -3$ in photometrically selected, low-luminosity galaxies, our LAEs are spectroscopically confirmed and luminous ($M_{\rm UV}<-20$ mag). We model their broadband spectral energy distributions (SEDs), and find that two $\beta\simeq-2.6\pm0.2$ galaxies can be well fitted with young and dust-free stellar populations. However, it becomes increasingly difficult to fit bluer galaxies. We explore further interpretations by including non-zero LyC escape fraction $f_{\rm esc}$, very low metallicities, and/or AGN contributions. Assuming $f_{\rm esc}\simeq0.2$, we achieve the bluest slopes $\beta\simeq-2.7$ when nebular emission is considered. This can nearly explain the SEDs of two galaxies with $\beta\simeq-2.8$ and --2.9 ($\sigma_{\beta}=0.15$). Larger $f_{\rm esc}$ values and very low metallicities are not favored by the strong nebular line emission (evidenced by the IRAC flux) or the observed (IRAC 1 - IRAC 2) color. Finally, we find that the $\beta\simeq-2.9$ galaxy can potentially be well explained by the combination of a very young population with a high $f_{\rm esc}$ ($\ge0.5$) and an old, dusty population. We are not able to produce two $\beta \simeq -3.4 \pm0.4$ galaxies. Future deep spectroscopic observations are needed to fully understand these galaxies.


INTRODUCTION
In recent years, a number of high-redshift galaxies have been discovered, thanks to the advances of instrumentation on the Hubble Space Telescope (HST) and large ground-based telescopes. These galaxies have played an important role in studies of galaxy formation and evolution at early epochs. Most known galaxies at z ≥ 6 are photometrically selected Lyman-break galaxies (LBGs) or candidates using the dropout technique (e.g., McLure et al. 2011;Yan et al. 2012;Ellis et al. 2013;Bouwens et al. 2015;Ono et al. 2018). Only a small fraction of them have been spectroscopically confirmed (e.g., Finkelstein et al. 2013;Oesch et al. 2015;Jiang et al. 2017;Laporte et al. 2017;Stark et al. 2017). The narrowband (or Lyα) technique provides a complementary method to find high-redshift galaxies, and narrowband observations were mostly done by large-area ground-based observations (e.g., Kashikawa et al. 2011;Ota et al. 2017;Zheng et al. 2017;Ouchi et al. 2018;Hu et al. 2019). While ground-based observations provide samples of luminous galaxies, the majority of faint galaxies have to come from deep HST observations.
The physical properties of high-redshift galaxies are also being investigated. The combination of the HST and Spitzer (IRAC 1 and IRAC 2) observations covers the rest-frame UV and optical ranges for these galaxies, and is thus critical to characterize these objects and measure their stellar populations. The UV continuum slope β (f λ ∝ λ β ) is of particular interest. It provides key information to constrain young stellar populations, yet it does not require Spitzer observations. Previous HST observations have tentatively found extremely steep slopes β ≃ −3 in photometrically selected, very faint galaxies at z ≥ 6 (e.g., Bouwens et al. 2010;Labbé et al. 2010). Such slopes have not been reported in lower-redshift galaxies, nor have they been predicted by cosmological simulations (e.g., Finlator et al. 2011;Ceverino et al. 2019). On the other hand, this finding is controversial, as others have demonstrated the likelihood that the extreme slopes could be caused by contamination and bias due to the nature of photometrically selected samples (e.g., McLure et al. 2011;Dunlop et al. 2012Dunlop et al. , 2013. As the debate goes on, the key is to construct a reliable sample of galaxies with β ≃ −3. We have carried out a detailed study of a sample of 67 spectroscopically confirmed galaxies at z ≃ 6 (Jiang et al. 2013a(Jiang et al. ,b, 2016, based on our Subaru optical, HST near-IR, and Spitzer mid-IR images. These galaxies have blue UV continua with a median slope β ≃ −2.3 at M UV < −19.5 mag. In particular, a non-negligible fraction of them exhibit extreme slopes β ≃ −3 (see Figure 3 in Jiang et al. 2013a). Due to their brightness and secure redshifts, the β measurements were not subject to the contamination or bias mentioned above. The uncertainties σ β are around 0.3−0.5, mainly due to the short UV-continuum baseline or large photometric errors. In order to reduce σ β , we have obtained new HST WFC3 data from HST Cycle 25 for seven relatively bright and promising candidates of extremely blue galaxies selected from the above sample.
In this paper, we present the new HST observations and confirm very steep slopes in six of the seven galaxies. The layout of the paper is as follows. In Section 2, we introduce our HST observations and data reduction. In Section 3, we present our main results of the β measurements. In Section 4, we model the spectral energy distributions (SEDs) these galaxies and discuss our results. We summarize our paper in Section 5. Throughout the paper, all magnitudes are expressed on the AB system. The HST F105W, F110W, F125W, F140W, and F160W bands are denoted as Y 105 , J 110 , J 125 , H 140 , and H 160 , respectively. We use a Λ-dominated flat cosmology with H 0 = 68 km s −1 Mpc −1 , Ω m = 0.3, and Ω Λ = 0.7.

OBSERVATIONS AND DATA REDUCTION
We first briefly review our previous HST observations of 67 galaxies at 5.7 ≤ z ≤ 6.6 (Jiang et al. 2013a). These galaxies were spectroscopically confirmed via Lyα emission lines. They represent the most luminous z ≃ 6 galaxies in terms of Lyα emission or UV-continuum emission. Most of them are located in the Subaru Deep Field (SDF; Kashikawa et al. 2004), and the remaining are located in the Subaru XMM-Newton Deep Survey (SXDS; Furusawa et al. 2008) field. They have deep optical images from the Subaru telescope. They also have deep near-IR and mid-IR images from our previous HST and Spitzer programs. The HST observations of the SDF galaxies were done with a mix of strategies. Most galaxies were observed in the WFC3 J 125 and H 160 bands, with a depth of two orbits per band. For the galaxies in the SXDS field, we have used the CANDELS (Grogin et al. 2011;Koekemoer et al. 2011) J 125 -and H 160 -band images. Based on these data, we previously found that our galaxies have blue UV slopes with a median value of β ≃ −2.3. In addition, a small fraction of them show extreme slopes of β ≃ −3, but the errors were relatively large (σ β = 0.3 ∼ 0.5).
In order to study these extreme galaxies, we have carried out new HST WFC3 observations of seven LAEs (Table 1) in HST Cycle 25 program 15137 (PI: L. Jiang). These were selected from the above sample to be relatively bright (> 10σ detection in J 125 ) with β ≤ −2.8. We took advantage of the existing data and improved the β measurements by extending the UV continuum baseline, increasing the spectral sampling, and/or improving the photometry in H 160 . We observed all objects except ID28 in Y 105 , which provides a longer UV continuum baseline. The depth was one HST orbit per object. Because of the very blue slopes, the galaxies were expected to be brighter in Y 105 , but much fainter in H 160 . Therefore, we obtained new H 160 -band images for four objects (ID07, ID28, ID30, and ID61), and the depth was also one orbit per object. ID28 already had images in J 110 and H 160 , so we obtained one-orbit H 140band images to improve the spectral sampling.
We used the software package DrizzlePac and followed the standard procedure to reduce our HST images. For each galaxy, all images were drizzled onto the same WCS system. The new Y 105 -band images of the two SXDS galaxies were drizzled to match the WCS in the CAN-DELS images. The 'pixfrac' parameter in DrizzlePac was set to be 0.8 and the output plate scale is 0. ′′ 06 pixel −1 . We combined the old and new images in the individual bands. The final depth is roughly one orbit in Y 105 , two orbits in J 125 , and three orbits in H 160 , corresponding to the 10σ detection limits (within a 0. ′′ 6 or 10-pixel circular aperture in diameter) of ∼ 26.3, 26.7, and 26.6 mag, respectively.
In addition to the HST images, we have used deep z ′band imaging data from Subaru Suprime-Cam to extend the UV-continuum baseline when calculating β. These data have been used to select our galaxies, and the details were described in Jiang et al. (2013a). We have slightly improved the z ′ -band photometry of the galaxies. Furthermore, we have also used deep ground-based K s -band data to further extend the UV baseline for two SXDS galaxies (see Section 3). Note-The ID numbers in Column 1 correspond to the numbers in Table 1 in Jiang et al. (2013a). For ID28, the photometry in Column 5 is for J110, and the photometry in Column 6 is for H140. The Ks-band photometry in Column 8 is from Galametz et al. (2013). The IRAC photometry in Columns 9 and 10 is from Jiang et al. (2016).

RESULTS
In this section, we will measure the photometry of the galaxies in multiple bands, and then calculate their UV-continuum slopes using the multi-band photometric data. We will show that our slopes are not estimated from a single UV color as many previous studies did for high-redshift galaxies. Instead, they are measured from 3-4 bands (Figures 1 and 2).

Photometry
We run SExtractor (Bertin & Arnouts 1996) to measure the photometry. The procedure is similar to those in the literature. Usually HST images in a longer wavelength band have a larger size (or FWHM) of the point-spread function (PSF). But this is not the case in our images. We performed sub-pixel dithering during the HST observations, which can improve the PSF. As mentioned earlier, the typical numbers of the HST orbits for our galaxies are 1, 2, and 3 in Y 105 , J 125 , and H 160 , respectively. This results in nearly consistent PSF sizes (within the errors) in the different bands (see e.g., Windhorst et al. 2011). For example, the PSF FWHM values in the three bands, measured from bright stars by SExtractor, are 0. ′′ 22 ± 0. ′′ 01, 0. ′′ 21 ± 0. ′′ 02, and 0. ′′ 22±0. ′′ 02, respectively. Therefore, we skip PSF matching here. To do photometry for a given object, we first stack (inverse-variance weighted average) all its HST images to make a combined image. We then run SExtractor in the dual image mode, using the combined image as the 'detection' image. We calculate the total magnitude within a MAG AUTO elliptical aperture with PHOT AUTOPARAMS values of 1.8 and 2.5. The aperture corrections are computed using standard (larger) PHOT AUTOPARAMS values of 2.5 and 3.5. The results are shown in Table 1.
The photometric uncertainties are also measured by SExtractor. We have converted the weight or variance maps (obtained from DrizzlePac) to rms images following the procedure of Dickinson et al. (2004) that includes correlated noise. The rms images are then fed to SExtractor. The average background noise is estimated using a 2. ′′ 4-wide annulus (larger than the default value by 67%), and the total errors are calculated within the adopted aperture. We evaluate our error measurements using mock data. For a given galaxy in a given image, we generate 100 objects with the same brightness and shape as the galaxy, and put them randomly on the image. We then do photometry using SExtractor as we did for real objects, and compute the mean and standard deviation of the measurements. The new results are consistent with the results above.
The aperture size that we have adopted is roughly 0. ′′ 6 in diameter, larger than those in the majority of the previous studies of z ≥ 6 galaxies. This results in relatively larger measurement errors. On the other hand, this conservative approach reduces the effect from any possible PSF mismatches between different bands. Our galaxies are isolated in the HST images, and thus no obvious neighbors are included in the photometry. The two galaxies ID63 and ID64 are bright enough to test if the magnitudes that we have measured are close to 'real' total magnitudes. We measure their photometry again within a very large aperture of 1 ′′ , and the results are well consistent with the values given in Table 1.
The z ′ -band photometry is done within a 2 ′′ circular aperture in diameter using SExtractor. Aperture cor- Figure 1. Transmission curves of the filters that have been used to measure the UV-continuum slopes of the galaxies. The magenta and cyan spectra represent two model LAEs at z = 6.0 and 6.5, respectively. The Lyα emission line enters the z ′ band for z ≃ 6 galaxies, and enters the Y105 and J110 bands for z ≃ 6.5 galaxies.
rections are derived from bright stars. All galaxies except ID61 are compact in the HST images, so they are point-like sources in the ground-based z ′ -band images. Therefore, the usage of a circular aperture and the aperture corrections from stars are robust. The results are shown in Column 6 of Table 1. ID61 is likely an interacting system (see also Jiang et al. 2013a), and its UV slope is not blue (Table 1).

UV continuum slopes
We measure UV-continuum slopes by fitting a powerlaw (f λ ∝ λ β ) to the photometric data above. The procedure is the same as what we did in Jiang et al. (2013a). Because AB magnitude m AB and log(λ) follow the linear relation m AB ∝ (β + 2) × log(λ), we actually perform a linear fit on m AB and log(λ) with the m AB errors included, i.e., we do a standard (weighted) least-square linear regression by minimizing the χ 2 error statistic. The errors propagate to the final parameters. The results are listed in Table 1. They are also illustrated in Figure 2. All galaxies except ID61 have steep slopes with β ≤ −2.6.
The Lyα emission line enters the z ′ band for z ≃ 6 galaxies, and enters the Y 105 and J 110 bands for z ≃ 6.5 galaxies (see Figure 1). Hence, the Lyα emission increases the relevant broadband flux. The Lyα line flux of our galaxies was provided in Jiang et al. (2013a). Note that we do not use the z ′ band for z ≃ 6.5 galaxies. On the other hand, the flux blueward of Lyα is nearly completely absorbed by the IGM at such high redshifts, which decreases the relevant broadband flux. These effects have been taken into account in this paper. The Figure 2. Measurement of the UV-continuum slopes for 6 galaxies with β ≤ −2.6. The color-coded symbols with 1σ error bars represent the photometric data points. The data points for ID63 and ID43 have been shifted for clarity. The first data points (open circles) have been corrected for the IGM absorption and Lyα emission. The dashed lines are the best fits to the data in 3-4 bands. data points shown in Figure 2 have also been corrected for the IGM absorption and Lyα emission. The details are as follows. To correct for the IGM absorption, we use the model in McGreer et al. (2013). Figure 1 shows two LAE model spectra with the IGM absorption at z = 6.0 and 6.5. Since the Lyα emission line of a z ≈ 6.0 (or 6.5) galaxy is at the edge of the z ′ (or Y 105 ) filter, the uncertainty of the correction for the IGM absorption is negligible.
We use our spectroscopic redshifts and reliable Lyα flux measurements to remove the Lyα contribution in the z ′ -or Y 105 -band photometry before we calculate UV continuum slopes or perform SED modeling. The system transmission, including the filter and detector, is considered. The typical rest-frame equivalent width of Lyα in our sample is about 50Å. The Lyα contribution to the z ′ -band photometry is roughly 20% (10% in Y 105 ), because the two filters are both very wide. Therefore, a 20% uncertainty on the Lyα line flux measurement causes a 4% uncertainty in the z ′ -band photometry, or a 2% uncertainty in Y 105 . We have added an additional error of 0.05 mag in quadrature to the relevant bands to include extra uncertainties introduced in this step.
The two SXDS galaxies ID63 and ID64 are the brightest in our sample, and can be well detected in the K band from the ground. We have adopted the VLT HAWKI-I K s -band photometry from Galametz et al. (2013). The total magnitudes are 25.72 ± 0.17 and 25.76 ± 0.17 mag, respectively. We also find that the Y 105 -band data points of the two objects largely deviate from the linear fits. The two galaxies are bright, compact, and isolated in the HST images, so their photometric measurements are robust. We check the four individual exposures (in four dither positions) for each galaxy, and find that the flux varies by ∼0.5 mag. Due to this unusually large variation, the actual errors (0.13-0.15 mag shown in Table 1) are much larger than the measurement errors (0.05 mag).
ID61 is the only galaxy that obviously has an elongated shape. It seems to be an interacting system, and its minor component is brighter at longer wavelengths. Our current HST images do not have enough resolution to reliably separate the two components. The magnitudes quoted in Table 1 are the total magnitudes, and the slope (β ≃ −1.41) is not steep. We will not discuss this galaxy later.
Among the six blue galaxies, two galaxies ID28 and ID30 have β ≈ −2.6 ± 0.2. Their slopes become flatter compared to our previous measurements (β ≈ −2.8 and -3.5 with σ β ≈ −0.4) in Jiang et al. (2013a). This is mainly due to the significant improvement of the photometry in H 160 . The new photometry shows that the two objects become much 'brighter' in H 160 , resulting in the flatter slopes. If we do not use their z ′ -band data, the slopes become β ≈ −2.9 ± 0.3. Two galaxies ID63 and ID64 respectively have β ≈ −2.9 and -2.8 with σ β ≈ 0.15, which are consistent with our previous measurements β ≃ −2.8 and -3.0, but with reduced errors. The main reason for the consistent results is that the two objects are bright, and the previous measurements were reliable. If we do not use their z ′ -band data, the slopes become β ≈ −2.8±0.2. And if we further exclude the Y 105 -band data (these data are poor as mentioned above), the slopes become β ≈ −3.0 ± 0.3. All these measurements are consistent within 1σ.
The other two galaxies ID07 and ID43 have very steep slopes β ≈ −3.4. But the errors σ β ≈ 0.4 are also larger, because they are relatively faint. For ID43, if we exclude its Y 105 -band data, the slope becomes β ≈ −3.2 ± 0.6. It would be important to improve their photometry or obtain their spectra in future observations. Figure 3 compares our results with previous measurements of β in photometrically selected z ≃ 6 LBGs from HST. Our six galaxies are luminous (by selection) with M 1500 between -22 and -20 mag. In this luminosity range, typical star-forming galaxies have slopes around β ≃ −2 or even redder (Jiang et al. 2013a;

Matthee et al. 2019). Our previous results show that
LAEs have an average slope of β ≃ −2.3 (black square in Figure 3) in the range of −22 < M 1500 < −19 mag. In addition, a small fraction of LAEs have β ≃ −3, but with larger errors. Such extreme slopes were only reported in very faint (J/H ≃ 28 mag or M 1500 ≃ −18.5 mag), photometrically selected LBGs from HST (e.g., Bouwens et al. 2010;Labbé et al. 2010). In this work, we have confirmed the existence of such extreme slopes in luminous and spectroscopically confirmed LAEs. We will discuss their nature in the next section.

DISCUSSION
Generally speaking, higher-redshift galaxies are bluer. For example, Figure 3 shows that z ≥ 6 galaxies have β ≤ −2 on average, while low-redshift galaxies typically have β > −1.5. For LBGs at z > 4, some studies also found that the same redshift-dependence is still valid, and lower-luminosity galaxies tend to be bluer (e.g., Bouwens et al. 2012). These results are controversial (e.g., Dunlop et al. 2013). In addition, the existence of a population of β ≃ −3 LBGs at z ≥ 6 are also in debate, as we introduced in Section 1. A slope of β ≤ −2.7 usually means a very young stellar population with low metallicity and no dust content. But UV colors cannot be arbitrarily blue for a given initial mass function (IMF). Theoretically, it is difficult to explain β ≤ −2.7 using standard stellar population models when nebular emission is included (e.g., Raiter et al. 2010). For example, the bluest slopes are roughly β ≃ −2.6 from the BEAGLE code (Chevallard & Charlot 2016;Williams et al. 2018). The implications of extreme slopes have been discussed in the literature (e.g., Bouwens et al. 2010;Wilkins et al. 2011).

SED modeling
In order to understand the nature of our galaxies with β < −2.6, we model their SEDs and explore the possible ranges of a few stellar population parameters. We include the IRAC data from Jiang et al. (2016). All galaxies except ID07 were detected in at least one IRAC band. We fit the broadband SEDs using the GALEV evolutionary synthesis models (Kotulla et al. 2009). We adopt a Kroupa (similar to Chabrier) IMF with a mass range of 0.1-100 M ⊙ , and assume a constant star formation rate (SFR). Metallicity is fixed to be 0.2 Z ⊙ , a typical value in high-redshift galaxies . We use the Calzetti reddening law (Calzetti et al. 2000) and include metallicity-dependent gaseous or nebular emission. The inclusion of nebular emission is important, because our wavelength range covers very strong emission lines (e.g., Hα, Hβ, or [O iii]) that affect the broadband SEDs. Note that our spectroscopic redshifts remove one critical free parameter 'redshift' and ensure that strong emission lines are covered by proper filters. We mainly constrain three quantities: dust reddening E(B-V), SED age, and stellar mass. The best fits of four example galaxies are illustrated in Figure 4.
Our results show that the dust reddening is consistent with zero and the typical age is ∼4 Myr (the minimum allowed age) for all five galaxies with β ≤ −2.6 and IRAC detections. Two galaxies ID28 and ID30 with β ≈ −2.6 are well fitted. However, it becomes more difficult to fit the bluer galaxies with this model. This is demonstrated in Figure 4, where the best-fitted UV colors are significantly redder than the observed values. The best-fitted slopes β are calculated using the same method as we did for the observed data, i.e., we first compute broadband photometry using the best-fitted spectra, and then estimate β from the broadband photometry (blue crosses in the figure). From the model spectra, we note that the continuum slopes at ≥2 µm (rest-frame ≥3000Å for z ≃ 6) are slightly flatter than those at shorter wavelengths. In addition, only ID63 and ID64 (both at z ≃ 6) have K s -band data. For consistency, we do not use the K s band when we calculate slopes β from the model spectra. Therefore, the restframe wavelength range used here is very similar to the originally defined range of 1250-2600Å (Calzetti et al. 1994). The observed slopes of ID63 and ID64 nearly remain the same without the K s -band data points. . SED modeling of four galaxies using GALEV (Kotulla et al. 2009). The red points with error bars are the observed photometric data points. The horizontal errors indicate the wavelength ranges of the filters. The black spectra represent the best models. The large blue crosses represent the photometric points predicted by the models. The Lyα contributions have been removed from the observed photometry (the first data points). The best-fitted UV slopes (in blue) are significantly redder than the observed values (in red) for ID43, ID63, and ID64.
Previous studies have shown that standard models cannot produce β < −2.7. We adopt a top-heavy IMF by increasing M up to 120 M ⊙ . But this barely changes our results, because very massive stars are efficient ionizing sources that produce large amounts of nebular emission, which results in redder colors. The presence of strong nebular line emission is strongly suggested by our IRAC detections (Figure 4).

Further exploration
We use the CIGALE code (Noll et al. 2009;Serra et al. 2011;Boquien et al. 2019) to explore further possibilities, including lower metallicities, non-zero escape fractions (f esc ) of Lyman continuum flux, and AGN contributions (f AGN ). We modify the code to remove Lyα emission lines since Lyα emission is not a free parameter for our galaxies. We adopt a Chabrier IMF with a mass range of 0.1-100 M ⊙ and assume a constant SFR. We allow age to vary between 1 and 50 Myr, in steps of 1 Myr. We focus on ID63 that has more observed data points than the others. The results are shown in Figure 5. The best-fitted dust extinction is all (nearly) zero, and ages are 1-4 Myr. In panel (a), we try two low metallicities 0.02 and 0.005 Z ⊙ , and find that the IRAC bands are poorly fitted. ID63 is much brighter in IRAC 1 than it is in IRAC 2. Many such galaxies have Figure 5. SED modeling of ID63 using CIGALE (Noll et al. 2009;Serra et al. 2011;Boquien et al. 2019). The meanings of the symbols are the same as Figure 4. For purpose of clarity, emission lines above 24 mag at λ < 2 µm in the model spectra have been cut. The combination of a young population and an old population can potentially explain β ≃ −2.9, shown in (d). The details of the four panels are explained in Section 4.2.
been found, and it is believed that the IRAC 1 flux is boosted by strong [O iii] line emission at this redshift. At a very low metallicity, however, the O/H ratio is reduced, which cannot produce the (IRAC 1 -IRAC 2) color seen in ID63. Bouwens et al. (2010) have used the Schaerer (2003) models to explain β ≃ −3. These models involve some extreme populations such as Population III and those with extremely low metallicities (< 0.001 Z ⊙ ) and IMFs with M up = 500 M ⊙ . When a metallicity is as low as ≤ 0.001 Z ⊙ , β can be as blue as -3.0 over a small age range of 10 ∼ 50 Myr. Bouwens et al. (2010) argue that this may explain some of the galaxies with β ≃ −3. However, it does not likely apply to ID63, since such a low metallicity cannot explain its (IRAC 1 -IRAC 2) color. It is worth pointing out that we have used the same metallicity for the continuum and nebular emission in our models. The gas-phase nebular abundances could be much higher than the stellar metallicity (e.g., Steidel et al. 2016). We have not tried to use difference abundances for the two components.
In the above models, we have assumed the Lyman continuum escape fraction f esc = 0, which is consistent with previous observations that the majority of the lowredshift galaxies observed so far have f esc < 5%. A much higher f esc (≥ 20%) is required to reionize the Universe at z > 6. Meanwhile, a higher f esc can appar-ently make a steeper slope, not only because a higher fraction of very blue stellar light directly leaks from the galaxy, but also because less UV radiation produces less nebular emission. We use f esc = 0.2 with a metallicity of 0.2 Z ⊙ (e.g., Smith et al. 2018), and we show the result in panel (b) of Figure 5. We obtain a better fit with β ≃ −2.7, which can nearly explain the observed β ≃ −2.89 ± 0.15. We further increase f esc to 0.5, and obtain a worse fit, shown in panel (c). This is because there are not sufficient ionizing photons that produce strong nebular line emission to match the IRAC flux. In the simplest case, a line intensity is proportional to n i × n e , where n i is the ion (atom) responsible for the line emission and n e is the electron density. As we mentioned earlier, the IRAC 1 and IRAC 2 flux is mainly boosted by [O iii] and Hα. An increased f esc or a decreased amount of ionizing photons means smaller n i and n e for Hα or smaller n e for [O iii]. Therefore, we cannot arbitrarily increase f esc .
We have not considered AGN contributions so far. Our galaxies are spectroscopically confirmed and do not show type 1 AGN features, such as broad line emission. However, the existing X-ray data are not deep enough to fully rule out AGN activities in such faint objects. We add a small AGN contribution f AGN = 0.1, with an assumption of f esc = 0.2. The AGN model has a central power source with a surrounding dust torus. The spectrum of the central source is a broken-power law, with slope β = 0.2, -1, and -1.5 at 10Å < λ < 300Å, 300Å < λ < 1250Å, and 1250Å < λ < 20 µm, respectively (Fritz et al. 2006). We take the nominal values for the torus parameters (Boquien et al. 2019). This is a coarse estimate, since we do not have observational constraints on the AGN contribution. The result is comparable to the case without the AGN contribution.
Finally, we try the combination of a young and an old stellar populations. This is motivated by the fact that a young population with a high f esc can produce a very steep slope, but with weak IRAC flux (e.g., in case of ID63). The IRAC flux can potentially be boosted by an old population. Recent studies have shown that two populations may exist at z ≥ 6, including an evolved population and a recent starburst likely induced by galaxy merging (e.g., Hashimoto et al. 2018). Given the limited number of the data points, we are not able to fully explore the parameter space. Instead, we make a simple case that includes a 1-Myr old, dust-free population with f esc = 0.6, and a 650-Myr old population with f esc = 0 and E(B-V)=0.5. In this case, the young population dominates the wavelength range blueward of the Balmer break, and the contribution of the old population is negligible in this range. Both populations con-tribute to the IRAC flux redward of the Balmer break. We obtain a reasonable result with β ≃ −2.9 for ID63, shown in Figure 5 (d). Although this is not a real multipopulation SED modeling, it provides a hint that the combination of different populations may be a solution to produce extreme slopes. A more comprehensive study is needed to physically explore this possibility.
For ID43 with β ≃ −3.4±0.4, we are not able to find a reasonable fit to explain its extreme slope. Based on the above models with single populations, the bluest slopes that we can achieve is β ≃ −2.7 when nebular emission is included. In the two-population model, we can obtain β ≃ −2.9. Steeper slopes cannot be well explained by these models. On the other hand, the measurement errors are still large. Deeper observations are needed to reduce the errors.

SUMMARY
We have carried out a study of six very blue LAEs at z ≃ 6. These galaxies are luminous and spectroscopically confirmed, selected from a large galaxy sample that has previous HST near-IR and Spitzer IRAC observations. In HST Cycle 25, we obtained new WFC3 images. The combination of the new data and existing data have allowed us to reliably measure the UV-continuum slopes using 3-4 photometric data points. The slopes are in the range of −3.4 < β < −2.6, with a typical error of 0.2-0.4. In particular, four galaxies have β ≤ −2.8. Such blue, luminous, spectroscopically confirmed galaxies have not been reported previously. We performed SED modeling using the multi-band data. The two β ≈ −2.6 galaxies can be well fitted using standard stellar populations with very young (a few Myr) ages and zero dust content. It is more difficult to fit bluer galaxies using single stellar populations. The bluest slopes that the models predict are β ≃ −2.7. We note that our IRAC detections have put strong constraints on f esc and metallicity, such that a larger f esc cannot produce strong nebular line emission as evidenced by the IRAC flux, and an extremely low metallicity cannot produce the observed (IRAC 1 -IRAC 2) color for ID63. The inclusion of a small AGN contribution does not change our conclusion. Finally, we found that the combination of a young and an old populations can potentially explain the observed SEDs of the β ≃ −2.9 galaxies. In this case, the young population has zero dust content and a high f esc = 0.6, which produces a very blue slope and dominates the β measurement. The old, dusty population with a strong Balmer break contributes little to β. Both young and old populations contribute to the IRAC flux. This model is not able to produce UV slopes as blue as β = −3.4 ± 0.4 in the remaining two galaxies.
In order to fully understand these galaxies, further spectroscopic observations are required. In particular, they are expected to have strong nebular emission lines in the mid-IR. Given the brightness of these galaxies, JWST will be able to directly detect their continuum and line emission in the near future.