Resolution criteria to avoid artificial clumping in Lagrangian hydrodynamic simulations with a multi-phase interstellar medium

Large-scale cosmological galaxy formation simulations typically prevent gas in the interstellar medium (ISM) from cooling below $\approx 10^4$ K. This has been motivated by the inability to resolve the Jeans mass in molecular gas (>>$10^5\,\mathrm{M}_{\odot}$) which would result in undesired artificial clumping. We show that the classical Jeans criteria derived for Newtonian gravity are not applicable in the simulated ISM if the spacing of resolution elements representing the dense ISM is below the gravitational force softening length and gravity is therefore softened and not Newtonian. We re-derive the Jeans criteria for softened gravity in Lagrangian codes and use them to analyse gravitational instabilities at and below the hydrodynamical resolution limit for simulations with adaptive and constant gravitational softening lengths. In addition, we define criteria for which a numerical runaway collapse of dense gas clumps can occur caused by over-smoothing of the hydrodynamical properties relative to the gravitational force resolution. This effect is illustrated using simulations of isolated disk galaxies with the smoothed particle hydrodynamics code Swift. We also demonstrate how to avoid the formation of artificial clumps in gas and stars by adjusting the gravitational and hydrodynamical force resolutions.


INTRODUCTION
Numerical simulations of cosmological structure formation are an invaluable tool to test theories of galaxy formation, galaxy evolution and cosmology.These simulations have become increasingly complex in the last decades -from pure N-body simulations (e.g.Springel et al. 2005b) to simulations that model the gas and stars in and around galaxies in great detail (e.g.Horizon-AGN: Dubois et al. 2014, Illustris: Vogelsberger et al. 2014, eagle: Schaye et al. 2015, IllustrisTNG: Pillepich et al. 2018a, Simba: Davé et al. 2019, FIREbox: Feldmann et al. 2023, see also the review from Crain & van de Voort 2023).Gravity is a key ingredient in the Universe on a wide range of scales, from the formation of individual stars to cosmological structure formation, and modelling gravity adequately is critical for a multitude of simulations.
An important approach to validate numerical methods is comparing the gravitational collapse of various objects in simulations to analytical expectations.For gaseous, self-gravitating structures the analytical framework for stability conditions is largely based on Jeans (1902).His analysis of the stability of spherical nebulae to "vibrations" led to various derivations of "Jeans lengths" which generally refer to a critical length scale above which a spherical, self-gravitating gas cloud of a given density and temperature collapses under its own gravity.This scale can be derived by (i) solving the linearized fluid equations and defining a wavenumber  J as the limiting wavenumber ★ E-mail: sylvia.ploeckinger@univie.ac.at between oscillatory solutions ( >  J , short wavelengths) and exponentially growing solutions ( <  J , long wavelengths); (ii) defining the Jeans length as the size of the gas cloud above which the free-fall time is smaller than the sound crossing time; or (iii) comparing the gravitational potential energy  to the internal energy  of a spherical gas cloud.Perturbations can grow for scales where  +  < 0 (for derivations and discussion see e.g.Binney & Tremaine 2008).Each of these derivations leads to a length scale  J above which density perturbations (or: gas clouds) of size , sound speed  s , and (unperturbed) density  are gravitationally unstable: where  is the gravitational constant and  is a dimensionless prefactor of order unity which depends on the exact derivation.Because the assumptions of perfect spherical symmetry and initial zero velocity are rarely fulfilled in real world applications, the Jeans criterion is an approximation and the pre-factors of order unity may vary.Unlike in Eulerian (i.e.grid-based) codes, in Lagrangian (i.e.particle-based) codes, resolution elements can have arbitrarily small separations.If the resolution elements represent an underlying smooth density distribution, gravity is softened below a given length scale in order to avoid close 2-body interactions.A classical prescription for softened gravity is the Plummer (1911) potential  P = −  ( 2 +  2 ) −1/2 for a point mass, , and the gravitational softening length, .The optimal choice of  is already non-trivial for a collisionless N-body system (e.g.Merritt 1996;Romeo 1998;Athanassoula et al. 2000;Dehnen 2001;Power et al. 2003;Rodionov & Sotnikova 2005;Romeo et al. 2008;Ludlow et al. 2019b) and is further complicated in hydrodynamical simulations (e.g.Price & Monaghan 2007;Ludlow et al. 2020).The softening length, , is a measure for the gravitational force resolution of a simulation, because gravity is non-Newtonian (i.e.softened) for  ≲ .
A common technique to solve the equations of motions for a collisional fluid is Smoothed Particle Hydrodynamics (SPH, originally from Lucy 1977;Gingold & Monaghan 1977 but see e.g.Price 2012;Hopkins 2013;Hu et al. 2014;Beck et al. 2016;Borrow et al. 2022 for modern examples) where gas properties, such as the gas density and energy, are smoothed over a kernel which consist of multiple particles, so-called neighbours.The softening length is related to the gravitational force resolution while the kernel size, with the smoothing length, ℎ, as the characteristic length scale, is related to the hydrodynamical force resolution.
The smoothing length, ℎ, generally decreases with increasing gas density, ℎ ∝  −1/3 , for a density-independent number of neighbours per kernel.Depending on the code, the gravitational and hydrodynamical force resolutions can be coupled and  ∝ ℎ ("adaptive softening") or  can be set to a constant value for all densities.Adaptive softening is implemented e.g. in the SPH codes Gadget4 (Springel et al. 2021), Phantom (Price et al. 2018), and the SPH and meshless finite mass / volume code Gizmo (Hopkins 2015).Their adaptive softening lengths are defined such that they are tied to the kernel length (≈ ℎ, the smoothing length of gas particles) or cell size and the hydrodynamic force resolution equals the gravitational force resolution.
In the moving mesh code Arepo (Weinberger et al. 2020), the softening length is adaptive and follows  Arepo =  h (3/(4)) 1/3 , where  h is an input parameter which defines the softening length relative to the cell radius, and  is the volume of the Voronoi cell.In the IllustrisTNG project (Weinberger et al. 2017;Pillepich et al. 2018a), which uses Arepo, the softening length,  is set to a minimum value  min much larger than their minimum cell size (TNG50:  min = 72 pc, minimum cell size: 6.5 pc, Pillepich et al. 2019), an example for a combination of adaptive (for low densities) and constant ( =  min , for high densities) softening lengths.In grid-based codes, such as Ramses (Teyssier 2002), gas resolution elements cannot have separations smaller than the minimum cell size and an additional gravitational softening parameter is unnecessary.
Whether gravitational instabilities are modelled correctly (i.e.follow the Jeans conditions for gravitational instabilities) has been studied mostly for adaptive softening.Bate & Burkert (1997) have defined resolution criteria based on the Jeans length for SPH codes.Their main test case is the isothermal collapse of a one solar mass gas cloud with solid body rotation and they follow both its collapse and fragmentation.Bate & Burkert (1997) advocate using an adaptive softening length equal to the kernel smoothing length ( = ℎ) and ensuring that both resolution parameters can get small enough to resolve the smallest Jeans mass cloud ( J = 4 3 J /3) in the simulation by 2 neigh particles (later reduced to 1.5 neigh by Bate et al. 2003), where  neigh is the number of neighbours in the SPH kernel.If these conditions are not fulfilled, the collapse or lack thereof of self-gravitating structures is determined by the resolution parameters (collapse / fragmentation artificially induced for  < ℎ and artificially inhibited for  > ℎ) and not by the physical conditions.Hubber et al. (2006) introduced the "Jeans test", a numerical test of a plane-wave perturbation within a homogeneous medium, and found that if the resolution criterion from Bate & Burkert (1997) is violated, the collapse of marginally unstable modes is suppressed but artificial fragmentation does not occur in SPH for their setup with an adaptive softening length  = ℎ, confirming the analytical results from Whitworth (1998).Yamamoto et al. (2021) repeated the Jeans test with various hydrodynamic schemes within the Gizmo code (meshless finite mass, meshless finite volume, density-based SPH, and pressure-based SPH, see Hopkins 2015 for details).They confirm the findings of Hubber et al. (2006) for SPH and extend them to other meshless methods, all for adaptive softening with  = ℎ.They found that for none of the hydrodynamic solvers artificial fragmentation is induced, for all solvers the collapse is slowed down if the resolution is lower than required by Bate & Burkert (1997), and the collapse is slowed down more for the meshless methods compared to the SPH methods.
Large-scale simulations of cosmological volumes cannot follow the collapse of individual gas clouds within the interstellar medium (ISM) down to the formation of individual stars and doing so will remain computationally too expensive for the foreseeable future.Utilizing appropriate subgrid prescriptions for star formation and stellar feedback nevertheless allows for realistic galaxy populations (see e.g. the review from Vogelsberger et al. 2020).As an example, the eagle simulation project (Schaye et al. 2015) used an SPH code and a constant gravitational softening length for all baryon particles of  = 700 pc for their (100 Mpc) 3 simulation with an initial baryon particle mass of 1.81×10 6 M ⊙ .This marginally fulfills the Jeans mass criterion set by Bate & Burkert (1997) (but not their recommendation to set  = ℎ) for the warm neutral medium (WNM) with a typical Jeans mass of a few times 10 7 M ⊙ but would violate it drastically for the cold gas phase with typical Jeans masses below 10 4 M ⊙ .An effective pressure floor, sometimes referred to as a polytropic equation of state (EOS),  eff ∝   eff H for gas with a density of  H > 0.1 cm −3 formally solves this issue.For an effective polytropic index of  eff = 4/3, the Jeans mass is constant for gas that is limited by the effective pressure floor (see the discussion in Schaye & Dalla Vecchia 2008).In eagle the normalization for the effective pressure,  eff , is chosen to be at an effective temperature  eff = 8000 K for  H = 0.1 cm −3 .These are typical conditions for the WNM and since the Jeans mass is (marginally) resolved for the WNM, it remains resolved for arbitrarily high densities at their effective pressure.
Therefore, if it is too computationally expensive to resolve the "real" Jeans mass, the "numerical" Jeans mass can be artificially increased.A polytropic EOS was already used for this purpose in the simulations by Bate & Burkert (1997).Richings & Schaye (2016) studied the effect of various pressure floor normalisations in simulations of a dwarf galaxy with a particle mass of 750 M ⊙ and showed that the lowest mass gaseous clumps that formed in their galaxy disks are less compact and less gravitationally bound with a higher artificial pressure floor.An artificially increased Jeans mass in the form of a polytropic EOS, a pressure or entropy floor, or the Springel & Hernquist (2003) subgrid model 1 is implemented in almost 2 all simulations of cosmologically representative volumes that reach  = 0.
The resolution of the gaseous component in a simulation is con-1 Hydrodynamically, the Springel & Hernquist (2003) model acts as a polytropic EOS with a density-dependent  eff , see their figure 1.
2 One recent exception is the FIREbox project (Feldmann et al. 2023).
Table 1.Overview of equations for Newtonian gravity, a Plummer softened gravitational potential, and a Wendland C2 softened potential.
Newtonian (no softening) Plummer softening with softening scale  =  Plummer

