Context

Hydrodynamical cosmological simulations are increasing their level of realism by considering more physical processes, having more resolution or larger statistics. However, one usually has to either sacrifice the statistical power of such simulations or the resolution reach within galaxies. Here, we introduce the NewHorizon project where a zoom-in region of  ~ (16  Mpc)3 embedded in a larger box is simulated at high resolution. A resolution of up to 40  pc, typical of individual zoom-in state-of-the-art resimulated halos is reached within galaxies, allowing such a simulation to capture the multi-phase nature of the interstellar medium.

Table of Content


1.1:  Initial conditions and resolution
1.2:  Radiative cooling and heating
1.3:  Star formation
1.4:  Feedback from massive stars
1.5:  MBHs and AGN
1.5.1:  MBH formation, growth and dynamics
1.5.2:  MBH spin evolution
1.5.3:  Radio and quasar modes of AGN feedback
1.6:  Identification of halos and galaxies

Basic properties

  1.1  Initial conditions and resolution

The NewHorizon  simulation is a zoom-in simulation from the 100 h − 1  Mpc size   Horizon − AGN simulation [22]. The Horizon  − AGN simulation initial conditions had a 10243 DM particles and initial grid resolution, with a ΛCDM cosmology. The total matter density is Ω m = 0.272, dark energy density ΩΛ  = 0.728, amplitude of the matter power spectrum σ8  = 0.81, baryon density Ω b = 0.045, Hubble constant H0  = 70.4  km s  − 1 Mpc − 1, and ns  = 0.967 compatible with the WMAP-7 data [40]. Within this large-scale box, we define a initial spherical patch of 10  Mpc radius with a 40963 effective resolution, i.e. a DM mass resolution of MDM,  hr = 1.2 × 106  M⊙. The high-resolution initial patch is embedded in buffered regions with decreasing mass resolution of 107  M⊙, 8  × 107  M⊙, 6 × 108  M⊙, 5 × 109  M⊙ for spheres of xxx radius. In order to follow the Lagrangian evolution of the initial patch, we fill this initial sub-volume with a passive color variable with values of 1 inside and zero outside, and we only allow for refinement when this passive color is above a value of 0.01. Within this ‘colored’ region, refinement is allowed in a quasi-Lagrangian manner down to a resolution of Δx = 34   pc at z = 0: if the total mass in a cell becomes greater than 8 times the initial mass resolution. The minimum cell size is kept roughly constant by adding an extra level of refinement every times the expansion factor is doubling, i.e. at a exp  = 0.1, 0.2, 0.4 and 0.8 (the minimum cell size is thus between Δx =  27 and 54 pc). We have also added a super-Lagrangian refinement criterion to enforce the refinement of the mesh if a cell has a size shorter than one Jean length wherever the gas number density is larger than 5  H cm − 3.
The   NewHorizon  simulation is run with the adaptive mesh refinement ramses code [60]. Gas is evolved with a second-order Godunov scheme and the approximate Harten-Lax-Van Leer-Contact [62] Riemann solver with linear interpolation of the cell-centered quantities at cell interfaces using a MinMod total variation diminishing scheme.

  1.2  Radiative cooling and heating

We adopt the equilibrium chemistry model for primordial species (H and He) assuming collisional ionization equilibrium in the presence of a homogeneous UV background. The primordial gas is allowed to cool down to  ≈  104  K through collisional ionization, excitation, recombination, Bremsstrahlung, and Compton cooling. Metal-enriched gas can cool further down to 0.1  K using rates tabulated by [58] above  ≈ 104  K and those from [18] below  ≈ 104  K. Heating of the gas from a uniform UV background takes place after redshift z reion = 10 following [27]. Motivated by the radiation-hydrodynamic simulation results that the UV background is self-shielded in optically thick regions (n H  ga0.01  H cm − 3) [52], we assume that UV photo-heating rates are reduced by a factor exp(  − n ⁄ n shield), where n shield  = 0.01  H cm  − 3.

  1.3  Star formation

