Early Low-Mass Galaxies and Star-Cluster Candidates at z~6-9 Identified by the Gravitational Lensing Technique and Deep Optical/Near-Infrared Imaging

TL;DR

This study identifies and analyzes very faint high-redshift galaxies with low stellar masses using gravitational lensing and deep imaging, revealing their properties, stellar mass functions, and potential star cluster candidates at z~6-9.

Contribution

It provides the first detailed stellar population analysis of ultra-faint galaxies at z~6-9, extending the galaxy stellar mass function down to 10^6 solar masses and identifying potential star cluster candidates.

Findings

Identified 357 faint galaxies with stellar masses 10^6-10^9 solar masses.

Derived galaxy stellar mass functions extending to 10^6 solar masses.

Found two compact objects possibly representing star clusters at high redshift.

Abstract

We present very faint dropout galaxies at z~6-9 with a stellar mass M* down to M*~10^6Mo that are found in deep optical/near-infrared (NIR) images of the full data sets of the Hubble Frontier Fields (HFF) program in conjunction with deep ground-based and Spitzer images and gravitational lensing magnification effects. We investigate stellar populations of the HFF dropout galaxies with the optical/NIR photometry and BEAGLE models made of self-consistent stellar population synthesis and photoionization models, carefully including strong nebular emission impacting on the photometry. We identify 357 galaxies with M*~10^6-10^9Mo, and find that a stellar mass to UV luminosity L_UV ratio M*/L_UV is nearly constant at M*~10^6-10^9Mo. Our best-estimate M*/L_UV function is comparable to a model of star-formation duration time of 100 Myr, but 2-7 times higher than the one of 10 Myr assumed in a previous study (at the 5sigma level) that would probably underestimate M* of faint galaxies. We derive the galaxy stellar mass functions (GSMFs) at z~6-9 that agree with those obtained by previous studies with no M*/L_UV assumptions at M*>~10^8Mo, and that extends to M*~10^6Mo. Estimating the stellar mass densities rho* with the GSMFs, we find that rho* smoothly increases from log(rho*/[Mo Mpc^(-3)])=5.91 +0.75/-0.65 at z~9 to 6.21 +0.39/-0.37 at z~6-7, which is consistent with the one estimated from star-formation rate density measurements. In conjunction with the estimates of the galaxy effective radii R_e on the source plane, we have pinpointed two objects with low stellar masses (M*<=10^7Mo) and very compact morphologies (R_e<=40 physical pc) that are comparable with those of globular clusters (GCs) in the Milky Way today. These objects are candidates of star clusters that should be a part or a dominant component of high-redshift low-mass galaxy, some of which may be related to GCs today.

Figures

Tables

Table 1. Table 1: Summary of the subsamples.

Subsample	Threshold	$N$
$z\sim 6-7,M_{\mathrm{UV}}=-21.3$	$M_{\mathrm{UV}}\le -21.0$	$4$
$-20.6$	$-21.0\le M_{\mathrm{UV}}\le -20.5$	$7$
$-20.2$	$-20.5\le M_{\mathrm{UV}}\le -20.0$	$19$
$-19.7$	$-20.0\le M_{\mathrm{UV}}\le -19.5$	$43$
$-19.2$	$-19.5\le M_{\mathrm{UV}}\le -19.0$	$58$
$-18.8$	$-19.0\le M_{\mathrm{UV}}\le -18.5$	$64$
$-18.3$	$-18.5\le M_{\mathrm{UV}}\le -18.0$	$30$
$-17.7$	$-18.0\le M_{\mathrm{UV}}\le -17.5$	$18$
$-17.3$	$-17.5\le M_{\mathrm{UV}}\le -17.0$	$15$
$-16.6$	$-17.0\le M_{\mathrm{UV}}\le -16.0$	$10$
$-15.1$	$-16.0\le M_{\mathrm{UV}}$	$5$
$z\sim 8,M_{\mathrm{UV}}=-20.4$	$M_{\mathrm{UV}}\le -20.0$	$13$
$-19.7$	$-20.0\le M_{\mathrm{UV}}\le -19.5$	$10$
$-19.3$	$-19.5\le M_{\mathrm{UV}}\le -19.25$	$10$
$-19.1$	$-19.25\le M_{\mathrm{UV}}\le -19.0$	$11$
$-18.7$	$-19.0\le M_{\mathrm{UV}}$	$10$
$z\sim 9,M_{\mathrm{UV}}=-20.3$	$M_{\mathrm{UV}}\le -19.9$	$8$
$-19.6$	$-19.9\le M_{\mathrm{UV}}\le -19.5$	$9$
$-19.3$	$-19.5\le M_{\mathrm{UV}}\le -19.1$	$9$
$-18.7$	$-19.1\le M_{\mathrm{UV}}$	$10$

Table 2. Table 2: Best-fit parameters of the M ⋆ subscript 𝑀 ⋆ M_{\star} - M UV subscript 𝑀 UV M_{\mathrm{UV}} relations.

$z$	$a_0$	$a_1$
$6-7$	${8.67}_{-0.11}^{+0.11}$	$-{0.47}_{-0.09}^{+0.09}$
$8$	${8.59}_{-0.19}^{+0.18}$	$-0.47$ (fixed)
$9$	${8.76}_{-0.15}^{+0.13}$	$-0.47$ (fixed)

Table 3. Table 3: Best-fit parameters of the Schechter function of the GSMFs.

Ref.	$M_{\star }$	$\alpha$	${\Psi }^{\ast }$
	$\left[\log M_{\odot }\right]$		$\left[{10}^{-5}{\mathrm{dex}}^{-1}{\mathrm{Mpc}}^{-3}\right]$