Potential
Free-fall time Jeans mass Wendland C2 / Swift softening with softening scale  = 3 =  Plummer and  = / sequently not only defined by one mass and two spatial (gravity and hydrodynamic) resolutions per particle type (i.e.dark matter, gas, and stars) but also affected by a polytropic EOS.For example, a simulation with a theoretically infinite mass and spatial resolution would still not capture the dynamics of the cold gas phase correctly in the presence of a polytropic EOS.An obvious step forward is to remove any artificial pressure floor and allow for a multi-phase ISM (i.e.hot ionized, warm ionized / neutral, cold neutral, as well as molecular gas) to form.Based on the Jeans mass arguments above, this would only be possible for simulations with a baryon mass resolution of ≪ 10 4 M ⊙ .We argue here that the Newtonian Jeans mass does not need to be resolved to avoid numerical clumping and given the success of modern cosmological simulations in reproducing general galaxy properties, it can be questioned how important it really is -in the context of the objectives of large-scale simulations -to model the collapse of individual gas clouds within the ISM correctly.
The aim of this work is to define criteria that help to avoid artificial fragmentation and collapse for simulations with both adaptive and constant values for the gravitational force softening in Lagrangian codes.This is particularly relevant for simulations that model the multi-phase ISM without an artificial effective pressure floor.This paper is organized as follows.In section 2 we introduce the "softened" Jeans criteria (softened Jeans length and softened Jeans mass) and explain why they are more appropriate for describing the conditions of gravitational instabilities in a simulation with softened gravity than the physical Jeans criteria that assume Newtonian gravity.We discuss in section 3 for which densities and temperatures gravitational instabilities are modelled physically correctly, or are numerically induced or suppressed, for both constant and adaptive softening.This section includes an independent criterion on the minimum smoothing length, ℎ min for which a too large value can lead to a numerical runaway collapse (section 3.1).The impact of different choices for the force resolution parameters  and ℎ min is illustrated using simulations of isolated disk galaxies with the modern hydrodynamics code Swift3 (Schaller et al. 2018(Schaller et al. , 2023) ) in section 4. The implications for (cosmological) galaxy formation simulations are discussed in section 5 and we summarize the results in section 6.
In the literature, the terms "smoothing" and "softening" are each sometimes used for either the hydrodynamical or gravitational force calculations involving a "smoothing" / "softening" kernel.Throughout this work we strictly use "smoothing" when referring to the calculation of hydrodynamical quantities and "softening" for calculations of gravitational forces.We use log as log 10 throughout this work.

SOFTENED JEANS CRITERIA FOR A CONSTANT SOFTENING LENGTH
In galaxy formation simulations, mass distributions are discretized by resolution elements such as particles, static grid cells, or moving mesh cells.In order to reduce the artificial 2-body scattering from individual particles that represent a continuous distribution, the gravitational forces are softened for close interactions at a distance  smaller than the softening length  in Lagrangian codes.The Newtonian gravitational acceleration,  N ∝  −2 is replaced by the softened gravitational acceleration, for example  P ∝  ( 2 +  2 ) −3/2 for a Plummer (1911) potential (see Table 1 for an overview).In the limit of large separations,  ≫ , the softened acceleration equals the Newtonian acceleration ( P →  N ∝  −2 ) but for  ≪ , the softened acceleration approaches  P ∝  and therefore diverges from the Newtonian gravity as  → 0. In simulations with a particle mass of  B ≳ 10 5 M ⊙ , gravitational instabilities for densities typical for the cold neutral and molecular gas phases in the ISM are unresolved (examples will be shown and discussed in section 3 and in particular in Fig. 4) and the classical (Newtonian) Jeans criteria provide an inaccurate estimate of the conditions for gravitational instability.
We re-derive the classical Jeans criteria in order to obtain a better description of gravitational instabilities that occur in numerical simulations with softened gravity 4 .We show the derivation for two different softened gravitational potentials: the classical Plummer (1911) potential as well as the modern Wendland (1995) C2 potential (as implemented in Swift) in appendix A. We do not repeat the derivations for other potential shapes, such as the cubic spline kernel (as e.g.used in Gadget4, Springel et al. 2021) but the difference between the different softening potentials is negligible compared to the difference between Newtonian and softened Jeans criteria.
Table 1 summarizes the equations for free-fall times, Jeans lengths and Jeans masses for Newtonian and softened (Plummer and Wendland C2) gravity.The free-fall time in softened gravity (here:  ff,s =  ff,W,fit from equation A34) is longer than the free-fall time in Newtonian gravity for an object of radius  and is proportional to for  ≫ .As example, doubling the softening length therefore slows down the gravitational free-fall in the simulations by a factor of 2.8 for constant initial .Fig. 1 shows an overview of the Jeans masses from equations (A14), (A22), and (A37) for the Newtonian (top panel), Plummer softened (solid line, bottom panel), and the Wendland C2 softened (dotted line, bottom panel) potential, respectively, for a Plummer equivalent softening length of  = 100 pc.At densities above (or at temperatures below) the dashed red line in the bottom panel, the Newtonian Jeans length is smaller than the softening length , and the softened Jeans masses (bottom panel) exceed the Newtonian Jeans mass (top panel).This means that the limit between growing 4 For the derivation of the softened Jeans criteria we assume that the hydrodynamic forces are calculated accurately.The impact of unresolved hydrodynamic forces on numerical instabilities is discussed in detail in section 3. and decaying density perturbations is described by the softened Jeans mass (length) rather than the Newtonian Jeans mass (length).In addition to the different fragmentation scale, the gravitational collapse follows the longer softened free-fall time, rather than the Newtonian free-fall time.
The difference between the Plummer (solid lines) and the Wendland C2 (dotted lines) softened potentials is small considering the approximate nature of the Jeans criterion.In contrast, the softened Jeans mass can be several orders of magnitude larger than the Newtonian Jeans mass for gas densities and temperatures that are typical for the cold ISM.For example, at a gas temperature of a few tens of K and a gas density of 100 cm −3 , the Newtonian Jeans mass is  J,N ≈ 100 M ⊙ while the softened Jeans masses for  = 100 pc is  J,s ≈ 10 5 M ⊙ (and for  = 1000 pc,  J,s ≈ 10 7 M ⊙ ) for both  J,s =  J,P,fit (Plummer) and  J,s =  J,W,fit (Wendland C2).
Throughout this work we consider for simplicity only thermal pressure as a stabilizing force against gravitational collapse in the derivation of the Newtonian Jeans length and mass.For applications where (unresolved) turbulent pressure or magnetic pressure is important, their contributions can be added to the Newtonian Jeans criteria by substituting the sound crossing time in equation (A1) with a more general form of a signal crossing time,  sig , that can also include a turbulent or a magnetic component (see Nobels et al. 2023 for an example of using the turbulent velocity dispersion,  turb , and the velocity dispersion of ions,  Alfven , as additional terms, when deriving the Jeans length).The softened Jeans mass and length are defined relative to the Newtonian Jeans mass and length.For example, .The density at the centre of each distribution is calculated with a Wendland C2 kernel and a smoothing length which is the maximum of ℎ correct (red solid circle) and ℎ min (blue dashed circle) and indicated in the top right of each panel ( H,SPH ).In the top panels the dense gas clump (black circles) consists of 100 particles (≈ number of neighbours in the kernel) and the bottom panels represent a more homogeneous medium without clumping on scales below ℎ min .If the particles clump on scales below ℎ min , then the density that the SPH solver calculates is too low (see top right panel) but if the gas distribution is smooth on scales smaller than ℎ min , then the SPH density estimate is still accurate, even if ℎ correct ≪ ℎ min (see bottom right panel).
for the Wendland C2 kernel (Table 1).The definitions of the softened Jeans mass and length in Table 1 are therefore unaffected by these additional components.
The softened Jeans masses derived in appendix A and listed in Table 1 differ from the Newtonian Jeans mass whenever the particle separations are below the softening length.This might be most prominent in simulations with a (large) constant value for  but is also applicable for simulations with adaptive softening when  is limited by a minimum value  min .

NUMERICAL INSTABILITIES RELATED TO LIMITED RESOLUTION
In large-scale simulations, for example of cosmological volumes, gas can reach gas densities and temperatures that are not formally resolved.Either because gravity or hydrodynamic forces are softened or smoothed on scales larger than the Jeans length or because realistic fragmentation masses are not resolved by enough resolution elements.This is often acceptable, for example if the averaged properties of the ISM within a galaxy are more of interest than correctly following the collapse of individual gas clouds, especially if the collapse of molecular clouds is approximated by subgrid prescriptions, e.g. for star formation.For code stability it is typically preferable to numerically suppress the collapse of gas clouds than to numerically induce collapse because the latter can lead to prohibitively small timesteps, and artificial clumps in the stellar distribution.
In this section we discuss two distinct instabilities that can occur in Lagrangian simulations and how they depend on the gravitational softening length (i.e. the gravitational force resolution), the hydrodynamic smoothing length (i.e. the hydrodynamical force resolution) and the mass resolution.In subsection 3.1 we identify the problematic behaviour of SPH-like simulation codes when the hydrodynamical smoothing length ℎ is limited by a too large minimum value, ℎ min .The resulting instability is caused by an underestimate of the gas density and leads to artificial clumps of closely spaced particles.Setting ℎ min to a very small non-zero value resolves this issue and we cannot identify any disadvantages of letting ℎ min → 0 (to be discussed in detail in section 4).
In subsection 3.2 we analyse gravitational instabilities in simulations with either constant or adaptive softening lengths at and below the hydrodynamical resolution limit, i.e. the size of an individual smoothing kernel.

