FOREVER22: the first bright galaxies with Population III stars at redshifts z ≃ 10-20 and comparisons with JWST data

概要

MNRAS 525, 4832–4839 (2023)

https://doi.org/10.1093/mnras/stad2497

Advance Access publication 2023 August 25

FOREVER22: the first bright galaxies with Population III stars at
redshifts z  10–20 and comparisons with JWST data
Hidenobu Yajima ,1 ‹ Makito Abe,1 Hajime Fukushima ,1 Yoshiaki Ono,2 Yuichi Harikane,2
Masami Ouchi,2,3,4 Takuya Hashimoto 5 and Sadegh Khochfar6
1 Center

for Computational Sciences, University of Tsukuba, 1-1-1 Tennodai, Tsukuba, Ibaraki 305-8577, Japan
for Cosmic Ray Research, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8582, Japan
3 Kavli IPMU (WPI), The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8583, Japan
4 National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan
5 Tomonaga Center for the History of the Universe (TCHoU), Faculty of Pure and Applied Sciences, University of Tsukuba, Tsukuba, Ibaraki 305-8571, Japan
6 Institute for Astronomy, University of Edinburgh, Royal Observatory, Edinburgh EH9 3HJ, UK
2 Institute

ABSTRACT

We study the formation of the first galaxies in overdense regions modelled by the FORmation and EVolution of galaxies in
Extremely overdense Regions motivated by SSA22 (FOREVER22) simulation project. Our simulations successfully reproduce
the star formation rates and the MUV –Mstar relations of candidate galaxies at z ∼ 10–14 observed by the JWST. We suggest that
the observed galaxies are hosted by dark matter haloes with Mh  1010 M and are in short-period starburst phases. On the
other hand, even simulated massive galaxies in overdense regions cannot reproduce the intense star formation rates and the large
stellar masses of observed candidates at z ∼ 16. Also, we show that the contribution of Population III stars to the ultraviolet (UV)
flux decreases as the stellar mass increases and it is a few per cent for galaxies with Mstar ∼ 107 M . Therefore, a part of the
observed flux by JWST could be the light from Population III stars. Our simulations suggest that the UV flux can be dominated
by Population III stars and the UV slope shows β ࣠ −3 if future observations would reach galaxies with Mstars ∼ 105 M at
z ∼ 20 of which the mass fraction of Population III stars can be greater than 10 per cent.
Key words: stars: Population III – galaxies: evolution – galaxies: formation – galaxies: high-redshift.

1 I N T RO D U C T I O N
Understanding galaxy formation is one of the central issues in current
astrophysics. In particular, the first galaxies at redshifts beyond z = 10
are the most likely drivers of cosmic reionization (Yajima et al. 2009,
2014; Yajima, Choi & Nagamine 2011; Paardekooper, Khochfar &
Dalla Vecchia 2013, 2015; Wise et al. 2014; Arata et al. 2019; Ma
et al. 2020; Rosdahl et al. 2022) and hosts of the first massive black
holes (Regan & Haehnelt 2009; Agarwal et al. 2014; Yajima &
Khochfar 2016; Wise et al. 2019; Latif et al. 2022a). Thus, revealing
the formation of the first galaxies is of great importance. Using
Lyman α lines, a lot of galaxies at z ࣠ 9 have been identified (e.g. Ono
et al. 2012; Shibuya et al. 2012; Finkelstein et al. 2013; Song et al.
2016; Ouchi et al. 2018). However, the transmission of the Lyman α
line is reduced significantly as the neutral degree of the intergalactic
medium increases, resulting in the difficulty of galaxy observation
beyond z ∼ 10 (Yajima, Sugimura & Hasegawa 2018). Recent
submillimetre observations have successfully detected high-redshift
galaxies at z ࣠ 9 via the detections of [C II] 158 μm and [O III] 88 μm
lines (e.g. Capak et al. 2015; Inoue et al. 2016; Hashimoto et al.
2018; Tamura et al. 2019). The metal-line observation is expected to
be difficult if target galaxies exceed z ∼ 10 because of insufficient



E-mail: yajima@ccs.tsukuba.ac.jp

