# Zooming on the internal structure of z 6 galaxies

A. Pallottini, A. Ferrara, S. Gallerani, L. Vallini, R. Maiolino, S. Salvadori

► To cite this version:

A. Pallottini, A. Ferrara, S. Gallerani, L. Vallini, R. Maiolino, et al. Zooming on the internal structure of z 6 galaxies. Monthly Notices of the Royal Astronomical Society, 2017, 465 (3), pp.2540-2558. 10.1093/mnras/stw2847. hal-01554208

HAL Id: hal-01554208

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

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>

# Zooming on the internal structure of $z \simeq 6$ galaxies

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

Accepted 2016 November 2. Received 2016 November 1; in original form 2016 July 25

## ABSTRACT

We present zoom-in, adaptive mesh refinement, high-resolution ($\simeq$30 pc) simulations of high-redshift ($z \simeq 6$) galaxies with the aim of characterizing their internal properties and interstellar medium. Among other features, we adopt a star formation model based on a physically sound molecular hydrogen prescription, and introduce a novel scheme for supernova feedback, stellar winds and dust-mediated radiation pressure. In the zoom-in simulation, the target halo hosts ‘Dahlia’, a galaxy with a stellar mass $M_* = 1.6 \times 10^{10} M_{\odot}$, representative of a typical $z \sim 6$ Lyman-break galaxy. Dahlia has a total $H_2$ mass of $10^{8.5} M_{\odot}$ that is mainly concentrated in a disc-like structure of effective radius $\simeq 0.6$ kpc and scale height $\simeq 200$ pc. Frequent mergers drive fresh gas towards the centre of the disc, sustaining a star formation rate per unit area of $\simeq 15 M_{\odot} \text{ yr}^{-1} \text{ kpc}^{-2}$. The disc is composed of dense ($n \gtrsim 25 \text{ cm}^{-3}$), metal-rich ($Z \simeq 0.5 Z_{\odot}$) gas that is pressure supported by radiation. We compute the 158 $\mu$m [C II] emission arising from Dahlia, and find that $\simeq 95$ per cent of the total [C II] luminosity ($L_{[C II]} \simeq 10^{7.5} L_{\odot}$) arises from the $H_2$ disc. Although 30 per cent of the C II mass is transported out of the disc by outflows, such gas negligibly contributes to [C II] emission, due to its low density ($n \lesssim 10 \text{ cm}^{-3}$) and metallicity ($Z \lesssim 10^{-1} Z_{\odot}$). Dahlia is underluminous with respect to the local [C II]–SFR relation; however, its luminosity is consistent with upper limits derived for most $z \sim 6$ galaxies.

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

## 1 INTRODUCTION

The discovery and characterization of primeval galaxies constitute some of the biggest challenges in current observational and theoretical cosmology.<fnref n="1" />

Deep optical/near-infrared (IR) surveys (Dunlop 2013; Madau & Dickinson 2014; Bouwens et al. 2015) have made impressive progresses in identifying galaxies well within the Epoch of Reionization ($z \simeq 6$). Such surveys yield key information about the star

formation (SF) of hundreds of galaxies in the early Universe. They also allow to statistically characterize galaxies in terms of their ultraviolet (UV) luminosity up to $z \sim 10$ (Bouwens et al. 2015). However, using these surveys with broad-band alone little can be learned about other properties, such as their gas and dust content, metallicity, interactions with the surrounding environment (e.g. Barnes, Garel & Kacprzak 2014), feedback (e.g. Dayal et al. 2014) and outflows (Gallerani et al. 2016).

To obtain a full picture of these systems, optical/IR surveys must be complemented with additional probes. Information on the metal content and energetic of the interstellar medium (ISM) can be obtained with observations of far-IR (FIR) fine structure lines, and in particular the [C II] ($^2P_{3/2} \rightarrow ^2P_{1/2}$) line at 157.74 $\mu$m. The [C II] line is the dominant coolant of the ISM being excited in different ISM phases, such as the diffuse cold neutral medium (CNM), warm neutral medium (WNM), high-density photodissociation regions (PDRs) and – to a lower extent – ionized gas (Tielens &

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

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

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

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

<footnote n="5">⁵Dipartimento di Fisica e Astronomia, Università di Bologna, viale Berti Pichat 6/2, I-40127 Bologna, Italy</footnote>

<footnote n="6">⁶INAF, Osservatorio Astronomico di Bologna, via Ranzani 1, I-40127 Bologna, Italy</footnote>

<footnote n="7">⁷GEPI, Observatoire de Paris, Place Jules Janssen, F-92195 Meudon, Paris, France</footnote>

<footnote>*E-mail: andrea.pallottini@sns.it, ap926@cam.ac.uk</footnote>

<footnote n="1">¹ In the following, we assume cosmological parameters compatible with Planck results, i.e. a $\Lambda$ cold dark matter ($\Lambda$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 h \text{ km s}^{-1} \text{ Mpc}^{-1}$ with $h = 0.678$, spectral index $n = 0.967$ and $\sigma_8 = 0.826$ (Planck Collaboration XVI 2014).</footnote>
Hollenbach 1985; Wolfire et al. 1995; Abel 2006; Vallini et al. 2013). As [C II] emission can be enhanced by shocks, it has been suggested as a good outflow tracer (e.g. Maiolino et al. 2012; Kreckel et al. 2014; Cicone et al. 2015; Janssen et al. 2016), and can thus in general be used to study feedback processes in galaxies.

Observationally, the [C II] line is a promising probe as it is often the brightest among FIR emission lines, accounting for up to ~1 per cent of the total IR luminosity of galaxies (e.g. Crawford et al. 1985; Madden et al. 1997). It has been successfully used to probe the low-z ISM (e.g. De Looze et al. 2014). The unprecedented sensitivity of the Atacama Large Millimeter/submillimeter Array (ALMA) makes it possible for the first time to use [C II] emission to characterize high-z galaxies. Before the ALMA advent, in fact, detections were limited to a handful of quasi-stellar object (QSO) host galaxies and rare galaxies with extreme SF rates (SFR $\simeq 10^3 M_\odot$ yr$^{-1}$; e.g. Maiolino et al. 2005; De Breuck et al. 2011; Gallerani et al. 2012; Carilli & Walter 2013; Cicone et al. 2015).

However, for 'normal' star-forming galaxies ($\lesssim 10^2 M_\odot$ yr$^{-1}$) at $z \sim 6-7$ early ALMA searches for [C II] lines have mostly yielded upper limits (e.g. Kanekar et al. 2013; Ouchi et al. 2013; Ota et al. 2014; Schaerer et al. 2015). The situation has changed recently with a number of robust [C II] detections (e.g. Capak et al. 2015; Maiolino et al. 2015; Willott et al. 2015; Knudsen et al. 2016).

In many cases, the high-z [C II] line luminosity is fainter than expected from the [C II]–SFR relation found in local galaxies (De Looze et al. 2014). To explain such [C II]–SFR deficit, some efforts have been devoted to model the [C II] emission from high-z galaxies (Nagamine, Wolfe & Hernquist 2006; Vallini et al. 2013, 2015; Muñoz & Furlanetto 2014; Olsen et al. 2015). In brief, these theoretical works show that the [C II]–SFR deficit can be ascribed to different effects.

(a) Lower metallicity of high-z galaxies (Vallini et al. 2013, 2015; Muñoz & Furlanetto 2014), in particular supported by observations of lensed galaxies (Knudsen et al. 2016).

(b) Suppression of the [C II] line around star-forming regions (Vallini et al. 2013), typically observed as a displacement of the [C II] with respect to the UV-emitting region, as seen e.g. in BDF 3299 (Maiolino et al. 2015) and in some of the Capak et al. (2015) galaxies. This would be a signature of stellar feedback heating/ionizing the putative [C II]-emitting gas.

(c) Suppression of the [C II] line by the increased cosmic microwave background (CMB) temperature in the WNM/CNM component (Pallottini et al. 2015; Vallini et al. 2015), similarly to what is observed for dust emission (da Cunha et al. 2013).

Simulating the ISM of early galaxies at sufficient resolution and including feedback effects might shed light on these questions. Feedback prescriptions are particularly important as such a process regulates the amount of (dense) gas likely to radiate most of the power detected with FIR lines. Several studies have explored optimal strategies to include feedback in galaxy simulations.

For some works, the interest is in the comparison between different kinds of stellar feedback prescription, as modelled via thermal and/or kinetic energy deposition in the gas from supernovae (SN), winds (Agertz et al. 2013; Hopkins et al. 2014; Agertz & Kravtsov 2015; Barai et al. 2015) and radiation pressure (Wise et al. 2012; Ceverino et al. 2014); other analyses focus on implementing complex chemical networks in simulations (Maio & Tescari 2015; Tomassetti et al. 2015; Bovino et al. 2016; Grassi et al. 2016; Richings & Schaye 2016) or radiative transfer effect (Petkova & Maio 2012; Roskar et al. 2014; Rosdahl et al. 2015;

Maio et al. 2016) or aim at removing tensions between different coding approaches (Kim et al. 2014).

Thus, we can improve galaxy simulations by providing theoretical expectations for [C II] that should be compared with state-of-the-art data. Such a synergy between theory and observations, in turn, can guide the interpretation of upcoming ALMA data and drive future experiments of large-scale [C II] mapping (Gong et al. 2012; Pallottini et al. 2015; Silva et al. 2015; Yue et al. 2015), which would lead to a statistical characterization of the high-z galaxy population. In this work, we simulate a $z \sim 6$ galaxy typically detected in [C II] with ALMA current observations.

The paper is structured as follows. In Section 2, we detail the numerical model used to set up the zoom-in simulation, and describe the adopted H$_2$ star formation prescription (Section 2.2), mass and energy inputs from the stellar populations (Section 2.3) and feedback (including SN, winds and radiation pressure, Section 2.4 – see also Appendices A and B). The results are discussed in Section 3, where we analyse star formation history and feedback effects in relation to ISM thermodynamics (Section 3.2) and its structural properties. The expected [C II] emission and other observational properties of high-z galaxies are discussed in Section 3.3. Conclusions are given in Section 4.

## 2 NUMERICAL SIMULATIONS

We carry out our simulation using a customized version of the adaptive mesh refinement (AMR) code RAMSES (Teyssier 2002). RAMSES is an octree-based code that uses a Particle Mesh N-body solver for the dark matter (DM) and an unsplit second-order MUSCL<fnref n="2" /> scheme for the baryons. Gravity is accounted by solving the Poisson equation on the AMR grid via a multigrid scheme with Dirichlet boundary conditions on arbitrary domains (Guillet & Teyssier 2011). For the present simulation, we choose a refinement based on a Lagrangian mass threshold-based criterion.

Chemistry and heating/cooling processes of the baryons are implemented with GRACKLE 2.1³ (Bryan et al. 2014), the standard library of thermochemical processes of the AGORA project (Kim et al. 2014). Via GRACKLE, we follow the H and He primordial network and tabulated metal cooling and photoheating rates calculated with CLOUDY (Ferland et al. 2013). Cooling also includes inverse Compton scattering off the CMB and heating from a redshift-dependent ionizing UV background (UVB; Haardt & Madau 2012). Since H$_2$ gas phase formation is not accounted for, we do not include the cooling contribution of such species.

Because of stellar feedback (Sections 2.3 and 2.4), the gas can acquire energy both in thermal and kinetic form. The distinction is considered by following the gas evolution of the standard thermal energy and a 'non-thermal' energy (Agertz et al. 2013). Such an approach is one of the possible schemes used to solve the overcooling problem that affects galaxy-scale simulations (see Dale 2015, and references therein). The non-thermal energy mimics turbulence, i.e. it is not affected by cooling. The non-thermal energy variation is due to gas advection ($v\nabla v$), work (PdV) and dissipation (Agertz & Kravtsov 2015). Following Mac Low (1999), we assume a dissipation time-scale proportional to the size of the cell (injection scale)

<footnote n="2">² MUSCL: Monotone Upstream-centred Scheme for Conservation Laws.</footnote>

<footnote n="3">³ See also https://grackle.readthedocs.org/.</footnote>
and inversely proportional to the Mach number.<fnref n="4" /> Since the dynamical time is essentially set by the free-fall time, the dissipation time can be written as $t_{\text{diss}} = 9.785(l_{\text{cell}}/100 \text{ pc})/(v_{\text{turb}}/10 \text{ km s}^{-1}) \text{ Myr}$. Then, the non-thermal energy loss due to dissipation can be written as $\dot{e}_{\text{nth}} = -e_{\text{nth}}/t_{\text{diss}}$ (Teyssier et al. 2013, see equation 2). As noted in Teyssier et al. (2013), such a scheme for non-thermal energy and its dissipation gives results qualitatively similar to a delayed cooling approach (Stinson et al. 2006).

## 2.1 Initial conditions

The initial conditions (IC) used for the suite are generated with MUSIC (Hahn & Abel 2011). MUSIC produces IC on a nested grid using a real-space convolution approach (cf. Bertschinger 1995). The adopted Lagrangian perturbation theory scheme is perfectly suited to produce IC for multimass simulations and – in particular – zoom simulations. To generate the ICs, the transfer functions are taken from Eisenstein & Hu (1998).

To set up the zoom-in simulation, we start by carrying out a cosmological DM-only run. The simulation evolves a volume $V^{\text{cosmo}} = (20 \text{ Mpc h}^{-1})^3$ from $z=100$ to 6 with DM mass resolution of $m_{\text{dm}}^{\text{cosmo}} = 3.4 \times 10^7 \text{ h}^{-1} \text{ M}_\odot$. The resolution of the coarse grid is $\Delta x^{\text{cosmo}} = 78.1 \text{ h}^{-1} \text{ kpc}$, and we do not include additional levels of refinement. Using HOP (Eisenstein & Hut 1998), we find the DM halo catalogue at $z=6$. The cumulative halo mass function extracted from the catalogue is in agreement with analytical expectations (e.g. Sheth & Tormen 1999), within the precision of halo-finder codes (e.g. Knebe et al. 2013).

From the catalogue, we select a halo with DM mass $M_h \simeq 10^{11} \text{ h}^{-1} \text{ M}_\odot$ (resolved by $\simeq 5 \times 10^4 \text{ DM particles}$), whose virial radius is $r_{\text{vir}} \simeq 15 \text{ kpc}$ at $z=6$. Using HOP, we select the minimum ellipsoid enveloping $10 r_{\text{vir}}$, and trace it back to $z=100$. As noted in Óñorbe et al. (2014), this is usually sufficient to avoid contamination.<fnref n="5" /> At $z=100$, the trace back volume is $V^{\text{zoom}} \simeq (2.1 \text{ Mpc h}^{-1})^3$. Using MUSIC we recalculate the ICs by generating three additional levels of refinement. For such a multimass set-up, the finer DM resolution is $m_{\text{dm}}^{\text{zoom}} = 6.7 \times 10^4 \text{ h}^{-1} \text{ M}_\odot$, which corresponds to a spatial resolution of $\Delta x^{\text{zoom}} = 9.7 \text{ h}^{-1} \text{ kpc}$. We note that because of the traced back volume, our simulation is expected to probe not only the target halo, but also its satellites and environment, similar to other works (e.g. Fiacconi, Feldmann & Mayer 2015, where the target halo is chosen at $z \simeq 3$).

In the zoom-in simulation, $\Delta x^{\text{zoom}}$ corresponds to our coarse grid resolution, and we allow for six additional refinement levels, based on a Lagrangian mass threshold-based criterion. At $z=6$, the baryonic component of the selected halo has a mass resolution of $m_b = 1.8 \times 10^4 \text{ M}_\odot$ and a physical resolution of $\Delta x^{\text{min}} = 31.9 \text{ pc}$. For convenience, a summary of the resolution outline can be found in Table 1. Note that the refined cell of our simulations has mass and size typical of molecular clouds (MC; e.g. Gorti & Hollenbach 2002; Federrath & Klessen 2013).

In this paper, we refer to metallicity (Z) as the sum of all the heavy element species without differentiating among them, and assume solar abundance ratios (Asplund et al. 2009). In the IC, the gas is characterized by a mean molecular weight $\mu = 0.59$ and has

<table><caption>**Table 1.** Resolution set-up for the cosmological run (COSMO) and subsequent zoom-in (ZOOM) simulation. $m_{\text{dm}}$ and $m_b$ are in units of $M_\odot \text{ h}^{-1}$ and indicate the dark matter (DM) and baryon mass resolution, respectively; $\Delta_x^{\text{max}}$ and $\Delta_x^{\text{min}}$ indicate the coarse grid and minimum available refinement scale, respectively. Both scales are reported in comoving kpc h$^{-1}$. For $\Delta_x^{\text{min}}$ we also report the physical pc scale at $z=6$. For the COSMO run, no refinement is used, and for the ZOOM, we indicate the increased resolution of the zoomed halo due to the multimass approach and the AMR.</caption><thead><tr><th></th><th><i>m</i><sub>dm</sub><br>(M<sub>☉</sub> h<sup>-1</sup>)</th><th><i>m</i><sub>b</sub></th><th>Δ<sub>x</sub><sup>max</sup><br>(kpc h<sup>-1</sup>)</th><th>Δ<sub>x</sub><sup>min</sup></th><th>Δ<sub>x</sub><sup>min</sup> at z = 6<br>(pc)</th></tr></thead><tbody><tr><td>COSMO</td><td>3.4 × 10<sup>7</sup></td><td>—</td><td>78.1</td><td>78.1</td><td>2.5 × 10<sup>3</sup></td></tr><tr><td>ZOOM</td><td>6.7 × 10<sup>4</sup></td><td>1.2 × 10<sup>4</sup></td><td>9.7</td><td>0.1</td><td>32.1</td></tr></tbody></table>

metallicity floor $Z = Z_{\text{floor}} > 0$. The metallicity floor mimics the pre-enrichment of the halo at high-$z$, when we do not have the resolution to follow precisely star formation and gas enrichment. We set $Z_{\text{floor}} = 10^{-3} Z_\odot$, a level that 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). Note that such low metallicity only marginally affects the gas cooling time but is above the critical metallicity for formation of Population III stars. Additionally, a posteriori, we have found that the metallicity floor contributes for only $\lesssim 0.2$ per cent of the total metal mass produced by stars by $z=6$ in the refined region.

