Deep into the structure of the first galaxies: SERRA views

TL;DR

This study uses high-resolution cosmological simulations to analyze the formation, structure, and emission properties of early galaxies during the Epoch of Reionization, revealing insights into their ISM, star formation, and emission line characteristics.

Contribution

It provides detailed simulation-based insights into the ISM, star formation, and emission line properties of first galaxies, including the effects of starbursts and metallicity.

Findings

Freesia has a stellar age of 409 Myr and a stellar mass of 4.2×10^9 M_sun.

[OIII] emission is concentrated in star-forming regions, while [CII] is more extended.

Freesia's ionising escape fraction is approximately 2%.

Abstract

We study the formation and evolution of a sample of Lyman Break Galaxies in the Epoch of Reionization by using high-resolution ($\sim 10 \,{\rm pc}$), cosmological zoom-in simulations part of the SERRA suite. In SERRA, we follow the interstellar medium (ISM) thermo-chemical non-equilibrium evolution, and perform on-the-fly radiative transfer of the interstellar radiation field (ISRF). The simulation outputs are post-processed to compute the emission of far infrared lines ([CII], [NII], and [OIII]). At $z=8$, the most massive galaxy, `Freesia', has an age $t_\star \simeq 409\,{\rm Myr}$, stellar mass $M_{\star} \simeq 4.2\times 10^9 {\rm M}_{\odot}$, and a star formation rate ${\rm SFR} \simeq 11.5\,{\rm M}_{\odot}{\rm yr}^{-1}$, due to a recent burst. Freesia has two stellar components (A and B) separated by $\simeq 2.5\, {\rm kpc}$; other 11 galaxies are found within $56.9 \pm 21.6 \, {\rm kpc}$. The mean ISRF in the Habing band is $G = 7.9\, G_0$ and is spatially uniform; in contrast, the ionisation parameter is $U = 2^{+20}_{-2} \times 10^{-3}$, and has a patchy distribution peaked at the location of star-forming sites. The resulting ionising escape fraction from Freesia is $f_{\rm esc}\simeq 2\%$. While [CII] emission is extended (radius 1.54 kpc), [OIII] is concentrated in Freesia-A (0.85 kpc), where the ratio $\Sigma_{\rm [OIII]}/\Sigma_{\rm [CII]} \simeq 10$. As many high-$z$ galaxies, Freesia lies below the local [CII]-SFR relation. We show that this is the general consequence of a starburst phase (pushing the galaxy above the Kennicutt-Schmidt relation) which disrupts/photodissociates the emitting molecular clouds around star-forming sites. Metallicity has a sub-dominant impact on the amplitude of [CII]-SFR deviations.

Tables

Table 1. Table 1: Summary of the FIR emission line properties of Freesia. Emission line maps are given in Fig. 7 , the spectra are plotted in Fig. 9 . Offset and radius are calculated from the location of the emission weighted mean and variance of the emission maps, respectively.

line	$L$	shift	width (${\sigma }_v$)	offset	radius
	$\left[L_{\odot }\right]$	[$\mathrm{km}{s}^{-1}$]	[$\mathrm{km}{s}^{-1}$]	[kpc]	[kpc]
[C $\mathrm{II}$]	$7.73\times {10}^7$	-	93.0	-	1.54
[N $\mathrm{II}$]	$5.33\times {10}^5$	138.1	163.0	0.52	0.50
[O $\mathrm{III}$]	$2.07\times {10}^7$	121.1	163.3	0.65	0.85

Equations

$$R_{H2-\mathrm{dust}}=3\times 10^{-17}n{n}_{\mathrm{H}}(\mathcal{D}/{\mathcal{D}}_{\odot }){\mathrm{cm}}^{-3}{\mathrm{s}}^{-1},$$

$${\dot{e}}_{\mathrm{nth}}=-\frac{e_{\mathrm{nth}}}{t_{\mathrm{diss}}},$$

$$t_{\mathrm{diss}}=9.785\left(\frac{l_{\mathrm{cell}}}{100\mathrm{pc}}\right){\left(\frac{{\sigma }_{\mathrm{turb}}}{10\mathrm{km}{\mathrm{s}}^{-1}}\right)}^{-1}\mathrm{Myr},$$

$${\sigma }_{\mathrm{i}}=\frac{{\int }_{{\nu }_i^{\mathrm{low}}}^{{\nu }_i^{\mathrm{up}}}{\sigma }_{\lambda }{L}_{\mathrm{h}\nu }\mathrm{d}\nu }{{\int }_{{\nu }_i^{\mathrm{low}}}^{{\nu }_i^{\mathrm{up}}}{L}_{\mathrm{h}\nu }\mathrm{d}\nu },$$

$${\dot{\rho }}_{\star }={\zeta }_{\mathrm{sf}}\frac{\mu m_{\mathrm{p}}{n}_{H2}}{t_{\mathrm{ff}}},$$

$$⟨N_{\star }⟩=\frac{\mu m_{\mathrm{p}}{n}_{H2}(\delta l{)}^3}{m_b}\frac{{\zeta }_{\mathrm{sf}}\delta t}{t_{\mathrm{ff}}}.$$

$${\zeta }_{\mathrm{cr}}=3\times 10^{-17}(\mathrm{SFR}/{\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1}){\mathrm{s}}^{-1}.$$

$$N_{\mathrm{H}+}\simeq N_d\log (1+59U\mathcal{D}/{\mathcal{D}}_{\odot }),$$

$$N_d\simeq 1.7\times 10^{21}({\mathcal{D}}_{\odot }/\mathcal{D}){\mathrm{cm}}^{-2}$$

$$\mathrm{d}P\equiv P_s\mathrm{d}s=(2\pi {\sigma }_s^2{)}^{-1/2}{\exp }^{-{\left(\frac{s-s_0}{\sqrt{2}{\sigma }_s}\right)}^2},$$

$$s=\ln (n/n_0),$$

$$s_0=-{\sigma }_s^2/2$$

$${\sigma }_s^2=\ln (1+(\mathcal{M}/2{)}^2),$$

$$d{N}_{cl}=(V/l_J^3)\mathrm{d}P.$$

$$L_k=\int F_k{l}_J^{2}d{N}_{cl}$$

$$M_y=\int x_y{l}_J^{2}d{N}_{cl},$$

$$l_{\mathrm{res}}>N_d/n\simeq 5.5(n/100{\mathrm{cm}}^{-3}{)}^{-1}(Z/{\mathrm{Z}}_{\odot }{)}^{-1}\mathrm{pc}.$$

$$[\mathrm{CII}]/{[\mathrm{CII}]}_{\mathrm{DL}}\equiv {\Sigma }_{[\mathrm{CII}]}/{\Sigma }_{[\mathrm{CII}],\mathrm{DL}},$$

$$\log ({\Sigma }_{[\mathrm{CII}],\mathrm{DL}})=1.075\log ({\dot{\Sigma }}_{\star })+7.51,$$

$${\kappa }_{\mathrm{s}}\equiv {\dot{\Sigma }}_{\star }/{\dot{\Sigma }}_{\star ,\mathrm{KS}}$$

$$\log ({\dot{\Sigma }}_{\star ,\mathrm{KS}})=1.4\log ({\Sigma }_g)-12.0,$$

Full text

††pagerange: Deep into the structure of the first galaxies: SERRA views–Deep into the structure of the first galaxies: SERRA views††pubyear: 2019

1 Introduction

Characterising the interstellar medium (ISM) properties of galaxies in the epoch of the reionisation (EoR) represents a key quest of modern cosmology.

Optical/near infrared (IR) surveys have been fundamental in identifying galaxies in the EoR, and further to give us an overview of their stellar masses, star formation rates and sizes up to redshift $z\sim 10$, well within the EoR (Dunlop, 2013; Madau & Dickinson, 2014; Bouwens et al., 2015; Oesch et al., 2018). In particular, lensing has enabled us to probe the faintest galaxies (Smit et al., 2014; Bouwens et al., 2017; Vanzella et al., 2018), that are likely the main responsible for the reionisation and metal enrichment of the intergalactic medium (Barkana & Loeb, 2001; Ciardi & Ferrara, 2005; Bromm & Yoshida, 2011; Pallottini et al., 2014a; Greig & Mesinger, 2017; Dayal & Ferrara, 2018; Maiolino & Mannucci, 2019).

However, to understand the properties of the ISM of such objects, spectral information is needed. In particular, far infrared (FIR) lines can give a wealth of diagnostics on the thermo-dynamical state of the gas and on the interstellar radiation field (ISRF). As these lines are the main coolants of the ISM (Dalgarno & McCray, 1972; Wolfire et al., 2003), they can be used to trace feedback processes responsible for the evolution of these systems. Additionally, since FIR lines are emitted by ions with low (${\rm C}^{+}$) and high (${\rm O}^{++}$) ionisation potential, their simultaneous detection can constrain the intensity and shape of the IRSF. Finally, detection of CO rotational transitions would help us to constrain the physical properties of molecular clouds (Solomon & Vanden Bout, 2005; Carilli & Walter, 2013), and thus understanding the processes of star formation in galaxies at the EoR.

The advent of the Atacama Large Millimeter/Submillmeter Array (ALMA) has made it possible to access FIR lines from “normal” star forming galaxies ({\rm SFR}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}100{\rm M}_{\odot}\,{\rm yr}^{-1}) in the EoR. In particular [C $\scriptstyle\rm II$] at  $158\mu$, being typically the strongest FIR line (Stacey et al., 1991), is now routinely observed at z\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}6 in both Lyman Alpha Emitters (LAE, Pentericci et al., 2016; Bradac et al., 2017; Matthee et al., 2017; Carniani et al., 2018b; Harikane et al., 2018) and Lyman Break Galaxies (LBG, Maiolino et al., 2015; Willott et al., 2015; Capak et al., 2015; Knudsen et al., 2016; Carniani et al., 2018a). Additionally, [C $\scriptstyle\rm II$] follow-up observations have enabled us to study the galaxy kinematics (Jones et al., 2017; Smit et al., 2018), albeit such observations have not yet the level of maturity as those concentrating on intermediate (z\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}3) redshift (e.g. De Breuck et al., 2014; Leung et al., 2019). Low surface brightness gas outside the target galaxy is possibly a tracer of outflows: probing such material would help us to constrain the feedback mechanism driving the evolution of EoR galaxies. However, so far only statistical evidence of its presence is currently available (Gallerani et al., 2018; Fujimoto et al., 2019). The presence of [O $\scriptstyle\rm III$] at $88\mu$ has been revealed in various observations (Inoue et al., 2016; Laporte et al., 2017; Hashimoto et al., 2018; Tamura et al., 2018); in a few cases both [O $\scriptstyle\rm III$] and [C $\scriptstyle\rm II$] have been simultaneously detected (Carniani et al., 2017; Hashimoto et al., 2018), thus hindering the possibility to constrain the ISRF. Regarding the detection of molecular lines, so far CO has been observed only in one normal star forming galaxy via a serendipitous detection (D’Odorico et al., 2018; Feruglio et al., 2018). Summarising, while a large progress has been made with respect to the first ALMA observation cycles, we currently do not have a complete picture of the FIR properties of these galaxies. It is still unclear whether the local relation between [C $\scriptstyle\rm II$] and star formation rate (De Looze et al., 2014; Herrera-Camus et al., 2015) is fulfilled by high redshift galaxies and which are the physical mechanisms responsible for its larger dispersion w.r.t. the local one (Carniani et al., 2018a). Finally, a convincing explanation of the spectral shifts and spatial offsets between different lines and/or UV continuum that are often observed in these objects is still missing (Capak et al., 2015; Carniani et al., 2017).

To address such issues, on the theoretical side, models of FIR emission from galaxies in the EoR have been developed; these models account for the typically lower metallicity of these systems, higher gas turbulence, and include the suppression of FIR emission by the CMB in low density gas (Vallini et al., 2013; Olsen et al., 2015; Vallini et al., 2015, 2018; Popping et al., 2019). Such models account for the observed ISM and ISRF properties of these objects by post-processing numerical hydrodynamical simulations aimed at describing the formation and evolution of high-redshift galaxies (see Olsen et al., 2018, for an extended discussion).

Cosmological simulations – and in particular zoom-in simulations – have been used in order to study such galaxies. Most works concentrate on the relative importance of different kinds of feedback (e.g. SN, winds from massive stars, radiation pressure) in shaping early galaxy evolution (Agertz & Kravtsov, 2015; Pallottini et al., 2017a; Hopkins et al., 2018b), the chemical evolution of these primeval systems (Maio et al., 2016; Smith et al., 2017; Pallottini et al., 2017b; Lupi et al., 2018; Capelo et al., 2018), the effect of radiation from local sources, the ISM ionisation state, and the consequences for the reionisation process (Katz et al., 2017; Trebitsch et al., 2017; Rosdahl et al., 2018; Hopkins et al., 2018a).

In the past few years we have developed serra, a set of zoom-in simulations of LBGs in the EoR. Starting with Pallottini et al. (2017a), we zoomed-in on the structure of high-z galaxies by studying the formation of few galactic systems and following their evolution down to tens of parsec scales. Then, in Pallottini et al. (2017b), we have analysed the impact of chemistry on the ISM by including thermo-chemical networks to follow the formation of ${\rm{H_{2}}}$, that ultimately led to the formation of stars. Complementing this numerical simulations with both line (Vallini et al., 2015, 2018; Behrens et al., 2019, [C $\scriptstyle\rm II$], CO, Ly$\alpha$) and continuum (Behrens et al., 2018, UV, IR) emission, we have been able to fairly compare our models with high-redshift observations. However, previous simulation were were lacking a consistent modelling of photoevaporation effects due to the ISRF, that can affect the emission properties of the FIR lines (Vallini et al., 2017) and the evolution of molecular clouds (Decataldo et al., 2017).

With the aim of further improving our models, in the present work we include on-the-fly radiative transfer in our hydrodynamical simulations. By also including all the main sources of feedback (radiative, mechanical, chemical), we are able to pinpoint the origin of the deviation from the [C $\scriptstyle\rm II$]-SFR relation that is observed for galaxies in the EoR. Our numerical model is presented in Sec. 2. An overview of the physical properties of our galaxy sample is given in Sec. 3. The FIR emission properties ([C $\scriptstyle\rm II$], [O $\scriptstyle\rm III$], [N $\scriptstyle\rm II$]) are covered in Sec. 4, while Sec. 5 focuses on the [C $\scriptstyle\rm II$]-SFR relation. Conclusions are given in Sec. 6.

2 Numerical simulations

Our simulation suite serra111Greenhouse in Italian. is focused on zooming-in on galaxies in the EoR. In this work, we present “Freesia”, a prototypical LBG galaxy that is hosted by a $M_{\rm h}\simeq 10^{11}{\rm M}_{\odot}$ dark matter (DM) halo at $z=6$. With respect to previous works (Pallottini et al., 2017a, b), here we explore the effect of local sources of radiation.

Gas and DM evolution is simulated with a customised version of the Adaptive Mesh Refinement (AMR) code ramses222https://bitbucket.org/rteyssie/ramses (Teyssier, 2002). In ramses, gas is tracked with a second-order Godunov scheme and particles evolution is computed with a particle-mesh solver (see also Guillet & Teyssier, 2011, for the gravity solver). Radiation coupling to hydrodynamics is performed with ramses-rt (Rosdahl et al., 2013), that solves photons advection within a momentum-based framework with the closure given by setting a M1 condition for the Eddington tensor (Aubert & Teyssier, 2008). Coupling between gas and photons is handled by implementing a non-equilibrium chemical network generated with krome333https://bitbucket.org/tgrassi/krome (Grassi et al., 2014). Metal ion abundances and emission lines are calculated in post-processing, by interpolating grids of models obtained from the photo-ionisation code cloudy V17444https://www.nublado.org (Ferland et al., 2017). The modelling for gas, radiation, stars and line emission is described in Sec.s 2.1, 2.3, 2.2, and 2.4 respectively.

Set-up

We generate cosmological initial conditions (IC)555We assume cosmological parameters compatible with Planck results: $\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 $\rm H_{0}=100\,h\,{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}$ with $h=0.678$, spectral index $n=0.967$, $\sigma_{8}=0.826$ (Planck Collaboration et al., 2014). at $z=100$ with music (Hahn & Abel, 2011). The cosmological volume is $(20\,{\rm Mpc}/h)^{3}$, and the base grid is resolved with a mass $m_{b}=6\times 10^{6}{\rm M}_{\odot}$ per gas resolution element. The Lagrangian volume of the target halo has a linear size of $2.1\,{\rm Mpc}/h$ and is progressively refined by 3 concentric layers with increasing mass resolution, reaching a gas mass of $m_{b}=1.2\times 10^{4}{\rm M}_{\odot}$ around the target halo. In this zoom-in region, we allow for 6 additional level of refinement by adopting a Lagrangian-like criterion. This enables us to reach scales of $l_{\rm res}\simeq 30\,{\rm pc}$ at $z=6$ in the densest regions, i.e. the most refined cells have mass and size typical of Galactic molecular clouds (MC, e.g. Federrath & Klessen, 2013). Note that the resolution and IC are the same used in Pallottini et al. (2017a, b) to allow a fair comparison.

2.1 Hydrodynamics

Chemical network

As in Pallottini et al. (2017b), we implement a non-equilibrium chemical network by using krome (Grassi et al., 2014). The selected network includes H, H*+, H-, He, He+, He++, H2*, H${}_{2}^{+}$ and electrons. The network follows a total of 48 reactions666The reactions, their rates, and corresponding references are listed in App. B of (Bovino et al., 2016): we use reactions from 1 to 31 and 53, 54, from 58 to 61, and from P1 to P9; the rates are reported in Tab. B.1, Tab. B.2, and Tab. 2 of Bovino et al. (2016), respectively., including photo-chemistry, dust processes and cosmic ray-induced reactions (see also Bovino et al., 2016, for the original implementation). Individual ICs for the various species and ions are computed accounting for the chemistry in a primordial Universe (Galli & Palla, 1998).

Metals and dust

Metallicity ($Z$) is tracked as the sum of heavy elements, and we assume solar abundance ratios of the different metal species (Asplund et al., 2009). Dust evolution is not explicitly tracked during simulation. We make the assumption that the dust-to-gas mass ratio scales with metallicity, i.e. $\mathcal{D}=\mathcal{D}_{\odot}(Z/{\rm Z}_{\odot})$, where $\mathcal{D}_{\odot}/{\rm Z}_{\odot}=0.3$ for the Milky Way (MW) (e.g. Hirashita & Ferrara, 2002). While in principle it is possible to incorporate the evolution of dust grains in galaxy formation simulations (e.g. Grassi et al., 2017; McKinnon et al., 2018, see also Asano et al. 2013; De Rossi & Bromm 2017 for semi-analytic models), this would bias the following comparison with Pallottini et al. (2017b) and will be explored in the future. The grain size distribution is important when modelling light extinction, and it is detailed in Sec. 2.2 (see in particular Fig. 2).

Dust provides a formation channel for molecular hydrogen: the formation rate of ${\rm{H_{2}}}$ on dust grains is approximated following Jura (1975):

[TABLE]

where $n$ and $n_{\rm H}$ are the total and Hydrogen gas densities, respectively. We note that for \mathcal{D}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10^{-2}\mathcal{D}_{\odot} the dust channel is dominant with respect to gas-phase formation.

We adopt an initial metallicity floor $Z_{\rm floor}=10^{-3}{\rm Z}_{\odot}$ since at z\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}40 our resolution does not allow us to reach a density high enough for efficient ${\rm{H_{2}}}$ formation in the pristine gas of mini-halos, and consequently recover the formation of first stars (e.g. O’Shea et al., 2015; Smith et al., 2018). Such floor only marginally affects the gas cooling time and it is compatible with the metallicity of diffuse enriched IGM in cosmological metal enrichment simulations (e.g. Pallottini et al., 2014a; Maio & Tescari, 2015; Jaacks et al., 2018).

To summarise, metals and dust are treated as passive scalars and we allow for metal enrichment by supernova (SN) explosions and by winds from massive stars (see Sec. 2.3).

Gas thermo-dynamics

We model both the evolution of thermal and turbulent energy content of the gas.

The thermal energy is evolved by the thermo-chemical framework set with krome (see Pallottini et al., 2017b, for details). Note that photo-chemical reaction rates in each gas cell depend on the local amount of radiation and its energy distribution (see Sec. 2.2). Since metal species are not followed individually, we use the equilibrium metal line cooling function calculated via cloudy (Ferland et al., 2013) with a Haardt & Madau (2012) UV background. Following cooling from individual metal species can change the thermodynamics of the low density ISM, but does not appreciably affect the star forming regions, as shown in Capelo et al. (2018, see also ). While such change in the thermodynamical state of the gas can be important to correctly compute emission lines, we recall that in the present work this is accounted for in post-processing (Sec. 2.4). Dust cooling is not explicitly included, however we note that it gives only a minor contribution to the gas temperature for $n<10^{4}{\rm cm}^{-3}$ (e.g. Bovino et al., 2016, in particular see their Fig. 3). We consider the contribution of cosmic microwave background (CMB), that effectively sets a temperature floor for the gas.

We model the turbulent energy content of the gas similarly to Agertz & Kravtsov (2015) (see also Pallottini et al., 2017a): turbulent (or non-thermal) energy density $e_{\rm nth}$ is injected in the gas by SN, winds and radiation pressure, and it is dissipated as (Teyssier et al., 2013, see eq. 2)

[TABLE]

where $t_{\rm diss}$ is the dissipation time scale, which can be written as in Mac Low (1999)

[TABLE]

where $\sigma_{\rm turb}=\sqrt{e_{\rm nth}}$ is the turbulent velocity dispersion. Note that we do not explicitly include a source term due to shear (see Maier et al., 2009; Scannapieco & Brüggen, 2010; Iapichino et al., 2017, for more refined turbulence models).

2.2 Radiation

In ramses-rt, photons are treated as a fluid that is spatially tracked by sharing the same AMR structure of the gas. Photons are separated in different energy bins, each one tracking an independent “fluid”.

For the present work we select 5 photon bins to cover both the energy range of the Galactic UV ISRF (Black, 1987; Draine, 1978) adopted in Pallottini et al. (2017b) and H ionising radiation777Photo-ionisation of He and He+ is not included, as the stellar SED are typically not hard enough to produce such photons; He and He+ ionisation is only due to collision in the simulation.. Fig. 1 shows the stellar energy distribution (SED) per unit mass and unit energy for a stellar population at different age (see Sec. 2.3 for details on the assumptions). The first two low energy bins cover the Habing (1968) band888In this paper, the Habing flux $G$ is indicated in unit of $G_{0}=1.6\times 10^{-3}{\rm erg}\,{\rm cm}^{-2}\,{\rm s}^{-1}$, the MW value. $6.0-13.6\,\rm eV$, which is fundamental in regulating the temperature of the ISM and photo-dissociation regions (PDR). The second bin included in the Habing band is specific for the Lyman-Werner radiation ($11.2-13.6\,\rm eV$), which photo-dissociates ${\rm{H_{2}}}$ via the two-step Solomon process (Stecher & Williams, 1967). The last three bins cover the H-ionizing photons up to the first ionisation level of He ($13.6-24.59\,\rm eV$). For H-ionizing photons, the energy width is chosen such that the bins have the same number of photons when the SED is averaged on a fiducial stellar population. Our fiducial stellar population has an age $10\,\rm Myr$ and $Z_{\star}={\rm Z}_{\odot}$: these stars are the main sources of the ISRF, since such young stars dominate the spectrum of a galaxy with an exponentially rising SFR, as expected for Freesia (Pallottini et al., 2017a, b). To calculate the SED-averaged quantities, we fix the SED to the fiducial one. Note that using three energy bins for H ionisation allows us to reasonably capture the temperature evolution of ${\rm H}^{+}$ regions, since photo-ionisation coupling with the gas is computed using the mean energy within each bin.

In ramses the time-step is determined by the Courant condition (Courant et al., 1928), i.e. $\delta t\sim l_{\rm cell}/v_{g}$, with $v_{g}$ being the gas velocity; taking $v_{g}$ of the order of the rotational velocity of the galaxy, this yields $\delta t\sim 30\,{\rm pc}/100\,{\rm km}\,{\rm s}^{-1}\sim 0.3\,\rm Myr$. ramses-rt adopts an explicit scheme for the time evolution of radiation. Since hydrodynamics and radiation are coupled, it follows that the time-step of the simulation is determined by the minimum between the sound and light crossing time. To limit the computational load, we consider the reduced speed of light approximation to propagate wave-fronts (Gnedin & Abel, 2001), adopting $c_{red}=10^{-2}c$. With such prescription, we expect $\delta t\sim l_{\rm cell}/c_{red}\sim 0.01\,\rm Myr$. Using a reduced speed of light approach yields artefacts when propagating light fronts very far from the sources, e.g. in IGM reionisation studies999Note that a more correct treatment on large scales is possible with a variable speed of light approach (Katz et al., 2017).; however, it well captures the radiation transfer in the ISM/CGM of galaxies (see Deparis et al., 2019, for a detailed study of the impact of a reduced speed of light), as such it is well suited for the present work.

Coupling with gas and dust

In the original implementation of ramses-rt, the thermo-chemical time step is performed simultaneously with the radiation propagation and absorption of photons by gas and dust. Such coupled step is sub-cycled in order to ensure simultaneous convergence for the absorbed photons, final ionisation state, and gas temperatures. This scheme is similar to the one adopted by Nickerson et al. (2018), which includes ${\rm{H_{2}}}$ formation in ramses-rt, albeit ${\rm H}^{-}$ and ${\rm H}_{2}^{+}$ are not explicitly followed.

Here we split the convergence steps: similarly to ramses-rt the absorption of photons is sub-cycled; for each of these sub-cycles, we obtain the ionisation state and gas temperatures with krome, that adaptively solves the thermo-chemical time evolution assuming a constant impinging flux.

For the gas absorption, photon cross sections are the same ones used in the chemical network (see Sec. 2.1). For dust we assume a MW-like grain composition from Weingartner & Draine (2001). Both for gas and dust, the cross sections $\sigma_{\lambda}$ are used to pre-compute the $i$-th cross section in the photon energy bin ${\rm h}\nu_{i}^{\rm low}-{\rm h}\nu_{i}^{\rm up}$ as

[TABLE]

i.e. flux-averaged cross sections, with a weight given by the selected impinging SED $L_{\rm h\nu}$. In Fig. 2 we plot the dust absorption cross section ($\sigma_{\lambda}$) per unit of $\mathcal{D}/\mathcal{D}_{\odot}$ for both the MW and Small Magellanic Cloud (SMC) -like dust composition, as a function of photon energy ($h\nu$). For both dust types, the $\sigma_{i}$ are overplotted with dashed lines. The difference in absorption between MW- and SMC-like is \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}1\%, except in the $6.0-11.2\,\rm eV$ band, where the difference is about $\simeq 30\%$.

The analysis of current data seems to favour a MW-like distribution for $\mbox{high-}z$  galaxies. Behrens et al. (2018) manage to explain the observation by Laporte et al. (2017) with a low amount of MW-like dust, that leads to a warm FIR SED. However, the situation is still unclear. De Rossi & Bromm (2017) analyses the role of silicate rich dust, that is expected from a Pop III dominated enrichment: the work shows that a warmer FIR SED can naturally arise because of the emission properties of silicates, that typically emit at higher wavelength with respect to the carbonaceous grains. Further, in De Rossi et al. (2018) it is shown that a combination of silicate and small amount of carbonaceous grains can reproduce the SED observed in Haro 11, a local low-metallicity starburst, thought to be an analogue of high-z galaxies: this entails that assuming a simplified dust models can possibly modify the properties inferred from high-z observations (cfr. Behrens et al., 2018). While dust composition (silicate vs carbonaceous grains) does change the FIR emission properties, in the present work dust is considered only with respect to its absorption properties, that are mainly dependent on the grain size distribution which is assumed to be time-independent (see Sec. 2.1). The analysis of the possible modification to the FIR emission due to a different dust composition is left for a future work.

Summarizing, while a different assumption on the dust distribution can in principle heavily affect the observed SED, it should produce only minor differences in the adopted model, since the only energy bin where $\sigma_{i}$ is different in the two cases is responsible for neither ${\rm{H_{2}}}$ dissociation nor H ionisation. Finally, the self-shielding of ${\rm{H_{2}}}$ from photo-dissociation is accounted for by using the Richings et al. (2014) prescription, given the ${\rm{H_{2}}}$ column density, temperature and turbulence (cfr. with Wolcott-Green et al., 2011), as detailed in Pallottini et al. (2017b).

While our scheme is different with respect to the one presented in Rosdahl et al. (2013) and Nickerson et al. (2018), the overall results are consistent; detailed tests of our adopted scheme are found in the Appendix of Decataldo et al., in prep. (2019), using PDR (Röllig et al., 2007; Nickerson et al., 2018) and H $\scriptstyle\rm II$ region (Iliev et al., 2009) benchmarks.

2.3 Stars

Formation

As in Pallottini et al. (2017b), stars form according to a Kennicutt-Schmidt-like relation (Schmidt, 1959; Kennicutt, 1998) that depends on the molecular hydrogen density ($n_{\rm H2}$):

[TABLE]

where $\dot{\rho}_{\star}$ is the local star formation rate density, $\zeta_{\rm sf}$ the star formation efficiency, $m_{\rm p}$ the proton mass, $\mu$ the mean molecular weight, and $t_{\rm ff}$ the free-fall time. The efficiency is set to $\zeta_{\rm sf}=10\%$, by adopting the average value observed for MCs (Murray, 2011, see also Agertz et al. 2013), while $n_{\rm H2}$ computation is included in the non-equilibrium chemical network. As done in Rasera & Teyssier (2006); Dubois & Teyssier (2008); Pallottini et al. (2014a), eq. 5 is solved stochastically at each time step $\delta t$ in each cell with size $\delta l$, by forming in each possible event a new star particle with mass $m_{\star}=m_{b}\,N_{\star}$, with $N_{\star}$ drawn from a Poisson distribution characterised by mean

[TABLE]

For numerical stability, no more than half of the cell mass is allowed to turn into a star particle. Additionally, we allow only star formation events that spawn stellar clusters with mass $m_{\star}\geq 1.2\times 10^{4}{\rm M}_{\odot}$, i.e. the gas mass resolution of the simulation.

Stellar populations

A single star particle in our simulations can be considered a stellar cluster, with metallicity $Z_{\star}$ set equal to that of the parent cell. For the stellar cluster, we assume a Kroupa (2001) initial mass function and, by using starburst99 (Leitherer et al., 1999), we adopt single population stellar evolutionary tracks given by the padova (Bertelli et al., 1994) library, that covers the $0.02\leq Z_{\star}/{\rm Z}_{\odot}\leq 1$ metallicity range. The stellar tracks are then used to calculate mechanical, chemical, and radiative feedback.

Mechanical and chemical feedback

As in Pallottini et al. (2017a), we account for stellar energy inputs and chemical yields that depend both on metallicity $Z_{\star}$ and age $t_{\star}$ of the stellar cluster. Stellar feedback includes SNe, winds from massive stars, and radiation pressure (see also Agertz et al., 2013).

Depending on the kind of feedback, stellar energy input can be both thermal and kinetic, and we account for the dissipation of energy in MCs for SN blastwaves (Ostriker & McKee, 1988) and OB/AGB stellar winds (Weaver et al., 1977), as detailed in Sec. 2.4 and Appendix A of Pallottini et al. (2017a).

In Pallottini et al. (2017a) radiation pressure was implemented by adding a source term to the turbulent (non-thermal) energy. Thus, to avoid double counting of such feedback, we turn off the original radiation pressure coupling of ramses-rt, that is done by following an extra infrared energy bin (see Rosdahl & Teyssier, 2015).

Radiative feedback

Stellar tracks are used to calculate photon production. At each time step, stars act as a source, dumping photons in each energy bin according to their the stellar age and metallicity (Sec. 2.2, in particular Fig. 1), then photons are advected and absorbed in the radiation step, contributing at the same time to the photo-chemistry.

We neglect the cosmic UVB, since the typical ISM densities are sufficiently large to ensure an efficient self-shielding (e.g. Gnedin, 2010). For example, Rahmati et al. (2013) have shown that at $z\simeq 5$ the hydrogen ionisation due to the UVB is negligible for n\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10^{-2}{\rm cm}^{-3}, the typical density of diffuse ISM.

We do not explicitly consider production of radiation from recombination, i.e. we assume that recombination photons are absorbed “on the spot”, which is a valid approximation in the optically thick regime (Rosdahl et al., 2013).

Cosmic-ray (CR) processes are not explicitly tracked during the simulation (see however Dubois & Commerçon, 2016; Pfrommer et al., 2017, for possible implementations). Similarly to Pallottini et al. (2017b), we assume a CR hydrogen ionisation rate proportional to the global SFR (Valle et al., 2002) and normalised to the MW value (Webber, 1998, see Ivlev et al. 2015 for the spectral dependence):

[TABLE]

Coulomb heating is accounted for by assuming that every CR ionisation releases an energy of $20$ eV (see Glassgold et al., 2012, for a more accurate treatment).

2.4 Ions and emission lines

We model metal ion abundances and line emissions using cloudy V17 (Ferland et al., 2017) in post-processing. However, there are some typical challenges and shortcomings to consider when combining emission line codes with simulations (see Olsen et al., 2018, for an overview).

In the following we elaborate on the fact that i) a direct post-processing is computationally unfeasible and not completely consistent, ii) resolution limits our possibility to recover the ionising radiation and internal structure of molecular clouds. We conclude with the description of the solution adopted for the present work.

Consistency of the post-processing

In the simulation, we solve non-equilibrium photo-chemistry by coupling ramses-rt and krome, while cloudy calculations assume photo-ionisation equilibrium. Moreover, cloudy does not account for dynamical effects – such as e.g. shocks – which might affect ion abundances and emission line intensities (cfr. Sutherland & Dopita, 2017; Sutherland et al., 2018). However, the typical PDR code (Röllig et al., 2007) includes – and models more accurately – a larger number of physical processes with respect to the ones typically considered in hydrodynamical simulations with radiative transfer.

One option for the post-processing would consist in using the ISRF resulting from the simulation, and apply a single cloudy model to each cell, given its ISM physical characteristics. On the one hand this is costly: a cloudy model is completed to convergence in $\sim 0.1-0.5$ CPU hours, depending on the chosen maximum column density. Since each snapshot typically contains $\sim 10^{7}-10^{8}$ cells, the cost in CPU hours to post-process a simulation snapshot would be comparable to the cost of the simulation itself (see also Katz et al., 2019, for a machine learning approach to the problem). On the other hand, it is not guaranteed that such approach would result in a more consistent result. For instance in the hydrodynamical simulation we adopt 5 spectral energy bins, which is a very sparse sampling of the SED when compared to photo-ionisation code calculations, where typically thousands of energy bins are included. Moreover, in a cloudy calculation, a cell is divided in optically thin slices, and the temperature and chemical structure is then calculated as a function of the optical depth. Also, cloudy integrates up to photo-ionisation equilibrium to resolve the internal structure of gas patches (single cells in the simulation). Thus, differences in ion abundances are expected with respect to the adopted krome scheme.

Limits given by the resolution

The resolution and refinement criterion of the hydrodinamical simulation do no guaranteed to resolve dense ${\rm H}^{+}$ regions. The column density of ${\rm H}^{+}$ in a slab of a dusty gas can be written as (Ferrara et al., in prep., 2019)

[TABLE]

is the column density due to dust where the optical depth to UV photons becomes unity (see Fig. 2). From Pallottini et al. (2017b), we expect the dense ISM of our galaxy to have $n=3\times 10^{2}{\rm cm}^{-3}$, $Z=0.5\,{\rm Z}_{\odot}$, and $N\simeq 10^{22}{\rm cm}^{-2}$, because of our $\simeq 10\,\rm pc$ resolution. Using these values and eq.s 8, in a typical photo-ionisation region ($U\simeq 10^{-2}$) we obtain a ionisation fraction $f_{\rm H+}\simeq 10\%$. A partial ionisation in a single cell implies that all the ionising photons are absorbed, since in ramses-rt photons are advected after the absorption step101010This is equivalent to state that ramses-rt keeps track of the absorbed flux and not the impinging one; the difference between is negligible only for optically thin cells. Thus, it is possible to find young star clusters embedded in dense gas patches that have $f_{\rm H+}>0$ and $U=0$.

Moreover, with our $\simeq 10\,{\rm pc}$ resolution we cannot resolve the internal structure of molecular clouds (MC), that are made of clumps and cores of size \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}0.1\rm pc. Accounting for such contribution is important to correctly compute the emission in high density ($n\sim 10^{2}{\rm cm}^{-3}$) regions illuminated by a strong ($G\sim 10^{3}G_{0}$) and ionising ($U\sim 10^{-2}$) ISRF (Vallini et al., 2017). These ISM regions are expected to be the main contributors of various FIR lines (i.e. [O $\scriptstyle\rm III$]) in $\mbox{high-}z$  galaxies (e.g. Carniani et al., 2017) and their lower $z$ analogues (e.g. Cormier et al., 2012).

Adopted model

To summarise, the limited resolution of a typical galaxy simulation does not allow us to recover the physical ion structure/line emission even if we would run a single cloudy model per cell, which 1) is computationally very expensive and 2) the assumptions are not completely consistent with the one adopted in the run. To overcome these problems we have adopted a different strategy, as described below.

Two distinct grids of models are calculated, with and without ionising radiation111111In practical terms, the cloudy models without ionising radiation are calculated by interposing a dust-free obscuring screen with column density of $N_{H}=10^{23}{\rm cm}^{-2}$ between the source and the gas.. Parameters for each grid are $n$, $G$, and $Z$, with the following ranges: $10^{-1}\leq n/{\rm cm}^{-3}\leq 10^{4.5}$, $10^{-1}\leq G/{\rm G}_{0}\leq 10^{4.5}$, $10^{-3}\leq Z/{\rm Z}_{\odot}\leq 10^{0.5}$; for each parameter the grid spacing is $0.5$ dex, thus there are a total of 1152 individual models per grid. As assumed in the simulation, dust is proportional to metallicity. We use an impinging SED taken from starburst99 (Leitherer et al., 1999) with age $t_{\star}=10\,\rm Myr$, metallicity $Z_{\star}={\rm Z}_{\odot}$, and rescaled with $G$. Additionally to the ISRF, we include the CMB at the appropriate redshift. Note that cloudy V17 explicitly considers the CMB suppression (da Cunha et al., 2013; Vallini et al., 2015; Pallottini et al., 2015) and does subtract isotropic backgrounds, similarly to what is done in an ALMA observation (Ferland et al., 2017). For each model we stop each calculation at $N=10^{23}{\rm cm}^{-2}$, after convergence has been reached.

Given $n$, $G$, $Z$, and $N$ in each cell, the ion abundances and emission lines can be interpolated from the values found in the computed grid. The grid that includes ionising radiation is selected for those gas patches that either have a ionisation parameter $U>U_{\rm th}\equiv 10^{-4}$ or contain young ($t_{\star}<10\,\rm Myr$) star clusters. The grid with non ionising radiation is chosen for all the other cells. For lines arising from high-ionisation state (i.e. [O $\scriptstyle\rm III$]), this method allows us to recover both the emission from the diffuse ionised medium and from possibly unresolved dense ionised regions. Changes in the selected $U_{\rm th}$ threshold do not yield large variations in the [O $\scriptstyle\rm III$] total luminosity, since (i) low $U$ entails low flux, while high radiation fields are needed to produce a substantial emission (G\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10^{2}G_{0}, see Vallini et al., 2017), (ii) regions containing young star clusters dominate the FIR emission of highly-ionised species (Cormier et al., 2012); this point is detailed in the results, i.e. Sec. 4 (in particular see Fig. 8).

To account for the internal structure of MCs, we use a model similar to Vallini et al. (2015) (see also Vallini et al., 2017, 2018). We assume that a MC with mean density $n_{0}$ and mach number $\mathcal{M}$ is characterised by a probability density function (PDF) given by a log-normal distribution (Padoan & Nordlund, 2011):

[TABLE]

In each MC, we assume that individual clumps have size given by the Jeans length $l_{J}=l_{J}(n,T)$, and volume $l_{J}^{3}$. Then, the differential number of clumps of a MC with total volume $V$ can be written as

[TABLE]

Case example

It is interesting to discuss a case example of such model for [C $\scriptstyle\rm II$] by looking at the efficiency $\eta=L_{\rm CII}/M$, i.e. the [C $\scriptstyle\rm II$] luminosity to total gas mass ratio of the MC.

We start by analysing single cloud models – i.e. without internal structure – for different $n_{0}$ with fixed metallicity ($Z=0.5\,{\rm Z}_{\odot}$) and a non-ionizing radiation field of intensity $G=10^{2}G_{0}$. The value of $\eta$ as a function of $N$ is plotted in the left panel of Fig. 3. At fixed $n_{0}$ for $N<N_{d}$, roughly $F_{\rm CII}\propto N=n_{0}l_{0}$, thus $\eta=L/M\propto F_{\rm CII}l_{0}^{2}/l_{0}^{3}$ is constant, with $l_{0}$ being the size of the cloud; for $N>N_{d}$, $\eta$ drops in all cases, because $F_{\rm CII}$ becomes constant, and thus $\eta\propto 1/l_{0}$.

An important corollary of this result is that, applying these arguments to a simulation with a fixed resolution $l_{\rm res}$ will likely underestimate [C $\scriptstyle\rm II$] emission for $N>N_{d}$, i.e. in cells with resolution (see eq. 8b)

[TABLE]

For typical current state-of-the-art galaxy simulations such value is very demanding, rarely reached, and not compatible with typical refinement criterion selected. Subgrid-models are needed to correctly recover and predict the emission arising from the internal structure of MCs.

To this aim, we use our model for the internal structure of MC (eq.s 9) and recompute the expected [C $\scriptstyle\rm II$] emission (eq. 10). The result is presented in the right panel of Fig. 3 as a function of $n$ and $G_{0}$ for $\mathcal{M}=10$ and $Z=0.5\,{\rm Z}_{\odot}$, i.e. the typical figures found in our simulated galaxies (Pallottini et al., 2017b; Vallini et al., 2018). A maximum value, $\eta\simeq 10^{-1}{\rm L}_{\odot}/{\rm M}_{\odot}$, fall approximately at $n\simeq 10^{3}{\rm cm}^{-3}$ and G\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}5\times 10^{2}G_{0}, while at high density ($n>10^{2}{\rm cm}^{-3}$) and low radiation field (G\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10G_{0}) the [C $\scriptstyle\rm II$] emission is inefficient (\eta\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10^{-2}{\rm L}_{\odot}/{\rm M}_{\odot}). Note that as at low densities (n\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10\,{\rm cm}^{-3}) MCs have little internal structure, the sub-grid model result coincide with single cloud ones.

Adopting such sub-grid model is not completely self-consistent, since cloudy uses an equilibrium approximation, different binning for ISRF, and more physical processes/chemical species. However, it heals some of the problems affecting the calculation of line emission from post-processing of galaxy simulations that cannot spatially resolve ${\rm H}^{+}$ regions and the internal structure of MCs.

3 Overview of the results

3.1 Galaxy formation histories

The formation history of the sample of galaxies in our simulation is plotted in Fig. 4, where we show the stellar mass build up ($M_{\star}$) and the star formation history (SFR, left panel) as a function of age ($t_{\star}$) and redshift ($z$, upper axis).

Freesia is the most massive galaxy in the sample. At $z=8$, it is hosted by a halo of mass $M_{\rm h}\simeq 10^{11}{\rm M}_{\odot}$; its age is $t_{\star}\simeq 409\,\rm Myr$, it has a stellar mass $M_{\star}\simeq 4.2\times 10^{9}{\rm M}_{\odot}$, and an instantaneous ${\rm SFR}\simeq(11.5\pm 1.8)\,{\rm M}_{\odot}\,{\rm yr}^{-1}$, where the error is given by the variance in the last $10\,\rm Myr$. The star formation shows variations on timescales of $\simeq 30\,\rm Myr$, peaking up to ${\rm SFR}\simeq 30\,{\rm M}_{\odot}\,{\rm yr}^{-1}$; overall the evolution is similar to our previous simulation without radiative transfer (Pallottini et al., 2017b), with stellar mass differences of the order $5\%$, and variations in the SFR mostly due to the stochasticity of the star formation prescription (Sec. 2.3, see eq. 6).

Along with Freesia, there is a sample of eleven more galaxies in the simulated region. They have distance from Freesia that has mean (variance) of $56.9\,{\rm kpc}$ ($21.6\,{\rm kpc}$); as they are are at least two virial radii away ($r_{\rm vir}\sim 12\,{\rm kpc}$) at this redshift, we do not label them as satellites. Two of these galaxies are relatively massive, with $M_{\star}\simeq 5\times 10^{8}{\rm M}_{\odot}$ and ${\rm SFR}\simeq 5{\rm M}_{\odot}\,{\rm yr}^{-1}$, and they are younger ($t_{\star}\simeq 300\,\rm Myr$) than Freesia. The other nine are smaller (5\times 10^{7}{\rm M}_{\odot}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}M_{\star}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10^{8}{\rm M}_{\odot}), with lower star formation rates ({\rm SFR}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}{\rm M}_{\odot}\,{\rm yr}^{-1}) and they are typically much younger (t_{\star}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}150\,\rm Myr). Such small galaxies are hosted in M_{\rm dm}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10^{9} dark matter haloes: their $M_{\rm dm}-M_{\star}$ relation has a large scatter, that is compatible with results from other theoretical works, which consider larger sample of galaxies (Xu et al., 2016) and zoom-in simulations focusing on smaller galaxies evolved with higher mass resolution (Jeon et al., 2015; Jeon & Bromm, 2019).

As we had already noted in Pallottini et al. (2017b), most of these objects were not present in our previous simulations without radiative transfer, as their star formation was suppressed; this was entailed by our assumption on the spatially uniform ISRF. In the present work, the ${\rm{H_{2}}}$ formation and hence star formation is not suppressed in objects that are located sufficiently far away from Freesia, as a consequence of flux dilution and attenuation.

In the remaining part of the Sec. along with Sec. 4 we focus on the properties of Freesia at $z=8$. In Sec. 5, the other systems are considered in the analysis.

3.2 Freesia structural properties

A face-on representation of the key structural properties of Freesia is shown in Fig. 5. Freesia has two stellar components – “A” and “B” – separated by $\simeq 4\,{\rm kpc}$, with Freesia-A containing about $\simeq 85\%$ of the total stellar mass and dominating the star formation rate ($90\%$). Both components are highly concentrated, with effective radius of about $\sim 200\,\rm pc$ in both cases; they show stellar surface density peaks with $\Sigma_{\star}\simeq 5\times 10^{11}{\rm M}_{\odot}\,{\rm kpc}^{-2}$, that are surrounded by a stellar halo with low surface density ($\Sigma_{\star}\simeq 10^{7}{\rm M}_{\odot}\,{\rm kpc}^{-2}$), that is likely due to the tidal interaction of the components. Note that some of these stellar clusters have formed recently; while they give a negligible contribution to the total SFR, they can be important in the ionising photon budget.

Looking at the gas density, Freesia-A has a spiral structure with arms characterised by a gas density 10^{2}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}n/{\rm cm}^{-3}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10^{3}; Freesia-B reaches similar densities, but it has a more uniform disk structure, because of its lower mass prevent the development of arms. The only other dense ($n\simeq 10^{2}{\rm cm}^{-3}$) structure is likely an unstable filament located $\simeq 2.5\,{\rm kpc}$ north-west of Freesia-A. These three components are embedded in a lower density medium ($n\simeq 5\,{\rm cm}^{-3}$), with very low density ($\simeq 10^{-2}{\rm cm}^{-3}$) shock-heated patches of gas121212Note that $n\sim 10^{-2}{\rm cm}^{-3}$ is a very low ISM density, but it is higher than the $\Delta=10^{2}$ baryon overdensity that is usually selected to mark the edge of halos in cosmological simulations..

The average131313Average values from the maps are typically quoted in the form ${\rm mean}\pm\,{\rm variance}$. radiation field is $G=(7.9\pm 23.1)\,G_{0}$ i.e. compatible with the assumption in Pallottini et al. (2017b), where $G=G_{0}({\rm SFR}/{\rm M}_{\odot}\,{\rm yr}^{-1}$) (see also Behrens et al., 2018). Note that the variance of the radiation field is $3$ times the mean, as in the MW (Habing, 1968; Wolfire et al., 2003). In analogy with $\Sigma_{\star}$ and $\dot{\Sigma}_{\star}$ peaks, the Habing field has two maxima located at the centre of Freesia-A and Freesia-B and decreases radially from these locations because of flux dilution, as well as gas and dust absorption. While the spatial variance of the Habing field is small, various peaks ($\simeq 2\times 10^{3}G_{0}$) are found in correspondence of stellar clusters, particularly in regions of recent star formation; pockets of gas with an enhanced local radiation can give an important contribution to line emission (Sec. 2.4, see Fig. 3).

The ionisation field has an asymmetric structure and shows a larger variation, with a mean $U=(2\times 10^{-3}\pm 2\times 10^{-2})$; high values (U\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10^{-1}) are co-located with recent star formation events, in particular in the region $2\,{\rm kpc}$ east of Freesia-A, where the gas density is low (n\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}5\,{\rm cm}^{-3}). In the same region we can see that there is a trail of ionising photons that is leaving the galaxy with a conical shape.

The temperature map shows that the spiral arms of Freesia-A and the dense gas in Freesia-B are cold (T\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}250\,\rm K) structures, blistered with warm ($T\sim 10^{4}\rm K$) spots, marking the presence of local radiative feedback. Shock heated regions ($T\sim 10^{5}\rm K$) that reach out from the two stellar components are caused by SN explosions, while the one west of Freesia-A is caused by accretion shocks.

The effect of local radiative feedback is more evident in the molecular and ionised hydrogen maps. The molecular material is concentrated in the spiral arms of Freesia-A and at the location of Freesia-B, with typical surface density peaks with $\Sigma_{\rm H2}\sim 10^{7.5}{\rm M}_{\odot}\,{\rm kpc}^{-2}$ and sizes about $\simeq 100-400\,\rm pc$141414See Leung et al., in prep. (2019) for a more complete analysis of individual MC properties found in our serra simulations.. ${\rm H}^{+}$ regions also show similar values in the peaks of the distribution, i.e. $\sigma_{\rm H+}\sim 10^{7.5}{\rm M}_{\odot}\,{\rm kpc}^{-2}$, but also enclose Freesia with a low surface density $\Sigma_{\rm H+}\sim 10^{6.5}{\rm M}_{\odot}\,{\rm kpc}^{-2}$ halo component. In Freesia-B the correspondence of ${\rm H}^{+}$ regions with spots of warm gas and local ionising field ($U\simeq 10^{-2}$) is particularly evident.

At this stage, Freesia is already mildly enriched, showing a metallicity of $Z\simeq 0.3\,{\rm Z}_{\odot}$ in the dense region, with central peaks up to $Z\simeq 2\,{\rm Z}_{\odot}$ in Freesia-A and Freesia-B. The surrounding gas is enriched at a mean $Z\simeq(0.02\pm 0.06)\,{\rm Z}_{\odot}$ up to $\sim 5\,{\rm kpc}$ from the stellar components. Freesia has a deeper potential well ($M_{\rm dm}\simeq 10^{11}{\rm M}_{\odot}$) with respect to the typical metal polluting galaxy (M_{\rm dm}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}\times 10^{9}{\rm M}_{\odot}, c.f.r Pallottini et al. 2014a): indeed a M_{\rm dm}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10^{8}{\rm M}_{\odot} halo is needed to have a galaxy with a mean Z\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}0.05\,{\rm Z}_{\odot} (Jeon et al., 2015, see in particular Fig. 4 therein); thus SN shocks originating from Freesia are less effective in enriching the intergalactic medium, as only the galaxy immediate surroundings can be easily accessed.

3.3 Thermo-chemical structure

To analyse the thermo-chemical structure of the gas we look at the density-temperature phase-diagrams, i.e. mass-weighted probability density functions (PDF) in the $n$-$T$ plane. In Fig. 6 we plot the phase-diagrams weighted by the total gas mass, the molecular gas mass, and the ionised component for the material within $\simeq 4.1\,{\rm kpc}$ from Freesia.

Considering the total mass (upper left panel of Fig. 6), we see that $\simeq 60\%$ of the gas is photoheated ($T\simeq 10^{4}\rm K$), due to photo-electric heating on dust grains illuminated by the radiation generated by local sources contributing to the Habing field. Lower temperatures can be reached by gas at density n\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10\,{\rm cm}^{-3}, accounting for $\simeq 35\%$ of the mass budget. The remaining \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}5\% is shock heated ($T\simeq 10^{6}\rm K$) by SN explosions and by accretion. The density has two small peaks around $n\simeq 1{\rm cm}^{-3}$ and $n\simeq 10^{2}{\rm cm}^{-3}$ that are superimposed to a roughly flat distribution. Overall, the total gas phase-diagram is similar to Pallottini et al. (2017b); this is expected since in the former work a uniform ISRF is assumed, and in Freesia we find that the Habing field is rather uniform few ${\rm kpc}$ away from the stellar component.

The double-peaked nature of the density distribution is due to the presence of the molecular and ionised components, as it is clear from the lower panels of Fig. 6. ${\rm{H_{2}}}$ is concentrated at high density (lower left panel of Fig. 6), with the peak at $n\simeq 5\times 10^{2}{\rm cm}^{-3}$. The diagram features a prominent peak at $T\simeq 3\times 10^{2}\rm K$ and a less pronounced one around $T\simeq 10^{3}\rm K$. The presence of the latter is linked to the formation of ${\rm{H_{2}}}$ in gas with $Z\sim 10^{-3}{\rm Z}_{\odot}$, i.e. the dust channel is disfavoured with respect to gas phase formation, which is enhanced at $T\simeq 10^{3}\rm K$. Note that there is a low (\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}1\%) amount of molecular hydrogen at even higher temperatures, albeit at lower (n\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10\,{\rm cm}^{-3}) densities. This gas is partially molecular and its presence is possible because of shielding from local LW sources. Recall that at this redshift ${\rm{H_{2}}}$ lines fall into the spectral window of SPICA (Spinoglio et al., 2017; Egami et al., 2018), thus in the future it will be possible to test the presence of such relatively high temperature ($T\sim 10^{3}$) ${\rm{H_{2}}}$, as the 0-0 S(1) at $17\,\mu\rm m$ and 0-0 S(0) at $28\,\mu\rm m$ are in the optimal sensitivity window of the instrument.

${\rm H}^{+}$ is responsible for the low density peak that is seen in the total gas diagram (lower right panel of Fig. 6). The bulk ($\simeq 70\%$) of ionised gas is centred at $n\simeq 5\times 10^{-2}{\rm cm}^{-3}$ and it is photoionised at a temperature $T\simeq 10^{4}\rm K$. Out of the total ${\rm H}^{+}$, $\simeq 25\%$ of the gas is in a shock heated state and mostly collisionally ionised. The remaining $\simeq 5\%$ of the ionised gas has typical densities n\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10\,{\rm cm}^{-3} and is only partially ionised.

The ionisation structure is better appreciated by looking at upper right panel of Fig. 6, that shows the phase-diagram of density vs ionised fraction, i.e. $f_{\rm ion}\equiv n_{\rm H+}/(n_{\rm H}+n_{\rm H+}+2\,n_{\rm H2})$. In Freesia, the ionised fraction decreases with density roughly as $f_{\rm ion}\propto n^{-0.5}$, with a dispersion of order of 0.3. The gas is found to be fully ionised ($f_{\rm ion}\simeq 1$) only in low density regions (n\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}1{\rm cm}^{-3}), while in potentially molecular regions ($n\sim 5\times 10^{2}{\rm cm}^{-3}$) $f_{\rm ion}\simeq 10^{-3}$. Such is likely a spurious result deriving from unresolved high density ${\rm H}^{+}$ regions (see Sec.s 2.2 and 2.4).

Radiation gives a non-negligible contribution to the ionised fraction both in high and low density regions. Assuming local thermodynamical equilibrium, the ionisation fraction due to collisions can be written as $f_{\rm ion}^{\rm coll}=\sqrt{(1+\xi)^{2}-1}-\xi$, where $\xi=\Gamma_{\rm H}/(2n\alpha_{\rm rec})$, with $\Gamma_{\rm H}=\Gamma_{\rm H}(T)$ and $\alpha_{\rm rec}=\alpha_{\rm rec}(T)$ being the collisional ionisation and recombination rate, respectively (see e.g. Dayal et al., 2008; Pallottini et al., 2014b). In Fig. 6 we plot $f_{\rm ion}^{\rm coll}$ for T\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10^{4}\rm K as a black solid line. $f_{\rm ion}$ is compatible with collisions only when n\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10^{-2}{\rm cm}^{-3}, while for progressively higher density the contributions from radiation becomes dominant, in particular considering that $f_{\rm ion}^{\rm coll}=0$ for T\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10^{3}\rm K, that is the typical temperature of the n\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10^{2}{\rm cm}^{-3} gas.

4 FIR emission properties

In this work we study the following FIR emission lines: [C $\scriptstyle\rm II$] $158\mu$ ${}^{2}P_{3/2}\rightarrow{}^{2}P_{1/2}$, [N $\scriptstyle\rm II$] $122\mu$ ${}^{3}P_{2}\rightarrow{}^{3}P_{1}$, and [O $\scriptstyle\rm III$] $88\mu$ ${}^{3}P_{1}\rightarrow{}^{3}P_{0}$, and we analyse the abundance and spatial distribution of their relative ions, ${\rm C}^{+}$ (ionisation potential $11.26\,\rm eV$), ${\rm O}^{++}$($35.11\,\rm eV$), and ${\rm N}^{+}$ ($14.53\,\rm eV$).

4.1 Imaging of the FIR emission and ions

In Fig. 7 we show the emission maps of Freesia in [C $\scriptstyle\rm II$], [N $\scriptstyle\rm II$], and [O $\scriptstyle\rm III$]. As a reference, the properties are summarised in Tab. 1.

The brightest line is [C $\scriptstyle\rm II$], with $L_{\rm CII}\simeq 7.7\times 10^{7}{\rm L}_{\odot}$. The two main peaks ($\Sigma_{\rm[CII]}\simeq 5\times 10^{7}{\rm L}_{\odot}\,{\rm kpc}^{-2}$) are spread over the Freesia-A and Freesia-B, with extents of the order of $1.5\,{\rm kpc}$. The only other prominent structure is the high density filament ($n\simeq 10^{2}{\rm cm}^{-3}$) North-West of Freesia, featuring a lower brightness – $\Sigma_{\rm[CII]}\simeq 5\times 10^{6}{\rm L}_{\odot}\,{\rm kpc}^{-2}$ – due to the lower metal content ($Z\simeq 5\times 10^{-2}{\rm Z}_{\odot}$, cfr. with Fig. 5). The bright spots are embedded in a faint – $\Sigma_{\rm[CII]}\simeq 5\times 10^{4}{\rm L}_{\odot}\,{\rm kpc}^{-2}$ – halo that marks the extent of region that has been metal enriched by Freesia. Note that such diffuse halo in Freesia gives only a small contribution to the [C $\scriptstyle\rm II$]: the $1.5\,{\rm kpc}$ extension of the emitting region extension is mainly determined by the presence of the two stellar components.

The emission from ions with higher ionisation state has a different morphology than [C $\scriptstyle\rm II$]. Both [N $\scriptstyle\rm II$] and [O $\scriptstyle\rm III$] show are less extended than [C $\scriptstyle\rm II$]. For [N $\scriptstyle\rm II$] and [O $\scriptstyle\rm III$], the luminosity of Freesia-A is a factor \lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10 larger than Freesia-B, while for [C $\scriptstyle\rm II$]  the factor is $\simeq 4$. Since they trace similar material (see later Fig. 8), the two lines have cospatial emission peaks, located in high density (n\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10^{2}{\rm cm}^{-3}) ionised regions that blister the disk of Freesia and enclose a total size \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}0.5{\rm kpc}. In particular, the [N $\scriptstyle\rm II$] line has a low luminosity, $L_{\rm NII}\simeq 5.3\times 10^{5}{\rm L}_{\odot}$. This is due to the smaller cooling efficiency with respect to the other lines (Dalgarno & McCray, 1972): the maximum surface brightness is $\Sigma_{\rm[NII]}\simeq 2\times 10^{7}{\rm L}_{\odot}\,{\rm kpc}^{-2}$, roughly one order of magnitude smaller than [C $\scriptstyle\rm II$] and [O $\scriptstyle\rm III$]. While the [O $\scriptstyle\rm III$] surface brightness is higher than [C $\scriptstyle\rm II$] ($\Sigma_{\rm[OIII]}\simeq 5\times 10^{8}{\rm L}_{\odot}\,{\rm kpc}^{-2}$), the smaller emitting region makes its total luminosity lower, i.e. $L_{\rm OIII}\simeq 2.1\times 10^{7}{\rm L}_{\odot}$. The [O $\scriptstyle\rm III$] shows a more diffuse halo with surface brightness $\Sigma_{\rm[OIII]}\simeq 5\times 10^{4}{\rm L}_{\odot}\,{\rm kpc}^{-2}$. However with respect to the [C $\scriptstyle\rm II$] halo the morphology is very different, since the [O $\scriptstyle\rm III$] halo is confined in the east direction, in correspondence of the low ionisation field ($U\sim 10^{-4}$).

Along with emission lines, the corresponding ion surface densities ($\Sigma_{\rm C+}$, $\Sigma_{\rm N+}$, and $\Sigma_{\rm O++}$) are plotted in Fig. 7. [C $\scriptstyle\rm II$] structure follows the ${\rm C}^{+}$ ion morphology. ${\rm C}^{+}$ is present in both diffuse and molecular material, and – without ionising radiation – we expect $\Sigma_{\rm C+}\propto nZ$ (Ferrara et al., in prep., 2019). Since the Habing field is almost constant on these scales, a rough proportionality between the luminosity and the ion abundance is expected. In both Freesia-A and Freesia-B the gas features a flat $\Sigma_{\rm C+}\simeq 10^{5}{\rm M}_{\odot}\,{\rm kpc}^{-2}$, that rapidly decreases to \Sigma_{\rm C+}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10^{2}{\rm M}_{\odot}\,{\rm kpc}^{-2} in the diffuse halo. As the unstable filament north-west of the galaxy has a lower metal enrichment than the two star forming components, its ${\rm C}^{+}$ abundance is similar to that of the halo. ${\rm N}^{+}$ and ${\rm O}^{++}$are similarly distributed and trace the ionised high density regions, as for the corresponding lines. In both cases, halos of low surface density material are present (\Sigma\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10{\rm M}_{\odot}\,{\rm kpc}^{-2}), however the ISRF is not high enough in order for the corresponding lines to be emitted efficiently. Note that the order of magnitude difference in the surface density is mostly due to the difference between the mass abundance of the elements.

4.2 FIR lines as a tracer of the ISM state

Using the phase-diagrams we can analyze the emission structure of the ISM in Freesia. In Fig. 8 we plot the phase-diagrams weighted by the luminosity of [C $\scriptstyle\rm II$] ($n$-$G$ plane), [N $\scriptstyle\rm II$] ($n$-$U$ plane), and [O $\scriptstyle\rm III$] ($n$-$U$ plane).

The [C $\scriptstyle\rm II$] diagram shows that most of the contribution to its emission comes from gas with $n\simeq 160\,{\rm cm}^{-3}$ illuminated by $G\simeq 20\,G_{0}$, i.e. dense gas in the two star forming components embedded in the average radiation field; this peak spans roughly an order of magnitude in both axis. Contribution from lower density gas (n\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10{\rm cm}^{-3}) is present but subdominant (\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10\%). Note that for a weak field (G\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}G_{0}) at low densities (n\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10{\rm cm}^{-3}), the contribution to the emission is also suppressed by the CMB (da Cunha et al., 2013; Vallini et al., 2015; Pallottini et al., 2015). Overall these results are consistent with our previous findings (Pallottini et al., 2017a, b), i.e. most of the [C $\scriptstyle\rm II$] emission is associated with material close to the molecular regions. A peak at high density is expected from our benchmark (Sec. 2.4, in particular see Fig. 3), which also shows that [C $\scriptstyle\rm II$] emission is favoured at relatively higher values of the Habing field (G\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}5\times 10^{2}G_{0}). These regions are relatively rare in Freesia, as they are associated with star forming regions, thus in our galaxy the peak contribution comes from regions with a milder radiation field.

Similar distributions are found for [N $\scriptstyle\rm II$] and [O $\scriptstyle\rm III$]. Two types of emitting regions are present: i) a stripe with roughly $U\propto n^{-2}$, which accounts for most of the luminosity of both lines, and ii) a \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}5\% contribution from diffuse ionised gas with n\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}1\,{\rm cm}^{-3} and $U\simeq 5\times 10^{-4}$.

In both cases the emission peaks in dense ${\rm H}^{+}$ regions. However, some differences are present: the [N $\scriptstyle\rm II$] peak is concentrated at $U\simeq(2\pm 0.1)\times 10^{-3}$ for gas with densities $n\simeq(95\pm 40)\,{\rm cm}^{-3}$, while [O $\scriptstyle\rm III$] shows a larger range for the ionisation parameter, $U\simeq(6\pm 2)\times 10^{-3}$, and arises from gas with lower densities, $n\simeq(53\pm 40)\,{\rm cm}^{-3}$. Thus, additionally to the higher typical cooling efficiency, the higher [O $\scriptstyle\rm III$] luminosity with respect to [N $\scriptstyle\rm II$] is also partially due to contribution from gas with $n\simeq 10\,cc$, which is more abundant then the $n\simeq 95\,{\rm cm}^{-3}$ gas that dominates the [N $\scriptstyle\rm II$] line emission.

To analyse the spectral hardness of the radiation field, in Fig. 8 we also plot phase-diagram in the $G$-$U$ plane weighted by photon number ($N_{\gamma}$). Overall, the bulk of the ISRF is relatively soft, with $\langle U\rangle=7\times 10^{-4}\pm 0.03$ and $\langle G\rangle=(15.0\pm 71)G_{0}$, and has a trend roughly given by $U\propto(G/G_{0})^{1/2}$ for G\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10^{2}G_{0}

The region in the phase-diagram corresponding to high $U$ and $G$ is characterised by young stellar clusters that have removed the gas from their surrounding, through their radiative and mechanical feedback (see also Fig. 5). This is motivated as follows: before absorption, a typical stellar SED would yield $G/G_{0}\sim U(n/{\rm cm}^{-3})10^{3}$ (see Fig. 1); as $G\sim 10^{3}G_{0}$ at most, the density in $U\sim 1$ regions must be n\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}1\,{\rm cm}^{-3}, i.e. lower than the original $n\simeq 3\times 10^{2}{\rm cm}^{-3}$ allowing the formation of stars.

Summarising, most of the radiation surrounding Freesia is non ionising, i.e. the escape fraction calculated on a sphere of radius of $4.1{\rm kpc}$ is of order $\simeq 2\%$, in agreement with estimate from observations of lower redshift galaxies with similar brightness (e.g. Bouwens et al., 2016; Grazian et al., 2017) and averaged values of simulated galaxy with similar masses (Xu et al., 2016). However, large variations are expected with different evolutionary stage and radii considered (Trebitsch et al., 2017); a detailed analysis is left for future works.

4.3 Spatial offsets

We have seen that while low (${\rm C}^{+}$) and high (${\rm O}^{++}$) ionisation species have a roughly similar spatial distribution, their corresponding emission structure is very different as a result of the modulation imposed by the IRSF, which is roughly spatially uniform in the Habing band, and very patchy in ionising radiation. This causes a spatial offsets between [C $\scriptstyle\rm II$] and [O $\scriptstyle\rm III$] emitting regions. In Freesia there are two different kinds of spatial offsets. In Freesia-A [C $\scriptstyle\rm II$] is extended and peaks at the edge of the disk, while [O $\scriptstyle\rm III$] is concentrated at the center. This configuration results in a $\simeq 250\rm pc$ offset for the [C $\scriptstyle\rm II$] and [O $\scriptstyle\rm III$] lines arising from Freesia-A. In Freesia-B  [C $\scriptstyle\rm II$] has luminosity $\sim 10^{7}{\rm L}_{\odot}$, but L_{\rm OIII}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10^{5}{\rm L}_{\odot}. While Freesia-B would not be detected in [O $\scriptstyle\rm III$] even with an extremely deep ALMA observation, its [C $\scriptstyle\rm II$] luminosity would move the center of the emission away from Freesia-A; thus a marginally resolved observation of the system would reveal an offset of $\simeq 2\rm kpc$ between [C $\scriptstyle\rm II$] and [O $\scriptstyle\rm III$]. This is similar to what is observed in BDF-3299 (Carniani et al., 2017, in particular see Fig. 4.)

We recall that in BDF-3299 $L_{\rm CII}/L_{\rm OIII}=0.27$ (see Tab. 4 of Carniani et al., 2017), while in Freesia we find $L_{\rm CII}/L_{\rm OIII}=3.73$. In Freesia, the [O $\scriptstyle\rm III$] emission predominantly comes from $Z\simeq 0.5\,{\rm Z}_{\odot}$ gas in Freesia-A; there, $\Sigma_{\rm[CII]}/\Sigma_{\rm[OIII]}\simeq 0.1$. The [O $\scriptstyle\rm III$] emission is limited to regions with an hard radiation field (U\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10^{-2}), while [C $\scriptstyle\rm II$] can excited by the mean UV ISRF ($G\simeq 7.9\,G_{0}$) in the diffuse medium surrounding Freesia-A and thus is emitted efficiently also from the material surrounding Freesia-A, yielding a ratio of total luminosity $L_{\rm CII}/L_{\rm OIII}<1$. This difference between Freesia and BDF-3299 can be explained if the latter has a larger ionised region; however it might be due to a different configuration of the system: further investigation is needed to have achieve a full classification.

Interestingly, neither in Freesia nor in the other simulated galaxies we find situations in which $L_{\rm CII}/L_{\rm OIII}\ll 1$, as shown in observations where [O $\scriptstyle\rm III$] is present but [C $\scriptstyle\rm II$] is undetected (Inoue et al., 2016; Laporte et al., 2017). We note that the three simulated galaxies at $9\leq z<12$ presented in Katz et al. (2019, see Fig. 11 therein) are also [C $\scriptstyle\rm II$]-dominated with typical values $L_{\rm CII}/L_{\rm OIII}\sim 10-50$. At present, simulations are apparently unable to reproduce the observed low $L_{\rm CII}/L_{\rm OIII}$ ratios. This issue requires further study and we leave it for future, more specific analysis.

4.4 Spectral shifts

Spectral shifts between different FIR lines are present in Freesia. To quantify this effect, we build spectra ($F_{\nu}$) as a function of the line of sight (l.o.s.) velocity ($v_{los}$). Each gas patch along the line of sight gives a contribution corresponding to its luminosity, which is kernel weighted by a Gaussian centred at the peculiar velocity of the gas and with a width that accounts for the thermal and turbulent broadening. Full detail of the model are given in Kohandel et al., in prep. (2019), along with an in-depth analysis of the kinematics of the [C $\scriptstyle\rm II$].

In Fig. 9 we plot the [C $\scriptstyle\rm II$], [N $\scriptstyle\rm II$], and [O $\scriptstyle\rm III$] spectra (see Tab. 1 for reference). The [C $\scriptstyle\rm II$] appears to be more peaked in velocity space, with a line width (second moment of the spectrum) $\sigma_{v}=93.0\,\,{\rm{km\,s^{-1}}}$. The [C $\scriptstyle\rm II$]  emission comes from both stellar components, with Freesia-A providing a broader turbulent disk component with |v_{\rm los}|\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}150\,\,{\rm{km\,s^{-1}}}, while Freesia-B provides a concentrated contribution at $v_{\rm los}=0$, that gives rise to the peak in the total spectrum. This is a rather common situation when multiple, possibly merging, components are present in the same system (Kohandel et al., in prep., 2019). Note that only \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10\% of the emission comes from high velocity (|v_{\rm los}|\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}150\,\,{\rm{km\,s^{-1}}}) gas, as most of the emission is concentrated in the dense ISM of the galaxy; only a small fraction possibly due to outflowing material, because of its lower metallicity (Gallerani et al., 2018). The [O $\scriptstyle\rm III$] line is shifted with respect to [C $\scriptstyle\rm II$] by $v_{\rm los}=138.14\,\,{\rm{km\,s^{-1}}}$, roughly one third of what is observed for BDF3299 (Carniani et al., 2017) as its emission is dominated by Freesia-A. In Freesia [O $\scriptstyle\rm III$] has a $\sigma_{v}$ that is $1.75$ larger than [C $\scriptstyle\rm II$], and it additionally features more prominent high velocity wings at |v_{\rm los}|\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}150\,\,{\rm{km\,s^{-1}}}, that are due to the diffuse low surface brightness [O $\scriptstyle\rm III$] halo in the West region (see Fig. 7). This possibly makes [O $\scriptstyle\rm III$]  a better tracer of outflowing gas. The spectrum of [N $\scriptstyle\rm II$] feature double peaks, due to both the contribution from Freesia-A and Freesia-B, that are almost coincident with the location of [C $\scriptstyle\rm II$] and [O $\scriptstyle\rm III$] peaks in the velocity space. A luminosity ratio [C $\scriptstyle\rm II$]/[N $\scriptstyle\rm II$]$\sim 1/100$ is expected also from theoretical works (Vallini et al., 2013). This makes [N $\scriptstyle\rm II$] observations very challenging.

Note that the spatial offsets and spectral shifts typically observed in high-$z$ galaxies are between UV/Ly-$\alpha$ and [C $\scriptstyle\rm II$]; only in a handful of cases [O $\scriptstyle\rm III$] is also available (Carniani et al., 2017, see Fig. 6.). Such analysis requires additional post-processing work to compute the UV continuum and Ly-$\alpha$ radiative transfer, as in e.g. Behrens et al. (2019). This is also deferred to a forthcoming paper (Kohandel et al., in prep., 2019).

5 On the [CII]-SFR relation

In Fig. 10 we plot the integrated SFR-[C $\scriptstyle\rm II$] relation for Freesia and the other galaxies that are on average within $(56.9\pm 21.6)\,{\rm kpc}$. For our sample of 11 galaxies we find that [C $\scriptstyle\rm II$] is increasing with SFR, however the slope is shallower with respect to the local De Looze et al. (2014) relation: galaxies with lower star formation ({\rm SFR}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}1{\rm M}_{\odot}\,{\rm yr}^{-1}) lie above or on top the local relation, while as we go to progressively high rates ({\rm SFR}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}5{\rm M}_{\odot}\,{\rm yr}^{-1}) galaxies fall below the relation. The trend for the simulated galaxies is not clear, as the dispersion is large and the sample is limited.

As noted in Carniani et al. (2018a), this trend is similar to what is observed at $\mbox{high-}z$ : it is unclear whether the local De Looze et al. (2014) relation holds at high redshift because of the low statistical significance of the observed sample, and there is evidence that the dispersion is larger by a factor $\times 1.8$ with respect to the local [C $\scriptstyle\rm II$]-SFR. In particular, Freesia is within $2\sigma$ from the De Looze et al. (2014) relation, similarly to what is observed in high-$z$ galaxies (Carniani et al., 2018a). Low metallicity alone cannot fully explain the tension from the local relation. The mean metallicity of the gas in Freesia and most of the other simulated galaxies is $Z=0.5\,{\rm Z}_{\odot}$. However, for $Z=0.5\,{\rm Z}_{\odot}$ the Vallini et al. (2015) model is consistent with the De Looze et al. (2014) relation.

Note that at redshift 12\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}z\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}9 the three simulated galaxies of Katz et al. (2019) have properties that are consistent with the local [C $\scriptstyle\rm II$]-SFR (De Looze et al., 2014). However, Katz et al. (2019) emission lines are computed starting from a library of cloudy models with constant gas temperature, i.e. not computing the temperature structure of the PDR. The temperature is taken from the hydrodynamical simulation and is computed with a simpler network with respect to cloudy or other simulations, e.g. not considering non-equilibrium metal cooling (Capelo et al., 2018). The resulting ion abundances and emission lines may be considerably different.

Instead, in the following Sec.s 5.1 and 5.2, we show that falling below the local [C $\scriptstyle\rm II$]-SFR relation is connected to the galaxy “star-burstiness”, i.e. the evolutionary stage in which the galaxy lies above the Kennicutt-Schmidt relation, therefore experiencing an enhanced stellar feedback. The theoretical background for this argument is worked out in detail in a companion work (Ferrara et al., in prep., 2019).

5.1 Connection with Kennicutt-Schmidt relation

The Kennicutt-Schmidt (KS) relation in our simulated sample is plotted in Fig. 11 (right panel) in the Krumholz et al. (2012) “formulation”, i.e. $\dot{\Sigma}_{\star}$ vs $\Sigma/t_{\rm ff}$. In this formulation the scatter is lower, and the relation holds for single MCs, local and moderate redshift unresolved galaxies.

Freesia and galaxy 1 are under the De Looze et al. (2014) relation and above the KS151515Recall that Freesia has recently experienced a burst of star formation (Fig. 4), this is the reason why it sits above the relation., while the opposite holds for galaxies 4 and 9, which are over the local [C $\scriptstyle\rm II$]-SFR and below the KS, while galaxies 6 and 8 are consistent with both. This is an indication that the two relations are intimately linked. However galaxies 2, 5, 7, and 10 are below both relations, while 3 is on the De Looze et al. (2014) but above the KS.

5.2 Spatially resolved relations

To clarify the connection between KS relation and [C $\scriptstyle\rm II$]-SFR relation, we analyse their spatially resolved versions, that are generally considered more fundamental in physical terms than the integrated ones (e.g. for [C $\scriptstyle\rm II$] see Herrera-Camus et al., 2015). Freesia is the only galaxy in our simulated sample with a SFR comparable to the one for currently observed galaxies, thus we focus again on that system.

As we are interested in the deviations from the local relations, it is convenient to parameterize them as follows. For the [C $\scriptstyle\rm II$]-SFR we use

[TABLE]

with $\Sigma_{\rm[CII]}$ and $\dot{\Sigma}_{\star}$ in units of ${{\rm L}_{\odot}\,{\rm kpc}^{-2}}$ and ${\rm M}_{\odot}\,{\rm yr}^{-1}\,{\rm kpc}^{-2}$, respectively (De Looze et al., 2014, data from Tab. 2). For the KS we can define

[TABLE]

where $\Sigma_{g}$ is in units of ${\rm M}_{\odot}\,{\rm kpc}^{-2}$ (Heiderman et al., 2010, see their eq. 2).

In Fig. 12 we plot the deviation from the relations along with the three quantities determining the local relations ($\Sigma_{\rm[CII]}$, $\dot{\Sigma}_{\star}$, $\Sigma_{g}$). To have a fair comparison with observations, the maps are smoothed at 100 pc, similar to the De Looze et al. (2014) resolution.

It is striking that in almost all the spots where Freesia is deficient in [C $\scriptstyle\rm II$] ($\log({\rm[CII]}/{\rm[CII]}_{\rm DL})<0$), it is locally star-bursting ($\log\kappa_{\rm s}>0$); this is particularly evident at the centres of both Freesia-A and Freesia-B, and at the location of the young stars stripped out in the halo of the system. Vice versa, being above the De Looze et al. (2014) relation implies a lack of star formation, as seen in the surroundings of Freesia-B and in the dense gas filament North-West of Freesia-A.

However, a location above (below) the KS relation is a necessary but not sufficient condition for a gas path to be under (over) luminous in [C $\scriptstyle\rm II$], as seen in the material around the filament and in the outer edges of the system. This happens because metallicity plays a secondary role in the link between the two relations (Fig. 12, bottom right). On the one hand, where the metallicity is very low (Z\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10^{-2}{\rm Z}_{\odot}), the [C $\scriptstyle\rm II$] is fainter than expected even for high $\Sigma_{g}$ (\Sigma_{g}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10^{8}{\rm M}_{\odot}\,{\rm kpc}^{-2}), as carbon abundance limits the emission. On the other hand, a low metallicity implies less dust, the main catalyst of ${\rm{H_{2}}}$ formation, thus consequently stars. This entails a lower $\dot{\Sigma}_{\star}$ for the patches of gas with similar $\Sigma_{g}$. Such low mean metallicity is not found for relatively massive (M_{\star}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10^{8}{\rm M}_{\odot}) and star forming ({\rm SFR}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}0.1{\rm M}_{\odot}\,{\rm yr}^{-1}) galaxies; low mean ISM metallicities (Z\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}0.05\,{\rm Z}_{\odot}) are typical of smaller (M_{\star}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10^{7}{\rm M}_{\odot}) systems with lower star formation ({\rm SFR}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}0.1\,{\rm M}_{\odot}\,{\rm yr}^{-1} Jeon et al. 2015; Jeon & Bromm 2019); for these galaxies, metallicity can play a role as relevant as burstiness, however currently no ALMA [C $\scriptstyle\rm II$] detection is available for galaxies with {\rm SFR}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}1\,{\rm M}_{\odot}\,{\rm yr}^{-1}.

The two resolved relations can be further analysed by extracting the PDFs of $\Sigma_{\rm[CII]}$-$\dot{\Sigma}_{\star}$ and $\dot{\Sigma}_{\star}$-$\Sigma_{g}$. From Fig. 13 we can see that $\Sigma_{\rm[CII]}$ is increasing with $\dot{\Sigma}_{\star}$ in Freesia: the slope is almost flat (\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}0.25) for low (\dot{\Sigma}_{\star}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}0.1{\rm M}_{\odot}\,{\rm yr}^{-1}\,{\rm kpc}^{-2}) and high (\dot{\Sigma}_{\star}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}5{\rm M}_{\odot}\,{\rm yr}^{-1}\,{\rm kpc}^{-2}) formation rates, while the trend is nearly linear for intermediate $\dot{\Sigma}_{\star}$. The scatter in the relation is decreasing with increasing $\dot{\Sigma}_{\star}$, because of the smaller spread in metallicity of the gas as we go to progressively higher $\Sigma_{\rm[CII]}$.

The local $\Sigma_{\rm[CII]}$-$\dot{\Sigma}_{\star}$ is fitted with an almost linear slope (eq. 12b) and the observed data range is 10^{5}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}\Sigma_{\rm[CII]}/{\rm L}_{\odot}\,{\rm kpc}^{-2}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10^{6.5} and 10^{-3}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}\dot{\Sigma}_{\star}/{\rm M}_{\odot}\,{\rm yr}^{-1}\,{\rm kpc}^{-2}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}1, with only few observed points (from the galaxy NGC1569) having 10^{-1}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}\dot{\Sigma}_{\star}/{\rm M}_{\odot}\,{\rm yr}^{-1}\,{\rm kpc}^{-2}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}1, where the data seem to indicate the presence of a deficit with respect to the linear trend (De Looze et al., 2014, see in Fig. 2 therein). The local relation is below the one found for Freesia at \dot{\Sigma}_{\star}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}0.1{\rm M}_{\odot}\,{\rm yr}^{-1}\,{\rm kpc}^{-2}, while it is above that at \dot{\Sigma}_{\star}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}1{\rm M}_{\odot}\,{\rm yr}^{-1}\,{\rm kpc}^{-2}. The overlap is not perfect in the intermediate region, because of the patch of low surface brightness regions (\Sigma_{\rm[CII]}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10^{5.5}{\rm L}_{\odot}\,{\rm kpc}^{-2}) at 0.01\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}\dot{\Sigma}_{\star}/{\rm M}_{\odot}\,{\rm yr}^{-1}\,{\rm kpc}^{-2}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}1 that drags down the mean relation from Freesia.

In Fig. 13 we see that most of the $\Sigma_{\rm[CII]}$-$\dot{\Sigma}_{\star}$ high-z data – obtained by integrating the various galaxy components and estimating their UV and IR size (see Carniani et al., 2018a, for details) – are nicely consistent with the average value extracted from Freesia, while they fall below the local relation. This is an indication those galaxy are dominated by \dot{\Sigma}_{\star}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}{\rm M}_{\odot}\,{\rm yr}^{-1}\,{\rm kpc}^{-2}, critical point where the local $\Sigma_{\rm[CII]}$-$\dot{\Sigma}_{\star}$ deviate from linear by saturating to an almost constant value (Ferrara et al., in prep., 2019).

To close the loop, we plot the PDF of the KS in Fig. 14. The local relation and the average found in Freesia are in good agreement up to \dot{\Sigma}_{\star}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10{\rm M}_{\odot}\,{\rm yr}^{-1}\,{\rm kpc}^{-2} (\Sigma_{g}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10^{9}{\rm M}_{\odot}\,{\rm kpc}^{-2}); for higher values, the intense and concentrated star formation and the consequent strong feedback cause these parts of the galaxy to deviate from the averaged KS. The majority of the regions of the galaxy is below the local $\Sigma_{\rm[CII]}$-$\dot{\Sigma}_{\star}$, in particular the patches of gas at high $\dot{\Sigma}_{\star}$, which are major contributors of its luminosity and SFR; thus, when spatially integrating $\Sigma_{\rm[CII]}$ and $\dot{\Sigma}_{\star}$, Freesia results to be below the local [C $\scriptstyle\rm II$]-SFR relation and consistent only within $2\sigma$ (Fig. 10).

6 Conclusions

We have studied the formation and evolution of a sample of Lyman Break galaxies in the Epoch of Reionisation by using crafted, cosmological zoom-in (spatial resolution $\sim 10$ pc) simulations, as part of the serra suite.

The serra simulations are based on a customised version of the Adaptive Mesh Refinement code ramses (Teyssier, 2002). The ISM thermo-chemical evolution is followed via a non-equilibrium network generated with krome (Grassi et al., 2014), that allows a precise tracking of the formation of ${\rm{H_{2}}}$, that can be converted into stars. With respect to previous works (Pallottini et al., 2017a, b), the present simulations perform a full on-the-fly radiative transfer of the interstellar radiation field (ISRF) thanks to ramses-rt (Rosdahl et al., 2013). In the post-processing phase we compute the intensities of several FIR lines ([C $\scriptstyle\rm II$], [N $\scriptstyle\rm II$], and [O $\scriptstyle\rm III$]) from a grid of models obtained with the photo-ionisation code cloudy V17 (Ferland et al., 2017). These calculations also account for the turbulent structure of the gas (Vallini et al., 2015). This procedure allows a fair comparison of the results with ALMA observations of high-$z$ galaxies.

At $z=8$, the most massive galaxy in the simulation is “Freesia”. It has a stellar age of $t_{\star}\simeq 409\,\rm Myr$, a stellar mass of $M_{\star}\simeq 4.2\times 10^{9}{\rm M}_{\odot}$, and a star formation rate of ${\rm SFR}\simeq(11.5\pm 1.8)\,{\rm M}_{\odot}\,{\rm yr}^{-1}$, due to a recent burst. The galaxy is composed by two concentrated ($\sim 200\,\rm pc$) stellar components (“A” and “B”) separated by $\simeq 2.5$ kpc. Freesia-A dominates both the mass ($\simeq 85\%$) and star formation ($90\%$) of the system. Around Freesia, other 11 galaxies are found within $56.9\pm 21.6$ kpc ($>2$ virial radii); while in our previous work such galaxies were present (Pallottini et al., 2017b), their SFR was likely suppressed in their formation stages, due to lack of proper treatment of the ISRF, which was spatially uniform and dominated in intensity by the most massive galaxy.

The properties of Freesia are overall similar to the ones found in previous work at the same evolutionary stage (Pallottini et al., 2017b), featuring comparable SFR and $M_{\star}$. The galaxy develops broken spiral disk of gas with radius $\simeq 0.5\,{\rm kpc}$ and average metallicity $Z\simeq 0.5\,{\rm Z}_{\odot}$. The molecular gas is primarily found at gas densities $n\simeq 200\,{\rm cm}^{-3}$, i.e. about 2$\times$ denser than previously found (Pallottini et al., 2017b) as a consequence of the higher local ISRF. This is because higher densities are required in order for the gas to self-shield from the stronger Lyman-Werner radiation. The bulk ($70\%$) of the ionised ISM has low densities ($n=5\times 10^{-2}{\rm cm}^{-3}$); $25\%$ is collisionally ionised by shocks; only $5\%$ is found in the high density regions (n\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10\,{\rm cm}^{-3}), where ionisation is due to the internal presence of young clusters, which gives large fluctuations in the ionising ISRF.

The ISRF of the galaxies features a mean (variance) Habing flux $G=7.9G_{0}$ ($G=23\,G_{0}$) and has a rather smooth spatial distribution. The ionisation radiation, parameterized by $U$, shows larger intensity variations in the range $U=2^{+20}_{-2}\times 10^{-3}$ and a patchy distribution peaked at the location of recent star formation events. The patchy ionisation structure yields an escape fraction of $f_{\rm esc}\simeq 2\%$, a value consistent with observations of lower redshift galaxies with similar brightness (e.g. Bouwens et al., 2016; Grazian et al., 2017), and averaged values of simulated galaxies with similar mass (Xu et al., 2016). We note, however, that large fluctuations are expected between different evolutionary stages (Trebitsch et al., 2017)

The [C $\scriptstyle\rm II$] line luminosity of Freesia is $L_{\rm CII}=7.73\times 10^{7}{\rm L}_{\odot}$. The emission extends on a few kpc scales, with peaks around Freesia-A, Freesia-B and dense non-starforming clumps found near the galaxy. The emission mostly comes from gas with $n\simeq 160\,{\rm cm}^{-3}$ and illuminated by an ISRF intensity $G\simeq 20\,G_{0}$. Instead, [O $\scriptstyle\rm III$] emission is concentrated around Freesia-A, the star-bursting component of the system, and its emission is dominated by gas with a wider range of properties, i.e. $n\simeq(53\pm 40)\,{\rm cm}^{-3}$ and $U\simeq(6\pm 2)\times 10^{-3}$. At the Freesia-A location the oxygen line is very bright, $\Sigma_{\rm[OIII]}/\Sigma_{\rm[CII]}\simeq 10$. However, the smaller extent of the [O $\scriptstyle\rm III$] emitting region implies a lower galaxy-integrated emission, i.e. $L_{\rm OIII}=2.07\times 10^{7}{\rm L}_{\odot}$. With respect to [C $\scriptstyle\rm II$], the [O $\scriptstyle\rm III$] show both a spatial offset ($\sim 2.5\,{\rm kpc}$) and a spectral shift of ($\simeq 120\,{\rm{km\,s^{-1}}}$), reminiscent of similar evidence found in some systems at high-$z$ (Carniani et al., 2017).

Freesia lies below the local [C $\scriptstyle\rm II$]-SFR relation (De Looze et al., 2014, within $\simeq 2\sigma$) as it is in a starburst phase, i.e. it sits above the Kennicutt-Schmidt (KS) relation (Schmidt, 1959; Kennicutt, 1998, within $\simeq 2\sigma$). Spatial analysis reveals that patches of the galaxy that are above the resolved local $\Sigma_{\rm[CII]}$-$\dot{\Sigma}_{\star}$ De Looze et al. (2014) relation are located below KS, and vice-versa.

In particular, due to the recent starburst, the dense center (\Sigma_{g}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10^{9}{\rm M}_{\odot}\,{\rm kpc}^{-2}) of Freesia lies above the KS (\dot{\Sigma}_{\star}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10{\rm M}_{\odot}\,{\rm yr}^{-1}\,{\rm kpc}^{-2}) relation. At such high star formation surface densities the $\Sigma_{\rm[CII]}$ increases with $\dot{\Sigma}_{\star}$ less rapidly than expected from the local relation De Looze et al. (2014), which however only covers the low-end of the $\dot{\Sigma}_{\star}$ range. Thus, the observed [C $\scriptstyle\rm II$]-SFR deficit can be primarily ascribed to negative stellar feedback during starburst phases, disrupting molecular clouds around star formation sites. This interpretation was originally proposed in Vallini et al. (2015); our results fully endorse it.

Metallicity effects have a weaker impact on the origin of [C $\scriptstyle\rm II$]-SFR deviations. Gas with extremely low (Z\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10^{-2}{\rm Z}_{\odot}) metallicities fails below both the local $\Sigma_{\rm[CII]}$-$\dot{\Sigma}_{\star}$ and the KS relations, because of lack of ${\rm C}^{+}$ ions and inefficient ${\rm{H_{2}}}$ formation, respectively. However, gas with Z\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10^{-2}{\rm Z}_{\odot} is scarce in Freesia and its combined contribution is subdominant when galaxy integrated relations are considered. This situation is common to most of the galaxies in our sample. A metallicity effect might be relevant for galaxies that cannot retain their metal production (Xu et al., 2016; Jeon & Bromm, 2019), that are expected to be smaller (M_{\star}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}5\times 10^{6}{\rm M}_{\odot}) and with a lower star formation ({\rm SFR}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}1{\rm M}_{\odot}\,{\rm yr}^{-1}), for which however there is currently no detection with ALMA.

Acknowledgments

This research was supported by the Munich Institute for Astro- and Particle Physics (MIAPP) of the DFG cluster of excellence “Origin and Structure of the Universe”. We thank A. Lupi, M. Trebitsch, J. Rosdahl, and the participants of the “The Interstellar Medium of High Redshift Galaxies” MIAPP conference for fruitful discussion. AF and SC acknowledge support from the ERC Advanced Grant INTERSTELLAR H2020/740120. LV acknowledges funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie Grant agreement No. 746119. We acknowledge use of the Python programming language (Van Rossum & de Boer, 1991), Astropy (Astropy Collaboration et al., 2013), Cython (Behnel et al., 2011), Matplotlib (Hunter, 2007), NumPy (van der Walt et al., 2011), pymses (Labadens et al., 2012), Pandas (McKinney, 2010), and SciPy (Jones et al., 2001).