$z\sim 6$
K19	${9.58}_{-0.15}^{+0.23}$	$-{1.85}_{-0.07}^{+0.07}$	${18.2}_{-9.3}^{+12.7}$
D14	${10.87}_{-1.06}^{+1.13}$	$-{2.00}_{-0.40}^{+0.57}$	${1.4}_{-1.4}^{+41.1}$
G15	$10.49\pm 0.32$	$-1.55\pm 0.19$	${6.19}_{-4.57}^{+13.50}$
S16	${10.72}_{-0.30}^{+0.29}$	$-{1.91}_{-0.09}^{+0.09}$	${1.35}_{-0.75}^{+1.66}$
B18	${10.29}_{-0.56}^{+0.48}$	$-{1.93}_{-0.07}^{+0.05}$	${5.63}_{-3.23}^{+7.12}$
$z\sim 7$
D14	$10.51$ (fixed)	$-{1.89}_{-0.61}^{+1.39}$	${3.6}_{-3.5}^{+30.1}$
G15	$10.69\pm 1.58$	$-1.88\pm 0.36$	${0.57}_{-0.57}^{+59.68}$
S16	${10.78}_{-0.28}^{+0.29}$	$-{1.95}_{-0.18}^{+0.18}$	${0.53}_{-0.38}^{+1.10}$
B18	${10.25}_{-0.49}^{+0.45}$	$-{1.95}_{-0.07}^{+0.06}$	${2.82}_{-1.88}^{+5.16}$
$z\sim 8$
K19	${8.93}_{-0.15}^{+0.23}$	$-{1.52}_{-0.26}^{+0.27}$	${45.7}_{-25.7}^{+33.7}$
S16	${10.72}_{-0.29}^{+0.29}$	$-{2.25}_{-0.35}^{+0.72}$	${0.035}_{-0.030}^{+0.246}$
B18	${10.48}_{-0.75}^{+0.90}$	$-{2.25}_{-0.29}^{+0.23}$	${0.089}_{-0.074}^{+0.36}$
$z\sim 9$
K19	${9.04}_{-0.18}^{+0.47}$	$-{1.55}_{-0.30}^{+0.29}$	${60.3}_{-46.8}^{+37.5}$
B18	${10.45}_{-0.85}^{+0.80}$	$-{2.33}_{-0.39}^{+0.30}$	${0.06}_{-0.06}^{+0.66}$

Table 4. Table 4: Galaxy stellar mass densities.

$z$	$\log {\rho }_{\star }\left[M_{\odot }{\mathrm{Mpc}}^{-3}\right]$
$6-7$	${6.21}_{-0.37}^{+0.39}$
$8$	${5.64}_{-0.51}^{+0.52}$
$9$	${5.91}_{-0.65}^{+0.75}$

Table 5. Table 5: Catalog of the early GC candidates.

ID	R.A. (J2000)	decl. (J2000)	$M_{\mathrm{UV}}$ [mag]	$M_{\star }\left[{10}^6{M}_{\odot }\right]$	$R_e$ [physical pc]	$\mu$
HFF5C-4260-1364	$342.177526$	$-44.526786$	$-15.7\pm 0.2$	$7.7\pm 7.1$	${40}_{-20}^{+30}$	19.12
HFF5C-4039-1566	$342.168321$	$-44.532412$	$-14.2\pm 0.1$	$1.6\pm 2.0$	${40}_{-20}^{+30}$	77.69

Equations

$$\text{fNMAD}\equiv 1.4826\times \mathrm{median}|{\mu }_{\mathrm{m}}(\mathbf{\theta })-\tilde{\mu }|/\tilde{\mu }.$$

$$M_{\mathrm{UV}}=m_{\mathrm{UV}}-5\log \left[\frac{d_{\mathrm{L}}(z)}{10\mathrm{pc}}\right]+2.5\log [\mu (z)]+2.5\log (1+z),$$

$$z_{\text{phot}}^{\text{median}}-0.1$$

$$1\mathrm{Myr}$$

$$-3$$

$$6$$

$$-2$$

$$0$$

$$P(\mathbf{\Theta }|\mathbf{D},\mathsf{H})=\frac{P(\mathbf{\Theta }|\mathsf{H})P(\mathbf{D}|\mathbf{\Theta },\mathsf{H})}{\int \mathrm{d}\mathbf{\Theta }P(\mathbf{\Theta }|\mathsf{H})P(\mathbf{D}|\mathbf{\Theta },\mathsf{H})}$$

$$\ln \mathcal{L}(\mathbf{\Theta })=-\frac{1}{2}\sum _k{\left[\frac{f_k-{\hat{f}}_k(\mathbf{\Theta })}{{\sigma }_k}\right]}^2.$$

$${\sigma }_k\equiv \sqrt{({\sigma }_k^{\text{obs}}{)}^2+({\sigma }_0{f}_k{)}^2},$$

$$\log (M_{\star }/{\mathrm{M}}_{\odot })=a_0+a_1(M_{\mathrm{UV}}+19.5),$$

$$\log (M_{\star }/{\mathrm{M}}_{\odot })=8.0-0.4(M_{\mathrm{UV}}+19.5)$$

$$\log (M_{\star }/{\mathrm{M}}_{\odot })=8.6-0.4(M_{\mathrm{UV}}+19.5),$$

$$\Psi (M_{\star })=\ln (10){\Psi }^{\ast }{\left(\frac{M_{\star }}{M^{\ast }}\right)}^{\alpha +1}\exp \left(-\frac{M_{\star }}{M^{\ast }}\right).$$

Full text

galaxies: formation — galaxies: high-redshift — galaxies: luminosity function, mass function — globular clusters: general — gravitational lensing: strong

††journal: ApJ

1 Introduction