Star formation occurs in regions with hydrogen gas number density above n0  = 10  H cm  − 3 (the stellar mass resolution is n0m pΔx3 = 1.3 × 104  M⊙) following a Schmidt law: ρ̇⋆  = ϵ⋆ρ g  ⁄ t ff, where ρ̇⋆ is the star formation rate mass density, ρ g the gas mass density, t ff = √3π  ⁄ (32Gρ g) the local free-fall time of the gas, G the gravitational constant, and ϵ⋆ is a varying star formation efficiency [37, 63, 64].
The theory of star formation provides a framework for working out the efficiency of the star formation where the gas density probability density function (PDF) is well approximated by a log-normal PDF [41, 49, 15, 24]. This PDF is related to the star forming cloud properties through the cloud turbulent Mach number ℳ = u rms ⁄ c s and virial parameter α vir = 2E kin  ⁄ E grav, and the efficiency is fully determined by integrating how much mass passes above a given density threshold using the multi-free fall approach of [29]:
( 1) ϵ⋆  = ϵ2φ texp⎛⎝38σ2s⎞⎠⎡⎣1  +  erf⎛⎝σ2s  − s crit√2σ2s⎞⎠⎤⎦,  
where s = ln(ρ ⁄ ρ0) is the logarithmic density contrast of the PDF with mean ρ0 and variance σ2 s = ln(1 + b2ℳ2). b = 0.4 conveys the fractional amount of solennoidal to compresionnal modes of the turbulence. The critical density contrast s crit is determined by [49]
( 2) s crit = ln⎛⎝0.067θ2α virℳ2⎞⎠.
In the   NewHorizon simulation, the turbulent Mach number is given by the local three-dimensional instantaneous velocity dispersion σ g, and the virial parameter takes also into account the thermal pressure support α vir, 0 = 5(σ2 g + c2 s) ⁄ (πρ gGΔx2). φ  − 1 t  = 0.57 and θ = 0.33 are empirical parameters of the model determined by the best-fit values between the theory and the numerical experiments [24]  [A] . We ignore the role of the magnetic field in this model despite the effect it has on the critical density and variance of the density PDF due to its large pressure with respect to the thermal pressure in the cold neutral medium [28, 17]. ϵ = 0.5 is a proto-stellar feedback parameter that controls the actual amount of gas above s crit that is able to form stars [44, 2, 4].
Such a star formation law shows a significantly different behaviour on galactic scales to simulations with constant (usually low) efficiencies since the efficiency can now vary by orders of magnitude between gravitationally bound α vir < 1 and highly turbulent regions ℳ > 1 where the efficiency can go well above 1 and regions that are marginally bound where the efficiency quickly drops to very low values. Star formation efficiency, in conjunction with stellar feedback, plays a key role in shaping galaxy properties [1, 47], and such potentially higher and more bursty star formation participates to drive stronger outflows and self-regulation of galaxy properties.

  1.4  Feedback from massive stars

