Forming early-type galaxies without AGN feedback: a combination of merger-driven outflows and inefficient star formation
Michael Kretschmer Romain Teyssier





firstpage–lastpage
Robust topological invariants of topological crystalline phases in the presence of impurities
    Ian Mondragon-Shem^1 and Taylor L. Hughes^2
    December 30, 2023
=============================================================================================




Regulating the available gas mass inside galaxies proceeds through a delicate balance between inflows and outflows, but also through the internal depletion of gas due to star formation. At the same time, stellar feedback is the internal engine that powers the strong outflows. Since star formation and stellar feedback are both small scale phenomena, we need a realistic and predictive subgrid model for both. 
We describe the implementation of supernova momentum feedback and star formation based on the turbulence of the gas in the RAMSES code.
For star formation, we adopt the so-called multi-freefall model.
The resulting star formation efficiencies can be significantly smaller or bigger than the traditionally chosen value of 1%.
We apply these new numerical models to a prototype cosmological simulation of a massive halo 
that features a major merger which results in the formation of an early-type galaxy without using AGN feedback.
We find that the feedback model provides the first order mechanism for regulating the stellar and baryonic content in our simulated galaxy.
At high redshift, the merger event pushes gas to large densities and large turbulent velocity dispersions, such that efficiencies come close to 10%, 
resulting in large SFR. We find small molecular gas depletion time during the starburst, in perfect agreement with observations. 
Furthermore, at late times, the galaxy becomes quiescent with efficiencies significantly smaller than 1%, resulting in small SFR and long molecular gas depletion time.



galaxies: formation –
galaxies: evolution –
galaxies: star formation –
stars: formation –
methods: numerical




§ INTRODUCTION


Although galaxy formation still remains one of the most important unsolved problems in astrophysics,
we have made spectacular progresses in the past decade. 
It became apparent that galaxy formation has to be very inefficient at forming stars, with a peak 
value around 20% of the available baryons for Milky Way sized haloes, down to much smaller efficiencies at lower and larger masses
<cit.>.
One of the key observables that triggered 
our recent advances comes from the abundance matching technique <cit.>,
from which we know relatively precisely the stellar mass of the central galaxy as a function of the parent halo mass.

The inefficiency of galaxy formation led our community to discover the importance of galactic winds to regulate star formation <cit.>, a theoretical idea corroborated
(some might say incepted) by observations of strong outflows at high redshift <cit.>.
In parallel, it was recognised that cosmological infall of fresh gas, in the form of cold streams, was also an important player in setting up the star formation rate
in high redshift galaxies <cit.>.
Galaxy mergers are now believed to play a minor role in the history of gas accretion, while major mergers (mass ratios close to 1:1)
trigger the extreme and rare star bursts we see nearby, as well as in the more distant universe <cit.>.

Regulating the available gas mass inside galaxies proceeds through a delicate balance between inflows and outflows <cit.>,
but also through the internal depletion of gas due to star formation. In the same time, stellar feedback is the internal
engine that powers the strong outflows, and since star formation and stellar feedback are both small scale phenomena, we need a realistic
and predictive subgrid model for both. Although it has been argued that strongly star forming galaxies have their star formation rate regulated by feedback <cit.>,
the transition to quenched galaxies, and their extremely low star formation rate still remains a puzzle. Most galaxies in the low redshift universe are quiescent, even disk
galaxies like the Milky Way. It is therefore of paramount importance to use the right model of star formation to get the correct 
star formation efficiency in this quiescent, quasi-quenched regime. Active Galactic Nuclei (AGN) feedback is invoked to explain the origin of massive, red and dead elliptical galaxies 
<cit.>. 
At intermediate masses, quiescent, spheroidal and dispersion dominated galaxies can be modelled without relying on AGN feedback, but on a 
combination of strong early feedback leading to gas removal and suppression of gas cooling due to a hotter halo temperature <cit.>.
Even in these extremely quiescent galaxies, star formation is still active, but in a quenched, very inefficient mode <cit.>. 

Quantitatively, the efficiency of star formation is characterised by the gas depletion time,
defined as the ratio of the star formation rate to the available gas mass. 
Observationally, the depletion time can be estimated using the so-called global Kennicutt-Schmidt (KS) relation
<cit.>. It turns out that this gas depletion time is very long in nearby galaxies, 
typical around 1 Gyr, with even longer depletion time around 3 Gyr for more massive quiescent galaxies <cit.>.

More information can be obtained using the resolved KS relation <cit.>. One can compute 
the local value of the depletion time in patches of roughly 1 kpc in size. This more local KS relation demonstrates
than the local star formation rate surface density scales with the local gas surface density, with a power law slope
around 1.5 <cit.>. This leads to the simple theoretical idea that the star formation rate
can be modelled locally using a so-called Schmidt law

    ρ̇_⋆=ϵ_ ffρ/t_ff    with   t_ff=√(3π/32 G ρ)

where ρ̇_⋆ is the local star formation rate density, ϵ_ ff is a fixed star formation efficiency
per freefall time, 
ρ is the local gas density and 
t_ff is the local gas freefall time. This simple model naturally explains the observed power law.
It is usually used only if the gas density lies above a fixed threshold corresponding to star forming gas. 
This fixed efficiency can be chosen to match the observed relation in nearby resolved galaxies <cit.>.
The value usually adopted in cosmological simulations of galaxy formation lies around ϵ∼ 0.01 <cit.>. 
Note that the star formation efficiency in these models is by construction uniform throughout the galaxy.

On smaller scales, close to individual Giant Molecular Clouds (GMC) of size between 1 and 10 pc, the situation is not so simple.
Although the average star formation efficiency is again between 1 and 2% <cit.>, it seems that ϵ_ ff can vary significantly 
from cloud to cloud <cit.>. Although this effect can be explained by a spurious
systematic effect in the observational protocol <cit.>, it could also have a physical origin,
related to different internal properties in the star forming molecular clouds. 

Current theories of star formation,
based on self-gravitating supersonic turbulence, favour a scenario where star formation is the result of a turbulent cascade.
In this framework, one can derive a local star formation efficiency by integrating over the log-normal distribution of 
the turbulent gas density <cit.>. In particular, the so-called multi-freefall model 
of star formation <cit.> give promising results, 
with ϵ_ ff depending on two important cloud structural parameters, namely the cloud mean density 
and the cloud mean Mach number (see below for more details). 

Recently, cosmological simulations of galaxy formation started exploring this new approach to model star formation.
For example, <cit.> used the GIZMO code together with a method that requires the gas to be self-gravitating, 
self-shielded and Jeans-unstable to form stars, very close in spirit to the multi-freefall approach.
Using the RAMSES code, several models based on the multi-freefall approach were presented by <cit.> and 
<cit.> with different implementation details.
In an isolated galaxy, <cit.> and <cit.> explored a similar model based on the local star formation efficiency models of <cit.> and <cit.> respectively.
All these different studies revealed that the local star formation efficiency in the simulated galaxies is very inhomogeneous,
and varies widely, both in time and space, similar to what we observe at small scales in molecular clouds.
Additionally, on the larger galactic scales, either globally or using 1 kpc patches, the observed KS relations were also successfully reproduced.

