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

# The impact of chemistry on the structure of high-z galaxies

A. Pallottini,<fnref n="1" />,<fnref n="2" />,<fnref n="3" />,<fnref n="4" />★ A. Ferrara,<fnref n="4" />,<fnref n="5" /> S. Bovino,<fnref n="6" /> L. Vallini,<fnref n="7" /> S. Gallerani,<fnref n="4" /> R. Maiolino<fnref n="2" />,<fnref n="3" /> and S. Salvadori<fnref n="8" />,<fnref n="9" />,<fnref n="10" />

Accepted 2017 July 13. Received 2017 July 12; in original form 2017 May 19

## ABSTRACT

To improve our understanding of high-z galaxies, we study the impact of H₂ chemistry on their evolution, morphology and observed properties. We compare two zoom-in high-resolution (30 pc) simulations of prototypical M\* ~ 10¹⁰ M☉ galaxies at z = 6. The first, ‘Dahlia’, adopts an equilibrium model for H₂ formation, while the second, ‘Althæa’, features an improved non-equilibrium chemistry network. The star formation rate (SFR) of the two galaxies is similar (within 50 per cent), and increases with time reaching values close to 100 M☉ yr⁻<fnref n="1" /> at z = 6. They both have SFR-stellar mass relation consistent with observations, and a specific SFR of ~5 Gyr⁻<fnref n="1" />. The main differences arise in the gas properties. The non-equilibrium chemistry determines the H → H₂ transition to occur at densities >300 cm⁻<fnref n="3" />, i.e. about 10 times larger than predicted by the equilibrium model used for Dahlia. As a result, Althæa features a more clumpy and fragmented morphology, in turn making SN feedback more effective. Also, because of the lower density and weaker feedback, Dahlia sits 3σ away from the Schmidt-Kennicutt relation; Althæa, instead nicely agrees with observations. The different gas properties result in widely different observables. Althæa outshines Dahlia by a factor of 7 (15) in [CII]157.74 µm (H₂17.03 µm) line emission. Yet, Althæa is underluminous with respect to the locally observed [CII]–SFR relation. Whether this relation does not apply at high-z or the line luminosity is reduced by cosmic microwave background and metallicity effects remain as an open question.

**Key words:** methods: numerical – galaxies: evolution – galaxies: formation – galaxies: high-redshift – galaxies: ISM – infrared: general.

## 1 INTRODUCTION

Understanding the properties of the interstellar medium (ISM) of primeval galaxies is a fundamental challenge of physical cosmology. The high sensitivity/spatial resolution allowed by current observations have dramatically improved our understanding of the ISM of local and moderate redshift (z = 2–3) galaxies (Osterbrock 1989; Stasińska 2007; Pérez-Montero 2017; Stanway 2017). We now have a clearer picture of the gas phases and thermodynamics (Daddi et al. 2010a; Carilli & Walter 2013), particularly for what concerns

the molecular component, representing the stellar birth environment (Krumholz 2015; Klessen & Glover 2016).

For galaxies located in the Epoch of Reionization (5 ≲ z ≲ 15) optical/near-infrared (IR) surveys have been very successful in their identification and characterization in terms of stellar mass and star formation rate (SFR; Dunlop 2013; Madau & Dickinson 2014; Bouwens et al. 2015). However, only recently we have started to probe the internal structure of such objects. With the advent of the Atacama Large Millimeter/Submillimeter Array (ALMA) it is now possible to access the far-infrared (FIR) band at high-z with an unprecedented resolution and sensitivity. Excitingly, this enables for the first time studies of ISM energetics, structure and composition in such pristine objects.

<footnote n="1">¹Centro Fermi, Museo Storico della Fisica e Centro Studi e Ricerche 'Enrico Fermi', Piazza delViminale 1, I-00184 Roma, Italy</footnote>

<footnote n="2">²Cavendish Laboratory, University of Cambridge, 19 J. J. Thomson Ave., Cambridge CB3 0HE, UK</footnote>

<footnote n="3">³Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="4">⁴Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy</footnote>

<footnote n="5">⁵Kavli IPMU, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa 277-8583, Japan</footnote>

<footnote n="6">⁶Hamburger Sternwarte, Universität Hamburg, Gojenbergsweg 112, D-21029 Hamburg, Germany</footnote>

<footnote n="7">⁷Nordita, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, SE-10691 Stockholm, Sweden</footnote>

<footnote n="8">⁸Dipartimento di Fisica e Astronomia, Universita' di Firenze, Via G. Sansone 1, I-50019 Sesto Fiorentino, Italy</footnote>

<footnote n="9">⁹INAF/Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125 Firenze, Italy</footnote>

<footnote n="10">¹⁰GEPI, Observatoire de Paris, PSL Research University, CNRS, Place Jule Janssen, F-92190 Meudon, France</footnote>

<footnote>* E-mail: andrea.pallottini@sns.it, andrea.pallottini@centrofermi.it</footnote>
Since C II is one of the major coolant of the ISM, [C II] ALMA detections (and upper limits) have so far mostly used this line for the above purposes (Capak et al. 2015; Maiolino et al. 2015; Willott et al. 2015) and to determine the sizes of early galaxies (Fujimoto et al. 2017). Line emission from different species (e.g. [O III]) has been used to derive the interstellar radiation field (ISRF) intensity (Inoue et al. 2016; Carniani et al. 2017), while continuum detections give us a measure of the dust content and properties (Watson et al. 2015; Laporte et al. 2017). Finally, some observations are beginning to resolve different ISM components and their dynamics by detecting spatial offsets and kinematic shifts between different emission lines, i.e. [C II] and optical–ultraviolet (UV) emission (Capak et al. 2015; Maiolino et al. 2015), [C II] and Ly α (Pentericci et al. 2016; Bradac et al. 2017) and [C II] and [O III] (Carniani et al. 2017).

In spite of these progresses, several pressing questions remain unanswered. A partial list includes the following: (a) What is the chemical composition and thermodynamic state of the ISM in high-z galaxies? (b) How does the molecular gas turns into stars and regulate the evolution of these systems? (c) What are the optimal observational strategies to better constrain the properties of these primeval objects?

Theoretically, cosmological numerical simulations have been used to attack some of these problems. The key idea is to produce a coherent physical framework within which the observed properties can be understood. Such learning strategy is also of fundamental importance to devise efficient observations from current (e.g. HST/ALMA), planned (JWST) and proposed (SPICA) instruments. Before this strategy can be implemented, though, it is necessary to develop reliable numerical schemes catching all the relevant physical processes. While the overall performances of the most widely used schemes have been extensively benchmarked (Scannapieco et al. 2012; Kim et al. 2014, 2016), high-resolution simulations of galaxy formation introduce a new challenge: they are very sensitive to the implemented physical models, particularly those acting on small scales.

