Equations of motion
StarSmasher does not use the textbook SPH equations. The formulation it integrates was derived in Appendix A of Gaburov, Lombardi & Portegies Zwart (2010), MNRAS 402, 105 (arXiv:0904.0997). The difference matters for exactly the problems this code is built for: collisions and mergers, where particles of very different mass end up mixed together in the same fluid.
The standard constraint
Standard SPH ties the smoothing length to the local density,
where \(C_i\) must carry dimensions of mass for \(h_i\) to come out as a length. Written the usual way, \(h_i = \eta\,(m_i/\rho_i)^{1/3}\), that mass is the particle’s own. In effect each particle keeps a fixed mass inside its kernel. That is unobjectionable when every particle has the same mass, and it is arguably the wrong constraint as soon as they do not.
In a stellar collision, particles from two parent stars mix. If the stars were
resolved differently, or if equalmass was used so that particle masses
already span orders of magnitude within a single star, then a particle
drifting into a region where the mean particle mass differs will find itself
with far too few neighbours, or far too many. Neither is recoverable: too few
and the density estimate is noise, too many and the resolution has quietly been
thrown away.
The constraint StarSmasher uses instead
Rather than fixing the mass inside the kernel, StarSmasher fixes an estimate of the number of neighbours. Each particle carries
and \(h_i\) is solved for so that \(N_i\) equals a target value. That
target is the nnopt of sph.input.
\(G\) is best understood by imagining a cruder version of itself. Suppose it were a step function, equal to 1 for a particle inside the kernel and 0 outside. Then \(N_i\) would be nothing more than a count of particle \(i\)’s neighbours, and the constraint would read “give every particle the same number of neighbours”.
\(G\) is that count made smooth. It falls off gradually instead of switching, so a particle drifting towards the kernel edge contributes a diminishing fraction rather than vanishing the instant it crosses. That is what makes the constraint differentiable, and so usable in a variational derivation, at the cost of \(N_i\) no longer being an integer count of anything.
Which \(G\) is used is set by gflag. gflag=0 selects the function
of equation (A2) and figure A1 of the paper. gflag=1, the default, is the
same except that \(G=1\) inside the inner half of the kernel, which works
better when massive compact objects or core particles are present. The
namelist default nnopt=22+gflag follows gflag, so changing one moves
the other.
gflag=2 selects \(G=(1-q^3)^2\) with \(q=r/2h\). Unlike the other
two, this \(G\) decreases strictly all the way out to \(r=2h\), so
\(N_i\) always increases with \(h_i\) whenever particle \(i\) has a
neighbour. With gflag 0 or 1, \(G\) is flat over part of the kernel and
nnopt is an integer, so \(N_i\) can be flat in \(h_i\) right at the
target. The smoothing length can then drift without the matching correction
term, and energy is not conserved. The error is small in ordinary stars but
large for gas near a massive black hole, because there a change in \(h_i\)
changes the gravitational softening a lot.
gflag=2 needs a smaller nnopt for the same number of neighbours. For
uniform density, the number of neighbours within \(2h\) is exactly
\(3\,n_{\rm opt}\), since \(3\int_0^1 (1-q^3)^2 q^2\,dq = 1/3\). In
relaxed stars it is about \(3.2\) to \(3.4\,n_{\rm opt}\), compared
with \(1.4\) to \(1.7\,n_{\rm opt}\) for gflag=1 (see below). So
gflag=2 needs roughly half the nnopt of gflag=1, and
nnopt should be set explicitly when using gflag=2.
Smoothing lengths are updated at every timestep, so this is a root find per particle per step. Standard SPH solves for \(h_i\) too, so the expense is not particular to this formulation.
How many neighbours you actually get
nnopt is the target of the constraint, not the number of neighbours. The
number of particles actually inside the kernel is larger, by a factor of roughly
1.4 to 1.7.
The factor is largest at small nnopt and settles to roughly 1.4 to 1.45
above about nnopt=60. Measurements across a polytrope, a low-mass MESA
model and a red giant span 1.38 to 1.66, and Gaburov et al. quote nnopt=22
giving 35 to 40 neighbours.
The density estimate is unchanged, and remains the ordinary kernel sum
What this costs
Because \(h_i\) now depends on particle positions through the neighbour constraint rather than through the density, the equation of motion acquires correction terms of a different form from the familiar \(\nabla h\) terms. Alongside the usual gradient of \(W\), each pair contributes a gradient of \(G\), weighted by factors (written \(\omega\) and \(\chi\) in the paper) that measure how the density and neighbour-number sums respond to a change in smoothing length. The result is still derived from a Lagrangian, so energy and momentum conservation are retained.
The price is smaller than that description suggests. Standard SPH already sweeps each particle’s neighbours to build \(\rho_i\), and already carries a \(\nabla h\) correction of its own. This formulation adds two more running totals to that same sweep, one for \(G\) and one for its derivative with respect to \(h\), together with a second table lookup for each pair because those two are evaluated on a slightly different length than the density is. There is no additional pass over pairs and no second neighbour search.
Why it is worth it
Particles of very different mass can mix freely while each retains a sensible
number of neighbours. That is precisely what a collision or merger demands, and
it is what makes equalmass a usable knob: a model can concentrate particles
where the mass is, spanning a wide range of particle masses, without any region
becoming under- or over-resolved as the star is disrupted.
Note
This is also why nnopt behaves the way Relaxing a star
describes. It is not a loose suggestion about neighbour counts but the
target of an equation solved for every particle at every step, which is why
leaving it unchanged while increasing N degrades the model.
Dynamic softening for compact objects
A compact object or core particle (a point particle, with u=0) feels and
exerts gravity only, softened over the length hco. By default that length is
fixed. With dynhco=1 it is instead solved from the same constraint as an
SPH smoothing length,
with the sum running over the SPH neighbours of the point particle only. The softening then shrinks when gas gathers around a black hole and grows when the gas leaves. Because \(h_i\) now depends on the positions of those neighbours, the gravity between the point particle and each of them gets the same kind of correction term as for an SPH particle (equations A14 and A15 of Gaburov et al. 2010), and energy stays conserved. There are no hydrodynamic forces on or from the point particle.
With dynhco=2 the solution \(h_{\rm dyn}\) is smoothly limited to
the range from hcomin to hcomax,
with \(p\) set by hcopnorm. The correction term includes the factor
\(dh/dh_{\rm dyn}\), so energy is still conserved. A hard clamp such as
\(h=\max(h_{\rm dyn}, h_{\rm comin})\) would not conserve energy, because
the correction force would switch on and off abruptly. Without a lower limit,
gas that collects around a black hole can drive the softening, and with it the
timestep, towards zero. hcomin or hcomax of zero or less removes that
limit, and neither is set by default. With neither limit set, dynhco=2
runs exactly as dynhco=1.
With dynhco=3 a smooth floor applies to every particle, SPH and point
particles alike. Each particle solves the constraint above for
\(\tilde h_i\), and its smoothing length is
with \(h_{\rm floor}\) set by hfloor and \(p\) by hcopnorm. A
point particle’s length is in addition capped at hcomax as above, so that
the gravity of a black hole is exact beyond \(2h_{\rm comax}\) while a star
is still far away. Because \(h_i\) depends on positions only through
\(\tilde h_i\), every correction term keeps the form of Gaburov et al.
(2010), with \(\nabla G\) and \(\chi_i\) evaluated at \(\tilde
h_i\), the kernel and the softened potential at \(h_i\), and \(\chi_i\)
divided by \(dh_i/d\tilde h_i\). That derivative goes to zero once
\(\tilde h_i \ll h_{\rm floor}\), so the correction terms fade out smoothly
when gas crowds into a black hole’s kernel. With the additive floor
\(h_i=\tilde h_i+h_{\rm floor}\) (hfloor with dynhco below 3) the
derivative stays 1, and the correction terms vary on the length
\(\tilde h_i\), which can become far smaller than any timestep resolves.
In a 4945-particle encounter of a 0.4 solar mass star with a 30 solar mass black
hole at periapse 0.1, gas collapsed into a clump of size \(10^{-5}\) at the
black hole and the energy jumped by more than 20 times its total with
dynhco=2, hcomin=0.05. With dynhco=3, hfloor=0.025 (with or
without hcomax=0.5) the energy changed by about \(2\times10^{-5}\), or
0.1 per cent, from periapse to five time units later, and the substeps
needed per dtmaxblk fell from 288000 to about 20000.
Blackollider
The secondary code in Blackollider/ relates \(h\) and \(\rho\)
differently. Instead of driving the neighbour-number estimate to a single
global nnopt, it solves a per-particle relation of the form
where \(a_i\), \(b_i\) and \(d_i\) are properties of the individual particle. Because \(d_i\) is added on, it acts as a floor: particle \(i\) can never have a smoothing length below it, whatever the neighbour count does.
That is what the code is for. In parallel_bleeding_edge it can be hard to
stop particles piling up around a massive compact object, and the very small
smoothing lengths that follow bring numerical trouble with them. A per-particle
lower limit on \(h\) is a convenient way to prevent that, and this relation
provides one.
For everything else, parallel_bleeding_edge and the neighbour-number
constraint above are the ones to use. They are what these pages document and
what the tutorials are written against.