We include feedback from Type II supernova (SN) assuming that each explosion initially releases the kinetic energy of 1051  erg. Because the minimum mass of a star particle is 104  M⊙, each particle is assumed to represent a simple stellar population with a Chabrier initial mass function (IMF) [13] where the lower (upper) mass cutoff is taken as M low = 0.1 (M upp = 150)  M⊙, respectively. We further assume that the minimum mass that explodes is 6  M⊙ in order to include electron-capture SNe [15]. The corresponding specific frequency of SN explosion is 0.015  M  − 1⊙. We increase this number by a factor of two (0.03  M  − 1⊙), because multiple SN explosions can increase the total radial momentum by making the ambient density into which subsequent SNe explode lower [33, 26]. SNe are assumed to explode instantaneously when a star particle becomes older than 5  Myr. The mass loss fraction from the explosions is 31%, with the metal yield (mass ratio of the newly formed metals over the total ejecta) of 0.05.
We employ the mechanical SN feedback scheme [35, 34] which ensures to transfer a correct amount of radial momentum to the surroundings. Specifically, the model examines whether the blast wave is in the Sedov-Taylor energy-conserving or momentum-conserving phase [14, 16, 9] by calculating the mass swept up by SN. If the SN explosion is still in the energy-conserving phase, the assumed specific energy is injected to the gas since hydrodynamics will naturally capture the expansion of the SN and imparts the correct amount of radial momentum. However, if the cooling length in the neighbouring regions is under-resolved due to a finite resolution, radiative cooling takes place rapidly, suppressing the expansion of the SN bubble. This leads to a under-estimation of the radial momentum, hence weaker feedback. In order to avoid this artificial cooling, the mechanical feedback model directly imparts the radial momentum expected during the momentum-conserving phase if the mass of the neighbouring cell exceeds some critical value. This is done by first measuring the local ratio of the swept-up gas mass over the ejecta mass, and examining if the ratio is greater than the critical ratio corresponding to the energy-to-momentum phase transition, i.e. 70 E − 2 ⁄ 17 51 n  − 4 ⁄ 17 1 Z’ − 0.28, where E 51 is the total energy released in units of 1051  erg, n1 is the hydrogen number density in units of  cm  − 3, and Z’ =  max[Z ⁄ Z⊙, 0.01] is the metallicity, normalised to the solar value (Z⊙  = 0.02). The final momentum in the snowplough phase per SN explosion is taken from [61] as:
( 3) q SN = 3 × 105  km s − 1M⊙ E16  ⁄ 1751n  − 2 ⁄ 17 1Z’  − 0.14 .
We further assume that the UV radiation from the young OB stars over-pressurize the ambient medium near young stars and increase the total momentum per SN to
( 4) q SN + PH = 5 × 105  km s − 1M⊙ E16  ⁄ 1751n  − 2 ⁄ 17 1Z’  − 0.14, 
following [25].
It is worth noting that the specific energy used for SN II explosion in this study is larger than previously assumed. A [12] IMF with a low-/high-mass cut-off of M low = 0.1 and M upp = 100  M⊙ and an intermediate-to-massive star transition mass at M IM = 8  M⊙ gives e SN  = 1.1 × 1049  erg M − 1⊙. However, e SN can be increased up to 3.6  × 1049  erg M − 1⊙ if a non-negligible fraction (f HN = 0.5) of hypernovae [32, 46] is taken into account, which is necessary to reproduce the abundance of heavy elements, such as Zinc [39], or if a lower transition mass M IM = 6  M⊙ and a shallower (Salpeter) slope of  − 2.1 at the high-mass end [65, 10, 43]. Furthermore, various sources of stellar feedback which would contribute to the overall formation of large-scale outflows including type Ia SNe, stellar winds, shock-accelerated cosmic rays [66, 55, 19], infrared multi-scattering on dust [30, 54, 53] or Lyman-α resonant line scattering [36, 57] are neglected. In addition runaway OB stars [11, 35, 3] or the unresolved porosity of the medium [31] are effects that are also ignored. In this regard, the   NewHorizon simulation is unlikely to over-estimate the effects of stellar feedback, as we discuss in more detail in Section 3.

  1.5  MBHs and AGN

  1.5.1  MBH formation, growth and dynamics

In   NewHorizon, massive black holes (MBHs) are assumed to form in cells that have gas and stellar densities above the threshold for star formation, a stellar velocity dispersion larger than 20  km s − 1 and that are located at a distance of at least 50 comoving kpc from any pre-existing MBH.
Once formed, the mass of MBHs grows at a rate Ṁ MBH = (1 − ϵ r)Ṁ Bondi, where ϵr is the spin-dependent radiative efficiency (see Eq:  7↓) and Ṁ Bondi is the Bondi-Hoyle-Lyttleton rate,
( 5) dM Bondidt  = 4πρ(GM MBH)2(u2  + c2s)3  ⁄ 2 , 
with G the gravitational constant, u the average MBH-to-gas relative velocity, c s the average gas sound speed, and ρ the average gas density. All average quantities are computed within 4Δx of the MBH, using mass weighting and a kernel weighting as specified in [21]. Note that we do not employ a boost factor in the formulation of the accretion rate, as is commonly done in cosmological simulations, as we have sufficient spatial resolution to model some of the multiphase structure of the interstellar medium of galaxies directly.
The Bondi-Hoyle-Lyttleton accretion rate is capped at the Eddington luminosity rate for the appropriate ϵ r
( 6) dM Edddt  = 4πGM MBHm pϵ rσTc , 
where σT is the Thompson cross-section, mp the proton mass and c the speed of light.
To avoid spurious motions of MBHs around high gas density regions due to finite force resolution effects, we include an explicit drag force of the gas onto the MBH, following [48]. This drag force term includes a boost factor with the functional form α = (n ⁄ n0)2 when n > n0, and α = 1 otherwise. We also enforce maximum refinement within a region of radius 4Δx around the MBH, which improves the accuracy of MBH motions [42]. MBHs are allowed to merge when they get closer than 4Δx ( ~ 150 pc) and when the relative velocity of the pair is smaller than the escape velocity of the binary. A detailed analysis of MBH mergers in   NewHorizon is presented in [67].

  1.5.2  MBH spin evolution

