Sink Particles

The Sink Particle algorithm is more or less taken from ​Federrath et al although Bondi type accretion is also implemented from ​Krumholz et al.

Before each step after overlaps have been applied - each processor searches through each info structure looking for new particles. This requires data out to 4 ghost cells. The processors then gather all of the new particles and the moments of the surrounding cells. Then each processor updates all of its grids by accreting material onto the sink particles. Then the change in particle mass from each processor is accumulated and the particles new properties (mass, mom, etc…) are then updated. The sink particles are then advanced redundantly using a leap-frog integration on every processor.

There are several criteria for creating sink particles. The Truelove criteria states the following must be satisfied:

\lambda_J=\left( \frac{\pi c_s^2}{G\rho} \right )^{1/2} > ~4\Delta x $

else artificial fragmentation can occur. Thus, the first criteria is that the density and sound speed in a cell must violate the Truelove criteria:

begin{equation} \lambda_J=\left( \frac{\pi c_s^2}{G\rho} \right )^{1/2} < ~4\Delta x \end{equation}

Regions that satisfy (1) are subject to potential artificial fragmentation.

A cell that violates the truelove criteria must also have converging flow and must be located at the minimum of a potential well within a radius of 4\Delta x

In addition to the central cell violating the truelove criteria - the surrounding region must be Jeans unstable or from the virial equation 2T+U < 0 where we have ignored external pressure. This implies that:

begin{equation} 2 \int{(c_V nKT) dV} - \frac{1}{2}G\int{\rho(r') \int{\frac{\rho(r)}{(r-r')}}} < 0 \end{equation}

For an ideal gas we have c_V nKT = \frac{P}{\gamma-1} = \frac{\rho c_s^2}{(\gamma-1) \gamma} or c_V nKT = \frac{3}{2} \rho c_s^2 for an isothermal gas.

If we assume a uniform sphere of radius R=4 \Delta x we get

begin{equation} 3\rho c_s^2 V - \frac{3 G \rho^2 V^2}{5R} < 0 \rightarrow c_s^2 - \frac{G\rho V}{5R} < 0 \end{equation}

or that

begin{equation} R^2 > \frac{15 c_s^2}{4 \pi G \rho} \rightarrow R > \left (\frac{15}{4 \pi^2} \right )^{1/2} \lambda_J \sim .6164 \lambda_J < \lambda_J < 4\Delta x \equiv R \end{equation}

So if a cell satisfies the Truelove criteria - then assuming that it is surrounded by cells of similar density within a radius of 4\Delta x it will also be Jeans unstable.

Magnetic fields modify the Jeans stability / virial equation to be 2T+U+M < 0 where M=\int{\frac{B^2}{2}}

Finally the region must be gravitationally bound

begin{equation} T+U+M+K<0 \end{equation}

Accretion

Federrath accretion removes just enough density from cells within 4 zones of the particle to keep them below the truelove limit.

Bondi accretions is a bit more involved. The host cell of the particle first estimates c_{\infty} and v_{\infty} using the host cell and calculates the bondi radius from the particles mass. It then calculates the surrounding mean density using \bar{\rho}=\frac{\sum{\rho wdV}}{\sum{wdV}}, and estimates \rho_{\infty} using the solution to bondi accretion \rho_{\infty}=\frac{\bar{\rho}}{\alpha(1.2\Delta x/r_{BH})}. It then calculates the accretion rate as

\dot{M}=4\pi\rho_{\infty}r_{BH}^2 \left ( \lambda^2 c_{\infty}^2 +v_{\infty}^2 \right )^{1/2}

The accretion rate can be written as \dot{M}=\frac{4\pi r_{BH}^2}{\alpha(1.2\Delta x/r_{BH})} \left ( \lambda^2 c_{\infty}^2 +v_{\infty}^2 \right )^{1/2}\frac{\sum{\dot{\rho} wdV}}{\sum{wdV}}=\frac{\sum{d\rho dV}}{\Delta t}. The Krumohlz paper suggests the amount of accreted material from each cell should be proportional to the weight factor d \rho \propto w but it makes more sense to make it proportional to the density times the weight factor d \rho = A \rho w \Delta t . If d\rho \propto \rho w then \sum{d\rho dV}= \sum{A \rho wdV \Delta t} = A \Delta t \sum{\rho wdV} which then gives A=\frac{4\pi r_{BH}^2}{\alpha(1.2\Delta x/r_{BH})\sum{wdV}} \left ( \lambda^2 c_{\infty}^2 +v_{\infty}^2 \right )^{1/2}

This makes implementation easier since each processor only needs to calculate A for local particles and then share them with other processors. This quantity is stored in particle%AccretionRate

Outflows

Particles can launch outflows following the description in ​Federrath et al 2014

The basic problem is to apply the following outflow function, but in a discrete way that is numerically symmetric and exact.

The change in rho has the following functional form d \rho=\rho_0 \mathcal{R}(r) \Theta(\theta)

as does the change in momentum d\mathbf{p}=d \rho \mathbf{\mathcal{V}}(\theta)

Since we want the total mass to equal some fraction of the accreted mass, as well as the top and bottom mass injection to be symmetric, we must have

\displaystyle{\sum_{top}{d \rho_i}}=\displaystyle{\sum_{bottom}{d \rho_ i}} = \frac{f \dot{M}}{2}

And since we also must have the total momentum be balanced, we need to have

\displaystyle{\sum_{top}{|dp^j_i|}}+\displaystyle{\sum_{bottom}{|dp^j_i|}} = f \dot{M}\mathcal{V}_0

as well as

\displaystyle{\sum_{top}{dp^j_i}}+\displaystyle{\sum_{bottom}{dp^j_i}} = \mathbf{0}

If we introduce scaling parameters for the density and momentum for the top and bottom…

d \rho_\pm=\alpha_{\pm} \rho_0 \mathcal{R}(r) \Theta(\theta)

d\mathbf{p}_\pm=d \rho_\pm \beta_\pm \mathbf{\mathcal{V}}(\theta)

and plug these into the above equations, we can solve for \alpha_\pm and \beta_\pm

\alpha_\pm = \frac{\displaystyle{\sum_\mp{\mathcal{R}_i \Theta_i}}}{\displaystyle{\sum_+{\mathcal{R}_i \Theta_i}\sum_-{\mathcal{R}_i \Theta_i}}} \frac{f \dot{M}}{2 \rho_0}

The momentum equations give us 4 constraints for two variables.

For reference on another way astrobear has implemented sinks in the past, see SinkParticle