Numerical instability caused by imposing a minimum hydrodynamical smoothing length
We have shown that depending on the gravitational softening scale, the softened Jeans mass can be orders of magnitude larger than the Newtonian Jeans mass and increases with density at constant temperature for  J,N <  (Fig. 1).Gas clumps would therefore be increasingly stabilised against gravitational collapse as their density gets larger.However, this assumes that the hydrodynamic forces are modelled accurately.For simulations with a minimum smoothing length, ℎ min ≠ 0, the hydrodynamic forces can be over-smoothed for high gas densities and we explain below how this can cause a numerical runaway collapse.
In this section, the smoothing length, ℎ, and its lower limit, ℎ min , are defined as in the hydrodynamics solver Sphenix (Borrow et al. 2022), implemented in the Swift code5 .In Swift, the smoothing length is based on the kernel standard deviation as in Dehnen & Aly (2012) and is calculated as where  res is a constant related to the number of neighbours in the kernel,  H is the hydrogen mass fraction,  B the baryon particle mass,  H the hydrogen particle mass,  H and the hydrogen number density.We use the typical values  res = 1.2348 (this corresponds to 58 neighbours with a Wendland C2 kernel) and  H = 0.74.
The minimum smoothing length in Swift is set by the dimensionless parameter ℎ min,ratio which is defined as the ratio between the kernel support, ℎ k , and the length scale above which gravity is Newtonian,  = 3.For the Wendland C2 kernel  k = 1.936492 (Dehnen & Aly 2012) and the minimum smoothing length is therefore ℎ min = 1.55  ℎ min,ratio .The density above which the hydrodynamic forces are over-smoothed is Note that the critical density  H,hmin can remain constant for simulations with different mass resolutions  B if ℎ min is adjusted accordingly.
The softened Jeans criteria derived in appendix A are therefore only valid for  H ≤  H,hmin .For higher densities the hydrodynamic forces can become inaccurate due to over-smoothing and this inaccuracy depends on the density contrast within ℎ min .The top panels of Fig. 2 demonstrate this for an extreme case with one clump of dense gas (represented by 100 resolution elements: circles in Fig. 2).The 3D particle positions are random but restricted to a volume that corresponds to a set input density  H which increases from log  H [cm −3 ] = 2 (leftmost panel) to log  H [cm −3 ] = 6 (rightmost panel).The red solid circle represents the smoothing length for this gas density and particle mass (equation 7) and the blue dashed circle represents the minimum smoothing length which is set to ℎ min = 15.5 pc (corresponding to e.g. = 100 pc and ℎ min,ratio = 0.1) in this example.
Each panel in Fig. 2 lists the input density as well as the density  H,SPH that the SPH solver calculates at the centre of each box using a Wendland C2 kernel.If ℎ becomes smaller than ℎ min and the gas is clumpy on scales smaller than ℎ min , then the estimated density (and pressure) is increasingly underestimated.This can lead to a runaway collapse if the underestimated pressure does not balance the (softened) gravitational forces.The artificial, unresolved collapse is difficult to spot in the gas properties because the SPH gas density does not correctly represent the particle distribution, but it can lead to the formation of clusters of stellar particles that are denser than expected based on the resolution of the simulation and the gas densities.The bottom panels of Fig. 2 show a similar situation, but here the dense gas is homogeneous on scales below ℎ min , in which case the output densities  H,SPH remain accurate, even for ℎ ≪ ℎ min (see also figure 1 in Borrow et al. 2021 for an illustration).
To avoid a potential run-away collapse, which is desirable for both numerical and physical reasons, we can define a problematic region in the density-temperature phase-space based on the conditions discussed above: i ℎ < ℎ min (or:  H >  H,hmin ) ii fragmentation on scales below ℎ min For an estimate of the clumping that is expected in the simulation, we use the softened Jeans length and (ii) then translates to  J,s < ℎ min .
Assuming that ℎ min <  (to avoid inducing numerical collapse, see e.g.Bate & Burkert 1997), the condition  J,s < ℎ min falls into the softened part of the phase-space where  ≡ 3 ≫  J,N and we therefore approximate the softened Jeans mass (here for the Wendland C2 kernel, equation A36) with  J,s ≈ 0.27 0.3  0.4 J,N  0.6 and the conditions (i) and (ii) for a potential run-away collapse correspond to the following regions in the phase-space diagram as This "runaway collapse zone" is indicated in Fig. 3 for various combinations of the particle mass  B , Plummer equivalent softening length , and the minimum smoothing length ℎ min .The default values (thick, black lines) with log  B [M ⊙ ] = 6.3,  = 700 pc, and ℎ min = 108.5 pc (i.e.ℎ min,ratio = 0.1) are those from the eagle (100 Mpc) 3 simulation and two of these parameters are constant in each panel while the third parameter is varied (top panel:  B , middle panel: , bottom panel: ℎ min ).
Each zone is limited through a combination of a vertical line, condition (i):  H =  H,hmin (equation 9) and a line of log  ∝ log  H , i.e. for a constant softened Jeans length of  J,s = ℎ min (condition ii, equation 10).Densities with higher values than the indicated lines are reported incorrectly by the SPH kernel density estimator if the gas is clumpy on scales below ℎ min and can lead to a runaway collapse.
Varying the baryon particle mass  B at constant  and ℎ min,ratio (Fig. 3, top panel) only changes  H,hmin but condition (ii) remains unaffected.A better mass resolution is even slightly counter-productive in resolving higher densities because  H,hmin ∝  B for constant ℎ min (equation 8).Similarly, reducing the gravitational softening length decreases the softened Jeans length and the runaway collapse zone includes higher temperatures (middle panel of Fig. 3).Decreasing the minimum smoothing length at constant  B and constant  on the other hand pushes the runaway collapse zone to significantly higher densities (from 10 2 cm −3 to > 106 cm −3 when reducing ℎ min by a factor of 10, bottom panel of Fig. 3).
If the highest (expected) density in a simulation is  sim,max , the minimum smoothing length should be set so that ℎ( sim,max ) ≥ ℎ min .With ℎ from equation ( 7), this means that the minimum smoothing length should follow for a criterion based only on the vertical lines in Fig. 3.As we will discuss in section 5, we cannot identify any benefit for larger values for ℎ min and therefore recommend to use a very small but non-zero 6 value for ℎ min that generously fulfills the conservative criterion from equation ( 11).Equation 11does not depend on  and applies to both simulations with adaptive and constant softening lengths.

Examples from the literature
The runaway collapse zone is indicated in Fig. 4 for three simulation projects that use the same code (a modified version of gadget3, Springel 2005) but span more than 2 orders of magnitude in mass resolution with baryon particle masses of  B = 1.81 × 10 6 M ⊙ (eagle),  B = 2.26 × 10 5 M ⊙ (eagle-high-res), and  B = 10 4 M ⊙ (apostle).
The resolution parameters (see Table 2) used in both eagle flagship runs (eagle, eagle-high-res, Schaye et al. 2015) and the highest resolution zoom-in simulations from the apostle simulation suite (Fattahi et al. 2016;Sawala et al. 2016) violate the condition from equation ( 11) for densities of  sim,max ≳ 100 cm −3 .However, for these simulations, the numerically induced runaway collapse was Table 2. Values for the gas particle mass ( B ), the gravitational force softening length ( ), and the minimum hydrodynamical smoothing length (ℎ min ) for a few simulations with constant softening from the literature and an alternative set of resolution parameter proposed in this work as shown in Fig. 4. avoided because gas with densities above 0.1 cm −3 follows a polytropic equation of state where the effective temperature7 is with the polytropic index  eff = 4/3.This effective temperature is indicated for reference as thick grey line in Figs. 3 and 4 and is identical for eagle and apostle.
The apostle zoom-in simulations have a particle mass of  B = 10 4 M ⊙ but the runaway collapse zone is barely affected by the higher mass resolution, compared to e.g.eagle.The effective temperature floor was therefore necessary for their values for ℎ min , despite the low particle mass.
The softened Jeans mass, as derived in appendix A for a Wendland C2 kernel and for the eagle and eagle-high-res softening parameters, is overlaid as contours in units of log  J,s [M ⊙ ] in Fig. 4. Artificial fragmentation that is caused by not resolving the Jeans mass with enough resolution elements (as in Bate & Burkert 1997) would not be expected as the softened Jeans mass increases with density (at a constant temperature) as soon as gravity is no longer Newtonian.While we argue that numerical runaway collapse would occur in these simulations without an effective temperature floor, the reason for this is not related to the Newtonian Jeans mass but to the adoption of a too large minimum smoothing length as discussed in section 3.1 and therefore not very sensitive to the particle mass.
The bottom panel of Fig. 4 illustrates an alternative set of resolution parameters for any mass resolution between those of eagle and apostle-L1 (the dependence on  B is very weak, see top panel of Fig. 3).A Plummer equivalent softening of  = 200 pc is small enough to model gravitational instabilities correctly in the WNM but still large enough for a minimum softened Jeans mass of ≈ 10 6 M ⊙ in all neutral phases of the ISM.Counter-intuitively, a smaller softening length would decrease the softened Jeans mass in the cold neutral medium (CNM) and molecular clouds (MCs), which means that the softened Jeans mass would be resolved by fewer particles in these phases.With a minimum smoothing length of ℎ min = 2 pc (instead of 54.25 pc in eagle-high-res), the runaway collapse zone only covers  H > 10 8 cm −3 and the molecular phase can be directly modelled without suffering from the numerical issues described above.Simulations of isolated galaxies at comparable mass resolutions and without an entropy floor are shown below in section 4.
Future simulations that include a multi-phase ISM need to carefully select the combination of mass (baryon particle mass), gravity  2).The contours indicate the softened Jeans mass in units of log  J,s [M ⊙ ].The red shaded region at high densities is the runaway collapse zone and SPH-estimated gas densities within this zone may be underestimated (see section 3.1).The grey thick line is the effective temperature floor used in eagle, eagle-high-res, and Apostle-L1.Typical densities and temperatures for the warm neutral medium, the cold neutral medium and molecular clouds (from low to high densities) are indicated with black patches for reference.An alternative set of parameters is shown in the bottom panel: for a mass resolution of order 10 5 to 10 6 M ⊙ , a softening scale of 200 pc and a minimum smoothing length of 2 pc (ℎ min,ratio = 0.006) push the runaway collapse zone to densities ≳ 10 8 cm −3 ).All phases in the ISM can therefore be modelled without triggering a numerical runaway collapse.
(softening scale), and hydrodynamic resolutions (minimum smoothing length) to avoid unwanted numerical artefacts.

Gravitational instabilities at and below the resolution limits
In subsection 3.1 we defined a numerical instability that is caused by an incorrect density estimate if the smoothing length is limited by a constant minimum value.In this subsection, we focus on the formation of gas clumps through gravitational fragmentation and collapse, considering both codes with adaptive and constant constant values for the gravitational softening length.
For a code-independent discussion of the conditions under which gravitational instabilities develop in simulations, we define the codeagnostic length scales  soft , over which gravitational forces are softened, and  smooth (with a minimum of  smooth,min ), over which the hydrodynamical forces are smoothed.We explain briefly how they connect to the smoothing lengths ℎ and the softening lengths  in some codes as examples: We use the kernel support as measure of  smooth .In Swift, ℎ is defined following Dehnen & Aly (2012) and based on the kernel standard deviation which is directly related to the numerical resolution of sound waves (equation 7).For the Wendland C2 kernel and the Swift definition of ℎ,  smooth =  k ℎ = 1.936492 ℎ.In practice, the definition of the smoothing length varies for each code.The SPH code Gadget4 (Springel et al. 2021) defines ℎ as the finite support of the kernel and therefore  smooth would be equal to ℎ Gadget4 .
The softening length is typically given as the Plummer-equivalent softening length, .Depending on the softening potential, the gravitational force is exactly Newtonian for particle separations of ≥ 2.8 (cubic spline kernel as in Gadget4, Springel et al. 2021) or ≥ 3 (Wendland C2 kernel as in Swift).Because the transition to softened gravity is very gradual, the deviations from Newtonian gravity are only significant for separations ≲ 1 − 1.5 (see e.g.figure 1 in Springel et al. 2021).We therefore use  soft = 1.5 as an estimate for the softening length scale.
For the following discussion, using a different measure for either  soft or  smooth would shift the lines slightly but does not change the general interpretation of the individual zones in temperature-density space.After concluding in section 3.1 that the minimum smoothing length should be very small to avoid a runaway collapse, we assume in this subsection  smooth,min → 0, if not explicitly mentioned otherwise.