metal enrichment (e.g. Bakx et al. 2023; Popping 2023; Yoon et al.
2023). Therefore, galaxies at z  10 have been investigated with
Lyman-break technique (e.g. Oesch et al. 2013, 2016; Bouwens et al.
2019). However, the number of samples has been limited and the
spectroscopic confirmations have been difficult for the sensitivities
of the telescopes with a reasonable integration time.
These situations are drastically changing with observations by
JWST. Using the data of the first cycle observation by JWST, highredshift galaxies have been identified. Donnan et al. (2023) found 44
new candidate galaxies and estimated the ultraviolet (UV) luminosity
functions at the redshifts z = 8–15. Harikane et al. (2023) found
candidate galaxies at z ∼ 16 with large stellar masses and star
formation rates (SFRs; see also Naidu et al. 2022). Furtak et al. (2023)
indicated that the candidate galaxies at z ∼ 9–16 had properties of
young galaxies with ages ∼10–100 Myr and very blue UV slope
down to β ∼ −3 (see also Topping et al. 2022; Cullen et al. 2023).
The properties of high-redshift galaxies at z > 6 have been
investigated in various simulation projects as CODA (Ocvirk et al.
2016), FLARES (Lovell et al. 2021), THESAN (Kannan et al.
2022), MILLENIUMTNG (Pakmor et al. 2023), UNIVERSEMACHINE (Behroozi et al. 2019), and Santa Cruz model (Gabrielpillai
et al. 2022). These simulations successfully reproduced statistical
properties like luminosity functions of observed galaxies at z ࣠ 8.
Also, some previous works provided theoretical predictions of galaxy
properties at z  10 from the simulation results (e.g. Behroozi et al.
© 2023 The Author(s)
Published by Oxford University Press on behalf of Royal Astronomical Society

Downloaded from https://academic.oup.com/mnras/article/525/4/4832/7251481 by Tsukuba Univ user on 02 November 2023

Accepted 2023 August 10. Received 2023 August 10; in original form 2022 November 22

FOREVER22: the first bright galaxies

2 METHODOLOGY
We use the results of our simulation project FOREVER22 (Yajima
et al. 2022) that focuses on protocluster regions in the cosmic volume
of (714 cMpc)3 . In this project, we use the GADGET-3 code (Springel
2005) with subgrid models developed in the OWLS project (Schaye
et al. 2010) and the FiBY project (Johnson et al. 2013). Besides,
we newly updated the code by adding the photoionization feedback,
the radiation pressure on dust, dust growth/destruction, black hole
growth, and its feedback (see more in Yajima et al. 2022). The
project consists of zoom-in simulations with three different levels
of the mass resolution and the size of zoom-in regions: protocluster
region [PCR; V = (28.6 cMpc)3 , smoothed particle hydrodynamics
(SPH) particle mass, mSPH = 4.1 × 106 M , and final redshift, zend =
2.0]; brightest protocluster galaxy [BCG; V ∼ (10 cMpc)3 , mSPH =
5.0 × 105 M , and zend = 4.0]; and First run [V ∼ (3 cMpc)3 , mSPH =
7.9 × 103 M , and zend = 9.5]. The PCR runs reproduce the observed
SFR densities of protoclusters at z ∼ 2–6. Also, we confirmed that the
mean density fields reproduced the observed stellar mass functions,
main sequences of star formation, gas fractions, and metallicities
of galaxies as a function of stellar mass well (Yajima et al. 2022).
In this work, we use First runs (First0 and First1 runs) in which
the most massive halo reaches Mh = 4.8 × 1011 M at z = 9.5.
The cosmological parameters are still under debate (Komatsu et al.

2011; Planck Collaboration VI 2020; Freedman 2021). Considering
the changing history of the parameter and Hubble parameter tension
(Freedman 2021), we adopt the cosmological parameters as M =
0.3, b = 0.045,  = 0.7, σ 8 = 0.82, and h = 0.7.
In this work, we consider both Pop II and Pop III stars. If the
metallicity is lower than a critical value, the initial mass function
(IMF) is likely to be a top-heavy (e.g. Chon, Omukai & Schneider
2021). Besides, the effective temperature of Pop III stars is high T ∼
105 K (Schaerer 2002). Therefore, Pop III stars can be strong sources
of radiative and SN feedback. We set the critical gas metallicity Z =
1.5 × 10−4 Z below which Pop III stars form (Omukai et al. 2005;
Frebel, Johnson & Bromm 2007; Chon et al. 2021). Although the
critical metallicity is still under debate, Abe et al. (2021) suggested
that the physical properties of the first galaxies did not depend on
it sensitively (see also Maio et al. 2010). We assume that the IMF
of Pop III stars is dn ∝ M−2.35 dM with the mass range 21–500 M ,
while that of Pop II is Chabrier IMF with the range 0.1–100 M .
Because of the expensive calculation costs for the first star formation
with radiative and magnetic feedback, the IMF of Pop III stars is still
under debate (Stacy & Bromm 2014; Susa, Hasegawa & Tominaga
2014; Hirano et al. 2015; Sugimura et al. 2020; Wollenberg et al.
2020; Latif, Whalen & Khochfar 2022b). Therefore, we adopt a
simple power-law function for the IMF of Pop III stars.
In evaluating the SFR, we consider the star formation model
based on the observed Kennicutt–Schmidt law that was developed
in Schaye & Dalla Vecchia (2008). The local SFR is measured as

−n  γ
(n−1)/2
˙ star ∝ mgas A 1 M pc−2
m
fP
, where mgas is the mass
G g
of a gas particle, γ = 5/3 is the ratio of specific heats, fg is the gas mass
fraction in the galactic disc, and P is the total interstellar medium
(ISM) pressure. Here, we set A = 1.5 × 10−4 M yr−1 kpc−2 and
γ = 1.4 for nH < 103 cm−3 and γ = 2.0 for nH ≥ 103 cm−3 . The
star formation model is the same as in EAGLE simulation project
(Schaye et al. 2015). Star formation occurs if local gas density

−0.64
exceeds nH = n0 cm−3 0.Z002
, where we set n0 = 10.0 for First
runs. In the estimation of the net cooling rate, we follow the nonequilibrium chemistry of primordial gas and the equilibrium state of
metals from pre-calculated tables with CLOUDY v07.02 code (Ferland
2000).
Once massive stars form, they give UV radiation feedback to
surrounding gas within their lifetime ∼107 yr. We take into account
the photoionization process of hydrogen and the dissociation of
hydrogen molecules. We estimate the volume of the ionized region
by taking the balance between the photon production rate and the
total recombination rate as (see the detail in Abe et al. 2021; Yajima
et al. 2022)
N˙ ion =

n

i=1

αB niH II nie

migas
i
ρgas

,

(1)

where N˙ ion is the photon production rate of a stellar particle, α B is
the case-B recombination coefficient, niH II and nie are the ionized
hydrogen and electron number densities of ith SPH particle. In
the ionized regions, the gas temperature is heated up to 3 × 104 K,
and star formation is prohibited. The dissociation rate of hydrogen
molecules is evaluated based on the contributions of stars in the
calculation box. First, we measure UV fluxes from stars with
distances to a target gas particle as

−2 

n

ri
m∗,i
JLW,21 =
fLW
(2)
,
1 kpc
103 M
i=1
where JLW, 21 is described in unit of 10−21 erg s−1 cm−2 Hz−1 sr−1 , ri
is the distance from ith stellar particle to a target gas particle and
MNRAS 525, 4832–4839 (2023)

Downloaded from https://academic.oup.com/mnras/article/525/4/4832/7251481 by Tsukuba Univ user on 02 November 2023

2020; Kannan et al. 2023). While these simulations allow us to study
the statistical natures of high-redshift galaxies with large cosmic
volumes, it is still difficult to study evolution from minihaloes hosting
the Population III (Pop III) stars to massive galaxies due to the limited
resolution.
In previous theoretical studies, galaxy formation at z  10
proceeds with the formation of Pop III stars, the radiative feedback,
and the metal enrichment via the first supernovae (SNe; e.g. Maio
et al. 2011; Wise et al. 2012; Johnson, Dalla & Khochfar 2013; Smith
et al. 2015; Xu et al. 2016; Chiaki & Wise 2019). Because of the metal
enrichment from Pop III stars, formation sites of new Pop III stars
move from higher to lower density regions in large-scale structure
(e.g. Tornatore, Ferrara & Schneider 2007; Pallottini et al. 2014;
Xu et al. 2016; Liu & Bromm 2020). Such numerical simulations
bridging from Pop III stars in minihaloes to first galaxies are still
challenging. Jeon & Bromm (2019) investigated the formation of
first galaxies with the halo mass of Mh ∼ 109 M at z = 9 and showed
their observational properties. Abe et al. (2021) studied the impact
of the initial mass function of Pop III stars on the physical properties
of first galaxies with Mh ∼ 108–9 M . They showed that inducing
frequent pair-instability SNe suppresses the gas mass fraction and
the SFRs of the first galaxies significantly for the top-heavy initial
mass function. The simulated halo masses in previous works have
been limited to ∼109 M . Therefore, the emergent UV fluxes were
too faint for the sensitivities of current telescopes.
Considering the brightness of observed candidates at z  10, they
can be hosted in massive haloes that likely form in overdense regions.
In this work, we investigate galaxy formation in overdense regions
in which the halo mass exceeds 1011 M at z ∼ 10. To study the
transition from Pop III to Population II (Pop II) stars, our simulations
resolve minihaloes and follow their growth up to the massive haloes.
Our paper is organized as follows. Section 2 shows our methodology and the information about the simulation set-up. In Section
3, we show the star formation histories and compare them with the
observational data by JWST. Also, we study the mass fraction of Pop
III stars with regard to the total stellar mass. Finally, we summarize
our results and discuss the limitations of our study in Section 4.