## 2.2 Star formation model

We model SF by assuming an H$_2$-dependent Schmidt-Kennicutt relation (Schmidt 1959; Kennicutt 1998)

$$ \dot{\rho}_{*} = f_{\text{H}_2} \rho / t_{\text{sf}}, \qquad (1a) $$

where $\dot{\rho}_{*}$ is the local SF rate (SFR) density, $f_{\text{H}_2}$ the molecular hydrogen fraction, $\rho$ the gas density and $t_{\text{sf}}$ the SF time-scale. In equation (1a), we assume the SF time-scale to be proportional to the free-fall time, i.e.

$$ t_{\text{sf}} = \zeta_{\text{sf}}^{-1} \sqrt{3\pi/(32 G \rho)}, \qquad (1b) $$

where $\zeta_{\text{sf}}$ describes the SF efficiency and it is treated as a parameter in this work (cf. Semenov, Kravtsov & Gnedin 2016, see discussion in Section 3.2.1). To calculate $f_{\text{H}_2}$ we adopt the KTM09 model (Krumholz, McKee & Tumlinson 2008, 2009; McKee & Krumholz 2010). Such a model considers H$_2$ formation on dust grains by computing radiative transfer on an idealized MC and assumes equilibrium between the formation and dissociation rate of H$_2$. The solution for $f_{\text{H}_2}$ can be approximated as

$$ f_{\text{H}_2} = (1 - 0.75 s/(1 + 0.25 s)) \Theta(2 - s), \qquad (2a) $$

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

$$ \chi = 71 \left( \sigma_{\text{d,21}} / \mathcal{R}_{-16.5} \right) \left( \frac{G}{G_0} \right) \left( \frac{n}{\text{cm}^{-3}} \right), \qquad (2c) $$

where $\Theta$ is the Heaviside function, $\tau_{\text{uv}}$ is the dust optical depth of the cloud, $\sigma_{\text{d}}^{-21} = \sigma_{\text{d}}/10^{-21} \text{ cm}^{-2}$ is the dust absorption cross-section (Li & Draine 2001), $\mathcal{R}/10^{-16.5} \text{ cm}^3 \text{ s}^{-1}$ is the formation rate coefficient of H$_2$ on dust grains (Wolfire et al. 2008), G is the FUV flux in the Habing band (6–13.6 eV) normalized to the average Milky Way (MW) value $G_0$ (Habing 1968; Draine 1978) and n is the hydrogen number density. As in KTM09, we calculate the dust optical depth by linearly rescaling the MW value, i.e. $\tau_{\text{uv}} = 10^{-21} \text{ cm}^{-2} \text{ N}_H Z/Z_\odot/\mu$, where N$_H$ is the hydrogen column

<footnote n="4">⁴ While the distinction in thermal and non-thermal is similar to previous works (e.g. Agertz & Kravtsov 2015), we note that usually the time-scale for dissipation is fixed to 10 Myr.</footnote>

<footnote n="5">⁵ A posteriori, we have checked that the haloes in the zoom-in region have a contamination level $\lesssim 0.1$ per cent.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** H₂ fraction ($f_{H_2}$) as a function of gas density (n) obtained using the KTM09 model (equations 2). Different solid lines correspond to different metallicity (Z) of the gas. Horizontal dotted grey lines mark $f_{H_2}$ values of 0.5 and 1. Vertical dashed lines indicate the critical density $n_c$ where $f_{H_2} = 0.5$ for different Z; these critical density values are obtained as a fit (equation 3b) to the KTM09 model (equation 2). See the text for details.</figcaption></figure>

density and $\mu$ the mean molecular weight. In the simulation, the column density is calculated as $N_H = n l_{cell}$; because of the mass threshold-based criterion used as a refinement in AMR, we expect $l_{cell} \propto n^{-1/3}$; thus, $N_H \propto n^{2/3}$.

Note that both $\sigma_d$ and $R$ are proportional to the dust mass that we assume to be proportional to the metallicity. Then the ratio between $\sigma_d$ and $R$ is independent of Z. Additionally, equation (2) can be simplified by assuming pressure equilibrium between the CNM and WNM. In this case, equation (2c) turns out to be independent on $G/G_0$ and can be written as (Krumholz et al. 2009)

$$ \chi = 0.75 (1 + 3.1 (Z/Z_{\odot})^{0.365}). \quad (2d) $$

As shown in Krumholz & Gnedin (2011), for $Z \gtrsim 10^{-2} Z_{\odot}$ such approximation gives H₂ fractions compatible with those resulting from full non-equilibrium radiative transfer calculations.

In Fig. 1, we plot $f_{H_2}$ from the KTM09 model as a function of the gas density. Different solid lines refer to different metallicity. At a fixed metallicity, the molecular fraction as a function of density vanishes for low values of n; it steeply rises up to $f_{H_2} \sim 0.8$ in one density dex and asymptotically reaches $f_{H_2} = 1$. The critical density where the gas can be considered molecular ($f_{H_2} = 0.5$) is roughly inversely proportional to the metallicity, i.e. $n_c \sim 25(Z/Z_{\odot})^{-1} \text{ cm}^{-3}$ (see also Agertz et al. 2013). We note that when detailed chemistry calculations are performed, such critical density depends on the chemical network and the assumptions regarding gas shielding from external radiation and clumpiness. As a consequence, the actual critical density can be higher than the one predicted by the KTM09 model (e.g. Bovino et al. 2016).

Because of the particular shape of the $f_{H_2}(n)$ relation, the adopted SF law (equations 1–2) is roughly equivalent to a prescription based on a density threshold criterion:

$$ \dot{\rho}_{*} = \Theta(n - n_c)m_p n/t_{sf}, \quad (3a) $$

where $m_p$ is the proton mass and the critical density

$$ n_c \simeq 26.45 (Z/Z_{\odot})^{-0.87} \text{ cm}^{-3} \quad (3b) $$

is calculated as a fit to the $f_{H_2}$ KTM09 model. In Fig. 1, we show $n_c$ for various metallicities (dashed vertical lines).

Equations (3) are not used to calculate the SFR in the simulation. However, being simpler, such formulation can be used to enhance our physical intuition of the adopted SF law<fnref n="6" /> in analysing the results. As noted in Hopkins, Narayanan & Murray (2013), the morphology of a galaxy is very sensitive to the minimum density of the cells that are able to form stars.

During the simulation, Equations (1) are 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). We impose that no more than half of a cell mass can be converted into a star particle in each event. This prescription ensures the numerical stability of the code (Dubois & Teyssier 2008). This is also consistent with the picture that nearly half of the mass in an MC is Jeans unstable (Federrath & Klessen 2013).

We allow SF only if the mass of a new star particle is at least equal to the baryon mass resolution. This avoids numerical errors for the star particle dynamics and enables us to treat the particle as a stellar population with a well-sampled initial mass function (IMF). Additionally, the SF law is driven by H₂ formation on dust grains; we do not allow gas to form stars if the dust temperature is larger than $\simeq 2 \times 10^3$, because of dust sublimation (see Section 2.4 and Appendix B for the details on the dust prescriptions).

For this work, we assume an SF efficiency $\zeta_{sf} = 10$ per cent, in accordance with the average values inferred from MC observations (Murray 2011, see also Agertz et al. 2013). Note that varying the parameters for the SF law should lead to a similar SFR once feedback is properly included, although the galaxy morphology can be different (Hopkins et al. 2013).

## 2.3 Mass and energy inputs from stars

Because of the finite mass resolution, it is necessary to introduce (according to equations 1–2d) ‘star particles’ to represent stellar populations. To this aim, we adopt a Kroupa (2001) IMF:

$$ \Phi(m) \propto \left[ m^{-\alpha_1} \Theta(m_1 - m) + m^{-\alpha_2} \Theta(m - m_1) m_1^{\alpha_2 - \alpha_1} \right], \quad (4a) $$

where $\alpha_1 = 1.3$, $\alpha_2 = 2.3$, $m_1 = 0.5 M_{\odot}$ and m is in the range $[10^{-1}-10^2] M_{\odot}$. The proportionality constant is chosen such that

$$ \int_{0.1 M_{\odot}}^{100 M_{\odot}} m \Phi dm = 1. \quad (4b) $$

Once formed, stars affect the environment with chemical, mechanical and radiative feedback. These stellar inputs are parametrized by the cumulative fraction of the returned gas mass, metals and energy (e.g. Salvadori, Ferrara & Schneider 2008; de Bennassuti et al. 2014; Salvadori, Skúladóttir & Tolstoy 2015). Mass and energy inputs are conveniently expressed per unit stellar mass formed ($M_*$).

Chemical feedback depends on the return fraction ($R$) and the yield ($Y$):

$$ R(t_{*}) = \int_{m(t_{*})}^{100 M_{\odot}} (m - w) \Phi dm, \quad (5a) $$

<footnote n="6">⁶ As a consequence of the rough equivalence, it is not necessary to manually prevent SF in underdense regions, by imposing that an overdensity $\Delta > 200$ is needed to form stars. At the start of the simulation ($z = 100$), the mean density of the gas is $\sim 0.1 m_p \text{ cm}^{-3}$, while the ‘effective’ SF threshold would be $n_c \sim 10^4 \text{ cm}^{-3}$ for the gas at $Z = Z_{\text{floor}}$.</footnote>
$$Y(t_*) = \int_{m(t_*)}^{100 M_{\odot}} m_Z \Phi \, dm, \qquad (5b)$$

where $w(m, Z_*)$ and $m_Z(m, Z_*)$ are the stellar remnant and the metal mass produced for a star of mass $m$ and metallicity $Z_*$ (e.g. Woosley & Weaver 1995; van den Hoek & Groenewegen 1997), and $m(t_*)$ is the minimum stellar mass with lifetime<fnref n="7" /> shorter than $t_*$, the time elapsed from the creation of the stellar particle (i.e. the ‘burst age’).

This approach is used both in zoom galaxy simulations (e.g. Kim et al. 2014) and cosmological simulations (e.g. Pallottini et al. 2014a, hereafter P14). Compared to cosmological simulations, though, zoom simulations have typically a better spatial and – consequently – time resolution (e.g. $\Delta t \sim 10^{-2}$ Myr versus $\Delta t \sim $ Myr). Thus, here we can follow the gradual release of both gas and metals in the ISM.

The mechanical energy input includes SN explosions and winds, either by OB or asymptotic giant branch (AGB) stars in young (<40 Myr) or evolved stellar populations:

$$\epsilon_{\text{sn}}(t_*) = \int_{m(t_*) > 8 M_{\odot}}^{40 M_{\odot}} e_{\text{sn}} \Phi \, dm, \qquad (5c)$$

$$\epsilon_w(t_*) = \int_{m(t_*)}^{100 M_{\odot}} e_w \Phi \, dm, \qquad (5d)$$

where $e_{\text{sn}} = e_{\text{sn}}(m, Z)$ and $e_w = e_w(m, Z)$ are the energy released by SN and stellar winds in units of $10^{51}$ erg $\equiv 1$ foe; we have further assumed that only stars with $8 \le m/M_{\odot} \le 40$ can explode as SN.

Radiative energy inputs can be treated within a similar formalism. The cumulative energy $\epsilon_{12}$ associated with the spectral range ($\lambda_1, \lambda_2$) can be written as

