Resolution criteria to avoid artificial clumping in Lagrangian hydrodynamic simulations with a multiphase interstellar medium

ABSTRACT

Large-scale cosmological galaxy formation simulations typically prevent gas in the interstellar medium (ISM) from cooling below |$\approx 10^4\, \mathrm{K}$|⁠. This has been motivated by the inability to resolve the Jeans mass in molecular gas (⁠|$\ll 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 oversmoothing of the hydrodynamical properties relative to the gravitational force resolution. This effect is illustrated using simulations of isolated disc 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.

1 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 k_J as the limiting wavenumber between oscillatory solutions (k > k_J, short wavelengths) and exponentially growing solutions (k < k_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 W to the internal energy U of a spherical gas cloud. Perturbations can grow for scales where W + U < 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 r, sound speed c_s, and (unperturbed) density ρ are gravitationally unstable:

where G is the gravitational constant and A is a dimensionless pre-factor 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 two-body interactions. A classical prescription for softened gravity is the Plummer (1911) potential ϕ_P = −GM(r² + ϵ²)^−1/2 for a point mass, M, 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, Schaye & Bower 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 r ≲ ϵ.

A common technique to solve the equations of motions for a collisional fluid is smoothed particle hydrodynamics (SPH, originally from Gingold & Monaghan 1977; Lucy 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, h, as the characteristic length-scale, is related to the hydrodynamical force resolution.

The smoothing length, h, generally decreases with increasing gas density, h ∝ ρ^−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 ϵ ∝ h (’adaptive softening’) or ϵ can be set to a constant value for all densities. Adaptive softening is implemented, for example, 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 (≈h, the smoothing length of gas particles) or cell size and the hydrodynamic force resolution equals the gravitational force resolution.

A constant parameter to define the Plummer-equivalent softening length, ϵ, in SPH is used for example in the eagle (Schaye et al. 2015, code: gadget3, based on Springel 2005), the astrid (Bird et al. 2022, code: mp-gadget, Feng et al. 2018), and the flamingo (Schaye et al. 2023, code: swift Schaller et al. 2016, 2018, 2023) simulations.

In the moving mesh code arepo (Weinberger, Springel & Pakmor 2020), the softening length is adaptive and follows ϵ_Arepo = f_h(3V/(4π))^1/3, where f_h is an input parameter which defines the softening length relative to the cell radius, and V is the volume of the Voronoi cell. In the i llustrisTNG 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: |$\epsilon _{\mathrm{min}} = 72\, \mathrm{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 (ϵ = h) and ensuring that both resolution parameters can get small enough to resolve the smallest Jeans mass cloud (⁠|$M_{\mathrm{J}} = 4\pi \lambda _{\mathrm{J}}^3 \rho /3$|⁠) in the simulation by 2N_neigh particles (later reduced to 1.5N_neigh by Bate, Bonnell & Bromm 2003), where N_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 ϵ < h and artificially inhibited for ϵ > h) and not by the physical conditions.

Hubber, Goodwin & Whitworth (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 ϵ = h, confirming the analytical results from Whitworth (1998). Yamamoto, Okamoto & Saitoh (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 ϵ = h. 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 |$\epsilon = 700\, \mathrm{pc}$| for their |$(100\, \mathrm{Mpc})^3$| simulation with an initial baryon particle mass of |$1.81 \times 10^6\, \mathrm{M}_{\odot }$|⁠. This marginally fulfills the Jeans mass criterion set by Bate & Burkert (1997, but not their recommendation to set ϵ = h) for the warm neutral medium (WNM) with a typical Jeans mass of a few times |$10^7\, \mathrm{M}_{\odot }$|⁠, but would violate it drastically for the cold gas phase with typical Jeans masses below |$10^4\, \mathrm{M}_{\odot }$|⁠. An effective pressure floor, sometimes referred to as a polytropic equation of state (EOS), |$P_{\mathrm{eff}} \propto n_{\mathrm{H}}^{\gamma _{\mathrm{eff}}}$| for gas with a density of |$n_{\mathrm{H}} \gt 0.1\, \mathrm{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, P_eff, is chosen to be at an effective temperature |$T_{\mathrm{eff}} = 8000\, \mathrm{K}$| for |$n_{\mathrm{H}} = 0.1\, \mathrm{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 normalizations in simulations of a dwarf galaxy with a particle mass of |$750\, \mathrm{M}_{\odot }$| and showed that the lowest mass gaseous clumps that formed in their galaxy discs 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,¹ is implemented in almost² all simulations of cosmologically representative volumes that reach z = 0.

The resolution of the gaseous component in a simulation is consequently 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 multiphase 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 |$\ll 10^4\, \mathrm{M}_{\odot }$|⁠. 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 multiphase 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, h_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 h_min is illustrated using simulations of isolated disc galaxies with the modern hydrodynamics code swift³ (Schaller et al. 2018, 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₁₀ throughout this work.

2 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 two-body scattering from individual particles that represent a continuous distribution, the gravitational forces are softened for close interactions at a distance r smaller than the softening length ϵ in Lagrangian codes. The Newtonian gravitational acceleration, a_N ∝ r⁻² is replaced by the softened gravitational acceleration, for example, a_P ∝ r(r² + ϵ²)^−3/2 for a Plummer (1911) potential (see Table 1 for an overview). In the limit of large separations, r ≫ ϵ, the softened acceleration equals the Newtonian acceleration (a_P → a_N ∝ r⁻²) but for r ≪ ϵ, the softened acceleration approaches a_P ∝ r and therefore diverges from the Newtonian gravity as r → 0.

In simulations with a particle mass of |$m_{\mathrm{B}} \gtrsim 10^5\, \mathrm{M}_{\odot }$|⁠, 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.⁴ 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: t_{ff, s} = t_{ff, W, fit} from equation A34))

is longer than the free-fall time in Newtonian gravity for an object of radius R and is proportional to

for ϵ ≫ R. 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 R.

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 |$\epsilon = 100\, \mathrm{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 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\, \mathrm{cm}^{-3}$|⁠, the Newtonian Jeans mass is |$M_{\mathrm{J,N}}\approx 100\, \mathrm{M}_{\odot }$|⁠, while the softened Jeans masses for |$\epsilon = 100\, \mathrm{pc}$| is |$M_{\mathrm{J,s}}\approx 10^5\, \mathrm{M}_{\odot }$| (and for |$\epsilon = 1000\, \mathrm{pc}$|⁠, |$M_{\mathrm{J,s}}\approx 10^7\, \mathrm{M}_{\odot }$|⁠) for both M_{J, s} = M_{J, P, fit} (Plummer) and M_{J, s} = M_{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, t_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,

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.

3 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 (MCs) is approximated by subgrid prescriptions, for example, 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 time-steps, 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 Section 3.1, we identify the problematic behaviour of SPH-like simulation codes when the hydrodynamical smoothing length h is limited by a too large minimum value, h_min. The resulting instability is caused by an underestimate of the gas density and leads to artificial clumps of closely spaced particles. Setting h_min to a very small non-zero value resolves this issue and we cannot identify any disadvantages of letting h_min → 0 (to be discussed in detail in Section 4).

In Section 3.2, we analyse gravitational instabilities in simulations with either constant or adaptive softening lengths at and below the hydrodynamical resolution limit, that is, the size of an individual smoothing kernel.

3.1 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 stabilized 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, h_min ≠ 0, the hydrodynamic forces can be oversmoothed for high gas densities and we explain below how this can cause a numerical runaway collapse.

In this section, the smoothing length, h, and its lower limit, h_min, are defined as in the hydrodynamics solver sphenix (Borrow et al. 2022), implemented in the swift code.⁵ 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, X_H is the hydrogen mass fraction, m_B the baryon particle mass, m_H the hydrogen particle mass, n_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 X_H = 0.74.

The minimum smoothing length in swift is set by the dimensionless parameter h_{min, ratio} which is defined as the ratio between the kernel support, hγ_k, and the length-scale above which gravity is Newtonian, H = 3ϵ. For the Wendland C2 kernel, γ_k = 1.936492 (Dehnen & Aly 2012) and the minimum smoothing length is therefore |$h_{\mathrm{min}} = 1.55 \, \epsilon \, h_{\mathrm{min,ratio}}$|⁠. The density above which the hydrodynamic forces are oversmoothed is

Note that the critical density n_{H, hmin} can remain constant for simulations with different mass resolutions m_B, if h_min is adjusted accordingly.

The softened Jeans criteria derived in Appendix  A are therefore only valid for n_H ≤ n_{H, hmin}. For higher densities, the hydrodynamic forces can become inaccurate due to oversmoothing and this inaccuracy depends on the density contrast within h_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 n_H which increases from log n_H[cm⁻³] = 2 (leftmost panel) to log n_H[cm⁻³] = 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 |$h_{\mathrm{min}} = 15.5\, \mathrm{pc}$| (corresponding to e.g. |$\epsilon = 100\, \mathrm{pc}$| and h_{min, ratio} = 0.1) in this example.

Each panel in Fig. 2 lists the input density as well as the density n_{H, SPH} that the SPH solver calculates at the centre of each box using a Wendland C2 kernel. If h becomes smaller than h_min and the gas is clumpy on scales smaller than h_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 h_min, in which case the output densities n_{H, SPH} remain accurate, even for h ≪ h_min (see also Fig. 1 in Borrow, Schaller & Bower 2021 for an illustration).

To avoid a potential runaway 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:

• h < h_min (or: n_H > n_{H, hmin})

• fragmentation on scales below h_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} < h_min.

Assuming that h_min < ϵ (to avoid inducing numerical collapse, see e.g. Bate & Burkert 1997), the condition λ_{J, s} < h_min falls into the softened part of the phase space where H ≡ 3ϵ ≫ λ_{J, N} and we therefore approximate the softened Jeans mass (here for the Wendland C2 kernel, equation A36) with |$\lambda _{\mathrm{J,s}}\approx 0.27^{0.3} \, \lambda _{\mathrm{J,N}}^{0.4} \, H^{0.6}$| and the conditions (i) and (ii) for a potential runaway 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 m_B, Plummer equivalent softening length ϵ, and the minimum smoothing length h_min. The default values (thick, black lines) with log m_B[M_⊙] = 6.3, |$\epsilon = 700\, \mathrm{pc}$|⁠, and |$h_{\mathrm{min}} = 108.5\, \mathrm{pc}$| (i.e. h_{min, ratio} = 0.1) are those from the eagle|$(100\, \mathrm{Mpc})^3$| simulation and two of these parameters are constant in each panel, while the third parameter is varied (top panel: m_B, middle panel: ϵ, and bottom panel: h_min).

Each zone is limited through a combination of a vertical line, condition (i): n_H = n_{H, hmin} (equation 9) and a line of log T ∝ log n_H, that is, for a constant softened Jeans length of λ_{J, s} = h_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 h_min and can lead to a runaway collapse.

Varying the baryon particle mass m_B at constant ϵ and h_{min, ratio} (Fig. 3, top panel) only changes n_{H, hmin}, but condition (ii) remains unaffected. A better mass resolution is even slightly counter-productive in resolving higher densities because n_{H, hmin} ∝ m_B for constant h_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 m_B and constant ϵ on the other hand pushes the runaway collapse zone to significantly higher densities (from |$10^2$| to |$\gt 10^6\, \mathrm{cm}^{-3}$| when reducing h_min by a factor of 10, bottom panel of Fig. 3).

If the highest (expected) density in a simulation is n_{sim, max}, the minimum smoothing length should be set so that h(n_{sim, max}) ≥ h_min. With h 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 h_min and therefore recommend to use a very small but non-zero⁶ value for h_min that generously fulfills the conservative criterion from equation (11). Equation 11 does not depend on ϵ and applies to both simulations with adaptive and constant softening lengths.

3.1.1 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 |$m_{\mathrm{B}} = 1.81\times 10^6\, \mathrm{M}_{\odot }$| (eagle), |$m_{\mathrm{B}} = 2.26\times 10^5\, \mathrm{M}_{\odot }$| (eagle-high-res), and |$m_{\mathrm{B}} = 10^4\, \mathrm{M}_{\odot }$| (apostle).

The resolution parameters (see Table 2) used in both eagle flagship runs (eagle and 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 |$n_{\mathrm{sim,max}}\gtrsim 100\, \mathrm{cm}^{-3}$|⁠. However, for these simulations, the numerically induced runaway collapse was avoided because gas with densities above |$0.1\, \mathrm{cm}^{-3}$| follows a polytropic EOS where the effective temperature⁷ 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 |$m_{\mathrm{B}} = 10^4\, \mathrm{M}_{\odot }$|⁠, but the runaway collapse zone is barely affected by the higher mass resolution, compared to, for example, eagle. The effective temperature floor was therefore necessary for their values for h_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 M_{\mathrm{J,s}} \, [\mathrm{M}_{\odot }]$| 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 m_B is very weak, see top panel of Fig. 3). A Plummer equivalent softening of |$\epsilon =200\, \mathrm{pc}$| is small enough to model gravitational instabilities correctly in the WNM but still large enough for a minimum softened Jeans mass of |$\approx 10^{6}\, \mathrm{M}_{\odot }$| in all neutral phases of the ISM. Counterintuitively, a smaller softening length would decrease the softened Jeans mass in the cold neutral medium (CNM) and MCs, which means that the softened Jeans mass would be resolved by fewer particles in these phases. With a minimum smoothing length of |$h_{\mathrm{min}} = 2\, \mathrm{pc}$| (instead of |$54.25\, \mathrm{pc}$| in eagle-high-res), the runaway collapse zone only covers |$n_\mathrm{H}\gt 10^8\, \mathrm{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 multiphase ISM need to carefully select the combination of mass (baryon particle mass), gravity (softening scale), and hydrodynamic resolutions (minimum smoothing length) to avoid unwanted numerical artefacts.

3.2 Gravitational instabilities at and below the resolution limits

In Section 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 code-agnostic length-scales l_soft, over which gravitational forces are softened, and l_smooth (with a minimum of l_{smooth, min}), over which the hydrodynamical forces are smoothed. We explain briefly how they connect to the smoothing lengths h and the softening lengths ϵ in some codes as examples:

We use the kernel support as measure of l_smooth. In swift, h 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 h, |$l_{\mathrm{smooth}} =\gamma _{\mathrm{k}}\, h = 1.936492\, h$|⁠. In practice, the definition of the smoothing length varies for each code. The SPH code gadget4 (Springel et al. 2021) defines h as the finite support of the kernel and therefore l_smooth would be equal to h_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. fig. 1 in Springel et al. 2021). We therefore use l_soft = 1.5ϵ as an estimate for the softening length-scale.

For the following discussion, using a different measure for either l_soft or l_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 l_{smooth, min} → 0, if not explicitly mentioned otherwise.

3.2.1 Gravitational instabilities at the hydro resolution limit l_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, l_smooth. The hydrodynamical forces are inaccurate on scales below l_smooth and fragmentation on scales of λ_{J, N} < l_smooth is unresolved. For the Newtonian Jeans length, the condition λ_{J, N} = l_smooth depends on the gas density and temperature as well as on the particle mass, m_B (⁠|$l_{\mathrm{smooth}}\propto m_{\mathrm{B}}^{1/3} n_{\mathrm{H}}^{-1/3}$|⁠), but it does not depend on the gravitational softening length, l_soft. The condition λ_{J, N} = l_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} > l_smooth.

For gravitational instabilities at the hydrodynamic resolution limit, l_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, l_smooth, for simulations with constant and adaptive softening lengths.

3.2.2 Constant softening length

In nature (or in simulations with m_B → 0, l_soft → 0, and l_{smooth, min} → 0), density perturbations with a length-scale of l_smooth grow if λ_{J, N} < l_smooth, and decay until they reach a new equilibrium for λ_{J, N} > l_smooth. The border⁸ (i.e. case λ_{J, N} = l_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} > l_smooth (unshaded area).

For softened gravity with a constant softening length (left and right panels: |$l_{\mathrm{soft}} = 100\, \mathrm{pc}$|⁠, middle panel: various values for l_soft, see the labels), the solid lines (λ_{J, s} = l_smooth) separate perturbations with length-scales l_smooth that are gravitationally (un)stable in softened gravity, that is, in the simulation. Densities and temperatures for which gas is gravitationally stable in nature for perturbations at the resolution limit l_smooth are also gravitationally stable in simulations with softened gravity (white area in the left panel of Fig. 5 for a simulation with |$m_{\mathrm{B}} = 10^5\, \mathrm{M}_{\odot }$|⁠, |$l_{\mathrm{soft}} = 100\, \mathrm{pc}$|⁠, and |$l_{\mathrm{smooth,min}} \rightarrow 0\, \mathrm{pc}$|⁠).

Gas with densities and temperatures within the shaded areas (λ_{J, N} < l_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 l_smooth because the softened Jeans length exceeds the kernel size (λ_{J, s} > l_smooth) and therefore perturbations on scales of l_smooth decay. The middle panel shows that the boundary λ_{J, s} = l_smooth depends on the value for the constant softening length, l_soft. The right panel shows that increasing (dashed–dotted lines) or decreasing (dotted lines) the particle mass by a factor of 8 from the fiducial value of |$m_{\mathrm{B}} = 10^5\, \mathrm{M}_{\odot }$| (dashed lines), affects both the thick green lines for λ_{J, N} = l_smooth as well as the thin black lines for the condition λ_{J, s} = l_smooth, because |$l_{\mathrm{smooth}}\propto m_{\mathrm{B}}^{1/3}$| (right panel of Fig. 5 for |$l_{\mathrm{soft}} = 100\, \mathrm{pc}$|⁠).

3.2.3 Adaptive softening length

For adaptive softening, typically l_soft = l_smooth and therefore λ_{J, s} ≈ λ_{J, N}. The instability criteria for perturbations at the resolution limit, that is, at the kernel size l_smooth, therefore follow the instability criteria for Newtonian gravity. The left panel of Fig. 6 (line styles as in Fig. 5) illustrates this for parameters representative for the firebox (Feldmann et al. 2023) simulation project with |$m_{\mathrm{B}} = 6\times 10^4\, \mathrm{M}_{\odot }$|⁠. Their adaptive softening length is equal to the average gas particle separation (l_soft = 1.5ϵ(m_B/ρ)^1/3, for gas density ρ), down to a minimum value of |$l_{\mathrm{soft,min}} = l_{\mathrm{smooth,min}} = 2.25\, \mathrm{pc}$| (⁠|$\epsilon _{\mathrm{min}} = 1.5\, \mathrm{pc}$|⁠). Density perturbations of length-scale l_smooth therefore follow the Newtonian Jeans criteria. The small offset between the dashed and solid lines is from the order of unity pre-factors related to the kernel shapes and exact definitions of l_smooth and l_soft. For simplicity, we use the same kernel shapes and pre-factors as in Fig. 5.

The slope change of the λ_{J, N} = l_smooth (dashed) and λ_{J, s} = l_smooth (solid) lines is caused by the minimum smoothing length l_{smooth, min}. The critical density n_{H, hmin} above which densities are oversmoothed, and hence n_{H, SPH} < n_H, is |$\approx 2\times 10^6\, \mathrm{cm}^{-3}$| [equation (8) for |$m_{\mathrm{B}} = 6\times 10^4\, \mathrm{M}_{\odot }$| and |$h_{\mathrm{min}}\approx 1.2\, \mathrm{pc}$|]. The small values for l_{smooth, min} and l_{soft, min} therefore prevent the instability described in Section 3.1 for densities below |$\approx 2\times 10^6\, \mathrm{cm}^{-3}$|⁠.

In the illustrisTNG project (Nelson et al. 2018; Pillepich et al. 2018b) 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: |$l_{\mathrm{soft,min}} = 1.5 \epsilon _{\mathrm{gas,min}} = 108\, \mathrm{pc}$|⁠, TNG100: |$l_{\mathrm{soft,min}} = 1.5 \epsilon _{\mathrm{gas,min}} = 285\, \mathrm{pc}$|⁠, and TNG300: |$l_{\mathrm{soft,min}} = 1.5\epsilon _{\mathrm{gas,min}} = 555\, \mathrm{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\, \mathrm{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} = l_smooth (dashed) and λ_{J, s} = l_smooth (solid) lines for values representative for TNG50: |$m_{\mathrm{B}} = 8.5\times 10^4\, \mathrm{M}_{\odot }$|⁠, |$l_{\mathrm{soft,min}} = 108\, \mathrm{pc}$|⁠, and an assumed negligible value for l_{smooth, min}.

Comparing the right panel of Fig. 6 with the left panel of Fig. 5, we see that gravitational instabilities on scales <l_{soft, min} in simulations with adaptive softening lengths effectively behave as in simulations with constant softening lengths.

3.2.4 Gravitational instabilities below the hydro resolution limit l_smooth

Density perturbations grow in nature on length-scales smaller than l_smooth in gas with temperatures and densities for which λ_{J, N} < l_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 subkernel instabilities (λ_{J, s} < l_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 (SFRs), and other density- or pressure-dependent subgrid models use the smoother SPH estimates, we briefly discuss the differences in the treatment of subkernel perturbations between codes with constant and adaptive softening lengths.

3.2.5 Constant softening length

Gas with λ_{J, s} < l_smooth is expected to be gravitationally unstable when exposed to fluctuations on length-scales between λ_{J, s} and l_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 Section 3.1). If l_smooth < l_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, l_soft can be (much) smaller than the size of a smoothing kernel, for a constant value of l_soft. In this case, dense particle configurations with sizes possibly even smaller than λ_{J, N} can form because the density and pressure estimates of subkernel 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} < l_smooth, as in Fig. 5) is further split into densities for which further fragmentation to even smaller scales is suppressed (l_soft > l_smooth) or potentially induced (l_soft < l_smooth), with the boundary, l_soft = l_smooth, indicated as vertical, dotted lines (left panel: |$l_{\mathrm{soft}} = 100\, \mathrm{pc}$| and right panel: various values for l_soft, see contour labels).

We conclude that for a constant softening length, gas in simulations with λ_{J, s} < l_smooth and l_soft < l_smooth (purple areas in Fig. 7) might be affected by subkernel 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, l_smooth. In contrast to simulations with adaptive softening, the constant value for l_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 l_smooth = l_soft (vertical lines in Fig. 7). For higher gas densities (i.e. l_smooth < l_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 subkernel clumping (‘subk’) is undesired in simulations with constant softening, the conditions λ_{J, s} < l_smooth and l_soft < l_smooth (purple areas in Fig. 7) should not cover regions in density-temperature space that are populated by many gas particles in the simulation. If we want to confine the area in density–temperature space for which l_soft < l_smooth to densities below n_{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, subkernel gravitational instabilities in large parts of the cold ISM, but at the expense of suppressing physical instabilities at the kernel scale, l_smooth (see the middle panel of Fig. 5).

We focus in this work on simulations with |$m_{\mathrm{B}}\gtrsim 10^5\, \mathrm{M}_{\odot }$|⁠, 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 |$m_{\mathrm{B}} = 4\, \mathrm{M}_{\odot }$|⁠, for which a constant softening length of at least |$\epsilon = 2\, \mathrm{pc}$| is needed to avoid subkernel clumping.

3.2.6 Adaptive softening length

In simulations with an adaptive softening length and l_soft ≥ l_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} < l_smooth are therefore artificially suppressed.

In order to avoid the runaway collapse described in Section 3.1, the minimum smoothing length, l_{smooth, min}, needs to be small (see equation 11). In simulations, such as FIREbox, where l_{soft, min} = l_{smooth, min}, the minimum softening length needs to have the same small value as the minimum smoothing length.

4 GALAXY SIMULATIONS

In Section 3.2, we used the smoothing and softening lengths l_smooth and l_soft for a code independent discussion on the individual zones. In this section, we use simulations of isolated galaxies with the public swift code⁹ (Schaller et al. 2016, 2018, 2023, www.swiftsim.com). The resolution parameters set by the user are ϵ and h_min which relate to l_soft = 1.5ϵ and l_{smooth, min} = 1.94h_min, respectively (see the text in Section 3.2 for details). We will demonstrate the numerical issues that can occur for different choices of ϵ and h_min, focusing on simulations of isolated galaxies with a constant softening length. The 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, Di Matteo & Hernquist (2005a) but modified to use the exact definition of the analytic dark matter halo mass M₂₀₀ (i.e. mass within R₂₀₀, 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\times 10^{10}\, \mathrm{M}_{\odot }$| with solar metallicity (Z_⊙ = 0.0134, Asplund et al. 2009) and an exponential stellar disc with a radial scale length of |$4.3\, \mathrm{kpc}$| and a mass of |$3.836\times 10^{10}\, \mathrm{M}_{\odot }$| The initial disc gas fraction is 30 per cent and the gas is initialized with a temperature of |$10^4\, \mathrm{K}$|⁠. The baryonic disc is in equilibrium with an analytic dark matter halo with a Hernquist (1990) profile of |$M_{\mathrm{200}} = 1.37 \times 10^{12}\, \mathrm{M}_{\odot }$| and a scale radius such that the central density profile matches that of a Navarro, Frenk & White (1997) profile with a concentration of c = 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 (CRs), an interstellar radiation field (ISRF) 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}\, \mathrm{erg}$| per SN is injected stochastically as thermal energy, following Dalla Vecchia & Schaye (2012), with a heating temperature of |$10^{7.5}\, \mathrm{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}\, \mathrm{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, for example, |$10^5\, \mathrm{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 |$n_{\mathrm{H}} \gt 0.1\, \mathrm{cm}^{-3}$| and cold gas (temperatures of |$T\lt 1000\, \mathrm{K}$|⁠) and the SFR for each gas particle with mass m_gas is given by the Schmidt (1959) relation

with the star formation efficiency e_sf and the Newtonian free-fall time t_{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 |$\epsilon = [250, 500]\, \mathrm{pc}$|⁠, minimum SPH smoothing length |$h_{\mathrm{min}} = [7.75, 24.5, 77.5]\, \mathrm{pc}$| (h_{min, ratio} = [0.02, 0.063, 0.2] for |$\epsilon = 250\, \mathrm{pc}$|⁠), and baryon particle mass|$m_{\mathrm{B}}= [10^5, 8\times 10^5]\, \mathrm{M}_{\odot }$|⁠; as well as the star formation efficiency |$e_{\mathrm{sf}} = [0.32, 1]\, \mathrm{per\, cent}$| in order to change the amount of high-density gas in a controlled way.

All simulations run until |$t_{\mathrm{end}}=1.2\, \mathrm{Gyr}$|⁠. We use swift to recalculate the gas densities based on the gas particle positions for h_min → 0 by restarting the original simulations at |$t = 1\, \mathrm{Gyr}$| for one very small time-step (⁠|$\Delta t_{\mathrm{max}} = 10\, \mathrm{yr}$|⁠) with a very small minimum smoothing length¹⁰ (h_{min, ratio} = 10⁻⁴). The snapshot file that is produced after this time-step 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 h_min → 0 (bottom row) for 3 simulations (first three columns; all with |$m_{\mathrm{B}}=10^5\, \mathrm{M}_{\odot }$|⁠, |$\epsilon = 250\, \mathrm{pc}$|⁠, and |$e_{\mathrm{sf}}=0.32\, \mathrm{per\, cent}$|⁠) with decreasing values for h_min (from left to right). The runaway collapse described in Section 3.1 is obvious in the simulation with the largest value for h_min (first 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 h_min (third column), the two density–temperature diagrams look identical.

The panels in the rightmost column show the cumulative mass fraction of particles above density n_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 n_{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 (⁠|$h_{\mathrm{min}} = 77.5\, \mathrm{pc}$|⁠) has much higher ‘real’ (but unresolved) densities, up to |$n_{\mathrm{H}}\gt 10^6\, \mathrm{cm}^{-3}$|⁠, than the formally highest resolution simulation (⁠|$h_{\mathrm{min}} = 7.75\, \mathrm{pc}$|⁠). If these dense clumps would appear through physical processes such as instabilities in the galaxy disc, 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 SFRs are underestimated as well (⁠|$\dot{m}_{\star }\propto n_{\mathrm{H}}^{1/2}$|⁠, 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 |$n_{\mathrm{fb}} \approx 10\, \mathrm{cm}^{-3}$|⁠. Their derivation assumes h > h_min which is fulfilled at densities of |$\approx 10\, \mathrm{cm}^{-3}$| in all simulations in this work. Comparing the values of n_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 |$t=1\, \mathrm{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 |$m_{\mathrm{B}}=10^5\, \mathrm{M}_{\odot }$|⁠, |$\epsilon = 250\, \mathrm{pc}$|⁠, and |$h_{\mathrm{min}} = 77.5\, \mathrm{pc}$| and a star formation efficiency of |$e_{\mathrm{sf}}=0.32\, \mathrm{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 standalone 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 |$m_{\mathrm{B}}=10^5\, \mathrm{M}_{\odot }$| and 18 for simulations with |$m_{\mathrm{B}}=8 \times 10^5\, \mathrm{M}_{\odot }$|⁠. 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 (⁠|$M_{\mathrm{clump}} \gt 10^8\, \mathrm{M}_{\odot }$|⁠) are highlighted by diamonds.

The galaxy with the fiducial resolution parameters (centre) has several stellar clumps due to the large amount of mass (log M(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 m_B; lower right). On the other hand, increasing the gravitational force softening length (increase ϵ; upper right) or increasing the star formation efficiency (increase e_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 different,¹¹ while the simulation with the higher star formation efficiency (⁠|$e_{\mathrm{sf}}=1\, \mathrm{per\, cent}$|⁠) has very little cold (⁠|$\lt 1000\, \mathrm{K}$|⁠) gas, the galaxy with the larger softening length (⁠|$\epsilon = 500\, \mathrm{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 disc.

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 h_min. For a particle mass of |$m_{\mathrm{B}} = 10^5\, \mathrm{M}_{\odot }$| (first two columns in Fig. 10), the alternative parameter set fulfills the conditions in equation (11) for h_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 h_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, |$m_{\mathrm{B}} = 10^5\, \mathrm{M}_{\odot }$|⁠, left two panels, and |$m_{\mathrm{B}} = 8 \times 10^5\, \mathrm{M}_{\odot }$|⁠, right two panels) and all simulations use the same value for the star formation efficiency (⁠|$e_{\mathrm{sf}} = 0.32\, \mathrm{per\, cent}$|⁠). The gravitational softening length is increased from 250 pc (fiducial) to 500 pc (alternative). The minimum smoothing length is reduced from |$h_{\mathrm{min}} = 77.5\, \mathrm{pc}$| (fiducial) to |$h_{\mathrm{min}} = 7.75\, \mathrm{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 (⁠|$\lt 10^2\, \mathrm{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 SFRs averaged over 10 Myr (large panel) and the relative computing time (small panel on the right) for simulations with a particle mass of |$m_{\mathrm{B}} = 10^5\, \mathrm{M}_{\odot }$| and a star formation efficiency of |$e_{\mathrm{sf}} = 0.32\, \mathrm{per\, cent}$|⁠. Variations in h_min (line colours) for isolated galaxies with a softening length of |$\epsilon =500\, \mathrm{pc}$| (solid lines) do not have a large impact on neither the star formation history (SFH) nor the computing time, because the large softening reduces the amount of high-density gas. For simulations with a softening length of |$\epsilon =250\, \mathrm{pc}$| (dashed lines), the computing time is reduced by a factor of ≈2 when reducing h_min by a factor of 10 (from |$h_{\mathrm{min}} = 77.5\, \mathrm{pc}$|⁠, red dashed line, to |$h_{\mathrm{min}} = 7.75\, \mathrm{pc}$|⁠, green dashed line).

The increased computing time for simulations with larger values of h_min and therefore more clumping from runaway collapse is caused by substantial amounts of shock-heated gas at high densities log n_H[cm⁻³] ≳ 2) and temperatures of a few hundred to a few thousand K (see first and third columns of Fig. 10). The temperature increase of a factor of 100 compared to simulations with smaller values of h_min (second and fourth columns of Fig. 10) reduces the time-step size |$\Delta t \propto m_{\mathrm{B}}^{1/3} n_{\mathrm{H}}^{-1/3} T^{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 time-step size because Δt is calculated from the (underestimated) SPH densities.

Summarizing, a too large value for h_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.

5 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 |$e_{\mathrm{sf}} = 0.32\, \mathrm{per\, cent}$| and |$e_{\mathrm{sf}} = 1\, \mathrm{per\, cent}$| to illustrate the discussed issues. Higher values of e_sf mean that gas particles are converted into star particles on shorter time-scales for a given density. The isolated galaxy simulations from this work do not produce artificial clumps for |$e_{\mathrm{sf}}\gt 1\, \mathrm{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 |$e_{\mathrm{sf}}\ge 1\, \mathrm{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} > l_smooth) only depends on the particle mass (i.e. the mass resolution). Gas with densities and temperatures such that λ_{J, N} < l_smooth, is formally unresolved independently of the softening length, which for a baryon particle mass of |$m_{\mathrm{B}} = 10^5\, \mathrm{M}_{\odot }$| 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 |$\epsilon _{\mathrm{min}} = 2\, \mathrm{pc}$| and a simulation with the same mass resolution of |$m_{\mathrm{B}} = 10^5\, \mathrm{M}_{\odot }$| but constant softening might use |$\epsilon = 100\, \mathrm{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 = h_min) and potentially a larger value for ϵ for constant softening if one wishes to suppress artificial subkernel 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 (⁠|$m_{\mathrm{B}} = 10^5$| and |$8 \times 10^5\, \mathrm{M}_{\odot }$|⁠), 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, and hydro) also need to take into account the effect of artificial heating of the stellar disc by dark matter and stellar halo particles (see e.g. Ludlow et al. 2019a, 2021, 2023; Wilkinson et al. 2023), 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 idealized simulations.

The SFR 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 SFR 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 z = 0 limit the pressure (or temperature) of the ISM by an effective pressure floor |$P\propto \rho ^{\gamma _{\mathrm{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 |$n_{\mathrm{H}}\gt 0.1\, \mathrm{cm}^{-3}$| and in the widely used Springel & Hernquist (2003) model, the polytropic index varies with density between γ_eff ≈ 2.5 at |$0.1\, \mathrm{cm}^{-3}$| to γ_eff ≈ 4/3 at |$100\, \mathrm{cm}^{-3}$| (for z = 0; see fig. 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 (t_{ff, s}, equation 3) is longer than the Newtonian free-fall time (t_{ff, N}, equation 2) in simulations with an adaptive softening length. In simulations 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.

6 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 time-scale to the sound crossing time-scale (see the derivation in Appendix  A). In a homogeneous medium with gas density n_H and temperature T, 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 two-body scattering. In parallel to the Newtonian Jeans criteria, the softened Jeans length (mass), λ_{J, s} (M_{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 better 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 M_{J, N} by several orders of magnitude for gas densities and temperatures typical of the cold ISM. For example, at a gas temperature of a few tens of K and a gas density of |$\approx 100\, \mathrm{cm}^{-3}$|⁠, the Newtonian Jeans mass is |$M_{\mathrm{J,N}}\approx 10^2 \, \mathrm{M}_{\odot }$|⁠, and the softened Jeans mass for a constant softening length of |$\epsilon =100\, \mathrm{pc}$| is |$M_{\mathrm{J,s}}\approx 10^5 \, \mathrm{M}_{\odot }$| (Fig. 1).

In simulations with particle masses of |$\approx 10^5\, \mathrm{M}_{\odot }$| that do not impose an artificial pressure or entropy floor in the form of an effective EOS, cold neutral and molecular ISM gas is formally unresolved because gravitational instabilities should form within an individual smoothing kernel (λ_{J, N} < l_smooth). Modelling a multiphase ISM 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 (Section 3.1) and (ii) pressure that is smoothed on length-scales larger than those on which gravity is softened (Section 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 subkernel instabilities is demonstrated in Section 4 in simulations of isolated galaxies using the swift code with a constant softening length, ϵ. As outlined in Section 3.1, the density of gas clumps that are smaller than the minimum allowed smoothing length, h_min, are underestimated. Recalculating the densities of all gas particles with a vanishing value for h_min reveals gas clumps with several orders of magnitude higher gas densities (compare the top and bottom rows 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, m_B and a maximum expected density, n_{sim, max}:

We recommend to generously fulfill this condition and could not identify a benefit of a larger value for h_min. If the minimum values for softening and smoothing lengths are identical (ϵ_min = h_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 and the softening lengths are approaching zero.¹² Assuming h_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 subkernel gravitational instabilities: an adaptive softening length that follows the smoothing length of the kernel leads to artificial suppression of subkernel instabilities, while a constant softening length can result in artificially inducing subkernel 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 subkernel 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, l_smooth, for larger regions in density–temperature space (see the middle panel in Fig. 5). The impact of subkernel instabilities on galaxy properties remains to be tested, because cooling rates, SFRs, and other density- or pressure-dependent subgrid 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 |$m_{\mathrm{B}} = 8\times 10^5\, \mathrm{M}_{\odot }$|⁠, but the dependence on m_B is weak) do not depend on an effective pressure floor for numerical stability.¹³ Simulations with comparable mass resolutions and without an artificial pressure floor do not suffer from the numerical issue discussed in Section 3.1 if h_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, for example, 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 MCs, 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.

ACKNOWLEDGEMENTS

The authors would like to thank the anonymous referee for a very constructive referee report. SP acknowledges helpful discussions, especially with Vadim Semenov, at the Computational Galaxy Formation meeting, organized by Thorsten Naab and Volker Springel in Ringberg, Germany in 2022 April. This research was funded in part by the Austrian Science Fund (FWF) grant no. V 982-N. This paper made use of the following python packages: astropy (Astropy Collaboration 2022), numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), matplotlib (Hunter 2007), unyt (Goldbaum et al. 2018), and swiftsimio (Borrow & Borrisov 2020; Borrow & Kelly 2021). This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K00042X/1, ST/P002293/1, ST/R002371/1, and ST/S002502/1, Durham University, and STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure. The computational results presented have been achieved in part using the Vienna Scientific Cluster (VSC).

DATA AVAILABILITY

The simulations were run with the public SPH code swift (Schaller et al. 2023, www.swiftsim.com) version: 0.9.0, revision v0.9.0-1182-g423e9dd8. The initial conditions are discussed in Nobels et al. (2023) and public within the ‘IsolatedGalaxy’ example within swift (here the ‘M5_disk.hdf5’ and ‘M6_disk.hdf5’ files are used). The routines to set up, run, and analyse the simulations, including the routines to create each figure, are public on gitlab (https://gitlab.phaidra.org/softenedjeanscriteria). The results are therefore fully reproducible, and only require limited computational resources.

Footnotes

References

APPENDIX A: DERIVATION OF SOFTENED JEANS CRITERIA

In the following, we re-derive the Jeans length for softened gravity. For simplicity, we focus on comparing the free-fall time, t_ff, to the sound crossing time, t_sc, of a self-gravitating structure with an initial radius, R. This derivation can easily be extended to general gravitational accelerations.

A system is in equilibrium if it has a radius R for which the free-fall time equals the sound-crossing time (R ≡ λ_J, the Jeans length). If the structure (here: a gas cloud) has a radius larger than the Jeans length, the cloud collapses or fragments, otherwise it expands. The sound crossing time is defined as

with the sound speed

and t_sc is independent of the shape of the gravitational potential. Here, γ is ratio of specific heats, k_B is the Boltzmann constant, T is the gas temperature, μ the mean particle mass, and m_H the hydrogen particle mass.

For a general definition of the free-fall time, we integrate the time it takes a test mass to fall from rest at a distance R to the centre of the potential:

Integrating the equation of motion

from R to r leads to

which is used in equation (A3) to calculate the free-fall time

A pressure-less collapse and therefore no shell-crossing is assumed, which means that the mass inside the radius of the test mass remains constant throughout the free-fall: M(< r) = M = 4πR³/3ρ. Finally, the Jeans length is calculated by solving the equation t_ff(R) = t_sc(R, c_s) for R.

A1 Newtonian (unsoftened) gravity

The Jeans length for Newtonian gravity can be found in many textbooks but depending on the derivation the pre-factors are slightly different. We repeat one of the textbook derivations here briefly to allow for a direct comparison of the individual steps with the derivation of the softened Jeans length and to demonstrate which pre-factors are included in this derivation.¹⁴

For Newtonian gravity,

and the velocity for the free-fall is

The free-fall time is

or

for the gas density ρ = n_Hm_H/X_H with the hydrogen number density, n_H, the hydrogen mass fraction X_H (here we use X_H = 0.74 for solar chemical composition, Asplund et al. 2009). The Jeans length for Newtonian gravity is

which corresponds to a Jeans mass of

where we use γ = 5/3 and μ = 1.28, values typical for a neutral atomic ideas gas with solar metallicity. For diatomic gas, such as molecular hydrogen, the ratio of specific heats would be γ = 1.4 and μ = 2, but these are comparable to other order of unity effects, such as the concrete particle distribution.

A2 Plummer softening

For a Plummer (1911) softened potential with

and a softening scale ϵ, the free-fall velocity, equation (A5), becomes

and the free-fall time (equation A6) is

This expression can be re-written as

with the Newtonian free fall time t_{ff, N} and the dimensionless parameters x = r/R and f_P = ϵ/R. This integral cannot be solved analytically but it can be numerically evaluated for different values of f_P. The following fit matches these data points to within a few per cent:

The free-fall in a Plummer softened potential is therefore slowed down by a factor of ≈2.0 if the softening scale equals the size of the clump (f_P = 1) and by a factor ≈45 for a softening scale that is 10 times the size of the clump (f_P = 10).

For the limiting case f_P ≫ 1 (i.e. softening is important), the Plummer free-fall time can be calculated analytically (see Nobels et al. 2023) and for f_P ≪ 1 the Plummer free-fall time approaches the Newtonian free-fall time. For the full expression of the Jeans length (solving t_{ff, P, fit} = R/c_s for R with λ_{J, N} = t_{ff, N}c_s, the Newtonian Jeans length), we use a root finding algorithm for

This dimensionless equation is solved numerically for values of ϵ/λ_{J, N} between 0.01 and 10⁵. The results are fit by

and the softened Jeans for the Plummer potential is fit by

The top left panel in Fig. A1 illustrates the Jeans length (black, left y-axis) and Jeans mass (red, right y-axis) as a function of the ratio between the gravitational softening length, ϵ, and the Newtonian Jeans length λ_{J, N}. For small values of ϵ/λ_{J, N} (i.e. the gas cloud is much larger than the softening scale), the softened Jeans length approaches the Newtonian Jeans length (log λ_{J, P}/λ_{J, N} → 0) but for large values of ϵ/λ_{J, N} the softened Jeans length and the softened Jeans mass are much larger than their Newtonian counterparts.

The bottom left panel shows the ratio between the exact (numerical) solution and the fits from equations (A21) and (A22) which differ by up to 20 per cent for the Jeans mass. This is small compared to other order of unity effects (i.e. different derivations, non-spherical particle distribution).

A3 Wendland C2 softening

For the Wendland (1995) C2 kernel, the acceleration |a| = |a_W| is defined in two parts, a_W(r < H, H) and a_W(r ≥ H) = a_N(r). This means that the free-fall velocity has distinct formulations depending on the relative size of the cloud, R, and the softening scale H = 3ϵ, where ϵ is the Plummer equivalent softening length.

If the size of the cloud is smaller than H (f_W ≡ H/R ≥ 1), the gravitational acceleration is softened throughout the free-fall and the velocity is calculated from

with the Wendland C2 kernel

where u = r/H. If the size of the cloud is larger than the softening length, the free-fall velocity is Newtonian for r > H and softened for r ≤ H (or: u ≡ r/H ≤ 1). The full integral for the velocity for r < H consists therefore of two parts:

and the free-fall velocity is Newtonian (equation A8) for u > 1. Combining all these cases results in the following expression for the free-fall velocity in a Wendland C2 softened gravitational potential:

These expressions have been numerically verified by an explicit leapfrog integration with a small, constant time-step of 0.1 Myr.

A3.1 Case 1: free-fall always softened (f_W ≥ 1 or R ≤ H)

The free-fall time is

which can be expressed relative to the Newtonian free-fall time t_{ff, N} and the dimensionless variables, x = r/R and f_W = H/R as

The dimensionless integral and therefore t_{ff, W}/f_{ff, N} is solved numerically for values of f_W between 1 and 10³.

A3.2 Case 2: free-fall partially softened (f_W < 1 or R > H)

If the size, R, of the self-gravitating structure is larger than the softening length, H, the free-fall time consists of two parts and is defined as

with the Newtonian free-fall velocity v_N. The first integral can be solved analytically

with C_W = (1/f_W − 1)^1/2. The second integral is rewritten as

The free-fall time for a partially softened trajectory in a Wendland C2 potential is

and is solved numerically for f_W between 10⁻³ and 1.

A3.3 Combined

The numerically calculated values for t_{ff, W}/t_{ff, N} are fit with

The Jeans length for a Wendland C2 softened gravitational potential is calculated by solving t_{ff, W, fit} = R/c_s for R with λ_{J, N} = t_{ff, N}c_s, the Newtonian Jeans length. We use a root finding algorithm for

and this dimensionless equation is fit with

which leads us the softened Jeans mass

The expressions from equations (A36) and (A37) as well as the accuracy of the fits are shown in the right panels of Fig. A1. The behaviour is very similar to the softened Jeans criteria for the Plummer potential (left panels of Fig. A1) with a slight offset.

APPENDIX B: APPLICATION TO HIGHER RESOLUTION SIMULATIONS

Throughout this work, we have discussed the derived instability criteria for large-scale simulations with baryon particle masses of |$\gtrsim 10^{5}\, \mathrm{M}_{\odot }$|⁠. Here, we briefly demonstrate how the discussion in Section 3 applies to simulations with much higher mass resolution.

Assuming a vanishing value for the minimum smoothing length, subkernel clumping may be of concern in simulations with a constant softening length. Fig. B1 repeats the analysis from Fig. 7 for simulations with a particle mass of |$m_{\mathrm{B}} = 4\, \mathrm{M}_{\odot }$| and constant softening lengths of l_soft = 0.15, 0.75, and |$1.5\, \mathrm{pc}$| (ϵ = 0.1, 0.5, and |$1\, \mathrm{pc}$|⁠). Gas with densities and temperatures below the black lines is gravitationally unstable in the respective simulation to perturbations with a length-scale of l_smooth, the size of an individual kernel and the spatial resolution limit of the hydrodynamic solver. In the low-density part (purple-shaded area, defined as l_soft < l_smooth) instabilities within the kernel may grow (see Section 3.2.1 and 3.2.4 for a detailed discussion). For |$m_{\mathrm{B}} = 4\, \mathrm{M}_{\odot }$| and |$l_{\mathrm{soft}} = 0.15\, \mathrm{pc}$| (e.g. Hislop et al. 2022; Steinwandel et al. 2023), subkernel instabilities may form over the full density and temperature range typical for molecular gas in the Milky Way (indicated as dark patch at 10² ≲ n_H[cm⁻³] ≲ 10⁶, |$T\approx 20\, \mathrm{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\, \mathrm{pc}$|⁠, leads to numerical issues caused by runaway stars with very high velocities. For the same particle mass (⁠|$m_{\mathrm{B}} = 4\, \mathrm{M}_{\odot }$|⁠), Hu et al. (2017) varied the softening length between |$\epsilon = 0.5$| and |$2\, \mathrm{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 SFR only has a weak dependence on ϵ. It would be interesting to see if the discussed subkernel 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, subkernel clumping can be avoided for simulations with |$m_{\mathrm{B}} = 4\, \mathrm{M}_{\odot }$|⁠, for softening lengths of |$\epsilon \ge 2\, \mathrm{pc}$|⁠, as, for example, used in Hu et al. (2016, 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 (Appendix C1) and the thermal state of the ISM (Appendix 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 |$E_{\mathrm{SN}} = 10^{51}\, \mathrm{erg}$|⁠. This energy is stochastically injected into the nearest gas particle by increasing its temperature by |$\Delta T_{\mathrm{heat}} = 10^{7.5}\, \mathrm{K}$|⁠. Dalla Vecchia & Schaye (2012) estimate the maximum gas density up to which thermal feedback is efficient by comparing the gas cooling time, t_c, to the sound-crossing time, t_sc, across a heated resolution element. If t_c ≫ t_sc (e.g. by a factor of f_t ≡ t_c/f_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 |$\approx 10\, \mathrm{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 |$\Delta T_{\mathrm{heat}} = 10^{7.5}\, \mathrm{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 |$E_{\mathrm{SN}} = 0.3\times 10^{51}$| and |$3 \times 10^{51}\, \mathrm{erg}$| depending on metallicity (higher E_SN for lower metallicity) and gas density at the time of star formation (higher E_SN for higher gas density).

Increasing E_SN (⁠|$10^{51}\, \mathrm{erg}$|⁠, |$2 \times 10^{51}\, \mathrm{erg}$|⁠, |$4 \times 10^{51}\, \mathrm{erg}$|⁠, columns 1–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: |$h_{\mathrm{min}} = 77.5\, \mathrm{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 E_SN reduces the highest density gas within the galaxy and for the largest energy per supernova (⁠|$E_{\mathrm{SN}} = 4\times 10^{51}\, \mathrm{erg}$|⁠, third column), barely any molecular gas (with |$n_{\mathrm{H}} \gt 100\, \mathrm{cm}^{-3}$|⁠) is left within the galaxy. On the other hand, a smaller value of h_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 ISRF and CR rate (see Ploeckinger & Schaye 2020 for details) and a redshift-dependent 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 (⁠|$m_{\mathrm{B}} = 10^5\, \mathrm{M}_{\odot }$| and |$\epsilon = 250\, \mathrm{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, h_min: (i) stellar clumps and (ii) shock-heated, high-density (⁠|$n_{\mathrm{H}}\gtrsim 100\, \mathrm{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 h_min is large (⁠|$77.5\, \mathrm{pc}$|⁠, columns 1–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 h_min to |$7.75\, \mathrm{pc}$| (columns 4–6), confirming the findings from Section 4.