4833

4834

H. Yajima et al.

m∗, i is the mass of ith stellar particle. Then, we take into account
the self-shielding effect with the local H2 density and Jeans length
(Johnson et al. 2013):



1/2
f H 2  n H 1 / 2
T
NH2 = 2 × 1015 cm−2
, (3)
10−6
10 cm−3
103 K
where fH2 is the fraction of H2 and nH is the hydrogen number
density. We consider the shielding factor derived in Wolcott-Green,
Haiman & Bryan (2011) as
0.965
0.035
+
(1 + x/b5 )1.1
(1 + x)0.5



×exp −8.5 × 10−4 (1 + x)0.5 ,

(4)

where x ≡ NH2 /5 × 1014 cm−2 and b5 ≡ b/105 cm s−1 . Here b is the
Doppler broadening parameter, b ≡ (kB T/mH )1/2 . Thus, we estimate
the H2 dissociation rate (κ diss ) by combining JLW, 21 and fshield as
κ diss ∝ fshield JLW, 21 . In addition, we consider the photodetachment
process of H− (Shang, Bryan & Haiman 2010). With the dissociation
and formation rates, we evaluate H2 abundance and its radiative
cooling rate that is a main factor in controlling the formation of Pop
III stars in minihaloes.
When the age of a stellar particle reaches 107 yr, SN feedback turns
on. Following the SN feedback model in Dalla Vecchia & Schaye
(2012), we stochastically select a neighbouring gas particle and heat
the temperature up to 107.5 K. This hot bubble rapidly expands and
induces galactic wind, resulting in the suppression of star formation.

3 R E S U LT S
Fig. 1 shows the column density maps of gas and stars. Stellar
distributions are smoothed with a point spread function of JWST. The
gas widely distributes within virial radii, while stellar distributions
are compact and concentrated at the galactic centres. The gas accretes
onto galaxies along the filamentary structures and the stellar feedback
disturbs the gas structure. Stellar distributions and sizes change with
time. At z = 12, stellar clumps distribute at 0.5 Rvir , which reflects the
minor merger phase. As the galaxy grows via baryon accretion, the
size of the stellar components increases, but becomes small rapidly
when major mergers happen (see also Ono et al. 2023). Note that the
size shrinkage after the merger process sensitively depends on the
gas fraction and the structure of progenitor galaxies (e.g. Dekel &
Cox 2006).
Fig. 2 presents the redshift evolution of halo, stellar masses, and
SFR. The halo masses of the main progenitors are ∼109 M at z ∼
20 and evolve to ∼1011.5 M at z ∼ 10. The fluctuations are due to
mergers and the ability of the friends-of-friends (FoF) group finder
to identify all member particles. The rarity of the halo with 1011.5 M
at z = 10 is dN /d ln Mh ∼ 2 × 10−7 cMpc−3 . The cosmic volume
to host such a massive halo in our simulations is ∼500 cMpc3 that
is similar to the volumes of photometric galaxy surveys with JWST
(e.g. Finkelstein et al. 2023). Therefore, it can be reasonable to
directly compare our simulations with JWST data. Note that the
rarity changes with time even for the most massive progenitors
in the same region because of the variety of the halo merger
history.
As the halo grows, the stellar mass increases. The stellar masses
of the main progenitors exceed 108 M at z ∼ 13(14) and finally
reach 4.3(6.8) × 109 M in First0 (First1) run. The stellar masses
are similar to observed galaxies at z ࣠ 12. On the other hand,
it is much lower than the observed ones at z ∼ 16. Suppose the
estimated stellar masses of observed candidates are accurate and the
MNRAS 525, 4832–4839 (2023)

Downloaded from https://academic.oup.com/mnras/article/525/4/4832/7251481 by Tsukuba Univ user on 02 November 2023

