21.4 Sink Particles

21.4.1 Basics of Sink Particles

Sink particles are required in collapse simulations to model dense core, star, or black hole formation and accretion. Using sink particles solves two main problems in collapse calculations:

1. The physical length scale associated with the collapse, the Jeans length,

   $$\lambda_\mathrm{J}=\sqrt{\frac{\pi c_\mathrm{s}^2}{G\rho}}$$ (21.31)

   
   
   decreases with increasing gas density $\rho$ (here, $G$ is the gravitational constant and $c_\mathrm{s}$ the sound speed). To avoid artificial fragmentation, the Jeans length must be resolved with at least 4 grid cells (Truelove et al. 1997, ApJ 489, L179). More recent calculations show that a minimum of 32 grid cells is required to resolve the kinetic energy in rotational motions inside the Jeans volume and to resolve minimum turbulent MHD dynamo action (Sur et al. 2010, ApJ 721, L134; Federrath et al. 2011, ApJ 731, 62). The FLASH sink particle module automatically refines the AMR grid based on the Jeans length (controlled by runtime parameters refineOnJeansLength, jeans_ncells_ref, and jeans_ncells_deref). However, once the highest level of the AMR hierarchy is reached, sink particles must be created to avoid artificial fragmentation when the Jeans length drops further during the collapse.
2. The typical timescale for collapsing objects is the freefall time, $t_\mathrm{ff}=\sqrt{3\pi/(32G\rho)}$. With increasing gas density, this timescale becomes shorter and shorter, which means that the code will literally grind to a halt during runaway collapse, as time steps become shorter and shorter. The first object going into freefall collapse thus determines the time step of the whole calculation, such that following the formation of a star cluster would be impossible. To avoid these problems, the gas in regions that are in a state of runaway collapse are replaced by sink particles.

21.4.2 Using the Sink Particle Unit

To include sink particles, specify REQUIRES Particles/ParticlesMain/active/Sink in the simulation Config file. This automatically includes all required subunits for the sink particles. Currently, sink particles require a Paramesh 4 Grid implementation and Poisson gravity. Currently, various subunits and implementation directories of the Particles unit and are also included automatically. Table 21.6 lists all runtime parameters of the sink particle module. The default values of these parameters should be ok in most simulations, except sink_density_thresh, sink_accretion_radius, and sink_softening_radius. Those three parameters have to be adopted to the resolution and physics of a given simulation.

First, the sink particle accretion radius, $r_\mathrm{acc}$ (sink_accretion_radius) should be calculated, based on the minimum cell size $\Delta x$ for a given simulation. The recommended value is $r_\mathrm{acc}=2.5\Delta x$. If sink particles are included, the code will print a message to the standard output, stating the number of cells that a chosen sink accretion radius corresponds to. The recommended value for the sink particle softening radius (sink_softening_radius) is to set it equal to the accretion radius.

In order to avoid artificial fragmentation, the gas density on the AMR grid must not exceed a critical density $\rho_\mathrm{thresh}$ in regions of gravitational collapse (Truelove et al. 1997). This density is related to the smallest resolvable Jeans length $\lambda_\mathrm{J}$ on the highest level of the AMR hierarchy. The density threshold $\rho_\mathrm{thresh}$ (sink_density_thresh) is obtained by solving the definition of the Jeans length, Equation (21.31) for $\rho$,

$$\rho_\mathrm{thresh} = \frac{\pi c_\mathrm{s}^2}{G \lambda_\mathrm{J}^2} = \frac{\pi c_\mathrm{s}^2}{4 G r_\mathrm{acc}^2}\;.$$ (21.32)

This equation relates the sink particle accretion radius to the sink particle density threshold and depends on the thermodynamics (sound speed) of a given simulation. Since the Jeans length should be resolved with at least 4 grid cells (Truelove et al. 1997), the sink particle accretion radius $r_\mathrm{acc}$ must not be smaller than 2 grid cells. The accretion radius is thus determined by the smallest linear cell size $\Delta x$ of the AMR grid. As recommended above, setting $r_\mathrm{acc}\simeq2.5\Delta x$ satisfies the Truelove criterion, because the Jeans length $\lambda_\mathrm{J}=2 r_\mathrm{acc}$ is thus resolved with 5 grid cells on the top level of the AMR hierarchy.

Please cite Federrath et al. (2010, ApJ 713, 269) if you are using this unit.

Table 21.6: Runtime parameters for the sink particle module.

Variable	Type	Default	Description
useSinkParticles	BOOLEAN	.false.	switch sinks on/off
sink_density_thresh	REAL	1.0e-14	density threshold for sink creation and accretion
sink_accretion_radius	REAL	1.0e14	creation and accretion radius
sink_softening_radius	REAL	1.0e14	gravitational softening radius
sink_softening_type_gas	STRING	"linear"	sink-gas softening type (options: "linear", "spline")
sink_softening_type_sinks	STRING	"spline"	sink-sink softening type (options: "linear", "spline")
sink_integrator	STRING	"leapfrog"	sink particle time integrator (options: "euler", "leapfrog", "leapfrog_cosmo")
sink_dt_factor	REAL	0.5	time step safety factor ($\leq0.5$)
sink_subdt_factor	REAL	0.01	time step safety factor for sink-sink subcycling ($\leq0.5$)
sink_convergingFlowCheck	BOOLEAN	.true.	creation check for converging gas flow
sink_potentialMinCheck	BOOLEAN	.true.	creation check for gravitational potential minimum
sink_jeansCheck	BOOLEAN	.true.	creation check for Jeans instability
sink_negativeEtotCheck	BOOLEAN	.true.	creation check for gravitationally bound gas
sink_GasAccretionChecks	BOOLEAN	.true.	check for bound and converging state before gas accretion
sink_merging	BOOLEAN	.false.	switch for sink particle merging
sink_offDomainSupport	BOOLEAN	.false.	support for sink particles to remain active when leaving the grid domain (in case of outflow boundary conditions)
sink_AdvanceSerialComputation	BOOLEAN	.true.	use the global sink particle array for time advancement (to greatly speed up computation of sink-sink interaction)
pt_maxSinksPerProc	INTEGER	100	number of sinks per processor
refineOnSinkParticles	BOOLEAN	.true.	sinks must be on highest AMR level
refineOnJeansLength	BOOLEAN	.true.	switch for refinement on Jeans length
jeans_ncells_ref	REAL	32.0	number of cells for Jeans length refinement
jeans_ncells_deref	REAL	64.0	number of cells for Jeans length de-refinement

