The internal structure of z\simeq 6 galaxies

# Zooming on the internal structure of z≃6 galaxies

## Abstract

We present zoom-in, AMR, high-resolution () simulations of high-redshift () 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 , representative of a typical Lyman Break Galaxy. Dahlia has a total  mass of , that is mainly concentrated in a disk-like structure of effective radius kpc and scale height pc. Frequent mergers drive fresh gas towards the centre of the disk, sustaining a star formation rate per unit area of . The disk is composed by dense (), metal-rich () gas, that is pressure-supported by radiation. We compute the m [C ] emission arising from Dahlia, and find that of the total [C ] luminosity () arises from the  disk. Although of the C  mass is transported out of the disk by outflows, such gas negligibly contributes to [C ] emission, due to its low density () and metallicity (). Dahlia is under-luminous with respect to the local [C ]- relation; however, its luminosity is consistent with upper limits derived for most galaxies.

###### keywords:
galaxies: high-redshift, formation, evolution, ISM – infrared: general – methods: numerical
12

## 1 Introduction

The discovery and characterization of primeval galaxies constitute some of the biggest challenges in current observational and theoretical cosmology3.

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 (). 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 (Bouwens et al., 2015). However – using these surveys broad band alone – little can be learned about other properties as their gas and dust content, metallicity, interactions with the surrounding environment (e.g. Barnes et al., 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 energetics of the interstellar medium (ISM) can be obtained with observations of Far IR (FIR) fine structure lines, and in particular the [C  line at 157.74 m. The [C ] line is the dominant coolant of the ISM being excited in different ISM phases, as the diffuse cold neutral medium (CNM), warm neutral medium (WNM), high density photodissociation regions (PDRs), and – to a lower extent – ionized gas (Tielens & Hollenbach, 1985; Wolfire et al., 1995; Abel, 2006; Vallini et al., 2013). As [C ] 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 ] line is a promising probe as it is often the brightest among FIR emission lines, accounting for up to 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- ISM (e.g. De Looze et al., 2014). The unprecedented sensitivity of the Atacama Large Millimeter/Submillmeter Array (ALMA) makes it possible for the first time to use [C ] emission to characterize high- galaxies. Before the ALMA advent, in fact, detections were limited to a handful of QSO host galaxies, and rare galaxies with extreme SF rates (, e.g. Maiolino et al., 2005; De Breuck et al., 2011; Carilli & Walter, 2013; Gallerani et al., 2012; Cicone et al., 2015).

However, for “normal” star forming galaxies () at early ALMA searches for [C ] lines have mostly yielded upper limits (e.g. Ouchi et al. 2013 Kanekar et al. 2013; Ota et al. 2014; Schaerer et al. 2015). The situation has changed recently with a number of robust [C ] detections (e.g. Maiolino et al. 2015; Capak et al. 2015; Willott et al. 2015; Knudsen et al. 2016).

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

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

• Suppression of [C ] line around star forming regions Vallini et al. (2013), typically observed as a displacement of the [C ]  with respect to the UV emitting region, as seen e.g. in BDF3299 (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 ]-emitting gas.

• Suppression of [C ] line by the increased CMB temperature in the WNM/CNM component (Pallottini et al., 2015; Vallini et al., 2015), similarly to what 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 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 kind 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; Barai et al., 2015; Agertz & Kravtsov, 2015), and radiation pressure (Wise et al., 2012; Ceverino et al., 2014); other analyses focus on implementing complex chemical networks in simulations (Tomassetti et al., 2015; Maio & Tescari, 2015; Bovino et al., 2016; Richings & Schaye, 2016; Grassi et al., 2016), 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 ] 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 ] mapping (Gong et al., 2012; Silva et al., 2015; Yue et al., 2015; Pallottini et al., 2015), which would led to a statistical characterization of the high- galaxy population. In the present work we simulate a galaxy typically detected in [C ] with ALMA current observations.