fshield (NH2 , T ) =

redshifts are actually ∼16. In that case, most gas is very efficiently
converted into stars even in the early Universe (Harikane et al. 2023).
Inayoshi et al. (2022) suggested that 0.1–0.3 of the gas should be
converted into stars by using the abundance-matching technique
with the observed UV luminosity functions. In our simulations, the
SN feedback efficiently works in the suppression of star formation.
As a result, the gas is gradually converted into stars, and the star
formation efficiency (Mstar /Mgas ) of main progenitors in First0 run is
1.1 × 10−1 , 2.0 × 10−2 , and 2.5 × 10−2 at z = 10, 14, and 17. Note
that the redshifts and the physical properties of the candidate galaxies
at z > 14 were estimated with the photometric data. It is difficult to
evaluate the impacts of emission lines only from the photometric
data (e.g. Schaerer & de Barros 2009). Therefore, their properties
can be changed with follow-up spectroscopy (e.g. Arrabal Haro et al.
2023).
The SFR also increases with the growth of halo mass. Because of
the cycle of suppression of star formation and the short recovery timescale of gas, the SFR fluctuates significantly with time (Yajima et al.
2017). Main progenitor galaxies can have 10 M yr−1 at z ࣠ 14 and
show starbursts with SFR = 50.4 (133.5) M yr−1 at z = 9.5 in First0
Mgas
(First1) run. If we consider the SFR as SFR =  τdyn
, where  is an
efficiency parameter and τ dyn is the dynamical time, the starbursts at
z ∼ 9.5 correspond to  ∼ 0.08–0.25. This value is much larger than
typical star-forming galaxies in the local Universe. In the last period
of 0.1 Gyr (z = 9.5–11.5), the halo mass of the main progenitors in
First0 run increases from 1.2 × 1011 M to 4.8 × 1011 M . The rapid
mass growth with major mergers can induce the starburst. The SFRs
at z ࣠ 14 nicely match with the observed ones by JWST (Naidu et al.
2022; Donnan et al. 2023; Harikane et al. 2023) and GN-z11 at z =
10.957 (Jiang et al. 2021). Also, the modelled galaxies reproduce the
observed MUV at those redshifts. As suggested in Yajima et al. (2017),
once the halo mass exceeds 1011 M , most gas can be trapped in the
deep gravitational potential against SN feedback. This can induce
the starburst with SFR  10 M yr−1 and make galaxies observable.
In the redshift range, black holes are still in the state of being initial
seeds with ∼105 M and the accretion rates are mostly much lower
than the Eddington limit. Therefore, active galactic nucleus (AGN)
feedback is negligible.
Fig. 3 presents SFRs as a function of stellar mass. The SFR
increases with the stellar mass. At Mstar ࣠ 108 M , the suppression
of SFR due to the feedback makes the large dispersion. On the
other hand, massive galaxies keep the star formation continuous. We
confirm that our results match the observed galaxies. The stellar mass
monotonically increases with the halo mass, although there is a large
dispersion. The ratios of stellar to halo mass are ∼10−3 and ∼10−2
for Mh = 1010 and 1011 M at z = 10. These values are similar to the
empirical models in UNIVERSEMACHINE project (Behroozi et al.
2020), while it is somewhat higher than MILLENIUMTNG (Kannan
et al. 2023). The conversion efficiency sensitively depends on the
resolution, the star formation, and the feedback models. In particular,
our simulations can resolve minihaloes and dwarf galaxies, and their
star formation. At high redshifts, stars formed in dwarf galaxies can
contribute to the stellar mass in more massive galaxies via frequent
merger processes.
We present the relationships between the UV flux MUV and the halo
and the stellar mass in Fig. 4. We estimate MUV by measuring the
mean UV flux densities at λ = 1500–2000 Å in modelled spectral
energy distributions (SEDs) that will be shown in Fig. 6. MUV is
tightly related to the SFR, although it somewhat changes depending
on the star formation history. As the halo mass increases, the SFR
becomes higher, resulting in the formation of bright galaxies. We find
that the brightness can exceed the observable level MUV  −18 mag,

FOREVER22: the first bright galaxies

4835

if the halo mass is larger than ∼1010 M . In low-mass haloes,
galaxies at higher redshifts form stars more efficiently because
they are compact and have higher gas density typically. Also, there
is a large dispersion. This can be due to the SN feedback that
induces the intermittent star formation history via the cycle of gas
inflow and outflow (Yajima et al. 2017). On the other hand, the
UV brightness is more tightly correlated with the stellar mass. Our
simulations reproduce the observed UV brightnesses nicely. Galaxies
with Mstar  108 M are likely to have observable UV brightness
MUV  −18 mag. As shown in Fig. 2, the SFR rapidly increases
as the halo mass increases. Therefore, the stellar masses in the
massive haloes can be contributed mainly by the current starburst
episode, resulting in the tight relation in the massive systems. Note
that some observed galaxies with Mstar ∼ 106–7 M are brighter than
our modelled galaxies, although they are within the error bars. As
one possibility, hidden faint AGNs might contribute to observe UV
fluxes (e.g. Bunker et al. 2023). Future deep spectroscopic studies
will allow us to investigate AGN activities.
The upper panel of Fig. 5 shows the mass fraction of young Pop
III stars to the total stellar mass. As the star formation proceeds,
the interstellar gas is metal enriched via Type II SNe. Therefore, the
fraction steeply decreases as the stellar mass increases. Also, some
fractions of galaxies have no young Pop III stars. This indicates Pop
III stars form only when primordial gas clouds accrete on a galaxy.
Once the stellar mass exceeds ∼107 M , the fraction becomes ࣠0.01.
Considering the sensitivities of current telescopes, only massive
galaxies with Mstar  107 M have been observed. Therefore, Pop
II stars mainly form in the observed candidate galaxies at z  10.
We find that the mass fraction is insensitive to the redshift in the
range of z = 10–20. Given that the metal production source is only
Type II SNe, the total metal mass released is simply proportional to
the stellar mass. Thus, the insensitive redshift dependence indicates
similar metal mixing with the interstellar gas in the redshift range.