Reproducing the KS relation on large scale could have nothing to do with the details of star formation 
at small scales. Star formation on large scales could be entirely self-regulated by stellar feedback, 
as advocated by <cit.> and <cit.>. 
In a seminal albeit recent paper, <cit.> have shown that once the star formation efficiency at small scale is large enough (typically 
larger than 1%), the global depletion time converges to a fixed value set by stellar feedback only (therefore independent of ϵ_ ff). 

However, for lower star formation
efficiencies (typically smaller than 1%), this is not true anymore and the global depletion time becomes longer and longer, 
inversely proportional to ϵ_ ff.
The transition between these two regimes occurs around 1%. 
This critical value depends on the strength of the adopted feedback model: A weaker stellar feedback model leads to shorter 
gas ejection phases and a shorter depletion time scale. This gives a larger critical value for ϵ_ ff. 
It is therefore of primary importance to use the correct local star formation model
and the correct feedback model to get the right transition between the quiescent regime and the feedback dominated regime.

Numerical recipe for modelling stellar feedback have seen tremendous progresses over the past decade.
Although several stellar evolution processes are believed to inject thermal and kinetic energy into the
surrounding interstellar medium (ISM) of a star or a star cluster <cit.>,
it is now believed that Type II supernovae explosion is the dominant mechanism for stellar feedback. The 
challenge is then to resolve the cooling radius that marks the transition from the energy conserving phase of 
the remnant to its momentum conserving phase. This allows to correctly describe the hot X-ray emitting gas in the ISM
and to accumulate enough momentum to get the right amount of terminal momentum in the final phase. 

The first numerical implementations
were based on suppressing cooling artificially and increasing the cooling radius to the spatial resolution of the simulation
<cit.>.
Although this model is qualitatively successful in modelling strong supernovae feedback, 
it overestimates the hot gas fraction and the terminal momentum of the explosion. It was later proposed
to directly inject the correct thermal energy and terminal momentum into neighbouring cells <cit.>.
This method is now widely adopted, as it seems to capture the right amount of kinetic energy, but as shown by <cit.>, 
it does not capture the hot phase properly.
These different methods were recently compared and cross-calibrated in <cit.>, who concluded that no recipe is truly satisfactory,
but direct injection of the correct terminal momentum seems to be the least bad method of all.

In this paper, we use the RAMSES code, for which several competing supernovae feedback implementations have been developed over the past years.
<cit.> injected directly the correct supernovae terminal momentum, 
using a as non-thermal energy variable and modifying the Riemann solver accordingly. 
<cit.> explored the effects of radiative feedback, within the framework of the <cit.> delayed cooling model. 
Finally, <cit.> implemented a mechanical feedback scheme that directly injects the correct terminal momentum in the surrounding cells, 
based on the earlier model of <cit.>.

In this paper, we combine these various new developments, both for star formation with a varying efficiency
and for supernovae momentum feedback, combining the best of each past implementation in a novel and unique way,
that we believe is superior to what has been done before in the RAMSES code.
We apply this new implementation of subgrid galaxy formation physics to a prototype cosmological simulation of a massive halo 
that features a major merger, and study specifically the interesting case of an early-type galaxy.
We are particularly interested in the effect of our local star formation recipe on the global galaxy properties.
For this, we explore two different local star formation models: a classical Schmidt law with uniform efficiency 
and the multi-freefall model proposed by <cit.>.
We finally discuss why this new model can explain the properties of quenched galaxies,
as a consequence of strong feedback combined with a very inefficient local star formation.


§ NUMERICAL METHODS


We use the RAMSES Adaptive Mesh Refinement (AMR) code to model the dynamical evolution of both dark matter and baryonic fluids
in a cosmological context <cit.> .
AMR allows to divide space into cubical cells that can be adaptively refined to increase the spatial resolution locally
according to some adopted refinement criterion. In this work, we use the so-called quasi-Lagrangian strategy,
for which cells are refined when the dark matter or baryonic mass exceeds a given threshold value, that corresponds to
the mass resolution adopted in the initial conditions (see below for details). For simulations that include baryonic physics, one must impose a maximum level of 
resolution (or a minimal cell size) because of limited computational resources. Current state-of-the-art cosmological galaxy formation
simulations reach a spatial resolution between tens to hundreds of parsecs, barely resolving the vertical thickness of the gaseous disks,
and still missing the scale of the largest molecular clouds in the Galaxy. Only very recently do we see a few examples of “cloud resolving”
galaxy simulations in the cosmological context <cit.>, while isolated, idealised simulations are already capable of resolving molecular clouds
down to sub-pc scales <cit.>. None of these simulations are however able to resolve the truly star-forming scale, namely 
individual subsonic (or mildly transonic) molecular cores below 0.1 pc. 

It is therefore mandatory to augment these numerical models with a subgrid description of the physics of the ISM. Subgrid models of galaxy formation
have been introduced more than 20 years ago to model star formation in cosmological simulations
using the previously discussed Schmidt law <cit.>. 
As explained in the introduction, our community has now realised that we need to go beyond this 
simple recipe to model the physics of star formation. Although a complete theory of star formation is still to be found, we have made tremendous 
progresses in the past couple of decades <cit.>, the key physical ingredient being supersonic turbulence <cit.>. 
Modern subgrid models of star formation and feedback all feature, with different levels of sophistication, a model for supersonic turbulence. 
We first describe our approach of the problem. We then explain how this subgrid turbulence model is coupled to a physically motivated recipe for star formation and stellar feedback.
 


 §.§ Subgrid model for turbulence



Turbulence can be modelled numerically using the Navier-Stokes equations while resolving the microscopic viscous diffusion scale.
This approach, called Direct Numerical Simulations (DNS) in engineering applications, is not realistic for galaxy formation,
as under normal ISM conditions the viscous scales are many orders of magnitude smaller than any achievable resolution. 
In the 60's, <cit.> proposed a model (later called Large Eddy Simulations or LES) for which the small scales are averaged over to
define the mean flow. The fluid equations are modified through this averaging procedure, introducing the turbulent pressure and several turbulent diffusion terms.
The LES approach allows to model turbulent flows without resolving the dissipative scale, at the expense of designing a Sub-Grid Scale (SGS) model to describe turbulent
effects at the macroscopic scale. Note that numerical dissipation already naturally provides an implicit SGS model for turbulence, sometimes referred to as implicit LES models.

In astrophysics, LES and SGS models have been introduced by <cit.> and <cit.> to study Type Ia supernovae turbulent flame combustion. 
The subgrid turbulence is used to compute the turbulent flame speed properly, accounting for unresolved eddies in the reacting flow. 
Later, the formalism was extended by <cit.> to supersonic turbulent flows and applied to the physics of the ISM. 
Finally, in the context of galaxy formation, <cit.> introduced a very similar SGS model, coupled to a subgrid star formation model, in the spirit of what we present here. 