Gravitational instabilities at the hydro resolution limit 𝑙 smooth
Gravitational instabilities cannot be modelled accurately in the simulation if the Newtonian Jeans length,  J,N , is smaller than the size of an individual smoothing kernel,  smooth .The hydrodynamical forces are inaccurate on scales below  smooth and fragmentation on scales of  J,N <  smooth is unresolved.For the Newtonian Jeans length, the condition  J,N =  smooth depends on the gas density and temperature as well as on the particle mass,  B ( smooth ∝ ), but it does not depend on the gravitational softening length,  soft .The condition  J,N =  smooth therefore applies in the same way for simulations with constant and adaptive softening lengths.In each case, gravitational instabilities are only resolved (i.e.modelled accurately) for  J,N >  smooth .
For gravitational instabilities at the hydrodynamic resolution limit,  smooth , differences between the formation and further evolution of gravitational instabilities emerge that depend on the assumed gravitational softening lengths.In this subsection we discuss the formation of gravitational instabilities at the hydrodynamical resolutions limit,  smooth , for simulations with constant and adaptive softening lengths.
Figure 5. Gravitational stability at the hydrodynamical resolution limit,  smooth , for simulations with a constant softening length  soft .Perturbations at length scale,  smooth , are expected to grow in Newtonian gravity for  J,N <  smooth ("gravitationally unstable in nature") and to decay for  J,N >  smooth ("gravitationally stable in nature"), separated by the thick dashed line with a slope of 1/3 ( ∝  1/3 H ) in each panel for  smooth,min = 0. Here, the gas density,  H , refers to the SPH-estimated density.In simulations with softened gravity, such perturbations are expected to grow for  J,s <  smooth ("gravitationally unstable in sim") and to decay for  J,s >  smooth ("gravitationally stable in sim"), separated by the solid line which diverges from the dashed line and follows a slope of -2/3 ( ∝   .Gravitational stability at the hydrodynamical resolution limit, max( smooth ,  smooth,min ), for simulations with an adaptive softening length,  soft .Line styles as in the left panel of Fig. 5 but for an adaptive softening length,  soft .In the left panel, the minimum value for the gravitational softening length equals the minimum value for the smoothing length.Here, the slope of both lines (dashed line:  J,N = max( smooth ,  smooth,min ), solid line:  J,s = max( smooth ,  smooth,min )) changes when  smooth is limited by  smooth,min .In the right panel, the softening length is adaptive down to a minimum value  soft,min , for which the softening length is effectively constant.In this panel, the slope of the solid line,  J,s = max( smooth ,  smooth,min ), changes for densities above which the softening length,  soft , is limited by a minimum value,  soft,min .Values for  B ,  soft,min , and  smooth,min were selected to represent FIREbox (left panel) and TNG50 (right panel, see text for details).As in Fig. 5, the gas density,  H , refers to the SPH-estimated density.Typical densities and temperatures for the WNM, CNM, and MCs are indicated with dark patches, as in Fig. 4.

Constant softening length:
In nature (or in simulations with  B → 0,  soft → 0, and  smooth,min → 0), density perturbations with a length scale of  smooth grow if  J,N <  smooth , and decay until they reach a new equilibrium for  J,N >  smooth .The border 8 (i.e.case  J,N =  smooth ) is indicated as a thick dashed line in all panels in Fig. 5.As discussed above, gravitational instabilities are only modelled accurately (i.e. are resolved) if  J,N >  smooth (unshaded area).
For softened gravity with a constant softening length (left and right panel:  soft = 100 pc, middle panel: various values for  soft , see the labels), the solid lines ( J,s =  smooth ) separate perturbations with length scales  smooth that are gravitationally (un)stable in softened gravity, i.e. in the simulation.Densities and temperatures for which gas is gravitationally stable in nature for perturbations at the resolution limit  smooth are also gravitationally stable in simulations with softened gravity (white area in the left panel of Fig. 5 for a simulation with  B = 10 5 M ⊙ ,  soft = 100 pc, and  smooth,min → 0 pc).
Gas with densities and temperatures within the shaded areas ( J,N <  smooth ), is gravitationally unstable for perturbations at the size of a smoothing kernel in Newtonian gravity.In softened gravity, gravitational instabilities in part of this area are suppressed at scales of  smooth because the softened Jeans length exceeds the kernel size ( J,s >  smooth ) and therefore perturbations on scales of  smooth decay.The middle panel shows that the boundary  J,s =  smooth depends on the value for the constant softening length,  soft .The right panel shows that increasing (dash-dotted lines) or decreasing (dotted lines) the particle mass by a factor of 8 from the fiducial value of  B = 10 5 M ⊙ (dashed lines), affects both the thick green lines for  J,N =  smooth as well as the thin black lines for the condition B (right panel of Fig. 5 for  soft = 100 pc).

Adaptive softening length:
For adaptive softening, typically  soft =  smooth and therefore  J,s ≈  J,N .The instability criteria for perturbations at the resolution limit, i.e. at the kernel size  smooth , therefore follow the instability criteria for Newtonian gravity.The left panel of Fig. 6 (linestyles as in Fig. 5) illustrates this for parameters representative for the FIREbox (Feldmann et al. 2023) simulation project with  B = 6 × 10 4 M ⊙ .Their adaptive softening length is equal to the average gas particle separation ( soft = 1.5 ( B /) 1/3 , for gas density ), down to a minimum value of  soft,min =  smooth,min = 2.25 pc ( min = 1.5 pc).Density perturbations of length-scale  smooth therefore follow the Newtonian Jeans criteria.The small offset between the dashed and solid lines is from the order of unity prefactors related to the kernel shapes and exact definitions of  smooth and  soft .For simplicity we use the same kernel shapes and pre-factors as in Fig. 5.
The slope change of the  J,N =  smooth (dashed) and  J,s =  smooth (solid) lines is caused by the minimum smoothing length  smooth,min .The critical density  H,hmin above which densities are over-smoothed, and hence  H,SPH <  H , is ≈ 2 × 10 6 cm −3 (equation 8 for  B = 6 × 10 4 M ⊙ and ℎ min ≈ 1.2 pc).The small values for  smooth,min and  soft,min therefore prevent the instability described in section 3.1 for densities below ≈ 2 × 10 6 cm −3 .
In the IllustrisTNG project (Pillepich et al. 2018b; Nelson et al.   8 The condition  J,N =  smooth can also be interpreted as the condition that the Jeans mass,  J,N = 4  3 J,N /3, is equal to the mass in the kernel  neigh  B ≈  3 smooth  (with a constant pre-factor of order unity).However, the number of neighbours,  neigh , is not well defined and we therefore use the better defined condition  J,N =  smooth instead.
2018) which uses the moving mesh code Arepo (Weinberger et al. 2020), the softening length is related to the adaptive sizes of the gaseous cells but here the minimum softening length (TNG50:  soft,min = 1.5 gas,min = 108 pc, TNG100:  soft,min = 1.5 gas,min = 285 pc, TNG300:  soft,min = 1.5 gas,min = 555 pc, see table 1 in Pillepich et al. 2019) differs from the minimum smoothing length (the minimum cell size in TNG50 is reported as 6.5 pc, Pillepich et al. 2019).
Gravitational instabilities in gas cells for which the gravitational softening is limited by a minimum value ( =  gas,min ) follow the softened Jeans criteria outlined in section 2. The right panel of Fig. 6 shows the  J,N =  smooth (dashed) and  J,s =  smooth (solid) lines for values representative for TNG50:  B = 8.5 × 10 4 M ⊙ ,  soft,min = 108 pc and an assumed negligible value for  smooth,min .
Comparing the right panel of Fig. 6 with the left panel of Fig. 5, we see that gravitational instabilities on scales <  soft,min in simulations with adaptive softening lengths effectively behave as in simulations with constant softening lengths.

Gravitational instabilities below the hydro resolution limit 𝑙 smooth
Density perturbations grow in nature on length-scales smaller than  smooth in gas with temperatures and densities for which  J,N <  smooth (shaded regions in Figs. 5 and 6).The gravitational collapse of these perturbations within an individual smoothing kernel is not resolved and therefore not expected to be modelled accurately with any simulation method.We can see from Figs. 5 and 6 that a large fraction of the neutral ISM in simulations might form sub-kernel instabilities ( J,s <  smooth , dark shaded regions).While the resulting clumpy particle configurations within a kernel might have a limited effect on the simulation because the cooling rates, star formation rates and other density-or pressure-dependent sub-grid models use the smoother SPH estimates, we briefly discuss the differences in the treatment of sub-kernel perturbations between codes with constant and adaptive softening lengths.

Constant softening length:
Gas with  J,s <  smooth is expected to be gravitationally unstable when exposed to fluctuations on length scales between  J,s and  smooth .In the derivation of the softened Jeans length,  J,s , we assume an accurate gas density and pressure estimate.If the gas pressure is underestimated, gravitational instabilities may form from perturbations on length scales below  J,s .Therefore, the softened Jeans criteria serve as an upper limit to the expected instabilities in the simulation because the hydrodynamic forces might be underestimated due to the smoothing over the full kernel (an extreme case is discussed in subsection 3.1).If  smooth <  soft , the gravitational forces are softened on larger scales than a smoothing kernel and further gravitational collapse and fragmentation within a kernel is suppressed.On the other hand,  soft can be (much) smaller than the size of a smoothing kernel, for a constant value of  soft .In this case, dense particle configurations with sizes possibly even smaller than  J,N can form because the density and pressure estimates of sub-kernel clumps can be inaccurate.Fig. 7 shows an overview of the behaviour of dense particle configurations within a smoothing kernel.Gas that is gravitationally unstable at the scale of an individual smoothing kernel ( J,s <  smooth , as in Fig. 5) is further split into densities for which further fragmentation to even smaller scales is suppressed ( soft >  smooth ) or potentially induced ( soft <  smooth ), with the boundary,  soft =  smooth , indi- We conclude that for a constant softening length, gas in simulations with  J,s <  smooth and  soft <  smooth (purple areas in Fig. 7) might be affected by sub-kernel clumping.For example: within a kernel of 100 particles, a subset of 10 particles has particle separations that are smaller than average for this kernel and hence smaller than the smoothing length,  smooth .In contrast to simulations with adaptive softening, the constant value for  soft can be smaller than the particle separation of these 10 particles.The gravitational forces between this subset of particles are therefore Newtonian, while the pressure is smoothed over the full kernel with 100 particles.An initial perturbation within the kernel may grow artificially or at a rate that is artificially high until  smooth =  soft (vertical lines in Fig. 7).For higher gas densities (i.e. smooth <  soft ), the further collapse is suppressed.The detailed impact of sub-kernel clumping on galaxy properties in large-scale simulations is beyond the scope of this work.
If sub-kernel clumping ("subk") is undesired in simulations with constant softening, the conditions  J,s <  smooth and  soft <  smooth (purple areas in Fig. 7) should not cover regions in densitytemperature space that are populated by many gas particles in the simulation.If we want to confine the area in density temperature space for which  soft <  smooth to densities below  max,subk , this condition translates into a minimum constant softening length for the smoothing length and softening length definition as above.
Unlike in simulations with adaptive softening, unresolved gravitational instabilities are not suppressed by design in simulations with constant softening lengths.Fulfilling equation ( 13) avoids potentially undesired, sub-kernel gravitational instabilities in large parts of the cold ISM, but at the expense of suppressing physical instabilities at the kernel scale,  smooth (see middle panel of Fig. 5).
We focus in this work on simulations with  B ≳ 10 5 M ⊙ but the analysis presented here is relevant for simulations of any mass resolution.In Appendix B (Fig. B1), the right panel of Fig. 7 is repeated for simulations with a particle mass of  B = 4 M ⊙ , for which a constant softening length of at least  = 2 pc is needed to avoid sub-kernel clumping.

Adaptive softening length:
In simulations with an adaptive softening length and  soft ≥  smooth any further fragmentation within a smoothing kernel is generally suppressed because gravitational forces are softened on scales larger than or equal to the smoothing kernel.Physical gravitational instabilities on scales  J,N <  smooth are therefore artificially suppressed.
In order to avoid the runaway collapse described in section 3.1, the minimum smoothing length,  smooth,min , needs to be small (see equation 11).In simulations, such as FIREbox, where  soft,min =  smooth,min , the minimum softening length needs to have the same small value as the minimum smoothing length.

GALAXY SIMULATIONS
In section 3.2, we used the smoothing and softening lengths  smooth and  soft for a code independent discussion on the individual zones.In this section we use simulations of isolated galaxies with the public Swift code9 (Schaller et al. 2016(Schaller et al. , 2018(Schaller et al. , 2023, www.swiftsim.com), www.swiftsim.com).The resolution parameters set by the user are  and ℎ min which relate to  soft = 1.5 and  smooth,min = 1.94ℎ min , respectively (see text in section 3.2 for details).We will demonstrate the numerical issues that can occur for different choices of  and ℎ min , focusing on simulations of isolated galaxies with a constant softening length.The with identical initial conditions.The densities (x-axis) in the top row are the SPH gas density estimates,  H,SPH , as used in the simulations and stored in the snapshots, while the bottom row shows the re-calculated densities for ℎ min → 0, based on the same particle positions.The colour of the temperature-density histograms is proportional to the log of gas mass per pixel and the minimum and maximum values for the colour maps are constant across all plots.From left to right the minimum smoothing length, ℎ min , decreases from 77.5 pc (ℎ min,ratio = 0.2; first column) and ℎ min = 24.5 pc (second column) to ℎ min = 7.75 pc (third column).The red shaded area is the runaway collapse zone as defined in section 3.1 and each panel includes a number for the total mass of particles within this zone as log  zoneC [M ⊙ ] if  zoneC > 0. For reference, the black lines indicate where the softened Jeans mass is resolved by 1 (dotted), 10 (dashed), and 100 (solid) particles.The panels in the rightmost column show the cumulative mass fraction for particles above a given gas density  H .The solid lines in the top panel show the density distributions for the phase-space diagrams in the top row and the dashed lines in the bottom panel for the density distribution from the bottom row (dashed lines are repeated for reference in the top panel).The short vertical dotted lines at the top of the figures in the rightmost row indicate the densities  H,hmin (equation 8) above which the smoothing length is limited by ℎ min .
individual simulations are based on the "IsolatedGalaxy-feedback" example in Swift and use the modern and open source SPH scheme SPHENIX (Borrow et al. 2022), implemented as the default hydrodynamic solver in Swift.Gravity is softened with a Wendland C2 kernel and a constant softening length, defined as the Plummer equivalent softening length, .
The initial conditions were created with MakeNewDisk which is based on the code used in Springel et al. (2005a) but modified to use the exact definition of the analytic dark matter halo mass  200 (i.e.mass within  200 , the radius within which the average density is 200 times the critical density of the Universe) instead of the halo mass integrated to infinity (see Nobels et al. 2023 for details).The isolated galaxy initially has a gas mass of 1.644×10 10 M ⊙ with solar metallicity ( ⊙ = 0.0134, Asplund et al. 2009) and an exponential stellar disk with a radial scale length of 4.3 kpc and a mass of 3.836 × 10 10 M ⊙ The initial disk gas fraction is 30 per cent and the gas is initialized with a temperature of 10 4 K.The baryonic disk is in equilibrium with an analytic dark matter halo with a Hernquist (1990) profile of  200 = 1.37 × 10 12 M ⊙ and a scale radius such that the central density profile matches that of a Navarro et al. (1997) profile with a concentration of  = 9.The initial conditions are available within the "IsolatedGalaxy" example in Swift.
We use the fiducial cooling tables from Ploeckinger & Schaye (2020) which include the effects of self-shielding, dust, cosmic rays, an interstellar radiation field and a UV background.In Appendix C2 we demonstrate that our conclusions remain unchanged if we instead use cooling tables appropriate for a weaker and stronger radiation field.No artificial pressure or entropy floor is included.A supernova energy of 10 51 erg per SN is injected stochastically as thermal energy, following Dalla Vecchia & Schaye (2012), with a heating temperature of 10 7.5 K.The high heating temperature of the stochastic feedback model increases the efficiency of thermal feedback because the cooling time of gas with 10 7.5 K is long and less energy is lost radiatively lost compared to using many smaller thermal energy injections that would heat up the gas to e.g. 10 5 K, the peak of the cooling curve (see Dalla Vecchia & Schaye 2012 for a detailed discussion).We inject the stellar feedback energy into the gas particle closest to the star, which has been shown to further increase the efficiency of supernova feedback, compared to selecting a random gas particle in the star's kernel (see "Min distance" model in Chaikin et al. 2022).Additional simulations with 2 and 4 times higher supernova energies are presented in Appendix C1.
Star formation is limited to densities of  H > 0.1 cm −3 and cold gas (temperatures of  < 1000 K) and the star formation rate for each gas particle with mass  gas is given by the Schmidt (1959) relation with the star formation efficiency  sf and the Newtonian free-fall time  ff,N (equation A10).A gas particle is converted into a star particle stochastically.The simulations analysed in this work vary the resolution parameters: gravitational softening  = [250, 500] pc, minimum SPH smoothing length ℎ min = [7.75, 24.5, 77.5] pc (ℎ min,ratio = [0.02,0.063, 0.2] for  = 250 pc), and baryon particle mass  B = [10 5 , 8 × 10 5 ] M ⊙ ; as well as the star formation efficiency  sf = [0.32,1] per cent in order to change the amount of high-density gas in a controlled way.
All simulations run until  end = 1.2 Gyr.We use Swift to recalculate the gas densities based on the gas particle positions for ℎ min → 0 by restarting the original simulations at  = 1 Gyr for one very small timestep (Δ max = 10 yr) with a very small minimum smoothing length10 (ℎ min,ratio = 10 −4 ).The snapshot file that is produced after this timestep contains the correct (i.e.not limited by a minimum smoothing length) gas densities based on the gas particle positions.
Fig. 8 shows the difference between the original SPH densities (top row) and the recalculated densities for ℎ min → 0 (bottom row) for 3 simulations (first three columns; all with  B = 10 5 M ⊙ ,  = 250 pc,  sf = 0.32 per cent) with decreasing values for ℎ min (from left to right).The runaway collapse described in section 3.1 is obvious in the simulation with the largest value for ℎ min (1st column).The distribution of recalculated, "real" densities (bottom row) extends to much higher values than the densities used during the simulation (top row) for many particles as they approach the runaway collapse zone (red shaded area).For the simulation with the smallest value of ℎ min (3rd column), the two density-temperature diagrams look identical.
The panels in the rightmost column show the cumulative mass fraction of particles above density  H (x-axis) for the SPH densities (solid lines, top panel) and the recalculated densities (dashed lines, bottom panel; repeated for reference in the top panel).As the densities exceed  H,hmin (equation 8; indicated as short vertical dotted lines) the SPH densities and the recalculated densities diverge.In this case, the formally lowest resolution simulation (ℎ min = 77.5 pc) has much higher "real" (but unresolved) densities, up to  H > 10 6 cm −3 , than the formally highest resolution simulation (ℎ min = 7.75 pc).If these dense clumps would appear through physical processes such as instabilities in the galaxy disk, they would persist in the simulation with better resolved hydrodynamic forces.We therefore conclude that the very high densities seen in the particle distributions are indeed artificial as outlined in section 3.1.
While the dense gas clumps form, they may become eligible to star formation and if their densities are underestimated, their star formation rates are underestimated as well (  ★ ∝  1/2 H , equation 14).As result, the gas clumps in the runaway collapse zone turn into star particles slower than expected from the value for the star formation efficiency and the particle positions (i.e."real densities").
After star particles form within dense gas clumps, stellar feedback is modelled by injecting thermal energy.Dalla Vecchia & Schaye (2012) derived a maximal gas density below which thermal feedback is efficient, which for the simulations used in this work means that thermal feedback is inefficient for densities above  fb ≈ 10 cm −3 .Their derivation assumes ℎ > ℎ min which is fulfilled at densities of ≈ 10 cm −3 in all simulations in this work.Comparing the values of  fb to the particle densities in Figs. 8 and 9 reveals that the artificially dense gas clumps are unlikely to be destroyed by stellar feedback and may thus result in artificially dense clusters of star particles.
In Fig. 9 we compare images of the stellar mass surface density of five simulations and their respective gas distributions in temperature -recalculated density phase-space (all at  = 1 Gyr) The lines and labels are as in Fig. 8.The central images for both gas and stars are for the fiducial simulation with resolution parameters  B = 10 5 M ⊙ ,  = 250 pc, and ℎ min = 77.5 pc and a star formation efficiency of  sf = 0.32 per cent.The four simulations in each corner vary one of these parameters at a time, keeping the other parameters at their fiducial values (see labels).
Some stellar mass surface density images (right side of Fig. 9) show dense star clumps.We identify the clump properties by running the Swift friends-of-friends algorithm as a stand-alone routine on the star particles within a snapshot output.We use a small linking length of 25 pc and a minimum number of particles of 50 for simulations with  B = 10 5 M ⊙ and 18 for simulations with  B = 8 × 10 5 M ⊙ .The minimum clump mass is a factor of ≈ 2.9 times larger for the low mass-resolution simulation compared to the simulation with the fiducial particle mass.This is a compromise between choosing the same number of particles and the same mass in a clump when comparing simulations with different mass resolutions.Slightly different parameters in the clump finding algorithm might affect the detailed analysis on the exact properties of the individual clumps but the qualitative conclusions that we focus on are insensitive to these choices.The positions of all clumps identified outside of the innermost 3 kpc are indicated by circles and the most massive clumps ( clump > 10 8 M ⊙ ) are highlighted by diamonds.
The galaxy with the fiducial resolution parameters (centre) has several stellar clumps due to the large amount of mass (log  (zoneC) [M ⊙ ] = 9.37) within the runaway collapse zone (zone C in section 3.2) in the gas phase.These clumps are mostly artificial as they largely disappear when increasing the hydrodynamic force resolution (shown in Fig. 8 and when comparing centre and lower left panels in Fig. 9) and therefore moving the runaway collapse zone to higher gas densities.
The number of clumps as well as their masses are not very sensitive to the mass resolution (increase  B ; lower right).On the other hand, increasing the gravitational force softening length (increase ; upper right) or increasing the star formation efficiency (increase  sf ; upper left) both reduce the number of stellar clumps drastically by reducing the amount of high-density gas and therefore the amount of gas within the runaway collapse zone.Interestingly, their gas distributions look very different11 : while the simulation with the higher star formation efficiency ( sf = 1 per cent) has very little cold (< 1000 K) gas, the galaxy with the larger softening length ( = 500 pc) has large amounts of cold gas at densities just outside the runaway collapse zone but its further collapse is delayed due to the larger softening scale.A tail towards very high densities it still noticeable but in this case it is limited to the central part of the galaxy and therefore does not result in clumps throughout the galactic disk.
In contrast to the comparisons shown in Figs. 8 and 9 for which only one parameter is varied at a time, we show in Fig. 10 an alternative combination of values for  and ℎ min .For a particle mass of  B = 10 5 M ⊙ (first two columns in Fig. 10) the alternative parameter set fulfills the conditions in equation ( 11) for ℎ min and equation ( 13) for .Note that these parameters are not necessarily the best choice for every application but illustrate the impact of choosing values for  and ℎ min that are informed by the conditions defined in sections 3.1 and 3.2.
Both simulations within each figure half in Fig. 10 use the same number of particles (i.e.baryon particle mass,  B = 10 5 M ⊙ , left two panels, and  B = 8 × 10 5 M ⊙ , right two panels) and all simulations use the same value for the star formation efficiency ( sf = 0.32 per cent).The gravitational softening length is increased from 250 pc (fiducial) to 500 pc (alternative).The minimum smooth-Figure 10.Stellar mass surface density (top row) and gas phase-space distribution (bottom row) as in Fig. 9 for the fiducial simulation parameters (first and third columns) and a set of alternative simulation parameters (as in Table 2 and Fig. 4; second and fourth columns) for a particle mass of  B = 10 5 M ⊙ (left two columns) and  B = 8 × 10 5 M ⊙ (right two columns), and a star formation efficiency of  sf = 0.32 per cent.The combination of a larger gravitational softening scale ( = 500 pc) and the smaller minimum smoothing length (ℎ min = 7.75 pc) result in a galaxy without artificial stellar clumps and without gas in the runaway collapse zone.The low star formation efficiency still results in a large amount of cold gas with  < 100 K but at better resolved gas densities.
ing length is reduced from ℎ min = 77.5 pc (fiducial) to ℎ min = 7.75 pc (alternative).These changes move the runaway collapse zone to very high densities.In addition, the zone of induced collapse (zone D in section 3.2, red shaded area at low densities in Fig. 10) is also slightly reduced.The slower gravitational collapse ensures that the majority of the cold gas is at more manageable densities (< 10 2 cm −3 ) and the density distribution does not extend into the runaway collapse zone.For both mass resolutions, the stellar mass surface density is smoother for the alternative resolution parameters and there are no dense star particle clumps, in contrast to the massive clumps present in the simulations with the fiducial resolution parameters.
Fig. 11 shows both the star formation rates averaged over 10 Myr (large panel) and the relative computing time (small panel on the right) for simulations with a particle mass of  B = 10 5 M ⊙ and a star formation efficiency of  sf = 0.32 per cent.Variations in ℎ min (line colours) for isolated galaxies with a softening length of  = 500 pc (solid lines) do not have a large impact on neither the star formation history nor the computing time, because the large softening reduces the amount of high-density gas.For simulations with a softening length of  = 250 pc (dashed lines), the computing time is reduced by a factor of ≈ 2 when reducing ℎ min by a factor of 10 (from ℎ min = 77.5 pc, red dashed line, to ℎ min = 7.75 pc, green dashed line).
The increased computing time for simulations with larger values of ℎ min and therefore more clumping from runaway collapse is caused by substantial amounts of shock-heated gas at high densities log  H [cm −3 ] ≳ 2) and temperatures of a few hundred to a few thousand K (see first and third column of Fig. 10).The temperature increase of a factor of 100 compared to simulations with smaller values of ℎ min (second and fourth column of Fig. 10) reduces the timestep size Δ ∝  1/3 B  −1/3 H  1/2 (see e.g.Borrow et al. 2022) by a factor of 10.The higher "real" densities in artificially dense gas clumps in the fiducial run do not directly affect the timestep size because Δ is calculated from the (underestimated) SPH densities.
Summarizing, a too large value for ℎ min does not only produce artificially dense clumps of gas and star particles, but also increases the computing time by a factor of ≈ 2 in our simulations.

DISCUSSION
The runaway collapse zone, as defined in section 3.1 and in particular equations ( 9) and (10), is an approximation because the clumping will depend on the exact particle configuration.Nevertheless, the isolated galaxy simulations show a clear correlation between the amount of gas in this zone and the presence of both artificial collapse (seen in the recalculated gas densities) as well as the number of dense stellar clumps, which is drastically reduced for smaller values of the minimum smoothing length (compare centre and bottom left panels in Fig. 9).
The star formation efficiencies are varied in the presented simulations with values of  sf = 0.32 per cent and  sf = 1 per cent to illustrate the discussed issues.Higher values of  sf mean that gas particles are converted into star particles on shorter timescales for a given density.The isolated galaxy simulations from this work do not produce artificial clumps for  sf > 1 per cent, because dense gas is quickly converted into stars and therefore very little dense gas is present.Cosmological simulations of galaxy formation include the evolution of a large variety of galaxies with diverse properties and mass accretion histories and can occupy different density-temperature regions at different times.It is therefore plausible that the discussed issues are present in cosmological simulations also for  sf ≥ 1 per cent.
One could expect that smaller values for the gravitational softening scale  always lead to more accurate simulations than larger values because gravitational forces would be calculated correctly up to higher densities, but we discussed in section 3.2 that the density -temperature zone within which the gravitational instabilities are modelled correctly ( J,N >  smooth ) only depends on the particle mass (i.e. the mass resolution).Gas with densities and temperatures such that  J,N <  smooth , is formally unresolved independently of the softening length, which for a baryon particle mass of  B = 10 5 M ⊙ includes most of the cold neutral and molecular gas in galaxies (see the left panel of Fig. 5).A simulation with adaptive softening might report using  min = 2 pc and a simulation with the same mass resolution of  B = 10 5 M ⊙ but constant softening might use  = 100 pc.Despite their very different values for  and  min , neither simulation models gravitational instabilities correctly in the CNM nor in the MCs but both need to choose values for  or  min that avoid numerical issues.As described in section 3.1, this means a very small value for  min to avoid numerical runaway collapse for adaptive softening (if  min = ℎ min ) and potentially a larger value for  for constant softening if one wishes to suppress artificial sub-kernel gravitational instabilities (see section 3.2).
In very high (mass and force) resolution simulations, clusters of stars can form for physical reasons, but for the mass resolution in the isolated galaxies examples ( B = 10 5 M ⊙ ,  B = 8 × 10 5 M ⊙ ), we showed that the high-density gas clumps and resulting dense star clusters are mainly numerical artefacts because their numbers and masses are drastically reduced in simulations with better hydrodynamic force resolution (compare the centre and bottom left panels in Fig. 9).
For cosmological simulations the choice of the resolution parameters (mass, gravity, hydro) also need to take into account the effect of artificial heating of the stellar disk by dark matter and stellar halo particles (see e.g.Ludlow et al. 2019aLudlow et al. , 2021Ludlow et al. , 2023;;Wilkinson et al. 2022), which is reduced for larger values of  (appendix D in Ludlow et al. 2021).This effect is not included here because we do not use a live dark matter halo for these idealised simulations.
The star formation rate for a gas particle in our simulation depends on the Newtonian free-fall time (see equation 14) but the collapse of gas clumps in the simulation follows the longer softened free-fall time (equation A34).While for numerical stability considerations the softened properties (free-fall time, Jeans mass and length) are more relevant, we use the Newtonian free-fall time for the star formation rate because we want to model this subgrid process physically rather than numerically (see also the discussion in Nobels et al. 2023).
Which resolution and star formation criteria to choose depends on the application of the simulation, but the softened Jeans criteria and the temperature -density zones defined in section 3, in particular the runaway collapse zone (section 3.1) provide a useful guideline to avoid undesired numerical behaviour that may introduce artefacts and can furthermore slow down the simulation.
To date most simulations of cosmological volumes that reach  = 0 limit the pressure (or temperature) of the ISM by an effective pressure floor  ∝   eff which is related to the gas density  through the effective polytropic index  eff .It has been argued that for  eff > 4/3, the Jeans mass does not increase with density and therefore prevents artificial fragmentation (e.g.Schaye & Dalla Vecchia 2008).In eagle (Schaye et al. 2015)  eff = 4/3 for  H > 0.1 cm −3 and in the widely used Springel & Hernquist (2003) model, the polytropic index varies with density between  eff ≈ 2.5 at 0.1 cm −3 to  eff ≈ 4/3 at 100 cm −3 (for  = 0; see figure 1 of Springel & Hernquist 2003 for details).We argue here that such pressure (or entropy) floors are unnecessary in modern Lagrangian simulations, even for relatively low resolutions, provided that appropriate values for the softening and minimum smoothing lengths are used to avoid artificially inducing gravitational instabilities and a numerical runaway collapse.
The free-fall time in simulations with a constant softening length ( ff,s , equation 3) is longer than the Newtonian free-fall time ( ff,N , equation 2) in simulations with an adaptive softening length.In simu- lations of galaxy formation, the stellar feedback processes might have to start earlier to efficiently counteract the gravitational collapse when the gravitational softening is adaptive.Early feedback processes, such as radiation and stellar winds from young stars, or a small delay between star formation and core-collapse supernova explosions, might therefore be more important in simulations with adaptive softening than in simulations with a constant softening length.

SUMMARY
The Jeans stability criteria based on the analysis by Jeans (1902) define a length scale,  J,N , the so called Jeans length, by comparing the free-fall timescale to the sound crossing timescale (see the derivation in appendix A).In a homogeneous medium with gas density  H and temperature , density perturbations on scales >  J,N are gravitationally unstable (i.e. the density perturbations grow) while density perturbations on scales <  J,N are gravitationally stable (i.e. the density perturbations decay with time until a new equilibrium is reached).
In this work we introduce "softened Jeans criteria" (section 2) for which we re-derive the Jeans length in softened gravity, as used in Lagrangian simulations to suppress 2-body scattering.In parallel to the Newtonian Jeans criteria, the softened Jeans length (mass),  J,s ( J,s ) describes the minimum length (mass) scale above which density perturbations grow and become gravitationally unstable in simulations with softened gravity.
For gas with densities and temperatures for which the Newtonian Jeans length is smaller that the gravitational softening length (densities above or temperatures below the red dashed line,  J,N = , in Fig. 1), gravitational fragmentation is described by the softened Jeans criteria instead of the Newtonian Jeans criteria.The further gravitational collapse is slowed down in softened gravity and bet-ter described by the softened free-fall time (equation 3) than the Newtonian free-fall time (equation 2).
Depending on the gravitational softening length, , the softened Jeans mass can exceed the Newtonian Jeans mass  J,N by several orders of magnitude for gas densities and temperatures typical of the cold interstellar medium.For example, at a gas temperature of a few tens of K and a gas density of ≈ 100 cm −3 , the Newtonian Jeans mass is  J,N ≈ 10 2 M ⊙ , and the softened Jeans mass for a constant softening length of  = 100 pc is  J,s ≈ 10 5 M ⊙ (Fig. 1).
In simulations with particle masses of ≈ 10 5 M ⊙ that do not impose an artificial pressure or entropy floor in the form of an effective equation of state, cold neutral and molecular ISM gas is formally unresolved because gravitational instabilities should form within an individual smoothing kernel ( J,N <  smooth ).Modelling a multi-phase interstellar medium at these mass resolutions therefore requires an understanding of how instabilities behave at and below the resolution limit (see section 3).
If the hydrodynamic forces are resolved by at least one smoothing kernel, gravitational instabilities would be modelled correctly (if  J,s =  J,N ), or be suppressed ( J,s >  J,N ) for all gas densities and temperatures in simulations, because softened gravity never exceeds Newtonian gravity and therefore  J,s ≥  J,N .Yet, we show in section 3 that perturbations can grow under specific conditions within a hydrodynamic smoothing kernel.Two distinct pathways are identified for numerically induced instabilities, both related to an inaccurate calculation of the hydrodynamic properties below the resolution limit: instabilities caused by (i) an underestimated gas density (subsection 3.1) and (ii) pressure that is smoothed on length scales larger than those on which gravity is softened (subsection 3.2).The former instability is relevant for simulations with both adaptive and constant softening lengths, while the latter only applies to simulations with a constant softening length (see discussion in section 3).
The effects of sub-kernel instabilities is demonstrated in section 4 in simulations of isolated galaxies using the Swift code with a constant softening length, .As outlined in subsection 3.1, the density of gas clumps that are smaller than the minimum allowed smoothing length, ℎ min , are under-estimated.Re-calculating the densities of all gas particles with a vanishing value for ℎ min reveals gas clumps with several orders of magnitude higher gas densities (compare top and bottom row of Fig. 8).These gas clumps result in dense star clusters that largely disappear for smaller minimum smoothing length values and are therefore artificial (Fig. 9).
Based on the analysis in section 3.1, we recommend for both simulations with constant and adaptive softening lengths to set the minimum smoothing lengths to a value small enough value that the smoothing lengths are not limited to a constant value in a simulation with particle mass,  B and a maximum expected density,  sim,max : We recommend to generously fulfill this condition and could not identify a benefit of a larger value for ℎ min .If the minimum values for softening and smoothing lengths are identical ( min = ℎ min ) in simulations with adaptive softening, a small value for  min is also required.Gravitational instabilities that would form within a smoothing kernel are unresolved by definition, even if both the minimum smoothing length and the softening length are approaching zero 12 .Assuming 12 A simulation with a given particle mass  B does not have a higher res-ℎ min → 0, the fundamental differences that remain between simulations with adaptive and constant softening lengths for instabilities within a kernel are described in section 3.2.Neither method gives the right solution for sub-kernel gravitational instabilities: an adaptive softening length that follows the smoothing length of the kernel leads to artificial suppression of sub-kernel instabilities while a constant softening length can result in artificially inducing sub-kernel instabilities.Because these length scales are by definition below the resolution limit of the simulation, their detailed impact on galaxy properties needs further investigation and is beyond the scope of this work.
Suppressing sub-kernel instabilities in simulations with constant softening lengths requires a softening length that exceeds the value given by equation 13.However, increasing the softening length artificially stabilizes gravitational instabilities at the hydrodynamical spatial resolution limit: the size of a smoothing kernel,  smooth , for larger regions in density-temperature space (see the middle panel in Fig. 5).The impact of sub-kernel instabilities on galaxy properties remains to be tested, because cooling rates, star formation rates and other density-or pressure-dependent sub-grid models use the smoother SPH estimates.The optimal value of a constant softening length, , will therefore depend on the application.
Finally, we argue that SPH simulations with relatively low baryon mass resolution (shown here up to  B = 8 × 10 5 M ⊙ but the dependence on  B is weak) do not depend on an effective pressure floor for numerical stability 13 .Simulations with comparable mass resolutions and without an artificial pressure floor do not suffer from the numerical issue discussed in section 3.1 if ℎ min satisfies equation (15).While the collapse of individual gas clouds is suppressed or slowed down when gravity is softened, this can be taken into account by subgrid prescriptions, e.g. for star formation, and is generally not a bottleneck if the aim of a simulation is to reproduce general galaxy properties.
The ideal simulation has both the mass and force resolution to accurately model the gravitational fragmentation and collapse of individual molecular clouds but in this work we showed a numerically stable alternative, both for adaptive and constant softening, for projects for which this is computationally too expensive.Milky Way (indicated as dark patch at 10 2 ≲  H [cm −3 ] ≲ 10 6 ,  ≈ 20 K).Steinwandel et al. (2023) mention that the global galactic properties are not very sensitive to variations in the softening length but that reducing  to below 0.1 pc, leads to numerical issues caused by runaway stars with very high velocities.For the same particle mass ( B = 4 M ⊙ ), Hu et al. (2017) varied the softening length between  = 0.5 pc and  = 2 pc in their appendix B2.They find that the gas density distribution in their simulations extends to higher densities for smaller values of , but that integrated galaxy properties, such as the total star formation rate only has a weak dependence on .It would be interesting to see if the discussed sub-kernel clumping is visible upon closer inspection in these simulations with very small values for , but such an analysis is beyond the scope of this work.
Based on Fig. B1, sub-kernel clumping can be avoided for simulations with  B = 4 M ⊙ , for softening lengths of  ≥ 2 pc, as e.g.used in Hu et al. (2016Hu et al. ( , 2017)).

APPENDIX C: ROBUSTNESS OF RESULTS
The simulations presented in section 4 illustrate the process of a numerical runaway collapse as outlined in section 3.1.In this section we vary the energy per supernova (C1) and the thermal state of the interstellar medium (C2) to demonstrate that our conclusions are robust to these changes.

C1 Supernova energy
The simulations in section 4 use an energy per supernova of  SN = 10 51 erg.This energy is stochastically injected into the nearest gas particle by increasing its temperature by Δ heat = 10 7.5 K. Dalla Vecchia & Schaye (2012) estimate the maximum gas density up to which thermal feedback is efficient by comparing the gas cooling time,  c , to the sound-crossing time,  sc , across a heated resolution element.If  c ≫  sc (e.g. by a factor of  t ≡  c /  sc = 10), the gas expands adiabatically before the thermal energy is lost through radiative cooling.
As discussed in section 4, the maximum density below which thermal feedback is efficient is ≈ 10 cm −3 (equation 18 from Dalla Vecchia & Schaye 2012) for the fiducial simulations.In Fig. C1, we show simulations with a 2 and 4 times higher energy per supernova, while keeping the heating temperature constant at Δ heat = 10 7.5 K.In the stochastic model, this means more frequent thermal energy injections.For reference, Schaye et al. (2015) used an energy per supernova that varies between  SN = 0.3 × 10 51 erg and 3 × 10 51 erg depending on metallicity (higher  SN for lower metallicity) and gas density at the time of star formation (higher  SN for higher gas density).
Increasing  SN (10 51 erg, 2 × 10 51 erg, 4 × 10 51 erg, columns 1 to 3) indeed reduces the clumping of the stars as well as the amount of shock-heated gas, even for large values of the minimum smoothing length (here: ℎ min = 77.5 pc).This agrees with the prediction from section 3.1 that the artificial clumping depends on the amount of gas in the runaway collapse zone.A larger value for  SN reduces the highest density gas within the galaxy and for the largest energy per supernova ( SN = 4 × 10 51 erg, 3rd column), barely any molecular gas (with  H > 100 cm −3 ) is left within the galaxy.On the other hand, a smaller value of ℎ min reduces the artificial clumping of stars also in the presence of cold gas.

C2 Photoheating and cosmic rays
The simulations presented in section 4 use the radiative cooling and heating rates from Ploeckinger & Schaye (2020).The rates are pre-calculated, assuming a pressure-dependent interstellar radiation field (ISRF) and cosmic ray (CR) rate (see Ploeckinger & Schaye 2020 for details) and a redshiftdependent radiation background from distant galaxies.In addition to the fiducial table ("UVB_dust1_CR1_G1_shield1"), Ploeckinger & Schaye (2020) also provide tables without an ISRF or CRs ("UVB_dust1_CR0_G0_shield1"), and with a 10 times higher normalization for the ISRF and CR rates ("UVB_dust1_CR2_G2_shield1").
The fiducial simulation ( B = 10 5 M ⊙ ,  = 250 pc) was rerun without an ISRF ("No ISRF") which also does not include CRs, and with the cooling and heating rates that correspond to the higher ISRF and CR normalization ("Strong ISRF").Fig. C2 demonstrates that the features that we have identified for large values of the minimum smoothing length, ℎ min : (i) stellar clumps and (ii) shock-heated, highdensity ( H ≳ 100 cm −3 ) gas with temperatures of a few hundred to a few thousand K, are insensitive to variations in the radiation field, if the value of ℎ min is large (77.5 pc, columns 1 to 3) and many gas particles are located within the runaway collapse zone.In turn, the presence of both features is drastically reduced, when reducing the value of ℎ min to 7.75 pc (columns 4 to 6), confirming the findings from section 4.   10.The drastically reduced clumping in stars, as well as the reduced amount of dense shock-heated gas is apparent for simulations with different cooling functions (no ISRF: 1st and 4th column, fiducial ISRF: 2nd and 5th column, strong ISRF: 3rd and 6th column) when reducing the minimum smoothing length, ℎ min from 77.5 pc (columns 1 to 3) to 7.75 pc (columns 4 to 6).Typical densities and temperatures for the WNM, CNM, and MCs are indicated with dark patches, as in Fig. 4.

Figure 1 .
Figure 1.The contours show the Jeans mass in units of log  J [M ⊙ ] for a Newtonian gravitational potential (top panel), a Plummer softened (bottom panel, solid lines) and a Wendland C2 softened (bottom panel, dotted lines) potential.The red dashed line in the bottom panel indicates where the Newtonian Jeans mass equals the Plummer softening scale  (here:  = 100 pc).

Figure 2 .
Figure2.Random distributions of particles (particle mass  B = 10 5 M ⊙ ) in 3D that represent input densities from log  H [cm −3 ] = 2 (leftmost panels) to log  H [cm −3 ] = 6 (rightmost panels).The density at the centre of each distribution is calculated with a Wendland C2 kernel and a smoothing length which is the maximum of ℎ correct (red solid circle) and ℎ min (blue dashed circle) and indicated in the top right of each panel ( H,SPH ).In the top panels the dense gas clump (black circles) consists of 100 particles (≈ number of neighbours in the kernel) and the bottom panels represent a more homogeneous medium without clumping on scales below ℎ min .If the particles clump on scales below ℎ min , then the density that the SPH solver calculates is too low (see top right panel) but if the gas distribution is smooth on scales smaller than ℎ min , then the SPH density estimate is still accurate, even if ℎ correct ≪ ℎ min (see bottom right panel).

Figure 3 .
Figure 3. Lines indicate the borders of the runaway collapse zone (see text)for various combinations of the baryon particle mass  B , the gravitational softening length  , and the minimum smoothing length ℎ min .At densities above the indicated limits (both the vertical and the inclined part), the SPH density estimate may be incorrect and a runaway collapse can occur.Each panel varies one resolution parameter (top:  B , middle:  , bottom: ℎ min ) while keeping the remaining two constant at the fiducial parameters of the eagle (100 Mpc) 3 simulation (black solid line, see labels).The thick grey line represents the effective temperature used in eagle for gas with densities > 0.1 cm −3 .

Figure 4 .
Figure 4. Overview of the density-temperature phase-space for the eagle, eagle-high-res, and apostle-L1 resolution parameters, as well as for an alternative set (from top to bottom, see also Table2).The contours indicate the softened Jeans mass in units of log  J,s [M ⊙ ].The red shaded region at high densities is the runaway collapse zone and SPH-estimated gas densities within this zone may be underestimated (see section 3.1).The grey thick line is the effective temperature floor used in eagle, eagle-high-res, and Apostle-L1.Typical densities and temperatures for the warm neutral medium, the cold neutral medium and molecular clouds (from low to high densities) are indicated with black patches for reference.An alternative set of parameters is shown in the bottom panel: for a mass resolution of order 10 5 to 10 6 M ⊙ , a softening scale of 200 pc and a minimum smoothing length of 2 pc (ℎ min,ratio = 0.006) push the runaway collapse zone to densities ≳ 10 8 cm −3 ).All phases in the ISM can therefore be modelled without triggering a numerical runaway collapse.
. The left panel shows the boundary  J,s =  smooth for  B = 10 5 M ⊙ and a constant softening length  soft = 100 pc (thick solid line).In the middle panel the solid lines represent  J,s =  smooth for different values for  soft (see contour labels) and  B = 10 5 M ⊙ .The right panel shows the dependence of both  J,N =  smooth (thick green lines) and  J,s =  smooth (thin black lines) on the particle mass  B for  soft = 100 pc and  B = 10 5 /8 , 10 5 , and 10 5 × 8 M ⊙ (dotted, dashed, dash-dotted lines, respectively).Typical densities and temperatures for the WNM, CNM, and MCs are indicated with dark patches, as in Fig.4.