Our simulations suggest that low-mass galaxies with Mstar ࣠ 105 M
host Pop III stars with non-negligible fraction 0.1. Recently, Riaz,
Hartwig & Latif (2022) showed the mass fraction of Pop III stars
to the total stellar masses based on their seminumerical models
(Hartwig et al. 2022). The mass fraction of low-mass galaxies with
Mstar ∼ 104–5 M is similar to their results. On the other hand,
our results for massive galaxies are much higher. Our cosmological
simulations indicate that the gas in minihaloes can survive as the
primordial state and contribute to the Pop III star formation in massive
galaxies.
The lower panel of Fig. 5 represents the contribution of Pop III
stars to the UV flux at λ = 1500 Å in a rest frame. The contribution
fraction also decreases with the mass of Pop III stars. However, since
the mass-to-light ratio of Pop III stars is large and their effective
temperature is high T ∼ 105 K (Schaerer 2002), the contribution is
moderately large even if the mass fraction of Pop III stars is low. In
the cases of Mstar ∼ 107 M , it can be a few per cent. Therefore,
a part of the observed fluxes by JWST could be contributed by
Pop III stars. The UV flux can be dominated by Pop III stars if
the stellar mass is lower than ∼105 M . However, the low-mass
systems are likely to be too faint for the sensitivities of current
telescopes. Therefore, next-generation telescopes or gravitationally
lensed galaxies by foreground sources might be required for direct
observations of Pop III star clusters. Very recently, Vanzella et al.
(2023) indicated a candidate of a Pop III star cluster with the mass of
࣠104 M with the gravitational lens effect. The total metallicity even
for low-mass galaxies with Mstar ࣠ 105 M exceeds ∼10−3 Z/Z .
Therefore, the formation of Pop III stars indicates inhomogeneous
metal enrichment in a galaxy. The formation sites of Pop III stars are
somewhat far from the high-density regions of Pop II stars where the
metal enrichment proceeds earlier. Also, we find that the number
fraction of galaxies hosting young Pop III stars increases from
∼0.3 for Mstar = 105 M to ∼0.7 for Mstar = 108 M . The low-

MNRAS 525, 4832–4839 (2023)

Downloaded from https://academic.oup.com/mnras/article/525/4/4832/7251481 by Tsukuba Univ user on 02 November 2023

Figure 1. Column density maps of gas and stars in the most massive progenitors at z = 12, 14, and 17. The stellar distributions are smoothed with the point
spread function of JWST. Halo (stellar) masses are 5.3 × 109 (1.6 × 107 ) M at z = 17, 3.0 × 1010 (7.0 × 107 ) M at z = 14, and 1.0 × 1011 (4.1 × 108 ) M
at z = 12. Solid and dashed circles represent virial radii and half-virial radii, respectively. The size of the boxes is 3 arcsec, corresponding to L = 11.0, 9.8, and
8.4 physical kpc at z = 12, 14, and 17.

4836

H. Yajima et al.

Figure 4. UV fluxes as a function of halo and stellar masses. Filled circles
are our simulation results. Different colours represent the different redshifts.
Open symbols show the observational data: squares (Harikane et al. 2023),
circles (Naidu et al. 2022), and pentagons (Furtak et al. 2023).

Figure 3. SFRs as a function of stellar mass. Filled circles are our simulation
results. Different colours represent the different redshifts. Open symbols show
the observational data: squares (Harikane et al. 2023), circles (Naidu et al.
2022), and inverted triangles (Bunker et al. 2023).

MNRAS 525, 4832–4839 (2023)

mass galaxies without young Pop III stars consist of two states: star
forming only with Pop II stars and quenching of star formation due
to the SN feedback. Massive haloes are likely to distribute near the
centre of overdense regions and primordial gas hosted by minihaloes
can accrete them frequently. Therefore, the massive haloes can host
young Pop III stars although the mass fraction is low.
Note that even for the Pop III stars, we model the star formation by
replacing a gas particle with a stellar particle with the uniform mass
∼8 × 103 M , which models a star cluster. Therefore, the stellar
particles release the same SN energy and metal mass. However,
if the total stellar mass is smaller than ∼103 M in low-mass
haloes, the IMF may not be universal, resulting in unequal SN
feedback and metal amount (Abe et al. 2021). This can enhance
spatial fluctuation of metal distribution in the large cosmic volume.
Therefore, the metal distributions in low-mass galaxies are likely to
change with the resolution and the model of Pop III stars. We will

Downloaded from https://academic.oup.com/mnras/article/525/4/4832/7251481 by Tsukuba Univ user on 02 November 2023

Figure 2. Redshift evolution of the halo, stellar mass, and SFR. The red and
blue lines represent the properties of the most massive progenitors in First0
and First1 runs, respectively. Open symbols show the observations: inverted
triangles (Bunker et al. 2023), squares (Harikane et al. 2023), circles (Naidu
et al. 2022), and triangles (Donnan et al. 2023).

FOREVER22: the first bright galaxies

4837

Figure 5. Mass fraction of young Pop III stars to the total stellar mass (upper
panel) and their contribution fraction to UV flux at λ = 1500 Å in the rest
frame (lower panel). Different colours represent the different redshifts. The
galaxies with mass fractions lower than 10−4 are plotted at −4.0.