The difficulty with supersonic turbulence is to model both velocity and density fluctuations. The mass density is decomposed into
the average, large scale density ρ, defined as the volume-averaged density field, smoothed at the scale of the cell Δ x
and the fluctuation ρ', so that the subgrid density writes as ρ = ρ + ρ'. The average velocity 
field and average temperature are however defined using a mass-weighted average (also called the Favre average) with

    v = ρ v/ρ    and   T = ρ T/ρ,
        with   v = v+v”    and   T=T+T”,

where the fluctuations are now defined relative to the Favre average with a double prime. The turbulent kinetic energy is finally defined as

    K_T=1/2ρ v”^2 = 1/2ρσ_ 3D^2,

where we introduce the turbulent three-dimensional velocity dispersion σ_ 3D=√(3)σ_ 1D.
Using these definitions, it is possible to derive new fluid equations for the mean flow variables ρ, v and T,
as well as a new equation for the turbulent kinetic energy. For details on the derivation and the complete form of these new equations, 
we refer to <cit.> and <cit.>. Note that in this paper, contrary to the strategy adopted in <cit.>,
we do not consider the modified form of the Euler equations. We follow the philosophy of <cit.> and rely on numerical diffusion
to provide an implicit LES model, without adding extra diffusion to an already too diffusive numerical approach. We only consider the extra equation
on the turbulent kinetic energy, that writes <cit.>

    tK_ T + x_j( K_ Tv_j) + P_ Tv_jx_j = C_ T - D_ T,

where the turbulent pressure P_ T=2/3 K_ T = ρσ_ 1D^2.
The creation term C_ T is prescribed in the eddy viscosity model 
(also called mixing length theory in different contexts) using the mean flow viscous stress as

    C_ T = 2 μ_ T∑_ij[ 1/2( v_ix_j+v_jx_i)
    - 1/3( ∇· v) δ_ij]^2 = 1/2μ_ T| S_ij|^2.

The destruction term D_ T represents the dissipation of turbulence in the subgrid turbulent cascade down to viscous scales.
It is modelled as

    D_ T = K_ T/τ_ diss.

The two important parameters in the SGS turbulence model are thus the turbulent viscosity parameter μ_ T
and the turbulent dissipation time scale τ_ diss. They are both related to the cell size (also called
in different contexts the smoothing length or the mixing length) by

    μ_ T = ρΔ x σ_ 1D    and   τ_ diss = Δ x/σ_ 1D.

Note the strong analogy with the Chapman-Enskog theory for viscous rarefied gases, except that individual colliding particles
are replaced here by viscous eddies.

Other recent implementations of subgrid models for star formation use an instantaneous estimate of the turbulent velocity dispersion 
<cit.>. This approach can be interpreted here as the stationary limit of 
our turbulent kinetic energy equation, for which

    C_ T = D_ T⟹σ_ 1D = Δ x/√(3)√(| S_ij|^2).

Solving the full turbulent kinetic energy equation allows to account for advection and work of turbulent pressure, as well as 
non-equilibrium dissipation of turbulent kinetic energy. There is however an important caveat in the SGS approach: The
turbulent creation term is modelling injection of turbulent kinetic energy through shear flows. In the absence of gravity,
shear flows are indeed always unstable and turbulent, owing to the famous Kelvin-Helmholtz instability. In presence of gravity, however, this not necessarily true anymore,
as gravity might stabilise the flow, owing for non-convective conditions. On the other hand, self-gravity in a marginally stable disk
might be an additional source of turbulence that we do not consider explicitly in our model.
Our SGS model might overestimate the amount of turbulence, especially in an equilibrium, centrifugally supported disk,
or underestimate it in a gravitationally unstable disk.
This is why we will explore alternative models of creation of turbulence in the Results section.



 §.§ Subgrid model for star formation



Once we know the turbulent kinetic energy in each cell, we can prolong the turbulent spectrum to smaller,
unresolved scales ℓ < Δ x using Burgers turbulence spectrum, for which

    σ(ℓ) = σ_ 1D( ℓ/Δ x)^1/2.

A critical unresolved scale is the sonic length, defined by σ(ℓ_s) = c_s and given be
ℓ_s = Δ x / ℳ^2,
where the cell Mach number is defined by ℳ = σ_ 1D/c_s.
These scales correspond to the subsonic (or mildly transonic) molecular cores,
in which stars form. Below the sonic scale, density fluctuations become very weak. 
Each molecular core can thus be seen as a quasi-homogeneous region of space that 
will eventually collapse and form a star. On scales larger than the sonic length, however,
density fluctuations are very large.
Following <cit.>, but adapting here slightly their methodology, we assume that
gas density distribution in this supersonic turbulent medium follows a log-normal PDF, with

    p(s)=1/√(2πσ_s^2)exp(s-s)^2/2σ_s^2

with the logarithmic density s=ln (ρ / ρ) where ρ is the local density and ρ the mean density of the cell.
The distribution is normalised such that 

    ∫_-∞^+∞ p(s) ds = 1 
        and   ∫ _-∞^+∞ρ p(s) ds = ρ.

The mean logarithmic density s=-1/2 σ_s^2 is related to the standard deviation σ_s which can be fitted, using non-magnetised,
isothermal turbulence simulations <cit.>, by

    σ_s^2 = ln( 1 + b^2 ℳ^2),

where b is a parameter related to the exact nature of the turbulence forcing (solenoidal or compressive).

Following the model of <cit.> for star formation and assuming that these homogeneous cores
are spherically symmetric with diameter ℓ_s, we can compute their virial parameter as

    α_s = 2 E_kin/|E_grav| = 15/πc_s^2+σ(ℓ_s)^2/G ρℓ_s^2 .

A molecular core will collapse and form stars if α_s < 1. This can be translated into a critical density threshold for star formation

    ρ > ρ_ crit = 15/π2 c_s^2 ℳ^4/G Δ x^2 = α_ virρ2 ℳ^4/1+ℳ^2,

where we defined the virial parameter of the whole cell as

    α_ vir = 15/πc_s^2+σ_ 1D^2/G ρΔ x^2 = 15/πc_s^2 /G ρΔ x^2 (1+ℳ^2).

Finally, we derive the lognormal critical density for star formation as

    s_crit = ln[ α_vir2 ℳ^4/1+ℳ^2].


In the multi-freefall models of <cit.>, the local star formation rate is expressed as

    ρ̇_* = ∫_s_crit^∞ρ/t_ff(ρ) p(s)  s = ϵ_ ffρ/t_ ff(ρ).

This formulation just states that each fluid element that satisfies the gravitational instability criterion collapses
in one free-fall time and converts all its mass into (one or several) stars. From this, we can deduce the local star formation efficiency

    ϵ_ff   = ∫_s_crit^∞t_ff(ρ)/t_ff(ρ)ρ/ρ p(s)  s 
       = 1/2exp( 3/8σ_s^2)[ 1 + erf( σ_s^2-s_crit/√(2 σ_s^2)) ].