The evolution of the spin parameter a is followed on-the-fly in the simulation, taking into account the effects of gas accretion and MBH-MBH mergers. The model of MBH spin evolution has been introduced in [23], to which we refer the reader for the technical details. The only change is that we now use a different MBH spin evolution model at low accretion rates: χ = Ṁ MBH ⁄ M Edd  < χ trans, where χ trans  = 0.01. At high accretion rates (χ ≥  χ trans), a thin accretion disc solution is assumed [56], as in [23], and the angular momentum direction of the accreted gas is used to decide whether the accreted gas feeds an aligned or misaligned Lense-Thirring disc precessing with the spin of the MBH [38], respectively spinning the MBH up or down. At low accretion rates (χ  < 0.01), we assume that jets are powered by energy extraction from MBH rotation [8], and that the MBH spin magnitude can only decrease. The change in the spin magnitude da ⁄ dM follows the results from [45], where we have fitted a fourth-order polynomial to their sampled values (from their table 7, sH for AaN100 runs, where a is the value of the MBH spin). The functional form of the spin evolution as a function of MBH spin at low accretion rates is represented in the top panel of Fig. ↓, where the dimensionless spin up parameter is shown: s  ≡ d(a ⁄ M) ⁄ dt, where if s and a have opposite sign the black hole spins down.
In addition, MBH spins change in magnitude and direction during MBH-MBH coalescences, with the spin of the remnant depending on the spins of the two merging MBHs and the orbital angular momentum of the binary, following analytical expressions from [51].
The evolution of the spin parameter is a key component of the AGN feedback model, as it controls the radiative efficiency of the accretion disc and the jet efficiency and. Therefore, the Eddington mass accretion rate, here used to cap the total accretion rate, and the AGN feedback efficiency in both jet and thermal mode will vary with spin values. The spin-dependent radiative efficiency is defined as
( 7) ϵ r  = f att(1  − E isco)  = f att(1  − √1  − 2 ⁄ (3r isco))
where r isco  = R isco ⁄ R g is the radius of the innermost stable circular orbit (ISCO) in reduced units and R g is half the Schwarzschild radius of the MBH. R isco depends on spin a. For the radio mode, the radiative efficiency used in the effective growth of the MBH is attenuated by a factor f att = min(χ ⁄ χ trans,  1) following [6]. MBHs seeds are initialised with a zero spin value.

  1.5.3  Radio and quasar modes of AGN feedback

AGN feedback is modelled in two different ways depending on the Eddington rate [21]: below χ < χ trans the MBH powers jets (a.k.a. radio mode) releasing mass, momentum and total energy into the gas [20], while above χ ≥ χ trans the MBH releases only thermal energy back into the gas (a.k.a. quasar mode, [59]). For the two different modes, the AGN releases a power that is a fraction of the rest-mass accreted luminosity on to the MBH, L AGN, R, Q = η R, QṀ MBHc2.
For the jet mode of AGN feedback, the efficiency η R is not a free parameter: it scales with the MBH spin, following the results from magnetically chocked accretion discs of [45], where we have fitted a fourth-order polynomial to their sampled values of jet plus wind efficiencies (from their table 5, ηj plus ηw, 0 for runs AaN100). When active in our simulation, the bipolar AGN jet deposits mass, momentum, and total energy within a cylinder of size Δx in radius and semi-height, centred on the MBH, whose axis is (anti)aligned with the MBH spin axis (zero opening angle). Jets are launched with a speed of 104  km s − 1.
The quasar mode of AGN feedback deposits internal energy into its surrounding within a sphere of radius Δx, within which the specific energy is uniformly deposited (uniform temperature increase). As only a fraction of the AGN driven wind is expected to thermalize, and only some of the multi-wavelength radiation emitted from the accretion disc will couple to the gas on ISM scales  [7], we scale the feedback efficiency in quasar mode by a coupling factor of ηc = 0.15 [21]. The effective feedback efficiency in quasar mode is therefore η Q = ϵrηc.

  1.6  Identification of halos and galaxies