The paper is structured as follows. In Sec. 2 we detail the numerical model used to set-up the zoom-in simulation, and describe the adopted  star formation prescription (Sec. 2.2), mass and energy inputs from the stellar populations (Sec. 2.3) and feedback (including SN, winds and radiation pressure Sec. 2.4 – see also App. B and App. A). The results are discussed in Sec. 3, where we analyze star formation history and feedback effects in relation to ISM thermodynamics (Sec. 3.2) and its structural properties. The expected [C ] emission and other observational properties of high- galaxies are discussed in Sec. 3.3. Conclusions are given in Sec. 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 Particle Mesh N-body solver for the dark matter (DM) and an unsplit 2nd-order MUSCL4 scheme for the baryons. Gravity is accounted by solving the Poisson equation on the AMR grid via a multi-grid 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.15 (Bryan et al., 2014), the standard library of thermo-chemical processes of the AGORA project (Kim et al., 2014). Via grackle, we follow the H and He primordial network and tabulated metal cooling and photo-heating rates calculated with cloudy (Ferland et al., 2013). Cooling includes also inverse Compton off the cosmic microwave background (CMB), and heating from a redshift-dependent ionizing UV background (Haardt & Madau, 2012, UVB). Since  gas phase formation is not accounted for, we do not include the cooling contribution of such species.

Because of stellar feedback (Sec 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 approach is one of the possible scheme used to solve the over-cooling problem that affect 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 (), work (), and dissipation (Agertz & Kravtsov, 2015). Following Mac Low (1999) we assume a dissipation time scale proportional to the size of the cell (injection scale) and inversely proportional to the Mach number6. Since the dynamical time is essentially set by the free-fall time, the dissipation time can be written as . Then, the non-thermal energy loss due to dissipation can be written as (Teyssier et al., 2013, see eq. 2). As noted in Teyssier et al. (2013), such 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 nested grid using a real-space convolution approach (cf. Bertschinger, 1995). The adopted Lagrangian perturbation theory scheme is perfectly suited to produce IC for multi-mass 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 from to with DM mass resolution of . The resolution of the coarse grid is , and we do not include additional levels of refinement. Using hop (Eisenstein & Hut, 1998) we find the DM halo catalogue at . 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 (resolved by DM particles), whose virial radius is at . Using hop we select the minimum ellipsoid enveloping , and trace it back to . As noted in Oñorbe et al. (2014), this is usually sufficient to avoid contamination7. At the trace back volume is . Using music we recalculate the ICs, by generating 3 additional level of refinement. For such multi-mass set-up, the finer DM resolution is , that corresponds to a spatial resolution of . 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 et al. 2015, where the target halo is chosen at ).

In the zoom-in simulation corresponds to our coarse grid resolution, and we allow for 6 additional refinement levels, based on a Lagrangian mass threshold-based criterion. At , the baryonic component of the selected halo has a mass resolution of and a physical resolution of . For convenience, a summary of the resolution outline can be found in Tab. 1. Note that the refined cell of our simulations have mass and size typical of molecular clouds (MC, e.g. Gorti & Hollenbach, 2002; Federrath & Klessen, 2013).