$$\epsilon_{12}(t_*) = \int_{0}^{t_*} \int_{m(t)}^{100 M_{\odot}} L_{12} \Phi \, dm \, dt, \qquad (5e)$$

$$L_{12}(t) = \int_{\lambda_1}^{\lambda_2} L_{\lambda} \, d\lambda, \qquad (5f)$$

where $L_{\lambda} = L_{\lambda}(m, Z_*)$ is the luminosity per unit wavelength and mass. For convenience, we express the radiation energy in units of foe, as for the mechanical energy (equation 5c). In the following, we specify $\epsilon_{12}$ in equation (5e), by separately considering ionizing radiation ($\lambda_1 = 0, \lambda_2 = 912$ Å) denoted by $\epsilon_{\text{ion}}$, and the soft UV band, $\epsilon_{\text{uv}}$, defined as the range ($\lambda_1 = 912$ Å, $\lambda_2 = 4000$ Å).

In equations (5), the quantities $w, m_Z, e_{\text{sn}}, e_w$ and $L_{\lambda}$ can be calculated from stellar evolutionary models. We adopt the PADOVA (Bertelli et al. 1994) stellar tracks for metallicities $Z_*/Z_{\odot} = 0.02, 0.2, 0.4$ and $1$ to compute the chemical,<fnref n="8" /> mechanical and radiative inputs using STARBURST99 (Leitherer et al. 1999, 2010).

In Fig. 2, we plot $R, Y, \epsilon_{\text{sn}}, \epsilon_w, \epsilon_{\text{ion}}$ and $\epsilon_{\text{uv}}$ as a function of $t_*$. For each curve, the shaded regions denote the $0.02 \le Z_*/Z_{\odot} \le 1$ metallicity range; single $Z_*$ tracks are indicated with dark lines. The time interval during which massive stars can explode as SN ($0.8 \lesssim \log t_*/\text{Myr} \lesssim 1.6$) is highlighted with vertical dashed lines, and the upper axis is labelled with the corresponding stellar mass.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Stellar inputs (cumulative fraction) as a function of stellar age ($t_*$). Shown are the return fraction ($R$), metal yield ($Y$), SN mechanical energy ($\epsilon_{\text{sn}}$), wind mechanical energy ($\epsilon_w$), ionizing radiation energy ($\epsilon_{\text{ion}}$) and UV radiation energy ($\epsilon_{\text{uv}}$). The fractions are given per unit stellar mass formed; energies are expressed in units of $10^{51}$ erg $\equiv$ foe. Cumulative fractions are indicated with different colours, as indicated in the legend: the shaded regions cover the $0.02 \le Z_*/Z_{\odot} \le 1$ metallicity range; dark lines denote single metallicity PADOVA stellar tracks (Bertelli et al. 1994). To guide the eye, the SN explosion period is bracketed by vertical dashed lines; in the upper axis we report the value of $m(t_*)$, the minimum stellar mass corresponding to the stellar lifetime $t_*$. For definitions, see Equations (5).</figcaption></figure>

Note that the OB stars’ contribution ($\log t_*/\text{Myr} \lesssim 0.8$) to $\epsilon_w, Y$ and $R$ is roughly proportional to $t_*$ and $Z_*$ (see also Agertz et al. 2013, in particular equations 4). As in the simulation, the metallicity floor is set to $Z_{\text{floor}} = 10^{-3} Z_{\odot}$, we slightly overestimate the wind contribution for low $Z_*$.

Finally, note that the change of behaviour of $\epsilon_{\text{uv}}$ at $\log t_*/\text{Myr} \lesssim 2$ is due to the ionizing ($\lambda \le 912$ Å) photon production suppression. At late times ($\log t_*/\text{Myr} \gtrsim 1.6$), AGB stars give a negligible mechanical energy contribution ($\epsilon_w \simeq \text{constant}$) but return mass and metals to the gas ($R, Y$).

## 2.4 Stellar feedback

Equations (5) provide us with the energy produced by stars in different forms. The next step is to understand what fraction of that energy is eventually deposited in the ISM. Consider a stellar population of initial mass $M_*$, metallicity $Z_*$ and age $t_*$ residing in a gas cell with volume $V_{\text{cell}}$. In our scheme, when the simulation evolves for a time $\Delta t$, the chemical feedback act as follows:

$$\rho = \rho + [R(t_* + \Delta t) - R(t_*)] M_*/V_{\text{cell}}, \qquad (6a)$$

$$Z = Z + [Y(t_* + \Delta t) - Y(t_*)] M_*/V_{\text{cell}}, \qquad (6b)$$

where $\rho$ and $Z$ are the gas density and metallicity, and $R$ and $Y$ are taken from equation (5a). Note that chemical enrichment is due to both the SN and AGB winds.

### 2.4.1 Supernova explosions

For the mechanical feedback, let us first consider the case of SNe. At each SN event, the specific energy of the gas changes as

$$e_{\text{th}} = e_{\text{th}} + f_{\text{th}} [\epsilon_{\text{sn}}(t_* + \Delta t) - \epsilon_{\text{sn}}(t_*)] M_*/V_{\text{cell}}, \qquad (6c)$$

<footnote n="7">⁷ Stellar lifetimes are roughly independent of metallicity for $Z_* > 10^{-4} Z_{\odot}$ (Raiteri, Villata & Navarro 1996, see equation 3).</footnote>

<footnote n="8">⁸ Similarly to Kim et al. (2014), when computing the yields in equation (5a), we assume that the metal mass is linked to the oxygen and iron masses via $m_Z = 2.09 m_O + 1.06 m_{\text{Fe}}$, as appropriate for Asplund et al. (2009) abundances.</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 3. Example of the adopted feedback model. Fractional energy evolution for a single SN explosion ($E_0 = 1$ foe) in a gas characterized by $n = 1 \text{ cm}^{-3}$ and $Z = 10^{-3} Z_{\odot}$ as a function of the time interval from the explosion $\Delta t$. We plot the total ($f = f_{th} + f_{kn}$), thermal ($f_{th}$) and kinetic ($f_{kn}$) energy fraction acquired by the gas (see equation 6c) with solid black, dashed red and dotted blue lines, respectively. Shaded regions indicate different stages of the SN evolution, i.e. the energy conserving Sedov–Taylor (ST) stage, shell formation (SF) stage and pressure-driven snowplow (PDS) stage. In the adopted formalism, the initial energy $E_0$ is a function of the stellar input (see Section 2.3 and equation 5c), e.g. $E_0 = [\epsilon_{sn}(t_* + \Delta t) - \epsilon_{sn}(t_*)]M_*$. The full model is presented in Appendix A.</figcaption></figure>

$$e_{\text{nth}} = e_{\text{nth}} + f_{\text{kn}} [\epsilon_{\text{sn}}(t_{*} + \Delta t) - \epsilon_{\text{sn}}(t_{*})] M_{*} / V_{\text{cell}}, \quad (6d)$$

where $e_{\text{th}}$ and $e_{\text{nth}}$ are the thermal and non-thermal energy densities, and $f_{\text{th}}$ and $f_{\text{kn}}$ are the fractions of thermal and kinetic energy deposited in the ISM. Thus, $e_{\text{nth}}$ accounts for the momentum injection by SN and $e_{\text{th}}$ for the thermal pressure part.

In this work, we have developed a novel method to compute such quantities. The method derives $f_{\text{th}}$ and $f_{\text{kn}}$ from a detailed modelling of the subgrid blast wave evolution produced by the SN explosion. We calculate $f_{\text{th}}$ and $f_{\text{kn}}$ by evaluating the shock evolution at time $\Delta t$, the time step of the simulation.<fnref n="9" />

The adopted blast wave model is based on Ostriker & McKee (1988, hereafter OM88), and it accounts for the evolution of the blast through its different evolutionary stages (energy conserving, momentum conserving, etc.). While each stage is self-similar, the passage from one stage to the next is determined by the cooling time. Thus, $f_{\text{th}}$ and $f_{\text{kn}}$ depend on the blast wave evolutionary stage. The latter, in turn, depends on the gas density, cooling time and initial energy of the blast ($E_0 = [\epsilon_{\text{sn}}(t_* + \Delta t) - \epsilon_{\text{sn}}(t_*)]M_*$, in equation 5c).

The model details are presented in Appendix A. As an example, in Fig. 3 we show the energy evolution for a single SN explosion ($E_0 = 1$ foe) in a gas characterized by $n = 1 \text{ cm}^{-3}$ and $Z = 10^{-3} Z_{\odot}$. The total energy $E(t)$ is constant in the Sedov–Taylor (ST) stage, it decreases to $0.5 E_0$ during the shell formation (SF) stage and it evolves as $\Delta t^{-2/7}$ in the pressure-driven snowplow (PDS) stage (see equation A1). In the ST stage, most of the energy is thermal, i.e. $f_{\text{kn}}/f_{\text{th}} \simeq 0.4$; however, in the SF stage $f_{\text{kn}}$ increases, since part of the thermal energy is radiated away and some is converted into kinetic

form (e.g. Cox 1972; Cioffi, McKee & Bertschinger 1988). Finally, during the PDS stage the ratio of thermal to kinetic energy is $f_{\text{kn}}/f_{\text{th}} \simeq 2$ (see equation 6.14 in OM88).

In this particular example – a 1 foe SN exploding in an $n = 1 \text{ cm}^{-3}$ cell – by assuming a simulation time step of $\Delta t \simeq 10^{-2}$ Myr, we find that the blast wave is in the PDS stage, and the gas receives (via equation 6c) a fraction of energy $f_{\text{th}} \simeq 8$ per cent and $f_{\text{kn}} \simeq 16$ per cent in thermal and kinetic form, respectively. During $\Delta t$, about $\simeq 75$ per cent of the initial SN energy has been either radiated away or lost to work done by the blast wave to exit the cell. The model is in broad agreement with other more specific numerical studies (e.g. Cioffi, McKee & Bertschinger 1988; Martizzi, Faucher-Giguère & Quataert 2015; Walch & Naab 2015).

### 2.4.2 Stellar winds

Stellar winds are implemented in a manner paralleling the above scheme for SNe. The energy variation can be calculated via equation (5c), where $\epsilon_{\text{sn}}$ is substituted with $\epsilon_w$, given in equation (6c). Then, $f_{\text{th}}$ and $f_{\text{kn}}$ for winds are calculated via a stage scheme similar to SN. The main difference in the efficiency factor calculation depends on the mode of energy production, i.e. impulsive for SNe, continuous for winds. The complete scheme is detailed in Appendix A.

The efficiency of SN is greatly increased when the gas is pre-processed by stellar winds (Walch & Naab 2015; Fierlinger et al. 2016), since the energy loss process is highly non-linear (Fierlinger et al. 2016, see fig. 8). For example, when an SN explodes in the lower density bubble produced by the stellar progenitor wind, the adiabatic phase lasts longer and consequently $f_{\text{kn}}$ and $f_{\text{th}}$ increase considerably.

### 2.4.3 Radiation pressure

Finally, we account for radiation pressure from stars. The coupling of the gas with the radiation can be expressed in terms of $\dot{p}_{\text{rad}}$, the rate of momentum injection (Krumholz & Matzner 2009; Hopkins, Quataert & Murray 2011; Krumholz & Thompson 2012; Wise et al. 2012; Agertz et al. 2013), and accounts for the contribution from ionization, and from dust UV heating and IR trapping

$$\dot{p}_{\text{rad}} = (L_{\text{ion}}/c)(1 - \exp(-\tau_{\text{ion}})) + (L_{\text{uv}}/c)((1 - \exp(-\tau_{\text{uv}})) + f_{\text{ir}}), \quad (7a)$$

where $c$ is the speed of light, $\tau_{\text{ion}}$ is the hydrogen optical depth to ionizing radiation and $f_{\text{ir}}$ is the term accounting for the IR trapping. $L_{\text{ion}}$ and $L_{\text{uv}}$ are calculated by integration of the stellar tracks (equation 5e). The calculation of $\tau_{\text{uv}}$ is modelled in Section 2.2 (equation 2b and related text). We compute $\tau_{\text{ion}}$ and $f_{\text{ir}}$ according to the physical properties of the gas, as detailed in Appendix B. Note that we do not assume, as sometimes done, $\tau_{\text{ion}} \sim \tau_{\text{uv}} \gg 1$, i.e. we allow for the possibility that some LyC photons can escape.

In smoothed particle hydrodynamics (SPH) codes, radiation pressure (equation 7a) can be implemented as a 'kick' (e.g. Hopkins et al. 2011; Barai et al. 2015). Namely, a velocity $\Delta v = \dot{p}_{\text{rad}}\Delta t/m_b$ is directly added to some of the SPH particles of mass $m_b$ near the photon source. The particles that receive kicks are statistically chosen according to a probability $P_{\text{kick}}$ and with kick direction $\hat{v}$ that is sampled from a random distribution. Considering the specific kinetic energy of the SPH particles, we would have

$$e_{\text{k}} = 0.5 \langle m_{\text{b}}(\mathbf{v} + \Delta v P_{\text{kick}} \hat{\mathbf{v}})^2 \rangle / V_{\text{cell}}, \quad (7b)$$

<footnote n="9">⁹ The underlying assumption is that the shock fronts exit the cell in $\lesssim \Delta t$. This is quite consistent because the shock is expected to be supersonic, and the sound crossing time is larger or comparable with the simulation time step $\Delta t$, dictated by the Courant–Friedrichs–Lewy conditions.</footnote>
where $v$ is the original particle velocity, the $\langle \rangle$ operator indicates the particles sum weighted by the SPH kernel and $V_{cell}$ is the kernel volume. Thus, because of the kick, the increase of energy density would be<fnref n="10" />

$$ \begin{aligned} \Delta e_k &= \langle m_b P_{\text{kick}}(\Delta v \mathbf{v} \hat{\mathbf{v}}) + 0.5(\Delta v)^2 \rangle / V_{\text{cell}} \\ &= 0.5 m_b (\Delta v)^2 / V_{\text{cell}} \\ &= 0.5 (\dot{p}_{\text{rad}} \Delta t)^2 / (m_b V_{\text{cell}}), \end{aligned} \quad (7c) $$

where $\dot{p}_{\text{rad}}$ can be calculated via equation (7a), and equation (7c) can be directly cast into the AMR formalism. Additionally, because of our approximate treatment of IR trapping (see Appendix B), we force energy conservation: $V_{cell}\Delta e_k \le \Delta t(L_{\text{ion}} + L_{\text{uv}})$, i.e. the deposited energy must not exceed the radiative input energy. Finally, we recall here that non-thermal energy is dissipated with a time-scale $t_{\text{diss}}$, as described in the beginning of Section 2.