Halos are identified with the AdaptaHOP halo finder [5]. The density field used in AdaptaHOP is smoothed over 20 particles. The minimum number of particles in a halo is 100 DM particles, while only halos with an average density larger than 80 times the critical density are considered. The center of the halo is recursively determined by seeking the centre of mass in a shrinking sphere, while decreasing its radius by 10 per cent recurrently down to a minimum radius of 0.5 kpc [50]. The maximum DM density in that radius is defined as the center of the halo. The shrinking sphere approach is used since strong feedback processes can significantly flatten the central DM density and smaller, but denser, substructures can be misidentified as being the center of the main halo.
The same identification technique, using either AdaptaHop or HOP, has been running onto stars to identify the galaxies in the simulation, except that we only consider galaxies with more than 50 star particles. The former solution removes substructures which include both in situ star forming clumps as well as satellites already connected to a galaxy. The latter keeps all substructures connected to the main structure. Both will be employed depending on context, as indicated in the corresponding text. For the centering the galaxies at the low mass end, particular attention has to be taken , since these galaxies tend to be extremely turbulent structures where bulges cannot be easily identified.
Since the   NewHorizon  simulation is a zoom simulation embedded in a larger cosmological volume filled with lower DM resolution particles, one also needs to remove halos of the zoom regions polluted with low resolution DM particles. To that end, we only consider halos, and their embedded galaxies and MBHs encompassed in their virial radius, which are found devoid of low resolution DM particles up to some threshold. With 100 % purity, there are 6500 uncontaminated galaxies in the whole galaxy sample with stellar mass larger than 5  × 105  M⊙ at z = 1.7, 404 and 103 with mass larger than 108  M⊙ and 109  M⊙ respectively.


References

[1] Oscar Agertz, Romain Teyssier, Ben Moore: “The formation of disc galaxies in a ΛCDM universe”,   mnras, pp. 1391-1408, 2011.

[2] J. Alves, M. Lombardi, C.~J. Lada: “The mass function of dense molecular cores and the origin of the IMF”,   aap, pp. L17-L21, 2007.

[3] Eric P. Andersson, Oscar Agertz, Florent Renaud: “How runaway stars boost galactic outflows”,   mnras, 2020.

[4] Ph. André, A. Men'shchikov, S. Bontemps, V. Könyves, F. Motte, N. Schneider, P. Didelon, V. Minier, P. Saraceno, D. Ward-Thompson, J. di Francesco, G. White, S. Molinari, L. Testi, A. Abergel, M. Griffin, Th. Henning, P. Royer, B. Merı́n, R. Vavrek, M. Attard, D. Arzoumanian, C.~D. Wilson, P. Ade, H. Aussel, J. -P. Baluteau, M. Benedettini, J. -Ph. Bernard, J.~A.~D.~L. Blommaert, L. Cambrésy, P. Cox, A. di Giorgio, P. Hargrave, M. Hennemann, M. Huang, J. Kirk, O. Krause, R. Launhardt, S. Leeks, J. Le Pennec, J.~Z. Li, P.~G. Martin, A. Maury, G. Olofsson, A. Omont, N. Peretto, S. Pezzuto, T. Prusti, H. Roussel, D. Russeil, M. Sauvage, B. Sibthorpe, A. Sicilia-Aguilar, L. Spinoglio, C. Waelkens, A. Woodcraft, A. Zavagno: “From filamentary clouds to prestellar cores to the stellar IMF: Initial highlights from the Herschel Gould Belt Survey”,   aap, pp. L102, 2010.

[5] D. Aubert, C. Pichon, S. Colombi: “The origin and implications of dark matter anisotropic cosmic infall on ̃L* haloes”,   mnras, pp. 376-398, 2004.

[6] A.~J. Benson, A. Babul: “Maximum spin of black holes driving jets”,   mnras, pp. 1302-1313, 2009.

[7] R. Bieri, Y. Dubois, J. Rosdahl, A. Wagner, J. Silk, G.~A. Mamon: “Outflows Driven by Quasars in High-Redshift Galaxies with Radiation Hydrodynamics”,   mnras, 2016.

[8] R.~D. Blandford, R.~L. Znajek: “Electromagnetic extraction of energy from Kerr black holes”, MNRAS, pp. 433-456, 1977.

[9] John M. Blondin, Eric B. Wright, Kazimierz J. Borkowski, Stephen P. Reynolds: “Transition to the Radiative Phase in Supernova Remnants”,   apj, pp. 342-354, 1998.

[10] Michele Cappellari, Richard M. McDermid, Katherine Alatalo, Leo Blitz, Maxime Bois, Frédéric Bournaud, M. Bureau, Alison F. Crocker, Roger L. Davies, Timothy A. Davis, P.~T. de Zeeuw, Pierre-Alain Duc, Eric Emsellem, Sadegh Khochfar, Davor Krajnović, Harald Kuntschner, Pierre-Yves Lablanche, Raffaella Morganti, Thorsten Naab, Tom Oosterloo, Marc Sarzi, Nicholas Scott, Paolo Serra, Anne-Marie Weijmans, Lisa M. Young: “Systematic variation of the stellar initial mass function in early-type galaxies”,   nat, pp. 485-488, 2012.