The stellar mass of a galaxy is an indicator of the masses that are aggregated through the previous star formation and merging processes. Various observations and simulations have found that the stellar masses correlate with basic properties of galaxies, such as the star formation rate (SFR; e.g., McLure et al., 2011; Speagle et al., 2014; Salmon et al., 2015; Santini et al., 2017; Iyer et al., 2018), age (e.g., Sparre et al., 2015), metallicity (e.g., Kojima et al., 2017; Barber et al., 2018), and size (Trujillo et al., 2004; van der Wel et al., 2014; Lange et al., 2015). These facts imply that the stellar mass plays a critical role in understanding the formation and evolution of galaxies in the early universe, which still are major open questions in astronomy today.

The number density of galaxies per stellar-mass ($M_{\star}$) interval, i.e., the galaxy stellar mass function (GSMF), is often used to study the evolution of the total stellar mass in the universe. The low-mass ends of the GSMFs are especially worth investigating at high redshift, because low-mass galaxies are expected to be dominant in the early universe according to the hierarchical cosmology. It is also worth mentioning that the low-mass galaxies ($M_{\star}\sim 10^{6}-10^{7}~{}{\mathrm{M}_{\odot}}$) at $z\sim 6-9$ are predicted to grow up to have stellar masses comparable to that of the Milky Way at $z\sim 0$ (e.g., Behroozi et al., 2013, 2018).

However, the low-mass ends of the GSMFs at high redshift are poorly constrained due to the lack of the sample of high-redshift low-mass galaxies caused by insufficient depths of observations. These observational limits can be resolved by the capability of the Wide Field Camera 3 (WFC3/IR) on board the Hubble Space Telescope (HST) together with the strong gravitational lensing effect caused by galaxy clusters. In this manner, the Hubble Frontier Fields (HFF; Coe et al., 2015; Lotz et al., 2017) spent more than 800 orbits of the HST to survey six galaxy cluster fields, extending the faint end of the rest-frame ultraviolet (UV) luminosity functions (LFs) to a intrinsic UV magnitude $M_{\mathrm{UV}}\approx-14~{}\mathrm{mag}$ at $z\gtrsim 6$ (Finkelstein et al., 2015; Laporte et al., 2016; McLeod et al., 2016; Livermore et al., 2017; Ishigaki et al., 2018; Atek et al., 2018; Bhatawdekar et al., 2018). Here, one should study the GSMFs complementary to the UV LFs that represent the fundamental physical quantity of the stellar mass.

The lensing effect has brought further benefit to studies of galaxy structures and morphologies. Lensed galaxies are stretched along critical curves, allowing the structures to be studied at high spatial resolution. Note that this can be fully exploited with the high resolving power of the WFC3/IR. Owing to magnification by the lensing effects, understanding the size evolution has been significantly advanced by recent studies (e.g., Shibuya et al., 2015; Bouwens et al., 2017; Kawamata et al., 2018). Bouwens et al. (2017) compare the effective radii $R_{\mathrm{e}}$ and the stellar masses $M_{\star}$ of the $z\sim 6-8$ galaxies to study the analogy of the high-redshift galaxies to the local stellar systems, reporting that some galaxies have the values of $R_{\mathrm{e}}$ and $M_{\star}$ comparable to those of the globular clusters (GCs), super star complexes, or star cluster complexes in the local universe. Here, Bouwens et al. (2017) assume $M_{\star}$-$M_{\mathrm{UV}}$ relations, i.e. stellar mass to UV luminosity $L_{\mathrm{UV}}$ ratios $M_{\star}/L_{\mathrm{UV}}$, to estimate $M_{\star}$ with $M_{\mathrm{UV}}$ measurements. In Bouwens et al. (2017), the assumed $M_{\star}$-$M_{\mathrm{UV}}$ relations are based on two representative models with star-formation duration time of 10 and 100 Myr (including the time difference by a factor of $\sim 4$) that have not been tested with observational data yet. The $M_{\star}$-$M_{\mathrm{UV}}$ relation should be determined by observations with no such assumptions, and the $R_{\mathrm{e}}$-$M_{\star}$ relation should be studied.

The structure of this paper is as follows. We describe the data and sample in detail in Section 2. The $M_{\star}$-$M_{\mathrm{UV}}$ relations are obtained by the staking analysis and the spectral energy distribution (SED) fitting technique in Section 3. In Section 4, we show the results and the discussions on the GSMFs and the $R_{\mathrm{e}}$-$M_{\star}$ relations of the galaxies. Finally, we give a summary in Section 5.

Throughout this paper, we adopt a cosmology with $\Omega_{\mathrm{m},0}=0.3$, $\Omega_{\Lambda,0}=0.7$, and $H_{0}=70~{}\mathrm{km}~{}\mathrm{s}^{-1}~{}\mathrm{Mpc}^{-3}$. Magnitudes are given in the AB system (Oke & Gunn, 1983). Densities and sizes of the galaxies are measured in comoving and physical scales, respectively. We adopt the Chabrier (2003) initial mass function (IMF) in a mass range of $0.1-100~{}{\mathrm{M}_{\odot}}$ to estimate stellar masses. In this paper, all of the stellar masses taken from the previous studies are converted to those estimated with the Chabrier (2003) IMF.

2 Data and Samples

2.1 HST Data and Samples