According to <cit.>, the model breaks down for ℳ≤ 2. 
Indeed, for ℳ = 1, the sonic length becomes equal to the cell size. 
In the simulations we present here, we have most of the time ℳ≥ 10,
but not all the time. Many cells can have subsonic turbulence, especially in the hot gas phase.
We need to provide a model that is valid also for low Mach numbers.
We thus modify the collapse criterion for ℳ < 1, requiring now the whole cell to be gravitationally unstable.

    α_s = 15/πc_s^2+σ_ 1D^2/G ρΔ x^2 < 1.

Note that we didn't use the mean cell density but the local density ρ in the previous equation, so now the density PDF 
plays the role of a probability for the entire cell to have a certain density, while ρ is just the expectancy. We finally have

    s_crit = ln[ α_vir]    for   ℳ≤ 1.

In order to account for both regime ℳ≤ 1 and ℳ≥ 1 in a single equation, we propose here to modify the original <cit.> model 
by combining them into one single critical density

    s_crit = ln[ α_vir( 1 + 2 ℳ^4/1+ℳ^2)],

that we use in conjunction with Equation <ref>. <ref> shows the resulting star-formation efficiency per freefall time ϵ_ff as a function of the two parameters α_vir and ℳ. 
For low Mach numbers, the efficiency for star formation goes sharply from 0 if α_ vir>1 to 1 otherwise. For an isothermal gas, the corresponding recipe is just a simple density threshold, allowing star formation
only when the Jeans length is not resolved anymore by the mesh. This method has been used for a long time in galaxy formation simulations <cit.>. For larger Mach numbers, however, the efficiency iso-contours become wider.
Gas with α_vir>1 is allowed to collapse, owing to the large density fluctuations caused by supersonic turbulence, although the efficiency can be significantly lower than 100%. Typical conditions for star forming regions
in the Milky Way are α_vir≃ 3-10 and ℳ≃ 10-20, and correspond to ϵ_ ff≃ 0.01-0.02, in good agreement with observations <cit.>.
Traditional models of galaxy formation often adopt a similar value for their fixed ϵ_ff=0.01 <cit.>.



We now discuss some of the main caveats of the present approach. First, we have assumed a very simple spherical geometry for the collapsing transonic cores.
This is not true in reality, where filamentary structures are very often seen in star forming regions <cit.>. This effect could be parametrised through a multiplicative fudge factor 
in front of α_ vir in the collapse criteria. 

Second, in the multi-freefall formulation, we have assumed that within each transonic core the efficiency of star formation is 100%, and no gas is left after stars are born.
In our view, stellar feedback, as implemented in the next section, is responsible for terminating star formation. It is however possible to parametrise the efficiency of star formation per transonic cores using another 
multiplicative fudge factor in front of the integral in Equation <ref>. 

Third, we have also considered collapsing fluid elements at the sonic scales, but larger, supersonic fluid elements could also collapse
at a smaller density, further fragmenting into individual protostars. Such a model has been proposed for example by <cit.> and <cit.>.

Fourth, our lognormal model for the density fluctuations is only strictly speaking applicable to isothermal supersonic turbulence. At the relatively large scales considered here with Δ x ≃ 100 pc,
the ISM is far from strictly isothermal. This can affect the density distribution in a non-trivial way <cit.>. Having these caveats in mind, we nevertheless believe that any modification of the 
proposed model will not affect our results strongly, and certainly not at a qualitative level.



 §.§ Subgrid model for stellar feedback


Stellar evolution is probably the most extreme subgrid aspect in galaxy formation. 
In principle, we need to follow the evolution of the internal structure of the stars, 
describe their evolution throughout their main sequence, until they die and explode,
at least for the most massive ones. Although it is reasonable to assume that stellar interiors are fully decoupled
from galactic scales, this is not true for the final explosive stages of stellar evolution. 



Since our stellar particles have a typical mass of ∼ 10^5 M_⊙, we have to account for individual supernovae using a subgrid model.
For this, we assume that supernovae explosions within this single stellar population are uniformly distributed between t_start = 3 and t_end=20 Myrs. 
At each simulation time step, we compute the supernova rate as

    Ṅ= η_SN M_ini/ M_SN (t_end - t_start),

where η_SN=0.2 is the mass fraction of the population that explodes in supernovae, M_ini is the initial mass of the stellar particle and M_SN=10 M_⊙ is the typical mass of a single progenitor.
Multiplying by the time step Δ t, we obtain the expected (average) number of supernovae ⟨ N ⟩ = ṄΔ t. The actual number of exploding supernova for a given stellar particle during the time step is drawn from a 
Poisson distribution <cit.>. This allows us to have supernovae explosions as discrete events 
realistically distributed in time, independent on the adopted mass resolution.

The other challenge of modelling supernovae explosions comes from the requirement of resolving the energy-conserving Sedov phase of the remnant.
This phase is absolutely crucial as it is responsible for the momentum build-up that will later efficiently accelerate the surrounding gas.
The scale that marks the transition from energy-conserving to momentum-conserving is the cooling radius R_cool.

If the cooling radius is unresolved by the grid, which usually is the case for high gas densities, the thermal energy will be spuriously radiated away, without having the chance to create enough momentum.
Therefore, to properly simulate the effects of supernova on the gas, we inject momentum additionally to the thermal injection if the cooling radius is not resolved. This technique,
usually referred to as kinetic feedback or momentum feedback, has progressively emerged as one of the best subgrid models for supernovae explosions <cit.>.

On the other hand, if the cooling radius is resolved by the grid, which typically happens for low density gas, the supernova energy is directly injected in the form of thermal energy. 
The energy-conserving blast wave that will follow will then properly accelerate the gas around it and deliver the correct amount of momentum.

The terminal momentum of the momentum-conserving blast wave and the corresponding cooling radius can be computed analytically for homogeneous background densities and simple cooling functions.
For more realistic cases, including explosions within a supersonic turbulent medium, we can use high-resolution simulations of individual explosions, such as in <cit.>. 
Using these numerical models, we use for the cooling radius

    R_cool = 3.0 pc(Z/Z_⊙)^-0.082( n_H/100 cm^-3)^-0.42,

where n_H is the gas density of the cell where the star is exploding.
If R_cool is unresolved by at least one grid cell, we inject for each individual supernova the terminal momentum P_SN

    P_SN = 1.42 × 10^5 km s^-1M_⊙(Z/Z_⊙)^-0.137(n_H/100 cm^-3)^-0.16

Note that in the previous two equations, we use a metallicity floor at 0.01Z_⊙ to model the effect of primordial cooling for a pristine gas. 
The total amount of injected scalar momentum from each star particle into its surrounding cell is

    P = P_SN N_SNmin( 1,(Δ x_ min/R_cool)^3/2)