# 3 RESULTS

At $z = 6$ ($t \simeq 920$ Myr), the simulated zoom-in region contains a group of 15 DM haloes that host galaxies. We target the most massive halo ($M_h = 1.8 \times 10^{11} M_{\odot}$) that hosts 'Dahlia', which is a galaxy characterized by a stellar mass of $M_* = 1.6 \times 10^{10} M_{\odot}$, therefore representative of a typical Lyman-break galaxy (LBG) at that epoch. Dahlia has 14 satellites located within $\simeq 100$ kpc from its centre. The six largest ones have a DM mass in the range $M_h = 2.5 \times 10^9 - 1.2 \times 10^{10} M_{\odot}$, and they host stars with total mass $M_* \lesssim 10^9 M_{\odot}$. Additionally, there are eight smaller satellites ($M_h \simeq 10^7 M_{\odot}$), with $M_* \simeq 10^5 M_{\odot}$.

## 3.1 Overview

We start by looking at the overall properties of Dahlia on decreasing scales. In the following, we refer to Fig. 4, which shows the simulated density ($n$), temperature ($T$), total (thermal+kinetic) pressure ($P$) and metallicity ($Z$) maps<fnref n="11" /> at $z = 6$.

### 3.1.1 Environment (scale $\simeq 160$ kpc)

Dahlia sits at the centre of a cosmic web knot and accretes mass from the intergalactic medium (IGM) mainly via three filaments of length $\simeq 100$ kpc, confirming previous findings (Dekel et al. 2009). These overdense filaments ($n \simeq 10^{-2} \text{ cm}^{-3}$) are slightly colder ($T \simeq 10^{3.5}$ K) than the IGM ($\langle T \rangle \simeq 10^{4.5}$ K) as a consequence of their shorter radiative cooling time ($t_{\text{cool}} \propto n^{-1}$). Along these cold streams, pockets of shock-heated ($T \gtrsim 10^{4.5}$ K) gas produced by both structure formation and feedback (SN and winds) are visible.

The galaxy locations can be pinpointed from the metallicity map, showing a dozen metal-polluted regions. The size of the metal bubbles ranges from $\simeq 20$ kpc in the case of Dahlia to a few kpc for the satellites. Bubble sizes increase with the total stellar mass (see P14, in particular fig. 13) and the age of the galaxy stellar population.

On these scales, the pressure is dominated by the thermal component ($P \simeq P_{\text{th}} \sim 10^4 \text{ K cm}^{-3}$); higher values of pressure, associated

with non-thermal feedback effects (e.g. gas bulk motion), are confined around star-forming regions, again traced by the metallicity distribution.

### 3.1.2 Circumgalactic medium (scale $\simeq 50$ kpc)

To investigate the circumgalactic medium (CGM), we zoom in a region within $\sim 3 r_{\text{vir}} = 47.5$ kpc from Dahlia's centre. On these scales, we can appreciate the presence of several Dahlia's satellites, i.e. extended (few kpc) structures that are $\sim 100$ times denser than the filament in which they reside. Two of these density structures are particularly noticeable. These are located at a distance of $\sim 10$ kpc from the centre in the upper left and lower left part of the map, respectively. By looking at the metallicity distribution, we find that both satellites reside within their own metal bubble, which is separated from Dahlia's one. This clearly indicates an in situ star formation activity.

Additionally, the density map shows about 20 smaller ($\sim 10-100$ pc) overdense clumps ($n \gtrsim 10 \text{ cm}^{-3}$). The ones within Dahlia's metal bubble are enriched to $Z \simeq Z_{\odot}$. This high $Z$ value is indicative of in situ self-pollution, which possibly follows an initial pre-enrichment phase from Dahlia. Clumps outside Dahlia's metal bubble have on average a higher density ($n \sim 10^2 \text{ cm}^{-3}$). Since these clumps are unpolluted, they have not yet formed stars, as the effective density threshold for star formation is $\sim 25/(Z/Z_{\odot}) \text{ cm}^{-3}$ (see equation 3b and Section 2.2). Such clumps represent molecular cloud complexes caught in the act of condensing as the gas streams through the CGM (Ceverino et al. 2016). Such clumps have gas mass in the range $10^5-10^6 M_{\odot}$ and are not DM-confined, as the DM density field is flat on their location.

Star-forming regions are surrounded by an envelope of hot ($T \simeq 10^{5.5}$ K), diffuse ($n \gtrsim 10^{-2} \text{ cm}^{-3}$) and mildly enriched ($z \sim 10^{-2} Z_{\odot}$) gas produced by SN explosions and winds. In the centre of star-forming regions, instead, the gas can cool very rapidly due to the high densities/metallicities. Nevertheless, these regions are highly pressurized due to bulk motions mostly driven by radiation pressure (see Fig. 8).

### 3.1.3 ISM (scale $\simeq 10$ kpc)

The structure of Dahlia's ISM emerges once we zoom in a region $\sim 0.5 r_{\text{vir}}$ from its centre. In the inner region ($\simeq 2$ kpc), a counterclockwise disc spiral pattern is visible, since the field of view is perpendicular to the rotation plane of the galaxy (see Gallerani et al. 2016 for the analysis of the velocity field of Dahlia). The presence of discs in these early systems has already been suggested by other studies. For example, Feng et al. (2015) show that already at $z \sim 8$ nearly 70 per cent of galaxies with $M_* \simeq 10^{10} M_{\odot}$ have discs (see also Section 3.3).

The spiral central region and the spiral arms are dense ($n \simeq 10^2 \text{ cm}^{-3}$) and cold ($T \simeq 10^3$ K), and the active SF produces a large in situ enrichment ($Z \simeq Z_{\odot}$). Winds and shocks from SN have no effect in the inner part of the galaxy, because of the high density and short cooling time of the gas; this implies that metals remain confined within $\sim 2$ kpc. Within spiral arms, radiation-pressure-induced bulk motions largely dominate the total pressure, which reaches values as high as $P \gtrsim 10^{6.5}$ K cm$^{-3}$. The imprint of SN shocks is evident in the temperature map in regions with $T \gtrsim 10^5$ K. Shock-driven outflows originated in spiral arms travel outward in the CGM, eventually reaching the IGM if outflow velocities exceed the escape velocity ($\sim 100$ km s$^{-1}$, see fig. 4 in Gallerani et al. 2016).

<footnote n="10">¹⁰ In equation (7c), when going from the first to the second line, note that first terms give a null contribution, as $v$ is ordered motion, while the kicks are randomly oriented via $\hat{v}$, and that, by definition, $\langle P_{\text{kick}} \rangle = 1$.</footnote>

<footnote n="11">¹¹ Most of the maps of this paper are obtained with a customized version of PYMSES (Labadens et al. 2012), a visualization software that implements optimized techniques for the AMR grid of RAMSES.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Maps of the simulated galaxy Dahlia $z=6$. From left to right, we plot subsequent zooms on the galaxy. From top to bottom, we plot the density ($n$), temperature ($T$), pressure ($P$) and metallicity ($Z$). Each map is obtained<sup>11</sup> by mass averaging the physical quantity along the line of sight piercing the field of view and centred on Dahlia. In all panels, the physical scale is indicated as an inset. Movies of Dahlia can be found at https://www.researchgate.net/profile/Andrea_Pallottini.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Density (n, blue), metallicity (Z, red) and molecular hydrogen density ($n_{H_2}$, yellow) radial profile (r) with respect to Dahlia centre. The profiles are spherically averaged, as indicated by the (>)$_V$ operator, and the upper axis shows the radial distance r as a function of the virial radius of Dahlia ($r_{vir}$).</figcaption></figure>

Outflows are either preferentially aligned with the galaxy rotation axis or they start at the edge of the disc. However, when spherically averaged, infall and outflow rates are nearly equal (~30 M$_{\odot}$ yr$^{-1}$ at z ~ 6; Gallerani et al. 2016), and the system seems to self-regulate (see also Dekel & Mandelker 2014).

Outside the disc, clumps with density $n \simeq 10^2$ cm$^{-3}$ are also present and are actively producing stars. These isolated star-forming MCs are located at a distance $\gtrsim 2$ kpc from the centre, and show up as spots of high pressure ($P \gtrsim 10^7$ K cm$^{-3}$); some of these MCs are completely disrupted by internal feedback, and they can be recognized by the low metallicity ($z \sim 10^{-3}$ Z$_{\odot}$): this is consistent with the outcome of numerical simulations of multiple SN explosions in a single MC (e.g. Körtgen et al. 2016).

### 3.1.4 Radial profiles

Fig. 5 shows spherically averaged density, metallicity and $H_2$ density profiles for the gas. The density profile rapidly decreases from $n \sim 30$ cm$^{-3}$ at $r \sim 0$ to $n \sim 0.1$ cm$^{-3}$ at $r \sim 6$ kpc (~0.5 $r_{vir}$) and then flattens at larger distances. Such a profile is consistent with the average profile of $z=4$ galaxies presented in P14. There we claimed that the density profile is universal once rescaled to the halo virial radius (see also Liang, Kravtsov & Agertz 2016). Superposed to the mean density profile, local peaks are clearly visible: they result from individual clumps/satellites, as discussed above.

The central metallicity is close to the solar value, but by $r \sim 12$ kpc ~ $r_{vir}$ it has already dropped to $Z = Z_{\text{floor}}$. Within $0 \lesssim r/kpc \lesssim 6$, the metallicity gradient closely tracks the density profile, while for $6 \lesssim r/kpc \lesssim 15$ the decrease is steeper. Pallottini, Gallerani & Ferrara (2014b) find that the metallicity profile is not universal; however, it usually extends up to a few virial radii, as for Dahlia. Further insights can be obtained by analysing the n-Z relation (Section 3.2.3).

In Fig. 5, we note that the Z gradient found in Dahlia at $z=6$ is slightly steeper than the one inferred from observations of $z \sim 3$ galaxies, i.e. we find $\Delta Z/r \sim -0.1$ dex kpc$^{-1}$ while the observed ones are ~0 dex kpc$^{-1}$ (Wuys et al. 2016) and ~+0.1 dex kpc$^{-1}$ (Troncoso et al. 2014). This suggests that the metallicity profile evolves with cosmic time and that the flattening is likely caused by stellar feedback, which in our Dahlia may occur in the following

Gyr of the evolution. However, to prove such a claim we should evolve the simulation to $z \sim 3$.

The $H_2$ profile is spiky, and each peak marks the presence of a distinct SF region.<fnref n="12" /> In Dahlia, $H_2$ is mainly concentrated within $r \lesssim 0.5$ kpc, and it is distributed in the disc-like structure seen in Fig. 4 (see Section 3.3). The location of the other peaks corresponds to the satellites, which are mostly colocated with metallicity peaks. With increasing metallicity, in fact, lower densities are needed to form $H_2$ (equation 3b).

## 3.2 Star formation and feedback history

We analyse the SF history of Dahlia and its major satellites by plotting in Fig. 6 the cumulative stellar mass ($M_*$) and star formation rate (SFR) versus time.<fnref n="13" />

For the whole galaxy sample, the time-averaged (±rms) specific star formation is $\langle \text{sSFR} \rangle = (16.6 \pm 32.8) \text{ Gyr}^{-1}$. This mean value is comparable to that obtained by previous simulations of high-z galaxies (Wise et al. 2012) and broadly in agreement with $z \sim 7$ observations (Stark et al. 2013). At early times, the sSFR reaches a maximum of ~100 Gyr$^{-1}$, while a minimum of 3.0 Gyr$^{-1}$ is found during the late-time evolution. Both the large SFR range and the maximum at early times are consistent with simulations by Shen et al. (2014). At late times, the SFR is in agreement with analytical calculation (Behroozi, Wechsler & Conroy 2013) and with $z=7$ observations (González et al. 2010), although we note Dahlia has a larger stellar mass with respect to the galaxies in the sample ($M_* \simeq 5 \times 10^9 M_{\odot}$).

At all times, Dahlia dominates both the stellar mass and the SFR, whose mean value is $\langle \text{SFR} \rangle \simeq (35.3 \pm 32.7) M_{\odot} \text{ yr}^{-1}$. Its stellar mass grows rapidly, and it reaches $M_* \sim 10^9 M_{\odot}$ by $t \simeq 400$ Myr ($z=11$), i.e. after $\simeq 120$ Myr from the first star formation event. Such rapid mass build-up is due to merger-induced SF, which plays a major role at high-z (Salvadori, Dayal & Ferrara 2010; Behroozi, Wechsler & Conroy 2013; Poole et al. 2016). The SFR is roughly constant from $z \sim 11$ to ~8.5 and reaches a maximum of $\simeq 130 M_{\odot} \text{ yr}^{-1}$ at $z \sim 6$. With respect to observations of $z \sim 6$ LBG galaxies (e.g. Stanway, Bunker & McMahon 2003; Stark et al. 2009), the SFR and $M_*$ of Dahlia are above the mean values but still consistent within $1\sigma$. Additionally, the combination of SFR, $M_*$ and $Z_*$ for Dahlia is compatible with the fundamental mass-metallicity relation observed in local galaxies (Mannucci et al. 2010).

The total stellar mass in satellites is $M_* \sim 10^9 M_{\odot}$. Typically, SF starts with a burst, generating $\sim 10^{7.5} M_{\odot}$ of stars during the first $\simeq 20$ Myr. Then, the SFR exponentially declines and becomes intermittent with a bursty duty cycle of ~100 Myr. This process can be explained as follows. As a halo forms, at its centre the density of the gas slowly rises. When the density is higher than the critical density of $H_2$ formation (equation 3b), the gas in the inner region is converted into stars in a few free-fall times. Then feedback, and in particular coherent SN explosions ($t_* \gtrsim 10$ Myr, see Fig. 2), quenches the SFR, and the star formation activity becomes self-regulated. As mergers supply fresh gas, the SFR suddenly goes out of equilibrium and becomes bursty again. Note that self-regulation is possible only for major satellites, since smaller ones ($M_h \lesssim 10^8 M_{\odot}$)

<footnote n="12">¹² We remind that the profiles are volume-weighted; thus, the plotted $n_{H_2}$ accounts for the fact that $H_2$ is present only in a fraction of the gas at a given radius.</footnote>

<footnote n="13">¹³ The SFR is averaged in steps of $\simeq 3$ Myr. We have checked that smaller steps do not alter the following analysis.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Left: Cumulative stellar masses ($M_*$, upper left-hand panel) and star formation rates (SFR, lower left-hand panel) of Dahlia and its satellites as a function of cosmic time ($t$). For each galaxy, individual $M_*$ and SFR are plotted with a solid line, coloured accordingly to the total dark matter mass ($M_h$) of the host halo at $z=6$. For both $M_*$ and SFR, Dahlia's tracks are plotted with a black solid line, and the totals (Dahlia+satellites) are in dash-dot black lines. Right: SFR as a function of cosmic time, with individual galaxies defined by the merger history up to $z \simeq 8.5$. Note the different $M_h$ colour bar scale with respect to the left-hand panel.</figcaption></figure>