investigate these impacts on the fraction of Pop III stars in future
work.
Here, we derive intrinsic SEDs of galaxies by using a stellar
synthesis code STARBURST99 (Leitherer et al. 1999). Fig. 6 shows the
SEDs of the most massive progenitors at z = 12, 14, 17, and 20. In this
work, we estimate SEDs with the optically thin approximation, i.e.
no dust attenuation, which can be reasonable for low-mass and lowmetallicity galaxies (Yajima et al. 2012, 2014; Cullen et al. 2017).
The observed blue UV slope (β ࣠ −2) supports the assumption
(Naidu et al. 2022; Furtak et al. 2023). The contribution from
Pop III stars is estimated with the assumption of the brightness
temperature of 105 K and the mass-to-luminosity ratio for 120 M
derived in Schaerer (2002). The mass fractions of Pop III stars of
the galaxies are 1.1 × 10−3 , 1.8 × 10−3 , 4.7 × 10−3 , and 0.2 at
z = 12, 14, 17, and 20, respectively. Their contributions to UV
luminosity densities at λ = 1500 Å are 1.4 × 10−2 , 1.6 × 10−2 ,
7.5 × 10−2 , and 0.57. We find that the contribution fraction can be



Mstar,Pop III
+ 0.20. At z = 20, the
fit by log10 fUV,Pop III = 0.66 log10
Mstar
light from Pop III stars dominates at UV wavelengths. At the lower
redshifts, the UV-continuum fluxes are dominated by Pop II stars.
We measure the UV slopes of the SEDs by using the flux densities
at λ = 1300, 2000, and 3000 Å. It shows β = −2.68, −2.76, −2.62,
and −3.33 at z = 12, 14, 17, and 20, respectively. These naturally
reproduce the very blue slopes of the observed galaxies (Atek et al.
2023; Furtak et al. 2023). Note that we do not consider the nebular
emission in the above SEDs. If the nebular continuum at the UV
wavelengths is added into the SEDs, the slopes are changed to
β = −2.55, −2.60, −2.57, and −2.98 at z = 12, 14, 17, and 20.
Furthermore, dust extinction can make the SEDs redder at lower
redshifts.
The galaxies at z ≤ 14 have UV flux densities of 10 nJy that
are observable by JWST with a reasonable integration time. We
suggest that a part of the observed fluxes of candidate galaxies at
z  10 could be contributed by Pop III stars. If the sensitivity
of future observations at 3 μm will reach ∼0.1 nJy, the UV light
dominated by Pop III stars can be observed directly. Also, note that
Lyman-continuum (LyC) fluxes can be dominated by Pop III stars
significantly even if the mass fraction of Pop III stars is low, ࣠0.1.
For example, the contribution of Pop III stars to the LyC flux at
900 Å for the galaxy at z = 17 is 0.7. Thus, these galaxies consisting
of both Pop II and Pop III stars can be strong ionizing sources.
In addition, they may have unique properties in SEDs with high
equivalent widths of doubly ionized oxygen, carbon, and helium
(see also Nakajima & Maiolino 2022). In this work, we do not take
into account the radiative transfer in the galaxies. However, since
the galaxies are low-mass and low-metallicity systems, non-ionizing
UV-continuum photons are expected to escape efficiently. Therefore,
the estimated UV slopes are unlikely to change significantly. On
the other hand, the escape fraction of LyC photons can change
with time depending on the inhomogeneous gas structure due to
the SN feedback (Yajima et al. 2014; Paardekooper, Khochfar &
Dalla Vecchia 2015; Kimm et al. 2016; Trebitsch et al. 2017). In
MNRAS 525, 4832–4839 (2023)

Downloaded from https://academic.oup.com/mnras/article/525/4/4832/7251481 by Tsukuba Univ user on 02 November 2023

Figure 6. SEDs of the most massive progenitors at z = 12, 14, 17, and 20.
Blue dashed and orange solid lines show the radiation only from Pop III and
Pop II stars, respectively. Black solid lines are SEDs considering both stellar
populations.

4838

H. Yajima et al.

the case of a high escape fraction, He II and other metal lines can
be faint. As shown in Ono et al. (2023), the sizes of our modelled
galaxies change with time significantly. In phases when dusty gas
compactly distributes star-forming regions at the galactic centre, UV
photons can be attenuated even at such high redshifts. In practice,
the gas at the galactic centre reaches a metallicity with Z  0.1 Z
at z ∼ 10 (see also Isobe et al. 2023). The escape fraction of photons
sensitively depends on the covering fraction of dusty gas clouds from
young stars. We will perform radiative transfer simulations in future
work.

We have investigated the star formation and physical properties
in the first galaxies formed in overdense regions modelled by
the FOREVER22 simulation project. Our simulations followed the
evolution from minihaloes hosting Pop III stars to massive galaxies
with Mh > 1011 M . Our findings are summarized as follows.
(i) SFR increases with the halo mass and changes in the shorttime period due to the SN feedback. Once the halo mass exceeds
∼1011 M , galaxies continuously form stars with SFR  10 M yr−1
and induce starbursts with ∼100 M yr−1 . Even massive galaxies in
overdense regions cannot reproduce the observed stellar masses and
SFRs of candidate galaxies at z ∼ 16 suggested by Donnan et al.
(2023), Harikane et al. (2023), and Naidu et al. (2022).
(ii) Our simulations reproduce the relation between MUV and Mstar
of the observed galaxies at z  10 nicely. The galaxies with Mstar
 108 M show UV brightness of −18 mag that is observable by