In our analysis, we make use of the image mosaics obtained in the HFF program, which targets six cluster fields—Abell 2744, MACS J0416.1$-$2403, MACS J0717.5$+$3745, MACS J1149.6$+$2223, Abell S1063, and Abell 370 (hereafter A2744, M0416, M0717, M1149, A1063, and A370, respectively)—and their accompanying six parallel fields. All of the 12 fields were observed with the three bands of the Advanced Camera for Surveys (ACS) and four bands of the WFC3/IR; F435W ($B_{435}$), F606W ($V_{606}$), F814W ($i_{814}$), F105W ($Y_{105}$), F125W ($J_{125}$), F140W ($J\!H_{140}$), and F160W ($H_{160}$). We utilize the drizzled and weight images that were produced by Shipley et al. (2018) in the manner summarized below. 111http://cosmos.phy.tufts.edu/~danilo/HFF/Download.html First, the HFF v1.0 images were downloaded from the MAST archive. 222http://www.stsci.edu/hst/campaigns/frontier-fields/ The point-spread functions (PSFs) of these images were homogenized to those of the $H_{160}$ images. The PSF FWHM of the homogenized images is $\approx 0\farcs 18$. Second, the bright cluster galaxies (bCGs) were modeled and subtracted from the images to avoid the diffuse intracluster light (ICL) in photometry of the background faint sources. All of the images have a pixel scale of $0\farcs 06$. We correct for the Galactic extinction using the values given by Shipley et al. (2018).

We divide each image into $3\times 3$ grid cells, and measure the limiting magnitude in each cell ($\sim 1~{}\mathrm{arcmin}^{2}$). This is because limiting magnitudes are not homogeneous due to the ICL (e.g., Montes & Trujillo, 2014; Ishigaki et al., 2015; Kawamata et al., 2016). The $5\sigma$ limiting magnitudes in the $H_{160}$ band images are $\approx 28.4-29.2~{}\mathrm{mag}$ in a $0\farcs 35$-diameter circular aperture.

We use the galaxy sample selected by Kawamata et al. (2018), which consists of 350 $i$-, 64 $Y$-, and 39 $Y\!J$-dropouts (or $z\sim 6-7,~{}z\sim 8,~{}z\sim 9$ Lyman break galaxies, respectively). Photometry is reperformed in the following way. First, we measure the aperture magnitude $m_{\text{ap}}$ with a diameter of $D_{\text{ap}}\equiv 0\farcs 35$ at the position of the dropouts, using the IRAF task phot (Tody, 1986, 1993). Second, we apply an aperture correction. To evaluate the aperture correction term, we create a median-stacked $J_{125}$-band image of the $i$-dropouts with PSF homogenization, and measure the aperture flux of the stacked dropout with changing the aperture diameter. We define $D_{99}$ by the diameter of the aperture that includes $\geq 99~{}\%$ of the total flux is included. The aperture correction term, $c_{\text{ap}}$, is defined as $c_{\text{ap}}\equiv m_{\text{ap}}(D_{\text{ap}})-m_{\text{ap}}(D_{99})$, where $m_{\text{ap}}(d)$ is a aperture magnitude for a given diameter $d$. For $B_{435}$ and $V_{606}$ bands, $c_{\text{ap}}$ is defined as the same value defined in the $i_{814}$ band, because $c_{\text{ap}}$ cannot be estimated due to low signal-to-noise (S/N) ratios in the $B_{435}$ and $V_{606}$ bands. The values of $c_{\text{ap}}$ are $0.98,~{}0.98,~{}0.98,~{}0.93,~{}0.95,~{}0.94$, and $0.90~{}\mathrm{magnitudes}$ in the $B_{435}$, $V_{606}$, $i_{814}$, $Y_{105}$, $J_{125}$, $J\!H_{140}$, and $H_{160}$ bands, respectively. Finally, we estimate the total magnitude $m_{\text{tot}}$ as $m_{\text{tot}}=m_{\text{ap}}(D_{\text{ap}})-c_{\text{ap}}$.

In our photometry, some dropouts may be detected at the bands blueward of the Lyman break. We remove the dropouts with a $\geq 2\sigma$ detection either in $B_{435}$ or $V_{606}$ band for $i$-dropouts, and in $B_{435}$, $V_{606}$, or $i_{814}$ band for $Y$- and $Y\!J$-dropouts. For our dropouts, we apply criteria of $\geq 3\sigma$ detections in $J_{125}$, $J\!H_{140}$, and $H_{160}$ bands for $i$-, $Y$-, and $Y\!J$-dropouts, respectively. Our final sample consists of 267 $i$-dropouts, 54 $Y$-dropouts, and 36 $Y\!J$-dropouts (357 in total).

2.2 VLT and Keck Data

We use the $K_{s}$ band images obtained by the K-band Imaging of the Frontier Fields (KIFF) program (Brammer et al., 2016). The deep $K_{s}$ images of the VLT/HAWK-I (Keck-I/MOSFIRE) are available in both the cluster and the parallel fields of A2744, M0416, A1063, and A370 (M0717 and M1149). The HAWK-I and MOSFIRE images have the PSF FWHM of $\approx 0\farcs 4-0\farcs 5$. We utilize the drizzled and weight images of Shipley et al. (2018), whose details are described in Section 2.1. The pixel scale of the images is $0\farcs 06$ that matches to that of the HST images. Photometry and limiting-magnitude measurements are conducted in the same manner as those in Section 2.1, but with $D_{\text{ap}}=0\farcs 6$ and $c_{\text{ap}}=0.98$. The $5\sigma$ limiting magnitudes are $\approx 25.3-26.2~{}\mathrm{mag}$.

2.3 Spitzer Data

We utilize the drizzled and weight images of ch1 ($3.6~{}\mu\mathrm{m}$) and ch2 ($4.5~{}\mu\mathrm{m}$) of the IRAC. The photometric data that were taken by 2016 December are combined by Shipley et al. (2018), who also applies bCG subtraction and Galactic extinction correction. In both ch1 and ch2, the PSF FWHM is $\approx 1\farcs 7-2\farcs 0$, and pixel scales are $0\farcs 3$. Photometry and limiting-magnitude measurements are conducted in the same manner as described in Section 2.1, but with $D_{\text{ap}}=3\farcs 0$ for both ch1 and ch2, and $c_{\text{ap}}=0.52,~{}0.55$ for ch1 and ch2, respectively (Ono et al., 2010). The $5\sigma$ limiting magnitudes are $\approx 24.5-25.5~{}\mathrm{mag}$.