cannot retain a large fraction of their gas following feedback events due to their shallow potential wells (see P14).

Note that the duty cycle and the amplitude of the burst are fairly in agreement with observations of $M_* \sim 10^8-10^{10} M_{\odot}$ galaxies at $z \lesssim 0.3$ (Kauffmann 2014). Furthermore, in our satellites we find that the typical behaviour of the burst phases – starburst–quiescent– post-starburst – is qualitatively similar to what is found by Read et al. (2016b), who simulate the evolution of a $M_* \simeq 10^9 M_{\odot}$ galaxy for $\simeq 1$ Gyr (see also Teyssier et al. 2013; Read, Agertz & Collins 2016a for further specific studies on the bursty nature of this kind of galaxies).

Since individual galaxies are defined as a group of star particles in the same DM halo at $z=6$, the SF history accounts for the sum of all the stars that formed in different progenitors of the considered halo. For comparison, in the right-hand panel of Fig. 6 we plot the SFR of individual haloes defined by their merger history at $z=8.7$. Galaxies with active SF at 300–550 Myr merge into Dahlia at a later time; thus, they do not appear individually in the left-hand panel of Fig. 6.

Superimposed to the global trend, the SF history of Dahlia and its satellites fluctuates on time-scales of $\sim 10$ Myr, corresponding to the time-scale of energy deposition by feedback (see e.g. Torrey et al. 2016).

### 3.2.1 Star formation efficiency

$\zeta_{sf} f_{H_2}$ represents the quantity of gas converted in stars within a free-fall time (see equation 1). In Fig. 7, we plot the effective star formation efficiency ($\zeta_{sf} f_{H_2}$) as a function of gas density, weighted by the $H_2$ mass fraction at $z=6$. Most of the $H_2$ is contained in the range $n = 10–100 \text{ cm}^{-3}$, and the effective efficiency $\zeta_{sf} f_{H_2}$ varies from $10^{-3}$ to $10^{-1}$. Since $\zeta_{sf} = \text{const.} = 0.1$, the spread is purely due to the dependence of $f_{H_2}$ on density and metallicity (see Fig. 1). Note that by construction $\zeta_{sf} f_{H_2} \le 0.1$, and the plot does not show values very close to such limit, since gas with higher effective efficiency is converted into stars within a few free-fall times (equation 1b).

Interestingly, our $H_2$-based star formation criterion is reminiscent of a density threshold one, as below $n \simeq 3 \text{ cm}^{-3}$ the efficiency drops

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Effective star formation efficiency ($\zeta_{sf} f_{H_2}$) versus density ($n$) at $z=6$. The distribution is $H_2$ mass weighted; we consider gas within $3 r_{vir} = 47.5$ kpc from Dahlia's centre.</figcaption></figure>

abruptly (equations 3). However, an important difference remains; i.e. in the present model at any given density the efficiency varies considerably as a result of the metallicity dependence. The relation between efficiency and density is also similar to that found by Semenov et al. (2016, hereafter S16). This is striking because these authors use a star formation efficiency that depends on the turbulent velocity dispersion of the gas, with no notion of the local metallicity. This comparison is discussed further in Section 4.

### 3.2.2 Feedback energy deposition

As discussed in Section 2.4, only a small fraction of the available energy produced by stars can couple to the gas. During the simulation, we find that the time average efficiency of the conversion is $f \sim 0.1$ per cent, regardless of the feedback type. These low efficiencies imply that energy is mostly dissipated within MCs where the stars reside and produce it. For SN and winds, such small efficiency is a consequence of the short cooling times in MCs (see also
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Rate of energy deposition in the gas, $dE/dt$, by feedback processes as a function of cosmic time. Different contributions (SN, wind and radiation) are plotted with a different colour, and we additionally distinguish between the kinetic (thick lines) and thermal (thin lines) energy variation. By definition, radiation pressure has no thermal contribution. Note the jump at $t \approx 500$ Myr due to the onset of radiation pressure by AGB stars. The upper axis indicates the corresponding redshift.</figcaption></figure>

Appendix A). For radiation pressure, the efficiency is limited by the relatively small dust optical depths (see also Appendix B).

Note that, typically in simulations (e.g. Wise et al. 2012; Agertz et al. 2013), energy from stars is directly deposited in the gas, and then dissipation (mostly by radiative losses) occurs during the hydrodynamical time step. Within our scheme, instead, the deposited energy is already dissipated within high-density cells, where cooling is important. Nevertheless, this does not appear to determine major differences in, e.g. SFR history and ISM thermodynamics, as discussed in Section 3.2.3.

In Fig. 8, we plot the energy deposition rate in the gas by various feedback processes as a function of time. Most evidently, radiation dominates the energy budget at all times: $\dot{E}_{rad} \approx 10^2 \dot{E}_{SN} \approx 10^3 \dot{E}_{w}$. The ratios of these energy rates somewhat reflect the stellar inputs shown in Fig. 5, although this is not a trivial finding, given that the interplay among different feedback types is a highly non-linear process.

As expected, the energy deposition rate behaves as $\dot{E} \propto \text{SFR}^q$, with $q \gtrsim 1$, apart from fluctuations and jumps as the one at $t \approx 500$ Myr. The scaling can be understood by simple dimensional arguments. Assume that most of the energy is deposited by radiation pressure. In the optically thick limit, we can combine equations (7c) and (7a) to write $\dot{E}_{rad} \Delta t \approx (L_{uv} \Delta t)^2 / (M_g c^2)$, where $M_g$ is the gas mass accelerated by radiation, and we neglect ionizing radiation. Then, using equation (5e), we can write $\dot{E}_{rad} \propto \text{SFR} (M_*/M_g)$. Initially, $M_g \approx M_*$, thus, $\dot{E}_{rad} \propto \text{SFR}$. Once the gas mass is expelled from the star-forming region or converted into stars, $M_g \ll M_*$, Thus, the deposition rate increases faster than the SFR and it is very sensitive to the amount of gas mass around the sources.

The previous argument holds until the gas remains optically thick. This is warranted by AGB metal/dust production which becomes important after for stellar ages $t_* \sim 100$ Myr (see Fig. 2). When combined with the parallel increase of UV photons by the same sources, it is easy to interpret the rapid increase of the radiative feedback efficiency at $t \approx 500$ Myr, i.e. after $\approx 200$ Myr from the first star formation events in Dahlia. We checked this interpretation by looking at the IR trapping recorded on the fly during the simulation.

<table><caption>**Table 2.** Summary of the gas masses for total metal, C II and H$_2$ within ~47.5 kpc ($3 r_{vir}$) from Dahlia's centre. In the table, we report also the fraction that is contained in different gas phases$^{14}$: rarefied (log (n cm$^{-3}$) $\le -1$), diffuse ($-1 < \log (n \text{ cm}^{-3}) \le 1$) and dense (log (n cm$^{-3}$) > 1). Discussion about gas and metal mass is found in Section 3.2.3; analysis of H$_2$ and C II is in Section 3.3 (see also Appendix C for C II calculation).</caption><thead><tr><th></th><th>Mass<br>(M<sub>☉</sub>)</th><th>Rarefied<br>(per cent)</th><th>Diffuse<br>(per cent)</th><th>Dense<br>(per cent)</th></tr></thead><tbody><tr><td>Gas</td><td>1.3 × 10<sup>10</sup></td><td>44</td><td>34</td><td>22</td></tr><tr><td>Metals</td><td>4.2 × 10<sup>5</sup></td><td>5</td><td>25</td><td>70</td></tr><tr><td>H<sub>2</sub></td><td>3.6 × 10<sup>8</sup></td><td>0</td><td>1</td><td>99</td></tr><tr><td>C II</td><td>2.2 × 10<sup>5</sup></td><td>4</td><td>22</td><td>74</td></tr></tbody></table>

We find that on average $f_{ir} \approx 10^{-2}$ for $t \lesssim 500$ Myr, and $f_{ir} \approx 0.1$ at later times, thus confirming our hypothesis.

The energy deposition rates for different feedback types are highly correlated in time (Pearson coefficients $\gtrsim 0.7$). This is partially due to the fact that the same stellar population inputs wind, radiation and supernova energy in the gas. Additionally, as we have just seen for the case of AGB stars, different types of feedback are mutually dependent. For example, radiation pressure is more effective when the gas is metal and dust enriched by SN and AGB stars; winds and SN can more efficiently couple with low-density gas (longer cooling time).

Note that short and intense peaks in energy deposition rate correspond to the complete disruption of multiple MCs. This occurs following strong SF events in small satellites ($M_h \sim 10^7 M_☉$) that cannot retain the gas and sustain a continuous star formation activity.

Finally, we remind that, when compared with observational/analytical constraints, the SFR and $M_*$ of Dahlia are higher than the mean but still consistent within $1\sigma$. We caution that this might imply a somewhat weak feedback prescription.

### 3.2.3 Feedback effects on ISM thermodynamics

Feedback leaves clear imprints in the ISM thermodynamics. For convenience, we classify ISM phases according to their density: we define the gas to be in the rarefied, diffuse and dense phase if $n \le 0.1 \text{ cm}^{-3}$, $0.1 \le n/\text{cm}^{-3} \le 10$ and $n > 10 \text{ cm}^{-3}$, respectively.$^{14}$

We focus at $z = 6$ and consider the gas in a region within $\approx 47.5$ kpc ($3 r_{vir}$) from Dahlia's centre, essentially the scale of the CGM described in Section 3.1.2. This region contains a total gas mass of $1.3 \times 10^{10} M_☉$ and metal mass of $4.2 \times 10^5 M_☉$ (additional data in Table 2).

Fig. 9 shows the equation of state (EOS; or phase diagram) of the gas. The fraction of gas in the rarefied, diffuse and dense phases is 44, 34 and 22 per cent; these phases contain 5, 25 and 70 per cent of the metals, respectively. Thus, while the gas mass is preferentially located in the lower density phases, metals are mostly found in dense gas, i.e. star-forming regions/MCs. Additionally, only $\sim 30$ per cent of the considered volume shows $Z > 10^{-3} Z_{\odot} = Z_{\text{floor}}$; i.e. it has been polluted by stars in the simulation. We note that the EOS in the $n-T$ plane is fairly consistent with the one found in other high-z galaxy simulations (e.g. see fig. 5 in Wise et al. 2012). Comparison

<footnote>$^{14}$ Compared to the definitions used in Klessen & Glover (2016), the rarefied phase corresponds to the warm and hot ionized medium, the diffuse phase to the cold and warm neutral medium and the dense phase to the molecular gas.</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Equation of state (EOS) of the gas within $z \simeq 47.5$ kpc (3 $r_{vir}$) from Dahlia centre at $z=6$. Each EOS consists in a mass- or metal-weighted probability distribution function (PDF) as specified by the colour bar. We plot the PDF in the $n-T$ plane (upper left-hand panel), in the $n-P$ plane (upper right-hand panel), the metal-mass weighted PDF in the $n-T$ plane (lower left-hand panel) and mass-weighted relation between gas $n$ and $Z$ (lower right-hand panel). Mean relations and rms dispersions are overplotted with solid black and dashed lines, respectively. In the upper horizontal axis of each panel, we indicate the overdensity ($\Delta$) corresponding to $n$. The density range of rarefied, diffuse and dense phases used in the text is indicated. For the panels on the left, the rarefied gas is additionally divided in photoionized ($T < 10^{4.5}$ K) and shock heated ($T \ge 10^{4.5}$ K). See Table 2 for a summary of the total values.</figcaption></figure>

between the EOS in the $n-T$ and $n-P$ plane highlights the relative importance of different feedback types.

The rarefied gas is characterized by long cooling times. Thus, once engulfed by shocks, such a phase becomes mildly enriched ($\langle Z \rangle \sim 10^{-2} Z_{\odot}$) and remains hot ($T \sim 10^6$ K). The enriched rarefied gas preferentially populates the $n \simeq 10^{-3} \text{ cm}^{-3}$ and $T \simeq 10^{6.5}$ K region of the phase diagram. However, part of the rarefied gas has $T \simeq 10^4$ K. This gas component has a temperature set by the equilibrium between adiabatic cooling and the photoheating by the UV background; it feeds the accretion on to Dahlia, but it is not affected by stellar feedback. As such, it is not central in the present analysis.

The dense gas is mostly unaffected by shocks, and it is concentrated in the disc. Typically, such gas has $n \sim 10^2 \text{ cm}^{-3}$ and $T \sim 10^2$ K; thus, a thermal pressure $P_{\text{th}}/k \sim 10^4 \text{ cm}^{-3}$ K is expected. However, the total gas pressure is $P/k \sim 10^7 \text{ cm}^{-3}$ K (see the $P-n$ EOS). The extra contribution is provided in kinetic form by radiation pressure, thanks to the strong coupling with the gas allowed by the high optical depth of this phase. This leads to the important implication that the central structure of Dahlia is radiation supported (see also Section 3.3).

The diffuse gas acts as an interface between the dense disc gas and the rarefied gas envelope. Diffuse gas is found both in hot ($T \sim 10^5$ K) and cold ($T \sim 10^3$ K) states. The cold part has a sufficiently high mean metallicity, $z \sim 0.1 Z_{\odot}$, to allow an efficient cooling of the gas. This is highlighted by the metal-weighted EOS, where we can see that most of the metals present in the diffuse phase are cold.

Note that the phase diagram also shows evidence for the classical two-phase medium shape for pressures around $P/k \sim 10^3 \text{ cm}^{-3}$ K, while at higher (and lower) pressures only one stable phase is allowed; nevertheless, at any given pressure a range of densities can be supported. Such a situation, though, is highly dynamic and does not correspond to a true thermal equilibrium.

A final remark is that by $z=6$ an $n-Z$ correlation is already in place, although considerable scatter is present. The relation gets steeper at large densities, and at the same time the scatter decreases. Such a relation arises from the superposition of the analogous relation for metal bubbles of individual galaxies (Dahlia and satellites). The scatter instead results from the fact that the slope of the $n-Z$ relation depends on the SFR history (for an in-depth analysis see P14). The average $n-Z$ relation found is consistent with the results from $z \simeq 3$ galaxies (Shen et al. 2014).
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Face-on (upper panels) and edge-on (lower panels) $z = 6$ Dahlia surface maps for H$_2$ density ($\Sigma_{H_2}/(\text{M}_\odot \text{kpc}^{-2})$) left-hand panels, C II density ($\Sigma_{\text{CII}}/(\text{M}_\odot \text{kpc}^{-2})$) middle panels) and [C II] brightness ($S_{[\text{CII}]}/(\text{L}_\odot \text{kpc}^{-2})$) right-hand panels). The scale is 10 kpc, as in the rightmost panels of Fig. 4 (Section 3.1.3). Note that lower limits for the maps are drawn for visualization purposes ($\log(\Sigma_{H_2}/(\text{M}_\odot \text{kpc}^{-2})) \simeq \log(S_{[\text{CII}]}/(\text{L}_\odot \text{kpc}^{-2})) \simeq 5$, $\log(\Sigma_{\text{CII}}/(\text{M}_\odot \text{kpc}^{-2})) \simeq 2$). Additionally, an average of the maps is plotted in Fig. 11.</figcaption></figure>