Figure 6
Figure 6.Gravitational stability at the hydrodynamical resolution limit, max( smooth ,  smooth,min ), for simulations with an adaptive softening length,  soft .Line styles as in the left panel of Fig.5but for an adaptive softening length,  soft .In the left panel, the minimum value for the gravitational softening length equals the minimum value for the smoothing length.Here, the slope of both lines (dashed line:  J,N = max( smooth ,  smooth,min ), solid line:  J,s = max( smooth ,  smooth,min )) changes when  smooth is limited by  smooth,min .In the right panel, the softening length is adaptive down to a minimum value  soft,min , for which the softening length is effectively constant.In this panel, the slope of the solid line,  J,s = max( smooth ,  smooth,min ), changes for densities above which the softening length,  soft , is limited by a minimum value,  soft,min .Values for  B ,  soft,min , and  smooth,min were selected to represent FIREbox (left panel) and TNG50 (right panel, see text for details).As in Fig.5, the gas density,  H , refers to the SPH-estimated density.Typical densities and temperatures for the WNM, CNM, and MCs are indicated with dark patches, as in Fig.4.

Figure 7 .
Figure 7.As Fig. 5 but focusing on densities and temperatures for which gravitation instabilities on sub-kernel scales are expected ( J,s <  smooth ).Left panel: The dotted vertical line indicates the boundary  soft =  smooth for  B = 10 5 M ⊙ and a constant softening length of  soft = 100 pc.If  soft <  smooth gas clumps within a smoothing kernel can fragment further, but this process is suppressed at higher densities where  soft >  smooth (see text for details).In the right panel all lines are repeated for different values of  soft (see contour labels).As in Fig. 5, the gas density,  H , refers to the SPH-estimated density.Typical densities and temperatures for the WNM, CNM, and MCs are indicated with dark patches, as in Fig. 4.

