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

Sylvia Ploeckinger,1,2 Folkert S.J. Nobels,3 Matthieu Schaller,2,3 Joop
Schaye3
1Department of Astrophysics, University of Vienna, Türkenschanzstrasse 17,
1180 Vienna, Austria
2Lorentz Institute for Theoretical Physics, Leiden University, PO Box 9506,
2300 RA Leiden, the Netherlands
3Leiden Observatory, Leiden University, PO Box 9513, NL-2300 RA Leiden, The
Netherlands E-mail<EMAIL_ADDRESS>0000-0002-1965-1650
0000-0002-2395-4902 0000-0002-0668-5560

(Accepted XXX. Received YYY; in original form ZZZ)

###### 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 over-
smoothing of the hydrodynamical properties relative to the gravitational force
resolution. This effect is illustrated using simulations of isolated disk
galaxies with the smoothed particle hydrodynamics code Swift. We also
demonstrate how to avoid the formation of artificial clumps in gas and stars
by adjusting the gravitational and hydrodynamical force resolutions.

###### keywords:

hydrodynamics – instabilities – methods: numerical – galaxies: ISM – galaxies:
general – galaxies: formation

††pubyear: 2023††pagerange: Resolution criteria to avoid artificial clumping
in Lagrangian hydrodynamic simulations with a multi-phase interstellar
medium–A.3