## 3.3 Additional ISM properties

We conclude our analysis by inspecting the distribution of two key ISM species, molecular hydrogen and C II, along with the expected surface brightness of the corresponding 158 µm [C II] line. The surface maps of these quantities in Dahlia ($z = 6$) are shown in Fig. 10 for the face-on and edge-on view cases. For reference, in Fig. 11 we additionally plot the radially averaged profiles of the same quantities, and in Table 3 we give their typical radial scales.

### 3.3.1 Molecular hydrogen

Dahlia has a total H$_2$ mass $M_{H_2} \simeq 3.6 \times 10^8 \text{ M}_\odot$ that is mainly concentrated in a disc-like structure of radius $\simeq 0.6$ kpc and scale height $\simeq 200$ pc, with a sharp cut-off beyond these scales.<fnref n="15" /> The disc has mean surface density ($\Sigma_{H_2}) \simeq 10^{7.5} \text{ M}_\odot \text{kpc}^{-2}$ that is approximately constant with radius and presents perturbed spiral arms along which the density is enhanced by a factor of $\simeq 3$. The spiral arms are less pronounced than in a more massive, MW-like galaxy (see Ceverino, Primack & Dekel 2015, S16). This trend with mass has already been pointed out by Ceverino, Dekel & Bournaud (2010).

The disc is composed of dense ($n \gtrsim 25 \text{ cm}^{-3}$), enriched ($Z \simeq 0.5 \text{ Z}_\odot$), radiation-pressure-supported gas, as already discussed. It

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Radially averaged profiles for face-on (full colour) and edge-on (transparent and hatched) views of Dahlia at $z = 6$. Upper left: H$_2$ surface density; upper right: stellar surface density; lower left: C II surface density; lower right: [C II] surface brightness.</figcaption></figure>

<footnote n="15">¹⁵ Such scales are calculated by using the principal component analysis of the H$_2$ distribution around the galaxy.</footnote>
<table><caption>**Table 3.** Summary of the effective radius ($r_{1/2}$) for the H$_2$, [C II], C II and stellar component in Dahlia at $z = 6$ for the face-on and edge-on case and corresponding values for the mass/luminosity. For each entry, $r_{1/2}$ is defined as the radius including half of the mass/luminosity. The quoted radii have a reference error of $\pm 0.015$ kpc. Note that the approximate values of masses/luminosity at $r_{1/2}$ are insensitive to the orientation of the projection (face-on/edge-on). The full profiles are shown in Fig. 11.</caption><thead><tr><th rowspan="2"></th><th colspan="2">r<sub>1/2</sub>/kpc</th><th rowspan="2">Approximate<br>value at r<sub>1/2</sub></th></tr><tr><th>Face-on</th><th>Edge-on</th></tr></thead><tbody><tr><td>H<sub>2</sub></td><td>0.59</td><td>0.36</td><td>10<sup>8.23</sup> M<sub>☉</sub></td></tr><tr><td>C II</td><td>0.64</td><td>0.38</td><td>10<sup>5.14</sup> M<sub>☉</sub></td></tr><tr><td>Stars</td><td>0.37</td><td>0.23</td><td>10<sup>9.89</sup> M<sub>☉</sub></td></tr><tr><td>[C II]</td><td>0.60</td><td>0.36</td><td>10<sup>7.25</sup> L<sub>☉</sub></td></tr></tbody></table>

is fed by frequent mergers driving fresh gas to the centre, and supports an SFR per unit area of $\simeq 15$ M$_{☉}$ yr$^{-1}$ kpc$^{-2}$, i.e. more than 1000 times the Milky Way value. Fragmentation of the disc is relatively weak (cf. Mayer et al. 2016), as indicated by a smooth surface density map, and also paralleling the flat metallicity profile in the inner $\simeq$ kpc. For the fragmentation of the H$_2$ component, we caution that this result has been obtained assuming a uniform UV interstellar field; stronger fragmentation in the H$_2$ distribution may occur when accounting for local radiation sources: Lyman-Werner photons from these sources might in fact locally dissociate the H$_2$ by generating pockets of H$_1$ in the distribution.

While most of the H$_2$ gas resides in the disc, we can clearly distinguish three clumps of molecular gas both in the face-on and edge-on maps. These clumps are located a few kpc away from the centre and are characterized by sizes of $\sim 150$ pc and $M_{H_2} \sim 5 \times 10^6$ M$_{☉}$. Such clumps are Jeans-unstable and form stars as they infall and stream through the CGM, as it can be appreciated by comparing the H$_2$ and stellar mass profiles (Fig. 11). The stellar mass profiles also highlight the presence of three stellar clumps at $t \sim 1$ kpc with no associated H$_2$. These ‘older’ clumps share the same nature of the previous ones, but the H$_2$ has been already consumed and/or dispersed by the star formation activity that produced the stars present at $z = 6$.

### 3.3.2 Singly ionized carbon

The C II abundance is calculated by post-processing the simulation outputs with the photoionization code CLOUDY (Ferland et al. 2013, and see Appendix C). The result is shown in Fig. 10. Dahlia contains a C II mass of $M_{C II} = 2.2 \times 10^5$ M$_{☉}$, accounting for $\sim 50$ per cent of the total metals produced. About 74 per cent of the C II mass is located in the dense phase, 22 per cent in the diffuse phase and 4 per cent in the rarefied phase. Note that the C II mass phase distribution differs only for $\lesssim 10$ per cent from the Z distribution (see Table 2). The difference arises because shock-heated gas can be collisionally excited to higher ionization states. Thus, to first order, we expect the C II spatial distribution to follow the metallicity one.

The face-on C II surface density has a central maximum ($\Sigma_{C II} \sim 10^5$ M$_{☉}$ kpc$^{-2}$), it gradually decreases to up to $\simeq 1.2$ kpc, and drastically drops to $\Sigma_{C II} \lesssim 10^2$ M$_{☉}$ kpc$^{-2}$ beyond that radius (see also Fig. 11). Thus, most of the C II is located into the disc, but a more extended envelope containing a sizable fraction of mass exists. On top of this smooth distribution, there are C II enhancements corresponding to the H$_2$ clumps described above.

The C II profile is similar for edge-on and face-on case. However, the edge-on has a higher C II central density and a steeper slope. While the higher central value is obviously due to the larger column density encountered along the disc, the sharp drop is related to metal

transport. As most of the star formation activity is located in the disc, metals above it can be only brought by outflows which become progressively weaker with distance. Metal outflows originating from the centre are preferentially aligned with the rotation axis, and the pollution region starting from the edge is stretched by the disc rotation and by tidal interaction with satellites.

### 3.3.3 Emission from singly ionized carbon

We finally compute the expected [C II] line emission using the same prescriptions of Vallini et al. (2013, 2015), as detailed in Appendix C. Note that for this work we assume uniform UV interstellar radiation. This approximation is valid in the MW, where variations around the mean field value are limited to a factor of 3. The results are plotted in Fig. 10.

Within 1 kpc from the centre the [C II] emission structure closely follows the C II distribution, and we find $S_{[C II]}/L_{\odot} \simeq 200 \Sigma_{C II}/M_{\odot}$. At larger radii, the [C II] surface brightness suddenly drops, although the peaks associated with H$_2$ clumps are preserved. This result holds both for the face-on and edge-on cases.

Such behaviour can be understood as follows. Take a typical MC with $n = 10^2$ cm$^{-3}$, $Z = Z_{\odot}$ and total mass M. Its [C II] luminosity is $L_{[C II]}/L_{\odot} \simeq 0.1(M/M_{\odot})$ (Goicoechea et al. 2015; Vallini et al. 2016). Also, the [C II] emission is $\propto Zn$ for $n \lesssim 10^3$, i.e. the critical density for C II collisional excitation by H atoms (Vallini et al. 2013). Then,

$$L_{[C II]} \simeq 0.1 \left( \frac{n}{100 \text{ cm}^{-3}} \right) \left( \frac{Z}{Z_{\odot}} \right) \left( \frac{M}{M_{\odot}} \right) L_{\odot}. \quad (8)$$

In the central kpc, where $n \simeq 10^2$ cm$^{-3}$ and $Z \simeq Z_{\odot}$, the luminosity depends only on the molecular mass contained in the disc, and the same holds even for H$_2$ clumps outside the disc. The envelope is instead more diffuse ($n \lesssim 10^4$ cm$^{-3}$) and only mildly enriched ($Z \lesssim 10^{-1}$ Z$_{\odot}$). As a result, its [C II] luminosity per unit mass is lower.

The emission from this diffuse component is further suppressed by the CMB (da Cunha et al. 2013; Pallottini et al. 2015; Vallini et al. 2015). Namely, for gas with $n \lesssim 0.1$ cm$^{-3}$, the upper levels of the [C II] transition cannot be efficiently populated through collisions, thus the spin temperature of the transition approaches the CMB one, and to a first order the gas cannot be observed in emission.

In summary, $\simeq 95$ per cent of Dahlia [C II] emission comes from dense gas located in the H$_2$ disc. Indeed, the [C II] half-light radius coincides with the H$_2$ half-mass radius, i.e. 0.59 kpc (0.36 kpc) in the face-on (edge-on) case (see also Table 3). Within such radius, the molecular gas has a mass $M_{H_2} \simeq 1.69 \times 10^8$ M$_{☉}$ and the luminosity is $L_{C II} \simeq 1.78 \times 10^7$ M$_{☉}$, i.e. with a [C II]-H$_2$ scaling ratio consistent within 15 per cent from the simple estimate in equation (8).

Dahlia has a total [C II] luminosity $L_{C II} \simeq 3.5 \times 10^7$ L$_{\odot}$; this is fainter than expected on the basis of the local [C II]-SFR relation ($L_{C II} \sim 10^8-10^9$ L$_{\odot}$; De Looze et al. 2014). However, at high-z, such relation seems to hold only for a small subset of the observed galaxies (Capak et al. 2015; Willott et al. 2015). The majority of the observed galaxies show a strong [C II]-SFR deficit, when considering both detections (e.g. BDF 3299, A383-5.1; Maiolino et al. 2015; Knudsen et al. 2016) and upper limits (e.g. Himiko, IOK1, MS0451-H; Ouchi et al. 2013; Ota et al. 2014; Knudsen et al. 2016).

For Dahlia, the [C II]-SFR deficit depends on multiple factors. The main contribution from [C II] emission is in the H$_2$ disc, that on average has $\langle Z \rangle \simeq 0.5$ Z$_{\odot}$, i.e. slightly lower than solar. Additionally, the gas in the disc is efficiently converted in stars (SFR $\simeq 100$ M$_{☉}$ yr$^{-1}$) and has $\langle n \rangle \simeq 25$ cm$^{-3}$, thus the [C II] emission
is hindered (equation 8). Finally, there is a marginal contribution to [C II] from the diffuse and rarefied phase: ≃30 per cent of C II is locked in the low density and metallicity gas that gives a negligible contribution to [C II] emission, particularly because of CMB suppression.

## 4 SUMMARY AND DISCUSSION

With the aim of characterizing the internal properties of high-z galaxies, we have performed an AMR zoom-in simulation of ‘Dahlia’, a z ≃ 6 galaxy with a stellar mass of $M_* = 1.6 \times 10^{10}$ M$_{\odot}$, therefore representative of LBGs at that epoch. We follow the zoom-in region with a gas mass resolution of $10^4$ M$_{\odot}$ and a spatial resolution of 30 pc.

The simulation contains a rich set of physical processes. We use a star formation prescription based on a H$_2$ dependent Schmidt–Kennicutt relation. The H$_2$ abundance is computed from the KTM09 model (Fig. 1). Using stellar evolutionary models (Bertelli et al. 1994; Leitherer et al. 1999), we include chemical, radiative and mechanical energy inputs, accounting for their time evolution and metallicity dependence on the stellar population properties (Fig. 2). We include feedback from SN, winds and radiation pressure with a novel, physically motivated coupling scheme between gas and stars. We also compute C II abundance and the 158 µm [C II] emission, by post-processing the outputs with CLOUDY (Ferland et al. 2013), and a FIR emission model drawn from radiative transfer numerical simulations (Vallini et al. 2013, 2015).

The main results can be summarized as follows.

(1) Dahlia sits at the centre of a cosmic web knot, and accretes mass from the IGM mainly via three filaments of length ≃100 kpc (Fig. 4). Dahlia has ~6 major satellites ($M_* \lesssim 10^9$ M$_{\odot}$) and is surrounded by ~10 minor ones ($M_* \sim 10^5$ M$_{\odot}$). The latter represent MC complexes caught in the act of condensing as the gas streams through the CGM (Fig. 5). Dahlia dominates both the stellar mass ($M_* \sim 10^{10}$ M$_{\odot}$) and the SFR of the galaxy ensemble (SFR ≃ 100 M$_{\odot}$ yr$^{-1}$, Fig. 6).

(2) Only a small fraction of the available energy produced by stars couples to the gas, as energy is mostly dissipated within MCs where the stars reside. Radiation dominates the feedback energy budget by a factor of >100 (Fig. 8).

(3) By $z=6$ Dahlia forms a H$_2$ disc of mass of $M_{H2} = 3.6 \times 10^8$ M$_{\odot}$, effective radius 0.6 kpc and scale height 200 pc (Fig. 10). The disc is dense ($n \gtrsim 25$ cm$^{-3}$), enriched ($Z \simeq 0.5$ Z$_{\odot}$) and it is fed by frequent mergers driving fresh gas to the centre, and supports a SFR per unit area of ≃15 M$_{\odot}$ yr$^{-1}$ kpc$^{-2}$.

(4) The disc is mostly unaffected by SN shocks, and it is pressure supported by radiation. SN/winds drive hot metal outflows (Fig. 9) that are either preferentially aligned with the galaxy rotation axis or start at the edge of the disc.