Among these, the role of feedback, i.e. how stars affect their own formation history via energy injection in the surrounding gas by supernova (SN) explosions, stellar winds and radiation, is far from being completely understood, despite considerable efforts to improve its modelling (Agertz & Kravtsov 2015; Martizzi, Faucher-Giguère & Quataert 2015) and understand its consequences on high-z galaxy evolution (Ceverino et al. 2014; Barai et al. 2015; O'Shea et al. 2015; Fiacconi et al. 2017; Hopkins et al. 2017; Pallottini et al. 2017).

Additionally, we are still lacking a completely self-consistent treatment of radiation transfer. This is an area in which intensive work is ongoing in terms of faster numerical schemes (Wise et al. 2012; Rosdahl et al. 2015; Katz et al. 2016), or improved physical modelling (Petkova & Maio 2012; Roskar et al. 2014; Maio et al. 2016).

A third aspect has received comparatively less attention so far in high-z galaxy formation studies, i.e. the implementation of adequate chemical networks. While various models have been proposed and tested (Krumholz, McKee & Tumlinson 2009; Bovino et al. 2016; Grassi et al. 2017), the galaxy-scale consequences of the different prescriptions are still largely unexplored (Tomassetti et al. 2015; Maio & Tescari 2015; Smith et al. 2017). Besides, there is no clear consensus on a minimal set of physical ingredients required to produce reliable simulations.

The purpose of this paper is to analyse the impact of H₂ chemistry on the internal structure of high-z galaxies. To this aim, we simulate

two prototypical $M_* \simeq 10^{10} M_\odot$ Lyman break galaxies (LBGs) at $z=6$, named ‘Dahlia’ and ‘Althæa’, respectively. The two simulations differ for the H₂ formation implementation, equilibrium versus non-equilibrium. We show how chemistry has a strong impact on the observed properties of early galaxies.

The paper is organized as follows. In Section 2, we describe the two simulations highlighting common features (Section 2.1), separately discussing the different chemical models used for Dahlia (Section 2.2) and Althæa (Section 2.3). Results are presented as follow. First, we perform a benchmark of the chemical models (Section 2.4), and compare the star formation (SF) and feedback history of the two galaxies (Section 3.1). Next, we characterize their differences in terms of morphology (Section 3.2), thermodynamical state of the ISM (Section 3.3), and predicted [C II] and H₂ (Section 4) emission line properties. Our conclusions are summarized in Section 5.

## 2 NUMERICAL SIMULATIONS

To assess the impact of H₂ chemistry on the internal structure of high-z galaxies, we compare two zoom-in simulations adopting different chemical models. Both simulations follow the evolution of a prototypical $z=6$ LBG galaxy hosted by a $M_h \simeq 10^{11} M_\odot$ dark matter (DM) halo (virial radius $r_{vir} \simeq 15$ kpc).

The first simulation has been presented in Pallottini et al. (2017, hereafter P17). The targeted galaxy (which includes also about 10 satellites) is called ‘Dahlia’ (see also Gallerani et al. 2016, for analysis of its infall/outflow structure). In such previous work, we showed that Dahlia’s specific SFR (sSFR) is in agreement with both analytical calculations (Behroozi, Wechsler & Conroy 2013), and with $z=7$ observations (González et al. 2010, see also Section 3.1).

In the second, new simulation we follow the evolution of ‘Althæa’, by using improved thermo-chemistry, but keeping everything else (initial conditions, resolution, SF and feedback prescriptions) unchanged with respect to the Dahlia simulation. We describe the implementation of these common processes in the following section. Next, we describe separately the chemical model used for Dahlia (Section 2.2) and Althæa (Section 2.3).

### 2.1 Common physical models

Both simulations are performed with a customized version of the Adaptive Mesh Refinement (AMR) code RAMSES (Teyssier 2002). Starting from cosmological IC<fnref n="1" /> generated with MUSIC (Hahn & Abel 2011), we zoom-in on the $z \simeq 6$ DM halo hosting the targeted galaxy. The total simulation volume is $(20 \text{ Mpc h}^{-1})^3$ that is evolved with a base grid with eight levels (gas mass $6 \times 10^6 M_\odot$); the zoom-in region has a volume of $(2.1 \text{ Mpc h}^{-1})^3$ and is resolved with three additional levels of refinement, thus yielding a gas mass resolution of $m_b = 1.2 \times 10^4 M_\odot$. In such region, we allow for six additional levels of refinement which allow us to follow the evolution of the gas down to scales of $l_{cell} \simeq 30$ pc at $z=6$, i.e. the refined cells have mass and size typical of Galactic molecular clouds (MCs; e.g. Federrath & Klessen 2013). The refinement is performed with a Lagrangian mass threshold-based criterion., i.e. a cell is refined

<footnote n="1">¹ We assume cosmological parameters compatible with Planck results: ΛCDM model with total matter, vacuum and baryonic densities in units of the critical density $\Omega_\Lambda = 0.692$, $\Omega_m = 0.308$, $\Omega_b = 0.0481$, Hubble constant $H_0 = 100 \text{ h km s}^{-1} \text{ Mpc}^{-1}$ with $h = 0.678$, spectral index $n = 0.967$, $\sigma_8 = 0.826$ (Planck Collaboration XVI 2014).</footnote>
if its total (DM+baryonic) mass exceed the mass resolution by a factor 8.

Metallicity (Z) is followed as the sum of heavy elements, assumed to have solar abundance ratios (Asplund et al. 2009). We impose an initial metallicity floor $Z_{floor} = 10^{-3} Z_{\odot}$ since at $Z \gtrsim 40$ our resolution is still insufficient to catch the metal enrichment by the first stars (e.g. O'Shea et al. 2015). Such floor is compatible with the metallicity found at high-$z$ in cosmological simulations for diffuse enriched gas (Davé, Finlator & Oppenheimer 2011; Pallottini et al. 2014a; Maio & Tescari 2015); it only marginally affects the gas cooling time.

Dust evolution is not explicitly tracked during simulations. However, we make the simple assumption that the dust-to-gas mass ratio scales with metallicity, i.e. $D = D_{\odot}(Z/Z_{\odot})$, where $D_{\odot}/Z_{\odot} = 0.3$ for the Milky Way (MW; e.g. Hirashita & Ferrara 2002; Asano et al. 2013).

### 2.1.1 Star formation

Stars form according to a linearly H₂-dependent Schmidt-Kennicutt (SK) relation (Schmidt 1959; Kennicutt 1998) i.e.

$$ \dot{\rho}_{*} = \zeta_{sf} f_{H2} \frac{\rho}{l_{ff}}, \qquad (7) $$

where $\dot{\rho}_{*}$ is the local SF rate density, $\zeta_{sf}$ is the SF efficiency, $f_{H2}$ is the H₂ mass fraction and $\rho = \mu m_p n$ is density of the gas of mean molecular weight $\mu$. Equation (1) is solved stochastically, by drawing the mass of the new star particles from a Poisson distribution (Rasera & Teyssier 2006; Dubois & Teyssier 2008; Pallottini et al. 2014a). In detail, in an SF event we create a star particle with mass $Nm_b$, with N an integer drawn from

$$ P(N) = \frac{\langle N \rangle}{N!} \exp(-\langle N \rangle), \qquad (2) $$

where the mean of the Poisson distribution is

$$ \langle N \rangle = \frac{f_{H2} \rho l_{cell}^3 \zeta_{sf} \delta t}{m_b t_{ff}}, \qquad (3) $$

with $\delta t$ the simulation time-step. For numerical stability, no more than half of the cell mass is allowed to turn into stars. Since we prevent formation of star particle with mass less than $m_b$, cells with density less than ~15 cm⁻<fnref n="3" /> (for $l_{cell} \simeq 30$ pc) are not allowed to form stars.

We set $\zeta_{sf} = 0.1$, in accordance with the average values inferred from MC observations (Murray 2011, see also Agertz et al. 2013); $f_{H2}$ depends on the adopted thermo-chemical model, as described later in Sections 2.2 and 2.3.

### 2.1.2 Feedback

Similarly to Kim et al. (2014), we account for stellar energy inputs and chemical yields that depend both on time and stellar populations by using STARBURST99 (Leitherer et al. 1999). Stellar tracks are taken from the padova (Bertelli et al. 1994) library with stellar metallicities in the range $0.02 \le Z_{*}/Z_{\odot} \le 1$, and we assume a Kroupa (2001) initial mass function. Stellar feedback includes SNs, winds from massive stars and radiation pressure (Agertz et al. 2013). We model the thermal and turbulent energy content of the gas according to the prescriptions by Agertz & Kravtsov (2015). The turbulent (or non-thermal) energy is dissipated as $\dot{e}_{nth} = -e_{nth}/t_{diss}$ (Teyssier et al. 2013, see equation 2), where, following Mac Low (1999), the dissipation time-scale can be written as

$$ t_{diss} = 9.785 \left( \frac{l_{cell}}{100 \text{ pc}} \right) \left( \frac{\sigma_{turb}}{10 \text{ km s}^{-1}} \right)^{-1} \text{Myr}, \qquad (4) $$

where $\sigma_{turb}$ is the turbulent velocity dispersion. Adopting the SN blastwave models and OB/AGB stellar winds from Ostriker & McKee (1988) and Weaver et al. (1977), respectively, we account for the dissipation of energy in MCs as detailed in section 2.4 and appendix A of P17.

## 2.2 Dahlia: equilibrium thermo-chemistry

In the Dahlia simulation, we compute $f_{H2}$ by adopting the KTM09 analytical prescription (Krumholz, McKee & Tumlinson 2008, 2009; McKee & Krumholz 2010). In KTM09, the H₂ abundance is derived by modelling the radiative transfer on an idealized MC and by assuming equilibrium between H₂ formation on dust grains and dissociation rates. For each gas cell, $f_{H2}$ can then be written as a function of n, Z and hydrogen column density ($N_H$). By further assuming pressure equilibrium between CNM and WNM (Krumholz et al. 2009), $f_{H2}$ turns out to be independent of the intensity of the ISRF, and can be written as

$$ f_{H2} = [1 - 0.75 s/(1 + 0.25 s)] \Theta(2 - s), \qquad (5a) $$

with

$$ s = \ln(1 + 0.6 \chi + 0.01 \chi^2) / (0.6 \tau_{UV}) \qquad (5b) $$

$$ \chi = 0.75 [1 + 3.1 (Z/Z_{\odot})^{0.365}], \qquad (5c) $$

and where $\Theta$ is the Heaviside function; $\tau_{UV}$ is the dust UV optical depth and it can be calculated by linearly rescaling the MW value,

$$ \tau_{UV} = \left( \frac{N_H}{1.6 \times 10^{21} \text{ cm}^{-2}} \right) \left( \frac{D}{D_{\odot}} \right). \qquad (6) $$

In Dahlia cooling/heating rates are computed using GRACKLE 2.1² (Bryan et al. 2014). We use an H and He primordial network, and tabulated metal cooling/photoheating rates from CLOUDY (Ferland et al. 2013). Inverse Compton cooling is also present, and we consider heating from a redshift-dependent ionizing UV background (UVB; Haardt & Madau 2012). Since H₂ is not explicitly included in the network, we do not include the corresponding cooling contribution.

## 2.3 Althæa: non-equilibrium thermo-chemistry

In Althæa, we implement a non-equilibrium chemical network by using KROME<fnref n="3" /> (Grassi et al. 2014). Given a set of species and their reactions, KROME can generate the code needed to solve the system of coupled ordinary differential equations that describe the gas thermochemical evolution.

### 2.3.1 Chemical network

Similarly to Bovino et al. (2016, hereafter B16), our network includes H, H⁺, H⁻, He, He⁺, He⁺⁺, H₂, H₂⁺ and electrons. Metal species are not followed individually in the network, as for instance done in Model IV from B16; therefore, we use an equilibrium metal line cooling calculated via CLOUDY tables.<fnref n="4" /> The adopted

<footnote n="2">² https://grackle.readthedocs.org/</footnote>

<footnote n="3">³ https://bitbucket.org/tgrassi/krome</footnote>

<footnote n="4">⁴ As a caveat, we point out that there is a formal inconsistency in the modelling. Metal line cooling tables are usually calculated with CLOUDY by assuming a Haardt & Madau (2012) UVB, while the ISRF SED we adopt is MW-like. To remove such inconsistency one should explicitly track metal</footnote>
network contains a total of 37 reactions, including photochemistry (Section 2.3.2), dust processes (Section 2.3.3) and cosmic rays (CR; Section 2.3.4). The reactions, their rates, and corresponding references are listed in appendix B of B16: specifically we use reactions from 1 to 31 (table B.1 in B16), 53, 54 and from 58 to 61 (table B.2 in B16).

### 2.3.2 Photochemistry

Photochemistry cross-sections are taken from Verner & Ferland (1996) and by using the SWRI<fnref n="5" /> and Leiden<fnref n="6" /> data bases. In the present simulation, the ISRF is not evolved self-consistently and it is approximated as follows. For the spectral energy density (SED), we assume an MW-like spectrum (Black 1987; Draine 1978), and we specify the SED using 10 energy bins from 0.75 to 14.16 eV. Beyond 13.6 eV the flux drops to zero, i.e. we do not include ionizing radiation.

We consider a spatially uniform ISRF whose intensity is rescaled with the SFR such that $G = G_0(SFR/M_\odot \text{yr}^{-1})$, where $G_0 = 1.6 \times 10^{-3} \text{ erg cm}^{-2} \text{ s}^{-1}$ is the far-UV (FUV) flux in the Habing band (6–13.6 eV) normalized to the average MW value (Habing 1968). Because of their sub-kpc sizes (Shibuya, Ouchi & Harikane 2015; Fujimoto et al. 2017), high $G_0$ values are expected in typical LBGs at $z \simeq 6$, as inferred also by Carniani et al. (2017). A similar situation is seen in some local dwarf galaxies (Cormier et al. 2015) that are generally considered as local counterparts of high-z galaxies. It is worth noting that the spatial variation of $G$ is very small in the MW, with an rms value $\simeq 3 G_0$ (Habing 1968; Wolfire et al. 2003). None the less, spatial fluctuations of the ISRF, if present, might play some role in the evolution of high-z galaxies (e.g. Katz et al. 2016). We will analyse this effect in future work.

On top of the ISRF, we consider the cosmic microwave background (CMB), which effectively sets a temperature floor for the gas. Additionally, we neglect the cosmic UVB, since the typical ISM densities are sufficiently large to ensure an efficient self-shielding (e.g. Gnedin 2010). For example, Rahmati et al. (2013) have shown that at $z \simeq 5$ the hydrogen ionization due to the UVB is negligible for $n \gtrsim 10^{-2} \text{ cm}^{-3}$, the typical density of diffuse ISM. The self-shielding of H$_2$ to photodissociation is accounted by using the Richings, Schaye & Oppenheimer (2014) prescription,<fnref n="7" /> thus in each gas cell the shielding can be expressed as an analytical function of its H$_2$ column density, temperature and turbulence (cf. with Wolcott-Green, Haiman & Bryan 2011).

### 2.3.3 Dust processes

As for Dahlia, the dust mass is proportional to the metal mass. Here, we also specify the dust size distribution to be the one appropriate for high-z galaxies, the Small Magellanic Cloud one, following Weingartner & Draine (2001). Dust grains can affect the chemistry

species, adopt a non-equilibrium metal line cooling and include radiative transfer. As noted in B16 (see their fig. 16), using non-equilibrium metal line cooling can typically change the cooling function by a factor $\lesssim 2$. This will be addressed in future work.

through cooling<fnref n="8" /> (Hollenbach & McKee 1979), photoelectric heating (Bakes & Tielens 1994), and by mediating the formation of molecules (Cazaux & Spaans 2009). In particular, the formation rate of H$_2$ on dust grains is approximated following Jura (1975)

$$R_{\text{H2-dust}} = 3 \times 10^{-17} n_{\text{H}} (D/D_{\odot}) \text{cm}^{-3} \text{s}^{-1}, \quad (7)$$

where $n_H$ is the hydrogen density. Note that for $D \gtrsim 10^{-2} D_\odot$, this dust channel is dominant with respect to gas-phase formation (e.g. reactions 6–7 and 9–10, table B.1 in B16).

### 2.3.4 Cosmic rays

CR ionization can become important in regions shielded from radiation, like MC interiors. We assume a CR hydrogen ionization rate $\propto$ SFR (Valle et al. 2002) and normalized to the MW value (Webber 1998):

$$\zeta_{\text{cr}} = 3 \times 10^{-17} (\text{SFR}/\text{M}_{\odot} \text{ yr}^{-1}) \text{ s}^{-1}. \quad (8)$$

The rate $\zeta_{\text{cr}}$ includes the flux of CR and secondary electrons (Richings et al. 2014). In the network, CR ionizations are proportional to $\zeta_{\text{cr}}$ and to coupling constants that depend on the specific ions; such couplings are taken from the KIDA data base (Wakelam et al. 2012). Additionally we account for Coulomb heating, by assuming that every CR ionization releases an energy<fnref n="9" /> of 20 eV.

### 2.3.5 Initial abundances of the species

Finally, following Galli & Palla (1998), we calculate IC for the various species by accounting for the primordial chemistry<fnref n="10" /> at $z \gtrsim 100$, for a density and temperature evolution corresponding to gas at the mean cosmic density.

## 2.4 Benchmark of H$_2$ formation models

As a benchmark for our simulations, we compare the formation of H$_2$ in different physical environments. For the Dahlia KTM09 model, we compute $f_{\text{H2}}$ from equation (5) as a function of $n$ and $Z$. We choose an expression for $N_H = n l_{\text{cell}} \mu \propto n^{2/3}$ resulting from the mass threshold-based AMR refinement criterion for which $l_{\text{cell}} \propto n^{-1/3}$. We restate that the equilibrium KTM09 model is independent of $G$ and the gas temperature $T$.

For the Althæa B16 model, we use KROME to perform single-zone tests varying $n$, $Z$ and $G$. In this case, we assume an initial temperature<fnref n="11" /> $T = 5 \times 10^3 \text{K}$, and we let the gas patch evolve at constant density until thermo-chemical equilibrium is reached. This typically takes 100 Myr.

The comparison between the two models is shown in Fig. 1 as a function of $n$ for different metallicities. For $G > 0$ and $Z < Z_\odot$, H$_2$ formation is hindered in B16 with respect to KTM09, i.e. higher $n$ are needed to reach similar $f_{\text{H2}}$ fractions. For $G = G_0$ and $Z = Z_\odot$, the two models are roughly in agreement: this is expected since KTM09 is calibrated on the MW environment. Finally, for $G > 0$ and at

<footnote n="5">⁵ http://phidrates.space.swri.edu</footnote>

<footnote n="6">⁶ http://home.strw.leidenuniv.nl/ewine/photo/</footnote>

<footnote n="7">⁷ The self-shielding formulation by Richings et al. (2014) does not account for a directional dependence as done in more computationally costly models (Hartwig et al. 2015).</footnote>

<footnote n="8">⁸ Dust cooling is not included in the current model, as it gives only a minor contribution for $n < 10^4 \text{cm}^{-3}$, i.e. see fig. 3 in B16.</footnote>

<footnote n="9">⁹ For a more accurate treatment of Coulomb heating refer to Glassgold, Galli & Padovani (2012).</footnote>

<footnote n="10">¹⁰ For a possible implementation of the Galli & Palla (1998) chemical network, see the 'early Universe' test contained in KROME.</footnote>

<footnote n="11">¹¹ The initial temperature corresponds to the virial temperature of the first star-forming haloes present in the zoomed region. The results depend very weakly on this assumption.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Benchmark of the formation of H₂ for the model used in Dahlia (KTM09, Section 2.2) and in Althæa (B16, Section 2.3). In each panel we plot the H₂ mass fraction $f_{H2}$ as a function of density ($n$), with different panels showing the results for different metallicities ($Z$). In each panel the dashed grey line indicates the KTM09 model, while the B16 models are plotted with solid lines, with different colours indicating a different impinging ISRF flux ($G$). In the upper axis, we indicate the free-fall times ($t_{ff}$) corresponding to $n$.</figcaption></figure>

$Z = 10^{-3} Z_{\odot}$ (the metallicity floor in our simulation set) the H₂ formation in the B16 model is strongly suppressed, e.g. $f_{H2} \simeq 10^{-3}$ for $n \simeq 10^3 \text{cm}^{-3}$. Note that these fractions are comparable to the ones expected for H₂ formation in a pristine environment where H₂ formation proceeds via gas-phase reactions.

As noted in P17, Dahlia's SF model (equations 1 and 5) is roughly equivalent to a density threshold criterion with metallicity-dependent critical density $n_c \simeq 26.45 (Z/Z_{\odot})^{-0.87} \text{cm}^{-3}$. Physically, this corresponds to the density at which $f_{H2} \ge 0.5$ (see also Agertz et al. 2013). Thus, Fig. 1 quantifies the density threshold required to spawn stars in the simulation.

Dahlia forms stars in gas with $n \simeq 30 \text{ cm}^{-3}$ and $Z \simeq 0.5 Z_{\odot}$ at a rate of about $10^2 \text{ M}_{\odot} \text{yr}^{-1}$ at $z = 6$. If Althæa has a similar SFR history (this is checked a posteriori, see Fig. 2), the resulting metallicity and ISRF intensity ($G \simeq 10^2 G_0$) should also be similar. Then, by inspecting Fig. 1 (middle-left panel) one can conclude that Althæa forms stars in much denser environments where $n > n_c \simeq 263 (Z/Z_{\odot})^{-1.19} \text{cm}^{-3}$ for $G = 10^2 G_0$. As noted in Hopkins, Narayanan & Murray (2013), although variations in the density threshold lead to similar total SFR, they might severely affect the galaxy morphology. We will return to this point in the next section.

We remind that in both simulations we use the cell radius to calculate column density, which is used e.g. to calculate the gas self-shielding. This is done to mainly ensure a fair comparison between the two simulations. In other simulations, e.g. MCs illuminated by an external radiation field, the prescriptions adopted accounts for the contribution of column density from nearby cells, i.e. by using Jeans or Sobolev-like length (see e.g. Hartwig et al. 2015 for a comparison between different prescriptions). However, in our simulation we expect stars to be very close or embedded in potential star-forming regions. Using the contribution to the column density from the surrounding gas would then overestimate the self-shielding effect. Such modelling uncertainty would be solved by including radiative transfer in the simulation. However, we note that at $z = 6$ the radius of our cells as a function of density can be approximated as $r_{\text{cell}} = 154.1 (n/\text{cm}^{-3})^{-1/3}$ pc, while the Jeans length is

<figure><img src="image_3.png" /><figcaption>**Figure 2.** SFR as a function of galaxy age ($t_*$) for Dahlia (black line and hatched region) and Althæa (orange line and transparent region). Also shown (grey dashed line) is an analytical approximation (within a factor 2 for both galaxies) to the average SFR trend. The redshift ($z$) corresponding to $t_*$ is plotted on the upper axis, and note that $t_* = 0$ corresponds to the first stellar formation event in Dahlia, and the plotted SFRs are averaged over 4 Myr.</figcaption></figure>

$l_j = 15.6 (n/\text{cm}^{-3})^{-1/2} (T/\text{K})^{1/2}$ pc. Thus, for typical values found for the molecular gas in Althæa ($n \simeq 300 \text{ cm}^{-3}$ and $T \simeq 100$ K, see later Fig. 8), the two prescriptions gives similar results, i.e. $r_{\text{cell}} \sim 20$ pc and $l_j \sim 10$ pc.

## 3 RESULTS

We now turn to a detailed analysis of the two zoomed galaxies, Dahlia and Althæa.<fnref n="12" /> We start by studying the SF and the build-up of the stellar mass from $z \simeq 15$ to $z = 6$ (Section 3.1). We then specialize at $z = 6$ to inspect the galaxy morphology (Section 3.2), the ISM multiphase structure (Section 3.3) and the predicted observable properties (Section 4). An overall summary of the properties of the two galaxies is given in Table 1.

### 3.1 Star formation history

In Fig. 2 we plot the SFR history as a function of ‘galaxy age’ ($t_*$) for Dahlia and Althæa; $t_* = 0$ marks the first SF event in Dahlia.<fnref n="13" /> For both Dahlia and Althæa the SFR has an increasing trend that can be approximated with good accuracy (within a factor of 2) as SFR = 1.5 log($t_*/(30 \text{ Myr})$). However, on average the SFR in Dahlia is larger by a factor $\simeq 1.5 \pm 0.6$ when averaged over the entire SFR history ($\simeq 700 \text{ Myr}$). Thus, in spite of very different chemical prescriptions, the SFR in the two galaxies shows very little variation. Stated differently, the higher critical density for SF arising from non-equilibrium chemistry does not alter significantly the rate at which stars form, as already noticed in Section 2.4. This

<footnote n="12">¹² We refrain from the analysis of the satellite population of the two galaxies due to the oversimplifying assumption of a spatially uniform ISRF artificially suppressing SF in environments with metallicity close to the floor value $Z_{\text{floor}} = 10^{-3} Z_{\odot}$.</footnote>

<footnote n="13">¹³ Note that even with the same modelling and IC, differences in the SFR may arise as a result of stochasticity in the SF prescription (equation 1). Such differences vanish once the SFR is averaged on time-scales longer than the typical free-fall time of the star-forming gas.</footnote>
<table><caption>**Table 1.** Physical properties of Dahlia and Althæa at z = 6. The values refer to gas and stars within 2.5 kpc from the galaxy centre (similar to the field of view in Figs 5 and 11). The effective radius, $r_d$, and gas scaleheight, $H$, are calculated from the principal component analysis of the density field. Values for $n$, $n_{H2}$, $Z$, $\Sigma$ and $\dot{\Sigma}_*$ represent mass-weighted averages.</caption><thead><tr><th>Property</th><th>Symbol</th><th>Dahlia</th><th>Althæa</th><th>(units)</th></tr></thead><tbody><tr><td>Star formation rate</td><td>SFR</td><td>156.19</td><td>136.50</td><td>M<sub>☉</sub> yr<sup>-1</sup></td></tr><tr><td>Specific SFR</td><td>sSFR</td><td>4.45</td><td>5.23</td><td>Gyr<sup>-1</sup></td></tr><tr><td>Stellar mass</td><td>M<sub>*</sub></td><td>3.51</td><td>2.61</td><td>10<sup>10</sup> M<sub>☉</sub></td></tr><tr><td>Metal mass in stars</td><td>M<sub>*</sub><sup>Z</sup></td><td>8.20</td><td>5.87</td><td>10<sup>8</sup> M<sub>☉</sub></td></tr><tr><td>Gas mass</td><td>M<sub>g</sub></td><td>1.23</td><td>2.72</td><td>10<sup>9</sup> M<sub>☉</sub></td></tr><tr><td>H<sub>2</sub> mass</td><td>M<sub>H2</sub></td><td>17.01</td><td>4.76</td><td>10<sup>7</sup> M<sub>☉</sub></td></tr><tr><td>Metal mass</td><td>M<sub>Z</sub></td><td>1.41</td><td>2.48</td><td>10<sup>7</sup> M<sub>☉</sub></td></tr><tr><td>Disc radius</td><td>r<sub>d</sub></td><td>610</td><td>504</td><td>pc</td></tr><tr><td>Disc scaleheight</td><td>H</td><td>224</td><td>191</td><td>pc</td></tr><tr><td>Gas density</td><td>⟨n⟩</td><td>23.89</td><td>164.41</td><td>cm<sup>-3</sup></td></tr><tr><td>H<sub>2</sub> density</td><td>⟨n<sub>H2</sub>⟩</td><td>6.62</td><td>4.95</td><td>cm<sup>-3</sup></td></tr><tr><td>Metallicity</td><td>⟨Z⟩</td><td>0.57</td><td>0.46</td><td>Z<sub>☉</sub></td></tr><tr><td>Gas surface density</td><td>⟨Σ⟩</td><td>37.89</td><td>222.02</td><td>M<sub>☉</sub> pc<sup>-2</sup></td></tr><tr><td>Star formation surface density</td><td>⟨Σ<sub>*</sub>⟩</td><td>0.40</td><td>0.83</td><td>M<sub>☉</sub> pc<sup>2</sup> Myr<sup>-1</sup></td></tr><tr><td>Luminosity [C II]</td><td>L<sub>C II</sub></td><td>3.39</td><td>21.08</td><td>10<sup>7</sup> L<sub>☉</sub></td></tr><tr><td>Luminosity H<sub>2</sub></td><td>L<sub>H2</sub></td><td>2.31</td><td>33.24</td><td>10<sup>5</sup> L<sub>☉</sub></td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** SFR versus stellar mass ($M_*$) for Althæa (circles) and Dahlia (squares), with symbols coloured accordingly to the age $t_*$, with crosses we overplot SFR and $M_*$ inferred from 27 galaxies observed at $z \simeq 6$ by J16. Following J16 analysis<fnref n="14" />, galaxies identified as young and old are plotted in blue and red, respectively. To guide the eye, the linear correlation between the data sets is also shown with a dashed line. See the text for more details.</figcaption></figure>

also entails a comparable metallicity, and we note that in both galaxy most of the metal mass is locked in stars (see Table 1), as they are typically formed from the most enriched regions.

It is interesting to check the evolutionary paths of Dahlia and Althæa (Fig. 3) in the standard SFR versus stellar mass ($M_*$) diagram, and compare them with data<fnref n="14" /> inferred from $z \simeq 6$ observations of 27 Lyman alpha emitters (LAEs) and LBGs (Jiang et al. 2016, hereafter J16). By using multiband data, precise redshift determinations, and an estimate of nebular emission from Ly$\alpha$, J16 were able to distinguish between a young ($t_* \lesssim 30$ Myr) and an old ($t_* \gtrsim 100$ Myr) subsample. Each subsample exhibits a linear correlation

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Ratio of mechanical ($\dot{E}_{sn+w}$) and radiative ($\dot{E}_{rad}$) energy deposition rates by stars as a function of galaxy age ($t_*$) for Dahlia (black line/hatched area) and Althæa (orange/transparent). Dashed lines indicate the $\simeq 700$ Myr time-averaged mean of the ratios for each galaxy. To guide the eye, we plot the unity value (dotted grey line). Similar to Fig. 2 the ratios are averaged over 4 Myr. The upper horizontal axis indicates redshift.</figcaption></figure>

in log SFR – log $M_*$, albeit with a different normalization: the young (old) subsample has an sSFR = SFR/$M_*$ = 39.7 Gyr$^{-1}$(4.1 Gyr$^{-1}$).

The SFR versus stellar mass of our simulated galaxies for $M_* \lesssim 10^{8.5} M_☉$ ($t_* \lesssim 100$ Myr) is fairly consistent with the young subsample relation (keeping in mind stochasticity effects at low stellar masses). At later evolutionary stages ($t_* \gtrsim 300$ Myr or $M_* \gtrsim 10^{9.5} M_☉$), Dahlia and Althæa nicely shift to the lower sSFR values characterizing the old J16 subsample data. This shift must be understood as a result of increasing stellar feedback: as galaxies grow, the larger energy input from the accumulated stellar populations hinders subsequent SFR events. Note that at late times ($t_* \gtrsim 300$ Myr), when $M_* = 5 \times 10^9 M_☉$, the sSFRs of Dahlia and Althæa are in agreement with analytical results by Behroozi et al. (2013), and with $z = 7$ observations by González et al. (2010).

As feedback clearly plays a major role in the overall evolution of early galaxies, we turn to a more in-depth analysis of its energetics. This can be quantified in terms of the stellar energy deposition rates in mechanical (SN explosions + OB/AGB winds,<fnref n="15" /> $\dot{E}_{sn+w}$) and radiative ($\dot{E}_{rad}$) forms. These are shown as a function of time in Fig. 4. The $\dot{E}_{sn+w}/\dot{E}_{rad}$ ratio shows short-term ($\simeq 20$ Myr) fluctuations corresponding to coherent burst of SF/SN activity.

Barring this time modulation, on average the mechanical/radiative energy ratio increases up to $\simeq 250$ Myr, when it suddenly drops and reaches an equilibrium value. This implies that radiation pressure dominates the energy input; consequently it represents the major factor in quenching SF. While this is true throughout the evolution, it becomes even more evident after $\simeq 250$ Myr, when the first stellar populations with $Z_* \gtrsim 10^{-1} Z_☉$ enter the AGB phase. At that time, winds from AGBs enrich their surroundings with metals and dust. As dust produced by AGBs remains more confined than SN dust around the production region, it provides a higher opacity, thus boosting radiation pressure via a more efficient dust–gas coupling (see also P17).

<footnote n="14">¹⁴ SFR and $M_*$ have been derived by assuming an exponentially increasing SFR, consistent with the history of both our simulated galaxies (Fig. 2).</footnote>

<footnote n="15">¹⁵ On average OB/AGB winds account only for $\lesssim 10$ percent of the SN power.</footnote>
For Dahlia the radiative energy input rate is about 20 times larger than the mechanical one, while for Althæa such ratio is on average eight times higher, although larger fluctuations are present. The latter are caused by the occurrence of more frequent and powerful bursts of SN events in Althæa. Why does this happen?

The answer has to do with the different gas morphology. As already noted discussing Fig. 1, the higher critical density for SF imposed by non-equilibrium chemistry has a number of consequences: (a) each formation event can produce a star cluster with an higher mass; (b) SF is more likely hosted in isolated high-density clumps (see later, particularly Fig. 6); (c) in a clumpier disc, SN explosions can easily break into more diffuse regions. The combination of (a) and (b) increases the probability of spatially coherent explosions having a stronger impact on the surrounding gas; due to (c), the blastwaves suffer highly reduced radiative losses (Gatto et al. 2015), and affect larger volumes. Similar effects have been also noted in the context of single giant MCs (~$10^6$ M$_{\odot}$), where unless the SNs explode coherently, their energy is quickly radiated away because of the very high gas densities (Rey-Raposo et al. 2017). For the reminder of the work we focus on $z = 6$, when the galaxies have an age of $t_* \simeq 700$ Myr.

## 3.2 Galaxy morphology

Dahlia and Althæa sit at the centre of a cosmic web knot and accrete mass from the intergalactic medium mainly via three filaments of length $\simeq 100$ kpc. In both simulations, the large-scale structure is similar, and we refer the reader to Section 3.1 of P17 for its analysis. Differences between the simulation are expected to arise on the ISM scale, whose structure is visible on $\simeq 7$ kpc scales. In Fig. 5, we show the gas density, temperature and H$_2$ density ($n_{H2} = f_{H2} n \mu$) fields for Dahlia and Althæa. The map$^{16}$ centres coincide with Dahlia’s stellar centre of mass.

### 3.2.1 Overview

Qualitatively, both galaxies show a clearly defined, albeit somewhat perturbed, spiral disc of radius $\simeq 0.5$ kpc, embedded in a lower density ($n \simeq 0.1$ cm$^{-3}$) medium. However, the mean disc gas density for Dahlia is $\langle n \rangle = 24$ cm$^{-3}$, while for Althæa $\langle n \rangle = 164$ cm$^{-3}$ (see Table 1). The temperature structure shows fewer differences, i.e. the inner disc is slightly hotter for Dahlia ($T \simeq 300$ K) than for Althæa ($T \simeq 100$ K), which features instead slightly more abundant and extended pockets of shock-heated gas ($T \gtrsim 10^6$). Such high-T regions are produced by both accretion shocks and SN explosions. In both cases the typical H$_2$ density is the same, i.e. $\langle n_{H2} \rangle = 5$ cm$^{-3}$, however, with respect to Dahlia, Althæa shows a slightly smaller disc that also seems more clumpy.

To summarize, the galaxies differ by an order of magnitude in atomic density, but have the same molecular density. In spite of this difference, the SFR are roughly similar. This can be explained as follows. To first order, in our model SFR $\propto n_{H2} n^{1/2} V$, where $V = 2\pi r_d^2 H$ is the galaxy volume (Table 1). It follows that the larger density is largely compensated by the smaller Althæa volume.

### 3.2.2 In-depth analysis

Fig. 5 visually illustrates the morphological differences between the two galaxies. The gas in Althæa appears clumpier than in Dahlia. To quantify this statement we start by introducing the H$_2$ clumping factor on the smoothing scale $r$, which is defined as$^{17}$

$$C(r) = \langle n_{H2}^2 \rangle / \langle n_{H2} \rangle^2. \qquad (9)$$

For Dahlia, $C(r)$ decreases from $10^3$ to 10 going from 30 pc to 1 kpc, while for Althæa $C(r)$ is $\gtrsim 2$ times larger on all scales.

A more in-depth analysis can be performed using the Minkowsky functionals (Schmalzing & Gorski 1998; Gleser et al. 2006; Yoshiura et al. 2017, Appendix A) that can give a complete description of the molecular gas morphological structure. For a 3D field, four independent Minkowsky functionals can be defined. Each of the functionals, $V_i(n_{H2})$ ($i = 0, \dots, 3$) characterizes a different morphological property of the excursion set with H$_2$ density $> n_{H2}$: $V_0$ gives the volume filling factor, $V_1$ measures the total enclosed surface, $V_2$ is the mean curvature, quantifying the sphericity/concavity of the set and $V_3$ estimates the Euler characteristic (i.e. multiple isolated components versus a single connected one. Appendix A gives more rigorous definitions with an illustrative application (Fig. A1).

In Fig. 6, we plot the Minkowsky functionals ($V_0$, $V_1$, $V_2$, $V_3$) calculated for the H$_2$ density field for Dahlia and Althæa. The $V_0$ functional analysis shows that Althæa is more compact, i.e. for each $n_{H2}$ value Dahlia’s excursion set volume is larger and it plummets rapidly at large densities. On the other hand, the set surface of Althæa is larger by about a factor of 5, implying that this galaxy is fragmented into multiple, disconnected components. This is confirmed also by Althæa’s larger ($10 \times$) Euler characteristic measure, $V_3$, an indication of the prevalence of isolated structures. This feature becomes more evident towards larger densities, as expected if H$_2$ is concentrated in MCs.$^{18}$

Further, in Dahlia most of the molecular gas resides in connected ($V_3 \lesssim 0$) disc regions, with a concave shape ($V_2 < 0$). For Althæa, there is a transition: for $\log(n_{H2}/\text{cm}^{-3}) \lesssim 1$ the gas has a concave ($V_2 < 0$), disjointed ($V_3 > 0$), filamentary structure, while for $\log(n_{H2}/\text{cm}^{-3}) \gtrsim 1$ the galaxy is composed by spherical clumps ($V_2 > 0$).

## 3.3 ISM thermodynamics

The thermodynamical state of the ISM can be analysed by studying the probability distribution function (PDF) of the gas in the density–temperature plane, i.e. the equation of state (EOS). In Fig. 7, we plot the mass-weighted EOS for Dahlia and Althæa at $z = 6$. We include gas within 30 kpc, or $\simeq 2 r_{vir}$, from the galaxy centre.

From the EOS, we can see that in both galaxies 70 per cent of the gas in a photoionized state ($T \sim 10^4$ K), which in Dahlia is induced by the Haardt & Madau (2012) UVB, while in Althæa is mainly due to photoelectric heating on dust grains illuminated from the uniform ISRF of intensity $G$. Only $\simeq 10$ per cent of the gas is in a hot $10^6$ K component produced by accretion shocks and SN explosions. A relatively minor difference descends from Althæa’s more effective mechanical feedback, already noted when discussing

<footnote>$^{16}$ The maps of this work are rendered by using a customized version of PYMSES (Labadens et al. 2012), a visualization software that implements optimized techniques for the AMR of RAMSES.</footnote>

<footnote>$^{17}$ To calculate the clumping factor, first we construct the 3D unigrid cube of the H$_2$ mass field, then we smooth it with a Gaussian kernel of scale $r$ and finally we calculate the mass-weighted average and variance of the smoothed H$_2$ density field.</footnote>

<footnote>$^{18}$ $V_3 > 0$ values at $\log(n_{H2}/\text{cm}^{-3}) \simeq 1.2$ in Dahlia are mainly due to the presence of the three satellites/clumps outside the disc.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Face-on maps¹⁶ of Dahlia (left-hand panels) and Althæa (right) at age $t_* \simeq 700$ Myr ($z=6$). Shown are line-of-sight mass-weighted average of the gas density (upper panels), temperature (middle) and H₂ density (lower) fields with amplitude given by the colourbar. The maps are 6.31 kpc on a side.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Morphological comparison of the molecular gas at $z = 6$. In the four panels, we plot the Minkowsky functionals ($V_0, V_1, V_2, V_3$) of the H$_2$ density field ($n_{H2}/\text{cm}^{-3}$). Functionals are plotted with black line and hatched regions for Dahlia, with orange line and transparent region for Althæa. Note that Minkowsky functionals are indicated in comoving units. For detail on the calculation of the Minkowsky functional, see Appendix A (in particular see Fig. A1).</figcaption></figure>

Fig. 4: small pockets of freshly produced very hot ($\ge 10^6$ K) and diffuse ($0.1 \text{ cm}^{-3}$) gas are twice more abundant in Althæa, as it can be appreciated from a visual inspection of the temperature maps in Fig. 5.

Fig. 7 (in particular compare the upper horizontal panels) shows that the density PDF is remarkably different in the two galaxies. In Dahlia the distribution peaks at $0.1 \text{ cm}^{-3}$; Althæa instead features a bimodal PDF with a second, similar amplitude peak at $n \sim 100 \text{ cm}^{-3}$. This entails the fact that the dense $\gtrsim 10 \text{ cm}^{-3}$ gas

is about two times more abundant in the latter system. In addition, the very dense gas ($n \gtrsim 300 \text{ cm}^{-3}$), only present in Althæa, can cool to temperatures of 30 K, not too far from the CMB one.

The high-density part of the PDF is worth some more insight as it describes the gas that ultimately regulates SF. This gas is largely in molecular form, and accounts (see Table 1) for 1.7 per cent (13.8 per cent) of the total gas mass in Dahlia (Althæa). Its H$_2$ density-weighted distribution in the $n-T$ plane is reported in Fig. 8. On average, the H$_2$ gas in Dahlia is 10 times less dense than in Althæa as a result of the new non-equilibrium prescription requiring higher gas densities to reach the same $f_{H2}$ fraction; at the same time the warm ($T \gtrsim 10^3$ K) H$_2$ fraction drops from 20 per cent (Dahlia) to an almost negligible value. Clearly, the warm component was a spurious result as (a) H$_2$ cooling was not included, and (b) $f_{H2}$ was considered to be independent of gas temperature (see equations 5). Note that in Althæa traces of warm H$_2$ are only found at large densities, in virtually metal-free gas in which H$_2$ production must proceed via much less efficient gas-phase reactions rather than on dust surfaces. This tiny fraction of molecular gas can survive only if densities large enough to provide a sufficient H$_2$ self-shielding against photodissociation are present.

Finally, the sharp EOS cut-off at $n \gtrsim 10^2 \text{ cm}^{-3}$ in Dahlia is caused by the density-threshold behaviour mimicked by the enforced chemical equilibrium: above $n_c \simeq 26.45 (Z/Z_\odot)^{-0.87} \text{ cm}^{-3}$ (Section 2.2) the gas is rapidly turned into stars. This spurious effect disappears in Althæa, implementing a full time-dependent chemical network.

## 4 OBSERVATIONAL PROPERTIES

As we already mentioned, the strongest impact of different chemistry implementations is on the gas properties, and consequently on ISM-related observables. In the following, we highlight the most important among these aspects.

<figure><img src="image_8.png" /><figcaption>**Figure 7.** EOS of the gas within 30 kpc for Dahlia (left-hand panel) and Althæa (right-hand panel) at $t_* \simeq 700$ Myr ($z=6$). EOS are shown as mass-weighted PDF in the density–temperature ($n-T$) plane, as specified by the colourbar. For both galaxies, the EOS projection on the $n(T)$ axis is additionally shown as a horizontal (vertical) inset. The 2D EOS are normalized such that the integral on the $n-T$ plane is unity; the projected EOS are normalized such that the sum of the bins is equal to 100 per cent.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** EOS of the molecular (H₂) gas for Dahlia (left-hand panel) and Althæa (right-hand panel) i.e. the H₂ mass-weighted PDF in the *n*-T plane. Notation is similar to Fig. 7, albeit a different region of *n*-T plane is shown.</figcaption></figure>

## 4.1 Schmidt-Kennicutt relation

We start by analysing the classical SK relation. This comparison should be interpreted as a consistency check of the balance between SF and feedback, since in the model we assume an SFR law that mimics an SK relation (equation 1).

The SK relation, in its most modern (Krumholz et al. 2012) formulation, links the SFR ($\dot{\Sigma}_{*}$) and total gas ($\Sigma$) surface density per unit free-fall time, $\dot{\Sigma}_{*} = \epsilon_{*}^{ff} \Sigma / t_{ff}$. The proportionality constant, often referred to as the efficiency per free-fall time following equation (1), is simply $\epsilon_{*}^{ff} = \zeta_{sf} f_{H2}$. Experimentally, Krumholz et al. (2012) find $\epsilon_{*}^{ff} = 0.015$ (see Krumholz 2015 for a complete review on the subject). This result is supported also by a larger set of observations including single MCs (Heiderman et al. 2010; Lada et al. 2010), local unresolved galaxies (Kennicutt 1998) and moderate redshift, unresolved galaxies (Bouché et al. 2007; Daddi et al. 2010a,b; Genzel et al. 2010; Tacconi et al. 2010). The SK relation is shown in Fig. 9, along with the location of Dahlia and Althæa at $z = 6$.

Dahlia appears to be overforming stars with respect to its gas mass, and therefore it is located about $3\sigma$ above the KS relation. As Althæa needs about 10 times higher density to sustain the same SFR, its location is closer to expectations from the SK. We have checked that the agreement is even better if we use only data relative to MC complexes (e.g. Heiderman et al. 2010; Murray 2011).

Dahlia's $\epsilon_{*}^{ff} = \zeta_{sf} f_{H2}$ is similar to the analogue values found by Semenov, Kravtsov & Gnedin (2016), who compute such efficiency using a turbulent eddy approach (Padoan, Haugbølle & Nordlund 2012), with no notion of molecular hydrogen fraction. The difference is that Dahlia misses the high-density gas. Althæa instead matches both the $\epsilon_{*}^{ff}$ and the amount of high-density gas found by Semenov et al. (2016). Also, its $\dot{\Sigma}_{*}-\Sigma$ relation is consistent with Torrey et al. (2017), who use an SF recipe involving self-gravitating gas with a local SK H₂ dependent relation.

From our simulations, it is also possible to perform a cell-by-cell analysis of the SK relation (Fig. 10). As expected, the results show the presence of a consistent spread in the local efficiency values

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Comparison of the observed and simulated SK relation expressed in terms of $\dot{\Sigma}_{*}-\Sigma / t_{ff}$. Observations are taken from single MCs (Heiderman et al. 2010; Lada, Lombardi & Alves 2010), local unresolved galaxies (Kennicutt 1998), and moderate redshift unresolved galaxies (Bouché et al. 2007; Daddi et al. 2010a,b; Genzel et al. 2010; Tacconi et al. 2010); the correlation (dispersion) for the observation found by Krumholz, Dekel & McKee (2012, see the text for details) is plotted with a black dashed line (grey shaded region). Dahlia and Althæa averaged value are plotted with black and orange stars, respectively (see Fig. 10 for the complete distribution in the simulated galaxies).</figcaption></figure>

which, however, has a different origin for Dahlia and Althæa. While in the former the variation is mostly due to a different enrichment level affecting H₂ abundance (equation 5), for Althæa the spread is larger because it results also from the individual evolutionary histories of the cells.

As noted by Rosdahl et al. (2017), for galaxy simulations with an SF model based on SK-like relation (equation 1), the resulting $\epsilon_{*} = \epsilon_{*}^{ff} / t_{ff}$ depends on how the feedback is implemented. However,
<figure><img src="image_11.png" /><figcaption>**Figure 10.** SK relation in Dahlia (left-hand panel) and Althæa (right-hand panel) at $t_* \simeq 700$ Myr ($z=6$). The relation is plotted using the H$_2$ mass weighted PDF of the instantaneous SFR surface density ($\dot{\Sigma}_*/M_{\odot} \text{ pc}^{-2} \text{ Myr}^{-1}$) versus the total gas surface density ($\Sigma/M_{\odot}/\text{pc}^2$). On both panels with dashed grey lines, we overplot the relation observed from Kennicutt & Evans (2012), i.e. $\dot{\Sigma}_* \propto \Sigma^{1.4}$, for several normalizations that written inline. Otherwise notation is similar to Figs 7 and 8.</figcaption></figure>

here we show that Althæa has a lower $\epsilon_*$ in spite of the fact that it implements exactly the same feedback prescription as Dahlia. The latter is qualitatively similar to a delayed cooling scheme used by Rosdahl et al. (2017) and others (Stinson et al. 2006; Teyssier et al. 2013). The lower efficiency $\epsilon_*$ is a consequence of chemistry. As under non-equilibrium conditions the gas must be denser to form H$_2$, the ISM becomes more clumpy (Fig. 6). These clumps can form massive clusters of OB stars which, acting coherently, yield stronger feedback and may disrupt completely the star-forming site.

## 4.2 Far- and mid-infrared emission

A meaningful way to compare the two galaxies is to predict their C II and H$_2$ line emission, which can be observable at high-z with ALMA, and possibly with SPICA (Spinoglio et al., in preparation; Egami et al., in preparation), respectively. Similarly to P17, we use a modified version of the [C II] emission model from Vallini et al. (2015, hereafter V15). Such model is based on temperature, density and metallicity grids built using CLOUDY (Ferland et al. 2013), as detailed in Appendix B. In Fig. 11, we plot the [C II] 157.74 µm and H$_2$ 17.03 µm surface brightness maps ($S/(L_{\odot} \text{ kpc}^{-2})$); the field of view is the same as in Fig. 5.

### 4.2.1 Far-infrared emission

Let us analyse first the C II emission. Dahlia has a [C II] luminosity of $\log(L_{CII}/L_{\odot}) \simeq 7.5$ which is about seven times smaller than Althæa, i.e. $\log(L_{CII}/L_{\odot}) \simeq 8.3$. Fig. 11 shows that the surface brightness morphology in the two galaxies is similar. Dahlia's emission is concentrated in the disc, featuring an average surface brightness of $\log S_{[CII]}/(L_{\odot} \text{ kpc}^{-2}) \simeq 6.4$ with peaks up to $\log S_{[CII]}/(L_{\odot} \text{ kpc}^{-2}) \simeq 7.4$ along the spiral arms. The analogous values for Althæa are 7.3 and 8.7, respectively.

This can be explained as follows. FIR emission from the warm ($\simeq 10^4$ K), low-density ($\lesssim 0.1 \text{ cm}^{-3}$) component of the ISM is suppressed at high-z by the CMB (Gong et al. 2012; da Cunha et al. 2013; Pallottini et al. 2015, V15; Appendix B), as the upper levels of the [C II] transition cannot be efficiently populated through collisions and the spin temperature of the transition approaches the CMB one (see Pallottini et al. 2015, for possibility of [C II] detection from low-density gas via CMB distortions). Thus, $\simeq 95$ per cent of the [C II] emission comes from dense ($\gtrsim 10 \text{ cm}^{-3}$, cold ($\simeq 100$ K), mostly molecular disc gas. As noted in V15 (see in particular their fig. 4) even when the CMB effect is neglected, the diffuse gas ($\lesssim 0.1 \text{ cm}^{-3}$) account only for $\lesssim 5$ per cent of the emission for galaxies with SFR $\sim 100 \text{ M}_{\odot} \text{ yr}^{-1}$ and $Z \sim Z_{\odot}$, while it can be important in smaller objects (Olsen et al. 2015). The emissivity (in $L_{\odot}/M_{\odot}$) of such gas can be written as in P17 (in equation 8, see also Vallini et al. 2013, 2017; Goicoechea et al. 2015):

$$ \epsilon_{[CII]} \simeq 0.1 \left( \frac{n}{10^2 \text{ cm}^{-3}} \right) \left( \frac{Z}{Z_{\odot}} \right), \qquad (10) $$

for $n \lesssim 10^3 \text{cm}^{-3}$, i.e. the critical density for [C II] emission.<fnref n="19" /> As the metallicity in the disc of the two galaxies is roughly similar ($\langle Z \rangle \simeq 0.5 Z_{\odot}$, see Table 1), difference in the luminosities is entirely explained by the larger density in Althæa. We stress once again that such density variation is a result of a more precise, non-equilibrium chemical network requiring to reach much higher densities before the gas is converted to stars. It is precisely that dense gas accounts for a larger FIR line emissivity from PDRs.

<footnote n="19">¹⁹ As the suppression of the CMB affects only the diffuse component ($\lesssim 0.1 \text{ cm}^{-3}$), no significant difference is expected in the emissivity from the discs of the two galaxies (equation 10), that is composed of much higher ($\gtrsim 20 \text{ cm}^{-3}$) density material.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Synthetic emission maps¹⁶ of the simulated galaxies Dahlia (left-hand panels) and Althæa (right-hand panels) at age $t_* \simeq 700$ Myr ($z=6$). Integrated surface brightness of [C II] ($S_{[CII]}/(L_\odot \text{ kpc}^{-2})$) and H₂ ($S_{H2}/(L_\odot \text{ kpc}^{-2})$) are shown in the upper and lower panels, respectively. The field of view is the same as in Fig. 5.</figcaption></figure>

We can also compare the calculated synthetic [C II] emission versus SFR with observations (Fig. 12) obtained for dwarf galaxies (De Looze et al. 2014), and available high-z detections or upper limits. The [C II] emission from Dahlia is lower than expected based on the local [C II]–SFR relation; its luminosity is also well below all upper limits for high-z galaxies. Although Althæa is ≃10 times more luminous, even this object lies below the local relation, albeit only by $1.3\sigma$. We believe that the reduced luminosity is caused by the combined effects of the CMB suppression and relatively lower Z. Note, however, that the predicted luminosity exceeds the upper limits derived for LAEs (e.g. Ouchi et al. 2013; Ota et al. 2014), but is broadly consistent with that of the handful of LBGs so far detected, like e.g. the four galaxies in the Pentericci et al. (2016). In general, observations are still rather sparse, with few [C II] detections with SFR comparable to Althæa (e.g. Capak et al. 2015). Also unclear is the amplitude of the scattering of the relation for high-z objects

compared with local ones. Improvements in the understanding of the ISM structure are expected from deeper observations and/or other ions (e.g. [O III]); Inoue et al. 2016; Carniani et al. 2017). Also helpful would be a larger catalogue of simulated galaxies (cf. Ceverino, Glover & Klessen 2017), to control environmental effects.

### 4.2.2 Mid-infrared emission

By inspecting the lower panel of Fig. 11 showing the predicted H₂ 17.03 µm line emission, we come to conclusions similar to those for the [C II]. Althæa outshines Dahlia by ≃15 × by delivering a total line luminosity of log (L<sub>H2</sub>/L<sub>☉</sub>) ≃ 6.5. Differently from the [C II] case, also the deviations from the mean are much more marked in Althæa, as appreciated from the figure.

Note that the H₂ 17.03 µm line emissivity is enhanced in high-density, high-temperature regions. Indeed, H₂ emission mostly
<figure><img src="image_13.png" /></figure>

**Figure 12.** The [C II]–SFR relation. Shown are Althæa (orange star) and Dahlia (black) at 700 Myr or $z=6$; the errors refer to rms variation in the last 50 Myr. Lines refer to results from the V15 model: constant metallicity models with $Z=Z_{\odot}$ (solid black), $Z=0.2Z_{\odot}$ (solid orange), $Z=0.05Z_{\odot}$ (pink dashed) and a model with mean $(Z/Z_{\odot}) = 0.05 + \text{density–metallicity relation}$ extracted from cosmological simulations (Pallottini, Gallerani & Ferrara 2014b, blue dot–dashed). Data for local dwarf galaxies (De Looze et al. 2014) are plotted with little circles and the grey hatched region gives the mean and rms variation in the sample. For high-$z$ galaxies, detections (upper-limits) are plotted with filled (empty) symbols, according to the inset legend. The high-$z$ sample include individual galaxies as BDF-3299 (Maiolino et al. 2015; Carniani et al. 2017), HCM6A (Kanekar et al. 2013), Himiko (Ouchi et al. 2013; Ota et al. 2014), IOK-1 (Ota et al. 2014) and data from Capak et al. (2015, $z \simeq 5.5$), Willott et al. (2015, $z \simeq 6$), Schaerer et al. (2015, $z \simeq 7$), (Pentericci et al. 2016, $z \simeq 7$), González-López et al. (2014, $z \simeq 8$) and lensed $z \simeq 6.5$ galaxies from Knudsen et al. (2016) and Bradac et al. (2017).

arises from shock-heated molecular gas, for which $100 \lesssim n/cm^{-3} \lesssim 10^5$ and for $10 \lesssim T/K \lesssim 3000$ (see Appendix B).

For Dahlia, the disc density is relatively low, $n \simeq 30 cm^{-3}$; in addition only 20 per cent of the gas is warm enough to allow some ($\epsilon_{H2} \simeq 0.01 L_{\odot} M_{\odot}^{-1}$) emission. In practice, such emission predominantly occurs along the outer spiral arms of the galaxy where these conditions are met due to the heating produced by SN explosions. In the denser Althæa disc, the gas emissivity can attain $\epsilon_{H2} \simeq 10^{-3}-0.01$ already at moderate $T = 200$ K. The brightness peaks are associated with a few ($\lesssim 1$ per cent) pockets of thousand-degree gas; they can be clearly identified in the map. This is particularly interesting because a galaxy-like Althæa might be detectable at very high-$z$ with SPICA, as suggested by Egami et al. (in preparation).

## 5 CONCLUSIONS

To improve our understanding of high-$z$ galaxies we have studied the impact of $H_2$ chemistry on their evolution, morphology and observed properties. To this end, we compare two zoom-in galaxy simulations implementing different chemical modelling. Both simulations start from the cosmological same initial conditions, and follow the evolution of a prototypical $M_* \simeq 10^{10} M_{\odot}$ galaxy at $z=6$ resolved at the scale of giant MCs (30 pc). Stars are formed according to an $H_2$ dependent SK relation. We also account for winds from massive stars, SN explosions and radiation pressure in

a stellar age/metallicity dependent fashion (see Section 2.1). The first galaxy is named Dahlia and $H_2$ formation is computed from the Krumholz et al. (2009) equilibrium model; Althæa instead implements a non-equilibrium chemistry network, following (B16). The key results can be summarized as follows:

(a) The SFR of the two galaxies is similar, and increases with time reaching values close to $100 M_{\odot} yr^{-1}$ at $z=6$ (see Fig. 2). However, Dahlia forms stars at a rate that is on average $1.5 \pm 0.6$ times larger than Althæa; it also shows a less prominent burst structure.

(b) Both galaxies at $z=6$ have an SFR-stellar mass relation compatible with J16 observations (Fig. 3). Moreover, they both show a continuous time evolution from specific SFR of $\text{sSFR} \simeq 40$ to $5 \text{ Gyr}^{-1}$. This is understood as an effect of the progressively increasing impact of stellar feedback hindering subsequent SF events.

(c) The non-equilibrium chemical model implemented in Althæa determines the atomic to molecular hydrogen transition to occur at densities exceeding $300 \text{ cm}^{-3}$, i.e. about 10 times larger than predicted by equilibrium model used for Dahlia (Fig. 1). As a result, Althæa features a more clumpy and fragmented morphology (Fig. 6). This configuration makes SN feedback more effective, as noted in point (a) above (Fig. 4).

(d) Because of the lower density and weaker feedback, Dahlia sits $3\sigma$ away from the SK relation; Althæa, instead nicely agrees with observations (Fig. 9). Note that although the SF efficiency is similar in the two galaxies and consistent with other simulations (Semenov et al. 2016), Dahlia is off the relation because of insufficient molecular gas content (Fig. 8).

(e) We confirm that most of the emission from the C II and $H_2$ is due to the dense gas forming the disc of the two galaxies. Because of Dahlia’s lower average density, Althæa outshines Dahlia by a factor of 7 (15) in [C II] 157.74 µm ($H_2$ 17.03 µm) line emission (Fig. 11). Yet, Althæa has a 10 times lower [C II] luminosity than expected from the locally observed [C II]–SFR relation (Fig. 12). Whether this relation does not apply at high-$z$ or the line luminosity is reduced by CMB and metallicity effects remains as an open question that can be investigated with future deeper observations.

To conclude, both Dahlia and Althæa follow the observed high-$z$ SFR-$M_*$ relation. However, many other observed properties (SK relation, C II and $H_2$ emission) are very different. This shows the importance of accurate, non-equilibrium implementation of chemical networks in early galaxy numerical studies.

## ACKNOWLEDGEMENTS

We are grateful to the participants of *The Cold Universe* programme held in 2016 at the KITP, UCSB, for discussions during the workshop. We thank P. Capelo, D. Celoria, E. Egami, D. Galli, T. Grassi, L. Mayer, S. Riolo and J. Wise for interesting and stimulating discussions. We thank the authors and the community of *RAMSES* and *PYMSES* for their work. AP acknowledges support from Centro Fermi via the project CORTES, ‘Cosmological Radiative Transfer in Early Structures’. AF acknowledges support from the ERC Advanced Grant INTERSTELLAR H2020/740120. RM acknowledges support from the ERC Advanced Grant 695671 ‘QUENCH’ and from the Science and Technology Facilities Council (STFC). SS acknowledges support from the European Commission through a Marie Skłodowska-Curie Fellowship, program PRIMORDIAL, Grant No. 700907. This research was supported in part by the National Science Foundation under Grant No. NSF PHY11-25915.
REFERENCES

Agertz O., Kravtsov A. V., 2015, ApJ, 804, 18

Agertz O., Kravtsov A. V., Leitner S. N., Gnedin N. Y., 2013, ApJ, 770, 25

Asano R. S., Takeuchi T. T., Hirashita H., Inoue A. K., 2013, Earth Planets
Space, 65, 213

Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481

Bakes E. L. O., Tielens A. G. G. M., 1994, ApJ, 427, 822

Barai P., Monaco P., Murante G., Ragagnin A., Viel M., 2015, MNRAS,
447, 266

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

Bertelli G., Bressan A., Chiosi C., Fagotto F., Nasi E., 1994, A&AS, 106, 275

Black J. H., 1987, in Hollenbach D. J., Thronson H. A. Jr, eds, Astrophysics and Space Science Library Vol. 134, Interstellar Processes. Reidel, Dordrecht, p. 731

Black J. H., Dalgarno A., 1976, ApJ, 203, 132

Bouché N. et al., 2007, ApJ, 671, 303

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

Bovino S., Grassi T., Capelo P. R., Schleicher D. R. G., Banerjee R., 2016, A&A, 590, A15 (B16)

Bradac M. et al., 2007, ApJ, 836, L2

Bryan G. L. et al., 2014, ApJS, 211, 19

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

Carilli C. L., Walter F., 2013, ARA&A, 51, 105

Carniani S. et al., 2014, preprint (arXiv:1701.03468)

Cazaux S., Spaans M., 2009, A&A, 496, 365

Ceverino D., Klypin A., Klimek E. S., Trujillo-Gomez S., Churchill C. W., Primack J., Dekel A., 2014, MNRAS, 442, 1545

Ceverino D., Glover S., Klessen R., 2017, MNRAS, 470, 2791

Choudhury T. R., Padmanabhan T., Srianand R., 2001, MNRAS, 322, 561

Ciardi B., Ferrara A., 2001, MNRAS, 324, 648

Coles P., Jones B., 1991, MNRAS, 248, 1

Cormier D. et al., 2015, A&A, 578, A53

da Cunha E. et al., 2013, ApJ, 766, 13

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

Daddi E. et al., 2010b, ApJ, 714, L118

Davé R., Finlator K., Oppenheimer B. D., 2011, MNRAS, 416, 1354

De Looze I. et al., 2014, A&A, 568, A62

Draine B. T., 1978, ApJS, 36, 595

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

Dunlop J. S., 2013, in Wiklind T., Mobasher B., Bromm V., eds, Astrophysics and Space Science Library Vol. 396, The First Galaxies. Springer-Verlag, Berlin, p. 223

Federrath C., Klessen R. S., 2013, ApJ, 763, 51

Ferland G. J. et al., 2013, Rev. Mex. Astron. Astrofis., 49, 137

Fiacconi D., Mayer L., Madau P., Lupi A., Dotti M., Haardt F., 2017, MNRAS, 467, 4080

Fujimoto S., Ouchi M., Shibuya T., Nagai H., 2017, preprint (arXiv:1703.02138)

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

Galli D., Palla F., 1998, A&A, 335, 403

Gatto A. et al., 2015, MNRAS, 449, 1057

Genzel R. et al., 2010, MNRAS, 407, 2091

Glassgold A. E., Galli D., Padovani M., 2012, ApJ, 756, 157

Gleser L., Nusser A., Ciardi B., Desjacques V., 2006, MNRAS, 370, 1329

Gnedin N. Y., 2010, ApJ, 721, L79

Goicoechea J. R. et al., 2015, ApJ, 812, 75

Gong Y., Cooray A., Silva M., Santos M. G., Bock J., Bradford C. M., Zemcov M., 2012, ApJ, 745, 49

González V., Labbé L., Bouwens R. J., Illingworth G., Franx M., Kriek M., Brammer G. B., 2010, ApJ, 713, 115

González-López J. et al., 2014, ApJ, 784, 99

Grassi T., Bovino S., Schleicher D. R. G., Prieto J., Seifried D., Simoncini E., Gianturco F. A., 2014, MNRAS, 439, 2386

Grassi T., Bovino S., Haugbølle T., Schleicher D. R. G., 2017, MNRAS, 466, 1259

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

Habing H. J., 1968, Bull. Astron. Inst. Neth., 19, 421

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

Hartwig T., Clark P. C., Glover S. C. O., Klessen R. S., Sasaki M., 2015, ApJ, 799, 114

Heiderman A., Evans II N. J., Allen L. E., Huard T., Heyer M., 2010, ApJ, 723, 1019

Hirashita H., Ferrara A., 2002, MNRAS, 337, 921

Hollenbach D., McKee C. F., 1979, ApJS, 41, 555

Hopkins P. F., Narayanan D., Murray N., 2013, MNRAS, 432, 2647

Hopkins P. F. et al., 2017, preprint (arXiv:1702.06148)

Inoue A. K. et al., 2016, Science, 352, 1559

Jiang L. et al., 2016, ApJ, 816, 16 (J16)

Jura M., 1975, ApJ, 197, 575

Kanekar N., Wagg J., Ram Chary R., Carilli C. L., 2013, ApJ, 771, L20

Katz H., Kimm T., Sijacki D., Haehnelt M., 2016, preprint (arXiv:1612.01786)

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

Kennicutt R. C., Evans N. J., 2012, ARA&A, 50, 531

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

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

Klessen R. S., Glover S. C. O., 2016, Star Formation in Galaxy Evolution: Connecting Numerical Models to Reality. Springer-Verlag, Berlin

Knudsen K. K., Richard J., Kneib J.-P., Jauzac M., Clément B., Drouart G., Egami E., Lindroos L., 2016, MNRAS, 462, L6

Kroupa P., 2001, MNRAS, 322, 231

Krumholz M. R., 2015, preprint (arXiv:1511.03457)

Krumholz M. R., McKee C. F., Tumlinson J., 2008, ApJ, 689, 865

Krumholz M. R., McKee C. F., Tumlinson J., 2009, ApJ, 693, 216

Krumholz M. R., Dekel A., McKee C. F., 2012, ApJ, 745, 69

Labadens M., Chapon D., Pomaréde D., Teyssier R., 2012, in Ballester P., Egret D., Lorente N. P. F., eds, ASP Conf. Ser. Vol. 461, Astronomical Data Analysis Software and Systems XXI. Astron. Soc. Pac., San Francisco, p. 837

Lada C. J., Lombardi M., Alves J. F., 2010, ApJ, 724, 687

Laporte N. et al., 2017, ApJ, 837, L21

Leitherer C. et al., 1999, ApJS, 123, 3

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

McKee C. F., Krumholz M. R., 2010, ApJ, 709, 308

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

Maio U., Tescari E., 2015, MNRAS, 453, 3798

Maio U., Petkova M., De Lucia G., Borgani S., 2016, MNRAS, 460, 3733

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

Martizzi D., Faucher-Giguère C.-A., Quataert E., 2015, MNRAS, 450, 504

Mecke K. R., Buchert T., Wagner H., 1994, A&A, 288, 697

Murray N., 2011, ApJ, 729, 133

Novikov D., Schmalzing J., Mukhanov V. F., 2000, A&A, 364, 17

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

Olsen K. P., Greve T. R., Narayanan D., Thompson R., Toft S., Brinch C., 2015, ApJ, 814, 76

Osterbrock D. E., 1989, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei. University Science Books, Mill Valley, CA

Ostriker J. P., McKee C. F., 1988, Rev. Mod. Phys., 60, 1

Ota K. et al., 2014, ApJ, 792, 34

Ouchi M. et al., 2013, ApJ, 778, 102

Padoan P., Haugbølle T., Nordlund Å., 2012, ApJ, 759, L27

Pallottini A., Ferrara A., Gallerani S., Salvadori S., D'Odorico V., 2014a, MNRAS, 440, 2498

Pallottini A., Gallerani S., Ferrara A., 2014b, MNRAS, 444, L105

Pallottini A., Gallerani S., Ferrara A., Yue B., Vallini L., Maiolino R., Feruglio C., 2015, MNRAS, 453, 1898

Pallottini A., Ferrara A., Gallerani S., Vallini L., Maiolino R., Salvadori S., 2017, MNRAS, 465, 2540 (P17)

Pentericci L. et al., 2016, ApJ, 829, L11

Pérez-Montero E., 2017, PASP, 129, 043001

Petkova M., Maio U., 2012, MNRAS, 422, 3067

Planck Collaboration XVI, 2014, A&A, 571, A16
Rahmati A., Pawlik A. H., Raicevic M., Schaye J., 2013, MNRAS, 430, 2427

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

Rey-Raposo R., Dobbs C., Agertz O., Alig C., 2017, MNRAS, 464, 3536

Richings A. J., Schaye J., Oppenheimer B. D., 2014, MNRAS, 442, 2780

Rosdahl J., Schaye J., Teyssier R., Agertz O., 2015, MNRAS, 451, 34

Rosdahl J., Schaye J., Dubois Y., Kimm T., Teyssier R., 2017, MNRAS, 466, 11

Roskar R., Teyssier R., Agertz O., Wetzstein M., Moore B., 2014, MNRAS, 444, 2837

Roussel H. et al., 2007, ApJ, 669, 959

Scannapieco C. et al., 2012, MNRAS, 423, 1726

Schaerer D., Boone F., Zamojski M., Staguhn J., Dessauges-Zavadsky M., Finkelstein S., Combes F., 2015, A&A, 574, A19

Schmalzing J., Buchert T., 1997, ApJ, 482, L1

Schmalzing J., Gorski K. M., 1998, MNRAS, 297, 355

Schmidt M., 1959, ApJ, 129, 243

Semenov V. A., Kravtsov A. V., Gnedin N. Y., 2016, ApJ, 826, 200

Shibuya T., Ouchi M., Harikane Y., 2015, ApJS, 219, 15

Smith B. D. et al., 2017, MNRAS, 466, 2217

Stanway E. R., 2017, preprint (arXiv:1702.07303)

Stasińska G., 2007, preprint (arXiv:0704.0348)

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

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

Teyssier R., Pontzen A., Dubois Y., Read J. I., 2013, MNRAS, 429, 3068

Timmermann R., Bertoldi F., Wright C. M., Drapatz S., Draine B. T., Haser L., Sternberg A., 1996, A&A, 315, L281

Tomassetti M., Porciani C., Romano-Díaz E., Ludlow A. D., 2015, MNRAS, 446, 3330

Torrey P., Hopkins P. F., Faucher-Giguère C.-A., Vogelsberger M., Quataert E., Keres D., Murray N., 2017, MNRAS, 467, 2301

Turner J., Kirby-Docken K., Dalgarno A., 1977, ApJS, 35, 281

Valle G., Ferrini F., Galli D., Shore S. N., 2002, ApJ, 566, 252

Vallini L., Gallerani S., Ferrara A., Baek S., 2013, MNRAS, 433, 1567

Vallini L., Gallerani S., Ferrara A., Pallottini A., Yue B., 2015, ApJ, 813, 36 (V15)

Vallini L., Ferrara A., Pallottini A., Gallerani S., 2017, MNRAS, 467, 1300

Verner D. A., Ferland G. J., 1996, ApJS, 103, 467

Wakelam V. et al., 2012, ApJS, 199, 21

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

Weaver R., McCray R., Castor J., Shapiro P., Moore R., 1977, ApJ, 218, 377

Webber W. R., 1998, ApJ, 506, 329

Weingartner J. C., Draine B. T., 2001, ApJ, 563, 842

Willott C. J., Carilli C. L., Wagg J., Wang R., 2015, ApJ, 807, 180

Wise J. H., Abel T., Turk M. J., Norman M. L., Smith B. D., 2012, MNRAS, 427, 311

Wolcott-Green J., Haiman Z., Bryan G. L., 2011, MNRAS, 418, 838

Wolfire M. G., McKee C. F., Hollenbach D., Tielens A. G. G. M., 2003, ApJ, 587, 278

Yoshiura S., Shimabukuro H., Takahashi K., Matsubara T., 2017, MNRAS, 465, 394

## APPENDIX A: MINKOWSKY FUNCTIONALS

In general, Minkowsky functionals are mathematical tools that give a complete characterization of the morphology of an IR$^n \mapsto$ IR field. In astrophysics they have been proposed as a mean to give a description of the large-scale structure (e.g. Mecke, Buchert & Wagner 1994), to study the topology of H\ddot{H} bubbles for reionization studies (e.g. Gleser et al. 2006; Yoshiura et al. 2017), and (for $n=2$ fields) analyse CMB anisotropies and non-Gaussianity (e.g. Schmalzing & Gorski 1998; Novikov, Schmalzing & Mukhanov 2000).

A formal definition can be given following Schmalzing & Buchert (1997). Let $u(\mathbf{x})$ denote a scalar field defined on a subset of $\mathbb{R}^3$ with volume $V$. Let us take $u$ such that it has zero mean $\langle\langle u(\mathbf{x}) \rangle\rangle = 0$) and variance $\langle\langle u^2(\mathbf{x}) \rangle\rangle = \sigma^2$. Then, we can define the excursion set $F_\nu(\mathbf{x})$ as the ensemble of regions in $V$ satisfying $u(\mathbf{x}) > \nu\sigma$. Then, the Minkowsky functionals can be defined in terms of volume ($d^3\mathbf{x}$) and surface ($d^2\mathbf{x}$) integrals as a function of the threshold $\nu$:

$$V_0(\nu) = (V)^{-1} \int_V \Theta(u - \nu\sigma) d^3\mathbf{x} \quad (A1a)$$

$$V_1(\nu) = (6V)^{-1} \int_{\partial F_\nu} d^2\mathbf{x} \quad (A1b)$$

$$V_2(\nu) = (6\pi V)^{-1} \int_{\partial F_\nu} (\kappa_1 + \kappa_2) d^2\mathbf{x} \quad (A1c)$$

$$V_3(\nu) = (4\pi V)^{-1} \int_{\partial F_\nu} \kappa_1 \kappa_2 d^2\mathbf{x}, \quad (A1d)$$

where $\Theta$ is the Heaviside function, $\partial F_\nu(\mathbf{x})$ is the surface beading the excursion set $F_\nu$ and $\kappa_1(\mathbf{x})$ and $\kappa_2(\mathbf{x})$ are the two principal curvatures of the surface. In practical terms, $V_0$ is a measure of the volume filling factor of the excursion set with threshold $\nu$, $V_1$ of the surface, $V_2$ of the mean curvature (sphericity/concavity) and $V_3$ of the Euler characteristic (shape of components).

The curvatures on the surface $\partial F_\nu$ can be expressed via the Koenderink invariant (Gleser et al. 2006, see appendix A and reference therein): by adopting the Einstein sum convention we can write

$$\kappa_1 + \kappa_2 = \epsilon^{ijk} \epsilon^{lmn} \delta_{kn}(\partial_i u)(\partial_j \partial_l u)(\partial_m u)/N_t^{3/2} \quad (A2a)$$

$$2\kappa_1\kappa_2 = \epsilon^{ijk} \epsilon^{lmn} (\partial_i u)(\partial_j u)(\partial_k u)(\partial_l u)/N_t^2 \quad (A2b)$$

$$N_t = (\partial_p u)(\partial^p u), \quad (A2c)$$

where $\epsilon^{ijk}$ is the Levi-Civita symbol, $\delta_{ij}$ is the Kronecker delta and $\partial_i$ is the $i$th component of the partial derivative operator. Finally, equations (A1) can be suitably expressed as integral over the volume by using the following relation:

$$\int_{\partial F_\nu} d^2\mathbf{x} = \int_V \delta(u - \nu\sigma) N_t^{1/2} d^3\mathbf{x}, \quad (A3)$$

where $\delta$ is the Dirac delta.

As an illustrative example, we can compute the Minkowsky functionals for a zero-mean Gaussian random field $u \equiv \log \Delta$ with variance $\sigma^2 = \langle\langle \log \Delta^2 \rangle\rangle$ and variance of the tangent field $\sigma_t^2 = \langle\langle (\partial_p \log \Delta)(\partial^p \log \Delta) \rangle\rangle$. For such Gaussian field, the Minkowsky functionals can be expressed using the following analytical expression (Schmalzing & Buchert 1997, see also Gleser et al. 2006)

$$V_0 = 1/2 - c_0 \int_0^\nu \exp(-x^2/2)dx \quad (A4a)$$

$$V_1 = c_1 \lambda \exp(-\nu^2/2) \quad (A4b)$$

$$V_2 = c_2 \lambda^2 \nu \exp(-\nu^2/2) \quad (A4c)$$

$$V_3 = c_3 \lambda^3 (\nu^2 - 1) \exp(-\nu^2/2), \quad (A4d)$$

where $\nu = \log \Delta/\sigma$, $\lambda = (6\pi)^{-1/2}\sigma_t/\sigma$, and $c_i$ are numerical constant with values $c_0 = c_3 = (2\pi)^{-1/2}$ and $c_1 = c_2 = 2/(3(2\pi)^{1/2})$.

We numerically compute the Minkowsky functionals for the log $\Delta$ field with $\sigma = 1$ on a $256^3$ unigrid box with volume (10 Mpc)$^3$, which thus result in a tangent field variance of $\sigma_t \simeq 5.6$ kpc. The
<figure><img src="image_14.png" /><figcaption>**Figure A1.** Example of Minkowsky functionals calculated for the zero-mean Gaussian random field log $\Delta$. The functionals are plotted with an orange line and transparent region and are normalized by $V_i/\lambda^i$, with $\lambda = (6\pi)^{-1/2}\sigma_t/\sigma$, where $\sigma = 1$ is the field variance and $\sigma_t \simeq 5.6/kpc$ is the variance of the tangent field. Analytical expected values for the functionals (equations A4a) are plotted with red crosses. To guide the eye, $\log \Delta = 0$ is marked with a dotted vertical black line, and both $V_2 = 0$ and $V_3 = 0$ are highlighted with dashed horizontal lines.</figcaption></figure>

resulting Minkowsky functionals are plotted in Fig. A1. We find a very good match with the analytical values.

Since the chosen log $\Delta$ Gaussian field is an approximation to the quasi-linear regime of the cosmic density field (e.g. Coles & Jones 1991; Choudhury, Padmanabhan & Srianand 2001), it is intuitive to analyse the properties of its Minkowsky functionals. In Fig. A1, the filling factor $V_0$ gives the probability of finding regions with increasing overdensity $\Delta$; note that at $\Delta = 1$ (mean density), $V_0 = 0.5$, i.e. it is equiprobable to find voids ($\log \Delta \lesssim -1$) and overdense regions ($\log \Delta \gtrsim 1$). Both voids and overdense regions are isolated $V_3(\gtrsim 1) = V_3(\lesssim -1) > 0$ and have a smaller area with respect to mean density regions ($V_1(\gtrsim 1) = V_1(\lesssim -1) \lesssim V_1(0)$). However, while overdensities have spherical shapes $V_2(\gtrsim 1) > 0$, voids are concave regions ($V_2(\lesssim 1) < 0$): both voids and overdense regions are delimited by connected ($V_2(\simeq 0) < 0$) mean density regions, which are almost flat ($V_2(\simeq 0) \simeq 0$) and have very large surface areas.

## APPENDIX B: EMISSION FROM C II AND H$_2$

To compute the emission from C II ions and H$_2$ molecules, we post-process the simulation outputs using the photoionization code CLOUDY (Ferland et al. 2013), similarly to what done in V15 and P17. We consider a grid of models based on the density (n), temperature (T) and metallicity (Z) of the gas in our simulation. We produce a total of $10^3$ models that are parametrized as a function of the column density (N). For each model, we adopt a plane-parallel geometry and assume a dust content proportional to the metallicity.

The radiation field includes the CMB background and an ISRF produced by stars, which is obtained by rescaling the MW spectrum (Black 1987) using the main galaxy (Dahlia or Althæa) SFR. At $z=6$, Dahlia has a star formation rate SFR = $156 M_{\odot} \text{ yr}^{-1}$ and Althæa has SFR = $136 M_{\odot} \text{ yr}^{-1}$, where the uncertainty is the rms in the last 50 Myr (Section 3.1, in particular see Fig. 2). For modelling convenience, in the CLOUDY calculation we set SFR = $100 M_{\odot} \text{ yr}^{-1}$. Note that a larger value for the rescaling does not yield a large variation of the expected [C II] in molecular gas (Vallini et al. 2017, with $G = G_0 \text{ SFR}/M_{\odot} \text{ yr}^{-1}$), and H$_2$ emission is relatively unaffected by the field, as the excitation is mostly due to shocks (e.g. Black & Dalgarno 1976; Ciardi & Ferrara 2001).

As noted in Section 2.3, accounting for the UVB is not relevant for the ionization state of the gas in the proximity of galaxies (Gnedin 2010). Thus, in our CLOUDY models we consider that the gas is shielded by a column density of $N \simeq 10^{20} \text{ cm}^{-2}$.

Regarding the [C II], we underline that the effect of CMB suppression of [C II] is included in the V15 model (see also da Cunha et al. 2013; Pallottini et al. 2015). Such effects suppress the emission where the spin temperature of the [C II] transition is close to the CMB one. This is relevant for low-density ($n \lesssim 10^{-1} \text{ cm}^{-3}$) medium that does not have enough collisions to decouple from the CMB. Here, we do not account for the photoevaporation effect on MC, which has an important impact on FIR emission (Vallini et al. 2017), particularly when including a spatially varying FUV field, not included in the present modelling.

Even though in this paper, we only show the H$_2$ line at $\lambda = 17.04 \text{ µm}$ (Section 4, in particular see Fig. 11), we used CLOUDY to compute the emission of the following H$_2$ roto-vibrational lines: 0–0 S(0), 0–0 S(1), 0–0 S(5) and 1–0 S(1), which correspond to transition of wavelength $\lambda/\mu\text{m} = 2.12, 6.91, 9.66 \text{ 17.04, and } 28.22$, respectively (see Spinoglio et al. in preparation; Egami et al. in preparation).

For the considered five transitions, the oscillator strength is a decreasing function of $\lambda$ (e.g. Turner, Kirby-Docken & Dalgarno 1977), going from $\simeq 3 \times 10^{-7} \text{ s}^{-1}$ for $\lambda = 2.12 \text{ µm}$ to $\simeq 3 \times 10^{-11} \text{ s}^{-1}$ for $\lambda = 28.22 \text{ µm}$. On the other hand, both the excitation temperature (Timmermann et al. 1996) and the critical density for collisional excitation (Roussel et al. 2007) decrease for decreasing $\lambda$ (see also Black & Dalgarno 1976), e.g. for $\lambda = 2.12 \text{ µm}$ we have an excitation temperature $T_{\text{ex}} \simeq 6 \times 10^4 \text{ K}$ and a critical density $n_{\text{cr}} \simeq 10^4 \text{ cm}^{-3}$, while for $\lambda = 28.22 \text{ µm}$ we have $T_{\text{ex}} \simeq 5 \times 10^3 \text{ K}$ and $n_{\text{cr}} \simeq 5 \text{ cm}^{-3}$. In the simulations, the bulk of the H$_2$ gas have $\langle T \rangle \sim 10^2 \text{ K}$ and density in the range $10^2 \lesssim n/\text{cm}^{-3} \lesssim 10^3$ (see Fig. 8); thus in both Dahlia and Althæa the most favoured H$_2$ transition is the $\lambda = 17.03 \text{ µm}$, followed by the $28.21 \text{ µm}$ line.

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