JWST.
(iii) Even when the galaxy is metal enriched and forms Pop II stars,
Pop III stars can form in zero-metallicity spots. The mass fraction of
Pop III stars decreases as the stellar mass increases, and it is <0.01
for galaxies with Mstar ∼ 107 M . Therefore, candidate galaxies at
z  10 by JWST can be dominated by Pop II stars. We suggest that
a part of galaxies with Mstar ࣠ 105 M can host Pop III stars with a
non-negligible fraction 0.1.
(iv) We model SEDs of galaxies with Pop II and Pop III stars. The
UV-continuum fluxes of massive galaxies at z ≤ 17 are dominated by
Pop II stars. However, a few per cent of UV fluxes can be from Pop
III stars because of their large mass-to-luminosity ratio. The galaxies
at z ≤ 14 have the brightness of 10 nJy at λ = 2 μm thath can be
observable by JWST with a reasonable integration time.
The estimated physical properties and redshifts of the observed
galaxies at z  10 are not robust. Future spectroscopic observations
would present more reliable data and constrain the physical properties
of galaxies. On the other hand, the physical properties of first galaxies
modelled by numerical simulations can depend on the resolution and
the models of star formation and feedback (Abe et al. 2021). Also, the
seeding of the first massive black holes is still under debate (Inayoshi,
Visbal & Haiman 2020), and it may change the physical properties of
galaxies and SEDs. The star formation efficiency of massive haloes
with Mh  1010 M at z ∼ 10 is similar to the results in Behroozi
et al. (2020) and somewhat higher than Kannan et al. (2023). In this
paper, we newly provide insights about the relationship between Pop
III stars and massive galaxies at z  10. The physical properties are
likely to sensitively depend on the resolution, the star formation, and
the feedback models. We will investigate the model and resolution
dependencies on the first galaxy formation in future work.
MNRAS 525, 4832–4839 (2023)

The numerical simulations were performed on the computer cluster,
XC50 in NAOJ, and Trinity at Center for Computational Sciences
in University of Tsukuba. This work is supported in part by
MEXT/JSPS KAKENHI Grant Numbers 17H04827, 20H04724,
21H04489 (HY), 20K22358, 22H01258 (TH), NAOJ ALMA Scientific Research Grant Number 2019–11A, JST FOREST Program,
Grant Number JP-MJFR202Z, and Astro Biology Center Project
research AB041008 (HY). For the purpose of open access, the author
has applied a Creative Commons Attribution (CC BY) licence to any
author accepted manuscript version arising from this submission.
DATA AVA I L A B I L I T Y
The data underlying this paper will be shared on reasonable request
to the corresponding author.
REFERENCES
Abe M., Yajima H., Khochfar S., Dalla Vecchia C., Omukai K., 2021,
MNRAS, 508, 3226
Agarwal B., Dalla Vecchia C., Johnson J. L., Khochfar S., Paardekooper J.-P.,
2014, MNRAS, 443, 648
Arata S., Yajima H., Nagamine K., Li Y., Khochfar S., 2019, MNRAS, 488,
2629
Arrabal Haro P. et al., 2023, preprint (arXiv:2303.15431)
Atek H. et al., 2023, MNRAS, 519, 1201
Bakx T. J. L. C. et al., 2023, MNRAS, 519, 5076
Behroozi P. et al., 2020, MNRAS, 499, 5702
Behroozi P., Wechsler R. H., Hearin A. P., Conroy C., 2019, MNRAS, 488,
3143
Bouwens R. J., Stefanon M., Oesch P. A., Illingworth G. D., Nanayakkara T.,
Roberts-Borsani G., Labbe´ I., Smit R., 2019, ApJ, 880, 25
Bunker A. J. et al., 2023, preprint (arXiv:2302.07256)
Capak P. L. et al., 2015, Nature, 522, 455
Chiaki G., Wise J. H., 2019, MNRAS, 482, 3933
Chon S., Omukai K., Schneider R., 2021, MNRAS, 508, 4175
Cullen F. et al., 2023, MNRAS, 520, 14
Cullen F., McLure R. J., Khochfar S., Dunlop J. S., Dalla Vecchia C., 2017,
MNRAS, 470, 3006
Dalla Vecchia C., Schaye J., 2012, MNRAS, 426, 140
Dekel A., Cox T. J., 2006, MNRAS, 370, 1445
Donnan C. T. et al., 2023, MNRAS, 518, 6011
Ferland G. J., 2000, Rev. Mex. Astron. Astrofis. Ser. Conf., 9, 153
Finkelstein S. L. et al., 2013, Nature, 502, 524
Finkelstein S. L. et al., 2023, ApJ, 946, L13
Frebel A., Johnson J. L., Bromm V., 2007, MNRAS, 380, L40
Freedman W. L., 2021, ApJ, 919, 16
Furtak L. J., Shuntov M., Atek H., Zitrin A., Richard J., Lehnert M. D.,
Chevallard J., 2023, MNRAS, 519, 3064
Gabrielpillai A., Somerville R. S., Genel S., Rodriguez-Gomez V., Pandya
V., Yung L. Y. A., Hernquist L., 2022, MNRAS, 517, 6091
Harikane Y. et al., 2023, ApJS, 265, 5
Hartwig T. et al., 2022, ApJ, 936, 45
Hashimoto T. et al., 2018, Nature, 557, 392
Hirano S., Hosokawa T., Yoshida N., Omukai K., Yorke H. W., 2015,
MNRAS, 448, 568
Inayoshi K., Harikane Y., Inoue A. K., Li W., Ho L. C., 2022, ApJ, 938, L10
Inayoshi K., Visbal E., Haiman Z., 2020, ARA&A, 58, 27
Inoue A. K. et al., 2016, Science, 352, 1559
Isobe Y. et al., 2023, preprint (arXiv:2307.00710)
Jeon M., Bromm V., 2019, MNRAS, 485, 5939
Jiang L. et al., 2021, Nat. Astron., 5, 256
Johnson J. L., Dalla V. C., Khochfar S., 2013, MNRAS, 428, 1857
Kannan R. et al., 2023, MNRAS, 524, 2594

Downloaded from https://academic.oup.com/mnras/article/525/4/4832/7251481 by Tsukuba Univ user on 02 November 2023

4 D I S C U S S I O N A N D S U M M A RY

AC K N OW L E D G E M E N T S

FOREVER22: the first bright galaxies

