# Formation and Morphology of the First Galaxies in the Cosmic Morning

Changbom Park, Jaehyun Lee, Juhan Kim, Donghui Jeong, Christophe Pichon, Brad K. Gibson, Owain N. Snaith, Jihye Shin, Yonghwi Kim, Yohan Dubois, et al.

► To cite this version:

Changbom Park, Jaehyun Lee, Juhan Kim, Donghui Jeong, Christophe Pichon, et al. Formation and Morphology of the First Galaxies in the Cosmic Morning. The Astrophysical Journal, 2022, 937, 10.3847/1538-4357/ac85b5. insu-03839648

HAL Id: insu-03839648

https://insu.hal.science/insu-03839648

Submitted on 4 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Formation and Morphology of the First Galaxies in the Cosmic Morning

Changbom Park<fnref n="1" />, Jaehyun Lee<fnref n="1" />, Juhan Kim<fnref n="1" />, Donghui Jeong<fnref n="1" />,<fnref n="2" />, Christophe Pichon<fnref n="1" />,<fnref n="3" />,<fnref n="4" />, Brad K. Gibson<fnref n="5" />, Owain N. Snaith<fnref n="6" />, Jihye Shin<fnref n="7" />, Yonghwi Kim<fnref n="1" />, Yohan Dubois<fnref n="3" />, and C. Gareth Few<fnref n="5" />

*Received 2022 February 24; revised 2022 July 22; accepted 2022 July 27; published 2022 September 16*

## Abstract

We investigate the formation and morphological evolution of the first galaxies in the cosmic morning ($10 \gtrsim z \gtrsim 4$) using the Horizon Run 5 (HR5) simulation. For galaxies above the stellar mass $M_{\star,\min} = 2 \times 10^9 M_{\odot}$, we classify them into disk, spheroid, and irregular types according to their asymmetry and stellar-mass morphology. We find that about two-thirds of the galaxies have a Sérsic index <1.5, reflecting the dominance of disk-type morphology in the cosmic morning. The rest are evenly distributed as incidental and transient irregulars or spheroids. These fractions are roughly independent of redshift and stellar mass up to $\sim 10^{10} M_{\odot}$. Almost all the first galaxies with $M_{\star} > M_{\star,\min}$ at $z > 4$ form at initial peaks of the matter-density field. Large-scale structures in the universe emerge and grow like cosmic rhizomes as the underlying matter-density fluctuations grow and form associations of galaxies in rare overdense regions and the realm of the galactic world is stretched into relatively lower-density regions along evolving filaments. The cosmic web of galaxies forms at lower redshifts when most rhizomes globally percolate. The primordial angular momentum produced by the induced tidal torques on protogalactic regions is correlated with the internal kinematics of galaxies and tightly aligned with the angular momentum of the total galaxy mass. The large-scale tidal field imprinted in the initial conditions seems responsible for the dominance of disk morphology and for the tendency of galaxies to reacquire a disk postdistortion.

*Unified Astronomy Thesaurus concepts: Galaxy evolution (594); Galaxy formation (595)*

## 1. Introduction

Which morphology do galaxies have when they form: Do the first galaxies have a single morphology or an initial morphological distribution? How do the initial conditions govern the formation and morphology of the first galaxies? How do subsequent evolutionary processes affect galaxy morphology? Galaxy morphology is one of the key physical properties of galaxies, but these questions regarding the galaxy morphology have not yet been fully answered. At the present epoch, the morphology of most bright galaxies can be categorized according to Hubble's tuning fork (Hubble 1926). In this system, galaxies are classified into irregular and regular types, which are in turn divided into ellipticals, spirals, and lenticulars.

Using a volume-limited sample of 1332 GOODS galaxies with rest-frame $M_B < -18$ and $0.4 \le z \le 1.0$, Hwang & Park (2009) have found that the fraction of early-type galaxies (ellipticals plus lenticulars) has been lower in the past, reaching as low as about 15% near $z=1$. They have also found the morphology-density relation (Oemler 1974; Dressler 1980; Park et al. 2007) observed at low redshifts almost disappears near $z=1$. Namely, near $z=1$, the early-type fraction stays at $\sim 15\%$ regardless of the density environment. Similarly, Tamburri et al. (2014) used 1302 GOODS-MUSIC galaxies brighter than $K_s(AB) = 22$ and with stellar mass roughly between $10^{10}$ and $3 \times 10^{11} M_{\odot}$ to find that at redshift $z \approx 1$, late

types are dominant and the early-type galaxy fraction is about 0.19. A natural question raised from these findings is whether or not galaxies at much higher redshifts ($z \gg 1$) are all late-type galaxies such as disks and irregulars.

There have been many observational efforts to study galaxy morphology at higher redshifts (see Elmegreen 2015; Marshall et al. 2019, for reviews). However, a consensus is yet to be reached. Some observations support the claim that high-redshift galaxies are more clumsy and irregular than low-redshift galaxies (Abraham & van den Bergh 2001; Papovich et al. 2005; Elmegreen & Elmegreen 2014). Cameron et al. (2011) have reported using the Hubble Ultra Deep Field data and the GOODS South data that irregulars are the most prominent type at $z > 2$. On the other hand, many reports do not support the above claims and find a larger fraction of spheroidals at high redshifts (Bundy et al. 2005; Franceschini et al. 2006; Lotz et al. 2006; Dahlen et al. 2007).

Ravindranath et al. (2006) have reported using GOODS images that, at $z \sim 3$, 40% of Lyman-break galaxies have exponential profiles, 30% have de Vaucouleurs profiles, and about 30% have disturbed morphology. For massive galaxies with $M_{\star} > 10^{11} M_{\odot}$ Buitrago et al. (2013) have found using a sample of HST galaxies that the fraction of peculiar-type galaxies is about the same as that of early- and late-type galaxies at $z \approx 2$. The recently launched James Webb Space Telescope will keep a close eye on this high-redshift universe and reveal the morphology of high-redshift galaxies.

Cosmological simulations have been invaluable tools to understand the formation and evolution of cosmic structures (e.g., Park 1990; Ocvirk et al. 2008; Park et al. 2012;

<figure><img src="image_2.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Korea Institute for Advanced Study, 85 Hoegi-ro, Dongdaemun-gu, Seoul 02455, Republic of Korea; syncphy@gmail.com</footnote>

<footnote n="2">² The Pennsylvania State University, University Park, PA, 16802, USA</footnote>

<footnote n="3">³ CNRS and Sorbonne Université, UMR 7095, Institut d'Astrophysique de Paris, 98 bis, Boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="4">⁴ IPhT, DRF-INP, UMR 3680, CEA, L'Orme des Merisiers, Bât 774, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="5">⁵ E.A. Milne Centre for Astrophysics, University of Hull, Hull, HU6 7RX, UK</footnote>

<footnote n="6">⁶ GEPI, Observatoire de Paris, PSL Research University, CNRS, 5 Place Jules Janssen, F-92190, Meudon, France</footnote>

<footnote n="7">⁷ Korea Astronomy and Space Science Institute, 776 Daedeokdae-ro, Yuseong-gu, Daejeon 34055, Republic of Korea</footnote>
Dubois et al. 2014b; Genel et al. 2014; Khandai et al. 2015; Schaye et al. 2015; Dolag et al. 2016; Tremmel et al. 2017; McCarthy et al. 2017; Marinacci et al. 2018; Davé et al. 2019; Lee et al. 2021). In recent years, several studies using large cosmological simulations accounting for gravity, hydrodynamics, and various subgrid physics have addressed galaxy morphology at high redshifts and the origin of Hubble's tuning fork classification system.

Dekel et al. (2009) have highlighted the specific importance of cold streams as the main mode of gas acquisition at high redshifts, while Pichon et al. (2011) and Stewart et al. (2013) have shown how steady angular momentum advected along these streams are key to the disk re-formation. The consistency of vorticity-rich inflow (Laigle et al. 2015) is due to the coherency of tides on larger scales (Danovich et al. 2015), which are themselves a result of the structure of the initial conditions (Bond et al. 1996). Because cold gas undergoes shocks during shell crossing (Kereš et al. 2005), it closely follows the ridges in the large-scale structure.

From the Horizon simulation, Dubois et al. (2016) highlights the interplay between cosmic gas accretion onto galaxies and galaxy mergers in driving the observed diversity in galaxy morphology at lower redshifts (see also Martin et al. 2018). They have noted the critical role of active galactic nuclei (AGNs) in setting up the correct galaxy morphology for the massive end of the population. AGNs accomplish this by freezing the morphological type of the postmerger remnant and quenching its quiescent star formation for a considerable duration.

Pillepich et al. (2019) have used the Illustris TNG50 simulations to focus on the evolution of star-forming galaxies and conclude that the dominant galaxy morphology at $z > 2$ is spheroidal for galaxies with $M_* < 10^{10} M_\odot$. They also claim that the fraction of disk-like galaxies increases at lower redshifts and for larger stellar mass. According to their results, disk-type galaxies become dominant at $z < 2$ for actively star-forming galaxies with $M_* > 10^{9.5} M_\odot$.

Trayford et al. (2019) have analyzed a sample of simulated galaxies with $M_* > 10^9 M_\odot$ in the EAGLE simulation. They assign morphology to galaxies using the fraction of the stellar mass that lies in the disk component and claimed that, at redshifts $z \gtrsim 1.5$, asymmetric galaxies dominate for all masses. They also reported that the Hubble sequence is established by $z \simeq 1.5$, and the spheroid type dominates after that epoch. The disk type is always subdominant in this work. It should be pointed out that their morphology classification of disk and spheroid types is based not on the shape of mass distribution but actually on kinematics parameterized by the circularity $\epsilon$.

Park et al. (2019) used the New Horizon simulation to find that thin disks in galaxies start to form only at $z \sim 1-2$ for those with $M_*(z=0.7) > 10^{10} M_\odot$ or at $z \approx 1$ for $M_*(z=0.7) = 10^9-10^{10} M_\odot$. The former era matches well with some measurements of the history of the Milky Way (e.g., Haywood et al. 2013).

Hence, neither observations nor simulations give a clear and consistent picture of the initial morphology of high-redshift galaxies. In this study, we use the recent Horizon Run 5 cosmological simulation (Lee et al. 2021), HR5, to inspect galaxies at high redshifts $10 \gtrsim z \gtrsim 4$, which we call cosmic morning.<fnref n="8" /> We will study the physical properties of the first

galaxies and large-scale structures in the cosmic morning focusing on their morphology. Relying on the unique volume and resolution of HR5, we shall clarify the origin and early evolution of the morphology of galaxies, respectively, by studying the relationship between morphology and the initial conditions and by closely examining the cosmic evolution of physical parameters.

## 2. Simulation Data