[11] Daniel Ceverino, Anatoly Klypin: “The Role of Stellar Feedback in the Formation of Galaxies”,   apj, pp. 292-309, 2009.

[12] Gilles Chabrier: “Galactic Stellar and Substellar Initial Mass Function”,   pasp, pp. 763-795, 2003.

[13] Gilles Chabrier: The Initial Mass Function: From Salpeter 1955 to 2005. 2005.

[14] Roger A. Chevalier: “The Evolution of Supernova Remnants. Spherically Symmetric Models”,   apj, pp. 501-516, 1974.

[15] Cesare Chiosi, Gianpaolo Bertelli, Alessandro Bressan: “New developments in understanding the HR diagram”,   araa, pp. 235-285, 1992.

[16] Denis F. Cioffi, Christopher F. McKee, Edmund Bertschinger: “Dynamics of Radiative Supernova Remnants”,   apj, pp. 252, 1988.

[17] Richard M. Crutcher: “Magnetic Fields in Molecular Clouds”,   araa, pp. 29-63, 2012.

[18] A. Dalgarno, R.~A. McCray: “Heating and Ionization of HI Regions”,   araa, pp. 375, 1972.

[19] Gohar Dashyan, Yohan Dubois: “Cosmic ray feedback from supernovae in dwarf galaxies”, arXiv e-prints, pp. arXiv:2003.09900, 2020.

[20] Y. Dubois, J. Devriendt, A. Slyz, R. Teyssier: “Jet-regulated cooling catastrophe”,   mnras, pp. 985-1001, 2010.

[21] Y. Dubois, J. Devriendt, A. Slyz, R. Teyssier: “Self-regulated growth of supermassive black holes by a dual jet-heating active galactic nucleus feedback mechanism: methods, tests and implications for cosmological simulations”,   mnras, pp. 2662-2683, 2012.

[22] Y. Dubois, C. Pichon, C. Welker, D. Le Borgne, J. Devriendt, et al.: “Dancing in the dark: galactic properties trace spin swings along the cosmic web”,   mnras, pp. 1453-1468, 2014.

[23] Y. Dubois, M. Volonteri, J. Silk, J. Devriendt, A. Slyz: “Black hole evolution - II. Spinning black holes in a supernova-driven turbulent interstellar medium”,   mnras, pp. 2333-2346, 2014.

[24] Christoph Federrath, Ralf S. Klessen: “The Star Formation Rate of Turbulent Magnetized Clouds: Comparing Theory, Simulations, and Observations”,   apj, pp. 156, 2012.

[25] S. Geen, J. Rosdahl, J. Blaizot, J. Devriendt, A. Slyz: “A detailed study of feedback from a massive star”,   mnras, pp. 3248-3264, 2015.

[26] Eric S. Gentry, Mark R. Krumholz, Piero Madau, Alessandro Lupi: “The momentum budget of clustered supernova feedback in a 3D, magnetized medium”,   mnras, pp. 3647-3658, 2019.

[27] F. Haardt, P. Madau: “Radiative Transfer in a Clumpy Universe. II. The Ultraviolet Extragalactic Background”,   apj, pp. 20-+, 1996.

[28] Carl Heiles, T.~H. Troland: “The Millennium Arecibo 21 Centimeter Absorption-Line Survey. IV. Statistics of Magnetic Field, Column Density, and Turbulence”,   apj, pp. 773-793, 2005.

[29] Patrick Hennebelle, Gilles Chabrier: “Analytical Star Formation Rate from Gravoturbulent Fragmentation”,   apjl, pp. L29, 2011.

[30] Philip F. Hopkins, Eliot Quataert, Norman Murray: “Self-regulated star formation in galaxies via momentum input from massive stars”,   mnras, pp. 950-973, 2011.

[31] Olivier Iffrig, Patrick Hennebelle: “Mutual influence of supernovae and molecular clouds”,   aap, pp. A95, 2015.

[32] K. Iwamoto, P.~A. Mazzali, K. Nomoto, H. Umeda, T. Nakamura, F. Patat, I.~J. Danziger, T.~R. Young, T. Suzuki, T. Shigeyama, T. Augusteijn, V. Doublier, J. -F. Gonzalez, H. Boehnhardt, J. Brewer, O.~R. Hainaut, C. Lidman, B. Leibundgut, E. Cappellaro, M. Turatto, T.~J. Galama, P.~M. Vreeswijk, C. Kouveliotou, J. van Paradijs, E. Pian, E. Palazzi, F. Frontera: “A hypernova model for the supernova associated with the γ-ray burst of 25 April 1998”,   nat, pp. 672-674, 1998.