2.4 Lens Models

We apply the best-fit magnification factors $\mu$ to estimate the intrinsic magnitudes and stellar masses of the dropouts. The magnification factors that we use are calculated with the glafic (Oguri, 2010) parametric models by Kawamata et al. (2016, 2018). 333https://archive.stsci.edu/prepds/frontier/lensmodels/ Priewe et al. (2017) estimate the uncertainties of the eight mass models of M0416, using a fractional normalized median absolute deviation, which is defined as

[TABLE]

Here $\mu_{\mathrm{m}}(\bm{\theta})$ is the magnification factor calculated with a given model m at a given position $\bm{\theta}$, and $\tilde{\mu}$ is the median of the magnification factors. It is found that $\text{fNMAD}\approx 0.3~{}(0.7)$ at $\tilde{\mu}\approx 2~{}(40)$ in M0416, which implies that magnification factors of different models differ only by a factor of $\approx 2$. The magnification uncertainties do not change our conclusions (see also Ishigaki et al. 2015; Meneghetti et al. 2017; Kawamata et al. 2018).

3 Stellar Population of the High-$z$ Galaxies

3.1 Stacking Analysis

To create images with high S/N ratios, we conduct a stacking analysis in the following way. We take $J_{125}$, $J\!H_{140}$, and $H_{160}$ band magnitudes as the rest-frame UV apparent magnitudes $m_{\mathrm{UV}}$ of the $i$-, $Y$-, and $Y\!J$-dropouts, respectively. The intrinsic absolute magnitude of a dropout is calculated via

[TABLE]

where $d_{\mathrm{L}}(z)$ is the luminosity distance to the dropout, and $\mu(z)$ is the best-fit magnification factor at the position of the dropout (see Section 2.4). We apply $z=6$, $8$, and $9$ for the $i$-, $Y$-, and $Y\!J$-dropouts, respectively, for consistency with the assumption in the size measurements (Kawamata et al., 2018). We divide our dropout samples into subsamples by the values of $M_{\mathrm{UV}}$ at $z\sim 6-7$, $8$, and $9$. A summary of the subsamples is shown in Table 1.

In each band, we cut out images centered at the positions of dropouts, and divide the pixel counts by $\mu$ of the dropouts. Here we assume that the values of $\mu$ do not change within sizes of the dropout galaxies, because the sizes of the dropout galaxies are too small to significantly impact on the dropout galaxy fluxes. We then median-stack the images of the dropouts for each subsample with iraf task imcombine.

Figures 1, 2, and 3 show the stacked images of the $z\sim 6-7$, $8$, and $9$ subsamples, respectively. Magnitudes for the subsamples are measured in the same manner as those for the individual dropouts in Section 2.1. To estimate uncertainties of the total fluxes (or the total magnitudes), we first follow the three steps; 1) randomly selecting positions in the sky area of the grid cell where the subsamples are located, 2) generating median-stacked sky noise images, and 3) performing aperture photometry on the sky noise images with the aperture correction to obtain the total flux $f_{i}$. The steps 1)-3) are repeated for 100 times. We make a histogram of $f_{i}$, and fit the histogram with a Gaussian profile. We regard the standard deviation divided by the median magnification factor of the subsample as the uncertainty of the total flux.

3.2 SED Modeling

To investigate the typical stellar mass $M_{\star}$ for a given $M_{\mathrm{UV}}$, we conduct the SED fitting method. We use the BEAGLE tool (Chevallard & Charlot, 2016), which is based on a recent version of the stellar population models of Bruzual & Charlot (2003) and the photoionization models of Gutkin et al. (2016) that are computed with CLOUDY (Ferland et al., 2013). The intergalactic medium (IGM) absorptions follow the models of Inoue et al. (2014). We adopt the Calzetti et al. (1994) law to the models for dust attenuation.

There are six free parameters in the fitting: i) redshift $z$ of the galaxy, ii) galaxy age $t_{\text{age}}$, iii) galaxy-wide ionization parameter $U_{\text{ion}}$, iv) total mass of the formed stars $M_{\star}$, v) stellar metallicity $Z$, and vi) $V_{606}$-band dust attenuation optical depth $\tau_{\mathrm{V}}$. We assume uniform prior probability distribution functions (PDFs) in the range of

[TABLE]

for the parameters i), ii), iii), iv), v), and vi), respectively. Here, $z_{\text{phot}}^{\text{median}}$ is the median of the photometric redshifts $z_{\text{phot}}$ of the subsample whose values are taken from Kawamata et al. (2018). We omit $z_{\text{phot}}$ values that are lower than $4$ to derive $z_{\text{phot}}^{\text{median}}$ in the same manner as the estimations of the median photometric redshifts Kawamata et al. (2018). The age of the universe at $z$ is represented by $t_{\text{univ}}(z)$. We assume a constant star formation history. The interstellar medium metallicity is assumed to equal the stellar metallicity. The dust-to-metal ratio is fixed to 0.3 (e.g., De Vis et al., 2017).

The posterior PDF $P(\bm{\Theta}|\bm{D},\mathsf{H})$ of a given parameter set $\bm{\Theta}=\{z,t_{\text{age}},U_{\text{ion}},M_{\star},Z,\tau_{\mathrm{V}}\}$ of a model $\mathsf{H}$ is calculated based on the Bayes’ theorem

[TABLE]

(e.g., Jeffreys, 1961), where $\bm{D}$ is the data set, i.e., the fluxes in the $B_{435}$, $V_{606}$, $i_{814}$, $Y_{105}$, $J_{125}$, $J\!H_{140}$, $H_{160}$, $K_{s}$, ch1, and ch2 bands. The likelihood function of $\bm{\Theta}$, $\mathcal{L}(\bm{\Theta})\equiv P(\bm{D}|\bm{\Theta},\mathsf{H})$, is defined via