where we use the Sedov solution when Δ x_ min < R_ cool < 4 Δ x_ min. 
Here N_SN is the number of supernovae for the current time step and for each star particle, while Δ x_ min is the minimum AMR cell size.
<ref> shows the dependence of the injected momentum on the gas density for one supernova and for an adopted resolution of Δ x_min= 150 pc.
From <cit.>, we designed two models shown in <ref>. The first one is valid for cells with a low Mach number, while the second one applies for high turbulent Mach number.
In this paper, we use the low Mach number solution, but ideally, one should use the Mach number as an additional parameter. This model is unfortunately unavailable at the present time.

Once we know the injected scalar momentum per star particle, we need to deposit isotropically this momentum into the grid.
For this, we use novel numerical techniques, inspired by the work of <cit.> and <cit.>. 

First, we deposit the individual particle momentum onto the grid using the cloud-in-cell interpolation technique. We obtain the grid scalar momentum density p = P/Δ x^3 and convert it into a momentum flux density using 

    Q = p Δ x/Δ t = P/Δ x^2 /Δ t

Note that Q is analogous to a pressure. This momentum flux is split equally among the 6 cell faces, defining a supernovae dynamical pressure noted P_⋆=Q/6. This pressure only accounts for the momentum flux due to the
supernovae explosions. It is directly added to the thermal pressure in the Riemann solver. This direct injection of momentum through the Riemann solver was first explored by <cit.>, using non-thermal energy and pressure variables
that are available in the RAMSES code. This method deposits the right amount of momentum to the gas, but also affects the thermal energy with undesirable effects (spurious heating or cooling of the gas). 

We follow here a different approach by adding only this new pressure term in the momentum equation 

    t( ρ v_i )+ x_j( ρ v_i v_j ) + x_i( P+P_⋆) = -ρϕx_i

and remove its corresponding p dV work from the energy equation, 
so that the internal energy equation remains unaffected by this momentum injection

    Et+ x_j( E+P+P_⋆) v_j = P_⋆v_jx_j - ρϕx_j v_j

This approach, a slight modification of the original method from <cit.>, will deliver the proper momentum flux to neighbouring cells, together with consistent mass and energy fluxes.
We also modify the time step stability condition, increasing the wave speed in the Courant condition to |v| + c_s + p/ρ, 
where p is here the scalar momentum density defined above.

Our new method is more efficient than the traditional recipe for kinetic feedback in the RAMSES code, as implemented in <cit.> and <cit.>.
These models inject the right solution at the relevant time of the deposition. In the limit of very poor resolution, they effectively inject the terminal momentum along with thermal energy.
However, in these earlier attempts, the momentum is deposited on the grid using an explicit spherical velocity profile spread over a sphere of 4 cells in radius. This is usually done every coarse step and is quite costly. This is incompatible with our objective of resolving discrete supernovae with possibly a very high frequency of explosions.

On the other hand, depositing the momentum in only the 6 direct neighbours of the supernovae triggers spurious grid-aligned effects, especially if multiple supernovae explode in the same cell <cit.>. 
To avoid this, we choose the location of each individual explosion at random among the 8 cloud-in-cell neighbours from the star particle, suppressing visible grid alignment effects.

Although momentum feedback has proven quite successful in delivering the right amount of kinetic energy to the surrounding gas, a significant improvement over previous delayed cooling recipe, it still suffers 
from the caveat of not resolving properly the hot diffuse phase filling the cavity bounded by the momentum-conserving dense gas shell. This could affect the wind properties. 
<cit.> have shown recently in the context of a dwarf galaxy that the main properties of the outflow are indeed captured by the momentum injection scheme,
but properly describing the multiphase temperature structure of the wind requires to resolve the supernovae cooling radius.

Additionally, we model HII regions around young stars using a simple recipe for which we maintain the gas temperature (only in the cell where the star sits) at 10^4 K for 20 Myr, until the last supernovae explodes.
Thermal feedback from each supernova is modelled by injecting E_SN = 10^51erg in the parent cell. Furthermore, we assume that each supernova ejects 1 M_⊙ of metal into the ISM.



§ NUMERICAL EXPERIMENTS


We now report on the numerical experiment we have performed to test these new implementations of the key subgrid models
for galaxy formation. For this, we model the cosmological evolution of a single dark matter halo, whose mass is comparable to the 
Milky Way mass, but whose accretion history is much more violent. In particular, we choose a halo with one major merger at relatively low redshift
(expansion factor a>0.5). This will give us the opportunity to study how star formation and feedback behave in extreme cases. As it turned out, this particular halo
features a strong starburst at a=0.65, followed by an extremely quiescent state, that we identify as a typical example of the population of 
early-type galaxies with long gas depletion time scales <cit.>. It is worth stressing 
that in our model we don't use AGN feedback, so that any quenching mechanism here must be a consequence of our adopted stellar feedback and/or star formation models.



 §.§ Simulation setup


We base our analysis on a cosmological zoom-in simulation that we performed with the AMR code RAMSES <cit.>.
We first ran a reference pure N-body simulation with 512^3 particles in a periodic box of size 25 h^-1Mpc. 
First, we have selected all halos with masses in the range of 5 × 10^11≤ M_vir/M_⊙≤ 1.5 × 10^12, where the virial mass is defined as the mass contained in a sphere of radius R_vir that encompasses an over-density of 200 times the critical density.
Secondly, we applied an isolation criterion at z=0 where halos are excluded if there is another halo within 5 R_vir with a mass of 10% of the target halo.
Finally, we looked at the mass-accretion histories of the available halos and selected one halo which experienced its last merger at a=0.65. This halo has a mass of M_vir=6.5 × 10^11M_⊙ at z=0 and was used for our analysis.

We then generated new initial conditions around this halo with the MUSIC code <cit.>, with an initial hierarchy of concentric grids from ℓ_ min=7, corresponding to a coarse grid resolution of 128^3 covering the entire periodic box, to ℓ_ max,ini=10, corresponding to an effective initial resolution of 1024^3. 
This gives us a dark matter particle mass of m_ dm=1.7 × 10^6M_⊙ and a baryonic initial mass resolution of m_ dm=2.6 × 10^5M_⊙.

We then performed 
several simulations including gas and galaxy formation physics. The maximum resolution was set to ℓ_ max=18 at z=0, while refinement levels were 
progressively released to enforce a constant physical resolution of Δ x_ min=100 h^-1pc. The adopted refinement criterion is the traditional quasi-Lagrangian 
approach, namely cells are individually refined when more than 8 dark matter particles are present or when the baryonic mass (gas and stars) exceed 8× m_ bar. 
Only the Lagrangian volume corresponding to twice the final virial radius of the halo was refined, the rest of the box being kept at a fixed, coarser resolution to provide
the proper tidal field. Our zoom simulation corresponds to a halo slightly smaller than the Milky Way, but it features two successive major mergers, one at a ≃ 0.35
with a stellar mass ratio of 1:1 and a stellar mass of each galaxy of M_* ≃ 1.5 × 10^9M_⊙, and one later at a ≃ 0.65 with a stellar mas ratio of 1:1 too, and a larger galaxy stellar mass of M_* ≃ 5 × 10^9M_⊙. This scenario turned out to be an ideal one to study extreme regime of 
star formation and analyse the impact of our star formation and feedback recipe. 

