# Young and turbulent: the early life of massive galaxy progenitors

Davide Fiacconi, Lucio Mayer, Piero Madau, Alessandro Lupi, Massimo Dotti, Francesco Haardt

► To cite this version:

Davide Fiacconi, Lucio Mayer, Piero Madau, Alessandro Lupi, Massimo Dotti, et al. Young and turbulent: the early life of massive galaxy progenitors. Monthly Notices of the Royal Astronomical Society, 2017, 467 (4), pp.4080-4100. 10.1093/mnras/stx335. hal-01554765

HAL Id: hal-01554765

https://hal.science/hal-01554765v1

Submitted on 29 Aug 2025

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
<figure><img src="image_2.png" /></figure>

# Young and turbulent: the early life of massive galaxy progenitors

Davide Fiacconi,<sup>1,2★</sup> Lucio Mayer,<sup>2</sup> Piero Madau,<sup>2,3,4</sup> Alessandro Lupi,<sup>4</sup> Massimo Dotti<sup>5,6</sup> and Francesco Haardt<sup>6,7</sup>

Accepted 2017 February 6. Received 2017 February 6; in original form 2016 September 26

## ABSTRACT

We present results from the ‘Ponos’ simulation suite on the early evolution of a massive, $M_{vir}(z=0) = 1.2 \times 10^{13}$ M$_{\odot}$ galaxy. At $z \gtrsim 6$, before feedback from a central supermassive black hole becomes dominant, the main galaxy has a stellar mass $\sim 2 \times 10^9$ M$_{\odot}$ and a star formation rate $\sim 20$ M$_{\odot}$ yr$^{-1}$. The galaxy sits near the expected main sequence of star-forming galaxies at those redshifts, and resembles moderately star-forming systems observed at $z > 5$. The high specific star formation rate results in vigorous heating and stirring of the gas by supernovae feedback, and the galaxy develops a thick and turbulent disc, with gas velocity dispersion $\sim 40$ km s$^{-1}$, rotation to dispersion ratio $\sim 2$, and with a significant amount of gas at $\sim 10^5$ K. The Toomre parameter always exceeds the critical value for gravito-turbulence, $Q \sim 1.5–2$, mainly due to the contribution of warm/hot gas inside the disc. Without feedback, a nearly gravito-turbulent regime establishes with similar gas velocity dispersion and lower $Q$. We propose that the ‘hot and turbulent’ disc regime seen in our simulations, unlike the ‘cold and turbulent’ gravito-turbulent regime of massive clumpy disc galaxies at $z \sim 1–2$, is a fundamental characterization of the main-sequence galaxies at $z \gtrsim 6$, as they can sustain star formation rates comparable to those of low-mass starbursts at $z=0$. This results in no sustained coherent gas inflows through the disc, and in fluctuating and anisotropic mass transport, possibly postponing the assembly of the bulge and causing the initial feeding of the central black hole to be highly intermittent.

**Key words:** turbulence – methods: numerical – galaxies: formation – galaxies: high-redshift – galaxies: ISM.

## 1 INTRODUCTION

Since the early 1990s, the growing number of galaxy surveys exploring different wavelengths and cosmological volumes have provided an immense body of data to reconstruct the formation and evolution of galaxies (e.g. York et al. 2000; Steidel et al. 2003; Le Fèvre et al. 2005; Lilly et al. 2007; Walter et al. 2008; Baldry et al. 2010; Grogin et al. 2011; Koekemoer et al. 2011; Kochanek et al. 2012). At the same time, the Λ cold dark matter (ΛCDM) cosmological paradigm of structure formation has produced detailed predictions on the statistical matter distribution in the Universe (e.g. Percival et al. 2001; Springel et al. 2005; Reed et al. 2007; Viel,

Bolton & Haehnelt 2009; Percival et al. 2010; Reid et al. 2010; Klypin, Trujillo-Gomez & Primack 2011), but the theory of galaxy formation within this framework is still developing, through analytical work (e.g. White & Rees 1978; Blumenthal et al. 1984; Mo, Mao & White 1998) and semi-analytic models (e.g. Kauffmann, White & Guideroni 1993; Cole et al. 1994; Somerville & Primack 1999; Springel et al. 2005; Bower et al. 2006; Croton et al. 2006; Guo et al. 2011; Henriques et al. 2013), as well as cosmological numerical simulations (e.g. Diemand et al. 2008; Crain et al. 2009; Schaye et al. 2010; Agertz, Teyssier & Moore 2011; Guedes et al. 2011; Dubois et al. 2014; Hopkins et al. 2014; Vogelsberger et al. 2014; Schaye et al. 2015; and references therein).

Most of the theoretical successes have been achieved at low to moderate redshifts ($z \sim 0–1.5$), where the majority and the highest quality data are available. The latter have been fundamental to

<footnote n="1"><sup>1</sup>Institute of Astronomy and Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="2"><sup>2</sup>Center for Theoretical Astrophysics and Cosmology, Institute for Computational Science, University of Zurich, Winterthurerstrasse 190, CH-8057 Zürich, Switzerland</footnote>

<footnote n="3"><sup>3</sup>Department of Astronomy and Astrophysics, University of California Santa Cruz, 1156 High St, Santa Cruz, CA 95064, USA</footnote>

<footnote n="4"><sup>4</sup>Institut d'Astrophysique de Paris, Sorbonne Universités, UPMC Univ Paris 6 et CNRS, UMR 7095, 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote n="5"><sup>5</sup>Dipartimento di Fisica, Università di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy</footnote>

<footnote n="6"><sup>6</sup>INFN, Milano-Bicocca, Pza della Scienza 3, I-20126 Milano, Italy</footnote>

<footnote n="7"><sup>7</sup>DiSAT, Università degli Studi dell'Insubria, Via Valleggio 11, I-22100 Como, Italy</footnote>

<footnote>*E-mail: fiacconi@ast.cam.ac.uk</footnote>
gauge the characteristics and the parameters of phenomenological, sub-resolution recipes introduced in numerical simulations (as well as in semi-analytic models) in order to model physical processes on scales that are currently inaccessible directly, such as gas radiative cooling, gas chemistry evolution, star formation, stellar feedback, black hole accretion and feedback (e.g. Stinson et al. 2006; Dalla Vecchia & Schaye 2008; Gnedin, Tassis & Kravtsov 2009; Wiersma et al. 2009; Vogelsberger et al. 2013; Keller et al. 2014). Large-box cosmological simulations, despite suffering from relatively low resolution, have been able to statistically reproduce the main properties of different galaxy populations, from small, star-forming irregular galaxies to massive quiescent ellipticals (Vogelsberger et al. 2014; Schaye et al. 2015). On the other hand, high-resolution zoom-in simulations focusing on the formation and growth of a limited number of systems have also satisfactory modelled the formation and growth of dwarf galaxies (e.g. Governato et al. 2010), Milky Way like hosts (e.g. Agertz et al. 2011; Guedes et al. 2011), as well as massive ellipticals at the centre of galaxy groups (e.g. Feldmann et al. 2010; Feldmann & Mayer 2015; Fiacconi, Feldmann & Mayer 2015).

At higher redshifts $z \sim 2$, recent observational campaigns have boosted the interest in the population of massive star-forming galaxies that show peculiar properties when compared to local counterparts of similar size. Such systems are often characterized by massive discs with baryonic mass $\sim 10^{11} M_{\odot}$ and star formation rates as high as $\sim 100 M_{\odot} \text{yr}^{-1}$, with a turbulent interstellar medium that accounts for $\gtrsim 30$ per cent of the baryonic mass (e.g. Genzel et al. 2006; Förster Schreiber et al. 2009; Daddi et al. 2010; Tacconi et al. 2010; Wisnioski et al. 2015). The new observations have triggered lively discussions in the theoretical community, also requiring the development of new theoretical ideas to explain the observations and to make new predictions (e.g. violent disc instability; Mandelker et al. 2014; Inoue et al. 2016).

Pushing it forward, the very high redshifts ($z \gtrsim 4$) still remain a partially unexplored territory, clearly because of the larger technical difficulties involved in new, cutting-edge observations. None the less, galaxies at $z > 4$ have been detected in a few different ways. Optical and near infrared observations have targeted star-forming galaxies by identifying them through the flux dropout in adjacent bands around the Lyman break (e.g. Madau et al. 1996; Steidel et al. 1999; Bouwens et al. 2003; Oesch et al. 2010a). Those have permitted to constrain the evolution of the cosmic star formation rate density as well as of the ultraviolet (UV) luminosity function of star-forming galaxies out to $z \gtrsim 8-10$, showing that the latter has a steep low luminosity tail, when compared to local samples (e.g. Bouwens et al. 2007, 2011; Oesch et al. 2012, 2014; Bouwens et al. 2015). Moreover, other wavelengths have been effective in providing information about the early galaxy population. Sub-millimetre galaxies are an example of massive (stellar mass $\sim 10^{11} M_{\odot}$), highly star-forming (star formation rates $\gtrsim 500 M_{\odot} \text{yr}^{-1}$), dusty galaxies that have been detected mostly at $z > 2.5-3$ (e.g. Chapman et al. 2005; Younger et al. 2009; Casey, Narayanan & Cooray 2014).

While those studies have been important to understand the global properties of the first galaxies in the Universe, they have mostly revealed the luminous tail of the galaxy population and they are still not able to characterize their structure (but see Oesch et al. 2010b; De Breuck et al. 2014). None the less, both available facilities, such as the Hubble Space Telescope (HST) or the Very Large Telescope (VLT), as well as new observatories that recently came online, such as the Atacama Large Millimeter Array (ALMA), are starting to discover smaller and possibly more typical galaxies at $z \ge 6$. For example, Bradley et al. (2012) have used HST imaging to identify a few Lyman break galaxy candidates at $z \approx 7$, lensed by the

foreground galaxy cluster A1703. The most luminous likely has a stellar mass $\sim 10^9 M_{\odot}$ and a star formation rate $\sim 8 M_{\odot} \text{yr}^{-1}$. Similarly, Watson et al. (2015) have combined HST, VLT and ALMA observations to constrain the stellar mass, star formation and gas content of another highly magnified Lyman break galaxy beyond the Bullet cluster. They also find a stellar mass $\sim 2 \times 10^9 M_{\odot}$, a star formation rate $\sim 10 M_{\odot} \text{yr}^{-1}$ and a gas fraction $\sim 40-50$ per cent, all confined within a physical surface $\sim 1.5 \text{ kpc}^2$. All these objects typically have specific star formation rates $\sim 5-10 \text{ Gyr}^{-1}$ (see also e.g. Tasca et al. 2015).