(5) The total [C II] luminosity of Dahlia is $10^{7.55}$ L$_{\odot}$, and ≃95 per cent of the emission is co-located with the H$_2$ disc (Fig. 11). The diffuse, enriched material surrounding Dahlia contains 30 per cent of the C II mass, but it negligibly contributes to the [C II] emission (Fig. 10) due to its low density ($n \simeq 10$ cm$^{-3}$) and metallicity ($Z \simeq 10^{-1}$ Z$_{\odot}$). Dahlia is underluminous with respect to the local [C II]–SFR relation; however, its luminosity is consistent with upper limits derived for most $z \sim 6$ galaxies.

We find clear indications that the SF subgrid prescription might considerably affect the [C II]–SFR relation and the ISM structure, as noted also by Hopkins et al. (2013). This is because stars form in gas of different densities depending on the chosen prescription. In our simulation, gas is converted into stars with an efficiency $\zeta_{sf} f_{H2}$,

where the H$_2$ fraction is computed from the KTM09 model and we set $\zeta_{sf} = 0.1$. In S16 the SF follows a total (i.e. not molecular) density Schmidt–Kennicutt relation. Further the SF efficiency depends on the free-fall time and the turbulent eddy turnover time. The SF relation is derived from an empirical fit to MC simulations (Padoan, Haugbølle & Nordlund 2012), with no notion of the local metallicity.

Interestingly, although the approaches are considerably different, the resulting efficiencies are compatible: in S16 the bulk of the star-forming gas has $n \sim 10^{1.5}$ cm$^{-3}$, as in Dahlia (Fig. 7). However, with respect to S16, Dahlia misses part of the very dense, star-forming gas, and its corresponding contribution to [C II] from $z \sim Z_{\odot}$ MCs with $n \sim 10^3$ cm$^{-3}$. These MC are expected to have high [C II] fluxes (see equation 8), but their abundance might be low (Padoan et al. 2012). Further investigation is needed before we draw any solid conclusion. To this aim, we plan to upgrade our simulations to a more sophisticated non-equilibrium H$_2$ evolution model. This is because the chemical equilibrium assumed in KTM09 does not hold in low-metallicity regimes.

Another important caveat is that we have assumed a uniform UV background. Instead, discrete sources (stellar clusters) might have a strong impact on star formation. For example, Lyman–Werner photons might locally dissociate the H$_2$ by generating pockets of H$_1$ in the gas distribution. Thus, unshielded (low dust column density) gas in the disc would contribute only marginally to the SFR.

Furthermore, a uniform UVB assumption likely leads to inaccurate computation of the ISM thermodynamic state. We find that $Z \simeq 10^{-3}$ Z$_{\odot}$ gas with $n \gtrsim 10^2$ cm$^{-3}$ has $T \simeq 10^4$ (Fig. 9), with the temperature been set by the UVB heating. However, such gas should be likely able to self-shield from the impinging UVB, whereas internal radiation sources could still play a role (e.g. Gnedin 2010).

Finally, local FUV flux variations can change the [C II] emission from individual regions of the galaxy. Also, very high FUV fluxes can photoevaporate MC on short time-scales ($\lesssim t_{\text{ff}}$ for gas with $z \sim 10^{-2}$ Z$_{\odot}$; Vallini et al. 2016). This effect are particularly important, as it might be responsible for the displacement between the [C II] and the UV emitting region observed in BDF 3299 (Maiolino et al. 2015), and in some of the Capak et al. (2015) galaxies. To solve these problems, a multifrequency radiative transfer computation must be coupled to the present simulations. This work is ongoing and will be presented elsewhere.

## ACKNOWLEDGEMENTS

We are grateful to the participants of *The Cold Universe* program held in 2016 at the KITP, UCSB, for discussions during the workshop. We acknowledge the AGORA project members and the DAVID group for stimulating discussion. We thank the authors and the community of PYMSES for their work. We thank B. Smith for support in implementing GRACKLE. This research was supported in part by the National Science Foundation under Grant No. NSF PHY11-25915. SS was supported by the European Research Council through a Marie–Skodolowska–Curie Fellowship, project PRIMORDIAL-700907. RM acknowledges support from the ERC Advanced Grant 695671 ‘QUENCH’ and from the Science and Technology Facilities Council (STFC).

## REFERENCES

* Abel N. P., 2006, MNRAS, 368, 1949

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

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

* Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481
Badjin D. A., Glazyrin S. I., Manukovskiy K. V., Blinnikov S. I., 2016, MNRAS, in press

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

Barnes L. A., Garel T., Kacprzak G. G., 2014, PASP, 126, 969

Bauer I., Finocchi F., Duschl W. J., Gail H.-P., Schloeder J. P., 1997, A&A, 317, 273

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

Bertschinger E., 1995, preprint (astro-ph/9506070)

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

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

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

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

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

Ceverino D., Primack J., Dekel A., 2015, MNRAS, 453, 408

Ceverino D., Sánchez Almeida J., Muñoz Tuñón C., Dekel A., Elmegreen B. G., Elmegreen D. M., Primack J., 2016, MNRAS, 457, 2605

Cicone C. et al., 2015, A&A, 574, A14

Cioffi D. F., McKee C. F., Bertschinger E., 1988, ApJ, 334, 252

Cox D. P., 1972, ApJ, 178, 159

Crawford M. K., Genzel R., Townes C. H., Watson D. M., 1985, ApJ, 291, 755

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

Dale J. E., 2015, New Astron. Rev., 68, 1

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

Dayal P., Hirashita H., Ferrara A., 2010, MNRAS, 403, 620

Dayal P., Ferrara A., Dunlop J. S., Pacucci F., 2014, MNRAS, 445, 2545

de Bennassuti M., Schneider R., Valiante R., Salvadori S., 2014, MNRAS, 445, 3039

De Breuck C., Maiolino R., Caselli P., Coppin K., Hailey-Dunsheath S., Nagao T., 2011, A&A, 530, L8

Dekel A., Mandelker N., 2014, MNRAS, 444, 2071

Dekel A. et al., 2009, Nature, 457, 451

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

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

Draine B. T., 2003, ARA&A, 41, 241

Draine B. T., 2011, ApJ, 732, 100

Draine B. T., Lee H. M., 1984, ApJ, 285, 89

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

Eisenstein D. J., Hu W., 1998, ApJ, 496, 605

Eisenstein D. J., Hut P., 1998, ApJ, 498, 137

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

Feng Y., Di Matteo T., Croft R., Tenneti A., Bird S., Battaglia N., Wilkins S., 2015, ApJ, 808, L17

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

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

Fierlinger K. M., Burkert A., Ntormousi E., Fierlinger P., Schartmann M., Ballone A., Krause M. G. H., Diehl R., 2016, MNRAS, 456, 710

Gallerani S. et al., 2012, A&A, 543, A114

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

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é I., Bouwens R. J., Illingworth G., Franx M., Kriek M., Brammer G. B., 2010, ApJ, 713, 115

Gorti U., Hollenbach D., 2002, ApJ, 573, 215

Grassi T., Bovino S., Haugboelle T., Schleicher D. R. G., 2016, MNRAS, in press (arXiv:1606.01229)

Guillet T., Teyssier R., 2011, J. Comput. Phys., 230, 4756

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

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

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

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

Hollenbach D. J., Tielens A. G. G. M., 1999, Rev. Modern Phys., 71, 173

Hopkins P. F., Quataert E., Murray N., 2011, MNRAS, 417, 950

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

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

Janssen A. W. et al., 2016, ApJ, 822, 43

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

Kauffmann G., 2014, MNRAS, 441, 2717

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

Kim C.-G., Ostriker E. C., 2015, ApJ, 802, 99

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

Klessen R. S., Glover S. C. O., 2016, Star Formation in Galaxy Evolution: Connecting Numerical Models to Reality. Saas-Fee Advanced Course, Vol. 43. Springer-Verlag, Berlin, p. 85

Knebe A. et al., 2013, MNRAS, 435, 1618

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

Körtgen B., Seifried D., Banerjee R., Vázquez-Semadeni E., Zamora-Avilés M., 2016, MNRAS, 459, 3460

Koyama H., Inutsuka S.-i., 2002, ApJ, 564, L97

Kreckel K. et al., 2014, ApJ, 790, 26

Kroupa P., 2001, MNRAS, 322, 231

Krumholz M. R., Gnedin N. Y., 2011, ApJ, 729, 36

Krumholz M. R., Matzner C. D., 2009, ApJ, 703, 1352

Krumholz M. R., Thompson T. A., 2012, ApJ, 760, 155

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

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

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

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

Leitherer C., Ortiz Otálvaro P. A., Bresolin F., Kudritzki R.-P., Lo Faro B., Pauldrach A. W. A., Pettini M., Rix S. A., 2010, ApJS, 189, 309

Li A., Draine B. T., 2001, ApJ, 554, 778

Liang C. J., Kravtsov A. V., Agertz O., 2016, MNRAS, 458, 1164

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

McLeod A. D., Whitworth A. P., 2013, MNRAS, 431, 710

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

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

Madau P., Ferrara A., Rees M. J., 2001, ApJ, 555, 92

Madden S. C., Poglitsch A., Geis N., Stacey G. J., Townes C. H., 1997, ApJ, 483, 200

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., 2005, A&A, 440, L51

Maiolino R. et al., 2012, MNRAS, 425, L66

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

Mannucci F., Cresci G., Maiolino R., Marconi A., Gnerucci A., 2010, MNRAS, 408, 2115

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

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

Muñoz J. A., Furlanetto S. R., 2014, MNRAS, 438, 2483

Murray N., 2011, ApJ, 729, 133

Nagamine K., Wolfe A. M., Hernquist L., 2006, ApJ, 647, 60

Nisikawa A., Habe A., Isibasi N., 1997, preprint (astro-ph/9711118)
Olsen K. P., Greve T. R., Narayanan D., Thompson R., Toft S., Brinch C., 2015, ApJ, 814, 76

Oñorbe J., Garrison-Kimmel S., Maller A. H., Bullock J. S., Rocha M., Hahn O., 2014, MNRAS, 437, 1894

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. Modern Phys., 60, 1 (OM88)

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'Olorico V., 2014a, MNRAS, 440, 2498 (P14)

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

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

Planck Collaboration XVI, 2014, A&A, 571, A16

Poole G. B., Angel P. W., Mutch S. J., Power C., Duffy A. R., Geil P. M., Mesinger A., Wyithe S. B., 2016, MNRAS, 459, 3025

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

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

Raymond J. C., Cox D. P., Smith B. W., 1976, ApJ, 204, 290

Read J. I., Agertz O., Collins M. L. M., 2016a, MNRAS, 459, 2573

Read J. I., Iorio G., Agertz O., Fraternali F., 2016b, MNRAS, 462, 3628

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

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

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

Salvadori S., Ferrara A., Schneider R., 2008, MNRAS, 386, 348

Salvadori S., Dayal P., Ferrara A., 2010, MNRAS, 407, L1

Salvadori S., Skúladóttir Á., Tolstoy E., 2015, MNRAS, 454, 1320

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

Schmidt M., 1959, ApJ, 129, 243

Semenov D., Henning T., Helling C., Ilgner M., Sedlmayr E., 2003, A&A, 410, 611

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

Shen S., Madau P., Conroy C., Governato F., Mayer L., 2014, ApJ, 792, 99

Sheth R. K., Tormen G., 1999, MNRAS, 308, 119

Silva M., Santos M. G., Cooray A., Gong Y., 2015, ApJ, 806, 209

Stanway E. R., Bunker A. J., McMahon R. G., 2003, MNRAS, 342, 439

Stark D. P., Ellis R. S., Bunker A., Bundy K., Targett T., Benson A., Lacy M., 2009, ApJ, 697, 1493

Stark D. P., Schenker M. A., Ellis R., Robertson B., McLure R., Dunlop J., 2013, ApJ, 763, 129

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

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

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

Tielens A. G. G. M., Hollenbach D., 1985, ApJ, 291, 747

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., 2016, MNRAS, in press (arXiv:1601.07186)

Troncoso P. et al., 2014, A&A, 563, A58

Valiante R., Schneider R., Bianchi S., Andersen A. C., 2009, MNRAS, 397, 1661

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

Vallini L., Ferrara A., Pallottini A., Gallerani S., 2016, MNRAS, in press (arXiv:1606.08464)

van den Hoek L. B., Groenewegen M. A. T., 1997, A&AS, 123, 305

Walch S., Naab T., 2015, MNRAS, 451, 2757

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

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

Wolfire M. G., Hollenbach D., McKee C. F., Tielens A. G. G. M., Bakes E. L. O., 1995, ApJ, 443, 152

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

Wolfire M. G., Tielens A. G. G. M., Hollenbach D., Kaufman M. J., 2008, ApJ, 680, 384

Woosley S. E., Weaver T. A., 1995, ApJS, 101, 181

Wuyts E. et al., 2016, ApJ, 827, 74

Yue B., Ferrara A., Pallottini A., Gallerani S., Vallini L., 2015, MNRAS, 450, 3829

## APPENDIX A: BLAST WAVE MODEL

In this section, we present the model used to calculate the kinetic ($f_{kn}$) and thermal ($f_{th}$) energy fractions that a gas cell acquires during an SN explosion/because of a stellar wind (see equation 6c).

First, we consider an SN explosion. Then, the picture is the following (e.g. Cioffi et al. 1988; Ostriker & McKee 1988; Walch & Naab 2015). In the first stage of the SN blast, the stellar ejecta follow a free expansion solution that ends when the blast has swept a mass of material roughly the mass of the ejecta, ~1 M⊙ per SN. Then, the shock enters in a ST stage. The ST stage is adiabatic, and the available energy is divided in thermal and kinetic parts that account for $f_{th} \simeq 0.7$ and $f_{kn} \simeq 0.3$ of the total,<fnref n="16" /> respectively. The ST stage ends when radiative losses becomes important, at $t \simeq t_{cool}$, the cooling time of the ambient gas. At the contact discontinuity, the gas cooling causes the formation of a thin shell (shell formation stage, SF), and after that the shock proceeds snowplowing through the ambient medium, driven by the pressure of the gas interior (pressure-driven snowplow, PDS). When all the thermal energy is radiated away, then the blast wave enters in the so-called momentum-conserving snowplow (MCS). In the MCS stage, the momentum is conserved, while the remaining energy, purely kinetic, is gradually lost because of the work done by the blast on the ambient material, and eventually the blast stops.

In each stage, the blast wave can be modelled by following the analysis by Ostriker & McKee (1988, hereafter OM88). A self-similar blast wave can be described by the evolution of the shock front at time $t$ as $r_s \propto t^{\eta}$, where $\eta$ is a constant that determines the type of blast (ST, MCS, ...). The velocity of propagation of the shock is $v_s = \eta r_s / t$, and, using the virial theorem (see equation 3.3 in OM88), the total energy of a blast wave can be written as

$$E(t) = (4\pi/3)\sigma_{\eta}\rho r_s^3 v_s^2 \propto t^{5\eta-2}, \quad (A1)$$