Shibuya T., Kashikawa N., Ota K., Iye M., Ouchi M., Furusawa H., Shimasaku
K., Hattori T., 2012, ApJ, 752, 114
Smith B. D., Wise J. H., O’Shea B. W., Norman M. L., Khochfar S., 2015,
MNRAS, 452, 2822
Song M., Finkelstein S. L., Livermore R. C., Capak P. L., Dickinson M.,
Fontana A., 2016, ApJ, 826, 113
Springel V., 2005, MNRAS, 364, 1105
Stacy A., Bromm V., 2014, ApJ, 785, 73
Sugimura K., Matsumoto T., Hosokawa T., Hirano S., Omukai K., 2020, ApJ,
892, L14
Susa H., Hasegawa K., Tominaga N., 2014, ApJ, 792, 32
Tamura Y. et al., 2019, ApJ, 874, 27
Topping M. W., Stark D. P., Endsley R., Plat A., Whitler L., Chen Z., Charlot
S., 2022, ApJ, 941, 153
Tornatore L., Ferrara A., Schneider R., 2007, MNRAS, 382, 945
Trebitsch M., Blaizot J., Rosdahl J., Devriendt J., Slyz A., 2017, MNRAS,
470, 224
Vanzella E. et al., 2023, preprint (arXiv:2305.14413)
Wise J. H., Demchenko V. G., Halicek M. T., Norman M. L., Turk M. J., Abel
T., Smith B. D., 2014, MNRAS, 442, 2560
Wise J. H., Regan J. A., O’Shea B. W., Norman M. L., Downes T. P., Xu H.,
2019, Nature, 566, 85
Wise J. H., Turk M. J., Norman M. L., Abel T., 2012, ApJ, 745, 50
Wolcott-Green J., Haiman Z., Bryan G. L., 2011, MNRAS, 418, 838
Wollenberg K. M. J., Glover S. C. O., Clark P. C., Klessen R. S., 2020,
MNRAS, 494, 1871
Xu H., Norman M. L., O’Shea B. W., Wise J. H., 2016, ApJ, 823, 140
Yajima H. et al., 2022, MNRAS, 509, 4037
Yajima H., Choi J.-H., Nagamine K., 2011, MNRAS, 412, 411
Yajima H., Khochfar S., 2016, MNRAS, 457, 2423
Yajima H., Li Y., Zhu Q., Abel T., 2012, MNRAS, 424, 884
Yajima H., Li Y., Zhu Q., Abel T., Gronwall C., Ciardullo R., 2014, MNRAS,
440, 776
Yajima H., Nagamine K., Zhu Q., Khochfar S., Dalla Vecchia C., 2017, ApJ,
846, 30
Yajima H., Sugimura K., Hasegawa K., 2018, MNRAS, 477, 5406
Yajima H., Umemura M., Mori M., Nakamoto T., 2009, MNRAS, 398, 715
Yoon I. et al., 2023, ApJ, 950, 61

This paper has been typeset from a TEX/LATEX file prepared by the author.

© 2023 The Author(s)
Published by Oxford University Press on behalf of Royal Astronomical Society

MNRAS 525, 4832–4839 (2023)

Downloaded from https://academic.oup.com/mnras/article/525/4/4832/7251481 by Tsukuba Univ user on 02 November 2023

Kannan R., Garaldi E., Smith A., Pakmor R., Springel V., Vogelsberger M.,
Hernquist L., 2022, MNRAS, 511, 4005
Kimm T., Cen R., Rosdahl J., Yi S., 2016, ApJ, 823, 52
Komatsu E. et al., 2011, ApJS, 192, 18
Latif M. A., Whalen D. J., Khochfar S., Herrington N. P., Woods T. E., 2022a,
Nature, 607, 48
Latif M. A., Whalen D., Khochfar S., 2022b, ApJ, 925, 28
Leitherer C. et al., 1999, ApJS, 123, 3
Liu B., Bromm V., 2020, MNRAS, 497, 2839
Lovell C. C., Vijayan A. P., Thomas P. A., Wilkins S. M., Barnes D. J.,
Irodotou D., Roper W., 2021, MNRAS, 500, 2127
Ma X., Quataert E., Wetzel A., Hopkins P. F., Faucher-Gigue` re C.-A., Keresˇ
D., 2020, MNRAS, 498, 2001
Maio U., Ciardi B., Dolag K., Tornatore L., Khochfar S., 2010, MNRAS,
407, 1003
Maio U., Khochfar S., Johnson J. L., Ciardi B., 2011, MNRAS, 414,
1145
Naidu R. P. et al., 2022, ApJ, 940, L14
Nakajima K., Maiolino R., 2022, MNRAS, 513, 5134
Ocvirk P. et al., 2016, MNRAS, 463, 1462
Oesch P. A. et al., 2013, ApJ, 773, 75
Oesch P. A. et al., 2016, ApJ, 819, 129
Omukai K., Tsuribe T., Schneider R., Ferrara A., 2005, ApJ, 626, 627
Ono Y. et al., 2012, ApJ, 744, 83
Ono Y. et al., 2023, ApJ, 951, 72
Ouchi M. et al., 2018, PASJ, 70, S13
Paardekooper J.-P., Khochfar S., Dalla Vecchia C., 2013, MNRAS, 429, L94
Paardekooper J.-P., Khochfar S., Dalla Vecchia C., 2015, MNRAS, 451, 2544
Pakmor R. et al., 2023, MNRAS, 524, 2539
Pallottini A., Ferrara A., Gallerani S., Salvadori S., D’Odorico V., 2014,
MNRAS, 440, 2498
Planck Collaboration VI, 2020, A&A, 641, A6
Popping G., 2023, A&A, 669, L8
Regan J. A., Haehnelt M. G., 2009, MNRAS, 396, 343
Riaz S., Hartwig T., Latif M. A., 2022, ApJ, 937, L6
Rosdahl J. et al., 2022, MNRAS, 515, 2386
Schaerer D., 2002, A&A, 382, 28
Schaerer D., de Barros S., 2009, A&A, 502, 423
Schaye J. et al., 2010, MNRAS, 402, 1536
Schaye J. et al., 2015, MNRAS, 446, 521
Schaye J., Dalla Vecchia C., 2008, MNRAS, 383, 1210
Shang C., Bryan G. L., Haiman Z., 2010, MNRAS, 402, 1249

4839