[33] Chang-Goo Kim, Eve C. Ostriker, Roberta Raileanu: “Superbubbles in the Multiphase ISM and the Loading of Galactic Winds”,   apj, pp. 25, 2017.

[34] T. Kimm, R. Cen, J. Devriendt, Y. Dubois, A. Slyz: “Towards simulating star formation in turbulent high-z galaxies with mechanical supernova feedback”,   mnras, pp. 2900-2921, 2015.

[35] Taysun Kimm, Renyue Cen: “Escape Fraction of Ionizing Photons during Reionization: Effects due to Supernova Feedback and Runaway OB Stars”,   apj, pp. 121, 2014.

[36] Taysun Kimm, Martin Haehnelt, Jérémy Blaizot, Harley Katz, Léo Michel-Dansac, Thibault Garel, Joakim Rosdahl, Romain Teyssier: “Impact of Lyman alpha pressure on metal-poor dwarf galaxies”,   mnras, pp. 4617-4635, 2018.

[37] Taysun Kimm, Harley Katz, Martin Haehnelt, Joakim Rosdahl, Julien Devriendt, Adrianne Slyz: “Feedback-regulated star formation and escape of LyC photons from mini-haloes during reionization”,   mnras, pp. 4826-4846, 2017.

[38] A.~R. King, S.~H. Lubow, G.~I. Ogilvie, J.~E. Pringle: “Aligning spinning black holes and accretion discs”,   mnras, pp. 49-56, 2005.

[39] Chiaki Kobayashi, Hideyuki Umeda, Ken'ichi Nomoto, Nozomu Tominaga, Takuya Ohkubo: “Galactic Chemical Evolution: Carbon through Zinc”,   apj, pp. 1145-1171, 2006.

[40] E. Komatsu, K.~M. Smith, J. Dunkley, C.~L. Bennett, B. Gold, G. Hinshaw, N. Jarosik, D. Larson, M.~R. Nolta, L. Page, D.~N. Spergel, M. Halpern, R.~S. Hill, A. Kogut, M. Limon, S.~S. Meyer, N. Odegard, G.~S. Tucker, J.~L. Weiland, E. Wollack, E.~L. Wright: “Seven-year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Cosmological Interpretation”,   apjs, pp. 18-+, 2011.

[41] Mark R. Krumholz, Christopher F. McKee: “A General Theory of Turbulence-regulated Star Formation, from Spirals to Ultraluminous Infrared Galaxies”,   apj, pp. 250-268, 2005.

[42] Alessandro Lupi, Francesco Haardt, Massimo Dotti: “Massive black hole and gas dynamics in galaxy nuclei mergers - I. Numerical implementation”,   mnras, pp. 1765-1774, 2015.

[43] Ignacio Martı́n-Navarro, Alexandre Vazdekis, Francesco La Barbera, Jesús Falcón-Barroso, Mariya Lyubenova, Glenn van de Ven, Ignacio Ferreras, S.~F. Sánchez, S.~C. Trager, R. Garcı́a-Benito, D. Mast, M.~A. Mendoza, P. Sánchez-Blázquez, R. González Delgado, C.~J. Walcher, The CALIFA Team: “IMF&ndashMetallicity: A Tight Local Relation Revealed by the CALIFA Survey”,   apjl, pp. L31, 2015.

[44] Christopher D. Matzner, Christopher F. McKee: “Efficiencies of Low-Mass Star and Star Cluster Formation”,   apj, pp. 364-378, 2000.

[45] J.~C. McKinney, A. Tchekhovskoy, R.~D. Blandford: “General relativistic magnetohydrodynamic simulations of magnetically choked accretion flows around black holes”,   mnras, pp. 3083-3117, 2012.

[46] Ken'ichi Nomoto, Nozomu Tominaga, Hideyuki Umeda, Chiaki Kobayashi, Keiichi Maeda: “Nucleosynthesis yields of core-collapse supernovae and hypernovae, and galactic chemical evolution”,   nphysa, pp. 424-458, 2006.

[47] Arturo Nuñez-Castiñeyra, Emmanuel Nezri, Julien Devriendt, Romain Teyssier: “Cosmological simulations of the same spiral galaxy: the impact of baryonic physics”, arXiv e-prints, pp. arXiv:2004.06008, 2020.

[48] E.~C. Ostriker: “Dynamical Friction in a Gaseous Medium”, ApJ, pp. 252-258, 1999.