In the present paper we refer to metallicity () 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 , and has metallicity floor . The metallicity floor mimics the pre-enrichment of the halo at high-, when we do not have the resolution to follow precisely star formation and gas enrichment. We set , a level that is compatible with the metallicity found at high- in cosmological simulations for diffuse enriched gas (Davé et al., 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 contribute for only of the total metal mass produced by stars by in the refined region.

### 2.2 Star formation model

We model star formation (SF) by assuming a  dependent Schmidt-Kennicutt relation (Schmidt, 1959; Kennicutt, 1998)

 ˙ρ⋆=fH2ρ/tsf (1a) where ˙ρ⋆ is the local SF rate (SFR) density, fH2 the molecular hydrogen fraction, ρ the gas density and tsf the SF time scale. In eq. 1a we assume the SF time scale to be proportional to the free-fall time, i.e. tsf=ζ−1sf√3π/(32Gρ), (1b)

where describes the SF efficiency and it is treated as a parameter in the present work (cf. Semenov et al., 2016, see discussion in Sec. 3.2.1). To calculate we adopt the KTM09 model (Krumholz et al., 2008, 2009; McKee & Krumholz, 2010). Such model considers  formation on dust grains by computing radiative transfer on a idealized MC and assumes equilibrium between formation and dissociation rate of . The solution for can be approximated as

 fH2 =(1−0.75s/(1+0.25s))Θ(2−s) (2a) s =ln(1+0.6χ+0.01χ2)/0.6τuv (2b) χ =71(σd,21/R−16.5)((G/G0)/(n/cm−3)), (2c) where Θ is the Heaviside function, τuv the dust optical depth of the cloud, σ−21d=σd/10−21cm−2 is the dust absorption cross section (Li & Draine, 2001), R/10−16.5cm3s−1 is the formation rate coefficient of H2 on dust grains (Wolfire et al., 2008), G is the FUV flux in the Habing band (6−13.6eV) normalized to the average Milky Way (MW) value G0 (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. τuv=10−21cm−2NHZ/Z⊙/μ, where NH is the hydrogen column density and μ the mean molecular weight. In the simulation, the column density is calculated as NH=nlcell; because of the mass threshold-based criterion used as a refinement in AMR, we expect lcell∝n−1/3, thus NH∝n2/3. Note that both σd and R are proportional to the dust mass, that we assume to be proportional to the metallicity. Then the ratio between σd and R is independent of Z. Additionally, eq. 2 can be simplified by assuming pressure equilibrium between the CNM and WNM. In this case, eq. 2c turns out to be independent on G/G0 and can be written as (Krumholz et al., 2009) χ=0.75(1+3.1(Z/Z⊙)0.365). (2d)

As shown in (Krumholz & Gnedin, 2011), for such approximation gives  fractions compatible with those resulting from a full non-equilibrium radiative transfer calculations.

In Fig. 1 we plot 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 ; it steeply rises up to in one density dex and asymptotically reaches . The critical density where the gas can be considered molecular () is roughly inversely proportional to the metallicity, i.e. (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 that the one predicted by the KTM09 model (e.g. Bovino et al., 2016).

Because of the particular shape of the relation, the adopted SF law (eqs. 12) is roughly equivalent to a prescription based on a density threshold criterion:

 ˙ρ⋆=Θ(n−nc)mpn/tsf, (3a) where mp is the proton mass and the critical density nc≃26.45(Z/Z⊙)−0.87cm−3 (3b)

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

Eqs. 3 are not used to calculate the in the simulation. However, being simpler, such formulation can be used to enhance our physical intuition of the adopted SF law8 in analyzing the results. As noted in Hopkins et al. (2013), the morphology of a galaxy is very sensitive to the minimum density of the cells that are able to form star.

During the simulation, eqs. 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 a 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  formation on dust grains, we do not allow gas to form stars if the dust temperature is larger than , because of dust sublimation (see Sec. 2.4 and App. B for the details on the dust prescriptions).

For the present work we assume a SF efficiency , 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 similar once feedback are 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 eqs. 12d) “star particles” to represent stellar populations. To this aim, we adopt a Kroupa (2001) IMF

 Φ(m)∝ [m−α1Θ(m1−m) (4a) + m−α2Θ(m−m1)mα2−α11], where α1=1.3, α2=2.3, m1=0.5M⊙, and m is in the range [10−1−102]M⊙. The proportionality constant is chosen such that ∫100M⊙0.1M⊙mΦdm=1. (4b)

Once formed, stars affect the environment with chemical, mechanical and radiative feedback. These stellar inputs are parameterized by the cumulative fraction of the returned gas mass, metals and energy (e.g. Salvadori et al., 2008; de Bennassuti et al., 2014; Salvadori et al., 2015). Mass and energy inputs are conveniently expressed per unit stellar mass formed ().

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

 R(t⋆)= ∫100M⊙m(t⋆)(m−w)Φdm (5a) Y(t⋆)= ∫100M⊙m(t⋆)mZΦdm, (5b) where w(m,Z⋆) and mZ(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 lifetime9 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. Δt∼10−2Myr vs Δt∼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 AGB stars in young (<40Myr) or evolved stellar populations: ϵsn(t⋆)= ∫40M⊙m(t⋆)>8M⊙esnΦdm, (5c) ϵw(t⋆)= ∫100M⊙m(t⋆)ewΦdm, (5d) where esn=esn(m,Z) and ew=ew(m,Z) are the energy released by SN and stellar winds in units of 1051erg≡1foe; we have further assumed that only stars with 8≤m/M⊙≤40 can explode as SN. Radiative energy inputs can be treated within a similar formalism. The cumulative energy ϵ12 associated to the spectral range (λ1,λ2) can be written as ϵ12(t⋆)= ∫t⋆0∫100M⊙m(t)L12Φdmdt (5e) L12(t)= ∫λ2λ1Lλdλ, (5f)

where is the luminosity per unit wavelength and mass. For convenience, we express the radiation energy in units of , as for the mechanical energy (eqs. 5c). In the following we specify in eq. 5e, by separately considering ionizing radiation (, ) denoted by , and the soft UV band, , defined as the range (, ).

In eqs. 5, the quantities , , , , and can be calculated from stellar evolutionary models. We adopt the padova (Bertelli et al., 1994) stellar tracks for metallicities to compute the chemical10, mechanical and radiative inputs using starburst99 (Leitherer et al., 1999, 2010).

In Fig. 2 we plot , , , , and as a function of . For each curve the shaded regions denote the metallicity range; single tracks are indicated with dark lines. The time interval during which massive stars can explode as SN () is highlighted with vertical dashed lines, and the upper axis is labelled with the corresponding stellar mass.

Note that the OB stars contribution () to , and is roughly proportional to and (see also Agertz et al. 2013, in particular eqs. 4). As in the simulation the metallicity floor is set to , we slightly overestimate the wind contribution for low .

Finally, note that the change of behavior of at is due to the ionizing () photon production suppression. At late times (), AGB stars give a negligible mechanical energy contribution () but return mass and metals to the gas (, ).

### 2.4 Stellar feedback

Eqs. 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 , metallicity and age residing in a gas cell with volume . In our scheme, when the simulation evolves for a time , the chemical feedback act as follows:

 ρ =ρ+[R(t⋆+Δt)−R(t⋆)]M⋆/Vcell (6a) Z =Z+[Y(t⋆+Δt)−Y(t⋆)]M⋆/Vcell, (6b) where ρ and Z are the the gas density and metallicity and R and Y are taken from eqs. 5a. Note that chemical enrichment is due both to the SN and AGB winds.

where and are the thermal and non-thermal energy densities, and and are the fractions of thermal and kinetic energy deposited in the ISM. Thus, accounts for the momentum injection by SN and for the thermal pressure part.

In the present work, we have developed a novel method to compute such quantities. The method derives and from a detailed modelling of the subgrid blastwave evolution produced by the SN explosion. We calculate and by evaluating the shock evolution at time , the time step of the simulation11.

The adopted blastwave 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, and depends on the blastwave evolutionary stage. The latter, in turn depends on the gas density, cooling time, and the initial energy of the blast (, in eq. 5c).

The model details are presented in App. A. As an example, in Fig. 3, we show the energy evolution for a single SN explosion () in a gas characterized by and . The total energy is constant in the Sedov-Taylor (ST) stage, it decrease down to during the shell formation (SF) stage, and it evolves as in the pressure driven snowplow (PDS) stage (see eq. 9). In the ST stage most of the energy is thermal, i.e. ; however, in the SF stage increases, since part of the thermal energy is radiated away and some is converted into kinetic form(e.g. Cox, 1972; Cioffi et al., 1988). Finally, during the PDS stage the ratio of thermal to kinetic is (see eqs. 6.14 in OM88).

In this particular example – a SN exploding in a cell – by assuming a simulation time step of , we find that the blastwave is in the PDS stage, and the gas receives (via eqs. 6c) a fraction of energy and in thermal and kinetic form, respectively. During , about % of the initial SN energy has been either radiated away or lost to work done by the blastwave to exit the cell. The model is in broad agreement with other more specific numerical studies (e.g. Cioffi et al., 1988; Walch & Naab, 2015; Martizzi et al., 2015).

#### Stellar winds

Stellar winds are implemented in a manner paralleling the above scheme for SNe. The energy variation can be calculated via eq. 5c, where is substituted with , given in eqs. 6c. Then, and for winds are calculated via a stage scheme similar to SN. The main difference in the efficiency factors calculation depends on the mode of energy production, i.e. impulsive for SNe, continuous for winds. The complete scheme is detailed in App. 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 a SN explodes in the lower density bubble produced by the stellar progenitor wind, the adiabatic phase lasts longer and consequently and increase considerably.

Finally, we account for radiation pressure from stars. The coupling of the gas with the radiation can be expressed in terms of , the rate of momentum injection (Krumholz & Matzner, 2009; Hopkins et al., 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

 ˙prad= (Lion/c)(1−exp(−τion)) (7a) + (Luv/c)((1−exp(−τuv))+fir), where c is the speed of light, τion the hydrogen optical depth to ionizing radiation, and fir is the term accounting for the IR-trapping. Lion and Luv are calculated by integration of the stellar tracks (eqs 5e). The calculation of τuv is modelled in Sec. 2.2 (eq. 2b and related text). We compute τion and fir according to the physical properties of the gas, as detailed in App. B. Note that we do not assume, as sometimes done, τion∼τuv≫1, i.e. we allow for the possibility that some LyC photons can escape. In smoothed particle hydrodynamics (SPH) codes, radiation pressure (eq. 7a) can be implemented as a “kick” (e.g. Hopkins et al., 2011; Barai et al., 2015). Namely, a velocity Δv=˙pradΔt/mb is directly added to some of the SPH particles of mass mb near the photon source. The particles that receive kicks are statistically chosen according to a probability Pkick, and with kick direction ^v that is sampled from a random distribution. Considering the specific kinetic energy of the SPH particles, we would have Missing \left or extra \right (7b) where v is the original particle velocity, the ⟨⟩ operator indicates the particles sum weighted by the SPH kernel, and Vcell is the kernel volume. Thus, because of the kick, the increase of energy density would be12 Δek =⟨mbPkick(Δvv^v+0.5(Δv)2⟩/Vcell =0.5mb(Δv)2/Vcell (7c) =0.5(˙pradΔt)2/(mbVcell),

where can be calculated via eq. 7a, and eq. 7 can be directly cast into the AMR formalism. Additionally, because of our approximate treatment of IR-trapping (see App. B), we force energy conservation: , 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 , as described in the beginning of Sec. 2.

## 3 Results

At (), the simulated zoom-in region contains a group of 15 DM haloes that host galaxies. We target the most massive halo () that hosts “Dahlia”, which is a galaxy characterized by a stellar mass of , therefore representative of a typical LBG galaxy at that epoch. Dahlia has 14 satellites located within from its centre. The six largest ones have a DM mass in the range , and they host stars with total mass . Additionally, there are eight smaller satellites (), with .

### 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 (), temperature (), total (thermal+kinetic) pressure (), and metallicity () maps13 at .

#### Environment (scale ≃160 kpc)

Dahlia sits at the centre of a cosmic web knot and accretes mass from the intergalactic medium (IGM) mainly via 3 filaments of length , confirming previous findings (Dekel et al., 2009). These overdense filaments () are slightly colder () than the IGM () as a consequence of their shorter radiative cooling time (). Along these cold streams, pockets of shock-heated () 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 of metal-polluted regions. The size of the metal bubbles ranges from kpc in the case of Dahlia to a few for the satellites. Bubble sizes increase with the total stellar mass (see P14, in particular Fig. 13), and age of the galaxy stellar population.

On these scales, the pressure is dominated by the thermal component (); higher values of pressure, associated to non-thermal feedback effects (e.g. gas bulk motion), are confined around star forming regions, again traced by the metallicity distribution.

#### Circumgalactic medium (scale ≃50 kpc)

To investigate the circumgalactic medium (CGM), we zoom in a region within kpc from Dahlia’s centre. On these scales, we can appreciate the presence of several Dahlia’s satellites, i.e. extended (few ) structures that are times denser than the filament in which they reside. Two of these density structures are particularly noticeable. These are located at a distance of 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 smaller () overdense clumps (). The ones within Dahlia’s metal bubble are enriched to . This high value is indicative of in-situ self-pollution, which possibly follows an initial pre-enrichment phase from Dahlia. Clumps outside Dahlia metal bubble have on average an higher density (). Since these clumps are unpolluted, they have not yet formed stars, as the effective density threshold for star formation is (see eq. 3b and Sec. 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 , 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 (), diffuse () and mildly enriched () 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).

#### ISM (scale ≃10 kpc)

The structure of Dahlia’s ISM emerges once we zoom in a region from its centre. In the inner region (), a counterclockwise disk 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 disks in these early systems has already been suggested by other studies. For example, Feng et al. (2015) show that already at nearly of galaxies with have disks (see also Sec. 3.3).

The spiral central region and the spiral arms are dense () and cold (), and the active SF produces a large in-situ enrichment (). 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 . Within spiral arms radiation pressure induced bulk motions largely dominate the total pressure, which reaches values as high as . The imprint of SN shocks is evident in the temperature map in regions with . Shock driven outflows originated in spiral arms travel outward in the CGM, eventually reaching the IGM if outflow velocities exceed the escape velocity (, see Fig. 4 in Gallerani et al. 2016).

Outflows are either preferentially aligned with the galaxy rotation axis, or they start at the edge of the disk. However, when spherically averaged, infall and outflow rates are nearly equal ( at , Gallerani et al. 2016), and the system seems to self-regulate (see also Dekel & Mandelker, 2014).

Outside the disk, clumps with density are also present and are actively producing stars. These isolated star forming MCs are located at a distance from the centre, and show up as spots of high pressure (); some of this MCs are completely disrupted by internal feedback and they can be recognized by the low metallicity (): this is consistent with the outcome of numerical simulations of multiple SN explosions in single MC (e.g. Körtgen et al., 2016).

Fig. 5 shows spherically averaged density, metallicity, and  density profiles for the gas. The density profile rapidly decreases from at to at , and then flattens at larger distances. Such profile is consistent with the average profile of galaxies presented in P14. There we claimed that the density profile is universal once rescaled to the halo virial radius (see also Liang et al. 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 it has already dropped to . Within , the metallicity gradient closely tracks the density profile, while for the decrease is steeper. Pallottini et al. (2014b) find that the metallicity profile is not universal, however it usually extend up to few virial radii, as for Dahlia; further insights can be obtained by analyzing the - relation (Sec. 3.2.3).

In Fig. 5 we note that the gradient found in Dahlia at is slightly steeper than the one inferred from observations of galaxies: i.e. we find while the observed ones are (Wuyts et al., 2016) and (Troncoso et al., 2014). This suggests that the metallicity profile evolve 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 claim we should evolve the simulation to .

The  profile is spiky, and each peak marks the presence of a distinct SF region14. In Dahlia  is mainly concentrated within and it is distributed in the disk-like structure seen in Fig. 4 (see Sec. 3.3). The location of the other peaks correspond to the satellites, which are mostly co-located with metallicity peaks. With increasing metallicity, in fact, lower densities are needed to form  (eq. 3b).

### 3.2 Star formation and feedback history

We analyze the SF history of Dahlia and its major satellites by plotting in Fig. 6 the cumulative stellar mass () and star formation rate () vs. time15.

For the whole galaxy sample, the time averaged ( r.m.s.) specific star formation is . This mean value is comparable to that obtained by previous simulations of high- galaxies (Wise et al., 2012) and broadly in agreement with observations (Stark et al., 2013). At early times the reaches a maximum of , while a minimum of is found during the late time evolution. Both the large range and maximum at early times are consistent with simulations by Shen et al. (2014). At late times, the is in agreement with analytical calculation (Behroozi et al., 2013), and with observations (González et al., 2010), although we note Dahlia has a larger stellar mass with respect to the galaxies in the sample ().

At all times, Dahlia dominates both the stellar mass and star formation rate, whose mean value is . Its stellar mass grows rapidly, and it reaches by (), i.e. after from the first star formation event. Such rapid mass build-up is due to merger-induced SF, that plays a major role at high- (Poole et al., 2016; Behroozi et al., 2013; Salvadori et al., 2010). The is roughly constant from to and reaches a maximum of at . With respect to observations of LBG galaxies (e.g. Stanway et al., 2003; Stark et al., 2009) the and of Dahlia are above the mean values, but still consistent within one sigma. Additionally, the combination of , , and for Dahlia are compatible with the fundamental mass metallicity relation observed in local galaxies (Mannucci et al., 2010).

The total stellar mass in satellites is . Typically, SF starts with a burst, generating of stars during the first . Then the exponentially declines and becomes intermittent with a bursty duty cycle of . This process can be explained as follows. As an halo forms, at its centre the density of the gas slowly rises. When the density is higher than the critical density of  formation (eq. 3b), the gas in the inner region is converted into stars in few free-fall times. Then feedback, and in particularly coherent SN explosions (, see Fig. 2), quenches , and the star formation activity becomes self-regulated. As mergers supply fresh gas, the suddenly goes out of equilibrium and becomes bursty again. Note that self-regulation is possible only for major satellites, since smaller ones () 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 galaxies at (Kauffmann, 2014). Furthermore, in our satellites we find that the typical behavior of the burst phases – starburst - quiescent - post-starburst – is qualitatively similar to what found by Read et al. (2016b), that simulate the evolution of a galaxy for (see also Teyssier et al. 2013; Read et al. 2016a for further specific studies on the bursty nature of this kind of galaxies).

Since individual galaxies are defined as group of star particles in the same DM halo at , 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 panel of Fig. 6 we plot the of individual halos defined by their merger history at . Galaxies with active SF at Myr merge into Dahlia at a later time, thus they do not appear individually in the left panel of Fig. 6.

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

#### Star formation efficiency

represents the quantity of gas converted in stars within a free-fall time (see eq. 1). In Fig. 7 we plot the effective star formation efficiency () as a function of gas density, weighted by the  mass fraction at . Most of the  is contained in the range , and the effective efficiency varies from to . Since , the spread is purely due to the dependence of on density and metallicity (see Fig. 1). Note that by construction , 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 (eq. 1b).

Interestingly, our -based star formation criterion is reminiscent of a density threshold one, as below the efficiency drops abruptly (eqs. 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) (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 Sec. 4.

#### Feedback energy deposition

As discussed in Sec. 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 , 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 App. A). For radiation pressure the efficiency is limited by the relatively small dust optical depths (see also App. 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., history and ISM thermodynamics, as discussed in Sec. 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: . 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 , with , apart from fluctuations and jumps as the one at . 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 eqs. 7 and 7a to write , where is the gas mass accelerated by radiation, and we neglect ionizing radiation. Then, using 5e, we can write . Initially, , thus . Once the gas mass is expelled from the star forming region or converted into stars, . Thus the deposition rate increases faster than the 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 (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 , i.e. after 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. We find that on average for , and at later times, thus confirming our hypothesis.

The energy deposition rates for different feedback types are highly correlated in time (Pearson coefficients ). 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 star, 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 () that cannot retain the gas and sustain a continuous star formation activity.

Finally, we remind that, when compared with observational/analytical constraints, the and of Dahlia are higher then the mean, but still consistent within one sigma. We caution that this might imply a somewhat weak feedback prescription.

#### 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 , , , respectively16.

We focus at and consider the gas in a region within from Dahlia’s centre, essentially the scale of the CGM described in Sec. 3.1.2. This region contains a total gas mass of , and metal mass of (additional data in Tab. 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 , and ; these phases contain , and 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/MC. Additionally only of the considered volume shows , i.e. it has been polluted by stars in the simulation. We note that the EOS in the - plane is fairly consistent with the one found in other high- galaxy simulations (e.g. see Fig. 5 in Wise et al., 2012). Comparison between the EOS in the - and - plane highlights the relative importance of different feedback types.

The rarefied gas is characterized by long cooling times. Thus, once engulfed by shocks, such phase becomes mildly enriched () and remains hot (). The enriched rarefied gas preferentially populates the and region of the phase diagram. However, part of the rarefied gas has . This gas component has a temperature set by the equilibrium between adiabatic cooling and the photo-heating by the UV background; it feeds the accretion onto 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 disk. Typically, such gas has and , thus a thermal pressure is expected. However, the total gas pressure is K (see the - 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 Sec. 3.3).

The diffuse gas acts as an interface between the dense disk gas and the rarefied gas envelope. Diffuse gas is found both in hot () and cold () states. The cold part has a sufficiently high mean metallicity, , 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 2-phase medium shape for pressures around , 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 situation, though, is highly dynamic and does not correspond to a true thermal equilibrium.

A final remark is that by a 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 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 relation depends on the history (for an in-depth analysis see P14). The average relation found is consistent with the results from galaxies (Shen et al., 2014).