[TABLE]

Here $f_{k}$, $\hat{f}_{k}(\bm{\Theta})$, and $\sigma_{k}$ are the observed flux, the flux predicted by the parameter set $\bm{\Theta}$, and the flux uncertainty, respectively. The subscript $k$ runs over all of the bands but the $i_{814}$, $Y_{105}$, and $J_{125}$ band for $z\sim 6-7$, $8$, and $9$ subsamples, respectively. We do not include the $i_{814}$, $Y_{105}$, and $J_{125}$ band for $z\sim 6-7$, $8$, and $9$ subsamples, respectively. This is because we need to avoid the broadband photometry contaminated by unknown Ly$\alpha$ emission and IGM absorption effects. The value of $\sigma_{k}$ is defined as

[TABLE]

where $\sigma_{k}^{\text{obs}}$ is the observational uncertainty estimated in Section 3.1, and $\sigma_{0}$ is the relative systematical uncertainty, e.g., errors in background subtraction, flux calibration, and model predictions (Brammer et al., 2008; Dahlen et al., 2013; Acquaviva et al., 2015; Chevallard & Charlot, 2016). We define $\sigma_{0}=0.04~{}(0.05)$ for the HST and $K_{s}$ (IRAC) bands, applying the values recommended in the user manual of $\tt BEAGLE$ (Chevallard & Charlot, 2016). If the S/N ratio of the stacked image in the band $k$ is less than $2$, we substitute $f_{k}$ values with [math]. The posterior PDFs are efficiently sampled by the Nested Sampling algorithm implemented in the MULTINEST tool (Feroz & Hobson, 2008; Feroz et al., 2009). We present the best-fit SEDs with the data photometries for $z\sim 6-7$, $8$, and $9$ subsamples in Figure 4, 5, and 6, respectively. An example of the posterior PDFs of the parameters is shown in Figure 7.

3.3 Stellar Mass to UV Luminosity Relations

Figure 8 shows the relations between the stellar mass $M_{\star}$ and the UV magnitude $M_{\mathrm{UV}}$. Here, $M_{\star}$ and $M_{\mathrm{UV}}$ represent the median value of the marginal posterior PDF and the median value of the UV magnitudes of the dropouts in the subsample, respectively. The vertical error bars represent 68 % confidence intervals, while the horizontal bars show the minimum and maximum values of $M_{\mathrm{UV}}$ of the dropouts in the subsamples. Figure 8 indicates that our results are broadly consistent with the previous results.

We fit a linear function

[TABLE]

where the intercept $a_{0}$ and the slope $a_{1}$ are set as free parameters, to our data in the magnitude range of $-21.0\leq M_{\mathrm{UV}}\leq-16.0$. At $z\sim 8$ and $9$, we fix $a_{1}$ to the best-fit value of $a_{1}$ of $z\sim 6-7$, because fitting is unstable due to the lack of data points. The black solid lines in Figure 8 show the best-fit relations whose parameters are listed in Table 2. We find that the best-fit value of $a_{1}$ at $z\sim 6-7$ is $-0.47\pm 0.09$ that is similar to those obtained in the previous studies within the $1\sigma$ uncertainties. This value of $a_{1}$ is consistent with $-0.4$ that is the slope for the constant $M_{\star}/L_{\mathrm{UV}}$ case, indicating that the average $M_{\star}/L_{\mathrm{UV}}$ of dropouts at $z\sim 6-7$ is nearly constant at $M_{\star}\sim 10^{6}-10^{9}~{}{\mathrm{M}_{\odot}}$. The best-fit values of intercept $a_{0}$ are $8.67_{-0.11}^{+0.11}$, $8.59_{-0.19}^{+0.18}$, and $8.76_{-0.15}^{+0.13}$ at $z\sim 6-7$, $8$, and $9$, respectively, showing no evolution beyond the errors. No evolution of $a_{0}$ indicates that stellar populations of the dropouts do not significantly change at $z\sim 6-9$ in the magnitude range.

In Figure 8, the blue dashed and dotted lines denote the $M_{\star}$-$M_{\mathrm{UV}}$ relations (corresponding to $M_{\star}/L_{\rm UV}$ ratios) that are assumed by Bouwens et al. (2017) for their $M_{\star}$ estimations from their $M_{\mathrm{UV}}$ measurements for $z\sim 6-8$ dropouts. Bouwens et al. (2017) assume two representative model stellar populations whose star-formation duration times of 10 and 100 Myr with no stellar mass/UV luminosity dependence (See Section 1 for more details). The assumptions of 10 and 100 Myr correspond to $(a_{0},a_{1})=(8.0,-0.4)$ and $(8.6,-0.4)$, i.e.

[TABLE]

and

[TABLE]

respectively. Figure 8 indicates that our best-fit $M_{\star}$-$M_{\mathrm{UV}}$ relations (the black solid lines) are comparable to the model of 100-Myr star-formation duration time (the blue dotted lines). However, our best-fit $M_{\star}$-$M_{\mathrm{UV}}$ relations fall above the model of 10-Myr star-formation duration time (the blue dashed lines). By the comparison with our best-fit parameters at $z\sim 6-7$, we find that the parameter set for the blue dashed line ($a_{0}=8.0$ and $a_{1}=-0.4$) is ruled out at the $5\sigma$ confidence level. The parameter set of $(a_{0},a_{1})=(8.0,-0.4)$ provides the $M_{\star}$-$M_{\mathrm{UV}}$ relation $2-7$ times lower than the one of our best-fit parameters at $z\sim 6-9$ in $M_{\star}$.

4 Results and Discussion

4.1 Galaxy Stellar Mass Functions

