---
abstract: 'Galaxy evolution depends strongly on the physics of the interstellar medium (ISM).  Motivated by the need to incorporate the properties of the ISM in cosmological simulations we construct a simple method to include the contribution of non-thermal components in the calculation of pressure of interstellar gas.  In our method we treat three non-thermal components - turbulence, magnetic fields and cosmic rays - and effectively parametrize their amplitude.  We assume that the three components settle into a quasi-steady-state that is governed by the star formation rate, and calibrate their magnitude and density dependence by the observed Radio-FIR correlation, relating synchrotron radiation to star formation rates of galaxies.  We implement our model in single cell numerical simulation of a parcel of gas with constant pressure boundary conditions and demonstrate its effect and potential.  Then, the non-thermal pressure model is incorporated into [[*RAMSES*]{} ]{}and hydrodynamic simulations of isolated galaxies with and without the non-thermal pressure model are presented and studied.  Specifically, we demonstrate that the inclusion of realistic non-thermal pressure reduces the star formation rate by an order of magnitude and increases the gas depletion time by as much.  We conclude that the non-thermal pressure can prolong the star formation epoch and achieve consistency with observations without invoking artificially strong stellar feedback.' 
author:
- |
    Yuval Birnboim$^\text{1}$, Shmuel Balberg$^\text{1}$, & Romain Teyssier$^2$\
    $^1$Racah Institute of Physics, The Hebrew University, Jerusalem 91904 Israel\
    $^2$Institute for Computational Science, University of Zürich, Winterthurerstrasse 190, CH-8057 Zürich, Switzerland\
date: 'Accepted —.  Received —; in original —'
title: 'Galaxy Evolution: Modeling the Role of Non-thermal Pressure in the Interstellar medium'
---

\[firstpage\]

galaxies: evolution — hydrodynamics — cosmic rays — ISM: general — ISM: magnetic fields —

Introduction
============

Both star formation and stellar feedback play crucial roles in galaxy evolution.  Star formation leads to stellar feedback, which in turn is assumed to regulate the star formation rate and prevent galaxies from turning all their gas into stars over less than a ${\,{\rm Gyr}}$; moderate star formation rates are implied from low and high redshift observations.  In contrast, pure hydrodynamic simulations and semi-analytic models of galaxy formation tend to predict high gas densities within galaxies.  These densities cause the gas to cool very efficiently and supersede the density threshold required for star formation [@schmidt59; @kennicutt98].  For purely hydrodynamic simulations, unrealistically strong stellar feedback is often necessary to regulate the star formation rate in galaxies [@Scannapieco2012]. 

Several approaches have been attempted to regulate simulations of star formation during galaxy evolution.  Within the framework of pure hydrodynamics, the most basic feedback mechanism is usually thermal feedback, or injection of some fraction of the supernova energy into the gas.  Since this energy is injected into dense, cold gas it cools efficiently and typically has small overall effect on galaxy evolution [@Scannapieco2012].  Momentum feedback is added by explicitly injecting momentum to the material that surrounds star forming regions [@Navarro1993; @springel03_sf; @Oppenheimer2006; @Dubois2008].  While the efficiency of such models is slightly better than thermal feedback [@Scannapieco2012], the ejected gas is almost always highly supersonic and kinteic energy is converted to thermal radiation very efficiently through shocks, unless the feedback is injected effectively over very large volumes [@Oppenheimer2006].  Attempts have also been made to inject the energy into a warm component of a two-phase gas, effectively delaying the cooling until energy is transferred from the diffuse gas to a denser gas which immediately cools [@springel03_sf; @Governato2007].  An additional technique sometime used has been to increase the efficiency of the feedback by releasing the feedback energy in bursts rather than spread out [@Crain2009].  However, in all these recipes, the cooling of a parcel of gas near the plane of the disc still occurs with isobaric boundary conditions set by the weight of the atmosphere on top of it.  Once the gas cools, it contracts over a crossing time to regain its pressure.  Since the cooling rate scales as $\rho^2$ (where $\rho$ is the gas density) this leads to runway cooling.  In order to increase the efficiency of supernova feedback, the cooling of injected energy is artificially delayed [@springel03_sf], and momentum feedback efficiency in enhanced by preventing it from interacting with its immediate, dense environment [@Oppenheimer2006]. 

In this paper we revisit the conjecture that non-thermal processes contribute to the total pressure.  With this additional pressure, gas can reach hydrostatic equilibrium with a considerably lower gas density that naturally predicts lower star formation rate and bypasses the need for unrealistic supernovae feedback.  The non-thermal pressure does not depend on the temperature of the gas, and the gas cools isochorically rather than isobarically, further stabilizing the gas.  We develop a simple, easy to use, parametric model which allows to study (analytically and in simulations) the effect of such non-thermal components on galaxy formation. 

The enhanced star formation problem is closely related to the general complication of modeling the interstellar medium (ISM) gas.  While the use of a standard, purely thermodynamic equation of state of an ideal gas is justified outside of galaxies in the IGM, it becomes less appropriate to use in haloes (haloes of galaxy clusters exhibit non-negligible magnetic fields) and even more so in the ISM of galaxies.  This gas is highly multiphased, and consists of cold, warm and hot gas arranged within atomic and molecular clouds, filaments, and bubbles [@mckee77; @ferriere01].  Complicated chemistry and dynamics, as well as radiation fields at multiple wavelengths affect the behaviour and interrelation between the different phases.  Moreover, the dynamics of the gas are strongly affected by non-thermal components, namely turbulence, magnetic fields and cosmic rays (CR).  Stars, through their formation, evolution and destruction pump energy into the ISM by stirring turbulence, emitting high energy particles (cosmic rays) and releasing radiation that heats and drives the gas [@deJong1985; @bell03_LIR].  Gravitational energy also powers turbulence [@dekel09] and heats the gas [see however @Hopkins2013].  Turbulence, could, in principle be accurately followed by pure hydrodynamics.  However, modeling turbulence requires high resolution and realistic driving of the turbulence which is still an open question [@Schmidt2009].  The detailed modeling of these effects is the subject of intensive ongoing efforts [e.g. @Korpi1999; @MacLow2004; @Elmegreen2004; @Dib2006; @Robertson2008; @Koyama2009; @Hopkins2012; @Kim2013].  All these physical phenomena are determind by the relatively small ($\sim\;$pc) scale of observed giant molecular clouds (GMCs), large eddies of turbulence [@Schmidt2010], and tangled magnetic fields. 

It is prohibitively challenging to include all the aforementioned effects and small scales in cosmological simulations.  Effective equations of state for star forming gas are thus constructed, directly pressurizing the thermal component of the gas [@springel03_sf; @schaye08].  An equation of state for subgrid turbulence has been proposed in @maier09.  @Joung2009 proposed an effective EoS directly related to star formation rates, and @Braun2012 incorporated turbulent pressure as a sub-grid model of the various phases motivated by @mckee77.  @Scannapieco2010 showed that sub-grid models for supersonic turbulence can have a large affect on dwarf galaxies.  While these attempts artificially pressurize the gas (as we propose below) they do not prevent the gas from over-cooling isobarically and still require the unrealistically strong feedback described above.  Recent attempts [@Salem2014; @Booth2013; @Hanasz2013] to simulate the effects of cosmic ray pressure on star formation and winds in galaxies have shown that they are able to drive significant outflows and could be efficient in regulating the star formation of star forming galaxies.  @Booth2013 [@Salem2014] separately implemented a two fluid approximation for cosmic rays and gas for single-galaxy simulations, and propagated the CR as a diffusive component with constant diffusion coefficient.  @Hanasz2013 used an MHD code to simulate anisotropic diffusion that preferentially diffuses CRs along magnetic fields.  In all three implementations the CRs are found to drive winds in some cases.  These results are encouraging but the simulations are of single, ideal galaxies and the suggested methods cannot be easily extended to the evolution phases of the galaxies and cosmological initial conditions.  Their physics is too detailed and the required resolutions are too fine to be practical in full hydrodynamic cosmological simulations.  To date, no simulations exists that include all the physical processes which are important for pressurizing the gas, and no cosmological-scale simulation will likely be able to simulate all this physics in the forthcoming future. 

In view of all these complications, here we follow a simpler, alternative avenue, and construct an effective EoS that mimics some of the main observed characteristics of the non-thermal pressure components.  First, we take advantage of the fact that a scale separation roughly exists between the parsec-scale phenomena discussed in the previous paragraph and that of typical observed scales for vertical scale heights of discs which range between $\approx 100 {\,{\rm pc}}$ for the Milky Way [@ferriere01] and quiescent star forming galaxies to $1{\,{\rm kpc}}$ for starbursting high redshift galaxies [@tacconi06].  This implies that in the context of galaxy formation simulations, it should be sufficient to resolve the ISM on scales much larger than $\sim$1 parsec in order to reproduce galactic discs with realistic characteristics.  Correspondingly we use a coarse grained, effective modeling of the ISM.  Our model bridges the gap between the parsec and kilo-parsec scales.  Second, we complement this scale separation with an effective equation of state that is straightforward to use.  In principle, one can attempt to apply a rigorous treatment of the physical processes that constitute the non-thermal physics as sub-grid models.  However, even if realistic such modules were constructed, there remains the problem of stipulating physically-consistant initial conditions: how to seed magnetic fields, how and when cosmic rays are generated and accelerated, and what drives turbulence and precisely on which scales.  In addition, one has to relate the initial conditions to star formation which most likely drives these effects.  Hence, we opt for a simple, easy to use, pressure-density relation for the non-thermal EoS which we develop below. 

A key feature of our implementation is the calibration of the EoS by the observed relation between the FIR radiation - a star formation indicator - and the radio radiation, which constrains the joint energy content of cosmic rays and magnetic fields.  This novel approach provides a quantitative relation between the gas density and the magnitude of the non-thermal components.  Additional physical assumptions which are required to complete the specific parametrization of the effective EoS are then limited to factors of order unity, rather than being arbitrary. 

The structure of this paper is as follows.  In [§\[sec:eos\]]{} we describe possible modifications for the equation of state of the gas that manifest some important aspects of the non-thermal pressure components.  In [§\[sec:toymodel\]]{} we demonstrate the effectivenes of such modified EoS and and the importance of the non-thermal conmponents in general, using simple calculations of a point (single-cell) model, focusing on the regaulation of the overall star formation rate. [§ \[sec:ramses\]]{} describes the incorporation of the model into the hydrodynamic code [[*RAMSES*]{} ]{}[@Teyssier2002] and the setup and results of isolated galaxy simulations with and without our model.  In [§\[sec:summary\]]{} we summarize and discuss our results. 