Gas cooling and heating is implemented using equilibrium chemistry for Hydrogen and Helium <cit.>, together with the metallicity dependent cooling function 
for metals from <cit.>. Additionally, self-shielded UV-heating of the gas is taken into account, where the self-shielding density is assumed to be 
n_H=0.01 H/cc <cit.>. We adopt the model of <cit.> for the UV background.

We ran three different simulations with different models for the galaxy formation physics. 


  * The first simulation was run without feedback but using our new multi-freefall star formation model. It is labelled “no feedback” in the figures. 

  * The second simulation was run with stellar feedback but adopting the old school constant efficiency model with ϵ_ ff=0.01
and a star formation density threshold of ρ_*=0.1H/cc. It is labelled “constant efficiency” in the figures.

  * The third simulation was run with both our new recipe, namely stellar feedback and our varying ϵ_ ff model. It is labelled “varying efficiency” in the figures.

Note that we could have run the second simulation with a higher density threshold, to mimic a criterion α_ vir<1 for some fixed temperature and Mach number.
This would give results intermediate between our two cases. Note that if one adopts a high density threshold, ρ_*=10H/cc for example, it is impossible to form star in lower density gas, so that quenching occurs automatically in this case, in a rather ad-hoc way. 
Our multi-freefall approach, on the other hand, allows for star formation even in gas with α_ vir>1 if the subgrid turbulence is strong enough.

We finally list our more technical RAMSES settings for these 3 simulations. We used the HLLC Riemann solver, modified to account for the supernovae dynamical pressure.
The star particle was set equal to the baryonic mass resolution m_*=m_ bar. We adopted the MinMod slope limiter, an important choice to ensure
the stability of the new feedback scheme. 
The adiabatic exponent of the gas was set to γ=5/3. 
We didn't use any polytropic pressure floor, as our star formation model, by construction, very quickly removes gas for which the Jeans length 
is not resolved. The gas metallicity was initialised to Z_ ini=10^-3Z_⊙ to account for early population III stars enrichment.





 §.§ Effect of the feedback model


We plot in <ref> the evolution of the stellar mass of the central galaxy in our simulated galaxy for the 3 different subgrid models. 
The central galaxy is defined here as everything within 0.1 × R_ vir from the centre. As a validation test, we also show as a shaded area 
the expected stellar mass of the central galaxy from abundance matching <cit.>, using the simulated halo mass as input at each epoch. 
The simulation without feedback completely overestimates the stellar mass (by a factor of 5), while our two models with feedback reproduce roughly the expected evolution.
We see a clear trend for the constant efficiency model to produce almost twice as many stars than the multi-freefall model. 
The two major mergers in the mass accretion history of this particular halo can be seen as two jumps in the stellar mass of the central galaxy at a ≃ 0.35 and a ≃ 0.65.

In <ref>, we show the evolution of the gas fraction within 0.1 × R_ vir, defined by f_gas = M_gas/(M_gas + M_⋆). 
Ironically, the lowest gas fraction is obtained for the no feedback run, as most of the gas is consumed and turned into stars. The constant efficiency model, on the other hand,
has the largest gas fraction, with f_ gas≃ 0.4 at late time, while the varying efficiency model shows larger fluctuations, with peaks of gas fraction
during star bursts, and a lower value of f_ gas≃ 0.2 at late time, in better agreement with observation (see Discussion section). 
Our interpretation is that stellar feedback in the varying efficiency model is stronger, owing to the higher efficiency in the star bursts mode (see next Section).
This is corroborated by <ref> that shows the baryon fraction with R_ vir. The no feedback case manages to accumulate more baryons 
in the virial radius than the universal fraction, while feedback maintains the baryon fraction below the universal value in all cases. The varying efficiency case,
however, only retained 40% of the baryons in the halo, meaning that 60% of the baryons are lurking outside the virial radius, probably never coming back.



 §.§ Effect of the star formation model


Although the stellar feedback subgrid model seems to play the main role in regulating the stellar mass of the galaxy, 
the star formation subgrid model seems to have a non-negligible impact on our results.
The top panel of <ref> shows the star formation rate (SFR) as a function of the scale factor for our 2 models with feedback.
It is defined as the instantaneous SFR within 0.1 × R_ vir at each time.
We ignore the model without feedback in what follows because we believe it is not realistic.
We clearly see the two major mergers as two strong starbursts with the SFR peaking at 10 M_⊙  yr^-1. 
We used a bin size of Δ a=0.01 corresponding roughly to 100 Myr at low redshifts. 

In the varying efficiency case, the SFR during the starbursts is slightly larger than for the constant efficiency model. But this is after the starbursts that 
the difference between the 2 models is striking. The SFR for the varying efficiency model drops significantly down to 0.1 M_⊙  yr^-1, while it remains
quite high, around 1 M_⊙  yr^-1 for the constant efficiency case.

This effect can be partly explained by the higher gas and baryon fraction in the constant efficiency case.
Indeed, we see in the top panels of Figure <ref> the morphological evolution of the corresponding central galaxy using true colour images with dust absorption
that allows us to render both stars and gas. The constant efficiency case, owing to its higher gas and baryon fraction, leads to the formation of a large gaseous
and star forming disk. This is a well-known outcome of gas-rich mergers <cit.>.
The varying efficiency case, on the other hand, leads to the formation of a spheroidal galaxy with a tiny nuclear gas disk (see bottom panels of Fig. <ref>). 
In this case, the gas fraction is much lower, leading to the formation of dispersion dominated systems as explained by the gas-poor, dry merger scenario 
<cit.>.

The origin of this difference can also be found in the star formation efficiency itself. 
In the lower panel of <ref>, we plot the mean efficiency of the star forming gas
within 0.1 × R_ vir, defined as 

    < ϵ_ ff> = ∫ϵ_ ffρ/t_ ff dV / ∫ρ/t_ ff dV

The constant efficiency with ϵ_ ff=0.01 is also shown for comparison. We clearly see that during the starbursts the average
efficiency peaks at up to 3%, while in the post-starburst quiescent phases, it drops down to 0.1%. The high efficiency we obtain during 
the mergers explains how the galaxy, although relatively gas poor compared to the other model, manages to have a SFR as high as 10 M_⊙  yr^-1,
like in the constant efficiency case. This helps maintaining the gas and baryon fractions low. 

Note that this high efficiency is significantly higher than the critical value
of 1% proposed by <cit.> that divides the feedback-dominated star formation regime from the efficiency-dominated star formation regime.
Our constant efficiency model with ϵ_ ff=0.01 sits just at the limit, so the resulting galaxy never reaches the feedback-dominated regime.
Our varying efficiency model, on the other hand, clearly explores deep into these two regimes: 
1- strong feedback-dominated starbursts, for which the exact value of ϵ_ ff
probably plays a minor role, as long as it is larger than 1% and 2- long, extended, post-starbursts quiescent phases, for which the galaxy is quenched and the efficiency drops down to 0.1%.