These recent observations, combined with those from the next generation of both space-based (James Webb Space Telescope; JWST) and earth-based (e.g. European Extremely Large Telescope; E-ELT) telescopes, require new interpretations and predictions on the theoretical side, where less has been done compared to low/medium redshifts, except regarding the most massive population of galaxies at $z > 6-8$ (e.g. Choi & Nagamine 2012; O'Shea et al. 2015; Paardekooper, Khochfar & Dalla Vecchia 2015; Ocvirk et al. 2016; Waters et al. 2016). Motivated by this, we investigate the early phases of the formation and the evolution of a galaxy that become a massive elliptical at $z=0$. We focus on the first burst of star formation, prior to the assembly of the central supermassive black hole and before active galactic nuclei feedback becomes dominant. We use a new high-resolution hydrodynamic run that is part of the recent Ponos programme of zoom-in cosmological simulations of massive galaxies (Fiacconi et al. 2016). Our new simulation reproduces the main features of recently observed star-forming galaxies at $z \sim 7$, and allows us to study in details the properties of the early interstellar medium that drive the star formation and may determine the early feeding habits of the central black hole.

The paper is organized as follows. In Section 2, we describe our numerical techniques and the features of the PonosHydro simulation. We describe our main results in Section 3, focusing on the properties of the interstellar medium and the early star formation history of the simulated system at $z \sim 6$. We present our conclusions in Section 4, cautioning the reader about the limitations of our results and discussing the possible implications of our findings. In the following, when not explicitly specified, all lengths and densities are given in physical units.

# 2 METHODS

## 2.1 Initial conditions

We perform and analyse a new simulation that complements the suite of cosmological, zoom-in simulations presented by Fiacconi et al. (2016) and named ‘Ponos’. This is a high-resolution version of the run focusing on the halo originally dubbed as ‘PonosV’. This new version includes hydrodynamics and self-consistent baryonic physics, and we refer to it as ‘PH’ or ‘PonosHydro’ in the following. The target halo evolves in a box of 85.5 comoving Mpc and it reaches a mass $\approx 1.2 \times 10^{13} M_{\odot}$ at $z=0$. We adopt a $\Lambda$CDM cosmology consistent with the results of Wilkinson Microwave Anisotropy Probe 7/9 yr, parametrized by $\Omega_{m,0} = 0.272$, $\Omega_{\Lambda,0} = 0.728$, $\Omega_{b,0} = 0.0455$, $\sigma_8 = 0.807$, $n_s = 0.961$ and $H_0 = 70.2 \text{ km s}^{-1} \text{ Mpc}^{-1}$ (Komatsu et al. 2011; Hinshaw et al. 2013). The original initial conditions (ICs) are part of the Assembling Galaxies of Resolved Anatomy (AGORA)<fnref n="1" /> code-comparison project (Kim et al. 2014).

<footnote n="1">¹ https://sites.google.com/site/santacruzcomparisonproject/home</footnote>
We generate new ICs of the same halo using the MUSIC<fnref n="2" /> code (Hahn & Abel 2011). Our ICs are optimized to follow the growth of the halo until $z = 6$. They consist of a base cube of $128^3$ particles per side starting at $z = 100$, with additional nested levels of refinement to increase the resolution within the Lagrangian region that maps the particles contained within $2.5R_{vir}$ at $z = 6$ on to the ICs. We define the virial radius $R_{vir}$ as the spherical radius that encompasses a mean matter density $\Delta(z)\rho_c(z)$, where $\rho_c(z)$ is the critical density that the Universe needs to be flat, and $\Delta(z)$ is the z-dependent virial overdensity defined by Bryan & Norman (1998). As a consequence, the virial mass is defined as $M_{vir} = 4\pi\Delta(z)\rho_c(z)R_{vir}^3/3$. At $z = 6$, the virial mass of the main galaxy is $\approx 1.5 \times 10^{11} M_{\odot}$. In fact, we use the following procedure to determine the high-resolution region on the ICs (Fiacconi et al. 2016): (i) we run a $128^3$, dark matter only, full-box simulation to identify the main halo at $z = 6$; (ii) we trace the particles within $2.5R_{vir}$ back on to the ICs; (iii) we locally increase the resolution of the ICs by adding one level of refinement within a rectangular box that contains all the identified particles; (iv) we evolve again the ICs until $z = 6$; (v) we repeat (ii) and we add two additional levels of refinement within the convex hull that contains all the identified particles. Every level of refinement increases the spatial and mass resolution by a factor 2 and 8, respectively. We iterate steps (iv) and (v) until we add seven additional levels of refinement above the base cube, introducing gas particles in the last level. We verified with dark matter only test runs and in the main run that this procedure allows us to maintain the fraction of contaminating particles<fnref n="3" /> in both mass and number well below 0.1 per cent within the virial volume at all $z \ge 6$. Finally, the highest resolution dark matter and gas particle have a mass $m_{dm} = 4397.6 M_{\odot}$ and $m_g = 883.4 M_{\odot}$, respectively. The force resolution is determined by the softening length, which is set to 1/60 of the mean particle separation at each level of refinement. This corresponds to a minimum dark matter softening $\epsilon_{dm} = 81.8$ physical pc and to a gas softening $\epsilon_g = 47.9$ physical pc. The softening is kept constant in physical units during evolution at $z < 9$, while it remains constant in comoving coordinates at higher redshifts. The total number of particles in the ICs is 118 694 002, while the total number of particles within the virial radius at $z = 6.5$ is 56 213 155.

## 2.2 Simulation code

We evolve the simulation using the GASOLINE code (Wadsley, Stadel & Quinn 2004). The code computes gravitational interactions using a fourth-order (i.e. hexadecapole) multipole expansion of the gravitational force on a k-D binary tree, following the original scheme adopted by the PKDGRAV code (Stadel 2001, 2013). GASOLINE models the gas dynamics using the smoothed particle hydrodynamics (SPH) algorithm (e.g. Gingold & Monaghan 1977; Lucy 1977; Hernquist & Katz 1989). In addition to the standard SPH formulation, the energy equation includes a term for thermal energy and metal diffusion (with a coefficient $C = 0.05$) as introduced and discussed by Wadsley, Veeravalli & Couchman (2008) and Shen, Wadsley & Stinson (2010). This approach reduces the artificial surface tension that arises close to strong density gradients where Kelvin-Helmholtz instability may develop (e.g. Agertz et al. 2007; Read & Hayfield 2012; Hopkins 2013; Keller et al. 2014).

GASOLINE includes several sub-resolution models to treat the radiative cooling of the gas, the formation of stars and the impact of their feedback in terms of supernovae-injected energy and mass released by winds. The gas is allowed to cool radiatively in the optically thin limit by solving the non-equilibrium reaction network of H I, H II, He I, He II and He III. We take into account the contribution from metal lines with a temperature floor $T_{floor} \approx 100$ K following Shen et al. (2013). The simulation also includes a uniform, redshift-dependent UV radiation background due to stellar and quasar reionization according to Haardt & Madau (2012). The implementation of star formation and stellar feedback mostly follows the prescriptions from Stinson et al. (2006). Stars form when: (i) the local density exceeds the threshold $n_{SF} = 10 H \text{ cm}^{-3}$; (ii) the local temperature goes below $T_{SF} = 10^4$ K; (iii) the local gas overdensity is $> 2.64$; and (iv) the flow is convergent and locally Jeans unstable. When the above criteria are fulfilled, stars form according to a Schmidt-like low, $\dot{\rho}_* = \epsilon_{SF}\rho_g/t_{dyn} \propto \rho_g^{3/2}$, where $\rho_*$ is the mass density of formed stars, $\rho_g$ is the local gas mass density and $t_{dyn} = 1/\sqrt{G\rho_g}$ is the local dynamical time. The efficiency $\epsilon_{SF} = 0.05$ is a phenomenological parameter meant to capture the average star formation efficiency. Each stellar particle has an initial mass $m_* = 0.4m_g = 353.4 M_{\odot}$ and represents a stellar population with a Kroupa (2001) initial mass function.

We choose $n_{SF}$ to be about the density reached when resolving the local Jeans mass with at least one kernel at the lowest temperature reached by the simulation, i.e. the temperature floor of the cooling function, namely

$$n_{SF} \lesssim 16.4 \left( \frac{T_{\text{floor}}}{100 \text{ K}} \right)^3 \left( \frac{N}{32} \right)^{-2} \left( \frac{m_g}{883.4 \text{ M}_{\odot}} \right)^{-2} \text{H cm}^{-3}, \quad (1)$$

where $N$ is the number of gas particles per kernel. This choice ensures that, even in the most extreme case (i.e. gas at the lowest temperature that has not formed stars yet), the local collapse is resolved. However, since $T_{SF} \gg T_{\text{floor}}$, the Jeans mass during collapses that lead to star formation episodes is likely resolved with a number of particles $\gg N$. Indeed, at $n_{SF}$ and $T_{SF}$, the local Jeans mass is resolved with $\sim 1220$ kernel masses, which corresponds to $\sim 39000$ particles in the regions at the verge of the collapse that might form stars.

During their lifetime, light stars with mass between 1 and $8 M_{\odot}$ release $\sim 40$ per cent of the particle mass and metals to the surrounding gas through winds following Weidemann (1987). On the other hand, massive stars between 8 and $40 M_{\odot}$ are responsible for Type II supernovae, which inject energy, mass and metals to the surrounding gas following the analytical blast wave solution of McKee & Ostriker (1977). In particular, each Type II supernova releases $10^{51}$ erg in thermal energy to the gas particles within the maximum radius that the blast wave can reach. The cooling of those gas particles is temporarily turned off for the time corresponding to the snowplough phase of the blast wave. Type Ia supernovae also contribute to the supernova feedback budget injecting the same energy as Type II supernovae, but without the shut off of the cooling. They releases a mass of $1.4 M_{\odot}$, where 0.63 and 0.13 $M_{\odot}$ are constituted by iron and oxygen, respectively. Their frequency is given by the binary fraction estimated by Raiteri, Villata & Navarro (1996).

We also adopt the pressure floor described by e.g. Roškar et al. (2015) in order to avoid spurious fragmentation. Specifically, the minimum pressure of each gas particle (at a given density and temperature) is $P_{min} = \alpha \max(\epsilon_g, h^2 G\rho^2)$, where $h$ is the local SPH smoothing length, $\rho$ is the gas mass density and $\alpha = 9$ is a safety factor. This choice of $\alpha$ is such that the local Jeans length

<footnote n="2">² http://www.phys.ethz.ch/hahn/MUSIC/</footnote>

<footnote n="3">³ Here, the term 'contaminating particles' refers to dark matter particles from coarser levels of refinement in the ICs.</footnote>
<table><caption>**Table 1.** List of performed simulations and of their main features. From left to right: label of the simulation, initial redshift of the simulation, final redshift of the simulation, $\alpha$ parameter of the pressure floor (see the text), initial mass of gas particles, mass of the lightest dark matter particles, gravitational softening of the gas and smallest gravitational softening of the dark matter.</caption><thead><tr><th>Label</th><th>z<sub>ini</sub></th><th>z<sub>end</sub></th><th>α</th><th>m<sub>g</sub><br>(M<sub>☉</sub>)</th><th>m<sub>dm</sub><br>(M<sub>☉</sub>)</th><th>ε<sub>g</sub><br>(pc)</th><th>ε<sub>dm</sub><br>(pc)</th></tr></thead><tbody><tr><td>PH</td><td>100</td><td>6.5</td><td>9</td><td>883.4</td><td>4397.6</td><td>47.9</td><td>81.8</td></tr><tr><td>PH_PF1</td><td>8</td><td>6.5</td><td>1</td><td>883.4</td><td>4397.6</td><td>47.9</td><td>81.8</td></tr><tr><td>PH_NF<sup>a</sup></td><td>8</td><td>7.1</td><td>9</td><td>883.4</td><td>4397.6</td><td>47.9</td><td>81.8</td></tr><tr><td>PH_LR</td><td>100</td><td>6.5</td><td>9</td><td>7067.2</td><td>35 180.8</td><td>95.8</td><td>163.6</td></tr></tbody></table>

Note. $^a$This run does not include either star formation or stellar feedback.

is always resolved with $N_λ \approx 6$ ($N_λ \approx 3$) resolution elements max ($\epsilon_g$, $h$), assuming that the Jeans length is defined as $\lambda_J = c_s\sqrt{\pi/(G\rho)}$ ($\lambda_J = \sqrt{15k_B T/(4\pi G\rho \mu m_p)}$). The local sound speed of the gas is $c_s = \sqrt{\gamma k_B T/(\mu m_p)}$, where $\gamma = 5/3$ is the adiabatic index, $k_B$ is the Boltzmann constant, $\mu \approx 0.6$ is the mean molecular weight (assuming that most of the gas is ionized) and $m_p$ is the proton mass. The criterion on the Jeans length translates in the local Jeans mass resolved with $\approx(\epsilon_g/h)^3 N_λ^3/8 \approx 15.6 N_λ^3$ kernel masses, where we substitute $\epsilon_g/h = 5$ in the last equality after noting that the distribution of the ratio $h/\epsilon_g$ within the main galaxy is peaked at about 0.2. We note that we adopt a pressure floor to avoid numerical fragmentation, while we leave the temperature of the gas free to evolve according to the energy equation and the effect of radiative cooling and heating (e.g. Richings & Schaye 2016). This means that the equation of state of the gas does not effectively follow an ideal gas law anymore when relating the pressure to the density and the temperature. The usage of the pressure floor also imply that the gravitational collapse is at least partially suppressed for structures initially collapsing on scales below $N_λ \epsilon_g \sim 150-300$ pc (Bate & Burkert 1997). However, we have checked that the gas phase-space region where the chosen pressure floor kicks in overlaps with the conditions for forming stars only at high densities and low temperatures. This implies that the initial phases of gravitational collapses that lead to star formation are physical and well resolved, while the pressure floor ensures that spurious fragmentation on smaller scales is avoided. None the less, this may have a dynamical role in the evolution of the interstellar medium of the simulated galaxy; therefore, we also run an additional run from $z \approx 8$ (after that the disc has re-formed; see Section 3) to $z = 6.5$ adopting a lower pressure floor with $\alpha = 1$. This additional run is dubbed ‘PH_PF1’. We also restart an additional version of run PH from $z \approx 8$ to $z = 7.1$, including radiative cooling but without star formation and stellar feedback, in order to test the impact of feedback on the properties of the interstellar medium (see Section 3.4). We refer to this run as ‘PH_NF’.

In addition to run PH, we have also simulated a lower resolution version, i.e. ‘PH_LR’, adopting the same parameters. This run has a mass and force resolution eight times and two times coarser than the main run, respectively. The initial total number of particles is 17 952 072, while the main halo contains about 5974 677 particles within the virial radius at $z = 6.5$. We use this run for resolution tests that we show in the Appendix. We summarize the labelling and the main features of all the simulations in Table 1.

## 2.3 Halo detection

We identify dark matter haloes (and then the contained galaxies) using the AMIGA HALO FINDER (Gill, Knebe & Gibson 2004; Knollmann & Knebe 2009). Every halo is defined as a gravitationally bound group of at least 100 particles within a virial radius $R_{vir}$. Then, we construct the merger tree of the main halo matching the dark matter particle IDs among every two snapshots from $z = 6.5$ to $z \approx 30$ (e.g. Fiacconi et al. 2015). The main progenitor branch is determined as the progenitor halo that maximizes the number $f_{\text{shared}} = N_{\text{shared}}/\sqrt{N_h N_{\text{prog}}}$ through each snapshot, where $N_{\text{shared}}$ is number of particles shared among the halo and its progenitor, while $N_h$ and $N_{\text{prog}}$ are the particles of the halo and the progenitor, respectively (Fiacconi et al. 2015).

# 3 RESULTS

## 3.1 Evolution of the main galaxy

The global evolution of the main halo of the run PH is summarized in Fig. 1. Specifically, we show the merger tree of the main halo, highlighting the major mergers with mass ratio $q < 4$ (where we define $q \equiv M_1/M_2 > 1$), and the growth of the virial mass $M_{vir}$ with redshift. The early growth of the main halo is characterized by a sequence of several major mergers with mass ratios between $q = 1$ and $q = 3$ from $z \approx 18.5$ to $z \approx 11.6$. During this phase, $M_{vir}$ quickly grows from $<10^8 M_{\odot}$ to $\sim 4 \times 10^9 M_{\odot}$ through the rapid accretion of both dark matter and gas, also triggered by the major mergers. By $z \approx 11$, the main halo has formed, i.e. it has reached about 5 per cent of its final mass at $z = 6.5$, which is $M_{vir} \approx 1.2 \times 10^{11} M_{\odot}$. Slightly later, at $z \approx 10$, the virial radius of the main halo ($R_{vir} \approx 7.7$ kpc in physical units) begins to overlap with the virial volume of a nearby halo with mass 3.7 times smaller than the main one. This is the last major merger of the main halo during the time span that we have simulated. It completes by $z \sim 9$, though the two central galaxies merge slightly after at $z \approx 8$, as shown by Fig. 1. Such a time lag between the merger of the two haloes and their central galaxies, $\sim 100$ Myr from $z = 9$ to $z = 8$, is consistent with the dynamical friction time-scale for the central galaxies to sink at the centre of the remnant halo after a $q = 4$ merger, namely $\sim 0.2 t_H(z) \sim 170$ Myr, where $t_H(z) = H^{-1}(z)$ is the Hubble time at $z$ and the last equality is obtained assuming $z = 9$ (Krivitsky & Kontorovich 1997; Boylan-Kolchin, Ma & Quataert 2008; Jiang et al. 2008; Hopkins et al. 2010).

After the last major merger, the galaxy quickly builds an extended stellar and gaseous disc at $z \approx 7.5$. This is visible e.g. in the mock UVJ image at $z = 7.1$ in Fig. 1. Those images are obtained similarly to Fiacconi et al. (2015): we consider each stellar particle as a stellar population with a Kroupa (2001) initial mass function of age $\tau$ and metallicity $Z$ and we determine its total luminosity by interpolating a table based on the synthetic stellar population models of the Padova group (Marigo et al. 2008; Girardi et al. 2010; Bressan et al. 2012). These tables<fnref n="4" /> span the stellar age interval from 4 Myr to 12.6 Gyr and the metallicity interval from $5 \times 10^{-3}$ to $1.6 Z_{\odot}$. We neglect the effect of dust attenuation.

Then, the galaxy evolves nearly in isolation until $z \approx 6.5$, when it starts to interact with another galaxy in a mildly minor $q \approx 5$ merger. During this period, the total gas fraction $f_{\text{gas}} = M_{\text{gas}}/(M_* + M_{\text{gas}})$ within $0.1 R_{vir}$ oscillates around 50 per cent (only mildly increasing

<footnote n="4">⁴ Similar tables can be obtained at http://stev.oapd.inaf.it/cgi-bin/cmd.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Upper row, from left to right: sequence of images at $z = 6.5, 7.1, 8.1$ and $9.1$, respectively, of the stellar component of the main galaxy in the rest-frame U, V and J bands. Middle row: merger tree of the main halo as a function of redshift. The size of each circle is proportional to the logarithm of the virial mass. Green circles denote the main branch of the merger tree, while thick, red circles connected by thick, red lines mark the occurrence of a major merger with mass ratio $q < 4$. Lower panel: evolution of the virial mass of the main halo as a function of redshift. Red dots and vertical dashed lines show major mergers together with their mass ratio $q$. The main halo reaches about $10^{11} M_{\odot}$ by $z \gtrsim 6$.</figcaption></figure>

with time), while the fraction of cold gas ($T < 10^4$ K) is lower and about 37 per cent. The stellar mass, measured as the mass within $0.1R_{vir}$, is $M_* \approx 2.5 \times 10^9 M_{\odot}$. By changing the filtering scale (e.g. assuming a fixed volume of 3 physical kpc or 20 comoving kpc), the difference on the determination of $M_*$ is less than 20 per cent. During this phase, the global properties of the main galaxy, namely stellar mass, gas fraction and star formation rate (see below), are quantitatively similar<fnref n="5" /> to those observed in a few, strongly lensed galaxies at $z \sim 7$ by Bradley et al. (2012) and Watson et al. (2015). Moreover, the stellar mass is a factor 2 higher than the halo mass–stellar mass relations determined by Behroozi, Wechsler & Conroy (2013), but yet consistent within the $2\sigma$ uncertainty.

Fig. 2 shows the evolution of the star formation rate of the main galaxy as a function of time. This is continuously increasing with time, with isolated peaks corresponding to the major mergers highlighted in Fig. 1. Around $z \approx 6.5–7$, the main galaxy is forming ~20 $M_{\odot}$ yr$^{-1}$ in new stars. On the other hand, the specific star formation rate decreases from $\gtrsim 10$ Gyr$^{-1}$ to $\sim 5–6$ Gyr$^{-1}$ at $z = 6.5$. These numbers are consistent with recent determinations of the typical star formation rate of galaxies in samples at $z > 4$. Smit et al. (2012) have computed the star formation rate function (normalized as the number of galaxies per unit volume and unit star formation

rate) starting from dust-corrected UV luminosity functions from Bouwens et al. (2007, 2011). They have found that between redshift 6 and 8 the star formation rate corresponding to the knee $L_*$ of a Schechter luminosity function is about $10 M_{\odot}$ yr$^{-1}$, suggesting that our galaxy is probably slightly larger than $L_*$. By correcting the value of $L_*$ at $z = 6.8$ determined by Bouwens et al. (2011) as described by Smit et al. (2012), we get a corresponding absolute magnitude after dust correction $M_{UV}^* = -20.6$, which is indeed slightly larger (i.e. less bright) than the rest-frame magnitude $M_U = -20.8$ of the stellar component within $0.1R_{vir}$ measured from run PH. We caution however that the observational measurements have been obtained at 1600 Å, while our estimate from the simulation comes from U rest-frame. We also compare the specific star formation rate as a function of mass with a population of slightly lower redshift galaxies ($4.5 < z < 5.5$, with spectroscopic redshifts) observed in the VIMOS Ultra-Deep Survey (VUDS) by Tasca et al. (2015). At the mass of a few $10^9 M_{\odot}$, our galaxy sits slightly above the bulk of the data at $z \sim 5$ (thought there is large scattering), which in turn are distributed around the relation inferred by Schreiber et al. (2015) from Herschel data. This is also consistent with the slowly growing trend of the average specific star formation rate with increasing redshift derived by Tasca et al. (2015), $\propto (1+z)^{1.2}$, which would predict $\sim 6$ Gyr$^{-1}$ at $z = 6.5$, although most of the galaxies used to get this determination have $M_* > 10^{10} M_{\odot}$. Overall, the global properties of the main galaxy in run PH seems to be consistent with typical star-forming galaxies recently observed at similar redshifts (Iye et al. 2006; Bradley et al. 2012; Watson et al. 2015).

<footnote n="5">⁵ Watson et al. (2015) have observed a galaxy at $z = 7.5 \pm 0.2$ with stellar mass $1.7_{-0.5}^{+0.7} \times 10^9 M_{\odot}$, star formation rate $9_{-2}^{+4} M_{\odot}$ yr$^{-1}$ (from infrared light) and gas fraction $55 \pm 25$ per cent.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Upper panel: time evolution of the star formation rate ($\dot{M}_{\star}$, blue continuous line, left y-axes) and of the specific star formation rate ($\dot{M}_{\star}/M_{\star}$, red dotted line, right y-axes). The latest two major mergers of the main halo are highlighted with their mass ratio $q$. Lower panel: specific star formation rate as a function of the stellar mass for our simulation at $z = 6.5$ (PH, blue square) and for the data from the VUDS survey in the redshift range $4.5 < z < 5.5$ (Tasca et al. 2015). The VUDS data are represented as red filled and empty circles when the redshift determination is ~100 per cent and ~70–75 per cent reliable, respectively. The thin dotted line shows the VUDS completeness limit in the determination of the star formation rates. The dashed and dot-dashed lines show the determination of the main sequence of star-forming galaxies at $z \approx 2$ from Daddi et al. (2007) and the analytical fit of the main sequence by Schreiber et al. (2015) extrapolated to $z = 5$, respectively.</figcaption></figure>

## 3.2 Structure of the galactic disc

Fig. 3 shows the evolution of the galactic disc in run PH between $6.5 < z < 7.6$, after the last major merger is completed. The galaxy is oriented face-on by determining the specific angular momentum of the gas within a sphere of radius $0.1R_{vir}$, after having centred the galaxy on the minimum of the potential and having removed the systemic velocity evaluated as the mass-weighted velocity of the particles within 500 physical pc from the centre. During this interval of time, the galaxy appears as a disc with tenuous nearly

axisymmetric structures in the stellar component. At $z = 7.6$, the two cores originally at the centre of the merging galaxies are finally coalescing at the centre of the remnant, forming a tiny bulge (also visible as a redder component in Fig. 1).

The surface density map of the gas reveals that the gaseous disc is highly inhomogeneous, with many overdense regions. This is particularly evident at $z = 6.5$, when the disc is less axisymmetric and more perturbed than at slightly larger redshifts, with more extended regions of dense gas at and above ~$10^8 M_{\odot}$ kpc$^{-2}$, possibly because of the interaction with the companion galaxy visible in Fig. 1. Indeed, the gas is highly multiphase, as also confirmed by the temperature maps. Local fluctuations can vary from temperatures <1000 K to ~$10^6$ K on scales of $\gtrsim 100$ pc. Those high temperatures are likely triggered by local injection of energy from stellar feedback.

Such a structure of the interstellar medium is potentially a signature of turbulence, as expected from previous work on high-z galaxies (e.g. Green et al. 2010; Bournaud et al. 2011; Hopkins et al. 2012b), though the main driver of turbulence is still debated. Indeed, we observe large fluctuations in the local velocity dispersion of the gas, calculated as SPH kernel-weighted average of the local standard deviation of the velocity, ranging from a few to almost 100 km s$^{-1}$, with a typical, average value $\sigma_g \sim 40$ km s$^{-1}$ across the disc at different redshifts. We discuss the turbulence within the disc in more detail in Section 3.3.

High-density regions are naturally associated with larger star formation rates, following a large-scale Kennicutt–Schmidt relation (Kennicutt 1998; Bigiel et al. 2008) with slope $n \approx 1–1.5$ at (total) surface densities $> 1 M_{\odot}$ pc$^{-2}$, consistent with the local volumetric Schmidt law (Gnedin, Tasker & Fujimoto 2014). The star formation within the disc is mostly distributed along the major features (i.e. spiral arms) in the stellar disc, with typical star formation rate surface densities $\gtrsim 1 M_{\odot}$ yr$^{-1}$ kpc$^{-2}$, while overdensities are associated with isolated regions of intense star formation as high as $\gtrsim 25 M_{\odot}$ yr$^{-1}$ kpc$^{-2}$.

Fig. 4 shows the rotation curve $V_{rot}$ and the surface brightness profiles of the main galaxy in run PH at $z = 7.1$. The rotation curve is obtained as the mean over 50 random line of sights through the disc mid-plane. We take into account only cold/warm gas with temperature <50 000 K to mimic the rotation curve that could be obtained through $H\alpha/H\beta$ observations. The galaxy has a flat rotation curve up to a few kpc from the centre with a nearly constant value about 100 km s$^{-1}$. The circular velocity $V_{circ} = \sqrt{GM(<r)/r}$ is also flat, but with an asymptotic large velocity about 140 km s$^{-1}$. This mismatch witnesses the degree of non-radial motion in the gas, also visible in the fluctuation of $V_{rot}$ at $R \ge 2$ kpc. The surface brightness profiles of the stellar component in the rest-frame U,V and K band show that the stellar disc is nearly exponential, with a disc scale radius (half-light radius) about 460 physical pc (760 physical pc) in K band, which is grossly consistent with the size of typical galaxies at $z \sim 7$ (Oesch et al. 2010b). The stellar disc drops slightly outside 2 kpc and it shows a tiny steepening for radii below ~200–300 pc owing to the presence of a tiny bulge. We perform a profile decomposition of the surface brightness in two exponential profiles, one for the disc and on for the bulge, finding a similar value of $B/T \approx 0.03$ per cent (i.e. fraction of the total light in the bulge component) across the three photometric bands.

Fig. 5 shows the edge-on evolution of the disc that forms in run PH. The disc mid-plane corresponds to $y = 0$, though it might look difficult to identify a smooth structure because the disc is indeed highly perturbed in the vertical direction, showing a clumpy and discontinuous structure. We calculate the vertical density profile of
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Evolution of the galactic disc in run PH between $z = 7.6$ and $z = 6.5$ (from bottom to top). From left to right: gas surface density, local gas temperature in the disc mid-plane, local gas velocity dispersion (computed locally as an SPH weighted average) in the disc mid-plane, stellar surface density, and star formation rate surface density. All quantities are measured within a 200 pc thick slice centred on the disc plane. The star formation rate surface density is determined from the surface density of stellar particles younger than 10 Myr. The circles in the lower-right corners of the first column have radii equal to the gravitational softening of the gas, 47.9 physical pc.</figcaption></figure>

the gas, azimuthally averaging the gas density in five linearly spaced bins between 20 pc and 2 kpc from the centre of the disc. We find that the gas density drops by almost an order of magnitude after a typical distance from the disc mid-plane between 100 and 300 pc, with the slight tendency to increase at larger cylindrical radii (i.e. the disc flares at large radii), as it is mildly hinted by Fig. 5. Indeed, we fitted the vertical profile of the gas density with an exponential law, finding typical vertical scales ≈300 pc. A similar result (even a slightly ticker disc, probably because of the different coupling between feedback energy and dense gas) is obtained analysing the PH_PF1 run with the reduced pressure floor, suggesting that the latter does not play a major role in setting the vertical structure of the disc. We repeat the same analysis for the stellar disc and we find that the latter is somewhat thinner, with exponential vertical scales typically ≈170–200 pc. This is consistent with stars forming preferentially in the denser and thinner cold gas disc.

The vertical temperature structure shows interesting features. The dense, clumpy gas that can be identified in the central disc by comparing with the density map is typically cold, with temperatures between ≳100 and ~8000 K owing to metal cooling, and it is embedded in an atmosphere of hot gas ≳10⁵ K. The hot gas is typically outflowing vertically from the central disc, funnelled through clumps of colder gas, likely suggesting that it originates from episodes of supernova feedback. This gas is indeed polluted to metallicity ≲Z⊙ and is injected into the circumgalactic medium.

Fig. 6 shows the probability distribution of the outflow velocities $V_{\text{out}}$ (i.e. the radial velocity $v_r > 0$ of outflowing gas) of the gas

particles within a spherical shell between 0.2 and $0.3R_{\text{vir}}$ from the central galaxy (e.g. Muratov et al. 2015). The particles are selected to have $v \cdot \hat{r} > 0$, where $v$ is the particle velocity (after removing the systemic velocity of the halo) and $\hat{r}$ is the versor along the direction from the centre of the halo to the particle position. The distribution is roughly exponential, i.e. $\propto \exp(-V_{\text{out}}/\tilde{V})$, where $\tilde{V} \approx 80 \text{ km s}^{-1}$ at $z = 7.1-6.5$ ($\approx 60 \text{ km s}^{-1}$ at $z = 7.6$) is a typical scale value for the positive radial velocity. In particular, $\tilde{V} \approx \sigma_r \equiv \sqrt{\langle V_{\text{out}}^2 \rangle - \langle V_{\text{out}} \rangle^2}$ (where $\langle \cdot \rangle$ is intended as mass-weighted average), suggesting that the genuinely outflowing material populates the extended tails up to velocities as high as 200–300 km s$^{-1}$. However, the fraction of the outflowing mass at high velocity, more specifically with $V_{\text{out}}$ larger than the escape velocity $V_{\text{esc}}$, is typically low, oscillating between 1 and 5 per cent, i.e. most of the gas expelled from the disc will be recycled through the halo and will eventually join the central galaxy. Here, we define the escape velocity $V_{\text{esc}}$ at the position $\mathbf{x}$ of a particle as $V_{\text{esc}} = \sqrt{2|\phi(\mathbf{x}) - \phi_0|}$, where $\phi$ is the local gravitational potential and $\phi_0$ is the reference potential at $R_{\text{vir}}$ calculated as the average gravitational potential of all particles in a thin spherical shell between $0.95R_{\text{vir}}$ and $R_{\text{vir}}$. We have repeated this analysis within two additional spherical shells, (i) with an outer radius $0.2R_{\text{vir}}$ and a thickness of 1 physical kpc, and (ii) between 2 and 2.5 physical kpc, and we have found very similar results, within a factor 2.

We have also computed the mass loading factor $\eta \equiv \dot{M}_{\text{out}}/\dot{M}_{\star}$, where $\dot{M}_{\text{out}} = \Delta r^{-1} \sum_i m_i v_i \cdot \hat{r}_i$ is the mass outflow rate, the sum involves only gas particles $m_i$ within the considered shell with
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Upper panel: rotation curve $V_{\text{rot}}$ at $z = 7.1$ from the cold gas ($T < 50\,000\,\text{K}$). The blue solid line is the average over 50 random line of sights through the galaxy mid-plane, while the shaded region marks the 68 per cent of the measurements. The red dashed line is the circular velocity $V_{\text{circ}} = \sqrt{GM(<r)/r}$ shown for comparison. Lower panel: rest frame surface brightness profiles of the stellar component in $U$ (blue solid), $V$ (green dashed) and $K$ (red dotted) bands at $z = 7.1$. The grey shaded region at the centre marks two gravitational softenings. The surface brightness profiles include cosmological dimming. The grey solid line shows the best fit of the profile decomposition of $\mu_*$ into two exponential profiles for the $K$ band as an example.</figcaption></figure>

$v_i \cdot \hat{r}_i > 0$, and $\Delta r$ is the radial thickness of the shell. The results between $z = 7.6$ and $z = 6.5$ are shown in the lower panel of Fig. 6, where we compare with the observational determinations of $\eta$ by Schroetter et al. (2015) and references therein. They are plotted against $V_{\text{max}}$, i.e. the asymptotic maximum of the rotation curve (Schroetter et al. 2015). We repeat the same analysis as for Fig. 4 to calculate $V_{\text{max}}$ from the simulation data, whose ‘uncertainty band’ in Fig. 6 is associated with the measurements within different shells. We also show the recent determination of $\eta \approx 2$ (where we add a generous uncertainty of 50 per cent) by Gallerani et al. (2016), which represents the average $\eta$ estimated for the sample of $z \sim 5$ and $M_* \sim 10^{10} \, \text{M}_\odot$ galaxies presented in Capak et al. (2015). In this case, we estimate $V_{\text{max}}$ by using the stellar mass–halo mass relation by Behroozi et al. (2013) to determine the virial velocity of the typical halo where such galaxies are supposed to live in at $z = 5$, with an uncertainty of 50 km s$^{-1}$. We find $\eta \approx 0.5–1$ for PH,

in fair agreement with the observational estimates. We have also repeated the analysis selecting particles with $v_i \cdot \hat{r}_i > \sigma_r$ and finding slightly lower values for $\eta$ by at most 40 per cent. We caution though that the observational data from Schroetter et al. (2015) are a collection of measurements at much lower redshift than our simulation, namely between 0.1 and 0.8. However, it is noteworthy that the observational estimate at $z \sim 5$ by Gallerani et al. (2016) is yet consistent with the lower redshift data and in particular with local starbursts (Heckman et al. 2015). Moreover, the agreement at face values is promising, given the scattering of order $> 10$ in $\eta$ among different successful theoretical/empirical models (see e.g. fig. 10 of Schroetter et al. 2015).

On the other hand, clumps and streams of cold gas are typically raining down radially on the galactic disc. Some of that gas is mildly polluted $Z \sim 0.2 \, \text{Z}_\odot$ and typically lives in overdense regions, compressed by the surrounding hot outflows, where the cooling time is short, likely because of the metals already seeded by previous stellar outflows (Costa, Sijacki & Haehnelt 2015). However, some of the inflowing material has low metallicity, possibly coming from direct feeding through cold flows. We select gas particles at $z = 6.5$ within a cubic box 6 physical kpc in size centred on the main galaxy (removing a slab 400 pc thick centred on the disc mid-plane) with temperature $T < 5 \times 10^4 \, \text{K}$, i.e. the gas that is (or will be soon) able to form stars. We trace them back among snapshots to $z = 8$, storing the time-evolution of the density, temperature and metallicity. We select the gas that is directly inflowing from cold flows as the gas whose temperature never exceeds $10^{5.5} \, \text{K}$, finding that it represents $\sim 50–60$ per cent of the mass of cold gas originally selected in the trial volume. A fraction $\sim 5–10$ per cent of that gas is also accreted at nearly primordial composition as its metallicity always remains below $0.1 \, \text{Z}_\odot$. We plan to devote a more detailed analysis to the gas inflow form larger scales as well as to the recycling of the gas within the halo in a forthcoming publication.

## 3.3 The properties of the disc turbulence

The gaseous disc of the main galaxy is turbulent and multiphase. We show the mass-weighted probability density function of the gas within the disc at the representative redshift $z = 7.1$ in Fig. 7. We select the gas within a cylinder with a radius of 2 physical kpc and 500 pc thick, centred on the disc mid-plane. The density distribution is not approximately a lognormal distribution, as generally found in simulations of isothermal turbulence (e.g. Padoan, Nordlund & Jones 1997; Price & Federrath 2010). Instead, it peaks at a few $M_\odot \text{ pc}^{-3}$ ($\sim 100 \, \text{H cm}^{-3}$), with an extended tail at low densities well below $10^{-2} M_\odot \text{ pc}^{-3}$ ($\sim 0.2 \, \text{H cm}^{-3}$) and a fast decline around $\sim 10 M_\odot \text{ pc}^{-3}$ ($\sim 250 \, \text{H cm}^{-3}$). This more complicated shape is a natural consequence of the mixing different phases of gas temperatures (Robertson & Kravtsov 2008). We can identify three main, nearly isothermal components<fnref n="6" /> characterized by an approximately lognormal distribution: a ‘cold’ component at $T < 10^3 \, \text{K}$, a ‘warm’ component with $5 \times 10^4 < T/K < 5 \times 10^5$ and a ‘hot’ component at temperatures larger than $5 \times 10^5 \, \text{K}$. They are associated with gas at average densities of 100, 10, and $1 \, \text{H cm}^{-3}$, respectively. As a whole, they represent 65 per cent of the total gas in the

<footnote n="6">⁶ In the following, we dub the different components within quotes in order to avoid confusion with the common phases of the interstellar medium (Ferrière 2001). In fact, their names are just labels for the temperature ranges used for the selection and, though similar to some of the phases of the interstellar medium, they do not refer directly to them.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** From left to right: edge-on view of the main galactic disc in run PH between $z = 7.6$ and $z = 6.5$. First row: gas temperature map averaged within a 200 pc thick slice perpendicular to the disc plane and centred on the disc centre. The arrows show the velocity field of the gas. Second row: the same of the top row, but for the gas metallicity. Third row: gas density map averaged within a 200 pc thick slice as for the temperature and the metallicity. Fourth row: stellar surface density projected within 4 kpc centred on the disc centre. The green lines mark equal surface density contours corresponding to $10^7$, $10^8$ and $10^9$ M$_\odot$ kpc$^{-2}$, from outside inward.</figcaption></figure>

disc, sub-divided between ‘cold’, ‘warm’ and ‘hot’ as ≈60 per cent, ≈29 per cent and ≈11 per cent, respectively. The lognormal fits of the ‘warm’ and ‘cold’ components do not agree very well with the data at large densities. This is possibly due to the effect of the pressure floor that limits the development of very high overdensities. Indeed, repeating the same analysis on the run PH_PF1 we find a lower peak of the ‘cold’ component, which instead extends up to densities a factor ~2–3 larger, up to about 20 M$_\odot$ pc$^{-3}$. In this case, the high-density tail is nicely described by the lognormal distribution and it does not show a power-law behaviour expected for highly self-gravitating turbulent flows (e.g. Scalo et al. 1998; Federrath et al. 2008; Kritsuk, Norman & Wagner 2011), at least at the densities of a few hundreds H cm$^{-3}$ that we are able to probe at the resolution of our simulation.

We analyse the properties of the turbulence by computing the power spectrum of velocity and density fluctuations at different scales. Given any quantity $w$, its two-points correlation function within a volume $V$ is defined as

$$ \xi_w(l) = \frac{1}{V} \int_V w(\mathbf{x} + l)w(\mathbf{x})d^3\mathbf{x}, \qquad (2) $$

which can be generalized for a vectorial quantity by means of the inner product $w(\mathbf{x}) \cdot w(\mathbf{x} + l)$. The power spectrum of $w$ is the Fourier transform of $\xi_w$, which can be rearranged as

$$ p_w(\mathbf{k}) = \frac{1}{(2\pi)^{3/2}V} \left| \int_V w(\mathbf{x})e^{-i\mathbf{k} \cdot \mathbf{x}} d^3\mathbf{x} \right|^2. \qquad (3) $$

We further assume that each mode is isotropic and it depends only on the module of the wavenumber $k = |\mathbf{k}|$. Then, the power spectrum becomes

$$ P_w(k) = k^2 \int_{4\pi} p_w(k, \theta, \phi) \sin\theta d\theta d\phi, \qquad (4) $$

where $(\theta, \phi)$ are the spherical coordinates in $\mathbf{k}$-space.

We perform this calculation for the gas surface density and for the velocity. In the first case, we build a two-dimensional surface density map $\Sigma_{i,j}$ by projecting the SPH gas density field within a cube of $L_{\text{box}} = 4$ physical kpc per side centred on the main galaxy (oriented face-on as described in Section 3.2) on a grid of 512 × 512 square grids $\Delta x \approx 7.8$ pc per side, similar to the typical smoothing length within the disc, i.e. ~10 pc. Then, we use the fast Fourier transform to numerically calculate the two-dimensional power spectrum $P_\Sigma$, after having padded with zeros the boundaries of the grid in order to
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Upper panel: mass-weighted probability distribution of the outflow velocity $V_{\text{out}}$ at $z = 6.5$ (blue solid line), $z = 7.1$ (red dashed line) and $z = 7.6$ (green dotted line). The arrows show the escape velocities at each $z$. The black dotted lines shows $dP/dV_{\text{out}} \propto \exp(-V_{\text{out}}/\bar{V})$, with $\bar{V} = 80 \text{ km s}^{-1}$. Lower panel: mass loading factor $\eta$ as a function of the asymptotic maximum rotation velocity. The red circles are the observational data from Schroetter et al. (2015) and references therein, the green square is the average determination of Gallerani et al. (2016) from the catalogue of $z \sim 5$ galaxies by Capak et al. (2015), and the blue line shows the results for run PH between $z = 7.6$ and $z = 6.5$. The blue shaded area spans the results associated with different shells used to measure $\dot{M}_{\text{out}}$ (see the text for details).</figcaption></figure>

reduce the aliasing in Fourier space produced by the non-periodic and finitely sampled content of the grid.

In the second case, we analyse the velocity fluctuations of the interstellar medium in the disc. First, we selected three boxes within the disc with side $L_{\text{box}} = 512$ physical pc at 1.5 kpc from the centre and at 120° angular separation from each other, in order to avoid any correlation among them and to minimize the effect of differential rotation since the rotation curve is roughly flat outside the central kpc (see Fig. 4). Each box typically contains from at least a few tens of thousands to a few hundreds of thousands gas particles. We subtract the systemic velocity of each cube as the mass-weighted gas velocity and we then interpolate the SPH-averaged gas velocity on a cubic grid of $128 \times 128 \times 128$ using the TIPGRID code.<fnref n="7" /> Then, we finally proceed as described above to calculate the velocity power spectrum $P_v$. The velocity power spectrum $P_v$ is normalized such

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Mass-weighted probability density function of the gas density within the galactic disc at $z = 7.1$. The black solid line shows the total distribution, while the blue, green and red solid lines show the distribution of the 'cold' ($T < 10^3$ K), 'warm' ($5 \times 10^4$ K $< T < 5 \times 10^5$ K) and 'hot' ($T > 10^5$ K) gas, respectively. The dashed lines of the same colour show the lognormal best fit of each component. The thin grey line shows the compositions of the lognormal best fits, renormalized by the fractional mass occupied by the three selected components (Robertson & Kravtsov 2008). The fits describe fairly accurately the total density distribution.</figcaption></figure>

as its integral over $k$ gives the total velocity dispersion $\sigma_g^2$ of the gas in each box, while $P_\Sigma$ is normalized such as its integral over $k$ gives 1. In both cases, we compute the power spectrum between $k_{\text{box}} = 2\pi/L_{\text{box}}$ and the Nyquist wavenumber $k_{\text{Nyq}} = \pi/\Delta x$.

The results of these calculations are shown in Fig. 8. The power spectrum of the surface density grossly shows a two power laws behaviour, with a break at $k \sim 30 \text{ kpc}^{-1}$, which corresponds to a scalelength $\sim 200$ pc. This is roughly consistent with our previous estimate of the disc scaleheight, as generally found in simulations (Bournaud et al. 2010) as well as observations (Elmegreen, Kim & Staveley-Smith 2001). This suggests a transition between more two- and three-dimensional like turbulence at low and high $k$ (or large and small length-scales), respectively. The power-law exponents oscillate significantly with redshift, ranging roughly between $-1.0$ and $-1.6$ for low $k$ and between $-2.3$ and $-2.7$ for high $k$, which in both cases is slightly less than what Bournaud et al. (2010) found. Comparing the results with the run PH_PF1, we find similar scatterings, but marginally different exponents (note however that the power spectra at $z = 7.1$ are fairly similar among the two runs). Typically, the large-scale modes have slightly shallower slopes between $-0.6$ and $-1.2$, while the power at small scales decays faster with $k$, typically having steeper slopes around $-3$. Moreover, the transition between the low- and high-$k$ branch and the latter as well are at somewhat higher normalization compared with the reference run. This different behaviours imply that there is slightly more power in large-scale modes and at the transition between two- and three-dimensional modes at smaller scales than in the reference simulation, likely because the lower pressure support promotes slightly larger overdensities to develop under the

<footnote n="7">⁷ TIPGRID was written by Joachim Stadel and it is available at http://astrosim.net/code/doku.php?id=home:code:analysistools:misctools.</footnote>
<figure><img src="image_10.png" /><figcaption>Figure 8. Upper panel: two-dimensional surface density power spectra of the gas in the galactic disc in a box $L_{\text{box}} = 4$ kpc. Blue solid, red dashed and green dotted lines refer to $z = 6.5, 7.1$ and $7.6$, respectively. Thick and thin lines are associated with run PH and PH\_PF1, respectively. Lower panel: three-dimensional velocity power spectra, obtained as average between three boxes $L_{\text{box}} = 512$ pc (see the text for details). The grey bands show the minimum–maximum range between the boxes at each redshift for run PH. The line style is the same as above. The dashed lines show $P_{\nu} \propto k^{-5/3}$ and $\propto k^{-2}$ for visual guidance. In both panels, the vertical dotted lines mark the $k$ associated with the gravitational softening, while the grey shaded regions refer to $k < k_{\text{box}}$ and $k > k_{\text{Nyq}}$.</figcaption></figure>

influence of self-gravity, though the gross structure of the disc (i.e. the disc thickness) roughly remains the same.

At globally smaller scales, the velocity power spectra (lower panel of Fig. 8) show much less differences in shape among the runs with different pressure floors, i.e. runs PH and PH\_PF1, though there are fluctuations of factors ~2–3 in the normalizations at different redshifts. This suggests that the development of turbulence is only mildly influenced by the presence of the pressure floor at any scale larger than the gravitational softening. The typical power-law slope measured between $k_{\text{box}}$ and $k = 100$ kpc$^{-1}$ (i.e. at scales larger than the gravitational softening) is about 1.6 with little variations among runs and redshifts. This looks consistent with a Kolmogorov-like

turbulence spectrum, i.e. for incompressible, subsonic turbulence; however, when we remove from the fits the first two modes, which are poorly sampled due to the box size, we obtain a slope close to −2, again with little scattering among runs and redshifts. Such a slope matches the Burgers law expected for compressible supersonic turbulence (e.g. Federrath 2013). However, it is difficult to unequivocally discriminate between the two descriptions; in fact, the gas is a mixture of phases characterized by different temperatures and Mach numbers since the density distribution of each box is similar to what shown in Fig. 7. We typically find a transonic or mildly supersonic turbulence with Mach number $\sigma_g/\langle c_s \rangle \gtrsim 1$, where $\langle c_s \rangle$ is the mean sound speed determined from the mass-weighted average temperature of the gas (i.e. averaged over the different phases) about $10^5$ K. None the less, we argue that the gas within the disc might be slightly better described by a Burgers-like law because of its compressible nature and the larger mass fraction in cold and supersonic gas.

### 3.4 What shapes the interstellar medium at high redshift?

Several physical processes may contribute to fostering turbulence in the interstellar medium. Among them, gravity and stellar feedback are often considered to be the dominant ones (e.g. Gómez & Cox 2002; Joung & Mac Low 2006; Brunt, Heyer & Mac Low 2009; Pan et al. 2016). Gravity can trigger turbulent motions in a disc – gravito-turbulence – when the local cooling time is factor of few to several the local orbital time. Despite that the exact transition is still matter of debate (Gammie 2001; Lodato & Clarke 2011; Meru & Bate 2011a,b; Paardekooper, Baruteau & Meru 2011), a gravito-turbulent state is typically characterized by a Toomre parameter $Q \approx 1.5–2$. Fig. 9 shows the Toomre parameter of the disc in run PH at different redshifts. We measure the Toomre parameter within a disc 2.5 physical kpc in radius and 500 physical pc thick. The general definition of the Toomre parameter is

$$Q = \frac{\kappa V}{AG\Sigma}, \qquad (5)$$

where $\kappa = \sqrt{2(V_{\phi}/R)^2(1 + d \log V_{\phi}/d \log R)}$ is the epicyclic frequency defined through the azimuthal velocity $V_{\phi}$ and the polar radius $R$, and $\Sigma$ is the surface density (either of the gas or the stars). The factor $A$ is $A_g = \pi$ and $A_\star = 3.36$ for gas and stars, respectively. The velocity $V$ is the radial velocity dispersion $V = \sigma_R$ for stars, and the sound speed $V = c_s$ for gas. If the gas is turbulent with a radial velocity dispersion $\sigma_{g,R}$, the ‘turbulent’ gas Toomre parameter $Q_{g,turb}$ adopts the corrected velocity $V = \sqrt{c_s^2 + \sigma_{g,R}^2}$.

We correct the Toomre parameter given by equation (5) for finite disc thickness effects following Romeo & Wiegert (2011, see also Romeo 1994; Romeo & Falstad 2013; Inoue et al. 2016); we multiply $Q$ by

$$T = \begin{cases} 1 + 0.6(\sigma_z/\sigma_R)^2 & \sigma_z/\sigma_R < 1/2 \\ 0.8 + 0.7(\sigma_z/\sigma_R) & \sigma_z/\sigma_R \ge 1/2 \end{cases}, \qquad (6)$$

where $\sigma_z$ and $\sigma_R$ are, respectively, the vertical and (polar) radial velocity dispersion, either of gas or stars. We generally measure the velocity dispersions as $\sigma = \sqrt{\langle v^2 \rangle - \langle v \rangle^2}$, where $\langle \cdot \rangle$ is the SPH average on the smoothing kernel of the appropriate velocity $v$.

The Toomre parameter describes the stability of rotating discs. However, the stellar component has a tiny bulge at the centre, which is dispersion dominated. Therefore, we calculate the ratio $\epsilon = j_z/j_{\text{circ}}$, where $j_z$ is the z component of the specific angular momentum of a particle computed in a reference frame centred
<figure><img src="image_11.png" /><figcaption>Figure 9. Radial profiles of Toomre parameter Q at z = 6.5, 7.1 and 7.6. Solid, dashed, dot-dashed and dotted lines refer to the total Q<sub>tot</sub> (according to Romeo & Wiegert 2011), the gas Q<sub>g</sub> without the contribution of the gas velocity dispersion, the gas Q<sub>g,turb</sub> including the gas velocity dispersion, and the stellar Q<sub>*</sub> for stars in the disc with j<sub>z</sub>/j<sub>circ</sub> > 0.5 (see the text for additional details). The vertical grey region marks the gravitational softening, while the horizontal ones show the regions of instability and marginal stability.</figcaption></figure>

on the galaxy, after having aligned the xy plane with the galac-
tic mid-plane. The maximum angular momentum of a particle at
distance r from the centre is the angular momentum on a circular
orbit, namely j<sub>circ</sub> = rV<sub>circ</sub>(r), where we estimate the circular veloc-
ity V<sub>circ</sub> &asymp; &radic;<span style="text-decoration: overline;">GM(<i>r</i>)/r</span>. In order to identify the stars that belong
to the rotating disc, we select those with &epsilon; > 0.5, i.e. the stars
whose angular momentum perpendicular to the disc plane is close
to maximal. Finally, we calculate the total Toomre parameter as

$$
Q_{\text{tot}}^{-1} = \begin{cases} W Q_{*}^{-1} + Q_{g}^{-1} & Q_{*} \ge Q_{g} \\ Q_{*}^{-1} + W Q_{g}^{-1} & Q_{*} < Q_{g} \end{cases} \tag{7}
$$

where $W = 2V_*V_g/(V_*^2 + V_g^2)$, where $V_g$ and $V_*$ stand for the ve-
locity adopted in equation (5) for the gas and the stars, respectively
(Romeo & Wiegert 2011).

Fig. 9 shows that the galactic disc is overall stable to gravitational perturbations, with the Toomre parameter typically ranging from ~4 to ~10 across the redshift interval z = 6.5–7.6, i.e. after the disc rebuilt. This is true for both the stellar and the gas components alone, and also when they are combined in Q<sub>tot</sub>, though the latter is slightly lower than the singular components. Such a high value for Q<sub>g,turb</sub> likely indicates that the gaseous disc is not in a gravito-turbulent state. The only exception is at z = 6.5, when the Toomre parameter of the gas component approaches 1.5–2 close to the centre, which in fact corresponds with a clumpy region in the gas, where overdensities might have been tidally enhanced by the close passage of the satellite galaxy that will eventually merge with the main one. We have checked that Q<sub>*</sub> goes to &le;1 within the inner 100–200 physical pc if we do not select particles with &epsilon; > 0.5, that is kinematically consistent with the presence of a central tiny bulge.

Fig. 9 also compares the Toomre parameter of the gas with and
without the contribution of turbulence (i.e. Q<sub>g,turb</sub> and Q<sub>g</sub>, respec-
tively). When the effect of turbulent motions is included, the Toomre
parameter is naturally larger. None the less, such effect accounts for
an increase in Q up to a factor ~2. This implies that the pres-
sure support alone would be enough to prevent widespread frag-
mentation in the disc, because the average temperature ~10<sup>5</sup> K
within the disc is such that the average sound speed of the gas
is &lang;c<sub>s</sub>&rang; ~ &sigma;<sub>g</sub> ~ 50 km s<sup>-1</sup>, i.e. pressure support and non-thermal
turbulence contributes similarly to the stability of the gaseous disc.

The high value of the gas Toomre parameter and the mild sensi-
tivity of the latter to turbulence suggest that (i) feedback may have a
major role in regulating the stability of the disc, and (ii) disc might
not be in a gravito-turbulent state, i.e. feedback is also the main
energy source for the gas turbulence. However, the latter point is
hard to demonstrate because gravity can still contribute by acceler-
ating free falling cold clouds from the surroundings that eventually
dissipate their kinetic energy by mixing with the denser gas of the
disc. In order to test the impact of feedback on the structure of
the interstellar medium, we have restarted the run PH without star
formation and feedback from z &asymp; 8, i.e. run PH_NF. Fig. 10 shows
the comparison between the gaseous disc at z = 7.1 in runs PH
and PH_NF. The two discs are visibly different. Under the effect of
feedback, the disc in run PH is overall slightly more extended in ra-
dius and it has a lower surface density on average. The gas in run PH
has larger density contrasts on small scales and a more flocculent
structure than in run PH_NF, which on the other hand shows a clear
two-armed spiral, initially triggered by the quasi-resonant tidal in-
teraction with a flying-by satellite (e.g. D’Onghia et al. 2010). The
disc in run PH_NF is also thinner than in run PH and it presents a
central bulge-like structure that is not dispersed by star formation
and feedback and that is also mildly elongated in a bar-like fashion.

Fig. 11 tries to quantify the morphological differences between the results of runs PH and PH_NF (comparing also run PH PF1) at z = 7.1. We compare the radial velocity dispersion of the gas σ<sub>R</sub>, the ratio V<sub>φ</sub>/σ<sub>g</sub> between the gas azimuthal velocity V<sub>φ</sub> and the three-dimensional velocity dispersion σ<sub>g</sub> (V<sub>φ</sub>/σ<sub>g</sub> quantifies the rotational support of the gaseous disc), and Q<sub>g,turb</sub>. We note that the gas in run PH_NF tends to have a slightly lower σ<sub>R</sub> ~ 25 km s<sup>-1</sup> than in both runs PH and PH PF1, except for a ‘bump’ at ~300–500 pc where the rotating disc joins the central spheroid. Below this radius, σ<sub>R</sub> drops to a few km s<sup>-1</sup> in run PH_NF, while it remains almost constant in the other cases, likely because the central gaseous bulge in run PH_NF is mostly pressure supported. This suggests an overall similar amount of turbulence with and without stellar feedback;<fnref n="8" />

<footnote n="8">⁸ The similarity of gas velocity dispersion achieved with and without feedback in presence of gravito-turbulence was pointed out in Agertz et al. (2009).</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Gas surface density maps of the gaseous disc of the main galaxy at $z = 7.1$ in runs PH (left column) and PH_NF (right column). The projections are face-on (upper row) and edge-on (lower row); sizes are in physical units.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Comparison between the disc properties in runs PH (blue solid line), PH PF1 (red dashed line) and PH_NF (green dotted line). From top to bottom: radial velocity dispersion of the gas, $V_\phi/\sigma_g$ ratio of the gas, and turbulent Toomre parameter $Q_{g,turb}$ of the gas. The vertical grey region common to all panels marks the gravitational softening of the gas.</figcaption></figure>

however, other dynamical properties of the disc are very different among the runs. Without stellar feedback, the disc is thinner and denser, and consistently more rotationally supported than in the runs including feedback, as shown by $V_\phi/\sigma_g \sim 3$, almost a factor of 2 larger than in the other cases. Moreover, while $Q_{g,turb} \sim 10$ is almost constant across the disc in run PH, it drops to $\lesssim 2$ within $R \sim 1.5$ kpc when the feedback is absent. This is due to the combined effect of the larger surface density of the disc and of the lower local support provided by both turbulence and thermal pressure, since the gas temperature across the disc typically corresponds to $c_s \sim 10 \text{ km s}^{-1}$ or lower when feedback is not included. Within the central 500 pc (i.e. the gaseous bulge), $Q_{g,turb}$ drops well below unity in run PH_NF, but it loses its physical significance because the central region is not rotationally supported.

The gas velocity dispersions are similar in the case with and without feedback but the stability properties of the discs are not. This suggests that a turbulent state with similar amplitude of turbulent motions can be achieved in different ways. When feedback is active, the disc settles in a ‘hot and turbulent’ configuration with a high Toomre parameter that results from efficient gas heating inside and beyond the disc. Alternatively, when radiative cooling effectively counterbalances heating, the disc remains in a marginally unstable ‘cold and turbulent’ state. While all these pieces of evidence do not unambiguously disentangle the relevance of stellar feedback and gravity as sources of turbulence, they do show that stellar feedback has a major impact in shaping the interstellar medium of such high-redshift galaxies. We argue that stellar feedback is possibly the main source of turbulent energy in this case. Indeed, only when we turn off stellar feedback the galaxy quickly readjusts to a new state that looks very similar to gravito-turbulence in terms of both morphology and structure of the disc, whereas when the feedback is active, the gaseous disc remains globally Toomre stable on a longer time-scale. This implies that supernova feedback affects the interstellar medium enough to prevent gravito-turbulence to take
place and in turn it suggests indirectly that feedback mainly powers the turbulent cascade in the gaseous disc.

However, we caution that the interpretation of the differences between the runs with and without feedback is complicated by other aspects: (i) the PH_NF run is restarted when the gas has already a notable amount of turbulent motions, hence its new dynamical state does not necessarily reflect only the onset of gravito-turbulence; (ii) the galaxy is still accreting both fresh material from larger scales and material recycled through the gaseous halo that likely contributes to additional turbulent motions when it joins the galactic disc (Klessen & Hennebelle 2010). Regarding (i), though, it is noteworthy that the gas cools on time-scales shorter than the orbital times when feedback is not active, which suggests that gravito-turbulence is needed to sustain the velocity dispersion in this new ‘cold’ state of the disc. In addition to that, the lack of stellar feedback makes PH_NF grows more than PH owing to the lower ‘resistance’ to gas inflows. As a consequence, the baryonic mass (in particular the gas component) within $0.1 R_{vir}$ of PH_NF is ~2 times larger than in PH, which leads to a steepening of the circular velocity curve within ~1.5 kpc; the latter peaks at ~190 km s⁻¹ and then asymptotically declines to ~160 km s⁻¹. This possibly stabilizes the disc as it increases $κ$ in the numerator of equation (5). At the same time, however, the lack of feedback increases the gas disc surface density through the enhancement of inflows and lowers the value of $V$, potentially decreasing $Q$. Since we observe a lower value of $Q$ in PH_NF, this suggests that feedback mostly influences the structure of the disc by means of both internal and external effects, namely the contribution to local turbulence and the indirect redistribution of mass through the enhancement of inflows, respectively. All these pieces of evidence are consistent with the possibility that feedback has a prominent role in shaping the disc turbulence in our simulations; none the less, this complex interplay between feedback and gravity prevents a clear identification of the ultimate cause of the turbulent cascade in the interstellar medium.

Finally, we mention that, as described in Section 3.3, we have also computed the power spectrum of the density-weighted velocity $w = \rho^{1/3}v$ in run PH in order to distinguish between a solenoidal or compressive mode of energy injection in turbulence. We find slopes typically ~-1.8 that would favour a dominant solenoidal mode, according to the results of Federrath (2013). This looks consistent with recent work showing that supernova-driven turbulence is only mildly compressive, at least on the molecular cloud scale (Pan et al. 2016), and it would support the interpretation that stellar feedback may be the main responsible for turbulence, though it is not clear whether and how this estimate is degenerate with and sensitive to the multiphase structure of the gas.

## 3.5 Mass flow through the disc

The turbulence within the disc can cause motions and mass transfer through the galactic disc. The upper panel of Fig. 12 shows the radial mass flow through the galactic disc. This is measured within a disc of 2 kpc radius and 400 pc thickness centred on the disc mid-plane. For each radial bin of width $\Delta R$ in polar coordinates, we measure the mass flow as $\dot{M} = \Delta R^{-1} \sum_j m_j v_{j,R}$, where $m_j$ and $v_{j,R}$ are the mass and the radial mid-plane velocity of the j-th gas particle, respectively, and the sum extends over the particles within each radial bin. Negative values of $\dot{M}$ are associated with inflows of gas through the galactic disc.

The figure shows that the gas is not steady inflowing towards the central region over the time span of ~160 Myr between $z = 7.6$ and $z = 6.5$. Instead, $\dot{M}$ fluctuates from negative to positive values

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Upper panel: radial mass flow through the gaseous disc. Blue solid, red dashed and green dotted lines refer to $z = 6.5, 7.1$ and $7.6$, respectively. The shaded regions with the same colours mark the symmetric inflow-outflow estimated by equation (8). Negative values of $\dot{M}$ are associated with inflows. Lower panel: radial profile of $V_\phi/\sigma_g$ at different redshifts.</figcaption></figure>

(i.e. from inflow to outflow) at different locations in the disc with typical absolute values within ~10 $M_{\odot}$ yr⁻¹. However, the gas flows in more strongly at $z = 6.5$, up to about $30 M_{\odot}$ yr⁻¹ between $R \approx 300$ pc and $\approx 2$ kpc, possibly because of the disturbance of the satellite galaxy visible in Fig. 1. We can fairly estimate the absolute value of $\dot{M}$ over the disc as

$$ \dot{M} \sim M(<R)/t_{\text{turb}} \sim M(<R)\sigma_R/R, \qquad (8) $$

where $M(<R)$ is the enclosed mass and $t_{\text{turb}} \sim R/\sigma_R$ is the turbulence crossing time, which also approximates the dissipation time of the turbulent kinetic energy if not continuously replenished (Mac Low 1999; Elmegreen 2000). However, it does not constrain the sign of $\dot{M}$, i.e. whether gas would preferentially inflow or outflow.

The overall behaviour seems qualitatively consistent with mass transport due to turbulence induced by feedback and not gravito-turbulence, as discussed above. Indeed, Goldbaum, Krumholz & Forbes (2015) have used controlled simulations of gravitationally unstable disc galaxies to show that gravito-turbulence would be able to sustain a net mass inflow over time through the disc plane because of the coherent torquing of the gas from persistent spiral arms. Specifically, they use models of Milky Way like galaxies and they find inflows ~1–2 $M_{\odot}$ yr⁻¹. The latter have lower absolute values than in our case, owing to the different conditions of the gas, i.e. different surface density, surface star formation rates, velocity dispersion, etc., but the main difference remains the steady inflows found by Goldbaum et al. (2015) in gravito-turbulent disc models. When we average $\dot{M}(R)$ over time between $z = 8.1$ and $z = 6.5$,
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Mollweide projections of the radial mass flow rate per unit steradian (negative and positive values represent inflow and outflow, respectively) at different scales and times. From top to bottom: redshift 6.5, 7.1 and 7.6, respectively. From left to right: mass flow rate in a spherical shell between 100 and 120 physical pc, between 500 and 550 physical pc, and between 1 and 1.2 physical kpc. The equator of the map corresponds with the disc mid-plane.</figcaption></figure>

we find large positive and negative fluctuations around zero, hinting against a net and continuous mass inflow due to gravito-turbulence, consistently with the previous arguments. However, we note that we do not have enough time resolution in the dumped snapshots to firmly assess the convergence of the timely averaged $\dot{M}(R)$ around zero across the entire disc.

Close to the centre, the gaseous disc becomes proportionally thicker, with an increasing aspect ratio $H/R \sim \sigma_g/V_\phi \gtrsim 1$ at radii $R \lesssim 300–500$ pc, as shown in the lower panel of Fig. 12. Therefore, the impact of three-dimensional turbulence in the mass flow close to the centre may also affect the angular distribution of moving matter and the direction of streaming gas that might eventually reach the galactic nucleus. We show that in Fig. 13 through Mollweide angular projections of the mass flow rate per unit of steradian, $d\dot{M}/d\Omega$, at different redshifts and at different radial distances from the centre. Specifically, we select gas particles within spherical test shells with radius $r$ and thickness $\Delta r$. Then, we tessellate the sphere by means of the HEALPIX algorithm<fnref n="9" /> and we calculate the mass flow through each angular tassel as $d\dot{M}/d\Omega = \Delta\Omega^{-1}\Delta r^{-1} \sum_j m_j v_{j,r}$, where $\Delta\Omega = 4\pi/N$ is the solid angle of the $N$ equal tassels, and the sum extends to the particles within each tassel.

The (instantaneous) mass flow is largely anisotropic on several scales. Most of the mass flow occurs through the disc plane, including both inflows and outflows at $\gtrsim 5$ M$_\odot$ yr$^{-1}$. On large scales $\sim 1$ kpc (the same order as the disc size), significant inflows and outflows proceed through ‘pockets’ of gas localized in solid angle from medium ($\sim 30^\circ$) to high ($\gtrsim 60^\circ$) latitudes both above and below the disc plane. This is consistent with the large-scale behaviour of the gas triggered by the interplay of stellar feedback and gravity discussed in Section 3.2; indeed, most of the inflowing gas at latitudes far from the disc plane typically has low mass-weighted mean temperature $\lesssim 10^4$ K, while the contrary is true for the outflowing gas, likely pushed away by supernova blast waves. At intermediate scales $\sim 500$ pc, the inflow-outflow episodes are mostly confined around the disc plane up to latitudes as high as $\sim 30^\circ–40^\circ$, consistently with the typical thickness of the gaseous disc $\sim 100–300$ pc as seen on the scale of the test shell, while fewer streams of gas cross the test shell at higher latitudes. On even smaller scales $\sim 100$ pc, the test shell is almost entirely embedded in the thick gaseous disc. At this scale, the dynamics of the gas within the galactic disc is dominated by turbulent motion and not by rotation, as hinted by the lower panel of Fig. 12. As therefore expected, the mass flow on small scales is highly anisotropic, with well-defined inflow-outflow regions representing cross-sections through the test shell of moving overdense gas clouds at every latitude. Consistently with Fig. 12, the integrated $d\dot{M}/d\Omega$'s shown in Fig. 13 effectively decrease from larger to smaller scales; however, the value of the net mass flow is

<footnote n="9">⁹For further information, see http://healpix.sourceforge.net/. We use the PYTHON implementation HEALPY, freely available at https://github.com/healpy/healpy.</footnote>
different because the former is integrated over a cylindrical shell selecting the flow along the disc plane, while the latter represents the radial flow through a spherical shell.

## 4 DISCUSSION AND CONCLUSIONS

In this paper, we present the results from PonosHydro, a high-resolution, zoom-in cosmological simulation meant to model the early evolution of a present-day massive galaxy down to $z \sim 6$, whose global properties appear to be consistent with the available data for galaxies of similar stellar masses at those redshifts (Iye et al. 2006; Bradley et al. 2012; Watson et al. 2015). Specifically, we study the assembly of the galaxy during its first starburst phase, before the supermassive black hole can exert significant feedback and quench star formation. We focus on the properties of the interstellar medium and the transport of mass across the galactic disc and study the conditions that determine the early evolution of the central regions of present day massive galaxies.

Before we discuss the potential implications of our findings, we briefly comment on the possible shortcomings of our calculations. Our results may be subject to the feedback model that we have used and this could in principle quantitatively affect our conclusions, at least to some extent. Various feedback schemes differently change the temperature of the gas locally after supernovae explosions. The delayed-cooling blast wave feedback produces gas at $\sim 10^5$ K and at densities $\sim 10–100$ H cm$^{-3}$ by injecting energy in the surrounding of the star-forming regions and preventing them from cooling. Before this gas expands adiabatically and gets blown away, it significantly contributes to the global stability of the galaxy against the onset of gravito-turbulence and fragmentation. While the contribution of this phase might be artificially enhanced by the feedback scheme, we argue that the physical conditions of run PH (in particular, the high specific star formation rate and the relatively small mass) are more prone to form a warm and dense gas phase than in local disc galaxies. Fig. 14 shows the density–temperature diagram of the gas within 3 and 10 physical kpc around the main galaxy in run PH at $z = 7.1$ and ErisMC$^{10}$ at $z = 3$, respectively. Despite that both simulations show some amount of gas around $10–100$ H cm$^{-3}$ and $\sim 10^5$ K, this accounts for an order of magnitude less mass fraction in the more quiescent ErisMC, that has almost 10 times smaller specific star formation rate (Shen et al. 2012). Similar results have been obtained with different feedback models (see e.g. fig. 11 of Hopkins, Quataert & Murray 2012a, where the warm phase is more sub-dominant compared to the cold phase at $\sim 100$ H cm$^{-3}$ in a Milky Way like galaxy than in a high-$z$-like galaxy). Moreover, the latest generation of stronger feedback models, designed to capture well the run between stellar mass and halo mass across cosmic scales and epochs, tend to produce lower density, thicker discs by driving more powerful and hot ($T \gtrsim 10^7$ K) gas outflow (Hopkins et al. 2014; Keller et al. 2014). Such galactic discs would likely be less gravitationally unstable and yet turbulent in the gas component (e.g. Hopkins et al. 2012a,b; Mayer et al. 2016). Therefore, we conclude that a different and stronger feedback model would likely lead

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Comparison of the phase diagram of run PH at $z = 7.1$ (upper panel) and ErisMC at $z = 3$ (lower panel). The colour bar shows the logarithm of the mass fraction per bin in the density–temperature plane. The red dotted line in the upper panel marks the density and star formation thresholds adopted in run PH.</figcaption></figure>

to qualitatively similar conclusions, namely that the gaseous discs of typical star-forming $z \sim 6$ galaxies would be maintained turbulent and stable against gravitational fragmentation by feedback, albeit future tests with different feedback schemes could better assess potential differences.

Two main processes have been often advocated in the literature to shape the global dynamics of the interstellar medium: gravitational instability and supernova feedback. In Section 3, we discussed the features of the gaseous disc of run PH in terms of star formation, outflows, turbulence, mass transport and gravitational stability, arguing that feedback likely plays a dominant role in the early evolution of a typical $z \sim 6–7$ galaxy. This seems to be somewhat different from what is usually expected both at low ($z \sim 0$) and intermediate ($z \sim 2$) redshift. Recently, Goldbaum et al. (2015) and Goldbaum, Krumholz & Forbes (2016) have thoroughly explored the relative role of gravity and stellar feedback with controlled simulations of present-day Milky Way like galaxies. Consistently with previous results (e.g. Agertz et al. 2009; Bournaud et al. 2010; Agertz, Romeo & Grisdale 2015), they find that stellar feedback is important to

<footnote>$^{10}$ ErisMC is a simulation of a Milky Way-like galaxy that adopts a sub-grid model very similar to ours and therefore it allows us to compare the thermodynamics of the gas in different physical conditions, though with the same feedback scheme. For further details about ErisMC, see Shen et al. (2012). Unfortunately, a snapshot at redshift lower than 3, which might have strengthened our considerations providing a more quiescent galaxy, is not available.</footnote>
locally regulate star formation and to create a multiphase interstellar medium, but the galaxy none the less settles to a gravito-turbulent state with a Toomre parameter $Q \sim 1$ that mostly controls the velocity dispersion and the mass transport through the disc.

Those results match the observations of nearby spiral galaxies with low star formation rates (e.g. Tamburro et al. 2009; Bagetakos et al. 2011). They are conceptually similar to what has been often argued for massive $\sim 10^{11}$ M$_{\odot}$, gas-rich galaxies at $z \approx 2$, i.e. star-forming galaxies that have not been quenched yet and are likely the progenitors of the most massive quiescent galaxies at $z=0$, undergoing so-called violent disc instability, i.e. gravitational instability that leads to the formation of massive star-forming clumps (e.g. Dekel, Sari & Ceverino 2009; Ceverino, Dekel & Bournaud 2010; Mandelker et al. 2014; Inoue et al. 2016; but see also Hopkins et al. 2012b; Tamburello et al. 2015). In this respect, massive discs at $z \approx 2$ would be the most extreme manifestation of the 'cold and turbulent' regime that eventually leads to gravitational instability and fragmentation in the gas, since they are already as massive as the most massive discs in the local Universe but proportionally more gas rich, perhaps because they appear near the peak of the cosmic star formation history (Madau & Dickinson 2014). Motivated by this analogy, Goldbaum et al. (2016) have proposed that gravitational instability is the dominant process setting mass transport and fuelling star formation over cosmic time. However, our results suggest that this might not be the case for typical $z \sim 6-7$ galaxies with stellar and gas mass $\gtrsim 10^9$ M$_{\odot}$, where the disc dynamics and mass transport seems to be significantly influenced by stellar feedback. Indeed, even in the favourable conditions of massive gaseous discs at $z \approx 2$, a phase of violent disc instability with its associated gravito-turbulent can or cannot occur depending on how effective is the feedback model at heating the gas and generating mass-loaded outflows. Recent work has shown that modern strong feedback models tend to suppress gravitational instability and fragmentation, and at the same time that blast wave feedback cannot suppress disc instability when conditions are favourable for its emergence (Mayer et al. 2016), corroborating at least the qualitative distinction between the two cases. Interestingly, Ceverino et al. (2017) recently reach analogous conclusions on the role of stellar feedback by looking at a galaxy with stellar mass $M_* \sim 10^9$ M$_{\odot}$ comparable to ours but at $z \sim 1$. They find a qualitatively similar evolution over time of the Toomre parameter and the $V_{\phi}/\sigma_g$ ratio despite they use a different approach to model stellar feedback, i.e. through non-thermal radiation pressure in addition to thermal dump without shut-off cooling.

We thus argue that the dominant role of stellar feedback in the early evolutionary phase of massive galaxy progenitors is likely controlled by the combination of the high specific star formation rate ($\gtrsim 5$ Gyr$^{-1}$) and of the relatively low mass at $z > 5$ ($\sim 10^9$ M$_{\odot}$). The first favours the impact of stellar feedback on the interstellar medium, while the latter proportionally weakens the dynamical role of gravity to lead to instabilities and eventually fragmentation. Those specific star formation rates are expected for galaxies on the main star-forming sequence at $z > 5$ (Fig. 2; Schreiber et al. 2015; Tasca et al. 2015); as our galaxy is consistent with the main sequence, this suggests that the 'hot and turbulent' regime that we characterize here could be typical of star-forming galaxies at $z > 5$ with baryonic/stellar masses comparable to ours. These should be fairly typical galaxies, as recent surveys begin to find (Bradley et al. 2012; Capak et al. 2015; Watson et al. 2015). In particular, recent ALMA observations by Maiolino et al. (2015) tend to qualitatively support the idea that stellar feedback has a dominant role in the early assembly of normal star-forming galaxies at $z \sim 6-7$. This

has immediate observational implications, as we would predict a significant amount of warm/hot gas with temperature $5 \times 10^4 \lesssim T/K \lesssim 5 \times 10^5$ inside and around the disc, possibly $\sim 0.1$ of the gas mass. Note that these are temperatures more akin to the circumgalactic medium distributed in the virial volume around galaxies at low redshift (Werk et al. 2013), but in our case it would be inside or surrounding the galactic disc.

Possible analogues in the local Universe may serve as preliminary test bed for our predictions. Those might be low-mass starburst galaxies, such as the prototypical M82, that have stellar mass and star formation rate rather similar to the main galaxy of run PH (with a factor 3–4 lower specific star formation rate due to the higher stellar mass; e.g. Förster Schreiber et al. 2003; Greco, Martini & Thompson 2012). While gas densities are expected to be lower at $z=0$, M82 might also host a significant fraction of warm/hot, turbulent gas in its disc, at least in the central kpc where the starburst is ongoing. This seems to be confirmed by observations (e.g. Griffiths et al. 2000) and at least in qualitative agreement with our results since gas temperature and phases found in our simulations would be somewhat dependent on the specific feedback model. Detailed characterization of the warm/hot interstellar medium in low-mass starburst galaxies could thus provide useful constraints to test our scenario.

In this 'hot and turbulent' regime, the mass transport is influenced by intense and clustered stellar feedback episodes. As a result, the gas flow through the disc is fluctuating and anisotropic, with no sustained coherent gas inflow within the disc. A coherent circumnuclear disc, which can be a way to funnel accretion towards the ultimate stage of the accretion disc, is not clearly seen to form at $\sim 100$ pc scales (though this would be barely resolved at our resolution). This might have some implications for the feeding of a central massive black hole (Gabor & Bournaud 2013; Dubois et al. 2015). On one hand, mean inflow rates could be small; however, episodic accretion events at high rates could occur through the infall of massive gas clouds, as we observe inflow rates that can occasionally peak at $\gtrsim 5$ M$_{\odot}$ yr$^{-1}$. None the less, if super-Eddington accretion is assumed, recent models show that even episodic accretion is enough for the rapid growth of central black holes (e.g. Lupi et al. 2016; Pezzulli, Valiante & Schneider 2016). Accretion may also occur in an anisotropic way, with large fluctuations in the angular momentum of the accreting matter, which would have implications for the nature of the accretion disc itself, if any, and for the evolution of the spin of the central black hole. However, we defer additional speculations on the evolution of a massive black hole in such environments to a forthcoming investigation.

As steady central gas inflows are not sustained, bulge/spheroid formation from dynamical and/or secular disc instabilities are unlikely to take place (e.g. Guedes et al. 2013). Indeed, the disc of PonosHydro remains nearly bulgeless for the whole simulation (see Fig. 4). However, we know from lower resolution runs going to $z=0$ that the galaxy will develop a dominant spheroid at lower redshifts as it grows to become a massive early-type galaxy (Fiacconi et al., in preparation). Since it has several mergers occurring at later times ($2 < z < 4$; Fiacconi et al. 2016), it is likely that such mergers will be the dominant driver of spheroid growth (Fiacconi et al. 2015). However, if the 'hot and turbulent' regimes characterize main-sequence galaxies at $z > 5$, this would imply that bulge formation may occur after the first billion year of evolution, possibly post-dating the growth of the massive black hole at the centre (see also Dubois et al. 2015; Habouzit, Volonteri & Dubois 2016). Therefore, we predict that gas-rich star-forming discs at $z > 5$ should not host a significant bulge.
The exploration described in this paper leads to interesting predictions about the early assembly of massive galaxies. However, our interpretations remain rather speculative, both from the theoretical and the observational point of view. On one hand, future simulations including different sub-grid models are necessary to quantitatively assess the different nature of galaxies at low and high redshifts, also comparing the results from codes with intrinsically different treatments of hydrodynamics (e.g. Kim et al. 2016). On the other hand, more detailed characterization of high-redshift galaxies are starting to be available from current observational facilities (e.g. ALMA). However, it is going to be in the next future that forthcoming observatories (e.g. JWST, E-ELT) will provide deep enough data to definitely test our predictions and at the same time to better guide the theoretical study of galaxies at the cosmic dawn.

## ACKNOWLEDGEMENTS

We thank the anonymous referee for constructive comments that helped us to improve the quality of the paper. We acknowledge useful discussions with Arif Babul, Rychard Bouwens, Nick Gnedin, Raffaella Schneider, Sijing Shen and Debora Sijacki. We thank Oliver Hahn and the AGORA collaboration for help with the ICs of the simulations. The simulations have been run on the ZBOX4 cluster at the University of Zurich and on the Pitz Dora cluster at CSCS, Lugano. We acknowledge the use of the PYTHON package PYNBODY (Pontzen et al. 2013, publicly available at https://github.com/pynbody/pynbody) in our analysis for this paper. DF is supported by the Swiss National Science Foundation under grant #No. 200021_140645. DF also acknowledges support by European Research Council (ERC) Starting Grant 638707 'Black holes and their host galaxies: coevolution across cosmic time'. Support for this work was provided to PM by the National Science Foundation (NSF) through grant AST-1229745, and by National Aeronautics and Space Administration (NASA) through grant NNX12AF87G. PM also acknowledges a NASA contract supporting the WFIRST-EXPO Science Investigation Team (15-WFIRST15-0004), administered by Goddard Space Flight Center (GSFC), and thanks the Préfecture of the Ile-de-France Region for the award of a Blaise Pascal International Research Chair, managed by the Fondation de l'Ecole Normale Supérieure.

## REFERENCES

* Agertz O. et al., 2007, MNRAS, 380, 963

* Agertz O., Lake G., Teyssier R., Moore B., Mayer L., Romeo A. B., 2009, MNRAS, 392, 294

* Agertz O., Teyssier R., Moore B., 2011, MNRAS, 410, 1391

* Agertz O., Romeo A. B., Grisdale K., 2015, MNRAS, 449, 2156

* Bagetakos I., Brinks E., Walter F., de Blok W. J. G., Usero A., Leroy A. K., Rich J. W., Kennicutt R. C., Jr, 2011, AJ, 141, 23

* Baldry I. K. et al., 2010, MNRAS, 404, 86

* Bate M. R., Burkert A., 1997, MNRAS, 288, 1060

* Behroozi P. S., Wechsler R. H., Conroy C., 2013, ApJ, 770, 57

* Bigiel F., Leroy A., Walter F., Brinks E., de Blok W. J. G., Madore B., Thornley M. D., 2008, AJ, 136, 2846

* Blumenthal G. R., Faber S. M., Primack J. R., Rees M. J., 1984, Nature, 311, 517

* Bournaud F., Elmegreen B. G., Teyssier R., Block D. L., Puerari I., 2010, MNRAS, 409, 1088

* Bournaud F. et al., 2011, ApJ, 730, 4

* Bouwens R. J. et al., 2003, ApJ, 595, 589

* Bouwens R. J., Illingworth G. D., Franx M., Ford H., 2007, ApJ, 670, 928

* Bouwens R. J. et al., 2011, ApJ, 737, 90

* Bouwens R. J. et al., 2015, ApJ, 803, 34

* Bower R. G., Benson A. J., Malbon R., Helly J. C., Frenk C. S., Baugh C. M., Cole S., Lacey C. G., 2006, MNRAS, 370, 645

* Boylan-Kolchin M., Ma C.-P., Quataert E., 2008, MNRAS, 383, 93

* Bradley L. D. et al., 2012, ApJ, 747, 3

* Bressan A., Marigo P., Girardi L., Salasnich B., Dal Cero C., Rubele S., Nanni A., 2012, MNRAS, 427, 127

* Brunt C. M., Heyer M. H., Mac Low M.-M. 2009, A&A, 504, 883

* Bryan G. L., Norman M. L., 1998, ApJ, 495, 80

* Capak P. L. et al., 2015, Nature, 522, 455

* Casey C. M., Narayanan D., Cooray A., 2014, Phys. Rep., 541, 45

* Ceverino D., Dekel A., Bournaud F., 2010, MNRAS, 404, 2151

* Ceverino D., Primack J., Dekel A., Kassin S. A., 2017, MNRAS, 467, 2664

* Chapman S. C., Blain A. W., Smail I., Ivison R. J., 2005, ApJ, 622, 772

* Choi J.-H., Nagamine K., 2012, MNRAS, 419, 1280

* Cole S., Aragon-Salamanca A., Frenk C. S., Navarro J. F., Zepf S. E., 1994, MNRAS, 271, 781

* Costa T., Sijacki D., Haehnelt M. G., 2015, MNRAS, 448, L30

* Crain R. A. et al., 2009, MNRAS, 399, 1773

* Croton D. J. et al., 2006, MNRAS, 365, 11

* D'Onghia E., Vogelsberger M., Faucher-Giguere C.-A., Hernquist L., 2010, ApJ, 725, 353

* Daddi E. et al., 2007, ApJ, 670, 156

* Daddi E. et al., 2010, ApJ, 713, 686

* Dalla Vecchia C., Schaye J., 2008, MNRAS, 387, 1431

* De Breuck C. et al., 2014, A&A, 565, A59

* Dekel A., Sari R., Ceverino D., 2009, ApJ, 703, 785

* Diemand J., Kuhlen M., Madau P., Zemp M., Moore B., Potter D., Stadel J., 2008, Nature, 454, 735

* Dubois Y. et al., 2014, MNRAS, 444, 1453

* Dubois Y., Volonteri M., Silk J., Devriendt J., Slyz A., Teyssier R., 2015, MNRAS, 452, 1502

* Elmegreen B. G., 2000, ApJ, 530, 277

* Elmegreen B. G., Kim S., Staveley-Smith L., 2001, ApJ, 548, 749

* Federrath C., 2013, MNRAS, 436, 1245

* Federrath C., Glover S. C. O., Klessen R. S., Schmidt W., 2008, Phys. Scr. T, 132, 014025

* Feldmann R., Mayer L., 2015, MNRAS, 446, 1939

* Feldmann R., Carollo C. M., Mayer L., Renzini A., Lake G., Quinn T., Stinson G. S., Yepes G., 2010, ApJ, 709, 218

* Ferrière K. M., 2001, Rev. Mod. Phys., 73, 1031

* Fiacconi D., Feldmann R., Mayer L., 2015, MNRAS, 446, 1957

* Fiacconi D., Madau P., Potter D., Stadel J., 2016, ApJ, 824, 144

* Förster Schreiber N. M., Genzel R., Lutz D., Sternberg A., 2003, ApJ, 599, 193

* Förster Schreiber N. M. et al., 2009, ApJ, 706, 1364

* Gabor J. M., Bournaud F., 2013, MNRAS, 434, 606

* Gallerani S., Pallottini A., Feruglio C., Ferrara A., Maiolino R., Vallini L., Riechers D. A., 2016, preprint (arXiv:1604.05714)

* Gammie C. F., 2001, ApJ, 553, 174

* Genzel R. et al., 2006, Nature, 442, 786

* Gill S. P. D., Knebe A., Gibson B. K., 2004, MNRAS, 351, 399

* Gingold R. A., Monaghan J. J., 1977, MNRAS, 181, 375

* Girardi L. et al., 2010, ApJ, 724, 1030

* Gnedin N. Y., Tassis K., Kravtsov A. V., 2009, ApJ, 697, 55

* Gnedin N. Y., Tasker E. J., Fujimoto Y., 2014, ApJ, 787, L7

* Goldbaum N. J., Krumholz M. R., Forbes J. C., 2015, ApJ, 814, 131

* Goldbaum N. J., Krumholz M. R., Forbes J. C., 2016, ApJ, 827, 28

* Gómez G. C., Cox D. P., 2002, ApJ, 580, 235

* Governato F. et al., 2010, Nature, 463, 203

* Greco J. P., Martini P., Thompson T. A., 2012, ApJ, 757, 24

* Green A. W. et al., 2010, Nature, 467, 684

* Griffiths R. E., Ptak A., Feigelson E. D., Garmire G., Townsley L., Brandt W. N., Sambruna R., Bregman J. N., 2000, Science, 290, 1325

* Grogin N. A. et al., 2011, ApJS, 197, 35

* Guedes J., Callegari S., Madau P., Mayer L., 2011, ApJ, 742, 76

* Guedes J., Mayer L., Carollo M., Madau P., 2013, ApJ, 772, 36
Guo Q. et al., 2011, MNRAS, 413, 101

Haardt F., Madau P., 2012, ApJ, 746, 125

Habouzit M., Volonteri M., Dubois Y., 2016, MNRAS, preprint (arXiv:1605.09394)

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

Heckman T. M., Alexandroff R. M., Borthakur S., Overzier R., Leitherer C., 2015, ApJ, 809, 147

Henriques B. M. B., White S. D. M., Thomas P. A., Angulo R. E., Guo Q., Lemson G., Springel V., 2013, MNRAS, 431, 3373

Hernquist L., Katz N., 1989, ApJS, 70, 419

Hinshaw G. et al., 2013, ApJS, 208, 19

Hopkins P. F., 2013, MNRAS, 428, 2840

Hopkins P. F. et al., 2010, ApJ, 724, 915

Hopkins P. F., Quataert E., Murray N., 2012a, MNRAS, 421, 3488

Hopkins P. F., Kereš D., Murray N., Quataert E., Hernquist L., 2012b, MNRAS, 427, 968

Hopkins P. F., Kereš D., Oñorbe J., Faucher-Giguère C.-A., Quataert E., Murray N., Bullock J. S., 2014, MNRAS, 445, 581

Inoue S., Dekel A., Mandelker N., Ceverino D., Bournaud F., Primack J., 2016, MNRAS, 456, 2052

Iye M. et al., 2006, Nature, 443, 186

Jiang C. Y., Jing Y. P., Faltenbacher A., Lin W. P., Li C., 2008, ApJ, 675, 1095

Joung M. K. R., Mac Low M.-M., 2006, ApJ, 653, 1266

Kauffmann G., White S. D. M., Guiderdoni B., 1993, MNRAS, 264, 201

Kaufmann T., Mayer L., Wadsley J., Stadel J., Moore B., 2007, MNRAS, 375, 53

Keller B. W., Wadsley J., Benincasa S. M., Couchman H. M. P., 2014, MNRAS, 442, 3013

Kennicutt R. C., Jr, 1998, ApJ, 498, 541

Kim J.-h. et al., 2014, ApJS, 210, 14

Kim J.-h. et al., 2016, ApJ, 833, 202

Klessen R. S., Hennebelle P., 2010, A&A, 520, A17

Klypin A. A., Trujillo-Gomez S., Primack J., 2011, ApJ, 740, 102

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

Kochanek C. S. et al., 2012, ApJS, 200, 8

Koekemoer A. M. et al., 2011, ApJS, 197, 36

Komatsu E. et al., 2011, ApJS, 192, 18

Kritsuk A. G., Norman M. L., Wagner R., 2011, ApJ, 727, L20

Krivitsky D. S., Kontorovich V. M., 1997, A&A, 327, 921

Kroupa P., 2001, MNRAS, 322, 231

Le Fèvre O. et al., 2005, A&A, 439, 845

Lilly S. J. et al., 2007, ApJS, 172, 70

Lodato G., Clarke C. J., 2011, MNRAS, 413, 2735

Lucy L. B., 1977, AJ, 82, 1013

Lupi A., Haardt F., Dotti M., Fiacconi D., Mayer L., Madau P., 2016, MNRAS, 456, 2993

McKee C. F., Ostriker J. P., 1977, ApJ, 218, 148

Mac Low M.-M., 1999, ApJ, 524, 169

Madau P., Dickinson M., 2014, ARA&A, 52, 415

Madau P., Ferguson H. C., Dickinson M. E., Giavalisco M., Steidel C. C., Fruchter A., 1996, MNRAS, 283, 1388

Maiolino R. et al., 2015, MNRAS, 452, 54

Mandelker N., Dekel A., Ceverino D., Tweed D., Moody C. E., Primack J., 2014, MNRAS, 443, 3675

Marigo P., Girardi L., Bressan A., Groenewegen M. A. T., Silva L., Granato G. L., 2008, A&A, 482, 883

Mayer L., Tamburello V., Lupi A., Keller B., Wadsley J., Madau P., 2016, ApJ, 830, L13

Meru F., Bate M. R., 2011a, MNRAS, 410, 559

Meru F., Bate M. R., 2011b, MNRAS, 411, L1

Mo H. J., Mao S., White S. D. M., 1998, MNRAS, 295, 319

Muratov A. L., Kereš D., Faucher-Giguère C.-A., Hopkins P. F., Quataert E., Murray N., 2015, MNRAS, 454, 2691

O'Shea B. W., Wise J. H., Xu H., Norman M. L., 2015, ApJ, 807, L12

Ocvirk P. et al., 2016, MNRAS, 463, 1462

Oesch P. A. et al., 2010a, ApJ, 709, L16

Oesch P. A. et al., 2010b, ApJ, 709, L21

Oesch P. A. et al., 2012, ApJ, 759, 135

Oesch P. A. et al., 2014, ApJ, 786, 108

Paardekooper S.-J., Baruteau C., Meru F., 2011, MNRAS, 416, L65

Paardekooper J.-P., Khochfar S., Dalla Vecchia C., 2015, MNRAS, 451, 2544

Padoan P., Nordlund A., Jones B. J. T., 1997, MNRAS, 288, 145

Pan L., Padoan P., Haugbolle T., Nordlund A., 2016, ApJ, 825, 30

Percival W. J. et al., 2001, MNRAS, 327, 1297

Percival W. J. et al., 2010, MNRAS, 401, 2148

Pezzulli E., Valiante R., Schneider R., 2016, MNRAS, 458, 3047

Pontzen A., Roškar R., Stinson G., Woods R., 2013, Astrophysics Source Code Library, record ascl:1305.002

Price D. J., Federrath C., 2010, MNRAS, 406, 1659

Raiteri C. M., Villata M., Navarro J. F., 1996, A&A, 315, 105

Read J. I., Hayfield T., 2012, MNRAS, 422, 3037

Reed D. S., Bower R., Frenk C. S., Jenkins A., Theuns T., 2007, MNRAS, 374, 2

Reid B. A. et al., 2010, MNRAS, 404, 60

Richings A. J., Schaye J., 2016, MNRAS, 458, 270

Robertson B. E., Kravtsov A. V., 2008, ApJ, 680, 1083

Romeo A. B., 1994, A&A, 286, 799

Romeo A. B., Falstad N., 2013, MNRAS, 433, 1389

Romeo A. B., Wiegert J., 2011, MNRAS, 416, 1191

Roškar R., Fiacconi D., Mayer L., Kazantzidis S., Quinn T. R., Wadsley J., 2015, MNRAS, 449, 494

Scalo J., Vázquez-Semadeni E., Chappell D., Passot T., 1998, ApJ, 504, 835

Schaye J. et al., 2010, MNRAS, 402, 1536

Schaye J. et al., 2015, MNRAS, 446, 521

Schreiber C. et al., 2015, A&A, 575, A74

Schroetter I., Bouché N., Péroux C., Murphy M. T., Contini T., Finley H., 2015, ApJ, 804, 83

Shen S., Wadsley J., Stinson G., 2010, MNRAS, 407, 1581

Shen S., Madau P., Aguirre A., Guedes J., Mayer L., Wadsley J., 2012, ApJ, 760, 50

Shen S., Madau P., Guedes J., Mayer L., Prochaska J. X., Wadsley J., 2013, ApJ, 765, 89

Smit R., Bouwens R. J., Franx M., Illingworth G. D., Labbé I., Oesch P. A., van Dokkum P. G., 2012, ApJ, 756, 14

Somerville R. S., Primack J. R., 1999, MNRAS, 310, 1087

Springel V. et al., 2005, Nature, 435, 629

Stadel J. G., 2001, PhD thesis, Univ. Washington

Stadel J., 2013, Astrophysics Source Code Library, record ascl:1305.005

Steidel C. C., Adelberger K. L., Giavalisco M., Dickinson M., Pettini M., 1999, ApJ, 519, 1

Steidel C. C., Adelberger K. L., Shapley A. E., Pettini M., Dickinson M., Giavalisco M., 2003, ApJ, 592, 728

Stinson G., Seth A., Katz N., Wadsley J., Governato F., Quinn T., 2006, MNRAS, 373, 1074

Tacconi L. J. et al., 2010, Nature, 463, 781

Tamburello V., Mayer L., Shen S., Wadsley J., 2015, MNRAS, 453, 2490

Tamburro D., Rix H.-W., Leroy A. K., Mac Low M.-M., Walter F., Kennicutt R. C., Brinks E., de Blok W. J. G., 2009, AJ, 137, 4424

Tasca L. A. M. et al., 2015, A&A, 581, A54

Viel M., Bolton J. S., Haehnelt M. G., 2009, MNRAS, 399, L39

Vogelsberger M., Genel S., Sijacki D., Torrey P., Springel V., Hernquist L., 2013, MNRAS, 436, 3031

Vogelsberger M. et al., 2014, MNRAS, 444, 1518

Wadsley J. W., Stadel J., Quinn T., 2004, New Astron., 9, 137

Wadsley J. W., Veeravalli G., Couchman H. M. P., 2008, MNRAS, 387, 427

Walter F., Brinks E., de Blok W. J. G., Bigiel F., Kennicutt R. C. Jr, Thornley M. D., Leroy A., 2008, AJ, 136, 2563

Waters D., Wilkins S. M., Di Matteo T., Feng Y., Croft R., Nagai D., 2016, MNRAS, 461, L51

Watson D., Christensen L., Knudsen K. K., Richard J., Gallazzi A., Michałowski M. J., 2015, Nature, 519, 327

Weidemann V., 1987, A&A, 188, 74
Werk J. K., Prochaska J. X., Thom C., Tumlinson J., Tripp T. M., O'Meara
J. M., Peebles M. S., 2013, ApJS, 204, 17

White S. D. M., Rees M. J., 1978, MNRAS, 183, 341

Wiersma R. P. C., Schaye J., Theuns T., Dalla Vecchia C., Tornatore L., 2009, MNRAS, 399, 574

Wisnioski E. et al., 2015, ApJ, 799, 209

York D. G. et al., 2000, AJ, 120, 1579

Younger J. D. et al., 2009, ApJ, 704, 803

APPENDIX A: RESOLUTION TESTS

We run a twin simulation at lower resolution dubbed ‘PH_LR’, as described in Section 2.2. Fig. A1 shows a comparison of the surface densities of the galactic disc in runs PH and PH_LR at z = 7.6, 7.1 and 6.5. The two runs produce galaxies that look qualitatively similar, at least between z = 7.6 and z = 7.1. However, we systematically observe that the gas in run PH_LR is slightly more compact and concentrated at the centre, in qualitative agreement with previous results from controlled simulations to study the angular momentum transport in cooling haloes as a function of resolution (Kaufmann et al. 2007). Moreover, the interstellar medium shows overall smoother features than in PH, as expected given the lower resolution. The major differences are at z = 6.5, when the two galaxies also react slightly differently to the perturbation induced by a nearby companion that is going to merge with the main galaxy.

More quantitatively, Fig. A2 shows the evolution with redshift of both the virial mass, $M_{vir}$, and the total stellar mass (main galaxy and satellites) within the virial radius, $M_*(<R_{vir})$. The growth of $M_{vir}$ with time is almost indistinguishable between the two runs, while the stellar mass shows mild differences at $z > 9$, with $M_*(<R_{vir})$ in run PH_LR being roughly a factor ≲2 smaller than in run PH, but this is also possibly related to the evolution of the satellites. Similarly, Fig. A3 shows the comparison of the (specific) star formation

<figure><img src="image_17.png" /><figcaption>**Figure A2.** Comparison of the redshift evolution of the virial mass $M_{vir}$ and of the total stellar mass $M_*$ within $R_{vir}$ (i.e. main galaxy and sub-haloes) for runs PH (blue solid line) and PH_LR (red dashed line).</figcaption></figure>

rate as a function of time for the main galaxy. The two runs show little differences at the two resolutions, suggesting that the combination of the sub-grid models with the evolution of the galaxies might have reached reasonable convergence.

Finally, Fig. A4 compares the surface density power spectra of run PH_LR with the fiducial one (see Section 3.3 for further details). The power spectra of the low-resolution case also shows a similar behaviour with two nearly power law for $k$ smaller and larger than ~20–30 kpc⁻¹; however, the exponents tend to be

<figure><img src="image_18.png" /><figcaption>**Figure A1.** Comparison of the surface density of the gas within the galactic disc between runs PH (upper row) and PH_LR (bottom row). From left to right, we show the data at $z = 6.5, 7.1$ and $7.6$, respectively. The white circle at the bottom-right corner of each panel shows the size of the softening in physical units.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure A3.** Comparison of the evolution of the star formation rate (thin lines, left y-axes) and of the specific star formation rate (thick lines, right y-axes) between runs PH (blue solid line) and PHLR (red dashed line).</figcaption></figure>

systematically slightly lower (i.e. steeper power law) for PHLR. This difference is more significant for large *k* modes, where the low-resolution run shows less power than the fiducial run. This is likely related to the lower mass and spatial resolution that tends to suppress overdensities at small scales. On the other hand, the largest scales ($k < 10$ kpc$^{-1}$) seem to reach reasonable convergence, at least at $z = 7.1$ and 7.6, while at $z = 6.5$ the difference is larger, possibly because of the different response of the main galaxy in the two runs to the perturbing satellite visible in Fig. 1. Despite these differences, the qualitative agreement in the shape of $P_Σ$ for the high- and low-resolution simulations suggests that the interstellar medium reaches

<figure><img src="image_20.png" /><figcaption>**Figure A4.** Comparison of the normalized surface density power spectra for runs PH (thick lines) and PHLR (thin lines). Blue solid, red dashed and green dotted lines refer to $z = 6.5, 7.1$ and $7.6$, respectively. The left and right grey shaded regions mark the scale *k* associated with the box size ($L_{box} = 4$ kpc) and to the Nyquist frequency associated with the grid spacing (see Section 3.3 for additional details).</figcaption></figure>

a similar configuration in the two cases. In addition, also the typical value of the gas velocity dispersion $\sigma_{gas} \sim 40$ km s$^{-1}$ is comparable among the two runs. However, due to the lower resolution (i.e. the lower number of gas particles in the galactic disc), we cannot compare the velocity power spectra in 512 pc cubic boxes directly.

This paper has been typeset from a TeX/LATeX file prepared by the author.