##  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_{\mathrm{J}}$ as the limiting wavenumber between oscillatory solutions
($k>k_{\mathrm{J}}$, short wavelengths) and exponentially growing solutions
($k<k_{\mathrm{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 $\lambda_{\mathrm{J}}$ above
which density perturbations (or: gas clouds) of size $r$, sound speed
$c_{\mathrm{s}}$, and (unperturbed) density $\rho$ are gravitationally
unstable:

$\lambda_{\mathrm{J}}=A\left(\frac{c_{\mathrm{s}}^{2}}{G\rho}\right)^{1/2}\;,$
(1)

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 2-body
interactions. A classical prescription for softened gravity is the Plummer
(1911) potential $\phi_{\mathrm{P}}=-GM(r^{2}+\epsilon^{2})^{-1/2}$ for a
point mass, $M$, and the gravitational softening length, $\epsilon$. The
optimal choice of $\epsilon$ 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; Ludlow et al., 2019b)
and is further complicated in hydrodynamical simulations (e.g. Price &
Monaghan, 2007; Ludlow et al., 2020). The softening length, $\epsilon$, is a
measure for the gravitational force resolution of a simulation, because
gravity is non-Newtonian (i.e. softened) for $r\lesssim\epsilon$.

A common technique to solve the equations of motions for a collisional fluid
is Smoothed Particle Hydrodynamics (SPH, originally from Lucy, 1977; Gingold &
Monaghan, 1977 but see e.g. Borrow et al., 2022 for a modern example) 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.

Targeting a constant number of neighbours in a kernel, $h$ decreases with
increasing gas density ($h\propto\rho^{-1/3}$). Depending on the code, the
gravitational and hydrodynamical force resolutions can be coupled and
$\epsilon\propto h$ (”adaptive softening") or $\epsilon$ can be set to a
constant value for all densities. Adaptive softening is implemented e.g. in
the SPH codes Gadget4 (Springel et al., 2021), Phantom (Price et al., 2018),
and the SPH and meshless finite mass / volume code Gizmo (Hopkins, 2015).
Their adaptive softening lengths are defined such that they are tied to the
kernel length ($\approx 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,
$\epsilon$, 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; Schaller et al., 2023) simulations.

In the moving mesh code Arepo (Weinberger et al., 2020), the softening length
is adaptive and follows
$\epsilon_{\mathrm{Arepo}}=f_{\mathrm{h}}(3V/(4\pi))^{1/3}$, where
$f_{\mathrm{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
IllustrisTNG project (Weinberger et al., 2017; Pillepich et al., 2018a), which
uses Arepo, the softening length, $\epsilon$ is set to a minimum value
$\epsilon_{\mathrm{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 ($\epsilon=\epsilon_{\mathrm{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 ($\epsilon=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_{\mathrm{neigh}}$ particles (later reduced to $1.5N_{\mathrm{neigh}}$ by
Bate et al., 2003), where $N_{\mathrm{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 $\epsilon<h$ and
artificially inhibited for $\epsilon>h$) and not by the physical conditions.

Hubber et al. (2006) introduced the “Jeans test”, a numerical test of a plane-
wave perturbation within a homogeneous medium, and found that if the
resolution criterion from Bate & Burkert (1997) is violated, the collapse of
marginally unstable modes is suppressed but artificial fragmentation does not
occur in SPH for their setup with an adaptive softening length $\epsilon=h$,
confirming the analytical results from Whitworth (1998). Yamamoto et al.
(2021) repeated the Jeans test with various hydrodynamic schemes within the
Gizmo code (Hopkins, 2015) (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 $\epsilon=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
$\epsilon=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}}>0.1\,\mathrm{cm}^{-3}$ formally solves this issue. For an
effective polytropic index of $\gamma_{\mathrm{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_{\mathrm{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
normalisations 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 disks are less compact and less gravitationally bound
with a higher artificial pressure floor. An artificially increased Jeans mass
in the form of a polytropic EOS, a pressure or entropy floor, or the Springel
& Hernquist (2003) subgrid model111Hydrodynamically, the Springel & Hernquist
(2003) model acts as a polytropic EOS with a density-dependent
$\gamma_{\mathrm{eff}}$, see their figure 1. is implemented in almost222One
recent exception is the FIREbox project (Feldmann et al., 2023). all
simulations of cosmologically representative volumes that reach $z=0$.

Table 1: Overview of equations for Newtonian gravity, a Plummer softened
gravitational potential, and a Wendland C2 softened potential. Newtonian (no
softening)
---
Potential | $\varphi_{\mathrm{N}}(r)$ | = | $-GMr^{-1}$
Acceleration | $|a_{\mathrm{N}}(r)|$ | = | $GMr^{-2}$
Free-fall time | $t_{\mathrm{ff,N}}$ | = | $\left(\frac{3\pi}{32G\rho}\right)^{1/2}=4.4\,\mathrm{Myr}\left(\frac{n_{\mathrm{H}}}{100\,\mathrm{cm}^{-3}}\right)^{-1/2}$
Jeans length | $\lambda_{\mathrm{J,N}}$ | = | $\left(\frac{3\pi\gamma X_{\mathrm{H}}k_{\mathrm{B}}T}{32Gm_{\mathrm{H}}^{2}n_{\mathrm{H}}}\right)^{1/2}=1.5\,\mathrm{pc}\left(\frac{T}{10\,\mathrm{K}}\right)^{1/2}\left(\frac{n_{\mathrm{H}}}{100\,\mathrm{cm}^{-3}}\right)^{-1/2}$
Jeans mass | $M_{\mathrm{J,N}}$ | = | $\frac{4\pi\rho}{3}\lambda_{\mathrm{J,N}}^{3}=46\,\mathrm{M}_{\odot}\,\left(\frac{T}{10\,\mathrm{K}}\right)^{3/2}\left(\frac{n_{\mathrm{H}}}{100\,\mathrm{cm}^{-3}}\right)^{-1/2}$
Plummer softening with softening scale $\epsilon=\epsilon_{\mathrm{Plummer}}$
Potential | $\varphi_{\mathrm{P}}(r,\epsilon)$ | = | $-GM\left(r^{2}+\epsilon^{2}\right)^{-1/2}$
Acceleration | $|a_{\mathrm{P}}(r,\epsilon)|$ | = | $GMr\left(r^{2}+\epsilon^{2}\right)^{-3/2}$
Free-fall time | $t_{\mathrm{ff,P,fit}}$ | = | $t_{\mathrm{ff,N}}\left(1+2^{2/3}\left(\frac{\epsilon}{R}\right)^{2}\right)^{3/4}$
Jeans length | $\lambda_{\mathrm{J,P,fit}}$ | = | $\lambda_{\mathrm{J,N}}\left(1+1.42\left(\frac{\epsilon}{\lambda_{\mathrm{J,N}}}\right)^{3/2}\right)^{2/5}$
Jeans mass | $M_{\mathrm{J,P,fit}}$ | = | $M_{\mathrm{J,N}}\left(1+1.42\left(\frac{\epsilon}{\lambda_{\mathrm{J,N}}}\right)^{3/2}\right)^{6/5}$
Wendland C2 / Swift softening with softening scale
$H=3\epsilon=\epsilon_{\mathrm{Plummer}}$ and $u=r/H$
Potential | $\varphi_{\mathrm{W}}(r<H,H)$ | = | $-GMH^{-1}W(u)$
|  with $W(u)$ | = | $\left(-3u^{7}+15u^{6}-28u^{5}+21u^{4}-7u^{2}+3\right)$
| $\varphi_{\mathrm{W}}(r\geq H)$ | = | $-GMr^{-1}$
Acceleration | $|a_{\mathrm{W}}(r<H,H)|$ | = | $GMrH^{-3}V(u)$
|  with $V(u)$ | = |  $-W^{\prime}(u)/u$ = $\left(21u^{5}-90u^{4}+140u^{3}-84u^{2}+14\right)$
| $|a_{\mathrm{W}}(r\geq H)|$ | = | $GMr^{-2}$
Free-fall time | $t_{\mathrm{ff,W,fit}}$ | = |  $t_{\mathrm{ff,N}}\left(1+\frac{1}{7}\left(\frac{H}{R}\right)^{3}\right)^{1/2}$ = $t_{\mathrm{ff,N}}\left(1+\frac{27}{7}\left(\frac{\epsilon}{R}\right)^{3}\right)^{1/2}$
Jeans length | $\lambda_{\mathrm{J,W,fit}}$ | = |  $\lambda_{\mathrm{J,N}}\left(1+0.27\left(\frac{H}{\lambda_{\mathrm{J,N}}}\right)^{2}\right)^{3/10}$ = $\lambda_{\mathrm{J,N}}\left(1+2.43\left(\frac{\epsilon}{\lambda_{\mathrm{J,N}}}\right)^{2}\right)^{3/10}$
Jeans mass | $M_{\mathrm{J,W,fit}}$ | = |  $M_{\mathrm{J,N}}\left(1+0.27\left(\frac{H}{\lambda_{\mathrm{J,N}}}\right)^{2}\right)^{9/10}$ = $M_{\mathrm{J,N}}\left(1+2.43\left(\frac{\epsilon}{\lambda_{\mathrm{J,N}}}\right)^{2}\right)^{9/10}$

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 multi-
phase ISM (i.e. hot ionized, warm ionized / neutral, cold neutral, as well as
molecular gas) to form. Based on the Jeans mass arguments above, this would
only be possible for simulations with a baryon mass resolution of $\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 multi-phase ISM without
an artificial effective pressure floor.

This paper is organized as follows. In section 2 we introduce the “softened”
Jeans criteria (softened Jeans length and softened Jeans mass) and explain why
they are more appropriate for describing the conditions of gravitational
instabilities in a simulation with softened gravity than the physical Jeans
criteria that assume Newtonian gravity. We discuss in section 3 for which
densities and temperatures gravitational instabilities are modelled physically
correctly, or are numerically induced or suppressed, for both constant and
adaptive softening. This section includes an independent criterion on the
minimum smoothing length, $h_{\mathrm{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 $\epsilon$ and $h_{\mathrm{min}}$
is illustrated using simulations of isolated disk galaxies with the modern
hydrodynamics code Swift333Swift is publicly available at: www.swiftsim.com
(Schaller et al., 2018; Schaller et al., 2023) in section 4. The implications
for (cosmological) galaxy formation simulations are discussed in section 5 and
we summarize the results in section 6.

In the literature, the terms “smoothing” and “softening” are each sometimes
used for either the hydrodynamical or gravitational force calculations
involving a “smoothing” / “softening” kernel. Throughout this work we strictly
use “smoothing” when referring to the calculation of hydrodynamical quantities
and “softening” for calculations of gravitational forces. We use $\log$ as
$\log_{\mathrm{10}}$ 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 2-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 $\epsilon$ in Lagrangian codes. The Newtonian gravitational
acceleration, $a_{\mathrm{N}}\propto r^{-2}$ is replaced by the softened
gravitational acceleration, for example $a_{\mathrm{P}}\propto
r(r^{2}+\epsilon^{2})^{-3/2}$ for a Plummer (1911) potential (see Table 1 for
an overview). In the limit of large separations, $r\gg\epsilon$, the softened
acceleration equals the Newtonian acceleration ($a_{\mathrm{P}}\rightarrow
a_{\mathrm{N}}\propto r^{-2}$) but for $r\ll\epsilon$, the softened
acceleration approaches $a_{\mathrm{P}}\propto r$ and therefore diverges from
the Newtonian gravity as $r\rightarrow 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 gravity444For the derivation of the softened Jeans criteria we
assume that the hydrodynamic forces are calculated accurately. The impact of
unresolved hydrodynamic forces on numerical instabilities is discussed in
detail in section 3.. 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_{\mathrm{ff,s}}=t_{\mathrm{ff,W,fit}}$
from equation 49)

$t_{\mathrm{ff,s}}=t_{\mathrm{ff,N}}\left(1+\frac{27}{7}\left(\frac{\epsilon}{R}\right)^{3}\right)^{1/2}$
(2)

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

$t_{\mathrm{ff,s}}\propto
t_{\mathrm{ff,N}}\left(\frac{\epsilon}{R}\right)^{3/2}$ (3)

for $\epsilon\gg 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$.

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

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

Fig. 1 shows an overview of the Jeans masses from equations (A.1), (37), and
(52) 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 $\epsilon$, 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_{\mathrm{J,s}}=M_{\mathrm{J,P,fit}}$
(Plummer) and $M_{\mathrm{J,s}}=M_{\mathrm{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 (16) with a more general form of a signal crossing time,
$t_{\mathrm{sig}}$, that can also include a turbulent or a magnetic component
(see Nobels et al., 2023). The softened Jeans mass and length are defined
relative to the Newtonian Jeans mass and length. For example,

$\displaystyle\lambda_{\mathrm{J,s}}$
$\displaystyle=\lambda_{\mathrm{J,N}}\left(1+2.43\left(\frac{\epsilon}{\lambda_{\mathrm{J,N}}}\right)^{2}\right)^{3/10}$
(4) $\displaystyle M_{\mathrm{J,s}}$
$\displaystyle=M_{\mathrm{J,N}}\left(1+2.43\left(\frac{\epsilon}{\lambda_{\mathrm{J,N}}}\right)^{2}\right)^{9/10}$
(5)

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 $\epsilon$ but is also applicable for simulations with
adaptive softening when $\epsilon$ is limited by a minimum value
$\epsilon_{\mathrm{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 is approximated by subgrid prescriptions, e.g.
for star formation. For code stability it is typically preferable to
numerically suppress the collapse of gas clouds than to numerically induce
collapse because the latter can lead to prohibitively small timesteps, and
artificial clumps in the stellar distribution.

In this section we discuss two distinct instabilities that can occur in
Lagrangian simulations and how they depend on the gravitational softening
length (i.e. the gravitational force resolution), the hydrodynamic smoothing
length (i.e. the hydrodynamical force resolution) and the mass resolution. In
subsection 3.1 we identify the problematic behaviour of SPH-like simulation
codes when the hydrodynamical smoothing length $h$ is limited by a too large
minimum value, $h_{\mathrm{min}}$. The resulting instability is caused by an
underestimate of the gas density and leads to artificial clumps of closely
spaced particles. Setting $h_{\mathrm{min}}$ to a very small non-zero value
resolves this issue and we cannot identify any disadvantages of letting
$h_{\mathrm{min}}\rightarrow 0$ (to be discussed in detail in section 4).

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

###  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
$\lambda_{\mathrm{J,N}}<\epsilon$ (Fig. 1). Gas clumps would therefore be
increasingly stabilised against gravitational collapse as their density gets
larger. However, this assumes that the hydrodynamic forces are modelled
accurately. For simulations with a minimum smoothing length,
$h_{\mathrm{min}}\neq 0$, the hydrodynamic forces can be over-smoothed for
high gas densities and we explain below how this can cause a numerical runaway
collapse.

In this section, the smoothing length, $h$, and its lower limit,
$h_{\mathrm{min}}$, are defined as in the hydrodynamics solver Sphenix (Borrow
et al., 2022), implemented in the Swift code555We discuss in section 3.2 how
this can be adapted to other simulations codes by defining a general smoothing
length scale, $l_{\mathrm{smooth}}$.. In Swift, the smoothing length is based
on the kernel standard deviation as in Dehnen & Aly (2012) and is calculated
as

$\displaystyle h$
$\displaystyle=\eta_{\mathrm{res}}\left(\frac{X_{\mathrm{H}}m_{\mathrm{B}}}{m_{\mathrm{H}}n_{\mathrm{H}}}\right)^{1/3}$
$\displaystyle=38.3\,\mathrm{pc}\left(\frac{m_{\mathrm{B}}}{10^{5}\,\mathrm{M}_{\odot}}\right)^{1/3}\left(\frac{n_{\mathrm{H}}}{100\,\mathrm{cm}^{-3}}\right)^{-1/3}$
(7)

where $\eta_{\mathrm{res}}$ is a constant related to the number of neighbours
in the kernel, $X_{\mathrm{H}}$ is the hydrogen mass fraction,
$m_{\mathrm{B}}$ the baryon particle mass, $m_{\mathrm{H}}$ the hydrogen
particle mass, $n_{\mathrm{H}}$ and the hydrogen number density. We use the
typical values $\eta_{\mathrm{res}}=1.2348$ (this corresponds to 58 neighbours
with a Wendland C2 kernel) and $X_{\mathrm{H}}=0.74$.

The minimum smoothing length in Swift is set by the dimensionless parameter
$h_{\mathrm{min,ratio}}$ which is defined as the ratio between the kernel
support, $h\gamma_{\mathrm{k}}$, and the length scale above which gravity is
Newtonian, $H=3\epsilon$. For the Wendland C2 kernel
$\gamma_{\mathrm{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 over-smoothed is

$\displaystyle n_{\mathrm{H,hmin}}$
$\displaystyle=\frac{X_{\mathrm{H}}}{m_{\mathrm{H}}}\eta_{\mathrm{res}}^{3}\left(\frac{m_{\mathrm{B}}}{h_{\mathrm{min}}^{3}}\right)$
$\displaystyle=5.6\,\mathrm{cm}^{-3}\left(\frac{m_{\mathrm{B}}}{10^{5}\,\mathrm{M}_{\odot}}\right)\left(\frac{h_{\mathrm{min}}}{100\,\mathrm{pc}}\right)^{-3}\;.$
(8)

Note that the critical density $n_{\mathrm{H,hmin}}$ can remain constant for
simulations with different mass resolutions $m_{\mathrm{B}}$ if
$h_{\mathrm{min}}$ is adjusted accordingly.

The softened Jeans criteria derived in appendix A are therefore only valid for
$n_{\mathrm{H}}\leq n_{\mathrm{H,hmin}}$. For higher densities the
hydrodynamic forces can become inaccurate due to over-smoothing and this
inaccuracy depends on the density contrast within $h_{\mathrm{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_{\mathrm{H}}$ which increases from $\log
n_{\mathrm{H}}[\mathrm{cm}^{-3}]=2$ (leftmost panel) to $\log
n_{\mathrm{H}}[\mathrm{cm}^{-3}]=6$ (rightmost panel). The red solid circle
represents the smoothing length for this gas density and particle mass
(equation 7) and the blue dashed circle represents the minimum smoothing
length which is set to $h_{\mathrm{min}}=15.5\,\mathrm{pc}$ (corresponding to
e.g. $\epsilon=100\,\mathrm{pc}$ and $h_{\mathrm{min,ratio}}=0.1$) in this
example.

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

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

To avoid a potential run-away collapse, which is desirable for both numerical
and physical reasons, we can define a problematic region in the density-
temperature phase-space based on the conditions discussed above:

  1. i

$h<h_{\mathrm{min}}$ (or: $n_{\mathrm{H}}>n_{\mathrm{H,hmin}}$)

  2. ii

fragmentation on scales below $h_{\mathrm{min}}$

For an estimate of the clumping that is expected in the simulation, we use the
softened Jeans length and (ii) then translates to
$\lambda_{\mathrm{J,s}}<h_{\mathrm{min}}$.

Assuming that $h_{\mathrm{min}}<\epsilon$ (to avoid inducing numerical
collapse, see e.g. Bate & Burkert, 1997), the condition
$\lambda_{\mathrm{J,s}}<h_{\mathrm{min}}$ falls into the softened part of the
phase-space where $H\equiv 3\epsilon\gg\lambda_{\mathrm{J,N}}$ and we
therefore approximate the softened Jeans mass (here for the Wendland C2
kernel, equation 51) 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 run-away collapse correspond to the following regions in
the phase-space diagram as

$\displaystyle n_{\mathrm{H}}$
$\displaystyle>5.6\,\mathrm{cm}^{-3}\left(\frac{m_{\mathrm{B}}}{10^{5}\,\mathrm{M}_{\odot}}\right)\left(\frac{h_{\mathrm{min}}}{100\,\mathrm{pc}}\right)^{-3}\;,$
(9) $\displaystyle T$
$\displaystyle<27\,\mathrm{K}\left(\frac{h_{\mathrm{min}}}{100\,\mathrm{pc}}\right)^{5}\left(\frac{\epsilon}{700\,\mathrm{pc}}\right)^{-3}\left(\frac{n_{\mathrm{H}}}{100\,\mathrm{cm}^{-3}}\right)\;.$
(10)

This “runaway collapse zone” is indicated in Fig. 3 for various combinations
of the particle mass $m_{\mathrm{B}}$, Plummer equivalent softening length
$\epsilon$, and the minimum smoothing length $h_{\mathrm{min}}$. The default
values (thick, black lines) with $\log
m_{\mathrm{B}}[\mathrm{M}_{\odot}]=6.3$, $\epsilon=700\,\mathrm{pc}$, and
$h_{\mathrm{min}}=108.5\,\mathrm{pc}$ (i.e. $h_{\mathrm{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_{\mathrm{B}}$, middle panel: $\epsilon$, bottom panel:
$h_{\mathrm{min}}$).

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

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

If the highest (expected) density in a simulation is $n_{\mathrm{sim,max}}$,
the minimum smoothing length should be set so that
$h(n_{\mathrm{sim,max}})\geq h_{\mathrm{min}}$. With $h$ from equation (7),
this means that the minimum smoothing length should follow

$h_{\mathrm{min}}<8.3\,\mathrm{pc}\left(\frac{m_{\mathrm{B}}}{10^{5}\,\mathrm{M}_{\odot}}\right)^{1/3}\left(\frac{n_{\mathrm{sim,max}}}{10^{4}\,\mathrm{cm}^{-3}}\right)^{-1/3}\;,$
(11)

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_{\mathrm{min}}$ and therefore recommend to use a very small but non-
zero666A non-zero value for $l_{\mathrm{smooth,min}}$ ensures that the
simulation does not break down in the highly unlikely, but not impossible,
case of a vanishing smoothing length for individual particles. value for
$h_{\mathrm{min}}$ that generously fulfills the conservative criterion from
equation (11).

Equation 11 does not depend on $\epsilon$ and applies to both simulations with
adaptive and constant softening lengths.

####  3.1.1 Examples from the literature

The resolution parameters used in both eagle flagship runs (eagle, eagle-high-
res, Schaye et al., 2015, see Table 2) 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 equation of state where the effective temperature777The polytropic
equation of state is in practice implemented as an entropy or internal energy
floor. The exact effective temperature therefore depends on the mean particle
mass, $\mu$, but this is irrelevant here as we show $T_{\mathrm{eff}}$ only
for reference. is

$T_{\mathrm{eff}}\approx
8000\,\mathrm{K}\left(\frac{n_{\mathrm{H}}}{0.1\,\mathrm{cm}^{-3}}\right)^{\gamma_{\mathrm{eff}}-1}$
(12)

with the polytropic index $\gamma_{\mathrm{eff}}=4/3$. This effective
temperature is indicated for reference as thick grey line in Fig. 3.

Fig. 4 illustrates why a pressure floor was necessary for a selection of
cosmological simulations with constant softening and which minor changes in
the resolution parameters, $\epsilon$ and $h_{\mathrm{min}}$, avoid numerical
runaway collapse in simulations with the same particle mass but without a
polytropic equation of state (bottom panel). The first and second panels shows
the runaway collapse zones for the resolution parameters of the two flaghship
simulations within the eagle project (“eagle”, “eagle-high-res”). The third
panel uses the resolution parameters from the highest resolution (“L1”)
simulations within the apostle project (Fattahi et al., 2016; Sawala et al.,
2016). 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. The effective temperature floor
was therefore necessary for their values for $h_{\mathrm{min}}$, despite the
low particle mass. The resolution parameters for all simulations are listed in
Table 2.

Table 2: Values for the gas particle mass ($m_{\mathrm{B}}$), the
gravitational force softening length ($\epsilon$), and the minimum
hydrodynamical smoothing length ($h_{\mathrm{min}}$) for a few simulations
with constant softening from the literature and an alternative set of
resolution parameter proposed in this work as shown in Fig. 4.

Simulation | $m_{\mathrm{B}}[\mathrm{M}_{\odot}]$ |  $\epsilon$ [pc] |  $h_{\mathrm{min}}$ [pc] | Reference
---|---|---|---|---
eagle | $1.81\times 10^{6}$ | 700 | 108.50 | (1)
eagle-high-res | $2.26\times 10^{5}$ | 350 | 54.25 | (1)
apostle-L1 | $1.00\times 10^{4}$ | 135 | 20.77 | (2,3)
Alternative |  $\approx 10^{5}$ to $10^{6}$ | 200 | 2.00 | -

(1) Schaye et al. (2015); (2) Fattahi et al. (2016); (3) Sawala et al. (2016)

Figure 4: Overview of the density-temperature phase-space for the eagle,
eagle-high-res, and apostle-L1 resolution parameters, as well as for an
alternative set (from top to bottom, see also Table 2). The contours indicate
the softened Jeans mass in units of $\log
M_{\mathrm{J,s}}[\mathrm{M}_{\odot}]$. The red shaded region at high densities
is the runaway collapse zone and SPH-estimated gas densities within this zone
may be underestimated (see section 3.1). The grey thick line is the effective
temperature floor used in eagle, eagle-high-res, and Apostle-L1. Typical
densities and temperatures for the warm neutral medium, the cold neutral
medium and molecular clouds (from low to high densities) are indicated with
black patches for reference. An alternative set of parameters is shown in the
bottom panel: for a mass resolution of order $10^{5}$ to
$10^{6}\,\mathrm{M}_{\odot}$, a softening scale of 200 pc and a minimum
smoothing length of 2 pc ($h_{\mathrm{min,ratio}}=0.006$) push the runaway
collapse zone to densities $\gtrsim 10^{8}\,\mathrm{cm}^{-3}$). All phases in
the ISM can therefore be modelled without triggering a numerical runaway
collapse.

The runaway collapse zone (equations 9 and 10) is indicated as a red shaded
area in each panel. Because the ISM is described by a single phase medium that
follows an effective temperature (thick grey line), the gas can enter the
runaway collapse zone only at very high densities
($n_{\mathrm{H}}>10^{6}\,\mathrm{cm}^{-3}$) and the eagle and apostle
simulations do not suffer from artificial runaway collapse and incorrect gas
density estimates. Without this effective temperature floor, these simulations
would however have been affected by a numerical runaway collapse at densities
$n_{\mathrm{H}}\gtrsim 100\,\mathrm{cm}^{-3}$. The adoption of a polytropic
equation of state with a polytropic index of $\gamma_{\mathrm{eff}}=4/3$ in
eagle was motivated by the numerical advantage of a constant Jeans mass for
the ISM and a constant ratio of the Jeans length to the smoothing length
(Schaye & Dalla Vecchia, 2008). Based on Bate & Burkert (1997), Schaye & Dalla
Vecchia (2008) argue that a constant Jeans mass prevents artificial
fragmentation because the (Newtonian) Jeans mass is always (marginally)
resolved by the same number of particles.

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.

Figure 5: Gravitational stability at the hydrodynamical resolution limit,
$l_{\mathrm{smooth}}$, for simulations with a constant softening length
$l_{\mathrm{soft}}$. Perturbations at length scale, $l_{\mathrm{smooth}}$, are
expected to grow in Newtonian gravity for
$\lambda_{\mathrm{J,N}}<l_{\mathrm{smooth}}$ (“gravitationally unstable in
nature”) and to decay for $\lambda_{\mathrm{J,N}}>l_{\mathrm{smooth}}$
(“gravitationally stable in nature”), separated by the thick dashed line with
a slope of 1/3 ($T\propto n_{\mathrm{H}}^{1/3}$) in each panel for
$l_{\mathrm{smooth,min}}=0$. Here, the gas density, $n_{\mathrm{H}}$, refers
to the SPH-estimated density. In simulations with softened gravity, such
perturbations are expected to grow for
$\lambda_{\mathrm{J,s}}<l_{\mathrm{smooth}}$ (“gravitationally unstable in
sim”) and to decay for $\lambda_{\mathrm{J,s}}>l_{\mathrm{smooth}}$
(“gravitationally stable in sim”), separated by the solid line which diverges
from the dashed line and follows a slope of -2/3 ($T\propto
n_{\mathrm{H}}^{-2/3}$) at high densities. The left panel shows the boundary
$\lambda_{\mathrm{J,s}}=l_{\mathrm{smooth}}$ for
$m_{\mathrm{B}}=10^{5}\,\mathrm{M}_{\odot}$ and a constant softening length
$l_{\mathrm{soft}}=100\,\mathrm{pc}$ (thick solid line). In the middle panel
the solid lines represent $\lambda_{\mathrm{J,s}}=l_{\mathrm{smooth}}$ for
different values for $l_{\mathrm{soft}}$ (see contour labels) and
$m_{\mathrm{B}}=10^{5}\,\mathrm{M}_{\odot}$. The right panel shows the
dependence of both $\lambda_{\mathrm{J,N}}=l_{\mathrm{smooth}}$ (thick green
lines) and $\lambda_{\mathrm{J,s}}=l_{\mathrm{smooth}}$ (thin black lines) on
the particle mass $m_{\mathrm{B}}$ for $l_{\mathrm{soft}}=100\,\mathrm{pc}$
and $m_{\mathrm{B}}=10^{5}/8\,,10^{5},\,\mathrm{and}\,10^{5}\times
8\,\mathrm{M}_{\odot}$ (dotted, dashed, dash-dotted lines, respectively).
Typical densities and temperatures for the WNM, CNM, and MCs are indicated
with dark patches, as in Fig. 4.

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_{\mathrm{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. Counter-intuitively, a _smaller_ softening
length would decrease the softened Jeans mass in the cold neutral medium (CNM)
and molecular clouds (MCs), which means that the softened Jeans mass would be
resolved by fewer particles in these phases. With a minimum smoothing length
of $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}}>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 multi-phase 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. If the ISM in the eagle simulations had contained the CNM
or the molecular cloud phase, i.e. if gas had not been limited by the
effective temperature floor, the choices for the gravitational softening and
the minimum smoothing length would have led to issues already for gas
densities of $n_{\mathrm{H}}>100\,\mathrm{cm}^{-3}$ (see top and middle panels
of Fig. 4), even for the high-resolution simulation. For this case, the
following options, or combinations thereof, are suggested to ensure numerical
stability by keeping gas particles out of the runaway collapse zone: (i)
_increasing_ the gravitational softening length, (ii) decreasing the values of
the minimum smoothing length; moderate changes already have a big impact (see
bottom panel in Fig. 3) or (iii) choosing a star formation condition that
converts gas particles to star particles on a very short timescale when they
enter the runaway collapse zone.

###  3.2 Gravitational instabilities at and below the resolution limits

In subsection 3.1 we defined a numerical instability that is caused by an
incorrect density estimate if the smoothing length is limited by a constant
minimum value. In this subsection, we focus on the formation of gas clumps
through gravitational fragmentation and collapse, considering both codes with
adaptive and constant constant values for the gravitational softening length.

For a code-independent discussion of the conditions under which gravitational
instabilities develop in simulations, we define the code-agnostic length
scales $l_{\mathrm{soft}}$, over which gravitational forces are softened, and
$l_{\mathrm{smooth}}$ (with a minimum of $l_{\mathrm{smooth,min}}$), over
which the hydrodynamical forces are smoothed. We explain briefly how they
connect to the smoothing lengths $h$ and the softening lengths $\epsilon$ in
some codes as examples:

We use the kernel support as measure of $l_{\mathrm{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_{\mathrm{smooth}}$ would be equal to $h_{\mathrm{Gadget4}}$.

The softening length is typically given as the Plummer-equivalent softening
length, $\epsilon$. Depending on the softening potential, the gravitational
force is exactly Newtonian for particle separations of $\geq 2.8\epsilon$
(cubic spline kernel as in Gadget4, Springel et al., 2021) or $\geq 3\epsilon$
(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 $\lesssim 1-1.5\epsilon$ (see e.g. figure 1 in Springel et
al., 2021). We therefore use $l_{\mathrm{soft}}=1.5\epsilon$ as an estimate
for the softening length scale.

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

####  3.2.1 Gravitational instabilities at the hydro resolution limit
$l_{\mathrm{smooth}}$

Gravitational instabilities cannot be modelled accurately in the simulation if
the Newtonian Jeans length, $\lambda_{\mathrm{J,N}}$, is smaller than the size
of an individual smoothing kernel, $l_{\mathrm{smooth}}$. The hydrodynamical
forces are inaccurate on scales below $l_{\mathrm{smooth}}$ and fragmentation
on scales of $\lambda_{\mathrm{J,N}}<l_{\mathrm{smooth}}$ is unresolved. For
the Newtonian Jeans length, the condition
$\lambda_{\mathrm{J,N}}=l_{\mathrm{smooth}}$ depends on the gas density and
temperature as well as on the particle mass, $m_{\mathrm{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_{\mathrm{soft}}$.
The condition $\lambda_{\mathrm{J,N}}=l_{\mathrm{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 $\lambda_{\mathrm{J,N}}>l_{\mathrm{smooth}}$.

For gravitational instabilities at the hydrodynamic resolution limit,
$l_{\mathrm{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_{\mathrm{smooth}}$, for simulations with constant and adaptive softening
lengths.

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

##### Constant softening length:

In nature (or in simulations with $m_{\mathrm{B}}\rightarrow 0$,
$l_{\mathrm{soft}}\rightarrow 0$, and $l_{\mathrm{smooth,min}}\rightarrow 0$),
density perturbations with a length scale of $l_{\mathrm{smooth}}$ grow if
$\lambda_{\mathrm{J,N}}<l_{\mathrm{smooth}}$, and decay until they reach a new
equilibrium for $\lambda_{\mathrm{J,N}}>l_{\mathrm{smooth}}$. The border888The
condition $\lambda_{\mathrm{J,N}}=l_{\mathrm{smooth}}$ can also be interpreted
as the condition that the Jeans mass,
$M_{\mathrm{J,N}}=4\pi\lambda_{\mathrm{J,N}}^{3}\rho/3$, is equal to the mass
in the kernel $N_{\mathrm{neigh}}m_{\mathrm{B}}\approx
l_{\mathrm{smooth}}^{3}\rho$ (with a constant pre-factor of order unity).
However, the number of neighbours, $N_{\mathrm{neigh}}$, is not well defined
and we therefore use the better defined condition
$\lambda_{\mathrm{J,N}}=l_{\mathrm{smooth}}$ instead. (i.e. case
$\lambda_{\mathrm{J,N}}=l_{\mathrm{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
$\lambda_{\mathrm{J,N}}>l_{\mathrm{smooth}}$ (unshaded area).

For softened gravity with a constant softening length (left and right panel:
$l_{\mathrm{soft}}=100\,\mathrm{pc}$, middle panel: various values for
$l_{\mathrm{soft}}$, see the labels), the solid lines
($\lambda_{\mathrm{J,s}}=l_{\mathrm{smooth}}$) separate perturbations with
length scales $l_{\mathrm{smooth}}$ that are gravitationally (un)stable in
softened gravity, i.e. in the simulation. Densities and temperatures for which
gas is gravitationally stable in nature for perturbations at the resolution
limit $l_{\mathrm{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
($\lambda_{\mathrm{J,N}}<l_{\mathrm{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_{\mathrm{smooth}}$ because the softened Jeans
length exceeds the kernel size ($\lambda_{\mathrm{J,s}}>l_{\mathrm{smooth}}$)
and therefore perturbations on scales of $l_{\mathrm{smooth}}$ decay. The
middle panel shows that the boundary
$\lambda_{\mathrm{J,s}}=l_{\mathrm{smooth}}$ depends on the value for the
constant softening length, $l_{\mathrm{soft}}$. The right panel shows that
increasing (dash-dotted lines) or decreasing (dotted lines) the particle mass
by a factor of 8 from the fiducial value of
$m_{\mathrm{B}}=10^{5}\,\mathrm{M}_{\odot}$ (dashed lines), affects both the
thick green lines for $\lambda_{\mathrm{J,N}}=l_{\mathrm{smooth}}$ as well as
the thin black lines for the condition
$\lambda_{\mathrm{J,s}}=l_{\mathrm{smooth}}$, because
$l_{\mathrm{smooth}}\propto m_{\mathrm{B}}^{1/3}$ (right panel of Fig. 5 for
$l_{\mathrm{soft}}=100\,\mathrm{pc}$).

##### Adaptive softening length:

For adaptive softening, typically $l_{\mathrm{soft}}=l_{\mathrm{smooth}}$ and
therefore $\lambda_{\mathrm{J,s}}\approx\lambda_{\mathrm{J,N}}$. The
instability criteria for perturbations at the resolution limit, i.e. at the
kernel size $l_{\mathrm{smooth}}$, therefore follow the instability criteria
for Newtonian gravity. The left panel of Fig. 6 (linestyles as in Fig. 5)
illustrates this for parameters representative for the FIREbox (Feldmann et
al., 2023) simulation project with $m_{\mathrm{B}}=6\times
10^{4}\,\mathrm{M}_{\odot}$. Their adaptive softening length is equal to the
average gas particle separation
($l_{\mathrm{soft}}=1.5\epsilon(m_{\mathrm{B}}/\rho)^{1/3}$, for gas density
$\rho$), 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_{\mathrm{smooth}}$ therefore follow the Newtonian Jeans criteria. The
small offset between the dashed and solid lines is from the order of unity
prefactors related to the kernel shapes and exact definitions of
$l_{\mathrm{smooth}}$ and $l_{\mathrm{soft}}$. For simplicity we use the same
kernel shapes and pre-factors as in Fig. 5.

The slope change of the $\lambda_{\mathrm{J,N}}=l_{\mathrm{smooth}}$ (dashed)
and $\lambda_{\mathrm{J,s}}=l_{\mathrm{smooth}}$ (solid) lines is caused by
the minimum smoothing length $l_{\mathrm{smooth,min}}$. The critical density
$n_{\mathrm{H,hmin}}$ above which densities are over-smoothed, and hence
$n_{\mathrm{H,SPH}}<n_{\mathrm{H}}$, is $\approx 2\times
10^{6}\,\mathrm{cm}^{-3}$ (equation 3.1 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_{\mathrm{smooth,min}}$ and $l_{\mathrm{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 (Marinacci et al., 2018; Naiman et al., 2018;
Nelson et al., 2018; Pillepich et al., 2018b; Springel et al., 2018; Pillepich
et al., 2019; Nelson et al., 2019) 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}$,
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 ($\epsilon=\epsilon_{\mathrm{gas,min}}$) follow
the softened Jeans criteria outlined in section 2. The right panel of Fig. 6
shows the $\lambda_{\mathrm{J,N}}=l_{\mathrm{smooth}}$ (dashed) and
$\lambda_{\mathrm{J,s}}=l_{\mathrm{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_{\mathrm{smooth,min}}$.

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

####  3.2.2 Gravitational instabilities below the hydro resolution limit
$l_{\mathrm{smooth}}$

Density perturbations grow in nature on length-scales smaller than
$l_{\mathrm{smooth}}$ in gas with temperatures and densities for which
$\lambda_{\mathrm{J,N}}<l_{\mathrm{smooth}}$ (shaded regions in Figs. 5 and
6). The gravitational collapse of these perturbations within an individual
smoothing kernel is not resolved and therefore not expected to be modelled
accurately with any simulation method.

We can see from Figs. 5 and 6 that a large fraction of the neutral ISM in
simulations might form sub-kernel instabilities
($\lambda_{\mathrm{J,s}}<l_{\mathrm{smooth}}$, dark shaded regions). While the
resulting clumpy particle configurations within a kernel might have a limited
effect on the simulation because the cooling rates, star formation rates and
other density- or pressure-dependent sub-grid models use the smoother SPH
estimates, we briefly discuss the differences in the treatment of sub-kernel
perturbations between codes with constant and adaptive softening lengths.

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

##### Constant softening length:

Gas with $\lambda_{\mathrm{J,s}}<l_{\mathrm{smooth}}$ is expected to be
gravitationally unstable when exposed to fluctuations on length scales between
$\lambda_{\mathrm{J,s}}$ and $l_{\mathrm{smooth}}$. In the derivation of the
softened Jeans length, $\lambda_{\mathrm{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
$\lambda_{\mathrm{J,s}}$. Therefore, the softened Jeans criteria serve as an
upper limit to the expected instabilities in the simulation because the
hydrodynamic forces might be underestimated due to the smoothing over the full
kernel (an extreme case is discussed in subsection 3.1). If
$l_{\mathrm{smooth}}<l_{\mathrm{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_{\mathrm{soft}}$ can be (much) smaller than the size of a smoothing kernel,
for a constant value of $l_{\mathrm{soft}}$. In this case, dense particle
configurations with sizes possibly even smaller than $\lambda_{\mathrm{J,N}}$
can form because the density and pressure estimates of sub-kernel clumps can
be inaccurate.

Fig. 7 shows an overview of the behaviour of dense particle configurations
within a smoothing kernel. Gas that is gravitationally unstable at the scale
of an individual smoothing kernel
($\lambda_{\mathrm{J,s}}<l_{\mathrm{smooth}}$, as in Fig. 5) is further split
into densities for which further fragmentation to even smaller scales is
suppressed ($l_{\mathrm{soft}}>l_{\mathrm{smooth}}$) or potentially induced
($l_{\mathrm{soft}}<l_{\mathrm{smooth}}$), with the boundary,
$l_{\mathrm{soft}}=l_{\mathrm{smooth}}$, indicated as vertical, dotted lines
(left panel: $l_{\mathrm{soft}}=100\,\mathrm{pc}$, right panel: various values
for $l_{\mathrm{soft}}$, see contour labels).

We conclude that for a constant softening length, gas in simulations with
$\lambda_{\mathrm{J,s}}<l_{\mathrm{smooth}}$ and
$l_{\mathrm{soft}}<l_{\mathrm{smooth}}$ (purple areas in Fig. 7) might be
affected by sub-kernel clumping. For example: within a kernel of 100
particles, a subset of 10 particles has particle separations that are smaller
than average for this kernel and hence smaller than the smoothing length,
$l_{\mathrm{smooth}}$. In contrast to simulations with adaptive softening, the
constant value for $l_{\mathrm{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_{\mathrm{smooth}}=l_{\mathrm{soft}}$ (vertical lines in Fig. 7). For higher
gas densities (i.e. $l_{\mathrm{smooth}}<l_{\mathrm{soft}}$), the further
collapse is suppressed. The detailed impact of sub-kernel clumping on galaxy
properties in large-scale simulations is beyond the scope of this work.

If sub-kernel clumping (“subk") is undesired in simulations with constant
softening, the conditions $\lambda_{\mathrm{J,s}}<l_{\mathrm{smooth}}$ and
$l_{\mathrm{soft}}<l_{\mathrm{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_{\mathrm{soft}}<l_{\mathrm{smooth}}$ to
densities below $n_{\mathrm{max,subk}}$, this condition translates into a
minimum constant softening length

$\displaystyle l_{\mathrm{min,subk}}=1.5\epsilon_{\mathrm{min,subk}}$
$\displaystyle>1.5\frac{\gamma_{\mathrm{k}}\eta_{\mathrm{res}}}{1.5}\left(\frac{X_{\mathrm{H}}m_{\mathrm{B}}}{m_{\mathrm{H}}n_{\mathrm{max,subk}}}\right)^{1/3}\;\mathrm{or}$
$\displaystyle l_{\mathrm{min,subk}}=1.5\epsilon_{\mathrm{min,subk}}$
$\displaystyle>160\,\mathrm{pc}\left(\frac{n_{\mathrm{max,subk}}}{10\,\mathrm{cm}^{-3}}\right)^{-1/3}\left(\frac{m_{\mathrm{B}}}{10^{5}\,\mathrm{M}_{\odot}}\right)^{1/3}$
(13)

for the smoothing length and softening length definition as above.

Unlike in simulations with adaptive softening, unresolved gravitational
instabilities are not suppressed by design in simulations with constant
softening lengths. Fulfilling equation (13) avoids potentially undesired, sub-
kernel gravitational instabilities in large parts of the cold ISM, but at the
expense of suppressing physical instabilities at the kernel scale,
$l_{\mathrm{smooth}}$ (see middle panel of Fig. 5).

##### Adaptive softening length:

In simulations with an adaptive softening length and $l_{\mathrm{soft}}\geq
l_{\mathrm{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 $\lambda_{\mathrm{J,N}}<l_{\mathrm{smooth}}$ are
therefore artificially suppressed.

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

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

##  4 Galaxy simulations

In section 3.2, we used the smoothing and softening lengths
$l_{\mathrm{smooth}}$ and $l_{\mathrm{soft}}$ for a code independent
discussion on the individual zones. In this section we use simulations of
isolated galaxies with the public Swift code999The simulations use version
0.9.0 and in particular revision v0.9.0-1182-g423e9dd8. (Schaller et al.,
2016, 2018; Schaller et al., 2023, www.swiftsim.com). The resolution
parameters set by the user are $\epsilon$ and $h_{\mathrm{min}}$ which relate
to $l_{\mathrm{soft}}=1.5\epsilon$ and
$l_{\mathrm{smooth,min}}=1.94h_{\mathrm{min}}$, respectively (see text in
section 3.2 for details). We will demonstrate the numerical issues that can
occur for different choices of $\epsilon$ and $h_{\mathrm{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, $\epsilon$.

The initial conditions were created with MakeNewDisk which is based on the
code used in Springel et al. (2005a) but modified to use the exact definition
of the analytic dark matter halo mass $M_{\mathrm{200}}$ (i.e. mass within
$R_{\mathrm{200}}$, the radius within which the average density is 200 times
the critical density of the Universe) instead of the halo mass integrated to
infinity (see Nobels et al., 2023 for details). The isolated galaxy initially
has a gas mass of $1.644\times 10^{10}\,\mathrm{M}_{\odot}$ with solar
metallicity ($Z_{\odot}=0.0134$, Asplund et al., 2009) and an exponential
stellar disk with a radial scale length of $4.3\,\mathrm{kpc}$ and a mass of
$3.836\times 10^{10}\,\mathrm{M}_{\odot}$ The initial disk gas fraction is 30
per cent and the gas is initialized with a temperature of
$10^{4}\,\mathrm{K}$. The baryonic disk 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 et al. (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, an interstellar
radiation field and a UV background. 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,
using a heating temperature of $10^{7.5}\,\mathrm{K}$) into the gas particle
closest to the star (“Min distance” model in Chaikin et al., 2022). Star
formation is limited to densities of $n_{\mathrm{H}}>0.1\,\mathrm{cm}^{-3}$
and cold gas (temperatures of $T<1000\,\mathrm{K}$) and the star formation
rate for each gas particle with mass $m_{\mathrm{gas}}$ is given by the
Schmidt (1959) law

$\displaystyle\dot{m}_{\star}$
$\displaystyle=\frac{e_{\mathrm{sf}}}{t_{\mathrm{ff,N}}}m_{\mathrm{gas}}$ (14)
$\displaystyle=7.2\times
10^{-5}\,\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}\left(\frac{e_{\mathrm{sf}}}{0.32\%}\right)\left(\frac{n_{\mathrm{H}}}{100\,\mathrm{cm}^{-3}}\right)^{1/2}\left(\frac{m_{\mathrm{gas}}}{10^{5}\,\mathrm{M}_{\odot}}\right)$

with the star formation efficiency $e_{\mathrm{sf}}$ and the Newtonian free-
fall time $t_{\mathrm{ff,N}}$ (equation 25). 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_{\mathrm{min,ratio}}=[0.02,0.063,0.2]$ for $\epsilon=250\,\mathrm{pc}$),
and baryon particle mass

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

$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 re-calculate the gas densities based on the gas particle positions for
$h_{\mathrm{min}}\rightarrow 0$ by restarting the original simulations at
$t=1\,\mathrm{Gyr}$ for one very small timestep ($\Delta
t_{\mathrm{max}}=10\,\mathrm{yr}$) with a very small minimum smoothing
length101010We use a non-zero value for $h_{\mathrm{min,ratio}}$ but the
selected value is small enough that the re-calculated smoothing lengths are
larger than $h_{\mathrm{min}}$ for all gas particles.
($h_{\mathrm{min,ratio}}=10^{-4}$). The snapshot file that is produced after
this timestep contains the correct (i.e. not limited by a minimum smoothing
length) gas densities based on the gas particle positions.

Fig. 8 shows the difference between the original SPH densities (top row) and
the recalculated densities for $h_{\mathrm{min}}\rightarrow 0$ (bottom row)
for 3 simulations (first three columns; all with
$m_{\mathrm{B}}=10^{5}\,\mathrm{M}_{\odot}$, $\epsilon=250\,\mathrm{pc}$,
$e_{\mathrm{sf}}=0.32\,\mathrm{per\,cent}$) with decreasing values for
$h_{\mathrm{min}}$ (from left to right). The runaway collapse described in
section 3.1 is obvious in the simulation with the largest value for
$h_{\mathrm{min}}$ (1st column). The distribution of recalculated, “real”
densities (bottom row) extends to much higher values than the densities used
during the simulation (top row) for many particles as they approach the
runaway collapse zone (red shaded area). For the simulation with the smallest
value of $h_{\mathrm{min}}$ (3rd column), the two density-temperature diagrams
look identical.

The panels in the rightmost column show the cumulative mass fraction of
particles above density $n_{\mathrm{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_{\mathrm{H,hmin}}$ (equation 3.1; 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}}>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 disk, they would persist in the simulation with better resolved
hydrodynamic forces. We therefore conclude that the very high densities seen
in the particle distributions are indeed artificial as outlined in section
3.1.

While the dense gas clumps form, they may become eligible to star formation
and if their densities are underestimated, their star formation rates are
underestimated as well ($\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_{\mathrm{min}}$ which is fulfilled at densities of
$\approx 10\,\mathrm{cm}^{-3}$ in all simulations in this work. Comparing the
values of $n_{\mathrm{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 stand-alone routine on the star particles within a
snapshot output. We use a small linking length of 25 pc and a minimum number
of particles of 50 for simulations with
$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 $\approx 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}}>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\mathrm{(zoneC)}[\mathrm{M}_{\odot}]=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.

Figure 10: Stellar mass surface density (top row) and gas phase-space
distribution (bottom row) as in Fig. 9 for the fiducial simulation parameters
(first and third columns) and a set of alternative simulation parameters (as
in Table 2 and Fig. 4; second and fourth columns) for a particle mass of
$m_{\mathrm{B}}=10^{5}\,\mathrm{M}_{\odot}$ (left two columns) and
$m_{\mathrm{B}}=8\times 10^{5}\,\mathrm{M}_{\odot}$ (right two columns), and a
star formation efficiency of $e_{\mathrm{sf}}=0.32\,\mathrm{per\,cent}$. The
combination of a larger gravitational softening scale
($\epsilon=500\,\mathrm{pc}$) and the smaller minimum smoothing length
($h_{\mathrm{min}}=7.75\,\mathrm{pc}$) result in a galaxy without artificial
stellar clumps and without gas in the runaway collapse zone. The low star
formation efficiency still results in a large amount of cold gas with
$T<100\,\mathrm{K}$ but at better resolved gas densities.

The number of clumps as well as their masses are not very sensitive to the
mass resolution (increase $m_{\mathrm{B}}$; lower right). On the other hand,
increasing the gravitational force softening length (increase $\epsilon$;
upper right) or increasing the star formation efficiency (increase
$e_{\mathrm{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 different111111The colourscales are the same across
all phase-space plots.: while the simulation with the higher star formation
efficiency ($e_{\mathrm{sf}}=1\,\mathrm{per\,cent}$) has very little cold
($<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 disk.

In contrast to the comparisons shown in Figs. 8 and 9 for which only one
parameter is varied at a time, we show in Fig. 10 an alternative combination
of values for $\epsilon$ and $h_{\mathrm{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_{\mathrm{min}}$ and equation (13) for $\epsilon$. Note that these
parameters are not necessarily the best choice for every application but
illustrate the impact of choosing values for $\epsilon$ and $h_{\mathrm{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 ($<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 star formation rates averaged over 10 Myr (large panel)
and the relative computing time (small panel on the right) for simulations
with a particle mass of $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_{\mathrm{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 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 $\approx 2$ when reducing $h_{\mathrm{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_{\mathrm{min}}$ and therefore more clumping from runaway collapse is caused
by substantial amounts of shock-heated gas at high densities $\log
n_{\mathrm{H}}[\mathrm{cm}^{-3}]\gtrsim 2$) and temperatures of a few hundred
to a few thousand K (see first and third column of Fig. 10). The temperature
increase of a factor of 100 compared to simulations with smaller values of
$h_{\mathrm{min}}$ (second and fourth column of Fig. 10) reduces the timestep
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
timestep size because $\Delta t$ is calculated from the (underestimated) SPH
densities.

Summarizing, a too large value for $h_{\mathrm{min}}$ does not only produce
artificially dense clumps of gas and star particles, but also increases the
computing time by a factor of $\approx 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_{\mathrm{sf}}$ mean that gas particles are converted into
star particles on shorter timescales for a given density. The isolated galaxy
simulations from this work do not produce artificial clumps for
$e_{\mathrm{sf}}>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}}\geq 1\,\mathrm{per\,cent}$.

One could expect that smaller values for the gravitational softening scale
$\epsilon$ 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
($\lambda_{\mathrm{J,N}}>l_{\mathrm{smooth}}$) only depends on the particle
mass (i.e. the mass resolution). Gas with densities and temperatures such that
$\lambda_{\mathrm{J,N}}<l_{\mathrm{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 $\epsilon$ and $\epsilon_{\mathrm{min}}$, neither simulation models
gravitational instabilities correctly in the CNM nor in the MCs but both need
to choose values for $\epsilon$ or $\epsilon_{\mathrm{min}}$ that avoid
numerical issues. As described in section 3.1, this means a very small value
for $\epsilon_{\mathrm{min}}$ to avoid numerical runaway collapse for adaptive
softening (if $\epsilon_{\mathrm{min}}=h_{\mathrm{min}}$) and potentially a
larger value for $\epsilon$ for constant softening if one wishes to suppress
artificial sub-kernel gravitational instabilities (see section 3.2).

In very high (mass and force) resolution simulations, clusters of stars can
form for physical reasons, but for the mass resolution in the isolated
galaxies examples ($m_{\mathrm{B}}=10^{5}\,\mathrm{M}_{\odot}$,
$m_{\mathrm{B}}=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, hydro) also need to take into account the effect of artificial
heating of the stellar disk by dark matter and stellar halo particles (see
e.g. Ludlow et al., 2019a, 2021, 2023; Wilkinson et al., 2022), which is
reduced for larger values of $\epsilon$ (appendix D in Ludlow et al., 2021).
This effect is not included here because we do not use a live dark matter halo
for these idealised simulations.

The star formation rate for a gas particle in our simulation depends on the
Newtonian free-fall time (see equation 14) but the collapse of gas clumps in
the simulation follows the longer softened free-fall time (equation 49). While
for numerical stability considerations the softened properties (free-fall
time, Jeans mass and length) are more relevant, we use the Newtonian free-fall
time for the star formation rate because we want to model this subgrid process
physically rather than numerically (see also the discussion in Nobels et al.,
2023).

Which resolution and star formation criteria to choose depends on the
application of the simulation, but the softened Jeans criteria and the
temperature - density zones defined in section 3, in particular the runaway
collapse zone (section 3.1) provide a useful guideline to avoid undesired
numerical behaviour that may introduce artefacts and can furthermore slow down
the simulation.

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

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
$\rho$ through the effective polytropic index $\gamma_{\mathrm{eff}}$. It has
been argued that for $\gamma_{\mathrm{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)
$\gamma_{\mathrm{eff}}=4/3$ for $n_{\mathrm{H}}>0.1\,\mathrm{cm}^{-3}$ and in
the widely used Springel & Hernquist (2003) model, the polytropic index varies
with density between $\gamma_{\mathrm{eff}}\approx 2.5$ at
$0.1\,\mathrm{cm}^{-3}$ to $\gamma_{\mathrm{eff}}\approx 4/3$ at
$100\,\mathrm{cm}^{-3}$ (for $z=0$; see figure 1 of Springel & Hernquist, 2003
for details). We argue here that such pressure (or entropy) floors are
unnecessary in modern Lagrangian simulations, even for relatively low
resolutions, provided that appropriate values for the softening and minimum
smoothing lengths are used to avoid artificially inducing gravitational
instabilities and a numerical runaway collapse.

The free-fall time in simulations with a constant softening length
($t_{\mathrm{ff,s}}$, equation 3) is longer than the Newtonian free-fall time
($t_{\mathrm{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, $\lambda_{\mathrm{J,N}}$, the so called Jeans length, by
comparing the free-fall timescale to the sound crossing timescale (see the
derivation in appendix A). In a homogeneous medium with gas density
$n_{\mathrm{H}}$ and temperature $T$, density perturbations on scales
$>\lambda_{\mathrm{J,N}}$ are gravitationally unstable (i.e. the density
perturbations grow) while density perturbations on scales
$<\lambda_{\mathrm{J,N}}$ are gravitationally stable (i.e. the density
perturbations decay with time until a new equilibrium is reached).

In this work we introduce “softened Jeans criteria” (section 2) for which we
re-derive the Jeans length in softened gravity, as used in Lagrangian
simulations to suppress 2-body scattering. In parallel to the Newtonian Jeans
criteria, the softened Jeans length (mass), $\lambda_{\mathrm{J,s}}$
($M_{\mathrm{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, $\lambda_{\mathrm{J,N}}=\epsilon$, 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, $\epsilon$, the softened
Jeans mass can exceed the Newtonian Jeans mass $M_{\mathrm{J,N}}$ by several
orders of magnitude for gas densities and temperatures typical of the cold
interstellar medium. For example, at a gas temperature of a few tens of K and
a gas density of $\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 equation of state, cold neutral and molecular ISM gas is formally
unresolved because gravitational instabilities should form within an
individual smoothing kernel ($\lambda_{\mathrm{J,N}}<l_{\mathrm{smooth}}$).
Modelling a multi-phase interstellar medium at these mass resolutions
therefore requires an understanding of how instabilities behave at and below
the resolution limit (see section 3).

If the hydrodynamic forces are resolved by at least one smoothing kernel,
gravitational instabilities would be modelled correctly (if
$\lambda_{\mathrm{J,s}}=\lambda_{\mathrm{J,N}}$), or be suppressed
($\lambda_{\mathrm{J,s}}>\lambda_{\mathrm{J,N}}$) for all gas densities and
temperatures in simulations, because softened gravity never exceeds Newtonian
gravity and therefore $\lambda_{\mathrm{J,s}}\geq\lambda_{\mathrm{J,N}}$. Yet,
we show in section 3 that perturbations can grow under specific conditions
within a hydrodynamic smoothing kernel. Two distinct pathways are identified
for numerically induced instabilities, both related to an inaccurate
calculation of the hydrodynamic properties below the resolution limit:
instabilities caused by (i) an underestimated gas density (subsection 3.1) and
(ii) pressure that is smoothed on length scales larger than those on which
gravity is softened (subsection 3.2). The former instability is relevant for
simulations with both adaptive and constant softening lengths, while the
latter only applies to simulations with a constant softening length (see
discussion in section 3).

The effects of sub-kernel instabilities is demonstrated in section 4 in
simulations of isolated galaxies using the Swift code with a constant
softening length, $\epsilon$. As outlined in subsection 3.1, the density of
gas clumps that are smaller than the minimum allowed smoothing length,
$h_{\mathrm{min}}$, are under-estimated. Re-calculating the densities of all
gas particles with a vanishing value for $h_{\mathrm{min}}$ reveals gas clumps
with several orders of magnitude higher gas densities (compare top and bottom
row of Fig. 8). These gas clumps result in dense star clusters that largely
disappear for smaller minimum smoothing length values and are therefore
artificial (Fig. 9).

Based on the analysis in section 3.1, we recommend for both simulations with
constant and adaptive softening lengths to set the minimum smoothing lengths
to a value small enough value that the smoothing lengths are not limited to a
constant value in a simulation with particle mass, $m_{\mathrm{B}}$ and a
maximum expected density, $n_{\mathrm{sim,max}}$:

$h_{\mathrm{min}}<8.3\,\mathrm{pc}\left(\frac{m_{\mathrm{B}}}{10^{5}\,\mathrm{M}_{\odot}}\right)^{1/3}\left(\frac{n_{\mathrm{sim,max}}}{10^{4}\,\mathrm{cm}^{-3}}\right)^{-1/3}\,.$
(15)

We recommend to generously fulfill this condition and could not identify a
benefit of a larger value for $h_{\mathrm{min}}$. If the minimum values for
softening and smoothing lengths are identical
($\epsilon_{\mathrm{min}}=h_{\mathrm{min}}$) in simulations with adaptive
softening, a small value for $\epsilon_{\mathrm{min}}$ is also required.

Gravitational instabilities that would form within a smoothing kernel are
unresolved by definition, even if both the minimum smoothing length and the
softening length are approaching zero121212A simulation with a given particle
mass $m_{\mathrm{B}}$ does not have a higher resolution just because it uses a
smaller value for $\epsilon$, see also discussion in section 5.. Assuming
$h_{\mathrm{min}}\rightarrow 0$, the fundamental differences that remain
between simulations with adaptive and constant softening lengths for
instabilities within a kernel are described in section 3.2. Neither method
gives the right solution for sub-kernel gravitational instabilities: an
adaptive softening length that follows the smoothing length of the kernel
leads to artificial suppression of sub-kernel instabilities while a constant
softening length can result in artificially inducing sub-kernel instabilities.
Because these length scales are by definition below the resolution limit of
the simulation, their detailed impact on galaxy properties needs further
investigation and is beyond the scope of this work.

Suppressing sub-kernel instabilities in simulations with constant softening
lengths requires a softening length that exceeds the value given by equation
13. However, increasing the softening length artificially stabilizes
gravitational instabilities at the hydrodynamical spatial resolution limit:
the size of a smoothing kernel, $l_{\mathrm{smooth}}$, for larger regions in
density-temperature space (see the middle panel in Fig. 5). The impact of sub-
kernel instabilities on galaxy properties remains to be tested, because
cooling rates, star formation rates and other density- or pressure-dependent
sub-grid models use the smoother SPH estimates. The optimal value of a
constant softening length, $\epsilon$, 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_{\mathrm{B}}$ is weak) do
not depend on an effective pressure floor for numerical
stability131313Simulations such as eagle (Schaye et al., 2015) needed the
implemented pressure floor to suppress the runaway collapse that could
otherwise occur due to their large adopted value for $h_{\mathrm{min}}$
($h_{\mathrm{min}}=108.5\,\mathrm{pc}$ for their $(100\,\mathrm{Mpc})^{3}$
simulation).. 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_{\mathrm{min}}$ satisfies equation (15). While the collapse
of individual gas clouds is suppressed or slowed down when gravity is
softened, this can be taken into account by subgrid prescriptions, e.g. for
star formation, and is generally not a bottleneck if the aim of a simulation
is to reproduce general galaxy properties.

The ideal simulation has both the mass and force resolution to accurately
model the gravitational fragmentation and collapse of individual molecular
clouds but in this work we showed a numerically stable alternative, both for
adaptive and constant softening, for projects for which this is
computationally too expensive.

## Acknowledgements

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 April 2022. This research was funded in part
by the Austrian Science Fund (FWF) grant number V 982-N. This paper made use
of the following python packages: astropy (Astropy Collaboration et al.,
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.

## References

  * Asplund et al. (2009) Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481
  * Astropy Collaboration et al. (2022) Astropy Collaboration et al., 2022, ApJ, 935, 167
  * Athanassoula et al. (2000) Athanassoula E., Fady E., Lambert J. C., Bosma A., 2000, MNRAS, 314, 475
  * Bate & Burkert (1997) Bate M. R., Burkert A., 1997, MNRAS, 288, 1060
  * Bate et al. (2003) Bate M. R., Bonnell I. A., Bromm V., 2003, MNRAS, 339, 577
  * Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. Princeton Univ. Press, Princeton, NJ
  * Bird et al. (2022) Bird S., Ni Y., Di Matteo T., Croft R., Feng Y., Chen N., 2022, MNRAS, 512, 3703
  * Borrow & Borrisov (2020) Borrow J., Borrisov A., 2020, Journal of Open Source Software, 5, 2430
  * Borrow & Kelly (2021) Borrow J., Kelly A. J., 2021, arXiv e-prints, p. arXiv:2106.05281
  * Borrow et al. (2021) Borrow J., Schaller M., Bower R. G., 2021, MNRAS, 505, 2316
  * Borrow et al. (2022) Borrow J., Schaller M., Bower R. G., Schaye J., 2022, MNRAS, 511, 2367
  * Chaikin et al. (2022) Chaikin E., Schaye J., Schaller M., Bahé Y. M., Nobels F. S. J., Ploeckinger S., 2022, MNRAS, 514, 249
  * Crain & van de Voort (2023) Crain R. A., van de Voort F., 2023, ARA&A, 61, 473
  * Dalla Vecchia & Schaye (2012) Dalla Vecchia C., Schaye J., 2012, MNRAS, 426, 140
  * Davé et al. (2019) Davé R., Anglés-Alcázar D., Narayanan D., Li Q., Rafieferantsoa M. H., Appleby S., 2019, MNRAS, 486, 2827
  * Dehnen (2001) Dehnen W., 2001, MNRAS, 324, 273
  * Dehnen & Aly (2012) Dehnen W., Aly H., 2012, MNRAS, 425, 1068
  * Dubois et al. (2014) Dubois Y., et al., 2014, MNRAS, 444, 1453
  * Fattahi et al. (2016) Fattahi A., et al., 2016, MNRAS, 457, 844
  * Feldmann et al. (2023) Feldmann R., et al., 2023, MNRAS, 522, 3831
  * Feng et al. (2018) Feng Y., Bird S., Anderson L., Font-Ribera A., Pedersen C., 2018, Mp-Gadget/Mp-Gadget: A Tag For Getting A Doi, Zenodo, doi:10.5281/zenodo.1451799
  * Gingold & Monaghan (1977) Gingold R. A., Monaghan J. J., 1977, MNRAS, 181, 375
  * Goldbaum et al. (2018) Goldbaum N. J., ZuHone J. A., Turk M. J., Kowalik K., Rosen A. L., 2018, Journal of Open Source Software, 3, 809
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Hernquist (1990) Hernquist L., 1990, ApJ, 356, 359
  * Hopkins (2015) Hopkins P. F., 2015, MNRAS, 450, 53
  * Hubber et al. (2006) Hubber D. A., Goodwin S. P., Whitworth A. P., 2006, A&A, 450, 881
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Jeans (1902) Jeans J. H., 1902, Philosophical Transactions of the Royal Society of London Series A, 199, 1
  * Lucy (1977) Lucy L. B., 1977, AJ, 82, 1013
  * Ludlow et al. (2019a) Ludlow A. D., Schaye J., Schaller M., Richings J., 2019a, MNRAS, 488, L123
  * Ludlow et al. (2019b) Ludlow A. D., Schaye J., Bower R., 2019b, MNRAS, 488, 3663
  * Ludlow et al. (2020) Ludlow A. D., Schaye J., Schaller M., Bower R., 2020, MNRAS, 493, 2926
  * Ludlow et al. (2021) Ludlow A. D., Fall S. M., Schaye J., Obreschkow D., 2021, MNRAS, 508, 5114
  * Ludlow et al. (2023) Ludlow A. D., Fall S. M., Wilkinson M. J., Schaye J., Obreschkow D., 2023, MNRAS, 525, 5614
  * Marinacci et al. (2018) Marinacci F., et al., 2018, MNRAS, 480, 5113
  * Merritt (1996) Merritt D., 1996, AJ, 111, 2462
  * Naiman et al. (2018) Naiman J. P., et al., 2018, MNRAS, 477, 1206
  * Navarro et al. (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
  * Nelson et al. (2018) Nelson D., et al., 2018, MNRAS, 475, 624
  * Nelson et al. (2019) Nelson D., et al., 2019, MNRAS, 490, 3234
  * Nobels et al. (2023) Nobels F. S. J., Schaye J., Schaller M., Ploeckinger S., Chaikin E., Richings A. J., 2023, arXiv e-prints, p. arXiv:2309.13750
  * Pillepich et al. (2018a) Pillepich A., et al., 2018a, MNRAS, 473, 4077
  * Pillepich et al. (2018b) Pillepich A., et al., 2018b, MNRAS, 475, 648
  * Pillepich et al. (2019) Pillepich A., et al., 2019, MNRAS, 490, 3196
  * Ploeckinger & Schaye (2020) Ploeckinger S., Schaye J., 2020, MNRAS, 497, 4857
  * Plummer (1911) Plummer H. C., 1911, MNRAS, 71, 460
  * Power et al. (2003) Power C., Navarro J. F., Jenkins A., Frenk C. S., White S. D. M., Springel V., Stadel J., Quinn T., 2003, MNRAS, 338, 14
  * Price & Monaghan (2007) Price D. J., Monaghan J. J., 2007, MNRAS, 374, 1347
  * Price et al. (2018) Price D. J., et al., 2018, Publ. Astron. Soc. Australia, 35, e031
  * Richings & Schaye (2016) Richings A. J., Schaye J., 2016, MNRAS, 460, 2297
  * Rodionov & Sotnikova (2005) Rodionov S. A., Sotnikova N. Y., 2005, Astronomy Reports, 49, 470
  * Romeo (1998) Romeo A. B., 1998, A&A, 335, 922
  * Sawala et al. (2016) Sawala T., et al., 2016, MNRAS, 457, 1931
  * Schaller et al. (2016) Schaller M., Gonnet P., Chalk A. B. G., Draper P. W., 2016, in Proceedings of the Platform for Advanced Scientific Computing Conference. p. 2 (arXiv:1606.02738), doi:10.1145/2929908.2929916
  * Schaller et al. (2018) Schaller M., Gonnet P., Draper P. W., Chalk A. B. G., Bower R. G., Willis J., Hausammann L., 2018, SWIFT: SPH With Inter-dependent Fine-grained Tasking, Astrophysics Source Code Library, record ascl:1805.020 (ascl:1805.020)
  * Schaller et al. (2023) Schaller M., et al., 2023, arXiv e-prints, p. arXiv:2305.13380
  * Schaye & Dalla Vecchia (2008) Schaye J., Dalla Vecchia C., 2008, MNRAS, 383, 1210
  * Schaye et al. (2015) Schaye J., et al., 2015, MNRAS, 446, 521
  * Schaye et al. (2023) Schaye J., et al., 2023, arXiv e-prints, p. arXiv:2306.04024
  * Schmidt (1959) Schmidt M., 1959, ApJ, 129, 243
  * Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  * Springel & Hernquist (2003) Springel V., Hernquist L., 2003, MNRAS, 339, 289
  * Springel et al. (2005a) Springel V., Di Matteo T., Hernquist L., 2005a, MNRAS, 361, 776
  * Springel et al. (2005b) Springel V., et al., 2005b, Nature, 435, 629
  * Springel et al. (2018) Springel V., et al., 2018, MNRAS, 475, 676
  * Springel et al. (2021) Springel V., Pakmor R., Zier O., Reinecke M., 2021, MNRAS, 506, 2871
  * Teyssier (2002) Teyssier R., 2002, A&A, 385, 337
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Vogelsberger et al. (2014) Vogelsberger M., et al., 2014, MNRAS, 444, 1518
  * Vogelsberger et al. (2020) Vogelsberger M., Marinacci F., Torrey P., Puchwein E., 2020, Nature Reviews Physics, 2, 42
  * Weinberger et al. (2017) Weinberger R., et al., 2017, MNRAS, 465, 3291
  * Weinberger et al. (2020) Weinberger R., Springel V., Pakmor R., 2020, ApJS, 248, 32
  * Wendland (1995) Wendland H., 1995, Advances in Computational Mathematics, 4, 389
  * Whitworth (1998) Whitworth A. P., 1998, MNRAS, 296, 442
  * Wilkinson et al. (2022) Wilkinson M. J., Ludlow A. D., Lagos C. d. P., Fall S. M., Schaye J., Obreschkow D., 2022, arXiv e-prints, p. arXiv:2208.07623
  * Yamamoto et al. (2021) Yamamoto Y., Okamoto T., Saitoh T. R., 2021, MNRAS, 504, 3986

##  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_{\mathrm{ff}}$, to the
sound crossing time, $t_{\mathrm{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\equiv\lambda_{\mathrm{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

$t_{\mathrm{sc}}=\frac{R}{c_{\mathrm{s}}}$ (16)

with the sound speed

$c_{\mathrm{s}}=\left(\frac{\gamma k_{\mathrm{B}}T}{\mu
m_{\mathrm{H}}}\right)^{1/2}$ (17)

and $t_{\mathrm{sc}}$ is independent of the shape of the gravitational
potential. Here, $\gamma$ is ratio of specific heats, $k_{\mathrm{B}}$ is the
Boltzmann constant, $T$ is the gas temperature, $\mu$ the mean particle mass,
and $m_{\mathrm{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:

$t_{\mathrm{ff}}=\int_{r=R}^{0}\mathrm{d}t=\int_{r=R}^{0}\frac{\mathrm{d}t}{\mathrm{d}r}\mathrm{d}r=\int_{r=R}^{0}\frac{1}{v(r)}\mathrm{d}r\,.$
(18)

Integrating the equation of motion

$\displaystyle v\frac{\mathrm{d}v}{\mathrm{d}r}$ $\displaystyle=-|a(r)|$ (19)

from $R$ to $r$ leads to

$v(r)=\left(2\int_{r^{\prime}=R}^{r}-|a(r^{\prime})|\mathrm{d}r^{\prime}\right)^{1/2}\,.$
(20)

which is used in equation 18 to calculate the free-fall time

$t_{\mathrm{ff}}(R)=\int_{r=R}^{0}v(r)^{-1}\mathrm{d}r=\int_{r=R}^{0}\left(2\int_{r^{\prime}=R}^{r}-|a(r^{\prime})|\;\mathrm{d}r^{\prime}\right)^{-1/2}\mathrm{d}r\,.$
(21)

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\pi R^{3}/3\rho$. Finally, the Jeans
length is calculated by solving the equation
$t_{\mathrm{ff}}(R)=t_{\mathrm{sc}}(R,c_{\mathrm{s}})$ for $R$.

###  A.1 Newtonian (unsoftened) gravity

The Jeans length for Newtonian gravity can be found in many textbooks but
depending on the derivation the prefactors 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 prefactors are included in this derivation141414A pre-
factor of order unity is not important due to the approximate nature of the
Jeans criterion in realistic applications but for comparisons to the
softenened Jeans criteria and for the sake of completeness, we keep all pre-
factors..

For Newtonian gravity,

$|a(r)|=|a_{\mathrm{N}}(r)|=\frac{GM}{r^{2}}$ (22)

and the velocity for the free-fall is

$v_{\mathrm{N}}(r)=\left(2GM\right)^{1/2}\left(\frac{1}{r}-\frac{1}{R}\right)^{1/2}\quad.$
(23)

The free-fall time is

$t_{\mathrm{ff,N}}=\left(\frac{\pi^{2}R^{3}}{8GM}\right)^{1/2}=\left(\frac{3\pi}{32G\rho}\right)^{1/2}$
(24)

or

$\displaystyle t_{\mathrm{ff,N}}$
$\displaystyle=4.4\,\mathrm{Myr}\left(\frac{n_{\mathrm{H}}}{100\,\mathrm{cm}^{-3}}\right)^{-1/2}$
(25)

for the gas density $\rho=n_{\mathrm{H}}m_{\mathrm{H}}/X_{\mathrm{H}}$ with
the hydrogen number density, $n_{\mathrm{H}}$, the hydrogen mass fraction
$X_{\mathrm{H}}$ (here we use $X_{\mathrm{H}}=0.74$ for solar chemical
composition, Asplund et al., 2009). The Jeans length for Newtonian gravity is

$\displaystyle\lambda_{\mathrm{J,N}}$ $\displaystyle=t_{\mathrm{ff,N}}\cdot
c_{\mathrm{s}}$ (26) $\displaystyle=\left(\frac{3\pi\gamma
X_{\mathrm{H}}k_{\mathrm{B}}T}{32G\mu
m_{\mathrm{H}}^{2}n_{\mathrm{H}}}\right)^{1/2}$ (27)
$\displaystyle=1.5\,\mathrm{pc}\,\left(\frac{T}{10\,\mathrm{K}}\right)^{1/2}\left(\frac{n_{\mathrm{H}}}{100\,\mathrm{cm}^{-3}}\right)^{-1/2}\;,$
(28)

which corresponds to a Jeans mass of

$\displaystyle M_{\mathrm{J,N}}$
$\displaystyle=\frac{4\pi\rho\lambda_{\mathrm{J,N}}^{3}}{3}$
$\displaystyle=46\,\mathrm{M}_{\odot}\,\left(\frac{T}{10\,\mathrm{K}}\right)^{3/2}\left(\frac{n_{\mathrm{H}}}{100\,\mathrm{cm}^{-3}}\right)^{-1/2}\;,$
(29)

where we use $\gamma=5/3$ and $\mu=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 $\gamma=1.4$ and $\mu=2$, but
these are comparable to other order of unity effects, such as the concrete
particle distribution.

Figure 12: Top panels: Ratios of gravitationally softened Jeans length
$\lambda_{\mathrm{J,s}}$ to Newtonian Jeans length $\lambda_{\mathrm{J,N}}$
(black lines and symbols, left y-axis) and softened Jeans mass to Newtonian
Jeans mass (red lines and symbols, right y-axis) for a Plummer potential with
a softening length $\epsilon$ (left panels) and a Wendland C2 potential with a
softening length $H=3\epsilon$ (right panels). The solid lines in the top
panels are the fits (equations 36, 37 for the Plummer potential and equations
51, 52 for the Wendland C2 potential) to the numerical root finding (points,
only selected points shown) and the bottom panels illustrate the accuracy of
the fit.

###  A.2 Plummer softening

For a Plummer (1911) softened potential with

$|a(r)|=|a_{\mathrm{P}}(r)|=\frac{GMr}{\left(r^{2}+\epsilon^{2}\right)^{3/2}}\,,$
(30)

and a softening scale $\epsilon$, the free-fall velocity, equation (20),
becomes

$v_{\mathrm{P}}(r)=\left(2GM\right)^{1/2}\left(\frac{1}{\left(r^{2}+\epsilon^{2}\right)^{1/2}}-\frac{1}{\left(R^{2}+\epsilon^{2}\right)^{1/2}}\right)^{1/2}\,,$
(31)

and the free-fall time (equation 21) is

$\displaystyle t_{\mathrm{ff,P}}=$ $\displaystyle\left(2GM\right)^{-1/2}\cdot$
$\displaystyle\int_{r=R}^{0}\left(\frac{1}{\left(r^{2}+\epsilon^{2}\right)^{1/2}}-\frac{1}{\left(R^{2}+\epsilon^{2}\right)^{1/2}}\right)^{-1/2}\mathrm{d}r\;.$
(32)

This expression can be re-written as

$t_{\mathrm{ff,P}}=t_{\mathrm{ff,N}}\int_{x=0}^{1}\frac{2}{\pi}\left(\left[x^{2}+f_{\mathrm{P}}^{2}\right]^{-1/2}-\left[1+f_{\mathrm{P}}^{2}\right]^{-1/2}\right)^{-1/2}\mathrm{d}x\;,$
(33)

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

$t_{\mathrm{ff,P,fit}}=t_{\mathrm{ff,N}}\left(1+2^{2/3}f_{\mathrm{P}}^{2}\right)^{3/4}\;.$
(34)

The free-fall in a Plummer softened potential is therefore slowed down by a
factor of $\approx 2.0$ if the softening scale equals the size of the clump
($f_{\mathrm{P}}=1$) and by a factor $\approx 45$ for a softening scale that
is 10 times the size of the clump ($f_{\mathrm{P}}=10$).

For the limiting case $f_{\mathrm{P}}\gg 1$ (i.e. softening is important) the
Plummer free-fall time can be calculated analytically (see Nobels et al.,
2023) and for $f_{\mathrm{P}}\ll 1$ the Plummer free-fall time approaches the
Newtonian free-fall time. For the full expression of the Jeans length (solving
$t_{\mathrm{ff,P,fit}}=R/c_{\mathrm{s}}$ for $R$ with
$\lambda_{\mathrm{J,N}}=t_{\mathrm{ff,N}}c_{\mathrm{s}}$, the Newtonian Jeans
length) we use a root finding algorithm for

$\left(\frac{\lambda_{\mathrm{J,P}}}{\lambda_{\mathrm{J,N}}}\right)^{10/3}-\left(\frac{\lambda_{\mathrm{J,P}}}{\lambda_{\mathrm{J,N}}}\right)^{2}-2^{2/3}\left(\frac{\epsilon}{\lambda_{\mathrm{J,N}}}\right)^{2}=0\;.$
(35)

This dimensionless equation is solved numerically for values of
$\epsilon/\lambda_{\mathrm{J,N}}$ between $0.01$ and $10^{5}$. The results are
fit by

$\lambda_{\mathrm{J,P,fit}}=\lambda_{\mathrm{J,N}}\left(1+1.42\left(\frac{\epsilon}{\lambda_{\mathrm{J,N}}}\right)^{3/2}\right)^{2/5}\;,$
(36)

and the softened Jeans for the Plummer potential is fit by

$M_{\mathrm{J,P,fit}}=M_{\mathrm{J,N}}\left(1+1.42\left(\frac{\epsilon}{\lambda_{\mathrm{J,N}}}\right)^{3/2}\right)^{6/5}\;.$
(37)

The top left panel in Fig. 12 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, $\epsilon$, and the Newtonian Jeans length
$\lambda_{\mathrm{J,N}}$. For small values of
$\epsilon/\lambda_{\mathrm{J,N}}$ (i.e. the gas cloud is much larger than the
softening scale), the softened Jeans length approaches the Newtonian Jeans
length ($\log\lambda_{\mathrm{J,P}}/\lambda_{\mathrm{J,N}}\rightarrow 0$) but
for large values of $\epsilon/\lambda_{\mathrm{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 (36) and (37) 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).

###  A.3 Wendland C2 softening

For the Wendland (1995) C2 kernel the acceleration $|a|=|a_{\mathrm{W}}|$ is
defined in two parts, $a_{\mathrm{W}}(r<H,H)$ and $a_{\mathrm{W}}(r\geq
H)=a_{\mathrm{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\epsilon$, where $\epsilon$ is the Plummer equivalent
softening length.

If the size of the cloud is smaller than $H$ ($f_{\mathrm{W}}\equiv H/R\geq
1$), the gravitational acceleration is softened throughout the free-fall and
the velocity is calculated from

$v(r)=\left[-2\int_{r^{\prime}=R}^{r}\frac{GMr^{\prime}}{H^{3}}V(u)\mathrm{d}r^{\prime}\right]^{1/2}$
(38)

with the Wendland C2 kernel

$\displaystyle W(u)$ $\displaystyle=-3u^{7}+15u^{6}-28u^{5}+21u^{4}-7u^{2}+3$
(39) $\displaystyle V(u)$
$\displaystyle=-W^{\prime}(u)/u=21u^{5}-90u^{4}+140u^{3}-84u^{2}+14$ (40)

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\leq H$ (or:
$u\equiv r/H\leq 1$). The full integral for the velocity for $r<H$ consists
therefore of two parts:

$\displaystyle v(r)$ $\displaystyle=$
$\displaystyle\sqrt{2}\left[\int_{r^{\prime}=R}^{H}-|a_{\mathrm{N}}(r)|\,\mathrm{d}r^{\prime}+\int_{r^{\prime}=H}^{0}-|a_{\mathrm{W}}(r<H,H)|\,\mathrm{d}r^{\prime}\right]^{1/2}\;,$
(41)

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

$\displaystyle v_{\mathrm{W}}(r)=$
$\displaystyle\sqrt{2GM}\begin{cases}H^{-1/2}\left[W(u)-W\left(\frac{1}{f_{\mathrm{W}}}\right)\right]^{1/2}&\text{if
$f_{\mathrm{W}}\geq 1$}\\\\[8.61108pt]
H^{-1/2}\left[W(u)-f_{\mathrm{W}}\right]^{1/2}&\text{if
$f_{\mathrm{W}}<1\wedge u\leq 1$}\\\\[8.61108pt]
\left[\frac{1}{r}-\frac{1}{R}\right]^{1/2}&\text{if $f_{\mathrm{W}}<1\wedge
u>1$}\end{cases}$ (42)

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

##### Case 1: free-fall always softened ($f_{\mathrm{W}}\geq 1$ or $R\leq H$):

The free-fall time is