In order to understand the origin of these wide variations in efficiency, we show in Figure <ref> a mass-weighted histogram of the star forming gas
within 0.1 × R_ vir in the α_ vir-ℳ plane. In the starburst case, at a=0.65, most of the star forming gas lies close
to the 10% efficiency line. This is a combination of two factors: First, the merger event triggers the fast migration of gas to much higher density and
smaller α_ vir <cit.>. The mean gas density during the starburst is roughly n_ H≃ 20 H/cc. 
In the same time, the merger favours a higher gas turbulent velocity dispersion, with a mass-weighted average value as high as σ_ 1D≃ 70 km/s. 
Both effects work in tandem to maintain the efficiency around 10% during the merger. This leads to extremely small depletion times, around 100 Myr,
but the gas is also quickly collapsing and this high efficiency state can be maintained during the duration of the merger.

In the quiescent phase, we see in Figure <ref> that most of the gas sits close to the 0.1% efficiency line. We are in the exact opposite situation than
the starburst: residual gas accretion maintains a non-negligible, albeit lower, level of turbulence, around 20 km/s, while in the same time the mean gas density dropped to
a rather low value with n_ H≃ 0.5 H/cc. This leads to extremely long depletion times, around 10 Gyr. This residual star formation occurs mostly
in a small nuclear disk a few kpc in size. 

Milky Way-like galaxies correspond to the intermediate regime between the starburst and the quenched galaxy. We have explored this regime
in a companion paper using a different halo (Kretschmer et al. in prep) and indeed, the mean efficiency settles naturally to 1% in the disk. Our varying efficiency model 
automatically adjusts its value to the conditions within the galaxy, each main galaxy type (starburst, disk, spheroid) giving rise to a different regime of star formation. 
This is due to a different structure of the ISM, with more or less gas being able to reach the limit of collapsing transonic cores.










 §.§ Effect of the subgrid turbulence model

The last important ingredient in our subgrid galaxy formation model is the SGS turbulence model.
As explained in the previous section, we use a source term for the turbulent kinetic energy based on the local shear tensor
and a sink term based on dissipation of kinetic energy over one turbulent crossing time. 
In order to test the robustness of our star formation model to the adopted model for turbulence, we decided to change the shear-based source term 
classically used in the SGS model by a supernova-based source term, in the spirit of the work of <cit.>. 
Note that in the classical SGS model, supernovae explosions will also inject turbulence indirectly though momentum deposition followed by the corresponding shear-induced source term. We nevertheless turn the shear source term off and replace it by direct turbulent kinetic energy 
being deposited in the cell where the supernova explodes, assuming 10% of the energy in turbulent form. 

The top panel of <ref> shows the evolution of the mass-weighted average turbulence resulting from these two models. 
The direct injection of turbulent energy from SN and the removal of the shear-based source term produces a level of turbulence 
roughly one order of magnitude smaller than the SGS model. 
Interestingly, this has a weak effect on the resulting star formation efficiency, as shown in the bottom panel of <ref>. 
It is apparent that the difference between the two models is quite small, even though the actual level of turbulence is very different. 
This is because the varying efficiency model adapts to the local conditions in the galaxy and will transfer gas to smaller α_ vir (or higher density) 
compared to the SGS model, maintaining a similar star formation efficiency.

There is a small difference in scalefactor for the starburst at a≃0.65 which originates from slightly different trajectories caused by numerical effects <cit.>.
Furthermore, the efficiencies for the SN induced turbulence model in the quiescent phases is larger compared to the efficiencies in the SGS run.
Consider for example the left panel of <ref>. In the SN turbulent case, many cells will have much smaller α_ vir and ℳ, bringing the individual gas cells into a regime where star formation is dominated by the collapse criterion through α_ vir. Therefore, during a quite phase, in the SN injected turbulence case the efficiencies are larger because the small values for the turbulence will bring cells into a regime where star formation is controlled by the virial collapse.









§ DISCUSSION


In this paper, we have introduced a new set of subgrid models for star formation and feedback, uniquely combined in this new version of RAMSES.
The feedback model, very similar qualitatively to earlier implementations of terminal momentum injection by supernovae, but quite different in the details, provides
the first order mechanism for regulating the stellar and baryonic content in our simulated galaxies. The multi-freefall star formation model, based on a varying 
efficiency based on a subgrid turbulence model, allows us to form star without invoking arbitrary values for both the efficiency and the density threshold. 
We have shown that in a starburst regime, the efficiency can be significantly higher than 1%, at odd with traditional models of galaxy formation, while in a
quenched regime, the efficiency can be as low as 0.1%, owing to the very different nature of the turbulent ISM in post-starburst, early-type galaxies.

Our star formation model, based on an SGS model for unresolved turbulence, and the subgrid model for supernovae momentum feedback 
are naturally designed in such a way that they are (relatively) resolution independent. When resolution is increased, the kinetic energy of turbulence
will decrease, following the scaling for Burger's (or Larson's) turbulence. Ultimately, when the sonic length is resolved, the SGS turbulence will become
subsonic and our star formation recipe with turn into a simple collapse criterion for thermally supported gas with 100% efficiency. 
At these scales, however, one usually adopts a sink particle formalism to describe star formation.

For supernovae feedback, increasing the resolution will allow us to resolve more the cooling radius, assuming the supernovae explode in a fixed density environment.
This won't be the case, unfortunately, as simulations with better and better resolution lead to denser and denser environments.
In the limit of sub-parsec resolution, we know from dedicated ISM studies <cit.> that a crucial ingredient is the inclusion of 
walk-away and run-away massive stars.
Overall, the galaxy formation simulations we present in this paper are far from resolving these scales, so that our subgrid models are crucial in providing a realistic
treatment of these unresolved phenomena. 

As explained in the previous section, our simulated halo allowed us to study two different regimes of star formation: a merger-induced starburst and a quenched early-type galaxy.
The SFR in these two extreme cases are widely different, with 10 M_⊙  yr^-1 for the starburst and 0.1 M_⊙  yr^-1 for the early-type. 
These variations can be explained partly by the different gas fraction, but also by the different typical density and Mach number in these galaxies. 
We now compare the properties of our simulated galaxy to observed ones using the global KS relation in Figure <ref>.
For this, we compute the total star formation rate in the galaxy by integrating the star formation rate in each gas cell within 0.1 × R_ vir
and the total gas mass within the same region. We then use the half-mass radius of the gas r_h to compute the global surface densities of both 
star formation and gas as <cit.>

    Σ̇_* = 0.5 Ṁ_*/π r^2_h    and   Σ_ gas = 0.5 M_ gas/π r^2_h.

In Figure <ref>, each symbol corresponds to the location of the galaxy in the global KS diagram at a different epoch. 
We also indicate 3 lines corresponding to a constant depletion time of 100 Myr, 1 Gyr and 10 Gyr. 
Although high redshift conditions lead in general to small depletion times, the galaxy is moving quite widely along the KS relation <cit.>.
The two diamond symbols show the most extreme conditions we have obtained, with the starburst at the rightmost tip of the KS relation,
and the quenched early-type at the leftmost end. 