Figure 8 .
Figure 8. Illustration of the runaway collapse described in section 3.1.Temperature-density histograms at time  = 1 Gyr for 3 simulations (columns 1 to 3) with identical initial conditions.The densities (x-axis) in the top row are the SPH gas density estimates,  H,SPH , as used in the simulations and stored in the snapshots, while the bottom row shows the re-calculated densities for ℎ min → 0, based on the same particle positions.The colour of the temperature-density histograms is proportional to the log of gas mass per pixel and the minimum and maximum values for the colour maps are constant across all plots.From left to right the minimum smoothing length, ℎ min , decreases from 77.5 pc (ℎ min,ratio = 0.2; first column) and ℎ min = 24.5 pc (second column) to ℎ min = 7.75 pc (third column).The red shaded area is the runaway collapse zone as defined in section 3.1 and each panel includes a number for the total mass of particles within this zone as log  zoneC [M ⊙ ] if  zoneC > 0. For reference, the black lines indicate where the softened Jeans mass is resolved by 1 (dotted), 10 (dashed), and 100 (solid) particles.The panels in the rightmost column show the cumulative mass fraction for particles above a given gas density  H .The solid lines in the top panel show the density distributions for the phase-space diagrams in the top row and the dashed lines in the bottom panel for the density distribution from the bottom row (dashed lines are repeated for reference in the top panel).The short vertical dotted lines at the top of the figures in the rightmost row indicate the densities  H,hmin (equation 8) above which the smoothing length is limited by ℎ min .