where $\sigma_{\eta}$ is a dimensionless constant and $\rho$ is the density of the material swept by the blast. Within the presented formalism (equation 6c), $f = E(t)/E_0$, where $E_0$ is the input energy, i.e. $E_0 = [\epsilon_{sn}(t_*) + \Delta t] - \epsilon_{sn}(t_*)]M_*$ (equation 5c).

The cooling time is critical in calculating the evolution of the blast between different stages. Here the cooling time $t_{cool}$ is defined as the time $t$ when $t = K_b T_s / (\rho_s / m_p \Lambda)$, where $K_b$ is the Boltzmann constant, $T_s = T_s(t)$ and $\rho_s$ are the temperature and density at the shock front, $m_p$ is the proton mass and $\Lambda = \Lambda(T_s, Z)$ is the cooling function.

Note that calculating the cooling function with GRACKLE or analytical approximations (e.g. Raymond, Cox & Smith 1976; Nisikawa, Habe & Isibasi 1997; Koyama & Inutsuka 2002) yields comparable results. This happens because the gas starts to cool at $T_s \gtrsim 10^7$ K,

<footnote n="16">¹⁶ Such values of $f_{th}$ and $f_{kn}$ are calculated by assuming $\gamma = 5/3$ for the gas adiabatic index. See table III in OM88 for the general solution.</footnote>
when Bremsstrahlung is the dominant cooling process, and it is in-
dependent of the metallicity, i.e. $\Lambda \propto T^{-1/2}$. Thus, as in Cioffi et al. (1988, see equation 3.10), the cooling time can be approximated as $t_{\text{cool}} = 3.61 \times 10^{-2}(E_0/\text{foe})^{3/14}(\rho/m_{\text{p}}\text{cm}^{-3})^{4/7}$ Myr (see also Kim & Ostriker 2015, in particular see equations 6 and 7).

The typical range of input energy is $1 \le E_0/\text{foe} \lesssim 10^3$ and ISM density is $10^{-1} \lesssim \rho/m_{\text{p}}\text{cm}^{-3} \lesssim 10^3$; thus, the cooling time is in the range $10^{-3} \lesssim t_{\text{cool}}/\text{Myr} \lesssim 1$. Since the expected simulation time step is $\Delta t \sim 10^{-2}$ Myr, we can further simplify the blast wave picture.<fnref n="17" />

The free expansion stage is shorter than our typical simulation time step ($t \lesssim 10^{-4}$ Myr; e.g. see equation 1 in Kim & Ostriker 2015); thus, we assume that the SN starts in the ST stage ($\eta = 2/5$). After $\simeq t_{\text{cool}}$, the energy of the shock is roughly half the initial value and the blast wave is in the PDS stage ($\eta = 2/7$). In both stages, the total energy is given by the sum of kinetic and thermal terms, and the relative fractions of kinetic ($f_{\text{kn}}$) and thermal ($f_{\text{th}}$) energies are constants (see equations 3.16 and 3.18 in OM88) that depends on $\eta$ and on the internal structure of the blast wave. In the intermediate SF stage, we approximate $f_{\text{kn}}$ and $f_{\text{th}}$ by linearly interpolating between the ST and initial PDS values. The time when the blast enters in the SF and PDS stages can be calculated following the analytical fit to the simulation presented in Cioffi et al. (1988), i.e. $0.14 t_{\text{cool}}$ and $0.4 t_{\text{cool}}$ for SF and PDS, respectively (see their equation 3.15). An example of such a model is presented in Section 2.4 (in particular, see Fig. 3). Note that the model is consistent with the result from SN exploding in an homogeneous medium; however, a blast wave propagating in an inhomogeneous medium behaves differently, since the blast travels unimpeded through the path of the lower density medium (Martizzi et al. 2015, see in particular the fit in equations 9 and 10 and the different scaling given in equations 11 and 12). We will explore this aspect in a future work.

The winds can be treated within the same blast wave formalism (see section VII in OM88). The stage evolution of the winds is similar to the SN one: first the wind is adiabatic, then the outer shock begins to radiate away and a thin shell is formed and eventually the shock becomes momentum conserving. However, in the wind case the energy injection is not impulsive, as the stars input a continuous energy injection with a (roughly) constant luminosity. Thus, we have a different blast wave structure and consequently $f_{\text{th}}$ and $f_{\text{kn}}$ in each stage.

The wind solution is given by Weaver et al. (1977, hereafter W77). The structure can be spatially divided in three parts (see fig. 1 in W77): near to the stars we have the stellar wind (a); the gas is contained in the region of shocked stellar wind (b) and in the shell of shocked gas (c).

In the adiabatic case, we have a (roughly) constant wind luminosity ($L$) and no radiative losses, thus, dimensional analysis implies $\eta = 3/5$ in equation (A1). In this stage, the gas in region (b) has 5/11 of the total energy, which is in a purely thermal form, while the gas in (c) contains energy in both thermal and kinetic form. Summing the contribution of both regions, W77 finds $f_{\text{th}} \simeq 0.78$ and $f_{\text{kn}} \simeq 0.22$ (see also equation 7.11 in OM88). Note that the thermal energy fraction is larger to the corresponding one in the ST stage of the SN-driven blast.

In the next stage, radiative losses become important ($t > t_{\text{cool}}$), and region (c) collapses to a thin shell. All the energy ($5/11 L t$) would be contained in the gas in (b) in thermal energy. The physical situation

is that the injected fluid is adiabatic, while in the ambient radiative stage losses would be dominant. This happens because the density of the injected fluid is much less than that of the ambient medium, and this consequently affects the cooling time-scales. However, the modellization of this stage cannot be readily implemented as a subgrid model in our simulation, because we cannot easily keep track of the internal structure of the cells in the simulation.

In our code, we opt to go directly from the adiabatic stage ($f_{\text{th}} = 0.78$, $f_{\text{kn}} = 0.22$, W77) to the momentum-conserving stage when $t > t_{\text{cool}}$. The latter is described by $f_{\text{kn}} = 1$, and $E(t)$ is given by equation (A1) with $\eta = 1/4$ (see equation 7.20 in OM88).

As noted in Agertz et al. (2013), Walch & Naab (2015), Fierlinger et al. (2016) and Körtgen et al. (2016), we expect wind injection to make the SN more efficient, since the SN blast sweeps through a gas with a lower density, thus with a longer cooling time. This point applies both to wind and radiation pressure injection.

Finally, note that in the blast wave modellization, we have neglected the instabilities that can arise during the thin-shell formation stage (Madau, Ferrara & Rees 2001; McLeod & Whitworth 2013; Badjin et al. 2016), and we do not explicitly consider the effect of multiple blast wave events (Walch & Naab 2015; Fierlinger et al. 2016). This issue will be addressed by a future work.

## APPENDIX B: RADIATION PRESSURE: OPTICAL DEPTH AND IR TRAPPING

The radiation pressure is dependent on the rate of momentum injection ($\dot{p}_{\text{rad}}$, equation 7a), which in turn depends on the optical depth to ionizing photons ($\tau_{\text{ion}}$) and the IR trapping of the UV radiation ($f_{\text{ir}}$).

For each gas cell, the optical depth $\tau_{\text{ion}}$ can be calculated by averaging the ionization cross-section on the stellar spectra assumed in our model (see equation 5e) which are taken from STARBURST99 (Leitherer et al. 1999, 2010). Namely, $\tau_{\text{ion}} = \sigma_{\text{ion}}N_{\text{H}}$, where ionization cross-section $\sigma_{\text{ion}}$ is obtained as a Rosseland mean:

$$ \sigma_{\text{ion}} = \int_{912 \AA}^{4000 \AA} L_{\lambda} \sigma_{\text{ion}}^{\lambda} d\lambda / L_{\text{ion}}, \qquad (\text{B1a}) $$

where the cross-section as a function of wavelength is given by Osterbrock (1989):

$$ \sigma_{\text{ion}}^{\lambda} = 6.3 \times 10^{-18} (1.34 (912 \AA / \lambda)^{-2.99} - 0.34 (912 \AA / \lambda)^{-3.99}). \qquad (\text{B1b}) $$

For the IR trapping, we assume $f_{\text{ir}} = \tau_{\text{ir}}$, the infrared dust optical depth, similarly to Agertz et al. (2013). As noted in Krumholz & Thompson (2012), the assumption $f_{\text{ir}} = \tau_{\text{ir}}$ is approximately correct for small values of the optical depth, i.e. $\tau_{\text{ir}} \lesssim 20$. This is our case, since in simulations with resolution scale larger than 10 pc, we expect $\tau_{\text{ir}} \lesssim 10$ (Rosdahl et al. 2015). As reported in Section 3.2.2, $\tau_{\text{ir}} \lesssim 10^{-1}$ in our simulation (see also Fig. 8).

Using the cross-section from Draine (2003, see also Semenov et al. 2003), we can write

$$ \tau_{\text{ir}} \simeq 1.79 \times 10^{-24} (N_{\text{H}}/\text{cm}^{-2})(Z/Z_{\odot})(T_{\text{dust}}/100 \text{ K})^2, \qquad (\text{B2a}) $$

where $T_{\text{dust}}$ is the dust temperature. $T_{\text{dust}}$ can be calculated as (Dayal, Hirashita & Ferrara 2010)

$$ T_{\text{dust}} = 6.73 (L_{\text{ir}}/L_{\odot})^{1/6} (M_{\text{dust}}/M_{\odot})^{-1/6}, \qquad (\text{B2b}) $$

where $L_{\text{ir}}$ is the reprocessed UV luminosity, i.e. $L_{\text{ir}} = L_{\text{uv}}(1 - \exp(-\tau_{\text{uv}}))$, and $M_{\text{dust}}$ is the dust mass, i.e.

$$ M_{\text{dust}} = D_{\odot} \rho V_{\text{cell}}(Z/Z_{\odot}), \qquad (\text{B2c}) $$

<footnote n="17">¹⁷ For the complete picture of blast wave evolution, we refer the reader to OM88, in particular fig. 1, table IV and references therein.</footnote>
with $D_⊙ = 6 \times 10^{-3}$ being the solar dust-to-gas ratio (e.g. Hirashita & Ferrara 2002).

We account for CMB heating on dust, which can be important at high redshift (e.g. da Cunha et al. 2013). The CMB heating is calculated by using a correction to the dust temperature (da Cunha et al. 2013, in particular see equation 12), i.e.

$$T_{\text{dust}}^{\text{corr}}(z) = (T_{\text{dust}}^{4+\beta} + T_{\text{CMB}}^{4+\beta}(z) - T_{\text{CMB}}^{4+\beta}(z=0))^{1/(4+\beta)}, \quad (\text{B2d})$$

where $T_{\text{CMB}}(z) = 2.725(1+z)$ is the average CMB temperature at $z$, and $\beta$ is the dust emissivity coefficient (Draine & Lee 1984). Similarly to da Cunha et al. (2013), we assume a fiducial dust emissivity of $\beta = 2$.

Note that dust can be destroyed by sublimation or evaporation. In the adopted chemical network (and the current GRACKLE version), dust abundance is not accounted self-consistently, and we calculate it by assuming that the $M_{\text{dust}}$ is proportional to the metal mass of the gas (equation B2c). To mimic dust destruction, we assume a negligible dust contribution to the optical depth when $T_{\text{dust}} > 2 \times 10^3$ (e.g. Bauer et al. 1997). We neglect dust sputtering driven by SNe (e.g. Valiante et al. 2009; Draine 2011), which might lead to dust destruction, and thus further reduce the radiation pressure efficiency.

Finally, we note that in our radiation feedback model (equations 7a and 7c), energy conservation is guaranteed by construction, if only the UV and ionizing contributions are present. However, the IR trapping is accounted by approximate formulas ($f_{\text{ir}} = \tau_{\text{ir}}$) that are consistently evaluated with the gas properties but do not account for radiative losses in the IR cascade. For this reason, we have added the additional energy conservation check presented in Section 2.4.3. Note that a posteriori we have found that IR trapping is a subdominant contribution to radiation pressure ($\tau_{\text{ir}} \lesssim 10^{-1}$, Section 3.2.2); thus, probably the check is not necessary.

## APPENDIX C: POST-PROCESSING IONIZATION STATE AND EMISSION

Radiative transfer is not followed during the evolution of the simulation. To compute the ionization state of the various atomic species, we post-process the simulation outputs using the photoionization code CLOUDY (Ferland et al. 2013). 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$).

The radiation fields include the UVB intensity at 912 Å (Haardt & Madau 2012), the CMB background and a galactic background that is obtained by rescaling with Dahlia SFR the Galaxy spectrum (Black 1987), in particular $G$, the FUV flux in the Habing band, which is usually normalized to the Galactic mean value $G_0$. Thus, we use $G = 130 G_0$ in our calculation. Note that a larger value of $G$ does not yield a large variation of the expected [C II] in molecular gas (Vallini et al. 2016).

Note that accounting for the UVB is not relevant for the ionization state of the gas in the proximity of galaxies (Gnedin 2010): because of the high-density environment, the gas is efficiently shielded and the galactic emission is the dominant radiation source. Thus, in our CLOUDY models we consider that the gas is shielded by a column density of $N \simeq 10^{20} \text{ cm}^{-2}$.

Spatial variation of the incident radiation is not considered in this work. It is to note that the variation of $G$ is very low in our Galaxy (Habing 1968; Wolfire et al. 2003), i.e. $\langle G \rangle = G_0$ and $((G/G_0)^2 - 1)^{1/2} \simeq 3$. However, in the close proximity of young OB associations, the flux can be very high, e.g. $10^6 G_0$ at 0.1 pc from a single OB star (Hollenbach & Tielens 1999), or – equivalently – $10^6 G_0$ at 10 pc from a starburst with SFR $\simeq M_⊙ \text{ yr}^{-1}$. As the flux propagates with $1/r^2$, such variation is not perceived in most of the volume. Such effect is not considered in this work.

Using the C II density obtained with CLOUDY, we compute the [C II] luminosity by using equation (3) in Vallini et al. (2013). The effect of CMB suppression of [C II] is included in the model (Pallottini et al. 2015; Vallini et al. 2015). Such effect suppresses the emission where the spin temperature of the [C II] transition is close to the CMB one. This is relevant for a low-density ($n \lesssim 10^{-1} \text{ cm}^{-3}$) medium that does not have enough collision to decouple from the CMB. We note that the result from the present [C II] modelling is consistent with the one obtained by directly using CLOUDY to compute the PDR emission.

Note that the current model does not account for the photoevaporation effect on MC, which has an important impact on FIR emission (Vallini et al. 2016), particularly when including a spatially varying FUV field. More detailed modelling will be accounted in a future work.

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