21.4.3 The Sink Particle Method

For technical details and tests of the FLASH sink particle unit, see Federrath et al. (2010, ApJ 713, 269) and Federrath et al. (2011, IAUS 270, 425). We only summarize here the most important aspects of the implementation.

Before a sink particle is created, we define a spherical control volume with the accretion radius $r_\mathrm{acc}$ around each cell that exceeds the resolution-dependent density threshold, $\rho_\mathrm{thresh}$ given by Equation (21.32), and check whether the gas in this control volume

• is on the highest level of refinement,
• is converging ( $\nabla\cdot\mathbf{v} < 0$; i.e., has negative radial velocity),
• has a central gravitational potential minimum,
• is gravitationally bound ( $\vert E_\mathrm{grav}\vert>E_\mathrm{int}+E_\mathrm{kin}+E_\mathrm{mag}$),
• is gravitationally unstable (Jeans-unstable),
• and is not within the accretion radius of an existing sink particle.

These checks are designed to avoid spurious sink particle formation and to trace only truly collapsing and star-forming gas. As soon as a sink particle is created, it can gain further mass by accretion from the AMR grid, if this gas is inside the sink particle accretion radius and is bound and collapsing towards the sink particle. If the accretion radii of multiple sink particles overlap, the gas of a given cell is accreted onto the sink particle to which the gas is most strongly bound.

We take all contributions to the gravitational interactions between the gas on the grid and the sink particles into account. Those interactions are

1. GAS-GAS (handled by either the multigrid or the tree Poisson solver)
2. GAS-SINKS (computed by direct summation over all particles and grid cells)
3. SINKS-GAS (computed by direct summation over all particles and grid cells)
4. SINKS-SINKS (computed by direct N-body summation over all sink particles).

We note that all interactions including sink particles (GAS-SINKS, SINKS-GAS, SINKS-SINKS) do not require the Poisson solver or any particle-grid mapping procedure. Since these are all computed by direct summation, computational costs can become quite expensive once the number of sink particles reaches thousands and more, but is significantly more accurate than mapping procedures. In fact, a previous implementation used interpolation of sink particle masses back onto the grid and employing the Poisson solver to compute all interactions. This resulted in very smooth gravitational interactions that made it impossible to follow close orbits and highly eccentric encounters of multiple (two and more) sink particles, in particular for the sink-sink interactions. Thus, the current implementation of sink particles is designed to follow a maximum of a few thousand particles. Linear and spline softening of the gravitational forces for extremely close encounters of multiple sink particles are implemented. By default, linear softening is used for gas-sinks and sinks-gas interactions (sink_softening_type_gas), motivated by an almost uniform gas density inside the sink particle radius. Spline softening is used for sink-sink interactions (sink_softening_type_sinks), as it is more suitable to follow N-body dynamics. Softening is only applied inside the sink softening radius (sink_softening_radius). A second-order accurate leapfrog integrator is used to advance the sink particles. For cosmological simulations, a cosmological version of the leapfrog integrator is available (see sink_integrator). The sink particles are fully integrated into the MPI parallelization of the FLASH code. Any split or unsplit solver for hydro or MHD (in which case the bound state check includes the contribution of the magnetic energy) is supported.

21.4.4 Sink Particle Unit Test

To invoke the sink particle unit test, use ./setup unitTest/SinkMomTest -auto -3d. This initializes a momentum conservation test. An initially uniform cloud with radius $0.025 \mathrm{pc}$ and density $10^{-18} \mathrm{g} \mathrm{cm}^{-3}$ at rest starts collapsing and forms a sink particle. An initial sink particle with $0.1 M_\odot$ and initial momentum $p_y=4\times10^{36} \mathrm{g} \mathrm{cm} \mathrm{s}^{-1}$ in positive $y$-direction is also present. The purpose of this test is to see how well the total $x$-momentum $p_x$ and the total $y$-momentum $p_y$ are conserved. This test is particularly suitable, because it involves gas collapse, sink particle creation, accretion, more than one sink particle, and thus all gravitational interactions (gas-gas, gas-sinks, sinks-gas, sinks-sinks), as well as close eccentric orbits of the two sink particles. Figure 21.3 shows $p_x$ and $p_y$ as a function of time for both sink particles and gas separately. The symmetry between sink and gas momenta shows that momentum is well conserved for several orbits. Figure 21.3 can be reproduced with the IDL and Python tools in /source/Simulation/SimulationMain/unitTest/SinkMomTest/utils/ by running the IDL script plot_mom_sinks.pro.

Figure 21.3: Sink particle unit test to check momentum conservation during sink particle creation, accretion, and eccentric orbits.