Equation of State {#sec:eos}
=================

To incorporate non-thermal pressure we need some typical scale for its density dependence and amplitude.  We will use the observed FIR-radio relation for a typical value of magnetic field and.  The model we describe in this paper assumes that (i) the non-thermal components are in equilibrium between themselves in the sense that energy can move quickly between them, and (ii) that this equilibrium does not depend on the magnitude of the energy.  The first assumption, of strong coupling between the components, is justified because the timescales for interactions between the components are eddy turnaround time for the turbulence, alfvénic crossing time for magnetic fields and diffusion times for cosmic rays, all on scales much smaller than galactic or cosmologic scales.  An eddy turnaround time for a $1{\,{\rm pc}}$ eddy rotating at a typical ISM speed of $5km/sec$ would take less than $1{\,{\rm Myr}}$ - much smaller than cosmological evolution timescales.  The alfvénic crossing time is even shorter: the velocity for $B=5\mu G$ and $n=10^{n-3}\text{cm}^{-3}$ is about $10km/sec.$ The diffusion coefficients for the cosmic rays indicate an even smaller timescale for equilibrium across a $1{\,{\rm pc}}.$ With a typical diffusion coefficient of $D\sim 10^{27}\text{cm}^2 sec^{-1}$ the diffusion over this scale-length will take $t\sim L^2/D\sim 300{\,{\rm yr}}$.

We stress that the second assumption is not necessary for the effect of the non-thermal pressure to be important, and that we use it below for the sake of simplicity.  There are, nonetheless, strong qualitative arguments which motivate such strong coupling.  Compelling physical arguments can be made in favor of the interrelation between turbulence energy and tangled magnetic field[^1] energy and cosmic rays.  Turbulence and magnetic fields are naturally related since turbulent flow can increase the energy in magnetic fields (by elongating and wraping the field lines), while large magnetic fields tend to rearrange and freeze the material in order to decrease the length of the flux tubes [see @Federrath2011 for a discussion including the dependence on magnetic field geometry and turbulence velocity] .  Either way, energy is naturally converted from one component to the other.  We assume that this qualitative argument holds, even though the dominant coupling process and the precise energy distribution between the components might vary somewhat [see discussion and references in @Lacki2013 who also assumes equal energies in turbulence and magnetic fields for starbursting galaxies].  Equipartition between magnetic and cosmic rays should also be natural [@longair94; @lisenfeld96; @bell03_LIR see however @stepanov12 [@lacki10a]].  Cosmic rays are relativistic electrons and protons (most likely accelerated during supernavae explosions), and travel along magnetic fields which confines them to the galaxy with some effective diffusion regulated by the fields.  Strong magnetic fields increase cosmic ray energy loses through synchotron radiation, and reduce the diffusion rate thus increasing the steady-state density of the cosmic rays in the galaxy. 

The simplest manifestation of such an approach is to assume equipartition between the three non-thermal components, so that the total energy density is three times that of each separate component.  Equipartition between magnetic fields and cosmic rays corresponds to a minimum of the total energy of cosmic rays and magnetic fields for a given measured synchrotron radio emission [@Lacki2013a].  In cases where the two are not in equipartition the combined non-thermal pressure due to the two components will be larger.  Finally, the cycle is completed by instabilities in the CR flows along magnetic fields giving rise to small scale turbulence [@Kulsrud1969].  It is worth mentioning that an equipartition between these three components are also observed in intra-cluster medium (ICM) where they each contribute about $5\%$ of the total pressure, while the bulk of the pressure comes from the thermal component of X-ray emitting gas.  We emphasize that the “equipartition anzatz” is basically just a scaling parameter for the effective non-thermal EoS; it is not an essential component of our model and different scalings can be used (see below). 

Once the distribution of energy among the non-thermal components is determined, we can naturally continue to develop an equation of state for them.  While a “proper” thermal equation of state relates all the thermodynamic variables to two independent variables (for example, the density and internal energy), a non-thermal component is, by nature, a single parameter function.  In the case of magnetic fields, for example, if a volume of space occupying magnetic fields is compressed, no heat is generated, and the process can be reversed [^2].  Hence, the magnetic pressure and energy are functions of only one variable (the field magnitude, $B$) which can in turn be related to only one thermodynamic variable.  This is completely analogous to the equation of state of cold matter, which is commonly used in the analysis of compact objects. 

It is important that the concept of entropy does exist in multi-component non-thermal system, through the requirement for equilibrium.  We note that the system is not closed, since energy is constantly pumped in by star formation and leaks out of the system by cosmic ray diffusion, electromagnetic radiation, reconnections, and dissipation of turbulence.  Hence the entropy of the non-thermal component can change while generating the equilibrium configuration.  In other words a cold component is not a zero entropy system (for example, we cannot use the adiabatic relation between the work done on an element and the internal energy within it).  Again, in analogy with compact objects, this is the basis for deteminning the composition of high density matter in neutron stars, while allowing for energy loss through the emission of neutrinos.  As mentioned above, in the absence of a first-princples model for the relation between magnetic fields, cosmic and turbulence, we simply assume that equipartition exists between these three components.  This simplification allows us to evaluate the entire non-thermal pressure based on a relation between one of these components and the star formation rate, and dictating the energy density in the other two components by equating it to the first. 

In accordance with this approach, we base our derivation on relations between magnetic fields and the star formation rate.  Specifically, we model the dependence of the magnitude of magnetic field, $B$, on the star formation rate, $\dot{\rho_*}$ (in mass per unit time per unit volume) as a power law: $B\propto \dot{\rho_*}^{\alpha_\text{1}}$. Combined with a [@schmidt59] power law relating the star formation rate to the gas density, $\rho$, i.e., $$\dot{\rho_*}=K\rho^\kappa,\label{eq:sfr}$$ we have a simple power law term of $$\label{eq:B_rho_alphaprime}
B\sim \rho^{\alpha_\text{2}},$$ for which $\kappa$, $\alpha_\text{1}$ and $\alpha_\text{2}=\kappa\alpha_\text{1}$ are all constants, which along with the proportionality factors must be constrained from observations.  Once this assumption has been made, the total non-thermal volumetric energy arising from these power laws takes the form: $E_\text{nt}=3 B^2/8\pi\propto\rho^\alpha$, $\alpha=2\alpha_\text{2}$ with the prefactor of $3$ originating from the contributions of the $3$ components in equipartition.  We reiterate that the equipartition is not a necessary assumption for this model.  Any constant distribution between the components is consistent with the assumptions and can be readily used.  The effective equation of state can easily be modified for further deviations from the“equipartition anzatz”, and any non-thermal pressure that is a monotonic function of density can be incorporated in an analogous way.  Essentially, the only requirement that cannot be simply generalized is that the non-thermal components are in quasi-steady state that depends on the gas density alone.  That is, that the non-thermal processes are related, and settle down on timescales which are short with respect to the evolution of galaxies. 

Stationary Non-thermal EoS {#sec:ssnteos}
--------------------------

We begin with the simplest model that incorporates the additional non-thermal pressure components.  In this model we assume that the energy in the non-thermal components is completely determined by the local instantaneous star formation rate.  Stipulating this assumption the energy in the non-thermal components becomes a simple function of $\rho$. The function is determined by various physical processes that contribute to the non-thermal pressure of the ISM, which depend differently on the density of the gas, but eventually materializes through equipartition.  At this point we will also assume that the gas is always in appropriate conditions so that star formation is turned on. 

Assuming a power-law dependence of $B$ (Equation \[eq:B\_rho\_alphaprime\]) and a Schmidt law for the star formation, the non-thermal pressure in equilibrium, $P_\text{nt}^\text{0}(\rho)$, is: $$P_\text{nt}^\text{0}(\rho)=
  A\rho^\alpha,\label{eq:Pnt}$$ where $A$ is a model dependent proportioanlity factor.  Here, and throughout the paper, we neglect the order-of-unity differences between pressure and energy density.  Magnetic field’s pressure depends on the field’s morphology and drops from $1$ for magnetic field in the disc’s plane, to $1/3$ for isotropically tangled field.  CR pressure equals $2/3$ of its energy density for non-relativistic particles.  We show below that the while the detailed analysis changes somewhat, the qualitative effect of the non-thermal pressure remains unaffected.  The total pressure at each point is the sum of the thermal and non-thermal pressure, $P_\text{tot}=P_\text{th}+P_\text{nt}^\text{0}$. As an aside we mention that the sound speed of the gas is striaghtforward to calculate in this effective EoS: For application in numerical codes, it is useful (for setting the timesteps according to the Courant conditions, for example) to calculate the numerical speed of sound of gas: $$\begin{aligned}
c_\text{s}^2&=\left(\frac{\partial P}{\partial \rho}\right)_\text{s}=\frac{\partial
  P_\text{nt}}{\partial \rho}+\left(\frac{\partial
  P_\text{th}}{\partial \rho}\right)_\text{s}\\
&=\alpha\frac{P_\text{nt}}{\rho}+\gamma\frac{P_\text{th}}{\rho}.\nonumber\end{aligned}$$ We note that it is not the physical speed of sound of the multi-component gas, which depends on the alfvénic velocity and largest eddy velocity in a non-trivial manner. 

Dynamic Non-thermal EoS {#sec:dnteos}
-----------------------

The EoS described in the previous subsection has the advantage of being extremely simple to implement since it requires a trivial addition to the ideal equation of state depending only on the gas density; There is no need to specifically trace the non-thermal component.  However, it suffers from undesired consequences that arise from the assumption that the non-thermal pressure traces the local instantaneuos star formation.  This means that if star formation was to suddenly begin (by passing the a threshold density, for example) or to suddenly end (perhaps if feedback blowing of the gas transfers it to a different thermodynamic regime where star formation is extinguished) a sudden jump in the pressure will follow, and potentially create spurious shocks and disturbances.  In addition, we know from observations that the magnetic field and cosmic ray vertical scale heights are considerably larger than those of the gas and star formation.  @ferriere01 estimates the magnetic scale height of the Galaxy at $\sim 1.4kpc$ based on rotation measures of pulsars [@Inoue1981] and the CR scale height of $\sim 2-4{\,{\rm kpc}}$ based on observed abundance of secondary particles such as $^{10}Be$ that are used to constrain various diffusion models for the Galaxy [@Garcia-Munoz1987; @Bloemen1993].  The magnetic field of distant galaxies is also observed via polarization and radio emission measurements at a vertical scale of a few kiloparsecs above the disc plane [see, for example @Krause2014].  The existence of non-thermal pressure based on radio observations and on a discrepancy between the vertical density and gravity profiles of the Galaxy perpendicular to the plane at the local neighborhood was advocated by [@boulares90; @cox05_ism].  The fact that magnetic fields and cosmic rays at regions that are far from star forming gas indicate that the coupling between the non-thermal components and star formation is more complicated than the simple assumptions in [§\[sec:ssnteos\]]{}.  Specifically, this suggests that the coupling is not instantaneous, but has a finite response time as energy convects with the gas or diffuses through it.  Alternatively, there may be additional sources for turbulence, magnetic fields and CR which are dominant at $\approx 1kpc$ altitudes above the disc (see, for example, @dekel09 for extra-galactic driving of turbulence and @Braun2012 for internal ISM instability-driven turbulence). 

To address this complication we introduce another independent variable into the equation of state so that the amount of non-thermal energy responds to star formation over a finite time.  We set $P_\text{nt}^\text{0}$ (Equation \[eq:Pnt\]) as the equilibrium value of non-thermal pressure for a given stellar density, and add a time-dependent form of the actual non-thermal pressure, $P_\text{nt}(t)$, which approaches $P_\text{nt}^\text{0}$ through some temporal dependence.  This adjustment expands the stationary EoS and includes a time integrated function of the star formation rate in the non-thermal pressure, thus ensuring both temporal and spatial continuity even if star formation flicks on and off.  Keeping with the spirit of our model, we do not attempt to describe the physics of the relaxation of the non-thermal component, and settle for a parametric description.  We note that this formulation also removes the numerical complications which arise from discontinuities (in the latter sense, this additional term has a stabilizing effect similar to the von-Neumann artificial viscosity which was introduced to help integrate over the non-smooth shock conditions).  We accomplish this by paramerterizing the non-thermal heating and cooling rate.  Since the two must cancel each other for a steady-state star formation, $P_\text{nt}^\text{0}$ must be an attractor of $P_\text{nt}(t)$ at a any given mass density: if $P_\text{nt}$ is too large then there should be a net cooling and vice versa. 

Non-thermal heating can be described by: $$\mathcal{H}_\text{nt}=f_\text{nt}~\mathsf{\epsilon}_\text{SN}~\eta_\text{SN}~\dot{\rho_*},\label{eq:heating}$$ where $\mathsf{\epsilon}_\text{SN}$ is the total energy injected into the gas per supernova, $\eta_\text{SN}$ is the number of supernovae per solar mass of stars that are created, and with $f_\text{nt}$ being the smaller than unity fraction of the supernova energy that ends up in the non-thermal components.  Non-thermal cooling is assumed to be responding to heating by the following characterization: $$\Lambda_\text{nt}=f_\text{nt}~\mathsf{\epsilon}_\text{SN}~\eta_\text{SN}
\left(\frac{P_\text{nt}}{P_\text{nt}^\text{0}}\right)^\beta~\dot{\rho_*},\label{eq:cooling1}$$ where $\beta$ is a free parameter which essentially controls the response time of the non-thermal components to changes in the star-formation rate.  The resemblance between the cooling term and the heating term arises from the requirement that the gas be in cooling/heating equilibrium at a star formation rate consistent with observations.  Combining these heating and cooling terms the time-dependent evolution of the non-thermal pressure at a constant gas density follows the simple form : $$\dot{P_\text{nt}}=\mathcal{H}_\text{nt}-\Lambda_\text{nt}=f_\text{nt}~\mathsf{\epsilon}_\text{SN}~\eta_\text{SN}
~\dot{\rho_*}\left[1-\left(\frac{P_\text{nt}}{P_\text{nt}^\text{0}}\right)^\beta\right].\label{eq:P_nt_oft}$$ Stability requires that: $$\label{eq:stability} \left.\frac{\partial\dot{P_\text{nt}}}{\partial
    P_\text{nt}}\right|_{P_\text{nt}^\text{0}}=-f_\text{nt}~\mathsf{\epsilon}_\text{SN}~\eta_\text{SN}\beta~\dot{\rho_*}<0,$$ so that $\beta>0$ ensures that the non-thermal pressure always approaches its asymptotic value for a steady star formation rate. 

The time dependent modification makes it possible to explicitly deal with a density threshold condition, as observed by @schmidt59 [@kennicutt98].  This condition cuts off star formation completely for gas densities below a critical value, $\rho_\text{c}$. We note that most numerical codes apply such a threshold (but for considerably lower threshold densities) also in order to prevent spurious star formation from occurring outside of galaxies.  Below this threshold, the steady-state non-thermal pressure is expected to vanish.  However, in reality, a non-star forming region can still maintain a steady state non-thermal pressure due to diffusion processes from neighboring regions [@Joung2009; @Scannapieco2012], or to various other sources [@Braun2012; @dekel09].  We partly account for that within our “single-cell” framework by setting a finite decay rate for non-star-forming regions.  In [eq.  (\[eq:cooling1\])]{}, we cannot set $P_\text{nt}^\text{0}=0$ for $\rho\leq \rho_\text{c}$, since the cooling rate then becomes ill defined.  We remedy this by formulating $P_\text{nt}^\text{0}$ as a function of $\rho$ (rather than $\dot{\rho_*}$) and redefining ([eq.  (\[eq:Pnt\])]{}) as follows:

$$\begin{aligned}
\label{eq:P_nt_0_new}
P_\text{nt}^\text{0} = A\rho^{\alpha '}&\text{ for } \rho>\rho_\text{c}\\
P_\text{nt}^\text{0} = A\rho_\text{c}^{\alpha '} &\text{ for } \rho<\rho_\text{c}.\nonumber\end{aligned}$$

Using [eq.  (\[eq:heating\])]{} and [eq.  (\[eq:cooling1\])]{} now assures that the heating turns off when no star formation occurs, and the cooling can proceed as the non-thermal pressure asymptotically approaches $0$. Combining [eq.  (\[eq:P\_nt\_oft\])]{} with [eq.  (\[eq:sfr\])]{} then yields: $$\dot{P_\text{nt}}=\mathcal{H}_\text{nt}-\Lambda_\text{nt}=f_\text{nt}~\mathsf{\epsilon}_\text{SN}~\eta_\text{SN}\left[
~\dot{\rho_*}-K\rho^\kappa\left(\frac{P_\text{nt}}{P_\text{nt}^\text{0}}\right)^\beta\right].\label{eq:P_nt_dot}$$

There is a noteworthy simplification in our model.  In order to achieve coarse-grained steady state of ISM gas two conditions must be met simultaneously: the total pressure must balance the external pressure, and the cooling must balance the heating.  For purely thermal pressure, these two equations are solved by varying two parameters - the density and temperature of the gas.  For gas with purely non-thermal pressure components of the type proposed in this work (which is a reasonable approximation for many external pressures, see [Fig.  \[fig:extp\]]{}), the pressure is a function of density alone, and the pressure equilibrium and heating/cooling equilibrium generally do not have a simultaneous solution.  We bypass this by relating the heating to the cooling in such a way that the observed relation is always achieved.  More advanced models which include physically motivated cooling and heating will introduce more dynamic parameters, allowing the gas to reach steady state more naturally. 

Dynamic Behaviour of the Non-Thermal EoS: A Quantitative Model {#sec:toymodel}
==============================================================

We now demonstrate the proporties and applicability of our effective EoS for non-thermal components with a point (zero-dimensional) model of the ISM.  In this model we evolve the conditions of a parcel of gas with isobaric boundary conditions, solving both thermal and non-thermal pressure components, in accordance with the models described in [§\[sec:eos\]]{}. 

Model Parameters
----------------

A quantitative implementation of our EoS requires the speficiation of the model’s free parameters.  For the stationary non-thermal pressure, these are the proportionality coeffecient and power which relate gas density to the pressure in magnetic fields ($P_\text{nt}^\text{0}=A\rho^\alpha$).  Even after applying our hypothesis of equaipartition among the non-thermal components, current uncertainties regarding the magnitude of magnetic fields in early galaxies, are quite large.  In essence, $(A,\alpha)$ may be treated as free parameters.  Some indication can, however, be gained from observed relations between star formation rates and synchrotron radiation [@Kennicutt83].  We chose to use the fits from equation A11 of @lacki10b of the form: $$B=B_\text{0}\Sigma_\text{gas}^a h^{-a},$$ where $\Sigma_\text{gas}$ and $h$ are galactic gas column densities and galactic scale-heights of galaxies.  These two parameters are fitted to observations using a one zone model for galaxies including the cosmic ray specra of primary and secondary rays, tracing self-consistently generation and evolution with effective diffusion coefficients [@lacki10a; @lacki10b; @Lacki2013].  The power law coefficient $a$, and normalization $B_\text{0}$, are observationally constrained by the FRC [far IR - radio correlation @Condon1991; @Yun2001] and by local measurements of CR and radio observations at $1.4GHz$ [see @lacki10a and reference within].  A model for the turbulent amplification rate of magnetic fields by star formation that recovers the FRC and predicts its breakdown was also suggested by @Schleicher2013.  Stipulating a typical vertical height for the magnetic fields of $1{\,{\rm kpc}}$ [@cox05_ism] one finds for the two fits suggested by @lacki10a [^3]: $$\begin{aligned}
  B&=6.65\left(\frac{\rho}{10^{-24}\text{gr}~\text{cm}^{-3}}\right)^{0.5}\mu G,\\
  B&=6.85\left(\frac{\rho}{10^{-24}\text{gr}~\text{cm}^{-3}}\right)^{0.6}\mu
  G.\nonumber\end{aligned}$$ These values are consistent with with observed Milky Way values [@ferriere01; @cox05_ism; @beck09] and with theoretical predictions [@lisenfeld96].  We convert those relations to non-thermal energy according to $P_\text{nt}=3~\frac{B^2}{8\pi}$ (again, recalling that the factor of 3 arises from the equipartition assumption) and find two similar (but not identical) realizations for the non-thermal EoS: $$\begin{aligned}
  P_\text{nt}^\text{0}(\rho)=5.3\times
  10^{-12}\left(\frac{\rho}{10^{-24}\text{gr}~\text{cm}^{-3}}\right)
  \text{erg}~\text{cm}^{-3}\label{eq:powerlaw1}\\
  P_\text{nt}^\text{0}(\rho)=5.6\times
  10^{-12}\left(\frac{\rho}{10^{-24}\text{gr}~\text{cm}^{-3}}\right)^{1.2}
  \text{erg}~\text{cm}^{-3}.\label{eq:powerlaw2}\end{aligned}$$ According to @lacki10a these fits reproduce the FIR-radio relations of galaxies equally well.  One could consider a case when the magnetic pressure scales linearly with density (corresponding to the first fit), or to the star formation (corresponding to a fit with $B\sim \rho^{0.75},$ or $B^2\sim \rho_*\sim \rho^{1.5}$) - but this parameter is not favored by the more detailed model there.  As we point below, in our model, the $B\sim rho^{0.5}$ fit is singular in the sense that its evolution can never derail it from static equilibrium once it has achieved.  For this reason, it is worthwhile to keep both fits at this stage. 

It is encouraging to note that these values are in rough agreement with the values needed to support the weight of the gas at the plane of the Galactic disc against its self gravity [@cox05_ism]. 

Other tracers of star formation could in principle be used to calibrate and constrain the non-thermal pressure terms.  X-rays are a good tracer for star formation as young OB stars emit X-rays [@ranalli03; @Mineo14] and at high redshifts [@Vattakunnel12].  However, X-rays are converted by neutral gas efficiently into UV and optical, and are never important as radiation pressure.  Further, X-rays are not a direct tracer for other sources of pressure of in the gas such as the amplitude of magnetic fields, turbulence or CR.  Measurements of turbulence by emission line broadening could, in principle, be used to calibrate the relation between star formation rate and the non-thermal components.  However, such measurements inevitably probe only some of the ISM gas phases, and measurements of high-redshift turbulence [@forster06] are few and might be driven by infalling gas rather than by star formation [@dekel09].  Alternative direct measurements of magnetic fields either trace the large scale magnetic fields of galaxies (polarization measurements) or the total line-of-sight magnetic fields (Faraday rotation measures) and are harder to correlate to total star formation than the synchrotron radiation used in this work. 

For the dynamic non-thermal EoS several additional parameters are required to define equations [eq.  (\[eq:heating\])]{} and [eq.  (\[eq:cooling1\])]{}.  The supernova energy $\mathsf{\epsilon}_\text{SN}$ and the supernova rate $\eta_\text{SN}$ can by taken from standard theories, but we do need to specify the parameters which control non-thermal heating and cooling, $f_\text{nt}$ and $\beta$. The fraction $f_\text{nt}$ sets the fraction of supernovae energy invested as non-thermal energy, while $\beta$ sets the power governing the rate at which the non-thermal energy reaches its equilibrium values.  These values are numerical by nature and should be set to allow the non-thermal energy to achieve equilibrium, while smoothing over pressure jumps arising from abrupt changes in the star formation rates.  As we show below, even for a low $f_\text{nt}$ of $0.1,$ the relaxation times for a wide choice of $\beta$ are shorter than a ${\,{\rm Myr}}.$ This value indicates that for smooth galactic histories the calculated state of the gas should be close to the asymptotic conditions constrained by observations. 

The Single Cell Isobaric Model
------------------------------

We incorporate our model for non-thermal pressure in a single cell model by tracing the evolution of gas under isobaric boundary conditions.  This represents a simplified behaviour of a single hydrodynamic cell embedded within a galaxy that evolves slowly and supports this cell with nearly constant external pressure.  For the general, dynamic case we solve the ordinary differential equations for the thermal internal energy of the gas and for the non-thermal pressure: $$\begin{aligned}
  \dot{e}_\text{th}=f_\text{th}\epsilon_\text{SN}\eta_\text{SN}\frac{\dot{\rho}_*}{\rho}-\Lambda_\text{t}(\rho,T)+P_\text{th}\frac{\dot{\rho}}{\rho^2},\label{eq:edot}\\
  \dot{P}_\text{nt}=f_\text{nt}~\mathsf{\epsilon}_\text{SN}~\eta_\text{SN}\left[
    ~\dot{\rho_*}-K\rho^\kappa\left(\frac{P_\text{nt}}{P_\text{nt}^\text{0}}\right)^\beta\right]+P_\text{nt}\frac{\dot{\rho}}{\rho}\label{eq:pntdot}.\end{aligned}$$ The last term in the right hand side of both equations is the contribution of the density change, $\dot{\rho}$, (a $PdV$ term for the energy equation).  The density is derived self-consistently by requiring that $$P_\text{nt}\left(\rho\right)+P_\text{th}\left(\rho,e_\text{th}\right)=P_\text{ext}.$$

We complement our model with an appropriate paramerterization of the star formation rate and the gas heating and cooling functions.  Star formation is modeled with a Schmidt law corresponding to a convention of $e_\text{eff}=5\%$ of the gas into stars every dynamic free-fall time of the gas: $$\dot{\rho_*}=e_\text{eff}\frac{rho}{t_\text{ff}}=e_\text{eff}\left(\frac{32G}{3\pi}\right)^{1/2}\rho^{3/2}.$$ For the sake of simplicity we begin with this star-formation rate with no density cutoff (we examine the implications of such a cutoff in section \[sec:cutoff\]).  The equation of state for the ideal gas implies $P_\text{th}=\left(\gamma-1\right)\rho e_\text{th}.$ In this work we allow the gas to cool according to the rates $\Lambda_\text{t}(\rho,T)$ from CLOUDY [version 96b4 @Ferland1998] by interpolating from tables described in @kravtsov03.  The cooling and heating of the gas includes Compton heating and cooling, redshift dependent UV heating and atomic and molecular cooling.  The tables provide the total cooling and heating and particle number as a function of the redshift, metalicity, density and temperature.  The temperature is related to the internal energy by integrating over the particle number dependent heat capacity as the number of particles changes by a factor of a few at recombination and at molecular formation.  For supernova heating we assume that one supernova occurs for every $160{M_\odot}$ of stars formed ($\eta_\text{SN}=1/160{M_\odot}$; see [eq.  (\[eq:heating\])]{}), corresponding to a Salpeter IMF between $0.1$ and $100{M_\odot}$ and supernovae occurring above $8{M_\odot}$ [@dobbs11].  We use a standard value for the average total energy released per supernova, ${\epsilon}_\text{SN}=10^{51}\text{erg}$. In most of the simulations described below we impose external pressure boundary conditions of $10^{-12}\text{erg}~\text{cm}^{-3}$ in rough correspondence with observed conditions in the plane of the Galactic disc [@cox05_ism]. 

For completeness we list the various definitions and default values for the coefficients in our equations in table \[tab:values\]. 

[ r r r | l ]{} Parameter & Units & Value & Definition\
\
\
\
$\gamma$ & & 5/3 & adiabatic constant\
$z$ & & 0 & redshift\
$Z$ & $Z_\odot$ & 1 & metalicity\
$\epsilon_\text{sfr}$ & & 0.05 & star formation efficiency\
$\eta_\text{SN}$ & & 1/160 & supernova per stellar mass formed\
$\epsilon_\text{SN}$ & erg & $10^{51}$ & supernova energy\
$f_\text{th}$ & & 0.1 & fraction of energy injected to thermal component\
\
\
\
$f_\text{nt}$ & & 0.2 & fraction of energy injected to thermal component\
$\alpha$ & & & density power law coefficient of non-thermal pressure ([eq.  (\[eq:Pnt\])]{})\
$A$ & see [eq.  (\[eq:Pnt\])]{} & & normalization of non-thermal pressure\
$\beta$ & & & cooling behaviour ([eq.  (\[eq:cooling1\])]{}\
\
\
\
$P_\text{ext}$ & $\text{erg} \, \text{cm}^{-3}$ & $10^{-12}$ & External pressure\
$P_\text{nt}^\text{i}$ & $\text{erg} \, \text{cm}^{-3}$ & & initial non thermal pressure\
$T^i$ & K & $10^5$ & initial gas temperature\

\[tab:values\]

Evolution of Gas with Purely Thermal Pressure {#sec:thermal}
---------------------------------------------

![ \[fig:thermal\_trho\] Time evolution of the temperature ([*top panel*]{}) and density ([*bottom panel*]{}) of the thermal-pressure-only models for varying supernovae efficiencies.  The pressure boundary conditions is $10^{-12}\text{dyn}~\text{cm}^{-2}$ and the initial temperature of the gas is $10^5$K. ](thermal_2t.eps){width="3.3in"}

![ \[fig:thermal\_rhodot\] Time evolution of the star formation rates and the depletion of gas into stars ($M_\text{gas}/M_\text{tot}$) for the model calculated in [Figure \[fig:thermal\_trho\]]{} . ]( thermal_1t.eps){width="3.3in"}

We start by demonstrating the properities of the single-cell simulation when only thermal pressure exists.  We set the initial conditions with a temperature of $T_\text{i}=10^5$K, which corresponds to a density of $\approx 10^{-25}\text{gr}~\text{cm}^{-3}.$ [Figure \[fig:thermal\_trho\]]{} demonstrates how cooling of the gas causes the temperature of the gas to decrease ([*top panel*]{}), and, correspondingly the isobaric boundary condition forces the density to scale as $1/T,$ forcing a density increase ([*bottom panel*]{}).  The bump at $\approx 0.1{\,{\rm Myr}}$ corresponds to the steep decrease in the cooling function at $T=10^4$K as gas becomes neutral and collisional excitation of lines becomes unimportant beyond this point. 

Greater densities enhance the star formation rate, as well as the resulting supernova feedback, and once the supernova feedback power balances the cooling rate the density and temperature of the gas become constant and the gas is converted into stars at a constant rate.  This is can be seen clearly in [Fig.  \[fig:thermal\_rhodot\]]{}, which depicts the specific star formation rate and the depletion of gas into stars.  This depletion ([Fig.  \[fig:thermal\_rhodot\]]{}, [*bottom panel*]{}) is calculated by noting that the mass in stars evolves as: $$M_\text{tot}=M_\text{gas}+M_*=M_\text{gas}+\int\dot{\rho_*}Vdt=\text{const},$$ where $M_\text{tot},M_\text{gas}$ and $M_*$ are the total, gas and stellar mass in our volume element, and $V$ and $\rho$ the time-dependent volume and density of the element.  In the single cell model $M_\text{tot}$ is fixed.  As gas is converted to stars the volume adjusts itself so the pressure corresponds to the external boundary condition.  This single cell assumption is self-consistently addressed in the full hydrodynamic implementation shown below ([§\[sec:ramses\]]{}).  Initially, $M_\text{tot}=M_\text{gas}=V_\text{0}\rho_\text{0}$ with $V_\text{0}$ and $\rho_\text{0}$ the initial volume and gas density, respectively.  $$\frac{V}{V_\text{0}}=\frac{M_\text{tot}-M_*}{V_\text{0}~\rho}=\frac{\rho_\text{0}}{\rho}-\frac{1}{\rho}\int\dot{\rho_*}\frac{V}{V_\text{0}}dt,\label{eq:vv0}$$ is an integral equation that can be evolved in time.  The depletion of gas is then shown as: $$\frac{M_\text{gas}}{M_\text{tot}}=1-\frac{M_*}{M_\text{tot}}=1-\frac{1}{\rho_\text{0}}\int
  \dot{\rho_*}\frac{V}{V_\text{0}}dt.\label{eq:depletion}$$ As is to be expected, once the density levels off at an equilibrium value, so does the specific star formation rate (the volume of the element continues to decrease over time in order to maintain a constant gas density with a decreasing mass). 

Our key observation is that without non-thermal components the gas achieves an equilibrium between the heating and the cooling after less than $1Myr$ and then converts most of the gas into stars quickly after that.  About $50\%$ of the gas is depleted during the first $100{\,{\rm Myr}}$ for low supernova efficiencies, and even when assuming perfect ($f_\text{SN}=1$) supernovae effeciencies, star formation has comsumed over one half of the gas by $400{\,{\rm Myr}}$. Absolute efficiency is certainly non-physical, since in reality most of the supernova energy gets converted into radiation that escapes the galaxy without contributing to pressure support of the gas.  In any case, we conclude that in a thermal-pressure-only model, boundary conditions corresponding to the pressure in the mid-plane of the Milky Way leads to gas being converted into stars over a few $100{\,{\rm Myr}}$ regardless of the efficiency of the thermal feedback. 

It is also noteworthy that since the gas achieves rate equilibrium during the first ${\,{\rm Myr}}$, the initial conditions of the gas do not affect the depletion time.  In [Figure \[fig:thermal2\_rhodot\]]{} we repeat the exercise with an initial temperature for the gas of $200$K, and the results are virtually unchanged, except that cooling/heating equilibrium is achieved after as little as $0.1{\,{\rm Myr}}$. We note that the unphysical $f_\text{SN}=1$ case reaches a different equilibrium point that exists on the molecular cooling branch at a lower temperature and leads to even faster gas depletion. 

![ \[fig:thermal2\_rhodot\] Same as [Figure \[fig:thermal\_rhodot\]]{}, but when the initial gas temperature is $200$K. ](thermal2_1t.eps){width="3.3in"}

Evolution of Gas with Thermal and Non-Thermal Pressure
------------------------------------------------------

We now proceed to examine the behaviour of a parcel of gas with similar boundary conditions as in [§\[sec:thermal\]]{}, but with additional non-thermal components, evolved according to [eq.  (\[eq:pntdot\])]{}. 

[Figure \[fig:ism\_trho\]]{} shows the evolution of a parcel of gas in terms of temperature, density and non-thermal pressure, again with a pressure boundary condition of $10^{-12}\text{dyn}~\text{cm}^{-2}$ and and initial temperature of $10^5$K. For the non-thermal pressure we use the parametrization described in [eq.  (\[eq:powerlaw1\])]{}. [ Fig.  \[fig:ism\_rhodot\]]{} describes the evolution of the specific star formation rate and gas depletion for the same model.  In all the simulations here the thermal supernova feedback is turned on with efficiency $f_\text{SN}=0.1$ as described in [§\[sec:thermal\]]{}, and the fraction of supernova energy that is injected into the non-thermal component here is $f_\text{nt}=0.2.$.

The green line shows the stationary ($P_\text{nt}=P_\text{nt}^\text{0}$) non-thermal EoS described in [§\[sec:ssnteos\]]{} (calculated by replacing [eq.  (\[eq:pntdot\])]{} with [eq.  (\[eq:Pnt\])]{}), and the blue, cyan and gray lines are for the dynamic non-thermal EoS ([§\[sec:dnteos\]]{}) with the relaxation power laws of $\beta=1$ and $5$ as indicated on the plot.  The blue and cyan lines correspond to models where we arbitrarily set a zero initial non-thermal pressure, $P_\text{nt}^\text{init}=0.$ This initial condition results in initial density of $\approx 10^{-25}\text{gr}~\text{cm}^{-3}$ as for the thermal case, whereas the gray line corresponds to an initial $P_\text{nt}^\text{init}$ which is in its steady-state value for an initial density, $\approx 6\times 10^{-26}\text{gr}~\text{cm}^{-3}$. Note that this value is only slightly below the initial value when non-thermal pressure is neglected. 

We note that the green and gray lines are identical: the power-law in [eq.  (\[eq:powerlaw1\])]{} is $\alpha=1,$ and when the gas is in equilibrium ($P_\text{nt}=P_\text{nt}^\text{0}$) its evolution according to [eq.  (\[eq:pntdot\])]{} is just $\dot{P}_\text{nt}=P_\text{nt}\frac{\dot{\rho}}{\rho}=A\dot{\rho}$ so it cannot evolve away from equilibrium once initially achieved.  We shall show below that the static evolution deviates from equilibrium initial conditions case when $\alpha\ne 1$ (Equation \[eq:powerlaw2\]). 

![ \[fig:ism\_trho\] Time evolution of the temperature (top panel) density (middle panel) and non-thermal pressure (bottom panel) of the gas.  Curves correspond to no non-thermal pressure (red), the steady-state non-thermal pressure (green), and dynamic non-thermal pressure (blue, cyan and gray); see text for detail - but note that in this case the green and gray lines overlap completely.  In all calculations the pressure boundary conditions is $10^{-12}\text{dyn}~\text{cm}^{-2}$ and the initial temperature of the gas is $10^5$K.](ism2_3t.eps){width="3.3in"}

![ \[fig:ism\_rhodot\] Time evolution of the star formation rates and the depletion of gas into stars ($M_\text{gas}/M_\text{tot}$) for the models calculated in [Figure \[fig:ism\_trho\]]{} . ]( ism2_1t.eps){width="3.3in"}

The two figures clearly demonstrate the distinct effect that non-thermal pressure has on the simulation.  For the initial conditions we set, the gas is initially supported (at least in part) by thermal pressure.  As the thermal energy is radiated away, temperature drops and the gas contracts, increasing the star formation rate.  However, the inclusion of non-thermal pressure removes the relation of $\rho
T\propto P_\text{tot}$, and introduces another degree of freedom.  The gas can then cool without a dramatic density increase, and so cooling does not necessarily lead to enhanced star formation.  In the calculations with the dynamic EoS the non-thermal component adjusts (increases) until a new stationary equilibrium is reached.  This equilibrium consists of a balance between supernovae feedback and cooling both for the gas and the non-thermal energy components, each separately.  Note that in all the simulations presented here the gas settles into this steady state in a few ${\,{\rm Myr}}$.

The shape of the the curves found with the dynamical EoS also deserves some elaboration.  Since the asymptotic non-thermal pressure is similar in all these cases, all trajectories with non-thermal pressure converge to the same values.  As gas contracts, its star formation and supernova rate increases, and, for the dynamic non-thermal EoS, it takes some time for the non-thermal reservoir to fill.  During this time the gas is actually under-pressurized with respect to its asymptotic values and the density is larger than its final value.  This overshoot is readily seen in the temperature and density of the gas of [Fig.  \[fig:ism\_trho\]]{} and in the specific star formation rate of [Fig.  \[fig:ism\_rhodot\]]{}.  The [*bottom panel*]{} of [Fig.  \[fig:ism\_trho\]]{} shows the gradual and monotonic increase of $P_\text{nt}$. The timescale for converging to the asymptotic value is set primarily by $f_\text{nt},$ and slightly depends on the relaxation power law $\beta.$ In all the runs here the gas quickly settles into a steady state for which the cooling is balanced by the thermal feedback and the total pressure is divided between the thermal component and the non-thermal component. 

The distinct effect of non-thermal pressure is easily seen by comparing the evolution in all of these calculations to the case in which non-thermal components are neglected, similar to [Fig.  \[fig:thermal\_trho\]]{} and [Fig.  \[fig:thermal\_rhodot\]]{} (shown for reference in [Figure \[fig:ism\_trho\]]{} and [Figure \[fig:ism\_rhodot\]]{} in red curves).  The outstanding feature is the dramatic difference in the asymptotic equilibrium between the two cases: the additional source of pressure allows the gas to cool without a dramatic density increase.  Hence, the equilibrium density is some $400$ times lower for the non-thermal case, and the temperature is $10$ times higher.  These different conditions lead to very different star formation rates and depletion times as can be observed in [Fig.  \[fig:ism\_rhodot\]]{}.  Once the non-thermal component is included, the equlibrium star formation rate is four order of magnitudes lower (corresponding to the Schmidt law used here that indicates $\dot{\rho_*}\propto\rho^{1.5}$).  Accordingly, the mass depletion time for the non-thermal EoS gas increases to $\approx
2{\,{\rm Gyr}}$ as opposed to about $100{\,{\rm Myr}}$ when only the thermal component in the pressure is included.  This increase in depletion times is related to the lower asymptotic density for this case.  Since the density approaches its asymptotic value much faster then the depletion time, most of the gas depletion occurs at the equilibrium density.  Hence the gas depletion time is $\tau_*\approx
M_\text{gas}/(V\dot{\rho}_*)=V\rho/(V\dot{\rho}_*)=\rho/\dot{\rho}_*$. For the Kennicutt Schmidt relation we use here, this leads to a $\tau_*\sim \rho^{-1/2}$ relation, so reducing the density by a factor of $400$ leads to a twenty fold increase of the depletion time. 

We also note the value of $\beta$ has a minor impact on the relaxation time scale in the dynamic models, which is ten to a few tens of Myr (as is to be expected, the model with $\beta=5$ has a shorter relaxation time than the one with $\beta=1$). 

In order to examine the robustness of the effects of the non-thermal EoS we repeat the calculations described here for (i) the parameters in [eq.  (\[eq:powerlaw2\])]{} and (ii) for a non-thermal component weaker by a factor of $3$ and $10$; $$P_\text{nt}=5.3\times 10^{-13}\left(\frac{\rho}{10^{-24}\text{gr}~\text{cm}^{-3}}\right) \text{erg}~\text{cm}^{-3},
\label{eq:powerlaw3}$$ instead of [eq.  (\[eq:powerlaw1\])]{}.  The results are described in [Fig.  \[fig:ism\_powerlaw2\]]{} and [Fig.  \[fig:ism\_equipartition\]]{}, respectively. 

![ \[fig:ism\_powerlaw2\] Same as [Figure \[fig:ism\_rhodot\]]{}, except that the non-thermal pressure is calculated with Equation ([eq.  (\[eq:powerlaw2\])]{}). ]( ism4_1t.eps){width="3.3in"}

![ \[fig:ism\_equipartition\] Time evolution of the star formation rates and the depletion of gas into stars ($M_\text{gas}/M_\text{tot}$) for various magnitudes of the non-thermal pressure based on the fit in [eq.  (\[eq:powerlaw1\])]{}) normalized to the full equipartition magnitude](ism10_1t.eps){width="3.3in"}

In contrast with the previous calculations, setting $\alpha\ne 1$ in [eq.  (\[eq:powerlaw2\])]{} causes the dynamic equation of state to slightly deviate from the static equation of state even when both calculations begin with initial conditions of $P_\text{nt}^\text{init}=P_\text{nt}^\text{0}.$ The actual difference in the evolution between these two cases is still small in comparison with the difference between them and calculations which begin with $P_\text{nt}^\text{init}=0.$, for which the gas density overshoots significantly, and peaks after about 0.5Myr (the dynamic calculation with a finite initial non-thermal pressure, shown in gray, does over shoot with the same time scale, but at a much smaller amplitude). 

Finally, [Fig.  \[fig:ism\_equipartition\]]{} shows that even when the non-thermal pressure is reduced by a factor of $10$, the gas depletion time is still a factor of five or so longer than the depletion time for any of pure thermal calculations with realistic thermal feedback effeciencies (see [§\[sec:thermal\]]{}).  Only a perfect thermal efficiency $\epsilon_\text{SN}=1.$ allows for a depletion time that is comparbale to the case when a weak non-thermal component is included.  This result emphasizes that non-thermal pressure is far more efficient in delaying gas depletion to star formation than enahncing thermal feedback from supernovae.  The time scale for relaxation in this weaked non-thermal pressure case is reduced, however, to one to a few Myr.  Our main conclusion is that any significant non-thermal pressure will inevitably lead to a large change in the gas depletion time when compared to pure-thermal pressure models.  We infer that this is a general consequence of non-thermal pressure, regardless of whether equipartition is assumed. 

Evolution with a Cutoff Density for Star Formation {#sec:cutoff}
--------------------------------------------------

![ \[fig:ism\_sfcutoff\] Same as [Figure \[fig:ism\_trho\]]{}, but with density dependent star formation cutoff introduced at $\rho=10^{-24}\text{gr}\;\text{cm}^{-3}$. ](ism6_3t.eps){width="3.3in"}

![ \[fig:extp\] The gas density ([*top panel*]{}) and the fraction of the non-thermal pressure from the total pressure ([ *bottom panel*]{}) as a function of the external pressure boundary condition for the parameters described in table \[tab:values\].  Below $10^-{11.4}~\text{erg}/\text{cm}^3$ the density is always the cutoff density for star formation ($10^{-24}\text{gr}/\text{cm}^3$ here).  Above this pressure the density increases as a power-law. ]( extp.eps){width="3.3in"}

The existence of a star formation density threshold is predicted by @Kennicutt83 and is present in essentially all numerical models of galaxy formation.  It is typically implemented by invoking a single numerical value, set to compensate for the inability to simulate star formation and reach the necessary (high) densities.  Moreover, the thershold is applied to prevent spurious star formation outside of galaxies.  In this subsection we study the response of our model to an inclusion of such a fiducial threshold.  We use a threshold value of $10^{-24}\text{gr}\;\text{cm}^{-3},$ which is typical in cosmological simulations.  We emphasize that this is a qualitative demonstration of the effect which is included self-consistently in any numerical simulation, such as those described below, in [§\[sec:ramses\]]{}. 

[Figure \[fig:ism\_sfcutoff\]]{} describes the temporal evolution of a gas parcel for a star formation law that includes a sharp cutoff for densities below $10^{-24}\text{gr}\;\text{cm}^{-3}.$ This threshold is higher than the equilibrium density calculated without the threshold ([Fig.  \[fig:ism\_trho\]]{}) so there is no strict equilibrium solution (i.e. static solution) for this case in which the cooling balances the heating at all times.  Instead, we find that the qualitative behaviour of the system is such, that the response time of the non-thermal components creates a cycle in which star formation flickers on and off and the time-averaged heating balances the continuous cooling.  By construction, the single cell model is clearly inadequate for a quantitative study of this duty cycle, because it coarse grains over the relevant spatial and temporal scales necessary .  We do confirm numerically that our integration does indeed flicker. 

It is noteworthy that applying the density cutoff does not imply that the star formation in a full simulation in a galactic ISM will occur at constant density.  Non-homogeneity in the ISM is expected [see for example @Ostriker2001] and implies that the external pressure boundary conditions should vary in space and time.  We demonstrate that this principle applies also in the case of non-thermal pressure by examining the dependence of the equilibrium density and of the star formation rate on the external pressure conditions.  We vary the latter and solve the equlibrium density and non-thermal pressure for the non-thermal relations in [eq.  (\[eq:powerlaw1\])]{} and [eq.  (\[eq:powerlaw2\])]{}.  Our results are presented in [Fig.  \[fig:extp\]]{}.  We find that the density cutoff imposes a transition that depends on the external pressures: for low external pressures, the equilibrium density settles at the cutoff density for star formation as described above.  However, for $P_\text{ext}\geq 10^{-11.4}~\text{erg/cm}^3$ star formation becomes possible and the equilibrium density is larger (see the top panel of [Fig.  \[fig:extp\]]{}).  We note that the non-thermal pressure dominates for practically any external pressure above this tansition value (lower panel in [Fig.  \[fig:extp\]]{}), so the density in this regime essentially scales as $\rho\propto P_\text{ext}^{1/\alpha}$. It is encouraging to note that for pressures that correspond to the plane of the Galaxy ($P\approx 3 \times 10^{12}\text{erg}/\text{cm}^3$) the relative contribution of the thermal component is a few per-cent,which is in agreement with observations [see fig.  2 of @cox05_ism and accompanying text]. 

The effect of non-thermal ISM equation of state on realistic galaxies {#sec:ramses}
=====================================================================

We test our model by implementing it on isolated spiral galaxy simulations ran on [[*RAMSES*]{} ]{}[@Teyssier2002].  In the following section we will describe in some detail the non-trivial aspects of our implementation ([§\[sec:implementation\]]{}), the simulations that were ran ([§\[sec:ramses\_ic\]]{}), and describe the effects of non-thermal feedback on the star formation history and on the morphology of the resulting galaxies ([§\[sec:ramses\_results\]]{}). 

Model implementation {#sec:implementation}
--------------------

We now describe the numerical methods we have used to solve for the Euler equations in presence of a non-thermal energy components.  The original equations have to be modified by adding to the total fluid energy the non-thermal energy and to the total pressure the non-thermal pressure.  The modified equations now read $$\frac{\partial \rho}{\partial t} + \nabla \cdot \left( \rho {\bf u}\right) = 0$$ $$\frac{\partial }{\partial t} \left( \rho {\bf u}\right) + \nabla \cdot \left( \rho {\bf u} \otimes {\bf u} + P_\text{tot} {\mathbb I}\right) = - \rho \nabla \Phi$$ $$\frac{\partial E_\text{tot}}{\partial t} + \nabla \cdot \left(  {\bf u} \left( E_\text{tot} + P_\text{tot }\right)\right) = - \rho {\bf u}\cdot \nabla \Phi$$ $$\frac{\partial e_\text{nt}}{\partial t} + \nabla \cdot \left(  {\bf u} e_\text{nt} \right)  + P_\text{nt } \nabla \cdot {\bf u} = 0$$ The total fluid energy is now defined as $$E_\text{tot}=\frac{1}{2}\rho u^2 + e + e_\text{nt}$$ and the total fluid pressure as $$P_\text{tot}=P + P_\text{nt}$$ where the thermal pressure is given by the EoS of the thermal component $$P=(\gamma-1)e,$$ and the non-thermal pressure by the EoS of the non-thermal component $$P_\text{nt}=(\gamma_\text{nt}-1)e_\text{nt},$$ For the thermal component only, we can also define the specific thermal energy $\epsilon$ as $$e=\rho \epsilon.$$ We see in the previous equation that the internal energy of the thermal component is obtained by subtracting from the total energy the other energy components, namely $$e=E_\text{tot}- e_\text{nt}-\frac{1}{2}\rho u^2.$$ In the hydrodynamics solver, we have to modify several components of the code to add this non-thermal energy variable.  First, the predictor step in our MUSCL scheme [@Teyssier2002; @Fromang2006] is augmented by an additional equation for the non-thermal pressure.  Second, in the same predictor step, the non-thermal pressure is added to the thermal pressure in the equation governing the velocity update.  The next important correction is for the Riemann solver, used to define the flux at cell interface, as a function of the left and right states interpolated with in space and time at the interface position.  We have modified our various approximate Riemann solver by just replacing the fluid energy and the fluid pressure by the total energy and the total pressure.  The sound-speed in the augmented hyperbolic system of (quasi-) conservation laws has to be modified as $$c^2_\text{s,tot}=\frac{\gamma P + \gamma_\text{nt}P_\text{nt}}{\rho}$$ We have tested successfully our new algorithm on simple shock tubes featuring the additional non-thermal energy.  An important point we would like to stress is that in our simple model, shock heating occurs only for the thermal component.  Since in the previous set of equations, there is no sources of non-thermal energy at shock fronts, and no coupling between the two energies, the evolution of the non-thermal component is strictly adiabatic. 

For our non-thermal pressure model described in this paper we allow for one additional scalar component of non-thermal energy with $\gamma_\text{nt}=2,$ and add a source term for it to the feedback routine according to the heating portion of [eq.  (\[eq:P\_nt\_dot\])]{} and a sink term to the baryonic cooling routine according to the cooling portion of that same equation. 

Simulation parameters {#sec:ramses_ic}
---------------------

Feedback in hydrodynamical simulations is typically ineffective in regulating star formation and reducing gas depletion times.  To overcome that, a combination of methods, all “pumped up” to be as efficient as physically possible, is used.  We apply here the ’standard’ tools used in the [[*RAMSES*]{} ]{}runs of the AGORA [@kim2014] isolated galaxy.  Our base-line simulation, “Standard sim” uses the standard tools used in RAMSES which include several ad-hoc measures calibrated to prevent reaching high density and short depletion times.  These methods include delayed cooling (preventing cooling for a period of time after feedback energy injection to account for the adiabatic phase of the Sedov-Taylor explosion).  Additional methods applied is to increase the stochasticity of the process by allowing feedback to operate according to a Poisson distribution with a typical mass scale of a giant molecular cloud (GMC) and, although not strictly motivated by feedback, incorporating a pressure floor for the ISM gas that prevents it from reaching extremely cold and dense states which would imply very large star formation rates.  We aim here to demonstrate that the parameters of these methods can be relaxed once our non-thermal feedback model is used.  For comparison we consider two more simulations.  In the first, delayed cooling is turned off and the simulation indeed exhibits an over-production of stars.  In the third simulation we introduce our non-thermal feedback and demonstrate its ability to reduce star formation without the delayed cooling model. 

The three simulations are defined in table \[tab:sims\].  The simulations were ran using the AGORA low resolution initial conditions for a Milky way-like galaxy.  The setup consists of a DM halo of $10^{12}{M_\odot}$, stellar disk of $3.4\times 10^{10}{M_\odot}$ and a gaseous disk of $8.5\times 10^9{M_\odot}.$ The maximal refinement level is $12$, with a box size of $400{\,{\rm kpc}}$ and a maximal resolution of $\sim
200{\,{\rm pc}}.$ Particle mass was $3\times 10^5{M_\odot}$ for the stellar component and $10^7{M_\odot}$ for the halo dark matter particles.  Execution time until a physical time of $2{\,{\rm Gyr}}$ on $240$ cores on the ICPL[^4] cluster took a few days. 

In all simulation the GMC mass is $6.4\times 10^6{M_\odot},$ and the pressure floor is defined as a minimal temperature of $T_\text{min}=10^4(n/0.1)^{2/3}$K. For the ’standard feedback’ simulation the delayed cooling timescale is $20{\,{\rm Myr}}.$ For the non-thermal feedback simulation the amount of energy that is injected into the non-thermal component is half of the energy that is injected into the thermal part, the asymptotic value for the non-thermal component was set by [eq.  (\[eq:powerlaw1\])]{} and the rate of energy dissipation was set to $\beta=1/4$ ([eq.  (\[eq:cooling1\])]{}). 

[ r r r ]{} Name & delayed cooling & non-thermal pressure\
\
Standard sim& on & off\
Weak sim & off & off\
Non-thermal sim & off & on\

\[tab:sims\]

![ \[fig:mw2a\] Density colormap of the central $40{\,{\rm kpc}}$ of our “Standard sim” simulation at $t=1100{\,{\rm Myr}}$. [*Top panel*]{} edge on, [*bottom panel*]{} face on. ]( mw.eps){width="3.3in"}

![ \[fig:histnt\] Density-Temperature histograms of “Non-thermal sim”. [* Panel (a)*]{} presents the thermal component, and [*panel (b)*]{} the non-thermal component at time $t=1100{\,{\rm Myr}}.$ The temperature of the ’non-thermal’ component is defined as the non-thermal pressure divided by the density.](histnt.eps){width="3.3in"}

![ \[fig:ppnt\] The spatial extent of the non-thermal pressure at $t=1100{\,{\rm Myr}}.$ The colormap represents the ratio between the non-thermal pressure to the total pressure.](ppnt.eps){width="3.3in"}

Results {#sec:ramses_results}
-------

Some results of the simulations described in [§\[sec:ramses\_ic\]]{} are presented in [Fig.  \[fig:mw2a\]]{}-[Fig.  \[fig:sfrcomp\]]{}. [ Fig.  \[fig:mw2a\]]{} shows a head-on and edge-on view of “Standard sim”.  The simulation produces a stable rotating disc, and strong outflowing winds that are noticeable as the rough “x” shaped overdensity extending diagonally to the edges of our $40{\,{\rm kpc}}$ box.  We note that in all these simulations there is no hot halo gas component, that is expected to strongly affect outflows. 

[Fig.  \[fig:histnt\]]{} presents a density-temperature histogram of our “Non-thermal sim”.  The top and bottom panels show the thermal and non-thermal components.  For the non-thermal component we define the effective temperature as the non-thermal pressure divided by the density.  In the histograms, isobaric lines would be represented by lines with slope of $-1$ ($P\sim n \times T=\textrm{const}$) and the pressure floor is the sharp cutoff diagonal through the lower right part of the plot with slope of $2/3.$ The non-thermal temperature is typically higher than the thermal pressure.  For the formulation of the non-thermal pressure used here ([eq.  (\[eq:powerlaw1\])]{}), this asymptotic temperature is independent of density and is simply $P/\rho
m_\text{p}/k_\text{B}=6.4\times 10^4$K, with $m_\text{p}$ the proton mass and $k_\text{B}$ Boltzmann’s constant.  For high densities, when the non-thermal cooling is efficient enough for gas to relax into the asymptotic value, this value roughly coincides with the effective temperature for the non-thermal component.  However, the scatter due to the stochasticity of the feedback process is very large. 

[Fig.  \[fig:ppnt\]]{} shows an edge-on view of the galaxy with the colormap representing the fraction of the non-thermal pressure to the total pressure.  It is evident that dynamically significant non-thermal pressure exists around the disk at a distance of a few kiloparsecs, consistent with magnetic field observations [@ferriere01]. 

![ \[fig:hists\] Density-Temperature histograms of the three models for feedback at time $t=1100{\,{\rm Myr}}.$ [*Panels (a), (b) and (c)* ]{} show results from “Standard sim”, “Weak sim” and “Non-thermal sim” respectively.  The black ellipses point to the star forming region on the n-T plot.](hist3pan.eps){width="3.3in"}

![ \[fig:sfrcomp\] Star formation histories ([*top panel*]{}) and star formation rates ([*bottom panel*]{}) for simulations with the three feedback models.](sfrcomp.eps){width="3.3in"}

[Fig.  \[fig:hists\]]{} compares the thermodynamic state of the gas of the various simulations.  Stars are formed through the high density branch with a slope of $2/3$, emphasized by the black ellipse in the plots.  While this line is truncated in “Standard sim” ([*panel (a)*]{}), at $n\sim 3\text{cm}^{-3},$ is continues beyond $10 \text{cm}^{-3}$ for “Weak sim” ([*panel (b)*]{}) which causes an increase in star formation (see [Fig.  \[fig:sfrcomp\]]{}).  The third panel is similar to the top panel of [Fig.  \[fig:histnt\]]{}.  When non-thermal feedback is added, the amount of high-density gas is reduced.  For the non-thermal case ([*panel (c)*]{}) it is evident that there is a considerable fraction of the gas at intermediate densities ($10^{-4}-10^{-1} \text{cm}^{-3}$) and $~10^4$K temperature that is absent from the two purely thermal simulations.  This is a result of the isochoric cooling process that occurs in that case (corresponding to purely vertical thermodynamic trajectories on the histogram plot) versus the more efficient isobaric cooling for the thermal cases (corresponding to diagonal trajectories going down and to the right, with slope of $-1$, as explained above).  The extra non-thermal pressure supports the gas and allows it to cool without contracting to the pressure floor.  In “Standard sim” ([*panel (a)*]{}), most of the gas has been blown out of the galaxy to $T=10^6-10^7$K and $n=10^{-7}-10^{-4}$ which effectively shut down star formation (see [Fig.  \[fig:sfrcomp\]]{}).  The non-thermal pressure ([*panel (c)*]{}) affects intermediate densities and keeps the gas pressurized with pressure corresponding to $\sim 10^5$K, without blowing it out of the galaxy altogether. 

[Fig.  \[fig:sfrcomp\]]{} shows the star formation history and star formation rates for the three simulations.  The most efficient quenching occurs for “Standard sim” and is a result of massive blowout of gas from the galaxy. “ Weak sim” depletes most of the gas rapidly, within the first $\sim 500{\,{\rm Myr}}$ of the simulation, demonstrating the overcooling catastrophe that occurs for insufficient feedback. “ Non-thermal sim” exhibits regulated star formation that produces stars at $\sim
1{M_\odot}/{\,{\rm yr}}$ throughout the simulation. 

Not surprisingly, our results indicate that non-thermal pressure can have a significant effect on star formation and prolong the depletion times of galaxies to the observed timescales.  Clearly further research will be required to quantify the relative importance of the non-thermal components to other feedback mechanisms that are applied.  Regardless, we emphasize that observational evidence, as well as robust theoretical motivation, point to the existence of this component and that, in one form or another, it should be incorporated in cosmological and galactic simulations. 

Summary and discussion {#sec:summary}
======================

On galactic scales, the interstellar medium exists at quasi-static pressure that is required to support the atmosphere above it.  In an equilibrium configuration, loss of pressure due to cooling processes is balanced by heating, which for typical disc galaxies at low redshifts, is dominated by stellar feedback.  Stellar feedback, through its dependence on star formation rate is related to the ISM gas density. 

It is a well-known result that when only thermal pressure is considered in simulations, the resulting ISM density constrained by the pressure and heating-cooling equilibria leads to relatively large star formation rates and short ($\sim 100$Myr) gas depletion time scales.  This is considerably faster than depletion times of $\approx
1{\,{\rm Gyr}}$ inferred from observations [@Kong2004; @Bauermeister2013; @Pflamm-Altenburg2009; @Tacconi2013].  In this work we consider the contribution of non-thermal pressure components to this picture.  Non-thermal pressure consists of turbulence, cosmic rays and magnetic fields, and we examine their impact in an effective model.  Current cosmological simulations generally do not include the later two, and do not always resolve turbulence.  We demonstrate that non-thermal pressure components can be instrumental in solving the depletion time discrepancy in two respects: they reduce the quasi-steady state density and the corresponding star formation rates and cooling times, and they stabilize the gas by adding longer relaxation times in cases where star formation flickers on and off.  The regulating effect has been shown previously for cosmic rays [@Salem2014; @Booth2013; @Hanasz2013] and turbulence [@Ostriker2001; @Braun2012] but the two were not considered together and in any case were not yet formulated in a way which is applicable to large scale cosmological simulations. 

To test our assumptions we construct a simplified physical model for which all the non-thermal components achieve a steady state that is solely a function of density.  While simplistic, the advantage of such an approach is that it is readily applicable in numerical simulations.  Furthermore, we calibrate this density dependence by using the observed relations between the star formation rate for various galaxy observations and the synchrotron radiation, so that the magnitude of the effect is reasonably constrained.  To study its effect we first implement it into a single-zone numerical model that traces the evolution of a parcel of star forming gas with varying physical conditions under isobaric boundary conditions that mimic the pressure confinement of the gas by the atmosphere around and above it.  Using this mode we find that for a given, realistic, thermal feedback the depletion times naturally grow from $\approx 100{\,{\rm Myr}}$ to $\approx 2{\,{\rm Gyr}}$ in better agreement with observations, and that the coarse grained density of the gas is reduced by several orders of magnitude.  Then, the model is implemented into the hydrodynamical code [[*RAMSES*]{} ]{}and we present three simulations of the same isolated Milky-way-like galaxy with three different physical models.  In the first we use some of the “standard” recipes generally used for feedback.  Using that model reduces star formation by blowing the gas to high temperatures and low densities, and expelling it from the galaxy.  Then, to demonstrate the problem we deliberately turn off one of the key feedback components - the delayed cooling and show that gas cools and accumulates at the numerical pressure floor with high densities that cause a depletion of the stars within $\sim 100{\,{\rm Myr}}$ - at odds with observed depletion times of $1-2{\,{\rm Gyr}}.$ In the third simulation we introduce the non-thermal model calibrated by the observed radio FIR relation.  For that model, the gas remains pressurized at intermediate temperatures and densities, reducing the non-physical low density gas of the first model, and the non-physical high density star forming gas of the second model.  This model is effective in regulating star formation for a long period of time ($\sim 1{\,{\rm Gyr}}$) without blowing the gas out of the galaxy altogether. 

This work is a natural first step in incorporating non-thermal pressure components in galactic scale simulation.  The next steps can, and should, pursue several avenues of research.  The first is to better model the various non-thermal components, their internal interaction and their interaction with the thermal component and star formation.  This would relax the assumption of equipartition, and replace the observational constraints with more physically motivated ones.  For this step, calibration against results from ISM-scale hydrodynamic simulations will be beneficial.  A different avenue to pursue, in tandem or separately, is to use our [[*RAMSES*]{} ]{}patch to run large scale numerical simulations and to demonstrate its effect and applicability on cosmic scales.  Cosmological simulations today generally do not include the magnetic fields and cosmic rays, and do not always resolve turbulence, and our approach allows to circumvent this difficulty by using a simple effective parametrization.  Such simulations will naturally include realistic boundary conditions for the ISM, namely the halo gas, and allow us to study its interactions with the ISM and its effect on winds.  Ultimately, once large scale cosmological simulations are possible with all the necessary physics, a toy sub-grid model can also be calibrated directly to those simulations and used as a cheaper approximation for them. 

Acknowledgements {#acknowledgements .unnumbered}
================

We thank Andrey Kravtsov for making the cooling tables available to us for the single cell calculations.  Computational resources were provided through ICPL (<http://ICPL.HUJI.AC.IL>). 

A., [Blitz]{} L., [Bolatto]{} A., [Bureau]{} M., [Leroy]{} A., [Ostriker]{} E., [Teuben]{} P., [Wong]{} T., [Wright]{} M., 2013, , 768, 132

R., 2009, IAU Symposium, 259, 3

R., [Brandenburg]{} A., [Moss]{} D., [Shukurov]{} A., [Sokoloff]{} D., 1996, , 34, 155

E. F., 2003, , 586, 794

J. B. G. M., [Dogiel]{} V. A., [Dorman]{} V. L., [Ptuskin]{} V. S., 1993, , 267, 372

C. M., [Agertz]{} O., [Kravtsov]{} A. V., [Gnedin]{} N. Y., 2013, , 777, L16

A., [Cox]{} D. P., 1990, , 365, 544

H., [Schmidt]{} W., 2012, , 421, 1838

J. J., [Huang]{} Z.-P., [Yin]{} Q. F., [Thuan]{} T. X., 1991, , 378, 65

D. P., 2005, , 43, 337

R. A., [Theuns]{} T., [Dalla Vecchia]{} C., [Eke]{} V. R., [Frenk]{} C. S., [Jenkins]{} A., [Kay]{} S. T., [Peacock]{} J. A., [Pearce]{} F. R., [Schaye]{} J., [Springel]{} V., [Thomas]{} P. A., [White]{} S. D. M., [Wiersma]{} R. P. C., 2009, , 399, 1773

M. A., [Breitschwerdt]{} D., 2005, , 436, 585

T., [Klein]{} U., [Wielebinski]{} R., [Wunderlich]{} E., 1985, , 147, L6

A., [Birnboim]{} Y., [Engel]{} G., [Freundlich]{} J., [Goerdt]{} T., [Mumcuoglu]{} M., [Neistein]{} E., [Pichon]{} C., [Teyssier]{} R., [Zinger]{} E., 2009, , 457, 451

S., [Bell]{} E., [Burkert]{} A., 2006, , 638, 797

C. L., [Burkert]{} A., [Pringle]{} J. E., 2011, , 417, 1318

Y., [Teyssier]{} R., 2008, , 477, 79

B. G., [Scalo]{} J., 2004, , 42, 211

C., [Chabrier]{} G., [Schober]{} J., [Banerjee]{} R., [Klessen]{} R. S., [Schleicher]{} D. R. G., 2011, Physical Review Letters, 107, 114504

G. J., [Korista]{} K. T., [Verner]{} D. A., [Ferguson]{} J. W., [Kingdon]{} J. B., [Verner]{} E. M., 1998, , 110, 761

K. M., 2001, Reviews of Modern Physics, 73, 1031

N. M., [Genzel]{} R., [Lehnert]{} M. D., [Bouch[é]{}]{} N., [Verma]{} A., [Erb]{} D. K., [Shapley]{} A. E., et al.,  2006, , 645, 1062

S., [Hennebelle]{} P., [Teyssier]{} R., 2006, , 457, 371

M., [Simpson]{} J. A., [Guzik]{} T. G., [Wefel]{} J. P., [Margolis]{} S. H., 1987, , 64, 269

F., [Willman]{} B., [Mayer]{} L., [Brooks]{} A., [Stinson]{} G., [Valenzuela]{} O., [Wadsley]{} J., [Quinn]{} T., 2007, , 374, 1479

M., [Lesch]{} H., [Naab]{} T., [Gawryszczak]{} A., [Kowalik]{} K., [W[ó]{}lta[ń]{}ski]{} D., 2013, , 777, L38

P. F., [Kere[š]{}]{} D., [Murray]{} N., 2013,

P. F., [Quataert]{} E., [Murray]{} N., 2012, , 421, 3488

M., [Tabara]{} H., 1981, , 33, 603

M. R., [Mac Low]{} M.-M., [Bryan]{} G. L., 2009, , 704, 137

R., 1983, , 120, 219

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

C.-G., [Ostriker]{} E. C., [Kim]{} W.-T., 2013, , 776, 1

J.-h., [ Abel]{} T., [Agertz]{} O., et al.  2014, , 210, 14

X., 2004, , 425, 417

M. J., [Brandenburg]{} A., [Shukurov]{} A., [Tuominen]{} I., [Nordlund]{} [Å]{}.,  1999, , 514, L99

H., [Ostriker]{} E. C., 2009, , 693, 1316

M., 2014, ArXiv:1401.1317

A. V., 2003, , 590, L1

R., [Pearce]{} W. P., 1969, , 156, 445

B. C., 2013, ArXiv:1308.5232

B. C., [Beck]{} R., 2013, , 430, 3171

B. C., [Thompson]{} T. A., 2010, , 717, 196

B. C., [Thompson]{} T. A., [Quataert]{} E., 2010, , 717, 1

U., [Voelk]{} H. J., [Xu]{} C., 1996, , 314, 745

M. S., 1994, [High energy astrophysics.  Vol.2: Stars, the galaxy and the interstellar medium]{}.  Cambridge University Press

M.-M., [Klessen]{} R. S., 2004, Reviews of Modern Physics, 76, 125

A., [Iapichino]{} L., [Schmidt]{} W., [Niemeyer]{} J. C., 2009, , 707, 40

C. F., [Ostriker]{} J. P., 1977, , 218, 148

S., [Gilfanov]{} M., [Lehmer]{} B. D., [Morrison]{} G. E., [Sunyaev]{} R., 2014, , 437, 1698

J. F., [White]{} S. D. M., 1993, , 265, 271

B. D., [Dav[é]{}]{} R., 2006, , 373, 1265

E. C., [Stone]{} J. M., [Gammie]{} C. F., 2001, , 546, 980

J., [Kroupa]{} P., 2009, , 706, 516

P., [Comastri]{} A., [Setti]{} G., 2003, , 399, 39

B. E., [Kravtsov]{} A. V., 2008, , 680, 1083

M., [Bryan]{} G. L., 2014, , 437, 3312

C., [Wadepuhl]{} M., [Parry]{} O. H., [Navarro]{} J. F., [Jenkins]{} A., [Springel]{} V., [Teyssier]{} R., [Carlson]{} E., [Couchman]{} H. M. P., [Crain]{} R. A., [Dalla Vecchia]{} C., [Frenk]{} C. S., et al.,  2012, , 423, 1726

E., [Br[ü]{}ggen]{} M., 2010, , 405, 1634

J., [Dalla Vecchia]{} C., 2008, , 383, 1210

D. R. G., [Beck]{} R., 2013, , 556, A142

M., 1959, , 129, 243

W., [Federrath]{} C., [Hupp]{} M., [Kern]{} S., [Niemeyer]{} J. C., 2009, , 494, 127

W., [Kern]{} S. A. W., [Federrath]{} C., [Klessen]{} R. S., 2010, , 516, A25

V., [Hernquist]{} L., 2003, , 339, 289

R., [Shukurov]{} A., [Fletcher]{} A., [Beck]{} R., [La Porta]{} L., [Tabatabaei]{} F., 2014, , 437, 2201

L. J., [Neri]{} R., [Chapman]{} S. C., [Genzel]{} R., [Smail]{} I., [Ivison]{} R. J., [Bertoldi]{} F., [Blain]{} A., [Cox]{} P., [Greve]{} T., [Omont]{} A., 2006, , 640, 228

L. J., [Neri]{} R., [Genzel]{} R., [Combes]{} F., [Bolatto]{} A., [Cooper]{} M. C., [Wuyts]{} S., [Bournaud]{} F., [Burkert]{} A., [Comerford]{} J., [Cox]{} P., [Davis]{} M., [F[ö]{}rster Schreiber]{} N. M., [Garc[í]{}a-Burillo]{} S., et al.,  2013, , 768, 74

R., 2002, , 385, 337

S., [Tozzi]{} P., [Matteucci]{} F., [Padovani]{} P., [Miller]{} N., [Bonzini]{} M., [Mainieri]{} V., [Paolillo]{} M., [Vincoletto]{} L., [Brandt]{} W. N., [Luo]{} B., [Kellermann]{} K. I., [Xue]{} Y. Q., 2012, , 420, 2190

M. S., [Reddy]{} N. A., [Condon]{} J. J., 2001, , 554, 803

\[lastpage\]

[^1]: We distinguish between ordered magnetic fields which slowly accumulate over the lifetime of a galactic discs due to galactic-scale dynamo effect, and tangled magnetic fields on scale of a few parsecs and below that is related to, and correlates with, the star formation [@beck96; @deavillez05].Throughout this paper we shall only be concerned with the latter. 

[^2]: Throughout the paper we assume the flux-frozen approximation for magnetic fields and neglect magnetic reconnection and ambipolar diffusion, as is appropriate for the densities and ionizations of the ISM gas. 

[^3]: We note that these scale-heights are higher than the ones used by @lacki10a.  Using a smaller scale-height would result in higher magnetic fields and even higher non-thermal pressure. 

[^4]: <http://icpl.huji.ac.il>