[49] Paolo Padoan, Åke Nordlund: “The Star Formation Rate of Supersonic Magnetohydrodynamic Turbulence”,   apj, pp. 40, 2011.

[50] C. Power, J.~F. Navarro, A. Jenkins, C.~S. Frenk, S.~D.~M. White, V. Springel, J. Stadel, T. Quinn: “The inner structure of ΛCDM haloes - I. A numerical convergence study”,   mnras, pp. 14-34, 2003.

[51] L. Rezzolla, E. Barausse, E.~N. Dorband, D. Pollney, C. Reisswig, J. Seiler, S. Husa: “Final spin from the coalescence of two black holes”,   prd, pp. 044002, 2008.

[52] J. Rosdahl, J. Blaizot: “Extended Lyα emission from cold accretion streams”,   mnras, pp. 344-366, 2012.

[53] J. Rosdahl, R. Teyssier: “A scheme for radiation pressure and photon diffusion with the M1 closure in RAMSES-RT”,   mnras, pp. 4380-4403, 2015.

[54] Rok Roškar, Romain Teyssier, Oscar Agertz, Markus Wetzstein, Ben Moore: “A systematic look at the effects of radiative feedback on disc galaxy formation”,   mnras, pp. 2837-2853, 2014.

[55] Munier Salem, Greg L. Bryan: “Cosmic ray driven outflows in global galaxy disc models”,   mnras, pp. 3312-3330, 2014.

[56] N.~I. Shakura, R.~A. Sunyaev: “Black holes in binary systems. Observational appearance.”,  aap, pp. 337-355, 1973.

[57] Aaron Smith, Volker Bromm, Abraham Loeb: “Lyman α radiation hydrodynamics of galactic winds before cosmic reionization”,  mnras, pp. 2963-2978, 2017.

[58] R.~S. Sutherland, M.~A. Dopita: “Cooling functions for low-density astrophysical plasmas”, apjs, pp. 253-327, 1993.

[59] R. Teyssier, B. Moore, D. Martizzi, Y. Dubois, L. Mayer: “Mass distribution in galaxy clusters: the role of Active Galactic Nuclei feedback”, mnras, pp. 195-208, 2011.

[60] R. Teyssier: “Cosmological hydrodynamics with adaptive mesh refinement. A new high resolution code called RAMSES”, aap, pp. 337-364, 2002.

[61] K. Thornton, M. Gaudlitz, H. -Th. Janka, M. Steinmetz: “Energy Input and Mass Redistribution by Supernovae in the Interstellar Medium”,  apj, pp. 95-119, 1998.

[62] E. Toro: Riemann Solvers and Numerical Methods for Fluid Dynamics. Springer-Verlag, 1999.

[63] Maxime Trebitsch, Jérémy Blaizot, Joakim Rosdahl, Julien Devriendt, Adrianne Slyz: “Fluctuating feedback-regulated escape fraction of ionizing radiation in low-mass, high-redshift galaxies”,  mnras, pp. 224-239, 2017.

[64] Maxime Trebitsch, Yohan Dubois, Marta Volonteri, Hugo Pfister, Corentin Cadiou, Harley Katz, Joakim Rosdahl, Taysun Kimm, Christophe Pichon, Ricarda S. Beckmann, Julien Devriendt, Adrianne Slyz: “The Obelisk simulation: galaxies contribute more than AGN to HI reionization of protoclusters”, arXiv e-prints, pp. arXiv:2002.04045, 2020.

[65] Tommaso Treu, Matthew W. Auger, Léon V.~E. Koopmans, Raphaël Gavazzi, Philip J. Marshall, Adam S. Bolton: “The Initial Mass Function of Early-Type Galaxies”,  apj, pp. 1195-1202, 2010.

[66] M. Uhlig, C. Pfrommer, M. Sharma, B.~B. Nath, T.~A. Enßlin, V. Springel: “Galactic winds driven by cosmic ray streaming”,  mnras, pp. 2374-2396, 2012.

[67] Marta Volonteri, Hugo Pfister, Ricarda S. Beckman, Yohan Dubois, Monica Colpi, Christopher J. Conselice, Massimo Dotti, Garreth Martin, Ryan Jackson, Katarina Kraljic, Christophe Pichon, Maxime Trebitsch, Sukyoung. K. Yi, Julien Devriendt, Sebastien Peirani: “Black hole mergers from dwarf to massive galaxies with the NewHorizon and Horizon-AGN simulations”, arXiv e-prints, pp. arXiv:2005.04902, 2020.