Figure 9 .
Figure 9.The distribution of gas (left) is represented as a (recalculated) density -temperature histogram (as in the bottom row of Fig. 8) and the distribution of stars (right) as a stellar mass surface density map for simulations with various resolution and star formation efficiency parameters at  = 1 Gyr.The central panel in each figure shows the results for the fiducial parameters  B = 10 5 M ⊙ ,  = 250 pc, ℎ min = 77.5 pc (ℎ min,ratio = 0.2), and  sf = 0.32 per cent and the panels in each corner vary one parameter while keeping the others constant: top right:  is increased by a factor of 2, bottom right:  B is increased by a factor of 8, bottom left: ℎ min is decreased by a factor of 10, top left:  sf is increased by a factor of 3. The phase-space plots on the left hand side show the zones where the numerical instability as described in section 3.1 can form as red-shaded areas.The total gas mass in this zone is shown as log  zone [M ⊙ ] in each panel, if  zone > 0. The right hand side consists of images of the stellar surface mass density within  = 25 kpc (excluding stars already present in the initial conditions).The clumps identified by the friends-of-friends algorithm are highlighted with circles (diamonds) if their mass is below (above) 10 8 M ⊙ .Both colour maps are logarithmic and span the same range across all simulations.

Figure 11 .
Figure 11.Main panel: Star formation histories (SFHs) for all simulations with a particle mass  B = 10 5 M ⊙ and a star formation efficiency of  sf = 0.32 per cent.Lines of different colours show the star formation rates (SFR, averaged over 10 Myr) of simulations with different minimum smoothing lengths, ℎ min , and different line styles indicate different gravitational softening lengths,  .The red dashed line is the SFH for the simulation with the fiducial parameters (centre of Fig. 9, first column in Fig. 10) and the green solid line shows the SFH for the simulation with the alternative set of parameters in the second column of Fig. 10.The small panel on the right gives an overview of the relative computing times normalized to the fastest simulation with the same linestyles as in the main panel.