The cosmological simulation HR5 (Lee et al. 2021) covers a volume of the universe in an 1.15 cGpc³ scale volume while achieving a resolution down to ~1 pkpc. We adopt the cosmological parameters of $\Omega_m = 0.3$, $\Omega_\Lambda = 0.7$, $\Omega_b = 0.047$, $\sigma_8 = 0.816$, and $h = 0.684$ that are compatible with the Planck data (Planck Collaboration et al. 2016). The linear power spectrum is computed using the CAMB package (Lewis et al. 2000) and the initial conditions are generated at $z = 200$ using the MUSIC package (Hahn & Abel 2011) with second-order Lagrangian perturbation theory (2LPT; Scoccimarro 1998; L'Huillier et al. 2014).

We use a modified version of the adaptive mesh refinement code RAMSES (Teyssier 2002) for HR5. The Poisson equations and Euler equations are solved based on the Particle-Mesh method (Guillet & Teyssier 2011) and the Harten–Lax–van Leer contact wave Riemann solver (Toro et al. 1994), respectively. All particles of dark matter, stars, and massive black holes (MBHs) are assigned into grids using a cloud-in-cell scheme to construct mass-density grids. We have implemented OpenMP parallelism into RAMSES on top of the Message Passing Interface originally equipped in RAMSES.

We set a high-resolution cuboid region of $1048.6 \times 119.0 \times 127.2$ cMpc³ at the center of the simulation box of $1048.6 \times 10^7$ cMpc³. The cosmological simulation box is set to have a background grid level 8, which corresponds to 256 coarse grids on a side with $\Delta x = 4.10$ cMpc. The high-resolution region is initially filled with $8192 \times 930 \times 994$ grids and dark-matter particles at a level 13, and the padding grids of levels from 9 to 12 are placed between the coarse- and high-resolution regions. A dark-matter particle has a mass resolution of $6.89 \times 10^7 M_\odot$ at the highest level, and the particle mass increases by a factor of 8 when the grid level decreases by 1. The grids are adaptively refined down to $\Delta x \sim 1$ pkpc when their mass is larger than eight times the minimum dark-matter particle mass. We have carried out the HR5 simulation down to $z = 0.625$.

The version of RAMSES used in HR5 computes gas cooling down to a temperature of $10^4$ K using the synthetic cooling functions of Sutherland & Dopita (1993), and metal-enriched gas can additionally cool down to ~750 K based on the cooling rates proposed by Dalgarno & McCray (1972). Cosmic reionization is approximated by adopting the uniform UV background of Haardt & Madau (1996). Star formation is computed using the statistical approach of Rasera & Teyssier (2006). The minimum mass of stellar particles is $2.56 \times 10^6 M_\odot$ at birth.

Supernova feedback operates in thermal and kinetic modes (Dubois & Teyssier 2008) and AGN feedback releases its energy in radio-jet and quasar modes, depending on the Eddington ratio (Dubois et al. 2012). Massive black holes (MHBs) are seeded with an initial mass of $10^4 M_\odot$ in grids when gas density is higher than the threshold of star formation and no other MBHs exist within 50 kpc (Dubois et al. 2014b). MBHs grow via accretion and coalescence, and their angular momentum obtained

<footnote n="8">⁸ We shall call this period “the cosmic morning,” reserving the name “cosmic dawn” for the period before $z \simeq 10$. Similarly, the cosmic noon is from $z \simeq 4$ to 1.5, and the cosmic afternoon is the period of $z \lesssim 1.5$.</footnote>
from the feeding processes is traced (Dubois et al. 2014a). We model the chemical evolution by adopting the method proposed by Few et al. (2012), and trace the abundance of H, O, and Fe based on a Chabrier initial mass function (IMF) (Chabrier 2003).

We calibrate the subgrid physics of HR5 to reproduce the star formation rate density, the galaxy stellar-mass function, and the mass–metallicity relation using test simulations. The performance of the parameter tuning was assessed by eyes with no formal fitting because of limited computing resources. The details of the calibration are given in Appendix A of Lee et al. (2021).

We compute the rest-frame luminosities of stellar particles in the Johnson and SDSS filter systems based on their ages and metallicities during postprocessing. We use the photometric predictions generated from E-MILES single stellar population SEDs (Vazdekis et al. 2012, 2016; Ricciardelli et al. 2012) in this modeling. For consistency with our chemical evolution modeling, we also assume a Chabrier IMF in this calculation. One can find further details of HR5 in Lee et al. (2021).

## 2.1. Galaxy Finding

RAMSES employs four types of astrophysical objects to study their mutual interactions in accordance with gravity and hydrodynamics. Dark matter, stars, and MBHs are represented by point-mass particles while information on gas is assigned to finite-volume cells. The particle-like representation has long been favored in most cosmological N-body simulations, and various structure-identification methods have successfully been developed using well-known particle-based algorithms, such as the friends of friends (FoF). On the other hand, a gas cell is stationary and has gas quantities being averaged over its volume. In the pipeline of structure finding, we assume the gas cells to be particles, like other matter components. Because the particles have varying mass, the linking length of a particle pair is derived by averaging the linking length of each particle.

To find substructures (or galaxies) in the HR5, we modified the original PSB subhalo finder (Kim & Park 2006), which was developed for N-body simulations (Kim et al. 2015, 2009). Like the PSB method, PGalf (Lee et al. 2021) is based on the coordinate-free density map and a network of neighboring particles. However, PGalf uses the stellar density map if there are stellar particles because stellar distributions are clumpier and more compact, making it easier to determine the boundaries of substructures (or galaxies). To identify star particle groups, we set the minimum number of star particles to be 10, which corresponds to ~2 × 10⁷ $M_\odot$ in the cosmic morning. When HR5 is calibrated, statistical comparison with observations is made for the “galaxies” with stellar mass above 1 × 10⁹ $M_\odot$, and we will use those having $M_* > 2 × 10^9 M_\odot$ in most parts of our galaxy morphology study. Therefore, galaxies in this study have more than 1,000 star particles.

An unbinding scheme is essential for subhalo identification to minimize contamination by host halo particles (Knebe et al. 2011). Thus, various unbinding algorithms have been introduced to halo finders, e.g., the local escape velocity scheme of AHF (Knollmann & Knebe 2009), the density gradient of SUBFIND (Springel et al. 2001), or the neighborhood link in six-dimensional configuration space of Rockstar (Behroozi et al. 2013) or VELOCIRaptor (Elahi et al. 2019). Following the scheme of Kim & Park (2006), we utilize two criteria to determine the membership of a particle with respect to a substructure. The first is the tidal radius limit and the other one

<table><caption>**Table 1**
Number of HR5 Galaxies above Given Stellar-mass Limits</caption><thead><tr><th>Redshift</th><th>$M_*$ > 2 × 10<sup>9</sup> $M_\odot$</th><th>$M_*$ > 5 × 10<sup>9</sup> $M_\odot$</th><th>Volume (10<sup>7</sup> cMpc<sup>3</sup>)</th></tr></thead><tbody><tr><td>9</td><td>14 (15)</td><td>– (–)</td><td>1.3508</td></tr><tr><td>8</td><td>71 (83)</td><td>5 (6)</td><td>1.3490</td></tr><tr><td>7</td><td>495 (593)</td><td>58 (67)</td><td>1.3459</td></tr><tr><td>6</td><td>2205 (2669)</td><td>1372 (446)</td><td>1.3425</td></tr><tr><td>5</td><td>8534 (10446)</td><td>1918 (2348)</td><td>1.3375</td></tr><tr><td>4</td><td>22136 (27441)</td><td>6095 (7590)</td><td>1.3317</td></tr></tbody></table>

**Note.** The numbers outside parentheses in the second and third columns are the counts of the galaxies located at least 3 cMpc away from low-level particles (the clean set) while those inside are for the entire set. The fourth column presents the volume uncontaminated by massive/low-level particles. In this study, we use the clean set only.

is the negative total energy criterion. We have used the tidal radius ($r_t$),

$$r_t = R \left( \frac{m/M}{\alpha + \beta} \right), \qquad (1)$$

where R is the distance between the host and satellite galaxies, M is the host galaxy mass contained within radius R,

$$\alpha(R) \equiv 2 - \left( \frac{d \ln M}{d \ln R} \right) |_{R}, \qquad (2)$$

and

$$\beta \equiv \frac{\Omega^2 R^3}{GM}, \qquad (3)$$

where Ω is the orbital angular speed of the satellite. In these equations, α is the contribution of the gravitational force by the host and β is the effect of centrifugal force during the orbital motion of the satellite.

A particle is first identified to be bound to a structure when it has negative total energy relative to the structure. The membership is probed if the particle is located inside the tidal radius of the structure. Therefore, even though a particle has negative total energy with respect to a structure, the particle does not belong to the structure if it is located outside the structure’s tidal radius. When a particle satisfies the membership conditions for multiple structures, we assign the particle to the least-massive galaxies. The tidal boundary scheme is implemented in PGalf to separate close pairs of galaxies strongly interacting with each other. In a compact group of galaxies, for example, galaxies are likely to be strongly tied to each other with complex tidal structures, and the tidal boundary condition may help us to draw border lines between contacting galaxies.

Table 1 shows the numbers of galaxies above a given stellar mass identified in the HR5 snapshots at $z=9-4$. Because stars form in high-density gas cells located only in the zoomed region, galaxies are also identified in the region filled with grids of level 13 or higher. However, low-level/massive particles initially placed outside the zoomed region gradually permeate into the high-resolution region due to the evolution of the density field as a function of cosmic time. In this study, we flag grids that contain low-level DM particles ($M_{DM} \ge 5.51 \times 10^8 M_\odot$) and only use galaxies at least 3 cMpc away from the contaminated regions to minimize the potential boundary
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Schematic diagram for the varying branches of merger trees constructed between snapshots *n* − 1, *n*, and *n* + 1. The gray shades are self-bound structures, the diamonds present the MBPs in the self-bound structures, and the colored ellipses are galaxies in them. The black solid arrows indicate main progenitor/descendant relations, dotted arrows mark the mass transfers between the structures with no main progenitor/descendant relations, and the dashed lines indicate mergers along the main descendant branches. The green solid arrows show the backtracking of the MBPs to repair the branches broken by structure misidentification, and the red solid arrow marks the main progenitor–descendant branch repaired by the MBP backtracking scheme.</figcaption></figure>

effects caused by massive particles. One can find the details of
the geometry evolution of the zoomed region in Appendix A of
Lee et al. (2021).

## 2.2. Merger Tree

We construct the merger trees of self-bound objects (SBOs)
throughout the entire snapshots of HR5 based on the most-
bound particles (MBPs) and stellar particles. The MBP is
defined as the dark-matter particle that has the least total
(kinetic+potential) energy in a SBO (Hong et al. 2016). We
assume that the evolutionary history of a SBO can be traced by
using the MBP, which, by definition, is most strongly bound to
the host SBO. Because dark-matter particles do not newly form
or disappear in HR5, we use dark-matter particles to identify
the MBPs of SBOs.

Figure 1 shows schematic diagrams of some examples of
various merger trees connecting SBOs in a series of three
snapshots. Case 1 in this figure illustrates how the MBP
scheme works for SBOs without stellar particles. The tree-
building algorithm first identifies the MBP of each SBO at an
initial time step and then finds the SBO that possesses the MBP
of interest at the next time step. We then assign the substructure
containing the MBP of the SBO at the next time step as the
main descendant of the SBO. An SBO can acquire multiple
MBPs from more than one SBO in the earlier snapshot. In that
case, we define the MBP donor with which the SBO shares its
largest dark-matter mass fraction as the main progenitor and
classify the rest of the MBP donors as merged ones. Thus, in

this algorithm, more than one SBOs at a given time step can
have a common descendant in the next snapshot. The multiple
MBPs coming from merged structures are tracked as the sub-
MBPs along the main branches of their host SBOs throughout
the entire snapshots. The main MBPs of SBOs are updated at
every time step because the internal kinematics of structures
continuously changes. The MBP scheme enables us to trace the
trajectories of SBOs that might have been missed in between
the time steps due to the structure mis-identifications that
frequently occur in dense environments (Srisawat et al. 2013).

In HR5, we define a lump of stellar particles in an SBO as a
galaxy. In principle, the MBP and the galaxy of an SBO should
be closely bound to each other, but, in practice, they are
sometimes decoupled due to complicated internal interactions,
branching into different structures in the next step (see Case 3
in Figure 1). In merger trees, this case can be seen, as if
galaxies suddenly appear or disappear with no legitimate
descendants or progenitors. To minimize such cases, we
additionally use the tree-building code ySAMtm (Jung et al.
2014; Lee et al. 2014) that forms merger tree branches based on
individual particle IDs. ySAMtm has been originally developed
to construct halo merger trees using dark-matter particles that
do not newly form or disappear throughout the duration of the
cosmological simulations. We have updated ySAMtm to be
able to deal with stellar particles that also have unique IDs but
newly formed in the simulation. ySAMtm first maps galaxies
from a snapshot to the galaxies in the next snapshot that share
at least one stellar particle. We define the main descendant/
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Flowchart showing the morphology classification scheme. We first classify galaxies into symmetric and asymmetric ones using the asymmetry parameter and divide the symmetric galaxies into disks, spheroids, and irregulars based on their radial stellar-mass-density profile parameterized by the Sérsic index $n_{\text{Sérsic}}$. Galaxies that are not fitted to the Sérsic profile are classified as irregulars together with the asymmetric ones.</figcaption></figure>

progenitor of a galaxy in a snapshot as the galaxy with which the galaxy of interest shares its largest stellar-mass fraction in the next/previous snapshot. If two galaxies in two adjacent snapshots are the main progenitor and descendant to each other, then the branch between them is classified as the main branch. A galaxy is regarded as merging into its main descendant when the main descendant has as the main progenitor a galaxy other than the galaxy of interest.

There are three cases that need to be addressed to clarify the evolutionary tracks of galaxies. (i) A galaxy is classified as new when all its stellar particles form between the current and previous time steps, and it has no progenitors (Case 2 in Figure 1). (ii) Some galaxies have stellar particles that are not bound to any galaxy in the previous snapshot but had already formed (Case 4 in Figure 1). These types of galaxies are classified as originating from accretion. They are mostly contained in SBOs that have particle numbers close to the minimum required to form FoF halos in our structure-finding scheme. Therefore, such galaxies make a practically negligible contribution at galaxy masses above the completion limit of $M_{*} \sim 10^{9} M_{\odot}$ in HR5 (see Lee et al. 2021). (iii) A galaxy can have a main progenitor that connects with a main descendant that differs from the galaxy of interest (Case 5 in Figure 1). Because the galaxies cannot have legitimate progenitors in this case, they are inevitably classified as new galaxies. This frequently occurs in fly-by galaxies in dense environments because their dynamical states are deeply entangled with central galaxies. Most of them soon escape and are subsequently identified as independent SBOs. To patch the broken branches, we backtrack the MBPs of the re-identified galaxies until the galaxies are explicitly separated from other SBOs in the past, as illustrated by the green arrows in Figure 1. In this case, the legitimate progenitors are initially identified as those merging into other galaxies. The merger branches of such progenitors are revised and connected to their legitimate descendants, instead of their massive neighbors.

The merger trees are constructed based on the MBP scheme first, and branches are revised based on stellar particles if galaxies are present in the structures. Thus, the branches formed based on stellar particles are prior to those made based on the MBP, except for Case 5 in Figure 1.

### 3. Morphology of Galaxies

We measure the galaxy morphology using the asymmetry and the radial profile of their stellar-mass-density distribution. Therefore, the morphology we study here may be different from the morphology defined in observations, which rely on the light distribution and take into account the locality of star formation within galaxies. As pointed out by Trayford et al.

(2019), albeit not a direct observable, the stellar-mass-density distribution provides us with the physical properties of galaxies, which must serve as a proxy to understand the origin and evolution of galaxy morphology. We assign morphological types to galaxies using the shape characteristics only and avoiding other indirect information such as kinematics (as done by e.g., Dubois et al. 2016; Trayford et al. 2019) or star formation activities. This is because the latter are not morphological properties, even though they can acquire some correlation with morphology in certain redshift ranges. Our morphology classification scheme is summarized in Figure 2.

We use the asymmetry parameter to classify galaxies into symmetric and asymmetric types. Radial stellar-mass-density profile is then used to measure the Sérsic index, according to which symmetric galaxies are divided into disks and spheroids. Irregulars include asymmetric galaxies and those for which the Sérsic model fit fails. We do not use the axis ratios as shape parameters because they usually depend on distance from the center, are sensitive to shot noise in the outskirts, and thus are not robust.

#### 3.1. Asymmetry Parameter

As shown in Figure 2, we first calculate the asymmetry parameter A for each galaxy to decide if the galaxy has a symmetric or asymmetric morphology. We measure A by adopting a methodology similar to that proposed by Conselice et al. (2000). Specifically, we use the stellar-mass density of the galaxy defined on a three-dimensional mesh covering four times the rms size of the galaxy. The stellar-mass density $\rho_{*,i}$ at pixel i is calculated by using a spline kernel with a variable size containing 100 adjacent stellar particles. Then, we measure the asymmetry parameter from

$$A = \frac{\sum_{i} w_{i} |\rho_{*,i} - \rho_{*,j}|}{2 \sum_{i} w_{i} \rho_{*,i}}, \quad (4)$$

where the pixel index i goes through all density pixels in the mesh. Here, pixels i and j are points reflected in the galactic center, and the weight $w_i$ is 1 when at least one of the pixels i and j contain more than 10 stellar particles and is 0 otherwise. We also exclude the central region with a diameter of 1.6 pkpc, as this region is not well resolved in HR5. The top-left panels of Figures 3 and 4 show the asymmetry parameter A as a function of stellar mass $M_{*}$ for all the galaxies identified at $z=7$ and 5, respectively. The galaxies with relatively high stellar mass seem to show a dichotomy in A at around $A=0.4$ to $-0.5$. Based on this finding, we classify galaxies with $A \ge 0.4$ as asymmetric galaxies and the rest as symmetric ones.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Relations between the asymmetry parameter A, i-band concentration parameter, Sérsic index $n_{Sérsic}$, ratio of rotational velocity to velocity dispersion $V_{rot}/\sigma$, and stellar mass $M_*$ of galaxies identified at $z=7$. The short blue line in the top-left panel marks the asymmetry parameter cut of A = 0.4 dividing symmetric and asymmetric galaxies. The short red lines in the bottom panels denote the Sérsic index cut demarcating disk-type and spheroid-type galaxies. The cyan curve in the bottom-right panel shows the median $V_{rot}/\sigma$ at a given $n_{Sérsic}$, and the error bars are the uncertainty limits of the medians.</figcaption></figure>

Numerical smoothing of the density profile over the simulation resolution scale is not expected to affect the asymmetry parameter much, as we find that galaxies become irregulars when they actually undergo structure-destructing events such as close encounters or mergers. Therefore, the fraction of irregulars is not much expected to be biased.

## 3.2. Sérsic Index

We calculate the Sérsic index, $n_{Sérsic}$, defined as

$$ \Sigma(R) = \Sigma_0 \exp(-kR^{1/n_{Sérsic}}), \qquad (5) $$

where $\Sigma(R)$ is the radial profile of the stellar-mass-density projected on the galactic plane defined by the major and intermediate axes, and R is the radius in the galactic plane (Sérsic 1963). We find the axes by diagonalizing the moment of inertia tensor of the stellar system. We fit the radial-density profile of each symmetric galaxy having A < 0.4 with the function in Equation (5) to find the best-fitting $n_{Sérsic}$. The fit is made over the radial interval from 0.8 pkpc to the radius $R_{90}$ where 90% of stellar particles are contained. Just like measuring A, we exclude the same central part to avoid spurious effects from the finite resolution of HR5 and also to exclude the influence of the galactic core, if any. For some galaxies, $n_{Sérsic}$ is undefined (that is, fitting fails to converge) because their radial profile is inconsistent with a Sérsic model. They are then added to the irregular galaxy category along with asymmetric galaxies (A > 0.4). Among the regular galaxies,

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Same as in Figure 3, but for galaxies at $z=5$. Only half the sample in $x < 0.5L_{box}$, where $L_{box}$ is the side length of the HR5 volume, is plotted here to make the figure less saturated.</figcaption></figure>

those with $n_{Sérsic} > 1.5$ are called “spheroids,” and those with $n_{Sérsic} < 1.5$ are called “disks.”

This choice is made based upon the Sérsic index of observed galaxies. We have measured the Sérsic indices of the SDSS galaxies in the training set of Park & Choi (2005) and find that $n_{Sérsic} = 1.5$ divides the galaxies into early and late types with a high reliability. Figure 21 in the Appendix shows the distribution of each type of SDSS galaxy in the isophotal axis ratio versus the Sérsic index plane.

$R_{90}$ containing 90% of total galaxy stellar mass is the outer limit of the radial range where the mass-density profile is calculated and used to estimate the Sérsic index. The medians of $R_{90}$ are 2.8, 3.2, 3.6, and 4.4 pkpc at $z=7, 6, 5$, and 4 for our least-massive galaxies with $M_* \sim 2-4 \times 10^9 M_\odot$, respectively. The 10% of stellar particles at the outermost region are dropped to minimize the effects of mass contribution of possible stray galaxies and tidal perturbation of neighboring galaxies. On the other hand, the inner limit is set to $r=0.8$ pkpc taking into account the force resolution of HR5. The finite simulation resolution is expected to make galaxies less concentrated at their center. Its impact on morphology classification is mitigated because the central part is excluded. The radial interval for density profile calculation is significantly larger than the simulation resolution scale for relatively more massive galaxies with $M_* \ge 5 \times 10^9 M_\odot$ and at lower redshifts $z \le 6$.

The bottom panels of Figures 3 and 4 show the Sérsic index as a function of stellar mass (left panel) and the rotational velocity to velocity dispersion ratio ($V_{rot}/\sigma$) as a function of Sérsic index (right panel), at $z=7$ and 5, respectively. The velocity components of the stellar particles are measured in a cylindrical coordinate system in which the direction of the
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Distributions of stars and gas of three galaxies identified at $z=6$. The first and second rows show edge-on and face-on views of stellar particle distributions, respectively, the third row shows stellar-mass-density distribution colored according to their ages (bluer is younger), and the bottom panels show gas-density distribution. Gas density increases from red to blue. The scale bar is in units of proper kpc. $M_*$, A, and $n_{\text{Sérsic}}$ are stellar mass, asymmetry parameter, and Sérsic index, respectively.</figcaption></figure>

angular-momentum vector of stellar mass is chosen to be the z-axis. The rotational velocity $V_{\text{rot}}$ is the mean tangential velocity of stellar particles inside the cylindrical shell of an inner radius $R_{1/3}$ and an outer radius $R_{2/3}$, where $R_{1/3}$ and $R_{2/3}$ are the radii containing one- and two-thirds the total stellar mass measured from the central axis, respectively. The velocity dispersion $\sigma$ is defined as $\sigma = \sqrt{(\sigma_R^2 + \sigma_{\text{rot}}^2 + \sigma_z^2)/3}$, where $\sigma_R$, $\sigma_{\text{rot}}$, and $\sigma_z$ are velocity dispersion in the radial, tangential, and vertical directions. The bottom-right panel of Figure 3 shows that galaxies with large $n_{\text{Sérsic}}$ have smaller $V_{\text{rot}}/\sigma$ at $z=7$, but the relationship becomes weaker at $z=5$ in the bottom-right panel of Figure 4. It is interesting to see that the relation is tighter at high redshifts and then nearly disappears at $z<5$ for these relatively small mass galaxies. This may be because of various evolutionary effects such as galaxy merger, as examined in the Discussion section. Figures 3 and 4 also show the relationships between various physical parameters including the i-band concentration index $L_i(<0.5R_h)/L_i(<2R_h)$, where $L_i(<R)$ is the total i-band luminosity within a radius R and $R_h$ is the half-light radius in a wavelength band.

Figure 5 shows edge-on and face-on views of representative disk (left), spheroid (middle), and irregular (right) type galaxies at $z=6$. The irregular galaxy is undergoing a merger, and has a large asymmetry parameter $A=0.723$. The spheroid galaxy in the middle panels ($n_{\text{Sérsic}}=3.26$) shows a much more compact

core than that of the disk galaxy ($n_{\text{Sérsic}}=0.764$) whose stellar-mass density is extended. Figure 5 demonstrates that we are able to conduct morphology classification of galaxies identified in the HR5 simulation following the scheme shown in Figure 2.

In this section we have shown how the asymmetry parameter and Sérsic index are related with other structural and kinematic parameters. Figures 3 and 4 show that the majority of galaxies are symmetric ones with A between 0.15 and 0.4, with an extended tail of asymmetric galaxies above $A \approx 0.4$. $n_{\text{Sérsic}}$ has a continuous distribution across 1.5 and is nearly independent of $M_*$ up to $\sim 5 \times 10^{10} M_{\odot}$. $V_{\text{rot}}/\sigma$ is a decreasing function of the concentration parameter or $n_{\text{Sérsic}}$ at high redshifts, which weakens at later epochs.

As shown below, the first galaxies in the cosmic morning are dominantly disks. In addition, we find that there are galaxies having spiral arms at this epoch. For these galaxies, the spiral arms are weak in the stellar-mass distribution but are often quite prominent in the gas-density distribution. An example is shown in Figure 6. It is interesting to note that the morphological characteristics of the present-day spiral galaxies, i.e., spiral arms or disks, have counterparts in the early cosmic history, even though the mechanism for their origin may not be the same (Dubois et al. 2016; Park et al. 2019).

# 4. Results

## 4.1. Sites of Galaxy Formation in the Cosmic Morning

By tracing the dark-matter particles belonging to galaxies back to their initial position, we find that almost all high-redshift galaxies form at high-density peaks in the initial density field. We also find that the galaxy-formation site in the initial density field migrates starting from the highest-density peaks to relatively lower-density peaks.

### 4.1.1. Galaxy-formation Sites in Initial Conditions

To find out where in the initial conditions galaxies would form later at a redshift $z$, we use the dark-matter particles belonging to the galaxies identified at that redshift. Specifically, we find the initial positions of these dark-matter particles and define the initial center, that is, the galaxy-formation site, for each galaxy by the center of mass of the dark-matter particles’ initial positions.

We then match the initial galaxy-formation sites to the peaks in the initial density field as follows. First, we smooth the initial density field with a Gaussian filter over the smoothing length $R_G = 0.35$ cMpc, which is 2.73 times the density pixel size. Then, density peaks are identified. Starting from the most-massive galaxies, we search for the highest unmatched peak within the $r=4.5R_G = 1.575$ cMpc ball around each initial galaxy-formation site. The search proceeds for the lighter galaxies after excluding the peaks already chosen. With this method, we have successfully identified, for example, the 99.6% of the galaxies with $M_* > 2 \times 10^9 M_{\odot}$ at $z=6$ to the initial density peaks.

Figure 7 shows the initial dark-matter-density field of HR5 in a 1.66 cMpc-thick slice. Contour levels delineate the dark-matter overdensity threshold $\nu = (\rho_{\text{DM}}/\bar{\rho}_{\text{DM}} - 1)/\sigma_{\text{DM}}$ in units of the standard deviation $\sigma_{\text{DM}}$. Symbols are the initial center of mass of dark-matter particles that belong to the dark-matter halos of galaxies identified at $z=7$ (circle), 6 (triangles), and 5 (crosses). The overlap of all the three symbols on top of a peak with $\nu_{\text{pk}}=5.34$ at the 7 o'clock position from the center of Figure 7
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Gas density (top) and stellar particle distribution of a small group at $z=6$. The group of galaxies is connected to each other via a gaseous filament. The central galaxy has a stellar mass of $M_{\star} = 1.76 \times 10^{10} M_{\odot}$, and its gas shows spiral arms. Small objects in the neighbor have a stellar mass smaller than $10^{9} M_{\odot}$.</figcaption></figure>

indicates that the galaxies identified at all three redshifts share the same initial density peak. We also find that the formation sites for the $z \ge 7$ galaxies correspond to the locally highest peaks residing on globally overdense regions. One typical example is shown in the Figure 7 as a circled region. This must be contrasted to the formation sites for $z=5$ galaxies (crosses) that correspond to relatively lower initial peaks that are sitting on filaments connecting the higher overdensity regions, showing the growth and emergence of large-scale structures of galaxies out of the underlying evolving matter distribution.

### 4.1.2. Initial Peak Height and Galaxy Properties

The top panels of Figure 8 show the relation between the stellar mass of the galaxies identified at $z=8, 7, 6$, and 5 (from left to right), and the peak height of the corresponding initial formation sites. Two evolutionary features are clearly seen here: the migration of the formation sites toward low-density regions and the increase in galaxy stellar mass. It should be noted that at a particular epoch, the galaxies with $M_{\star} \approx 2 \times 10^{9} M_{\odot}$ have originated from the initial density peaks with a narrow range ($\Delta v_{\text{pk}} = 0.23, 0.32$, and 0.53 from the mean $\bar{v}_{\text{pk}} = 4.60, 4.11, 3.58$ at $z=7, 6$, and 5, respectively) in their initial peak height.

To quantify the evolution, we show the mean and standard deviation of the peak heights as a function of redshifts for two narrow stellar-mass ranges in the left panel of Figure 9. The plot shows that galaxies with a fixed stellar mass form at lower and lower initial density peaks as a redshift decreases. The galaxies having $M_{\star} = 2.0-2.4 \times 10^{9} M_{\odot}$ at an epoch $z$ have formed typically out of the initial density peaks with mean height $\bar{v}_{\text{pk}} = 0.472z + 1.32 \pm 0.01$. Massive galaxies are just those that have formed earlier (at higher redshifts) with smaller mass, then acquired stellar mass through internal star formation, accretion, or merger. The galaxies with $M_{\star} = 5-6 \times 10^{9} M_{\odot}$ at $z=5$, for example, are roughly those with $M_{\star} = 2-2.4 \times 10^{9} M_{\odot}$ at $z=6$, and have formed out of the initial density peaks with $v_{\text{pk}} \approx 4.1$.

A close inspection of the distribution of initial peak heights supports the above scenario. The three solid curves in the right panel are the distributions of $v_{\text{pk}}$ of the peaks corresponding to the galaxies with $M_{\star} = 2-2.4 \times 10^{9} M_{\odot}$ at $z=8$ (top), 7 (middle), and 6 (bottom). The bottom curve for $z=6$ is quite similar to the dashed curve, which is the distribution function of galaxies with $M_{\star} = 5-6 \times 10^{9} M_{\odot}$ at $z=5$, supporting again the above observation. The initial density peaks with $v_{\text{pk}} \approx 4.1$ collapse to form $M_{\star} = 2-2.4 \times 10^{9} M_{\odot}$ galaxies by $z=6$, and grow to reach $M_{\star} = 5-6 \times 10^{9} M_{\odot}$ without losing their identity until $z=5$. On the other hand, the distribution function for those with $M_{\star} = 5-6 \times 10^{9} M_{\odot}$ at $z=6$ (dotted curve) is somewhat skewed to high $v_{\text{pk}}$ and is rather different from the middle solid curve, even though the mean peak heights are close ($\approx 4.6$).

The second and third rows of Figure 8 show that the Sérsic index and $V_{\text{rot}}/\sigma$ of galaxies do not clearly depend on the initial peak density. We will return to this question of the effects of the initial conditions on galaxy morphology and internal kinematics in Section 4.4.

How is galaxy bias related with the peak height of the initial formation site? Figure 10 shows the galaxy two-point correlation functions at three different redshifts ($z=7, 6, 5$) for a fixed stellar-mass cut ($M_{\star} > 2 \times 10^{9} M_{\odot}$, top panel) and for a fixed total number of galaxies ($N_{\text{gal}}=439$, bottom panel). As shown in Figure 9, again, for a fixed stellar-mass limit, the high-$z$ galaxies have formed at the initial density peaks with higher $v_{\text{pk}}$, and the galaxies formed at higher $v_{\text{pk}}$ (at the rarer peak) are more biased and strongly clustered. We can see this trend most clearly for $r < 20$ cMpc in the upper panel of Figure 10, where the clustering amplitude is larger for higher-redshift galaxies. On the other hand, when fixing the total number of galaxies, the significance $v_{\text{pk}}$ remains approximately the same. That is, the same galaxies maintaining their identity are sampled at the three redshifts. In this case, we may approximate that the galaxy bias evolves passively (Fry 1996; Tegmark & Peebles 1998), in which case the density contrast of highly biased ($b \gg 1$) galaxies stays constant. As a result, the two-point correlation function of these galaxies does not evolve in time, as shown in the bottom panel of Figure 10.

### 4.1.3. Emergence of the First Large-scale Structures

We note that the large-scale structures in the universe as traced by galaxies emerge in the cosmic morning in a way similar to the process that rhizomes spread out roots and shoot out culms from their nodes on the way. The burst of galaxy formation starts first in the highest-density regions, the locations of the first cosmic rhizomes, in which the first galaxies appear at the rare highest peaks. As the underlying dark matter and gas collapse into filaments, the formation sites
<figure><img src="image_9.png" /><figcaption>**Figure 7.** A 1.66 cMpc-thick slice of HR5's initial dark-matter-density field. Gray-scale map: the highest density within the slice along the direction perpendicular to the slice, after being smoothed over 0.35 cMpc with the Gaussian filter. Color lines: contour levels correspond to the overdensity in units of the rms overdensity, $\nu_{\text{DM}} \equiv (\rho_{\text{DM}}/\bar{\rho}_{\text{DM}} - 1)/\sigma_{\text{DM}} = 1.5$ (blue), 2.0 (green), 2.5 (orange), 3.0 (magenta), and 3.5 (red). Symbols: the formation sites of the first galaxies, defined as the initial center of mass of dark-matter particles that form the dark-matter halos of galaxies identified at $z=7$ (the circle), 6 (triangles), and 5 (crosses). The large-scale structures in the universe emerge at the highest-density regions and grow in time as galaxy-formation sites move from high- to relatively lower-density peaks.</figcaption></figure>

of the next new galaxies propagate away from the highest-density regions along attached filaments, increasing the realm of each rhizome; see Figure 7.

Galaxies at the outer boundary of rhizomes are in a different formation and evolution environment from the earlier ones and therefore may have different physical properties. The initial birth places of the first rhizomes are likely to evolve to protoclusters, as galaxies flow into them along filaments or as filaments collapse onto them. As time passes and matter fluctuations grow, initially lower-density peaks also become nonlinear and new less prominent rhizomes appear. And the first-generation rhizomes become multiply connected and continue to spread out toward less dense regions. The full cosmic web of galaxies appears when most cosmic rhizomes percolate with one another and become globally connected (see also Gott et al. 1986; Shandarin & Zeldovich 1989; Park 1990; Vogeley et al. 1994; Bond et al. 1996). The picture depends strongly on the mass of galaxies used for defining large-scale structures. For example, the most-massive galaxies live in disconnected regions even today and never form a web structure.

The results that the initial $\nu_{pk}$ of newly forming galaxies decreases in time (Figures 8 and 9) and the monotonic decrease of the correlation function of galaxies with fixed stellar mass (Figure 10) quantitatively support the cosmic rhizome picture

of the growth of large-scale structures in the universe. Such a description of cosmic-web formation is broadly consistent with the theories based upon the threshold density (Kaiser 1984), such as peak theory (Bardeen et al. 1986) or excursion set theory (Bond et al. 1991), although the detailed mathematical modeling requires further understanding of galaxy-formation physics (following, e.g., Musso et al. 2018, or more generally a formulation in terms of redshift-dependent galactic bias). Yet, it does not fully capture the dynamical top-down scale coupling for the gas, which funnels galaxy morphology through cold streams (Pichon et al. 2011), as discussed below.

## 4.2. Morphological Types of Galaxies

As described in Section 3, we have classified galaxy morphology based on the asymmetry parameter and the Sérsic index measured from the radial stellar-mass-density profile. We measure the fraction of disks ($n_{\text{Sérsic}} < 1.5$), spheroids ($n_{\text{Sérsic}} > 1.5$), and irregulars ($A > 0.4$ or undefined $n_{\text{Sérsic}}$) in each mass bin at three redshifts. The dotted, dashed, and solid lines in Figure 11 are the fractions at $z=7, 6, \text{and } 5$, respectively. It is clear that the dominant morphological type of the galaxies in the cosmic morning is the disk type. The fraction of the disk type is between 0.63 and 0.69 for the galaxies with $M_{\star} > 2 \times 10^9 M_{\odot}$ at the three epochs. The
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Stellar mass (top), Sérsic index (middle), and rotational velocity to velocity dispersion ratio (bottom) of galaxies at redshifts 8, 7, 6, and 5 as a function of the initial dark-matter overdensity peak height $v_{\text{pk}} = (\rho_{\text{DM,pk}} - \bar{\rho}_{\text{DM}})/\sigma_{\text{DM}}$ matched to the initial center of mass of the dark matter in galaxies at each epoch.</figcaption></figure>

fraction of the spheroid type is also quite stable and is between 0.19 and 0.14 for the same mass range. Interestingly, the fraction of irregulars is not high and is only a little higher than that of spheroids. However, the irregular fraction increases at the highest-mass end at each redshift, which indicates that the relatively more massive galaxies are likely to be undergoing structure-destroying interactions.

From Figure 11, we conclude that the fraction of disk plus irregular morphological types does not depend much on stellar mass over the range we have explored. For example, the fraction is ~0.87 within the range of $M_* = 2-30 \times 10^9 M_{\odot}$ at $z=5$. On the other hand, the observed fraction decreases sharply for local galaxies from ~0.85 to ~0.25 over the same mass range (Figure 4 of Kelvin et al. 2014). This fact may imply that at high redshifts, the mechanism determining galaxy morphology is different from the local one. The fraction of spheroid-type galaxies (~0.14 at $z=5$) found in this study is in good agreement with the fraction of early-type galaxies observed at $z \approx 1$ as measured by Hwang & Park (2009), who have found using the GOODS data that the fraction converges to ~0.15 in both high- and low-density regions near $z=1$ (Figure 11 in Hwang & Park 2009). We note that the observed early-type fraction of about 0.15 at $z \approx 1$ has already been set up in HR5 at high redshifts even though the early types at $z>4$ are most probably not those at $z \approx 1$ according to our finding of frequent morphology transformation of high-redshift galaxies presented in the next section.

Using the EAGLE simulation, Trayford et al. (2019) has reported a distribution of high-redshift galaxy morphology very different from ours. For example, Figure 6 of their paper shows that irregulars (i.e., asymmetric galaxies) dominate at high redshifts and all the way down to $z \approx 1.4$. The difference comes from their low choice of asymmetry parameter cut of $A = 0.2$ for irregulars. The choice is made from a visual inspection of galaxy images in their Figure 3. In the current study, the distribution of the measured $A$ is used to make the choice and the chosen threshold separates the disturbed galaxies from the large group of normal galaxies as seen in Figures 3 and 4. To some extent, the division of irregular galaxies is a matter of subjective choice. In Trayford et al., the fractions of disks, spheroids, and irregulars are about 14%, 24%, and 63% at $z=6$, which are very different from our results. It should be pointed out that their morphology classification of disk and spheroid types is actually based not on the galaxy stellar-mass distribution but on kinematics parameterized by the circularity $\epsilon$. As we show throughout this paper, the relation between morphology and internal kinematics is far from tight in the cosmic morning.

## 4.3. Effects of Evolution

### 4.3.1. Evolution of the Sérsic Index of Individual Galaxies

By closely examining the galaxies in HR5, we find that the spheroidal and irregular morphologies of galaxies in the cosmic morning are incidental and transient phenomena triggered by
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Left: initial dark-matter overdensity peak height $\nu_{pk}$ corresponding to the center of mass of the dark matter in the galaxies having a given stellar mass. The solid and dashed lines are for the galaxies with $M_{*} = 2-2.4 \times 10^{9} M_{\odot}$ and $M_{*} = 5-6 \times 10^{9} M_{\odot}$, respectively. The error bars show the 16th–84th percentile distributions. The open circles mark the median $\nu_{pk}$ of spheroids, and the asterisks are for disk-type galaxies. Right: histograms of galaxy number count as a function of $\nu_{pk}$. The cyan, blue, and magenta solid curves correspond to the galaxies with $M_{*} = 2-2.4 \times 10^{9} M_{\odot}$ at $z = 8, 7, $ and $6$. The magenta dotted and red dashed curves show the distributions of the galaxies with $M_{*} = 5-6 \times 10^{9} M_{\odot}$ at $z = 6$ and $5$, respectively.</figcaption></figure>

galaxy–galaxy interactions. To illustrate this, we show in Figure 12 the evolution of the Sérsic index of a few galaxies as a function of their stellar mass through $z=4$. The red and blue lines in the top panel are the tracks of two most-massive galaxies identified at $z=4$ after their stellar mass reached $2 \times 10^{9} M_{\odot}$. The dotted lines indicate the periods when the Sérsic index is not measurable or the asymmetry parameter is greater than 0.4. For the most-massive galaxy (red line), some specific redshifts are marked. The middle and bottom panels show the tracks of two randomly selected galaxies with $M_{*} \approx 5$ and $2.5 \times 10^{10} M_{\odot}$ at $z=4$, respectively. The gray horizontal lines ($n_{\text{Sérsic}} = 1.5$) are the demarcation dividing spheroids and disks in our study.

It can be immediately seen that galaxy morphology in the cosmic morning is not a fixed property but fluctuates in time. Galaxies spend most of their lives as disk types, but frequently change their morphology to spheroid or irregular type. One interesting phenomenon to note is that galaxies tend to recover their disk morphology after they become spheroids or irregulars. This tendency is likely to be due to reacquisition of angular momentum from the vorticity-rich surrounding material (see, e.g., Laigle et al. 2015; Danovich et al. 2015), as we discuss below.

The Sérsic index of relatively low-mass galaxies with $M_{*} \le 5 \times 10^{10} M_{\odot}$ at $z=4$ stays under the division line most of the time and only occasionally crosses the line. The tracks of the more massive galaxies often have dotted lines indicating the periods of significantly perturbed or destructed structure. Some of them are associated with merger or breakup events resulting in abrupt mass changes, but not always so. When galaxies

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Top: the two-point correlation function for the galaxies with $M_{*} \ge 2 \times 10^{9} M_{\odot}$ at three different redshifts, $z=7$ (Magenta), $6$ (Green), and $5$ (Blue). Bottom: the two-point correlation function at the same three redshifts but with different stellar-mass cuts to sample the same number of the most-massive galaxies, giving the fixed number of galaxies of $N_{\text{gal}} = 439$. For each panel, we also show the linear matter two-point correlation function at the same redshifts as a solid line ($z=5$), dotted line ($z=6$) and dashed line ($z=7$). The stronger evolution of the small-scale nonlinear bias shown for the fixed stellar-mass-limit case (top) has diminished when fixing the total number of galaxies.</figcaption></figure>

undergo destructive interactions resulting in irregular morphology, the Sérsic index tends to increase. Among the galaxies with $M_{*} = 5-25 \times 10^{9} M_{\odot}$ at $z=6$ that have suffered from structure-destructing events $56\% \pm 4\%$ have increased their Sérsic indices after each event while $44\% \pm 3.7\%$ have decreased it. On the other hand, it is nearly equally likely for the Sérsic index to increase or decrease when all major stellar-mass jumps, $\Delta M_{*}/M_{*} > 0.25$ per simulation output time step, are considered regardless of morphology transformation to irregular type.

Figure 13 shows the change in the Sérsic index of galaxies and fractional stellar-mass change that occur at each output time step of $\Delta t \approx 2.5 \times 10^{7}$ yr. The time steps between $z=4$ and the epoch when a galaxy first reaches $M_{*} = 2 \times 10^{9} M_{\odot}$ are used for the plot. The upper redshift limit thus varies among
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Fractions of disk (blue), spheroid (red), and irregular (black) morphological types, at redshifts 5 (solid), 6 (dashed), and 7 (dotted) as a function of galaxy stellar mass.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Evolution tracks of galaxies in Sérsic index vs. stellar-mass plane. The top panel shows the evolution of two most-massive galaxies at $z=4$. The middle and bottom panels exhibit the tracks of the galaxies with $M_{*} \approx 5$ or $2.5 \times 10^{10} M_{\odot}$ randomly chosen at $z=4$. All the galaxies are traced after they reach $M_{*} = 2 \times 10^{9} M_{\odot}$. The dotted lines mark the periods when galaxy morphology is irregular. Also shown are the time steps when the mass of the most-massive supermassive black hole among the black holes belonging to a galaxy increases more than 50% (filled circles) or 25% (open circles) per $2.5 \times 10^{7}$ yr.</figcaption></figure>

galaxies. The galaxies that have $M_{*} = 5-6 \times 10^{9} M_{\odot}$ (blue points) or $M_{*} > 6 \times 10^{9} M_{\odot}$ (red points) at $z=6$ are selected and traced along the merger tree. Cyan and magenta lines are the mean values for the two mass intervals, respectively, and error bars are the uncertainties of the means. The Sérsic index

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Sérsic index change $\Delta n_{\text{Sérsic}}$ vs. fractional stellar-mass change $\Delta M_{*}/M_{\text{star}}$ per unit output time step of $\Delta t \approx 2.5 \times 10^{7}$ yr. Each dot corresponds to a particular time step between $z=4$ and the epoch when each galaxy first reaches $M_{*} = 2 \times 10^{9} M_{\odot}$. The blue dots indicate galaxies with $M_{*} = 5-6 \times 10^{9} M_{\odot}$, and the red dots are for those with $M_{*} > 6 \times 10^{9} M_{\odot}$ at $z=6$. The cyan and magenta lines show the mean values of $\Delta n_{\text{Sérsic}}$ for blue and red points, respectively. The error bars exhibit the $1\sigma$ distribution of the mean curves.</figcaption></figure>

shows no trend with mass change when $\Delta M_{*}/M_{*} > 0.15$. We conclude that the destructive interactions resulting in spheroids or irregulars are not always identical to merger events. On the other hand, there is a very weak but statistically significant positive relation between $\Delta M_{*}/M_{*} = 0.05$ and 0.15 (more than $6\sigma$). This means that when $\Delta M_{*}/M_{*} \lesssim 0.15$, galaxies on average tend to have a more positive Sérsic index change for a larger fractional stellar-mass increase even though the change at each step appears stochastic. This is the regime of internal stellar-mass growth due to star formation within galaxies or accretion of low-mass objects. The star formation rate of the galaxies with $M_{*} = 5 \times 10^{9} M_{\odot}$ at $z=6$ is typically about $20 M_{\odot}/\text{yr}$ and they increases stellar mass of about $5 \times 10^{8} M_{\odot}$ or $\Delta M_{*}/M_{*} \approx 0.1$ during $2.5 \times 10^{7}$ yr.

In Figure 12, we mark the time steps when the mass of the most-massive supermassive black hole (SMBH) among the black holes belonging to a galaxy increases more than 50% (filled circles) or 25% (open circles) per $2.5 \times 10^{7}$ yr. Studies with much larger samples need to be made to draw any conclusions on whether or not the sudden growth of SMBH mass is correlated with morphological transformation at high redshifts.

### 4.3.2. Effective Morphology

Because galaxy morphology at a particular epoch cannot represent the shape of galaxies over an extended period, we instead introduce the fraction of the time period during which a galaxy has spent its life as a spheroid, to assign a representative galaxy morphology. That is, we define

$$f_{\text{age,spheroid}}(z) = \frac{\int_{t_i}^{t(z)} \Theta(n_{\text{Sérsic}}(t) - 1.5)dt}{t(z) - t_i}, \quad (6)$$

where $t_i$ is the time when the galaxy reached stellar mass of $2 \times 10^{9} M_{\odot}$, and $\Theta$ is the step function. Similarly, we define the
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Fraction of age that galaxies have spent their lives as spheroids by $z=6$ (top) and 4 (bottom) after their stellar mass reaches $M_{\star} = 2 \times 10^{9} M_{\odot}$. The abscissa is the fraction of stellar mass increased through mass jumps typically accompanied by mergers (see the text). The blue and red dots indicate galaxies having $M_{\star} = 5 - 6 \times 10^{9} M_{\odot}$ and $M_{\star} > 6 \times 10^{9} M_{\odot}$ at $z=6$, respectively. There are 44 (20) blue and 38 (18) red points at the origin out of a total 278 (146) points in the upper (lower) panel. The dashed curves present the medians at given $f_{\text{mass,merger}}$, and error bars are the uncertainties of the medians.</figcaption></figure>

fraction of stellar mass that has increased through mergers by a given redshift $z$,

$$f_{\text{mass,merger}}(z) = \frac{\int_{z_i}^{z} \Theta(\Delta M_{\star}/M_{\star} - 0.25) dM_{\star}(z')}{M_{\star}(z) - M_{\star}(z_i)}, \quad (7)$$

where the mass change $\Delta M_{\star}$ is again measured for our output time steps. Here, the stellar-mass increase $\Delta M_{\star}$ not only includes the stellar mass directly carried to the major galaxy by the infalling minor one but also that immediately formed after merger events. Therefore, the criterion $\Delta M_{\star}/M_{\star} > 0.25$ per time step does not always find the conventional major mergers used to construct the merger tree but also includes some minor mergers that trigger strong starburst.

Figure 14 shows the fraction of time that galaxies live as spheroids from the time when they acquired stellar mass of $2 \times 10^{9} M_{\odot}$ through $z=6$ (upper panel) or 4 (bottom panel). The abscissa is $f_{\text{mass,merger}}$ integrated through the corresponding redshifts. The blue and red points are again the galaxies having stellar mass of $5-6 \times 10^{9} M_{\odot}$ or greater than $6 \times 10^{9} M_{\odot}$ at $z=6$, respectively. They indicate that both effective morphology and merger-driven stellar-mass fraction do not depend on the stellar mass of galaxies. Note that there are many points at (0, 0) overlapping one another. They are the galaxies that have never had spheroidal morphology and never experienced major stellar-mass jumps. Dashed curves are the median $f_{\text{age,spheroid}}$.

$f_{\text{mass,merger}}$ is on average a little higher in the upper panel, because at the beginning of cosmic morning, galaxies more frequently merge with other galaxies with similar mass. As can

be seen from the median curves, the majority of galaxies spend their lives as late types (i.e., disks or irregulars) and the probability for them to live as spheroid is essentially independent of $f_{\text{mass,merger}}$ as long as $f_{\text{mass,merger}}$ is not zero. Interestingly, the galaxies that have never experienced major mass jumps have $f_{\text{age,spheroid}} = 0.077_{-0.005}^{+0.011}$ at $z=6$, which is somewhat lower than $f_{\text{age,spheroid}} = 0.169_{-0.019}^{+0.030}$ for those that have ever experienced major mass jumps.

In most cases, our statistical study includes the galaxies with $M_{\star} = 2-10 \times 10^{9} M_{\odot}$ in the cosmic morning ($z>4$). Our results suggest that these galaxies are typically disks and temporally become irregulars or spheroids. The mechanism transforming galaxies from disks to spheroids on average is not simply mass-jumping merger events, as such mergers do not preferentially make galaxies more spheroidal in this mass range. Instead, a statistical trend for transformation to spheroidal type is observed when galaxies go through an irregular morphology phase. Various kinds of cataclysmic events like close encounters, major and minor mergers, and episodic strong accretion that manage to affect galaxy structure significantly seem responsible for the trend. Welker et al. (2014) has showed that the magnitude of the spin of simulated galaxies increases steadily and aligns itself preferentially with the nearest filament when no significant merger occurs. Strong internal star formation gives a weak tendency toward spheroidal type, too (Martin et al. 2018).

### 4.4. Effects of the Initial Conditions

It is expected that the monolithic collapse of mass around an initially very-high-density region ends up producing galaxies morphologically different from those formed in filaments or walls. It is also sensible to expect for kinematically different regions that are in different tidal force environment (Codis et al. 2015, their Figure 12) to produce morphologically different galaxies. In this section, we examine how the initial conditions impact galaxy properties in the cosmic morning.

#### 4.4.1. Effects of the Initial Peak Height

We inspect in Figure 15 if the initial density peak height affects galaxy mass, morphology, merger history, and kinematic property of galaxies. The abscissa is the height of the initial density peak matched with the initial center of mass of the dark-matter particles composing the dark-matter halo of each galaxy identified at $z=6$. The initial dark-matter-density field represented on a 0.128 cMpc pixel mesh and Gaussian smoothed over $R_G = 0.35$ cMpc is used to identify the initial density peaks and to measure their height $\nu_{pk} \equiv (\rho_{\text{DM}}/\bar{\rho}_{\text{DM}} - 1)/\sigma_{\text{DM}}$, where $\sigma_{\text{DM}}$ is the rms density fluctuation. Plotted points are for the galaxies that have $M_{\star} = 5-6 \times 10^{9} M_{\odot}$ (left panels) or $M_{\star} > 6 \times 10^{9} M_{\odot}$ (right panels) at $z=6$. The blue and red points in the top and middle panels are the fractions calculated through $z=6$ and 4, and the red/blue lines with error bars are the medians of the red/blue dots, respectively.

In the middle panels, $f_{\text{mass,merge}}$ is the fraction of galaxy stellar mass that has increased during mass jumps. Here we lower the threshold for the mass jump to $\Delta M_{\star}/M_{\star} = 0.1$ per time step to include more minor mergers, so it is defined somewhat differently from that in Figure 14. The threshold is lowered because the majority of galaxies have zero amount of mass acquired through major mass jumps, and it is difficult to
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Fraction of age that galaxies spend their lives as spheroids (top), fraction of stellar mass increased during mergers (middle), and rotational velocity to velocity dispersion ratio at z = 6 (bottom) as a function of the height $v_{pk}$ of the initial dark-matter overdensity peak matched to galaxies. Density peaks are defined on the “galaxy” scale of $R_G = 0.35$ cMpc. The blue and red dots in the top and middle panels are the fractions integrated through z = 6 and 4, respectively, after their stellar mass reaches $2 \times 10^9 M_\odot$. The left panels are for the galaxies with $M_* = 5-6 \times 10^9 M_\odot$, and the right panels for those with $M_* > 6 \times 10^9 M_\odot$ at z = 6. Median curves and their uncertainty limits are given.</figcaption></figure>

see any trend with $v_{pk}$. The stellar mass increase of 10% during 0.025 Gyr is the mass that galaxies at this epoch typically increase just through internal star formation.

The top panels of Figure 15 show that galaxy morphology is independent of the height of the initial density peaks that end up forming the first galaxies. This contrasts with the trends that merger fraction and $V_{rot}/\sigma$ do have some dependence on $v_{pk}$ (middle and bottom panels). Higher-initial-density peaks tend to suffer from less amount of mergers by z = 6 or 4. Higher-density peaks tend to have higher $V_{rot}/\sigma$ with a marginal statistical significance of $2.6\sigma$. However, these trends in merger history and kinematic property are not related to galaxy morphology. We also do not see any clear dependence of morphology on age (blue versus red points) or mass (left versus right panel). Note that no galaxy in our sample has a stellar mass higher than $3 \times 10^{10} M_\odot$ before z = 6 as can be seen in Figure 8.

### 4.4.2. Effects of the Large-scale Background Density

To look into the dependency of galaxy morphology on large-scale density environment, in Figure 16, we plot $f_{age,spheroid}$ and $f_{mass,merger}$ as a function of the initial overdensity contrast ($\rho_{DM}/\bar{\rho}_{DM} - 1)/\sigma_{DM}$ smoothed with a top-hat filter over the “galaxy protocluster” scale of 4 cMpc. The corresponding smoothing mass is $1.0 \times 10^{13} M_\odot$. We chose a sample of galaxies identified at z = 6 having stellar mass within a narrow range of $5-6 \times 10^9 M_\odot$. The fractions are integrated from the epoch when their stellar mass reached $2 \times 10^9 M_\odot$ to z = 6 (blue dots) or 4 (red dots).

<figure><img src="image_18.png" /><figcaption>**Figure 16.** Fraction of age that galaxies spend their lives as spheroids (top) and stellar-mass fraction increased during mergers (bottom) integrated through z = 6 (blue) and 4 (red) because galaxy mass reaches $M_* = 2 \times 10^9 M_\odot$. The abscissa is the initial background overdensity contrast ($\delta_{initial, bg} \equiv (\rho_{DM}/\bar{\rho}_{DM} - 1)/\sigma_{DM}$) smoothed with a top-hat filter over the “protocluster” scale of $R_{TH} = 4$ cMpc.</figcaption></figure>

It can be seen in the upper panel that galaxy morphology is basically independent of the initial large-scale background density and that spheroidal shape is a minority everywhere and late-type (disk plus irregular) morphology dominates both in high- and low-density regions in these early periods. The stellar-mass fraction of galaxies acquired during mass jumps is also nearly independent of the initial background density by z = 6. Galaxies grow mainly via mergers by that epoch as presented by the blue crosses in the lower panel. However, merger events decrease in time at different rates in different density environments, and by z = 4 (stellar mass), accretion and internal star formation rather than mergers become the main channel of stellar-mass growth of the galaxies in low-background-density regions (red dots). This environmental dependence of the merger/interaction rate is expected to be responsible for the environmental galaxy mass and morphology dependence at later epochs.

### 4.4.3. Impact of the Initial Tides

Next we examine if galaxy morphology, merger fraction, and kinematics depend on the primordial angular momentum of galaxies. According to tidal torque theory (Doroshkevich 1970; White 1984), the large-scale matter distribution can exercise tidal torque on a protogalactic region, producing the angular momentum (hereafter AM) of the galaxy forming there. In the Zel'dovich approximation of the Lagrangian perturbation theory (Zel'dovich 1970b), the AM of a mass distribution is given as follows (Sugerman et al. 2000; Schäfer & Merkel 2012):

$$ \mathcal{L}_{\alpha}(\mathbf{q}_c) = a^3 H(a) \frac{dD}{da} \epsilon_{\alpha\beta\gamma} \sum_{\sigma} I_{\beta\sigma} \Psi_{\sigma\gamma}(\mathbf{q}_c), \quad (8) $$
where a is the scale factor, H is the Hubble expansion rate, D is the linear growth function that encodes the time evolution of the growing mode solution of linear density perturbation, $\epsilon_{\alpha\beta\gamma}$ is the Levi-Civita symbol, and $q_c$ is the center of mass. The moment-of-inertia tensor in terms of the initial matter position $q$ is

$$I_{\beta\sigma} = \rho_0 a^3 \int_V d^3 q (q - q_c)_\beta (q - q_c)_\sigma, \quad (9)$$

where the integration is made over the protogalactic volume V. The initial shear tensor is given by

$$\Psi_{\sigma\gamma}(q_c) = \partial_{\sigma}\partial_{\gamma}\Psi(q_c), \quad (10)$$

where $\Psi(q_c)$ is the displacement potential in the Zel'dovich approximation and we have assumed that the displacement field is irrotational. As we are interested only in the relative difference of AM, all the position-independent terms will be omitted in our evaluation of angular momentum. As the scalar displacement field is proportional to the gravitational potential in the linear approximation, we adopt a quantity proportional to the primordial AM as

$$L_{\alpha} \equiv \epsilon_{\alpha\beta\gamma} \sum_{\sigma} I_{\beta\sigma} \partial_{\sigma} \partial_{\gamma} \Phi(q_c), \quad (11)$$

where $\Phi$ is the gravitational potential field. We calculate the AM in the following way. We identify galaxies at $z=6$ and make a list of all dark-matter particles belonging to the dark-matter halo of each galaxy. We then use their initial positions to calculate the center of mass $q_c$ and the moment of inertia tensor $I_{\beta\sigma}$. The gravitational potential field is smoothed with a Gaussian filter over $R_G = 0.35$ cMpc, and the shear tensor $\partial_{\sigma}\partial_{\gamma}\Phi$ is calculated at $q_c$ by using the finite difference method. The moment of inertia tensor and the gravitational shear tensor are used to calculate the AM in Equation (11).

The top panels of Figure 17 show that galaxy morphology is basically independent of the initial AM of protogalactic regions. That is, disk morphology dominates regardless of the magnitude of L. Here $f_{\text{age,spheroid}}$ is nearly independent of L both for relatively low- and high-mass galaxies, even though the fraction seems to develop a moderately increasing trend with L between $z=6$ and 4 for massive galaxies (i.e., red line in the top-right panel). Even though the galaxy morphology does not show a dependence on the initial AM, $f_{\text{mass,merge}}$ and $V_{\text{rot}}/\sigma$ of stars are weakly increasing function of L as shown in the middle and bottom panels.

Evidence for a positive correlation between the initial AM and internal kinematics and for the negligible correlation between L and galaxy morphology can also be seen in Figure 18. The rightmost panels of Figure 18 presents histograms showing the alignment of the initial AM with the AM of total mass ($L_{\text{tot}}$, top panel), stellar component ($L_*$, middle panel) and with the minor axis of the moment of inertia of the stellar-mass distribution ($\hat{e}_3$, bottom panel) of galaxies at $z=6$. It can be seen that the primordial AM is tightly aligned with the total AM and rather well aligned with the AM of the stellar component. $\hat{L} \cdot \hat{L}_{\text{galaxy}}(z=6)$ is greater than 0 for 88.3% or 73.6% of galaxies when total (dominantly dark matter) or stellar mass is used for AM calculation, respectively. However, $\hat{L}$ shows only a very weak alignment with the short axis of the stellar-mass distribution. It should be pointed out again that the initial L is tightly aligned with galaxy internal motion but not much aligned with the mass shape.

<figure><img src="image_19.png" /><figcaption>**Figure 17.** Fraction of age that galaxies spend their lives as spheroids (top), stellar-mass fraction increased during mergers (middle), and rotational velocity to velocity dispersion ratio of stars at $z=6$ (bottom) as a function of the initial angular momentum (in arbitrary units). Blue and red points in the top and middle panels are the fractions integrated through $z=6$ and 4, respectively, after their stellar mass reaches $2 \times 10^9 M_{\odot}$. The left panels are for the galaxies with $M_* = 5-6 \times 10^9 M_{\odot}$, and the right panels for those with $M_* > 6 \times 10^9 M_{\odot}$ at $z=6$.</figcaption></figure>

This suggests that one should not use internal kinematics to define galaxy morphology at high redshifts.

The spin parameter $\lambda'$ of galaxies at $z=6$ in Figure 18 is defined as

$$\lambda' = \frac{\mathcal{L}}{\sqrt{2} M v R}, \quad (12)$$

where $\mathcal{L}$ is the total galaxy AM, M is virial mass, R is the virial radius, and v is circular velocity (Bullock et al. 2001; Aubert et al. 2004; Knebe & Power 2008; Schäfer & Merkel 2012). The alignment is somewhat stronger for galaxies with higher L or $\lambda'$, as shown by the red dashed median lines.

A possible interpretation is that the short axis of the stellar mass basically shows nearly no alignment with the initial AM because the morphology of the stellar component of galaxies in the cosmic morning can be easily and significantly affected by perturbations, and hence has a limited memory of the initial conditions. On the other hand, the AM of the stellar mass of galaxies does show some alignment, and the AM of the total mass has even stronger alignment with the initial AM. As the (gas-dominated) mass reservoir of galaxies remembers the initial AM produced by the tidal torque of the large-scale dark-matter distribution, later infalling gas onto the galaxies tends to reset the original disk morphology (e.g., Pichon et al. 2011; Stewart et al. 2013).

## 5. Discussion

In this work, we have found that the majority (~2/3) of galaxies in the cosmic morning before $z=4$ have a disk morphology, and galaxies transiently take a spheroidal or irregular morphology when they experience destructive interactions.
<figure><img src="image_20.png" /><figcaption>**Figure 18.** Alignment between the initial angular momentum and angular momentum of the total mass (top panels), angular momentum of stellar particles (middle) and the minor axis of the moment of inertia of the stellar-mass distribution (bottom) as functions of the magnitude of the initial angular momentum (left column) and spin parameter $\lambda'$ (middle column). Histograms are shown on the right. Galaxies with $M_* > 2 \times 10^9 M_\odot$ at $z=6$ are used. The initial angular momentum is strongly aligned with the total angular momentum of galaxies.</figcaption></figure>

This finding, in broad agreement with previous studies, raises many questions: (i) Why should the default initial morphology be a disk? (ii) What makes the fraction of the transient morphology (i.e., spheroid or irregular type) a particular value of about 1/3? (iii) Why does the fraction not change much with time or stellar mass? And (iv) why do merger events not always preferentially result in spheroidal morphology, as found for the massive galaxies at low redshifts? The inspection of the relations between the galaxy morphology and the initial conditions leads us to raise more questions such as (v) when does the environmental dependence of galaxy morphology appear and (vi) when is the correlation between morphology and kinematics (re-)established?

While we cannot answer all of these questions in view of this preliminary work on HR5, let us address specifically a few of them here. One key ingredient seems to be the consistency of gas-rich inflow, set up by the larger cosmic scales, via a top-down causation. Galaxy morphology is disk dominated but is poorly correlated with the initial angular momentum (AM) as seen in Figures 17 and 18. The initial AM produced by the tidal torque on protogalactic regions tends to make newly born galaxies disk like (Navarro et al. 2004; Bailin & Steinmetz 2005; Pichon et al. 2011; Codis et al. 2012), making disks the dominant morphology in the cosmic morning. Indeed, the dynamics within neighboring rhizomes drives cold gas out of the surrounding voids into accreting filaments. As boundaries of asymmetric voids, they acquire a net transverse motion, which explains the angular-momentum-rich nature of infall. This large-scale driven consistency explains why cold streams are so efficient at building up cosmic morning disks. The tight correlations of the AM and spin parameter of the

stellar component of galaxies with the initial AM found in this simulation support this interpretation.

It is also expected that structure formation and evolution on a galactic scale have been largely controlled by the large-scale power of the tidal torque, which has played an important role in setting up the preferred direction of the inflow (Sugerman et al. 2000; Lee & Pen 2000; Aubert et al. 2004; Bailin & Steinmetz 2005). The total mass of galaxies having $M_* = 1-2 \times 10^9 M_\odot$ at $z=6$ is typically $M_{\text{tot}} = 1 \times 10^{11} M_\odot$. The radius of top-hat sphere containing this much of mass in the initial condition is $R_{\text{TH}} = 0.849 \, \text{cMpc}$ in the cosmology adopted in this paper. The slope of our $\Lambda$CDM power spectrum is $d \ln P/d \ln k = -2.63$ at $k = 2\pi/R_{\text{TH}}$. When the slope of power spectrum is close to $d \ln P/d \ln k = -3$, the structure formation proceeds almost at the same time from subgalactic to supergalactic scales, and the influence of large-scale modes on galaxy formation should therefore be significant, leading to a top-down impact of the larger scale on galactic infall (Zel'dovich 1970a) (see also Codis et al. 2012, 2015; Laigle et al. 2015, for a complementary geometric formulation in terms of the size of vorticity and AM caustics and the mass of nonlinearity at that redshift). The situation is quite different for galaxy clusters where the slope of power spectrum at the relevant scale is about -1. The first nonlinear “galactic-scale” structures appearing at cosmic dawn ($z \gtrsim 10$) are isolated spots corresponding to the highest peaks in the initial density field, smoothed over a few hundred ckpc. Nearly simultaneously to the collapse of the peaks, the surrounding matter also collapses to form a few nonlinear filaments attached to them. Rapid mass accretion occurs onto peaks and filaments (along the cold flows; Kereš et al. 2005; Dekel & Birnboim 2006; Ocvirk et al. 2008).

Figure 19 shows zoomed snapshots of a relatively high-density region with 1.5 cMpc thickness in HR5 where some first galaxies and filaments are forming. The left column shows gas density and stars, and right column shows dark-matter particles. Velocity vectors projected on the slice are shown where gas density is relatively high. The spin direction of galaxy stellar-mass projected on the slice is shown with thick arrows, and the sense of the spin perpendicular to the slice is marked with encircled arrows. At most places, mass-flow patterns are very similar between gas and dark matter, except in the interior of dense gas clumps (see also Codis et al. 2012, their Figure 10). On large scales, the velocity field is determined by dark matter, the most-massive component of the universe, and gas is just following the dark-matter velocity field. Gas and stellar-mass flows deviate from the dark-matter flow only within galaxies.

In the particular region of HR5 shown in Figure 19, a few nonlinear overdensity peaks and filament segments of gas and dark matter appear at $z \gtrsim 8$. When they first form, the protogalactic density peaks are elongated in alignment with the filaments that they are embedded in. During the period $8 \gtrsim z \gtrsim 6$, however, they turn around and their disks become roughly perpendicular to the filaments. As a result the filament is wrapped around the forming galaxies. Meanwhile, the spin axes of galaxies are roughly aligned with their filament (see the spin directions within and perpendicular to the slice). This phenomenon seems due to the tidal torque exerted to the protogalaxies by the large-scale gravitational shear field and the corresponding cold gas flow (Kimm et al. 2011; Tillson et al. 2015). This indicates that the galaxies with $M_* \sim 10^9 M_\odot$ in the cosmic morning at $z=6$ may have acquired disk morphology through asymmetric mass infall onto filaments carrying the initial AM.
<figure><img src="image_21.png" /><figcaption>**Figure 19.** An 1.5 cMpc-thick slice of a relatively high-density region in HR5 showing gas density and stars (left) and dark matter (right) from $z = 10$ to 4. Gas density is higher in redder shades. Whitish gray dots in the left illustrate stellar-mass distribution. The white arrows show the velocity field of gas in dense areas. The black arrows at dense cores mark the spin direction of galaxy stellar-mass projected on the slice. The spin direction perpendicular to the slice (clockwise or counterclockwise rotation) is also shown on top of the spin arrows. Numbers are logarithmic stellar mass (left) and dark-matter halo mass (right) in $M_{\odot}$ units. The areas where dark-matter velocity field is shown are the same as in the gas-density plot.</figcaption></figure>
After $z=6$, these first galaxies can undergo mass accretion and mergers along filaments and their spins can swing (e.g., Bailin & Steinmetz 2005; Aragón-Calvo et al. 2007; Hahn et al. 2007), which weakens the alignment of their spin with the initial AM (Bett & Frenk 2012; Dubois et al. 2014b; González et al. 2017; Wang & Kang 2018). Filamentary segments then merge together (Cadiou et al. 2020) to build thicker filaments or collapse by themselves and form kinematically and morphologically more complex systems of filaments beaded with galaxies ($6 \gtrsim z \gtrsim 4$).

The fact that the total AM of galaxies has a very tight correlation with their initial AM supports the above scenario and can also explain why irregulars or spheroids tend to transform back to disks (Welker et al. 2014). Even though their stellar component is temporarily perturbed due to encounters and mergers, the matter surrounding them still traces the initial AM and can be accreted consistently along vorticity-rich cold flows, in such a way that the initial disk morphology is recovered. The correlation, however, gradually weakens with decreasing redshifts. At $z=4.5$, the fraction of galaxies in $\hat{L} \cdot \hat{L}_{\text{tot}} > 0.8$ is only one-third the fraction at $z=6$.

In the meantime, the morphology of the stellar component of galaxies seems quite susceptible to transformation into irregular or spheroidal morphology. This makes the initial correlation between internal kinematics and morphology weaker with cosmic time, and the correlation between the initial AM and galaxy orientation becomes almost negligible by $z=6$ (see also Danovich et al. 2015; Cadiou et al. 2022). It remains for us to explain the fraction of temporary morphology (i.e., irregular plus spheroid types) of about 1/3 at mass scales $M_* \lesssim 10^{10} M_\odot$ in the cosmological model adopted here.

Let us return to the near independence of morphological type fractions from stellar mass. This finding is true only for the stellar-mass range between 2 and $\sim 10 \times 10^9 M_\odot$ and before the epoch $z \sim 4$. Statistics is not good at higher masses. Even though we focus on galaxy morphology during the cosmic morning before $z=4$, we have also calculated the fraction of morphological types near the cosmic high noon at $z \approx 2$ to get a hint on what happens to galaxy morphology at later epochs.

Figure 20 shows that the fraction of each morphology is nearly the same as in the cosmic morning for galaxies with $M_* \lesssim 10^{10} M_\odot$ but is radically different at higher mass. In particular, the fraction of spheroids below $M_* \approx 5 \times 10^{10} M_\odot$ is consistent with the fraction of the GOODS early-type galaxies at $z \approx 1$ reported by Hwang & Park (2009). The absolute magnitude limit of $M_B = -18.0$ of the GOODS galaxies used in the analysis corresponds to the stellar mass of $M_* \sim 3 \times 10^{10} M_\odot$. This agreement on galaxy morphology fraction at relatively low redshifts between the observation and our simulation provides the basis for our work on morphology of galaxies in the cosmic morning.

A notable feature is a sharp rise of the irregular fraction from $\sim 1$ to $3 \times 10^{10} M_\odot$ and about 35% of massive galaxies above $3 \times 10^{10} M_\odot$ are irregulars. So, we expect massive irregular types appear and become frequent as we approach cosmic noon from cosmic morning (see also, e.g., Dubois et al. 2016).

Another dramatic feature in Figure 20 is another steep rise of the spheroid fraction at $M_* \sim 5 \times 10^{10} M_\odot$. This is probably related to the phenomenon that the irregular fraction starts to rise at smaller mass of $M_* \sim 1 \times 10^{10} M_\odot$. One can interpret the mass scale difference as a mass-dependent morphology transformation process from disk morphology via irregular

<figure><img src="image_22.png" /><figcaption>**Figure 20.** Morphological type fractions of galaxies as a function of galaxy stellar mass at $z \approx 2$. The blue, black, and red solid lines are for disks, irregulars, and spheroids, respectively.</figcaption></figure>

type to spheroid type as galaxies grow in mass. Above the mass scale of $1 \times 10^{11} M_\odot$, disk morphology is now a minor type with its fraction of only about 26%. In a forthcoming paper, we will inspect how the fraction of each galaxy morphological type as a function of stellar mass evolves during the cosmic noon (from $z=4$ to $1.5$ as defined by this paper) and afternoon ($z < 1.5$).

## 6. Conclusions

Using our new cosmological simulation, HR5, we find that the initial physical morphology of galaxies in the cosmic morning, before $z=4$, is disk dominant with the Sérsic index of stellar-mass distribution less than 1.5 for about two-thirds of galaxies. The fraction of irregulars or spheroids is about one-sixth each. The fractions are nearly independent of redshift and also of stellar mass up to $\sim 10^{10} M_\odot$. Stellar-mass jump events in the cosmic morning do not preferentially drive galaxy morphology toward spheroid shape at the mass scales $M_* \lesssim 10^{10} M_\odot$. However, the galaxies passing through irregular-type phase have a weak tendency to end up having spheroid morphology.

Almost all the first galaxies in the cosmic morning form at the initial matter-density peaks. Higher peaks form galaxies at higher redshifts. For example, the galaxies having $M_* = 2.0-2.4 \times 10^9 M_\odot$ at a redshift $z$ correspond to the peaks with mean height $\bar{v}_{pk} = 0.472z + 1.32 \pm 0.01$ in the initial overdensity field smoothed over $R_G = 0.35$ cMpc. That is, galaxies with $M_* \approx 2 \times 10^9 M_\odot$ at $z=6$ have formed at initially $\sim 4.1\sigma$ density peaks. More massive galaxies at the same epoch are the ones that have formed earlier at higher peaks and grown in mass.

Large-scale structures in the universe emerge and grow like cosmic rhizomes as the underlying matter-density fluctuations become nonlinear and form associations of galaxies in globally overdense regions and stretch the realm of the world of galaxies into lower-density regions along multiply connected filaments. The cosmic-web structure of the large-scale distribution traced
by galaxies forms at lower redshifts when most rhizomes percolate.

The initial angular momentum produced by the tidal torque on protogalactic regions is tightly correlated with the internal kinematics and spin parameters of galaxies, confirming the validity of tidal torque theory. However, the actual shape of the stellar component of galaxies has almost no correlation with the initial angular momentum. These findings can be reconciled if the initial morphology of all galaxies is disk like, while the morphology of the stellar component of the first galaxies is set up by the anisotropy of the larger-scale tides, which induce coherent cosmic gas inflow.

The validity of our predictions, of course, depends on how accurately HR5 matches the observed universe. Because we have studied the physical morphology of galaxies using the stellar-mass density, our results cannot be directly compared with the apparent morphology from observations that rely on light distribution. Nonetheless, we hope that many of our findings be judged in view of upcoming observations of the deep universe by the James Webb Space Telescope.

We thank Katarina Kraljic for helpful comments. C.P. and J.K. are supported by KIAS Individual Grants (PG016903, KG039603) at Korea Institute for Advanced Study. J.L. is supported by the National Research Foundation of Korea (NRF-2021R1C1C2011626). J.S. acknowledges support from the National Research Foundation of Korea grant (2021R1C1C1003785) funded by the Ministry of Science, ICT & Future Planning. B.K.G. and C.G.F. acknowledge the support of STFC through the University of Hull Consolidated Grant ST/ R000840/1, access to VIPER, the University of Hull High Performance Computing Facility, and the European Unions Horizon 2020 research and innovation program (ChETEC-INFRA—Project No. 101008324). This work benefited from the outstanding support provided by the KISTI National Supercomputing Center and its Nurion Supercomputer through the Grand Challenge Program (KSC-2018-CHA-0003). This research was also partially supported by the ANR-19-CE31-0017 (http://www.secular-evolution.org). Large data transfer was supported by KREONET, which is managed and operated by KISTI. This work is supported by the Center for Advanced Computation at Korea Institute for Advanced Study.

# Appendix
## Sérsic Index and Morphology of Galaxies

We classify galaxy morphology into disk and spheroid types for symmetric galaxies according to the Sérsic index $n_{Sérsic}$. We determine the Sérsic index criterion from the morphology training set of SDSS galaxies given by Park & Choi (2005). The sample contains visually identified morphology information of 1982 nearby bright late-type (spirals and irregulars) and early-type (ellipticals and lenticulars) galaxies. We measure the radial surface brightness profile using an ellipse-fitting method to account for the inclination effects and find the best-fit $n_{Sérsic}$ between 0.2 and 2 times the Petrosian radius (Park & Choi 2005; Park et al. 2007). The central region has been discarded to mitigate the effects of the finite point-spread function of the SDSS survey (about 1"/4). Figure 21 shows the isophotal axis ratios and Sérsic indices of late-type galaxies (green circles and triangles) and early-type SDSS galaxies (red/magenta circles and crosses) in the set. The distribution of $n_{Sérsic}$ is continuous, and the transition from late to early type is

<figure><img src="image_23.png" /><figcaption>**Figure 21.** Relation between the isophotal axis ratio and Sérsic index of late-type (green circles), early-type (red or magenta circles), abnormal late-type (triangles), and abnormal early-type (crosses) SDSS galaxies in the training set of Park & Choi (2005).</figcaption></figure>

not abrupt. However, it can be seen that when galaxies are classified by using $n_{Sérsic}$ alone, the threshold of $n_{Sérsic} = 1.5$ is a good division criterion.

## ORCID iDs

Changbom Park https://orcid.org/0000-0001-9521-6397

Jaehyun Lee https://orcid.org/0000-0002-6810-1778

Juhan Kim https://orcid.org/0000-0002-4391-2275

Donghui Jeong https://orcid.org/0000-0002-8434-979X

Christophe Pichon https://orcid.org/0000-0003-0695-6735

Brad K. Gibson https://orcid.org/0000-0003-4446-3130

Jihye Shin https://orcid.org/0000-0001-5135-1693

Yohan Dubois https://orcid.org/0000-0003-0225-6387

C. Gareth Few https://orcid.org/0000-0003-3940-7687

## References

* Abraham, R. G., & van den Bergh, S. 2001, *Sci*, 293, 1273

* Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, *A&A*, 594, A13

* Aragón-Calvo, M. A., van de Weygaert, R., Jones, B. J. T., & van der Hulst, J. M. 2007, *ApJ*, 655, L5

* Aubert, D., Pichon, C., & Colombi, S. 2004, *MNRAS*, 352, 376

* Bailin, J., & Steinmetz, M. 2005, *ApJ*, 627, 647

* Bardeen, J. M., Bond, J. R., Kaiser, N., & Szalay, A. S. 1986, *ApJ*, 304, 15

* Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013, *ApJ*, 762, 109

* Bett, P. E., & Frenk, C. S. 2012, *MNRAS*, 420, 3324

* Bond, J. R., Cole, S., Efstathiou, G., & Kaiser, N. 1991, *ApJ*, 379, 440

* Bond, J. R., Kofman, L., & Pogosyan, D. 1996, *Natur*, 380, 603

* Buitrago, F., Trujillo, I., Conselice, C. J., & Häußler, B. 2013, *MNRAS*, 428, 1460

* Bullock, J. S., Kolatt, T. S., Sigad, Y., et al. 2001, *MNRAS*, 321, 559

* Bundy, K., Ellis, R. S., & Conselice, C. J. 2005, *ApJ*, 625, 621

* Cadiou, C., Dubois, Y., & Pichon, C. 2022, *MNRAS*, 514, 5429

* Cadiou, C., Pichon, C., Codis, S., et al. 2020, *MNRAS*, 496, 4787

* Cameron, E., Carollo, C. M., Oesch, P. A., et al. 2011, *ApJ*, 743, 146

* Chabrier, G. 2003, *ApJL*, 586, L133

* Codis, S., Pichon, C., Devriendt, J., et al. 2012, *MNRAS*, 427, 3320

* Codis, S., Pichon, C., & Pogosyan, D. 2015, *MNRAS*, 452, 3369

* Conselice, C. J., Bershady, M. A., & Jangren, A. 2000, *ApJ*, 529, 886

* Dahlen, T., Mobasher, B., Dickinson, M., et al. 2007, *ApJ*, 654, 172

* Dalgarno, A., & McCray, R. A. 1972, *ARA&A*, 10, 375
Danovich, M., Dekel, A., Hahn, O., Ceverino, D., & Primack, J. 2015, MNRAS, 449, 2087

Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827

Dekel, A., & Birnboim, Y. 2006, MNRAS, 368, 2

Dekel, A., Birnboim, Y., Engel, G., et al. 2009, Natur, 457, 451

Dolag, K., Komatsu, E., & Sunyaev, R. 2016, MNRAS, 463, 1797

Doroshkevich, A. G. 1970, ApJ, 6, 320

Dressler, A. 1980, ApJ, 236, 351

Dubois, Y., Devriendt, J., Slyz, A., & Teyssier, R. 2012, MNRAS, 420, 2662

Dubois, Y., Peirani, S., Pichon, C., et al. 2016, MNRAS, 463, 3948

Dubois, Y., & Teyssier, R. 2008, A&A, 477, 79

Dubois, Y., Volonteri, M., & Silk, J. 2014a, MNRAS, 440, 1590

Dubois, Y., Pichon, C., Welker, C., et al. 2014b, MNRAS, 444, 1453

Elahi, P. J., Cañas, R., Poulton, R. J. J., et al. 2019, PASA, 36, e021

Elmegreen, D. M. 2015, Lessons from the Local Group: A Conf. in Honor of David Block and Bruce Elmegreen (Cham: Springer), 455

Elmegreen, D. M., & Elmegreen, B. G. 2014, ApJ, 781, 11

Few, C. G., Courty, S., Gibson, B. K., et al. 2012, MNRAS, 424, L11

Franceschini, A., Rodighiero, G., Cassata, P., et al. 2006, A&A, 453, 397

Fry, J. N. 1996, ApJL, 461, L65

Genel, S., Vogelsberger, M., Springel, V., et al. 2014, MNRAS, 445, 175

González, R. E., Prieto, J., Padilla, N., & Jimenez, R. 2017, MNRAS, 464, 4666

Gott, J. R., Melott, A. L., & Dickinson, M. 1986, ApJ, 306, 341

Guillet, T., & Teyssier, R. 2011, JCoPh, 230, 4756

Haardt, F., & Madau, P. 1996, ApJ, 461, 20

Hahn, O., & Abel, T. 2011, MNRAS, 415, 2101

Hahn, O., Porciani, C., Carollo, C. M., & Dekel, A. 2007, MNRAS, 375, 489

Haywood, M., Di Matteo, P., Lehnert, M. D., Katz, D., & Gómez, A. 2013, A&A, 560, A109

Hong, S. E., Park, C., & Kim, J. 2016, ApJ, 823, 103

Hubble, E. P. 1926, ApJ, 64, 321

Hwang, H. S., & Park, C. 2009, ApJ, 700, 791

Jung, I., Lee, J., & Yi, S. K. 2014, ApJ, 794, 74

Kaiser, N. 1984, ApJL, 284, L9

Kelvin, L. S., Driver, S. P., Robotham, A. S. G., et al. 2014, MNRAS, 444, 1647

Kereš, D., Katz, N., Weinberg, D. H., & Davé, R. 2005, MNRAS, 363, 2

Khandai, N., Di Matteo, T., Croft, R., et al. 2015, MNRAS, 450, 1349

Kim, J., & Park, C. 2006, ApJ, 639, 600

Kim, J., Park, C., Gott, J. R., & Dubinski, J. 2009, ApJ, 701, 1547

Kim, J., Park, C., L'Huillier, B., & Hong, S. E. 2015, JKAS, 48, 213

Kimm, T., Devriendt, J., Slyz, A., et al. 2011, arXiv:1106.0538

Knebe, A., Knollmann, S. R., Muldrew, S. I., et al. 2011, MNRAS, 415, 2293

Knebe, A., & Power, C. 2008, ApJ, 678, 621

Knollmann, S. R., & Knebe, A. 2009, ApJS, 182, 608

Laigle, C., Pichon, C., Codis, S., et al. 2015, MNRAS, 446, 2744

L'Huillier, B., Park, C., & Kim, J. 2014, NewA, 30, 79

Lee, J., & Pen, U.-L. 2000, ApJL, 532, L5

Lee, J., Shin, J., Snaith, O. N., et al. 2021, ApJ, 908, 11

Lee, J., Yi, S. K., Elahi, P. J., et al. 2014, MNRAS, 445, 4197

Lewis, A., Challinor, A., & Lasenby, A. 2000, ApJ, 538, 473

Lotz, J. M., Madau, P., Giavalisco, M., Primack, J., & Ferguson, H. C. 2006, ApJ, 636, 592

Marinacci, F., Vogelsberger, M., Pakmor, R., et al. 2018, MNRAS, 480, 5113

Marshall, M. A., Mutch, S. J., Qin, Y., Poole, G. B., & Wyithe, J. S. B. 2019, MNRAS, 488, 1941

Martin, G., Kaviraj, S., Devriendt, J. E. G., Dubois, Y., & Pichon, C. 2018, MNRAS, 480, 2266

McCarthy, I. G., Schaye, J., Bird, S., & Le Brun, A. M. C. 2017, MNRAS, 465, 2936

Musso, M., Cadiou, C., Pichon, C., et al. 2018, MNRAS, 476, 4877

Navarro, J. F., Abadi, M. G., & Steinmetz, M. 2004, ApJL, 613, L41

Ocvirk, P., Pichon, C., & Teyssier, R. 2008, MNRAS, 390, 1326

Oemler, A. J. 1974, ApJ, 194, 1

Papovich, C., Dickinson, M., Giavalisco, M., Conzelice, C. J., & Ferguson, H. C. 2005, ApJ, 631, 101

Park, C. 1990, MNRAS, 242, 59P

Park, C., & Choi, Y.-Y. 2005, ApJL, 635, L29

Park, C., Choi, Y.-Y., Kim, J., et al. 2012, ApJL, 759, L7

Park, C., Choi, Y.-Y., Vogeley, M. S., et al. 2007, ApJ, 658, 898

Park, M.-J., Yi, S. K., Dubois, Y., et al. 2019, ApJ, 883, 25

Pichon, C., Pogosyan, D., Kimm, T., et al. 2011, MNRAS, 418, 2493

Pillepich, A., Nelson, D., Springel, V., et al. 2019, MNRAS, 490, 3196

Rasera, Y., & Teyssier, R. 2006, A&A, 445, 1

Ravindranath, S., Giavalisco, M., Ferguson, H. C., et al. 2006, ApJ, 652, 963

Ricciardelli, E., Vazdekis, A., Cenarro, A. J., & Falcón-Barroso, J. 2012, MNRAS, 424, 172

Schäfer, B. M., & Merkel, P. M. 2012, MNRAS, 421, 2751

Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521

Scoccimarro, R. 1998, MNRAS, 299, 1097

Sérsic, J. L. 1963, BAAA, 6, 41

Shandarin, S. F., & Zeldovich, Y. B. 1989, RvMP, 61, 185

Springel, V., White, S. D. M., Tormen, G., & Kauffmann, G. 2001, MNRAS, 328, 726

Srisawat, C., Knebe, A., Pearce, F. R., et al. 2013, MNRAS, 436, 150

Stewart, K. R., Brooks, A. M., Bullock, J. S., et al. 2013, ApJ, 769, 74

Sugerman, B., Summers, F. J., & Kamionkowski, M. 2000, MNRAS, 311, 762

Sutherland, R. S., & Dopita, M. A. 1993, ApJS, 88, 253

Tamburri, S., Saracco, P., Longhetti, M., et al. 2014, A&A, 570, A102

Tegmark, M., & Peebles, P. J. E. 1998, ApJL, 500, L79

Teyssier, R. 2002, A&A, 385, 337

Tillson, H., Devriendt, J., Slyz, A., Miller, L., & Pichon, C. 2015, MNRAS, 449, 4363

Toro, E. F., Spruce, M., & Speares, W. 1994, ShWav, 4, 25

Trayford, J. W., Frenk, C. S., Theuns, T., Schaye, J., & Correa, C. 2019, MNRAS, 483, 744

Tremmel, M., Karcher, M., Governato, F., et al. 2017, MNRAS, 470, 1121

Vazdekis, A., Koleva, M., Ricciardelli, E., Röck, B., & Falcón-Barroso, J. 2016, MNRAS, 463, 3409

Vazdekis, A., Ricciardelli, E., Cenarro, A. J., et al. 2012, MNRAS, 424, 157

Vogeley, M. S., Park, C., Geller, M. J., Huchra, J. P., & Gott, J. R. 1994, ApJ, 420, 525

Wang, P., & Kang, X. 2018, MNRAS, 473, 1562

Welker, C., Devriendt, J., Dubois, Y., Pichon, C., & Peirani, S. 2014, MNRAS, 445, L46

White, S. D. M. 1984, ApJ, 286, 38

Zel'dovich, Y. B. 1970a, A&A, 5, 84

Zel'dovich, Y. B. 1970b, A&A, 500, 13