To derive the GSMFs, we apply the best-fit $M_{\star}$-$M_{\mathrm{UV}}$ relations (Section 3.3) to the LFs of Ishigaki et al. (2018). We show our GSMFs in Figure 9 with black dots. One can find that our GSMFs are broadly consistent with those of previous studies in the range of $M_{\star}\sim 10^{7}-10^{9}~{}{\mathrm{M}_{\odot}}$ at $z\sim 6-7$. In the low-mass end, our GSMFs reach $M_{\star}\sim 10^{6}~{}{\mathrm{M}_{\odot}}$, which is lower than those of the previous studies by $\gtrsim 1~{}\mathrm{dex}$.

We parametrize our GSMFs with the Schechter (1976) function,

[TABLE]

The characteristic stellar mass $M^{*}$, the low-mass end slope $\alpha$, and the normalization $\Psi^{*}$ are free parameters whose best-fit values are obtained by the least-squares fitting method. To improve the statistical accuracy, we fit Schechter functions simultaneously to our GSMFs and the GSMFs derived by Song et al. (2016) who do not use the HFF data but the CANDELS/GOODS (Grogin et al., 2011; Koekemoer et al., 2011) and the HUDF (Beckwith et al., 2006; Illingworth et al., 2013) data. We combine our GSMFs at $z\sim 6-7$ and the GSMFs of Song et al. (2016) at $z\sim 6$, because the median value of $z_{\text{phot}}$ of our dropouts is $z=6.3$ that is regarded as $z\sim 6$. Table 3 summarizes the best-fit Schechter parameters at each redshift.

4.2 Galaxy Stellar Mass Density

We derive the galaxy stellar mass densities (GSMDs) at $z\sim 6-9$. We integrate our best-fit Schechter functions of the GSMFs over the stellar mass range of $10^{8}\leq M_{\star}/{\mathrm{M}_{\odot}}\leq 10^{13}$ that is adopted in previous studies. Table 4 lists our GSMDs at $z\sim 6-9$. The uncertainties of our GSMDs include the statistical errors and the cosmic variance uncertainties, the latter of which are added in quadrature with the values of $0.3$, $0.4$ and $0.5~{}\mathrm{dex}$ at $z\sim 6-7$, $8$ and $9$, respectively (Robertson et al., 2014; Ishigaki et al., 2015; Bhatawdekar et al., 2018).

We show our GSMDs in Figure 10 together with the results obtained in the previous studies. Our GSMDs are consistent with those of the previous studies within $1\sigma$ uncertainties.

If the GSMDs and the star-formation rate densities (SFRDs) are accurately measured, the time integral of the SFRDs should be consistent with the GSMDs. In Figure 10, we show the evolution of the GSMDs that is derived as the time integration of the SFRDs obtained in Madau & Dickinson (2014) with the black solid curve. The SFRDs of Madau & Dickinson (2014) declines smoothly by $(1+z)^{-3.9}$ at $z>8$. Similar trends of the SFRD evolution are found in the studies of Finkelstein et al. (2015), McLeod et al. (2016) and Bhatawdekar et al. (2018; see also Ellis et al. 2013 and Madau 2018). However, Oesch et al. (2018) claim that the SFRDs evolve more rapidly by $(1+z)^{-10.9}$ at $z>8$ that is shown with the blue solid curve in Figure 10. There are two scenarios of the SFRD evolutionary trends, the smooth evolution of $(1+z)^{-3.9}$ and the rapid evolution of $(1+z)^{-10.9}$. We calculate $\chi^{2}$ values with our GSMDs for the two functions of the smooth evolution and the rapid evolution. The $\chi^{2}$ value is the smaller for the smooth evolution than for the rapid evolution by a factor of $\approx 2.8$, implying that our results support the smooth evolution of the SFRDs. The smooth evolution of the SFRDs suggests that the star-formation rate per dark-matter mass increases at $z>8$ (Harikane et al., 2018; Oesch et al., 2018).

4.3 Size-Mass Relations

We compare the distribution of the local GCs and our dropouts in the parameter space of $R_{\mathrm{e}}$ and $M_{\star}$. The GCs in the Milky Way today have $R_{\mathrm{e}}$ values up to $\approx 40$ physical pc (e.g. Pal 5 and Pal 14) and $M_{\star}$ up to $\sim 10^{7}~{}{\mathrm{M}_{\odot}}$ (e.g. NGC 5139) that are shown in the catalog compiled by Baumgardt & Hilker (2018). We thus regard the dropouts with $R_{\mathrm{e}}\leq 40$ physical pc and $M_{\star}\leq 10^{7}~{}{\mathrm{M}_{\odot}}$ as early GC candidates.

We plot the $R_{\mathrm{e}}$-$M_{\star}$ values of our dropouts at $z\sim 6-7$, $8$, and $9$ in Figure 11, 12, and 13, respectively. The values of $R_{\mathrm{e}}$ are taken from Kawamata et al. (2018), who obtain the best-fit value by the fit of ellipsoidal Sérsic profiles to the galaxy profiles that are observed with the lensing effects. We measure $M_{\mathrm{UV}}$ of the dropouts in the same manner as those for the stacked images. The values of $M_{\mathrm{UV}}$ are converted to $M_{\star}$ with the best-fit $M_{\star}$-$M_{\mathrm{UV}}$ relations obtained in Section 3.3.

We compare the $R_{\mathrm{e}}$-$M_{\star}$ distributions with those of the stellar systems in the local universe (Norris et al. 2014 and the references therein). We find that two of our dropouts, HFF5C-4260-1364 and HFF5C-4039-1566, 444These IDs are identical to the ones by Kawamata et al. (2018) and Ishigaki et al. (2018). meet the criteria of the early GC candidates. The parameters of the candidates are listed in Table 5.