To understand better the origin of these two extreme cases, we show in Figure <ref> the distribution of the gas mass
in a resolved KS diagram. For this, we compute the star formation and gas surface densities using square patches of 200 pc size
and compare it to the observed KS relation as well as constant depletion time models. In each case (starburst and early-type) we also show
as a single diamond symbol the global values. We clearly see that the gas and star formation surface densities is the starburst case are both very large,
with very small depletion times, around 1 Gyr in average but as small as 100 Myr in some regions. 
The early-type galaxy at the final epoch, on the other hand, has an average depletion time around 20 Gyr, with some regions reaching barely 1 Gyr.

In order to compare to molecular observations, we use the model of <cit.> to compute the molecular H_2 mass within 
each cell, assuming for the interstellar radiation field a global uniform value given by G=G_0 (SFR / M_⊙ yr^-1) at each epoch.
For the merger-induced starburst at a ≃ 0.65, we find M_ gas≃ 2 × 10^9 M_⊙ and M_H_2≃ 9 × 10^8  M_⊙,
while for the early-type quenched galaxy at a=1 we have M_ gas≃ 4 × 10^9  M_⊙ and M_ H_2≃ 8 × 10^8  M_⊙.
Note that our model for H_2 formation is based on exactly the same turbulence-based subgrid model for star formation, so that there is a built-in correlation
between star formation and molecular gas formation, but no direct causal relation.

In the early-type case, the long depletion time coincides with a low H_2 fraction, so that the molecular gas depletion time becomes shorter by a factor of ∼5, around 20 Gyr. 
This is in perfect agreement with the COLD GASS sample observed at IRAM <cit.>, with our early-type case corresponding to a low specific star formation rate sSFR≃ 10^-12  yr^-1 and long molecular depletion times t_ dep( H_2) ≃ 20 Gyr, and our starburst case corresponding to LIRGs or ULIRGs <cit.>. 

Interestingly, post-starburst galaxies show many similarities with our quenched mode of star formation. 
In several recent papers, galaxies showing signs of recent starburst activity are completely quenched,
although they contain enough gas to form stars at a rate much higher than observed <cit.>.
Residual turbulence in the ISM or high levels of radiation were proposed as possible explanations to quench star formation in such systems. 
These systems are probably ideal examples of galaxies with a very peculiar ISM structure, leading to an exceptionally low star formation efficiency.


At late times the baryon fraction in the halo stays constant. At the same time the gas fraction in the disk increases by ∼ 10% due to gas cooling from the halo onto the disk. The star formation rate however does not increase with the gas fraction but rather decreases. This is because the efficiency depends on the shear-induced turbulence in a non-linear way. Namely the level of turbulence in the disk from accretion or shear will prevent star formation.
In contrast, the star formation rate given by the constant efficiency model does increase with the gas fraction of that run.



Another interesting class of objects are circumnuclear disks at the centre of early-type galaxies. Our simulation with varying efficiency 
shows at the final epoch such a nuclear disk (see Fig. <ref>). Nuclear disks in early-type galaxies are compact, relatively dense in their centre 
with Σ_ gas≃ 100  M_⊙  pc^-2 but very inefficient in forming stars <cit.>,
when compared to normal galaxies with the same gas surface densities. The strong shear observed in these systems was proposed as 
an explanation for their surprisingly low star formation efficiency. One particularly striking example is NGC 4429 from the same ATLAS3D sample:
In a recent attempt to model it, Naab and co-workers (private communication) had to use a constant efficiency of ϵ_ ff≃ 0.2%
to reproduce various properties such as gas content and star formation rate. This is very close to the efficiency we obtained using the 
varying efficiency model in our final early-type galaxy. 



§ SUMMARY AND CONCLUSION

In this paper, we have introduced a new implementation of recently proposed subgrid models for star formation and feedback in the RAMSES code. 
The stellar feedback model is based on individual, discrete Type II supernova events. If the radius marking the transition from the energy-conserving phase 
to the momentum-conserving phase is unresolved by our computational cell, we inject directly momentum in the neighbourhood of the exploding star. 
This is done by modifying the Riemann solver in our Godunov scheme, allowing us to deliver the right amount of momentum, together with consistent mass and energy fluxes.

The multi-freefall star formation model is based on a subgrid turbulence model and provides a varying star formation efficiency that depends on the local conditions in each 
computational cell, through the cell's virial parameter and its turbulent Mach number. 
The star formation efficiency can be very large for large density fluctuations caused by supersonic turbulence, as expected, for example, during a merger. 
On the other hand, the star formation efficiency can be very low for subsonic and marginally stable disks, as expected during more quiescent phases or in early-type galaxies. 
This new versatile model allows us to form star without invoking arbitrary values for both the efficiency and the density threshold.

We apply this new implementation of subgrid galaxy formation physics to a prototype cosmological simulation of a massive halo that features a major merger, 
and study specifically the interesting case of the formation of an early-type galaxy.
We find that the feedback model provides the first order mechanism for regulating the stellar and baryonic content in our simulated galaxies.
Together with the multi-freefall star formation model, feedback is strong during starbursts. At redshift zero, only 40% of baryons are retained in the halo and the gas fraction is ∼ 20%. 
The resulting stellar masses are in good agreement with those expected from abundance matching results.

In a starburst regime, the efficiencies can be significantly higher than 1%, at odd with traditional models of galaxy formation, while in a quenched regime, the efficiencies can be as low as 0.1%, 
owing to the very different nature of the turbulent ISM in post-starburst, early-type galaxies.
Indeed, the merger driven event pushes gas to large densities and large turbulent velocity dispersions, which causes larger ℳ and small α_vir. 
As a result, efficiencies can reach locally 10% and the sSFR is large. 
The small value for the molecular gas depletion time during the starburst is in perfect agreement with observations. 
Additionally, at late times, when the galaxy is quiescent, the sSFR is small and the molecular gas depletion time is long, also in perfect agreement with observations.

In summary, at high redshifts, very efficient star formation together with strong feedback regulates the baryonic content. 
At low redshift, local star formation becomes very inefficient explaining the properties of quenched galaxies.



§ ACKNOWLEDGEMENTS

The authors thank the referee for their constructive comments that improved the quality of the paper.
We acknowledge stimulating discussions with Robert Feldmann, Lucio Mayer, Robbert Verbeke, Renyue Cen and Sandro Tacchella. 
This work was supported by the Swiss National Supercomputing Center (CSCS) project s890 - “Predictive models for galaxy formation” 
and the Swiss National Science Foundation (SNSF) project 172535 - “Multi-scale multi-physics models of galaxy formation”.
The simulations in this work were performed on Piz Daint at the Swiss Supercomputing Center (CSCS) in Lugano,
and the analysis was performed with equipment maintained by the Service and Support for Science IT, University of Zurich.
We also made use of the pynbody package <cit.>.

mnras