Figure B1 .
Figure B1.Lines and colours are as in the right panel of Fig. 7 but for simulations with much higher mass resolution ( B = 4 M ⊙ ) and constant softening lengths of  soft = 0.15, 0.75, and 1.5 pc ( = 0.1, 0.5, and 1 pc).Typical densities and temperatures for the CNM and MCs are indicated with dark patches, as in Fig. 4.

Figure C1 .
Figure C1.Stellar mass surface density maps (top row) and gas densitytemperature histograms for ℎ min = 77.5 pc, as in Fig. C2.Increasing the energy deposited per supernova from the canonical value of  SN = 10 51 erg (1st column) by a factor of 2 (2nd column) and by a factor of 4 (3rd column) reduces the amount of gas in the runaway collapse zone (red shaded area), resulting in reduced stellar clumping.Typical densities and temperatures for the WNM, CNM, and MCs are indicated with dark patches, as in Fig. 4.

Figure C2 .
Figure C2.Stellar mass surface density maps (top row) and gas density-temperature histograms, as in Fig.10.The drastically reduced clumping in stars, as well as the reduced amount of dense shock-heated gas is apparent for simulations with different cooling functions (no ISRF: 1st and 4th column, fiducial ISRF: 2nd and 5th column, strong ISRF: 3rd and 6th column) when reducing the minimum smoothing length, ℎ min from 77.5 pc (columns 1 to 3) to 7.75 pc (columns 4 to 6).Typical densities and temperatures for the WNM, CNM, and MCs are indicated with dark patches, as in Fig.4.