As stated in the introduction, Bouwens et al. (2017) study the $R_{\mathrm{e}}$-$M_{\star}$ distributions of the $z\sim 6-8$ dropouts. However, Bouwens et al. (2017) use only the first four HFF cluster (A2744, M0416, M0717, and M1149) data that were available at the time of their publication. No studies of $R_{\mathrm{e}}$-$M_{\star}$ distribution for $z\gtrsim 6$ are conducted in the rest of two clusters, A1063 and A370. Our early GC candidates of HFF5C-4260-1364 and HFF5C-4039-1566 are placed in the A1063 field. These early GC candidates are identified for the first time.

In Bouwens et al. (2017), there are 18 $z\sim 6-8$ dropouts meeting the criteria of the early GC candidates with their stellar masses obtained with the assumption of the 100-Myr star-formation duration time (Equation 14). One out of the 18 dropouts, M0416I-6115434445, is also selected in our study that refers to the dropout as HFF2C-1156-3446. This dropout is also reported as a compact star-forming region GC1 in Vanzella et al. (2017). The stellar masses and effective radii obtained in our study, Bouwens et al. (2017), and Vanzella et al. (2017) are broadly consistent within the $1-2\sigma$ uncertainties. Note that another dropout from Bouwens et al. (2017), M0416I-6118103480, is selected as the dropout HFF2C-1181-3480 in the original sample of Kawamata et al. (2018) and Ishigaki et al. (2018), but removed from our sample in Section 2.1, due to the possible detection of the blue continuum flux. Because our detection threshold would be conservative, the rest of the 16 dropouts shown in Bouwens et al. (2017) are not securely detected in our analysis.

The previous studies suggest that the GCs in the Milky Way form almost instantaneously at $z\gtrsim 5$, arguing that the stars inside the GCs have high number densities and tight isochrones on the HR diagram (e.g., Vandenberg et al., 1996; Forbes & Bridges, 2010). The dropouts that we select are thus candidates of star clusters that should be a part or a dominant component of high-redshift low-mass galaxies, some of which may be related to GCs today. Our calculations for observational feasibilities suggest that these candidates can be spectroscopically confirmed with the James Webb Space Telescope and the Extremely Large Telescopes such as the Thirty Meter Telescope.

5 Summary

We present stellar populations of low-mass dropout galaxies at $z\sim 6-9$ down to a stellar mass of $M_{\star}\sim 10^{6}~{}{\mathrm{M}_{\odot}}$ identified with the full data sets of the HFF. The stellar populations are studied with the HST and Spitzer deep optical to NIR images and the cluster-lensing magnification effects. The major results of our study are summarized below.

We investigate the stellar populations of the dropout galaxies with the optical/NIR photometric measurements and the BEAGLE SED modeling tool. We estimate $M_{\star}$ as a function of the UV magnitude $M_{\mathrm{UV}}$ at $M_{\star}\sim 10^{6}-10^{9}~{}{\mathrm{M}_{\odot}}$. We find that our best-estimate $M_{\star}/L_{\mathrm{UV}}$ is consistent with the constant $M_{\star}/L_{\mathrm{UV}}$ within the $1\sigma$ uncertainties. Our best-estimate $M_{\star}/L_{\rm UV}$ function is comparable to a model of star-formation duration time of 100 Myr, but $2-7$ times higher than the one of 10 Myr assumed in a previous study at the $5\sigma$ level. 2. 2.

We apply the best-estimate $M_{\star}$-$M_{\mathrm{UV}}$ relation (corresponding to $M_{\star}/L_{\mathrm{UV}}$) to the LFs of Ishigaki et al. (2018) to derive the GSMFs at $z\sim 6-9$. The GSMFs agree with those obtained by previous studies broadly at $M_{\star}\gtrsim 10^{8}~{}{\mathrm{M}_{\odot}}$, extending to the low-mass limit of $\sim 10^{6}~{}{\mathrm{M}_{\odot}}$. We reach this low-mass limit, exploiting the deep HFF full data sets and the gravitational-lensing magnification effects. 3. 3.

Integrating the GSMFs, we find that the GSMDs $\rho_{\star}$ smoothly increases from $\log(\rho_{\star}/[{\mathrm{M}_{\odot}}~{}\mathrm{Mpc}^{-3}])=5.91_{-0.65}^{+0.75}$ at $z\sim 9$ to $6.21_{-0.37}^{+0.39}$ at $z\sim 6-7$. This trend is consistent with the one estimated from time integration of the SFRD evolution compiled by Madau & Dickinson (2014). 4. 4.

The stellar masses of the dropouts are presented as a function of effective radius $R_{\mathrm{e}}$ with our best-fit $M_{\star}$-$M_{\mathrm{UV}}$ relations. We find that the two dropouts have low stellar masses ($M_{\star}\leq 10^{7}~{}{\mathrm{M}_{\odot}}$) and compact morphologies ($R_{\mathrm{e}}\leq 40$ physical pc), which are comparable with those of GCs in the Milky Way today. These dropouts are candidates of star clusters that should be a part or a dominant component of high-redshift low-mass galaxies. These dropouts may be related to the present GCs, given the fact that the old stellar populations are found in GCs of the Milky Way today.

We are grateful to Michiko Fujii, Seiji Fujimoto, Ryo Higuchi, Ryohei Itoh, Nobunari Kashikawa, Shiro Mukae, and Haibin Zhang for useful comments and discussions. We particularly thank Masafumi Ishigaki and Ryota Kawamata for their sample catalogs and mass models. GB acknowledges financial support through PAPIIT projects IG100319 from DGAPA-UNAM. This work is supported by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan, and KAKENHI (15H02064, 17H01110, 17H01114, 15H05892, and 18K03693) Grant-in-Aid for Scientific Research (A) through Japan Society for the Promotion of Science. This work is based on data and catalog products from HFF-DeepSpace, funded by the National Science Foundation and Space Telescope Science Institute (operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555).