doc/RTD/source/SubgridModels/GEAR/sinks/sink_timesteps.rst

+.. Sink particles in GEAR model
+   Darwin Roduit, 24 November 2024
+
+.. _sink_GEAR_timesteps:
+
+Sink timesteps
+~~~~~~~~~~~~~~
+
+Sink particles interact with the surrounding gas through accretion. To accurately follow the local gas dynamics and assess the stability of gas hydrodynamics, we must impose timestep constraints on the sink particles.
+
+First, we want to ensure sinks know the *local dynamics* and thus they obey a Courant-Friedrichs-Lewy (CFL)-like timestep constraint:
+
+.. math::
+   \Delta t_s \leq \Delta t_\text{CFL} =  C_\text{CFL} \frac{r_{\text{cut, min}}}{\sqrt{c_{s,s}^2 + \| \Delta \mathbf{v}_s \|^2}} \; ,
+
+where :math:`r_\text{cut, min} = \min(r_{\text{cut}, s}, \gamma_k \min_j(h_j))` is the minimal cut-off radius between the sink :math:`s` and the gas neighbours:math:`j`, :math:`c_{s, s}` and :math:`\Delta \mathbf{v}_s` are the gas sound-speed and the relative gas-sink velocity at the sink location. The latter two are reconstructed at the sink location with SPH interpolation. The value of :math:`C_\text{CFL}` is given in the YAML parameter file with ``GEARSink:CFL_condition``.
+
+Since sink particles accrete gas, they must anticipate the surrounding gas infall. We achieve this with a gas free-fall time criterion similar to `Grudic et al. (2021) <https://academic.oup.com/mnras/article/506/2/2199/6276745>`_:
+
+.. math::
+   \Delta t_s \leq \Delta t_\text{ff} = \sqrt{ \frac{3 \pi}{32 G \rho_s} } \quad \text{with} \quad \rho_s = \frac{3 m_s}{4 \pi {r_{\text{cut, min}}}} \; ,
+
+with :math:`m_s` the sink mass.
+
+These constraints ensure smooth gas accretion by removing a few particles per timestep. This is important since, in SPH, gas particles serve as interpolation points. By preventing the removal of a large number of particles, we ensure the stability of the hydrodynamic computations.
+
+Sink 2-body encounters
+++++++++++++++++++++++
+
+To accurately follow *sink mergers*, we implemented `Grudic et al. (2021) <https://academic.oup.com/mnras/article/506/2/2199/6276745>`_ two body timestep constraints between all sink particles:
+
+.. math::
+   \Delta t_s \leq \Delta t_\text{2-body} = \frac{ t_\text{c, min} t_\text{dyn, min}}{t_\text{c, min} + t_\text{dyn, min}} \; ,
+
+with
+
+.. math::
+  \quad t_\text{c, min} = \min_{s \neq s'} \frac{ |\varphi^{-1}(r_{ss'}, \, H_\text{sink})| }{v_{ss'}} \quad \text{and} \quad t_\text{dyn, min} = \min_{s \neq s'} \sqrt{ \frac{ |\varphi'(r_{ss'}, \, H_\text{sink})|^{-1}} { G (m_s + m_{s'})}    } \; ,
+
+where :math:`r_{ss'}` is the sinks relative separation, :math:`v_{ss'}` the relative velocity, :math:`m_{s}` and :math:`m_{s'}` their masses and :math:`H_\text{sink}` is the sink (fixed) gravitational softening. The function :math:`\varphi(r, H)` is the potential corresponding to the Wendland C2 kernel density field (see `Schaller et al. (2024) <https://doi.org/10.1093/mnras/stae922>`_ section 4.1) and  :math:`\varphi'(r, H) \equiv \frac{\mathrm{d} \varphi(r, H)}{\mathrm{d} r}` its derivative.
+
+Timesteps per sink's age categories
++++++++++++++++++++++++++++++++++++
+
+We also implemented maximal timesteps sizes depending on the sink age; :math:`\Delta t_\text{max,s}^\text{age}`. A sink can be young, old or dead. In the first two cases, the sink's timestep is :math:`\min(\Delta t_\text{max,s}^\text{age}, \Delta t_s)`. In the last case, we impose :math:`\Delta t_\text{2-body}` only if a dead sink is involved in a two-boy encounter with an alive sink. Otherwise, the sink has no timestep constraint (apart from gravity). The parameters controlling the transition between the young and old is ``GEARSink:timestep_age_threshold``, and the one between old and dead is ``GEARSink:timestep_age_threshold_unlimited_Myr``. The maximal timesteps are given by  ``GEARSink:max_timestep_young_Myr`` and  ``GEARSink:max_timestep_old_Myr``.
+
+Notice that sink particles also satisfy a gravity timestep constraint, as do all gravitational particles in Swift.
+
+Star formation constraint
++++++++++++++++++++++++++
+
+Although the stars' masses are sampled from the IMF, the stars' metallicities are not. If a sink accretes mass, it can create many star particles simultaneously. However, these stars will all have the same metallicity, which does not represent the actual metals' evolution during the accretion. In such a situation, the galaxies' properties are affected and do not represent the underlying physics.
+
+Another problem is that we can spawn many stars simultaneously, and the code may complain. Such a situation could be better. Although the constraints in the previous section will help, more is needed. Our solution is to introduce a new accretion criterion using the IMF properties. However, since our politics is that accretion should be feedback-regulated and not based on an arbitrary accretion rate, we reduce the sink time step to avoid limiting the star formation rate to an arbitrary value.
+
+The new accretion criterion is the following. The swallowed gas and sink mass does not exceed ``n_IMF`` times the IMF mass (see the IMF sampling section), but make sure to swallow at least one particle: :math:`M_\text{swallowed} \leq n_\text{IMF} M_\text{IMF} \text{ or } M_\text{swallowed} = 0`.
+
+Since we artificially restrict mass accretion, we keep track of the mass :math:`M_\text{eligible}` that would be swallowed without this criterion. Then, we compute the error :math:`\Delta M` between the restricted and unrestricted swallow. The absolute error is :math:`\Delta M = M_\text{swallowed} - M_\text{eligible}` and the relative error is :math:`| \Delta M | / M_\text{eligible}`.
+
+When :math:`\Delta M < 0` (i.e. :math:`M_\text{swallowed} \neq M_\text{eligible}`), we know the accretion was restricted and we can apply another time-step contraint. To compute a timestep, we convert :math:`\Delta M` to accretion rate by dividing by :math:`\Delta t_\text{s, estimated} = \min(\Delta t_\text{CFL}, \, \Delta t_\text{ff}, \Delta  t_\text{2-body})`. Hence, we have the constraint:
+
+.. math::
+   \Delta t_s \leq \Delta t_\text{SF} = \eta \cfrac{M_\text{eligible} \Delta t_\text{s, estimated}}{\Delta M} \text{ if } \Delta M < 0 \; ,
+
+where :math:`\eta` is a tolerance parameter. This parameter corresponds to ``GEARSink:tolerance_SF_timestep`` in the code.
+

doc/RTD/source/SubgridModels/GEAR/sinks/theory.rst

+.. Sink particles in GEAR model
+   Darwin Roduit, 24 November 2024
+
+.. _sink_GEAR_model_summary:
+
+Model summary
+-------------
+
+.. figure:: sink_scheme.png
+    :width: 400px
+    :align: center
+    :figclass: align-center
+    :alt: Illustration of the sink scheme.
+
+    This figure illustrates the sink scheme. Eligible gas particles (blue) are converted to sink particles (orange). Then, the sink searches for eligible gas/sink particles (red edges) to swallow and finally accretes them. The final step is to spawn star particles (yellow) by sampling an IMF. These stars represent a continuous portion of the IMF or an individual star (the figure does not distinguish the two types of stars).
+
+Here, we provide a comprehensive summary of the model. Sink particles are an alternative to the current model of star formation that transforms gas particles into sink particles under some criteria explained below. Then, the sink can accrete gas and spawn stars. Sink particles are collisionless particles, i.e. they interact with other particles only through gravity. They can be seen as particles representing unresolved regions of collapse. 
+
+We sample an IMF to draw the stars' mass and spawn them stochastically. Below, we provide a detailed explanation of the IMF sampling. In short, we split the IMF into two parts. In the lower part, star particles represent a continuous stellar population, similar to what is currently implemented in standard models. In the second upper part, star particles represent individual stars. Then, the feedback is improved to take into account both types of stars. Currently, only supernovae feedback is implemented. Thus, the sink particle method allows us to track the effects of individual stars' supernovae in the simulation. 
+
+The current model includes sink formation, gas accretion, sink merging, IMF sampling, star spawning and finally supernovae feedback (type Ia and II). The figures below illustrates the scheme and the associated tasks.
+
+Our main references are the following papers: `Bate et al. <https://ui.adsabs.harvard.edu/abs/1995MNRAS.277..362B/abstract>`_, `Price et al. <https://ui.adsabs.harvard.edu/abs/2018PASA...35...31P/abstract>`_ and `Federrath et al. <https://ui.adsabs.harvard.edu/abs/2010ApJ...713..269F/abstract>`_
+
+.. note::
+   Sink examples are available in ``examples/SinkParticles/``. They include self-gravity, cooling and feedback effects.
+
+.. figure:: ../../../Task/sink.png
+    :width: 400px
+    :align: center
+    :figclass: align-center
+    :alt: Task dependencies for the sink scheme.
+
+    This figure shows the task dependencies for the sink scheme.
+    The first rectangle groups the tasks that determine if sink particles will swallow other
+    sink particles or gas particles.
+    In the second one, the gas particles tagged as "to be swallowed" are effectively swallowed.
+    In the third one, the sink particles tagged as "to be swallowed" are effectively swallowed.
+    This was done with SWIFT v0.9.0.
+
+
+Conversion from comoving to physical space
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+In the following, we always refer to physical quantities. In non-cosmological simulations, there is no ambiguity between comoving and physical quantities since the universe is not expanding, and thus, the scale factor is :math:`a(t)=1`. However, in cosmological simulation, we need to convert from comoving quantities to physical ones when needed, e.g. to compute energies. We denote physical quantities by the subscript `p` and comoving ones by `c`. Here is a recap:
+
+* :math:`\mathbf{x}_p = \mathbf{x}_c a`
+* :math:`\mathbf{v}_p = \mathbf{v}_c/a + a H \mathbf{x}_c`
+* :math:`\rho_p = \rho_c/a^3`
+* :math:`\Phi_p = \Phi_c/a + c(a)`
+* :math:`u_p = u_c/a^{3(\gamma -1)}`
+
+
+Here, :math:`H` is the Hubble constant at any redshift, :math:`c(a)` is the potential normalization constant and :math:`\gamma` the gas adiabatic index. Notice that the potential normalization constant has been chosen to be :math:`c(a) = 0`.
+
+
+Sink formation
+~~~~~~~~~~~~~~
+
+.. figure:: sink_accretion_radius.png
+    :width: 400px
+    :align: center
+    :figclass: align-center
+    :alt: GEAR sink accretion radius representation
+
+    This figure shows a sink particle (in orange) newly formed among other gas particles (in blue). The accretion radius is :math:`r_{\text{acc}}`. It is the one used for sink formation. There is also an inner accretion radius :math:`f_{\text{acc}} r_{\text{acc}}` (:math:`0 \leq f_{\text{acc}} \leq 1`) that is used for gas swallowing. Particles within this inner radius are eaten without passing any other check, while particles between the two radii pass some check before being swallowed. 
+
+At the core of the sink particle method is the sink formation algorithm. This is critical to form sinks in regions adequate for star formation. Failing to can produce spurious sinks and stars, which is not desirable. However, there is no easy answer to the question. We chose to implement a simple and efficient algorithm.
+The primary criteria required to transform a gas particle into a sink are:
+
+1. the density of a given particle :math:`i` is exceeds a user-defined threshold density: :math:`\rho_i > \rho_{\text{threshold}}` ;
+2. if the particle's density lies between the threshold density and a user-defined maximal density: :math:`\rho_{\text{threshold}} \leq \rho_i \leq \rho_{\text{maximal}}`, the particle's temperature must also be below a user-defined threshold: :math:`T_i < T_{\text{threshold}}`;
+3. if the particle’s density exceeds the maximal density: :math:`\rho_i > \rho_{\text{threshold}}`, no temperature check is performed.
+
+The first criterion is common, but not the second one. We check the latter to ensure that sink particles, and thus stars, are not generated in hot regions. The third one ensures that if, for some reason, the cooling of the gas is not efficient, but the density gets very high, then we can form a sink. The parameters for those threshold quantities are respectively called ``density_threshold_Hpcm3``, ``maximal_density_threshold_Hpcm3`` and ``temperature_threshold_K``.
+
+Then, further criteria are checked. They are always checked for gas particles within the accretion radius :math:`r_{\text{acc}}` (called the ``cut_off_radius`` in the parameter file) of a given gas particle :math:`i`. Such gas particles are called *neighbours*.
+
+.. note::
+   Notice that in the current implementation, the accretion radius is kept *fixed and the same* for all sinks. However, for the sake of generality, the mathematical expressions are given as if the accretion radii could be different. 
+
+So, the other criteria are the following:
+
+3. The gas particle is at a local potential minimum: :math:`\Phi_i = \min_j \Phi_j`.
+4. Gas surrounding the particle is at rest or collapsing: :math:`\nabla \cdot \mathbf{v}_{i, p} \leq 0`. (Optional)
+5. The smoothing kernel's edge of the particle is less than the accretion radius: :math:`\gamma_k h_i < r_{\text{acc}}`, where :math:`\gamma_k` is kernel dependent. (Optional)
+6. All neighbours are currently active.
+7. The thermal energy of the neighbours satisfies: :math:`E_{\text{therm}} < |E_{\text{pot}}|/2`. (Optional, together with criterion 8.)
+8. The sum of thermal energy and rotational energy satisfies: :math:`E_{\text{therm}} + E_{\text{rot}} < | E_{\text{pot}}|`. (Optional, together with criterion 7.)
+9. The total energy of the neighbours is negative, i.e. the clump is bound to the sink: :math:`E_{\text{tot}} < 0`. (Optional)
+10. Forming a sink here will not overlap an existing sink :math:`s`: :math:`\left| \mathbf{x}_i - \mathbf{x}_s \right| > r_{\text{acc}, i} + r_{\text{acc}, s}`. (Optional)
+
+Some criteria are *optional* and can be *deactivated*. By default, they are all enabled. The different energies are computed as follows:
+
+* :math:`E_{\text{therm}} = \displaystyle \sum_j m_j u_{j, p}`
+* :math:`E_{\text{kin}} = \displaystyle \frac{1}{2} \sum_j m_j (\mathbf{v}_{i, p} - \mathbf{v}_{j, p})^2`
+* :math:`E_{\text{pot}} = \displaystyle \frac{G_N}{2} \sum_j m_i m_j \Phi_{j, p}`
+* :math:`E_{\text{rot}} = \displaystyle \sqrt{E_{\text{rot}, x}^2 + E_{\text{rot}, y}^2 + E_{\text{rot}, z}^2}`
+* :math:`E_{\text{rot}, x} = \displaystyle \frac{1}{2} \sum_j m_j \frac{L_{ij, x}^2}{\sqrt{(y_{i, p} - y_{j, p})^2 + (z_{i,p} - z_{j, p})^2}}`
+* :math:`E_{\text{rot}, y} = \displaystyle \frac{1}{2} \sum_j m_j \frac{L_{ij, y}^2}{\sqrt{(x_{i,p} - x_{j,p})^2 + (z_{i,p} - z_{j,p})^2}}`
+* :math:`E_{\text{rot}, z} = \displaystyle \frac{1}{2} \sum_j m_j \frac{L_{ij, z}^2}{\sqrt{(x_{i, p} - x_{j, p})^2 + (y_{i,p} - y_{j,p})^2}}`
+* The  (physical) specific angular momentum: :math:`\mathbf{L}_{ij} = ( \mathbf{x}_{i, p} - \mathbf{x}_{j, p}) \times ( \mathbf{v}_{i, p} - \mathbf{x}_{j, p})`
+* :math:`E_{\text{mag}} = \displaystyle \sum_j E_{\text{mag}, j}`
+* :math:`E_{\text{tot}} = E_{\text{kin}} + E_{\text{pot}} +  E_{\text{therm}} + E_{\text{mag}}`
+
+.. note::
+   Currently, magnetic energy is not included in the total energy, since the MHD scheme is in progress. However, the necessary modifications have already been taken care of.
+
+   The :math:`p` subscript is to recall that we are using physical quantities to compute energies.
+
+   Here, the potential is retrieved from the gravity solver. 
+
+
+Some comments about the criteria:
+
+The third criterion is mainly here to prevent two sink particles from forming at a distance smaller than the sink accretion radius. Since we allow sinks to merge, such a situation raises the question of which sink should swallow the other. This can depend on the order of the tasks, which is not a desirable property. As a result, this criterion is enforced.
+
+The tenth criterion prevents the formation of spurious sinks. Experiences have shown that removing gas within the accretion radius biases the hydro density estimates: the gas feels a force toward the sink. At some point, there is an equilibrium and gas particles accumulate at the edge of the accretion radius, which can then spawn sink particles that do not fall onto the primary sink and never merge. Moreover, the physical reason behind this criterion is that a sink represents a region of collapse. As a result, there is no need to have many sinks occupying the same space volume. They would compete for gas accretion without necessarily merging. This criterion is particularly meaningful in cosmological simulations to ensure proper sampling of the IMF. *This criterion can be disabled*.
+
+Once a sink is formed, we record it birth time (or scale factor in cosmological runs). This information is used to put the sink into three categories: young, old and dead. If a sink is dead, it cannot accrete gas or sink anymore. However, a dead sink can still be swallowed by a young/old sink. Young and old sink only differ by their maximal allowed timestep. Details are provided in :ref:`sink_GEAR_timesteps`.
+
+.. note::
+  However, notice that contrary to  `Bate et al. <https://ui.adsabs.harvard.edu/abs/1995MNRAS.277..362B/abstract>`_, no boundary conditions for sink particles are introduced in the hydrodynamics calculations.
+
+.. note::
+   Note that sink formation can be disabled. It can be useful, for example if you already have sinks in your initial conditions. 
+
+
+Gas accretion
+~~~~~~~~~~~~~
+
+Now that sink particles can populate the simulation, they need to swallow gas particles. To be accreted, gas particles need to pass a series of criteria. In the following, :math:`s` denotes a sink particle and :math:`i` is a gas particle. The criteria are the following:
+
+#. The sink is not dead. If it is dead, it does not accrete gas. A sink is considered dead if it is older than ``timestep_age_threshold_unlimited_Myr``.
+#. If the gas falls within :math:`f_{\text{acc}} r_{\text{acc}}` (:math:`0 \leq f_{\text{acc}} \leq 1`), the gas is accreted without further check.
+#. In the region  :math:`f_{\text{acc}} r_{\text{acc}} \leq |\mathbf{x}_i| \leq r_{\text{acc}}`, then, we check:
+   
+   #. The specific angular momentum is smaller than the one of a Keplerian orbit at :math:`r_{\text{acc}}`: :math:`|\mathbf{L}_{si}| \leq |\mathbf{L}_{\text{Kepler}}|`.
+   #. The gas is gravitationally bound to the sink particle: :math:`E_{\text{tot}} < 0`.
+   #. The gas size is smaller or equal to the sink size: :math:`\gamma_k h_i \leq r_{\text{acc}}`.
+   #. Out of all pairs of sink-gas, the gas is the most bound to this one. This case is illustrated in the figure below.
+   #. The total swallowed mass does not exceed ``n_IMF`` times the IMF mass (see the IMF sampling section), but make sure to swallow at least one particle: :math:`M_\text{swallowed} \leq n_\text{IMF} M_\text{IMF} \text{ or } M_\text{swallowed} = 0`.
+
+The physical specific angular momenta and the total energy are given by:
+
+* :math:`\mathbf{L}_{si} = ( \mathbf{x}_{s, p} - \mathbf{x}_{i, p}) \times ( \mathbf{v}_{s, p} - \mathbf{x}_{i, p})`,
+* :math:`|\mathbf{L}_{\text{Kepler}}| = r_{\text{acc}, p} \cdot \sqrt{G_N m_s / |\mathbf{x}_{s, p} - \mathbf{x}_{i, p}|^3}`.
+* :math:`E_{\text{tot}} = \frac{1}{2}  (\mathbf{v}_{s, p} - \mathbf{x}_{i, p})^2 - G_N \Phi(|\mathbf{x}_{s, p} - \mathbf{x}_{i, p}|) + m_i u_{i, p}`.
+
+.. note::
+   Here the potential is the softened potential of Swift.
+
+Those criteria are similar to `Price et al. <https://ui.adsabs.harvard.edu/abs/2018PASA...35...31P/abstract>`_ and `Grudic et al. (2021) <https://academic.oup.com/mnras/article/506/2/2199/6276745>`_, with the addition of the internal energy. This term ensures that the gas is cold enough to be accreted. Its main purpose is to avoid gas accretion and star spawning in hot regions far from sink/star-forming regions, which can happen, e.g., if a sink leaves a galaxy.
+
+Let's comment on the fourth criterion, specific to our star formation scheme. This criterion restricts the swallowed mass to avoid spawning too many stars in a single time step. Swallowing too many gas particles in a time-step can lead to instabilities in the hydrodynamics, given that gas particles act as interpolation points. Also, creating many stars at once is prejudicial for two reasons. First, the stars' mass samples an IMF, but the star's metallicities do not. So, all stars end up with the same metal content. This situation does not reflect the history of metal accretion and will lead to poor galaxy properties. Second, we need to specify at runtime the maximal number of memory allocated for extra stars until the next tree rebuild. If we create more stars than this limit, the code will stop and send an error.
+
+Since our politics is not to arbitrarily restrict the accretion using some arbitrary mass accretion rate (in fact, the accretion must be feedback-regulated), we then lower the sink time step to swallow the remaining gas particles soon. So, instead of eating a considerable amount of mass and spawning many stars in a big time step, we swallow smaller amounts of gas/sink and create fewer stars in smaller time steps. Details about the how we reduce the timestep are given in :ref:`sink_GEAR_timesteps`.
+
+Once a gas is eligible for accretion, its properties are assigned to the sink. The sink accretes the *entire* gas particle mass and its properties are updated in the following way:
+
+* :math:`\displaystyle \mathbf{v}_{s, c} = \frac{m_s \mathbf{v}_{s, c} + m_i \mathbf{v}_{i, c}}{m_s + m_i}`,
+* Swallowed physical angular momentum:  :math:`\mathbf{L}_{\text{acc}} = \mathbf{L}_{\text{acc}} + m_i( \mathbf{x}_{s, p} - \mathbf{x}_{i, p}) \times ( \mathbf{v}_{s, p} - \mathbf{x}_{i, p})`,
+* :math:`X_{Z, s} = \dfrac{X_{Z,i} m_i + X_{Z,s} m_s}{m_s + m_i}`, the metal mass fraction for each element,
+* :math:`m_s = m_s + m_i`.
+
+.. figure:: sink_overlapping.png
+    :width: 400px
+    :align: center
+    :figclass: align-center
+    :alt: Example of two sinks overlapping
+
+    This figure shows two sink particles (in orange) with gas particles (in blue) falling in the accretion radii of both sinks. In such cases, the gas particles in the overlapping regions are swallowed by the sink they are the most bound to. 
+
+Sink merging
+~~~~~~~~~~~~
+
+Sinks are allowed to merge if they enter one's accretion radius. We merge two sink particles if they respect a set of criteria. The criteria are similar to the gas particles, namely:
+
+#. At least one of the sinks is not dead. A sink is considered dead if it is older than ``timestep_age_threshold_unlimited_Myr``.
+#. If one of the sinks falls within the other's inner accretion radius, :math:`f_{\text{acc}} r_{\text{acc}}` (:math:`0 \leq f_{\text{acc}} \leq 1`), the sinks are merged without further check.
+#. In the region  :math:`f_{\text{acc}} r_{\text{acc}} \leq |\mathbf{x}_i| \leq r_{\text{acc}}`, then, we check:
+
+   #. The specific angular momentum is smaller than the one of a Keplerian orbit at :math:`r_{\text{acc}}`: :math:`|\mathbf{L}_{ss'}| \leq |\mathbf{L}_{\text{Kepler}}|`.
+   #. One sink is gravitationally bound to the other: :math:`E_{\text{mec}, ss'} < 0` or  :math:`E_{\text{mec}, s's} < 0`.
+   #. The total swallowed mass does not exceed ``n_IMF`` times the IMF mass (see the IMF sampling section), but make sure to swallow at least one particle: :math:`M_\text{swallowed} \leq n_\text{IMF} M_\text{IMF} \text{ or } M_\text{swallowed} = 0`.
+
+We compute the angular momenta and total energies in the same manner as gas particles, with the difference that we do not use internal energy. Notice that we have two energies: each sink has a different potential energy since their mass can differ.
+
+When sinks merge, the sink with the smallest mass merges with the sink with the largest. If the two sinks have the same mass, we check the sink ID number and add the smallest ID to the biggest one.
+
+IMF sampling
+~~~~~~~~~~~~
+
+.. figure:: sink_imf.png
+    :width: 400px
+    :align: center
+    :figclass: align-center
+    :alt: Initial mass function split into the continuous and discrete part.
+
+    This figure shows an IMF split into two parts by :math:`m_t`: the continuous (orange) and the discrete (blue) part. The IMF mass is :math:`M_\text{IMF} = M_c + M_d`.
+
+Now remains one critical question: how are stars formed in this scheme? Simply, by sampling an IMF. 
+In our scheme, population III stars and population II have two different IMFs. For the sake of simplicity, in the following presentation, we consider only the case of population II stars. However, this can be easily generalized to population III. 
+
+Consider an IMF such as the one above. We split it into two parts at ``minimal_discrete_mass_Msun`` (called :math:`m_t` on the illustration). The reason behind this is that we want to spawn star particles that represent *individual* (massive) stars, i.e. they are "discrete". However, for computational reasons, we cannot afford to spawn every star of the IMF as a single particle. Since the IMF is dominated by low-mass stars (< 8 :math:`M_\odot` and even smaller) that do not end up in supernovae, we would have lots of "passive" stars.
+
+.. note::
+   Recall that currently (July 2024), GEAR only implements SNIa and SNII as stellar feedback. Stars that do not undergo supernovae phases are "passive" in the current implementation.
+
+As a result, we group all those low-mass stars in one stellar particle of mass ``stellar_particle_mass_Msun``. Such star particles are called "continuous", contrary to the "discrete" individual stars.  With all that information, we can compute the number of stars in the continuous part of the IMF (called :math:`N_c`) and in the discrete part (called :math:`N_d`). Finally, we can compute the probabilities of each part, respectively called :math:`P_c` and :math:`P_d`. Notice that the mathematical derivation is given in the theory latex files.
+
+Thus, the algorithm to sample the IMF and determine the sink's ``target_mass`` is the following :
+
+* draw a random number :math:`\chi`  from a uniform distribution in the interval :math:`(0 , \; 1 ]`;
+* if  :math:`\chi < P_c`: ``sink.target_mass = stellar_particle_mass``;
+* else: ``sink_target_mass = sample_IMF_high()``.
+
+We have assumed that we have a function ``sample_IMF_high()`` that correctly samples the IMF in the discrete part.
+
+Now, what happens to the sink? After a sink forms, we give it a target mass with the abovementioned algorithm. The sink then swallows gas particles (see the task graph at the top of the page) and spawns stars.  *While the sink possesses enough mass*, we can continue to choose a new target mass. When the sink does have enough mass, the algorithm stops for this timestep. The next timestep, the sink may accrete gas and spawn stars again. The sink cannot spawn stars if it never reaches the target mass. In practice, sink particles could accumulate enough pass to spawn individual (Pop III) stars with masses 240 :math:`M_\odot` and more!
+
+For low-resolution simulations (:math:`m_\text{gas} > 100 \; M_\odot`), we also add a minimal sink mass constraint: the sink can spawn a star if ``m_sink > target_mass`` *and* ``m_sink - target_mass >= minimal_mass``. In low-resolution simulation, when a gas particle turns into a sink, the latter can have enough mass to spawn stars, depending on the sink stars and IMF parameters. As a result, the sink spawns the stars and then ends up with :math:`m_\text{sink} \ll m_\text{gas}`. Such a situation is detrimental for two reasons: 1) the sink mass is so low that gas can seldom be bound to it and thus stops spawning stars and 2) the sink can get kicked away by gravitational interactions due to the high mass difference. The parameter controlling the sink's minimal mass is ``GEARSink:sink_minimal_mass_Msun``.
+
+As explained at the beginning of this section, GEAR uses two IMFs for the population of II and III stars. The latter are called the first stars in the code. How does a sink decide which IMF to draw the target mass from? We define a threshold metallicity, ``GEARFeedback:imf_transition_metallicity`` that determines the first stars' maximal metallicity. When the sink particle's metallicity exceeds this threshold, it uses the population II IMF, defined in ``GEARFeedback:yields_table``.
+
+Star spawning
+~~~~~~~~~~~~~
+
+Once the sink spawns a star particle, we need to give properties to the star. From the sink, the star inherits the chemistry properties. The star is placed randomly within the sink's accretion radius. We draw the star's velocity components from a Gaussian distribution with mean :math:`\mu = 0` and standard deviation :math:`\sigma` determined as follows:
+
+.. math::
+   \sigma = f \cdot \sqrt{\frac{G_N M_s}{r_{\text{acc}}}} \; ,
+
+where :math:`G_N` is Newton's gravitational constant, math:`M_s` is the sink's mass before starting to spawn stars, and :math:`f` is a user-defined scaling factor. The latter corresponds to the ``star_spawning_sigma_factor`` parameter.
+
+
+Stellar feedback
+~~~~~~~~~~~~~~~~
+
+Stellar feedback *per se* is not in the sink module but in the feedback one. However, if one uses sink particles with individual stars, the feedback implementation must be adapted. Here is a recap of the GEAR feedback with sink particles. 
+
+All details and explanations about GEAR stellar feedback are provided in the GEAR :ref:`gear_stellar_evolution_and_feedback` section. Here, we only provide the changes from the previous model.
+
+In the previous model, star particles represented a population of stars with a defined IMF. Now, we have two kinds of star particles: particles representing a *continuous* portion of the IMF (see the image above) and particles representing a *single* (discrete) star. This new model requires updating the feedback model so that stars eligible for SN feedback can realise this feedback.
+
+**Discrete star particles:** Since we now have individual star particles, we can easily track SNII feedback for stars with a mass larger than 8 :math:`M_\odot`. When a star's age reaches its lifetime, it undergoes SNII feedback.
+
+**Continuous star particles**: In this case, we implemented SNII and SNIa as in the previous model. At each timestep, we determine the number of SN explosions occurring. In practice, this means that we can set the ``minimal_discrete_masss`` to any value, and the code takes care of the rest.

doc/RTD/source/SubgridModels/GEAR/supernova_feedback.rst

+.. Supernova feedback in GEAR model
+   Darwin Roduit, 30 March 2025
+
+.. gear_sn_feedback_models:
+
+.. _gear_sn_feedback_models:
+
+GEAR supernova feedback
+=======================
+
+When a star goes into a supernova, we compute the amount of internal energy, mass and metals the explosion transfers to the star's neighbouring gas particles. We will group all these in the “fluxes” term.  
+We have two models for the distribution of these fluxes and the subgrid modelling of the supernovae: GEAR model and GEAR mechanical model.
+
+.. note::
+   We may sometimes refer to GEAR feedback as GEAR thermal feedback to clearly distinguish it from GEAR mechanical feedback.
+
+
+.. _gear_sn_feedback_gear_thermal:
+
+GEAR model
+----------
+
+In the GEAR (thermal) model, the fluxes are distributed by weighing with the SPH kernel:
+
+.. math::
+
+   w_{{sj}} = W_i(\| \vec{{x}}_{{sj}} \|, \, h_s) \frac{{m_j}}{{\rho_s}}
+
+for :math:`s` the star and :math:`j` the gas (`Revaz and Jablonka 2012 <https://ui.adsabs.harvard.edu/abs/2012A%26A...538A..82R/abstract>`_).
+
+In the GEAR model, we do not inject momentum, only *internal energy*. Then, internal energy conversion to kinetic energy is left to the hydrodynamic solver, which will compute appropriately the gas density, temperature and velocity.  
+However, if the cooling radius :math:`R_{\text{cool}}` of the explosion is unresolved, i.e. the cooling radius is smaller than our simulation resolution, the cooling radiates away the internal energy.
+
+To understand why this happens, let us remind the main phases of an SN explosion in a homogeneous medium. We provide a simple picture that is more complicated than the one explained here. See `Haid et al. 2016 <https://ui.adsabs.harvard.edu/abs/2016MNRAS.460.2962H/abstract>`_ or `Thornton et al. 1998 <https://iopscience.iop.org/article/10.1086/305704>`_ for further details.
+
+* The first stage of the SN explosion is the **free expansion**. In this momentum-conserving regime, the ejecta of the stars sweeps the ISM. At the end of this phase, 72% of the initial SN energy has been converted to thermal energy.
+* Once the SN ejecta has swept an ISM mass of comparable mass, the blast wave enters the **energy-conserving Sedov-Taylor phase**. It continues with an adiabatic expansion, performing some :math:`P \, \mathrm{d}V` work on the gas. In this phase, the internal energy is converted into kinetic energy as the ejecta continues sweeping the ISM gas. This phase continues until radiative losses become significant after some radius :math:`R_{\text{cool}}`.
+* At this point, the blast wave enters the **momentum-conserving snowplough phase** and forms a thin shell. In this regime, efficient cooling radiates away the internal energy, and thus, the blast wave slows down.
+
+Now, we better understand why the internal energy is radiated away. It is a consequence of efficient cooling in the snowplough phase. When this happens, the feedback is unresolved and its energy does not affect the ISM, apart from the mass and metal injection. To circumvent this problem, GEAR thermal feedback implements a **fixed delayed cooling mechanism**. The cooling of the particles affected by feedback is deactivated during some mega year, usually 5 Myr in our simulations. The time is controlled by the ``GrackleCooling:thermal_time_myr`` parameter. This mechanism allows the internal energy to transform into kinetic energy without immediately being radiated away. However, such an approach poses the question of the time required to prevent gas from cooling in the simulations.

doc/RTD/source/SubgridModels/QuickLymanAlpha/index.rst

+.. Quick Lyman-alpha sub-grid model
+   Matthieu Schaller, 7th March 2020
+
+
+Quick Lyman-alpha (QLA) model
+=============================
+
+This section of the documentation gives a brief description of the different
+components of the quick Lyman-alpha sub-grid model. We mostly focus on the
+parameters and values output in the snapshots.
+
+Given the nature of the model, no feedback or black holes are used. The star
+formation model is minimalist and the chemistry/cooling models are limited to
+primordial abundances.
+
+.. _QLA_entropy_floors:
+
+Gas entropy floor
+~~~~~~~~~~~~~~~~~
+
+The gas particles in the QLA model are prevented from cooling below a certain
+temperature. The temperature limit depends on the density of the particles. The
+floor is implemented as a polytropic "equation of state":math:`P = P_c
+\left(\rho/\rho_c\right)^\gamma` (all done in physical coordinates), with the
+constants derived from the user input given in terms of temperature and Hydrogen
+number density. We use :math:`gamma=1` in this model. The code computing the
+entropy floor is located in the directory ``src/entropy_floor/QLA/`` and the
+floor is applied in the drift and kick operations of the hydro scheme.
+
+An additional over-density criterion above the mean baryonic density is applied
+to prevent gas not collapsed into structures from being affected. To be precise,
+this criterion demands that the floor is applied only if :math:`\rho_{\rm com} >
+\Delta_{\rm floor}\bar{\rho_b} = \Delta_{\rm floor} \Omega_b \rho_{\rm crit,0}`,
+with :math:`\Delta_{\rm floor}` specified by the user, :math:`\rho_{\rm crit,0}
+= 3H_0/8\pi G` the critical density at redshift zero [#f1]_, and
+:math:`\rho_{\rm com}` the gas co-moving density. Typical values for
+:math:`\Delta_{\rm floor}` are of order 10.
+
+The model is governed by 3 parameters for each of the two limits. These are
+given in the ``QLAEntropyFloor`` section of the YAML file. The parameters are
+the Hydrogen number density (in :math:`cm^{-3}`) and temperature (in :math:`K`)
+of the anchor point of the floor as well as the minimal over-density required to
+apply the limit. To simplify things, all constants are converted to the internal
+system of units upon reading the parameter file.
+
+For a normal quick Lyman-alpha run, that section of the parameter file reads:
+
+.. code:: YAML
+
+  QLAEntropyFloor:
+    density_threshold_H_p_cm3: 0.1       # Physical density above which the entropy floor kicks in expressed in Hydrogen atoms per cm^3.
+    over_density_threshold:    10.       # Over-density above which the entropy floor can kick in.
+    temperature_norm_K:        8000      # Temperature of the entropy floor at the density threshold expressed in Kelvin.
+
+
+SWIFT will convert the temperature normalisations and Hydrogen number density
+thresholds into internal energies and densities respectively assuming a neutral
+gas with primordial abundance pattern. This implies that the floor may not be
+exactly at the position given in the YAML file if the gas has different
+properties. This is especially the case for the temperature limit which will
+often be lower than the imposed floor by a factor :math:`\frac{\mu_{\rm
+neutral}}{\mu_{ionised}} \approx \frac{1.22}{0.59} \approx 2` due to the
+different ionisation states of the gas.
+
+Recall that we additionally impose an absolute minimum temperature at all
+densities with a value provided in the :ref:`Parameters_SPH` section of the
+parameter file. This minimal temperature is typically set to 100 Kelvin.
+
+
+.. _QLA_cooling:
+     
+Gas cooling: Wiersma+2009a with fixed primordial metallicity
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The gas cooling is based on the redshift-dependent tables of `Wiersma et
+al. (2009)a <http://adsabs.harvard.edu/abs/2009MNRAS.393...99W>`_ that include
+element-by-element cooling rates for the 11 elements (`H`, `He`, `C`, `N`, `O`,
+`Ne`, `Mg`, `Si`, `S`, `Ca` and `Fe`) that dominate the total rates. The tables
+assume that the gas is in ionization equilibrium with the cosmic microwave
+background (CMB) as well as with the evolving X-ray and UV background from
+galaxies and quasars described by the model of `Haardt & Madau (2001)
+<http://adsabs.harvard.edu/abs/2001cghr.confE..64H>`_. Note that this model
+ignores *local* sources of ionization, self-shielding and non-equilibrium
+cooling/heating. The tables can be obtained from this `link
+<http://virgodb.cosma.dur.ac.uk/swift-webstorage/CoolingTables/EAGLE/coolingtables.tar.gz>`_
+which is a re-packaged version of the `original tables
+<http://www.strw.leidenuniv.nl/WSS08/>`_. The code reading and interpolating the
+table is located in the directory ``src/cooling/QLA/``.
+
+The Wiersma tables containing the cooling rates as a function of redshift,
+Hydrogen number density, Helium fraction (:math:`X_{He} / (X_{He} + X_{H})`) and
+element abundance relative to the solar abundance pattern assumed by the tables
+(see equation 4 in the original paper). Since the quick Lyman-alpha model is
+only of interest for gas outside of haloes, we can make use of primordial gas
+only. This means that the particles do not need to carry a metallicity array or
+any individual element arrays. Another optimization is to ignore the cooling
+rates of the metals in the tables.
+
+Above the redshift of Hydrogen re-ionization we use the extra table containing
+net cooling rates for gas exposed to the CMB and a UV + X-ray background at
+redshift nine truncated above 1 Rydberg. At the redshift or re-ionization, we
+additionally inject a fixed user-defined amount of energy per unit mass to all
+the gas particles.
+
+In addition to the tables we inject extra energy from Helium II re-ionization
+using a Gaussian model with a user-defined redshift for the centre, width and
+total amount of energy injected per unit mass. Additional energy is also
+injected instantaneously for Hydrogen re-ionisation to all particles (active and
+inactive) to make sure the whole Universe reaches the expected temperature
+quickly (i.e not just via the interaction with the now much stronger UV
+background).
+
+The cooling itself is performed using an implicit scheme (see the theory
+documents) which for small values of the cooling rates is solved explicitly. For
+larger values we use a bisection scheme.  The cooling rate is added to the
+calculated change in energy over time from the other dynamical equations. This
+is different from other commonly used codes in the literature where the cooling
+is done instantaneously.
+
+We note that the QLA cooling model does not impose any restriction on the
+particles' individual time-steps. The cooling takes place over the time span
+given by the other conditions (e.g the Courant condition).
+
+Finally, the cooling module also provides a function to compute the temperature
+of a given gas particle based on its density, internal energy, abundances and
+the current redshift. This temperature is the one used to compute the cooling
+rate from the tables and similarly to the cooling rates, they assume that the
+gas is in collisional equilibrium with the background radiation. The
+temperatures are, in particular, computed every time a snapshot is written and
+they are listed for every gas particle:
+
++---------------------+-------------------------------------+-----------+-------------------------------------+
+| Name                | Description                         | Units     | Comments                            |
++=====================+=====================================+===========+=====================================+
+| ``Temperatures``    | | Temperature of the gas as         | [U_T]     | | The calculation is performed      |
+|                     | | computed from the tables.         |           | | using quantities at the last      |
+|                     |                                     |           | | time-step the particle was active |
++---------------------+-------------------------------------+-----------+-------------------------------------+
+
+Note that if one is running without cooling switched on at runtime, the
+temperatures can be computed by passing the ``--temperature`` runtime flag (see
+:ref:`cmdline-options`). Note that the tables then have to be available as in
+the case with cooling switched on. 
+
+The cooling model is driven by a small number of parameter files in the
+`QLACooling` section of the YAML file. These are the re-ionization parameters
+and the path to the tables. A valid section of the YAML file looks like:
+
+.. code:: YAML
+
+   QLACooling:
+     dir_name:     /path/to/the/Wiersma/tables/directory # Absolute or relative path
+     H_reion_z:            11.5      # Redshift of Hydrogen re-ionization
+     H_reion_ev_p_H:        2.0      # Energy injected in eV per Hydrogen atom for Hydrogen re-ionization.
+     He_reion_z_centre:     3.5      # Centre of the Gaussian used for Helium re-ionization
+     He_reion_z_sigma:      0.5      # Width of the Gaussian used for Helium re-ionization
+     He_reion_ev_p_H:       2.0      # Energy injected in eV per Hydrogen atom for Helium II re-ionization.
+
+.. _QLA_star_formation:
+
+Star formation
+~~~~~~~~~~~~~~
+
+The star formation in the Quick Lyman-alpha model is very simple. Any gas
+particle with a density larger than a multiple of the critical density for
+closure is directly turned into a star. The idea is to rapidly eliminate any gas
+that is found within bound structures since we are only interested in what
+happens in the inter-galactic medium. The over-density multiple is the only
+parameter of this model.
+
+The code applying this star formation law is located in the directory
+``src/star_formation/QLA/``. 
+
+For a normal Quick Lyman-alpha run, that section of the parameter file reads:
+
+.. code:: YAML
+
+   # Quick Lyman-alpha star formation parameters
+   QLAStarFormation:
+     over_density:              1000      # The over-density above which gas particles turn into stars.                                                                                          
+
+
+.. [#f1] Recall that in a non-cosmological run the critical density is
+	 set to 0, effectively removing the over-density
+	 constraint of the floors.

doc/RTD/source/SubgridModels/index.rst

@@ -13,6 +13,9 @@ be use as an empty canvas to be copied to create additional models.
    :maxdepth: 2
    :caption: Available models:
 	      
-   Basic/index	      
    EAGLE/index
+   QuickLymanAlpha/index
    GEAR/index
+   Basic/index	     
+   AGNSpinJets/index
+   AGORA/index

doc/RTD/source/Task/adding_to_dependency_plotting_tool.rst

+.. Dependency Plotting Additions for Tasks
+   Mladen Ivkovic, Sep 2020
+
+
+.. _task_adding_to_plotting_tool:
+.. highlight:: python
+
+
+
+Adding your Task to the Dependency Plotting Tool
+================================================
+
+How to create a task dependency graph is described in :ref:`Analysis_Tools`.
+By default, it should pick up and plot the new task. However, you might want to
+customize the plotting a bit, e.g. if you're introducing a new 'family' of tasks
+or you'd like them plotted as a cluster instead of individually. In both cases, we
+need to modify the ``tools/plot_task_dependencies.py`` script.
+
+
+
+Colouring in the Task Nodes
+---------------------------
+
+First, decide on a colour. If you want to use the same colour as already existing
+task types, find its key in the ``task_colours`` dict, which is defined at the 
+top of the file. If you want to add a new colour, add it to the ``task_colours``
+dict with a new key.
+
+Then the script needs to identify the task types with which it is working.
+To do so, it will check the task names, which are generated following the scheme
+``taskname_subtaskname``, where ``taskname`` is defined in ``taskID_names`` and
+``subtasknbame`` is defined in ``subtaskID_names`` in ``task.c``. In 
+``tools/plot_task_dependencies.py``, you'll have to write a function that recognizes your 
+task by its name, like is done for example for gravity::
+
+    def task_is_gravity(name):
+        """
+        Does the task concern the gravity?
+
+        Parameters
+        ----------
+
+        name: str
+            Task name
+        """
+        if "gpart" in name:
+            return True
+        if "grav" in name:
+            return True
+    return False
+
+You'll need to add the check to the function ``get_task_colour()``::
+
+    if taskIsGravity(name):
+        colour = task_colours["gravity"]
+
+Feel free to pick out a `nice color <http://graphviz.org/doc/info/colors.html>`_ for it :)
+
+
+
+
+
+
+
+
+Adding Clusters
+---------------
+
+In certain cases it makes sense to group some tasks together, for example the self 
+and pair tasks when computing hydro densities, gradients, or forces. To do this, 
+you'll need to modify the function ``task_get_group_name`` in ``src/task.c``. The group
+is determined by the task subtype, e.g.
+
+.. code-block:: c
+
+    case task_subtype_grav:
+      strcpy(cluster, "Gravity");
+      break;
+
+But since the task type itself is also passed to the function, you could use that
+as well if you really really need to. And that's it!

doc/RTD/source/Task/adding_your_own.rst

@@ -8,16 +8,20 @@ Adding a Task
 =============
 
 First you will need to understand the dependencies between tasks
-using the file ``dependency_graph.dot`` generated by swift at the beginning of any simulation and then decide where it will fit (see :ref:`task`).
+using the file ``dependency_graph.dot`` generated by swift at the beginning of 
+any simulation and then decide where it will fit (see :ref:`task`).
 
-For the next paragraphs, let's assume that we want to implement the existing task ``cooling``.
+For the next paragraphs, let's assume that we want to implement the existing 
+task ``cooling``.
 
 Adding it to the Task List
 --------------------------
-First you will need to add it to the task list situated in ``task.h`` and ``task.c``.
+First you will need to add it to the task list situated in ``task.h`` and 
+``task.c``.
 
-In ``task.h``, you need to provide an additional entry to the enum ``task_types`` (e.g. ``task_type_cooling``).
-The last entry ``task_type_count`` should always stay at the end as it is a counter of the number of elements.
+In ``task.h``, you need to provide an additional entry to the enum 
+``task_types`` (e.g. ``task_type_cooling``). The last entry ``task_type_count`` 
+should always stay at the end as it is a counter of the number of elements.
 For example::
 
     enum task_types {
@@ -43,8 +47,9 @@ For example::
     } __attribute__((packed));
 
 
-In ``task.c``, you will find an array containing the name of each task and need to add your own (e.g. ``cooling``).
-Be careful with the order that should be the same than in the previous list.
+In ``task.c``, you will find an array containing the name of each task and need
+to add your own (e.g. ``cooling``).  Be careful with the order that should be the
+same than in the previous list.
 For example::
 
   /* Task type names. */
@@ -59,13 +64,15 @@ For example::
 Adding it to the Cells
 ----------------------
 
-Each cell contains a list to its tasks and therefore you need to provide a link for it.
+Each cell contains a list to its tasks and therefore you need to provide a link
+for it.
 
-In ``cell.h``, add a pointer to a task in the structure.
-In order to stay clean, please put the new task in the same group than the other tasks.
+In ``cell_<particle_type>.h``, add a pointer to a task in the structure. For 
+example, cooling couples to the hydro particles, so we'll be adding our task
+to ``cell_hydro.h``.
 For example::
 
-  struct cell {
+  struct cell_hydro {
     /* Lot of stuff before. */
     
     /*! Task for the cooling */
@@ -183,11 +190,16 @@ Adding your Task to the System
 ------------------------------
 
 Now the tricky part happens.
-SWIFT is able to deal automatically with the conflicts between tasks, but unfortunately cannot understand the dependencies.
+SWIFT is able to deal automatically with the conflicts between tasks, but 
+unfortunately cannot understand the dependencies.
 
-To implement your new task in the task system, you will need to modify a few functions in ``engine.c``.
+To implement your new task in the task system, you will need to modify a 
+few functions in ``engine_maketasks.c``.
 
-First, you will need to add mainly two functions: ``scheduler_addtask`` and ``scheduler_addunlocks`` in the ``engine_make_hierarchical_tasks_*`` functions (depending on the type of task you implement, you will need to write it to a different function).
+First, you will need to add mainly two functions: ``scheduler_addtask`` and 
+``scheduler_addunlocks`` in the ``engine_make_hierarchical_tasks_*`` functions
+(depending on the type of task you implement, you will need to write it to a 
+different function).
 
 In ``engine_make_hierarchical_tasks_hydro``,
 we add the task through the following call::
@@ -205,21 +217,68 @@ and the second kick cannot be done before the cooling::
 
 
 The next step is to activate your task
-in ``engine_marktasks_mapper``::
+in the relevant section of ``cell_unskip.c`` (things are split
+by type of particles the tasks act on)::
 
   else if (t->type == task_type_cooling || t->type == task_type_sourceterms) {
     if (cell_is_active_hydro(t->ci, e)) scheduler_activate(s, t);
   }
 
-Then you will need to update the estimate for the number of tasks in ``engine_estimate_nr_tasks`` by modifying ``n1`` or ``n2``.
+Then you will need to update the estimate for the number of tasks in 
+``engine_estimate_nr_tasks`` in ``engine.c`` by modifying ``n1`` or ``n2``,
+and give the task an estimate of the computational cost that it will have in 
+``scheduler_reweight`` in  ``scheduler.c``::
+
+      case task_type_cooling:
+        cost = wscale * count_i;
+        break;
+
+This activates your tasks once they've been created.
+
+
+If your task has some computational weight, i.e. does some actual computation
+on particles, you'll also need to add it to the list of task types checked for
+weights in ``partition.c:partition_gather_weights(...)``::
+
+        /* Get the cell IDs. */
+        int cid = ci - cells;
+
+        /* Different weights for different tasks. */
+        if (t->type == task_type_init_grav || t->type == task_type_ghost ||
+            ... long list of task types ...
+            add your new task type here )
+
+            do stuff
+
+And the same needs to be done in the ``check_weights(...)`` function further
+down in the same file ``partition.c``,  where the same list of task types
+is being checked for.
+
 
 Initially, the engine will need to skip the task that updates the particles.
-It is the case for the cooling, therefore you will need to add it in ``engine_skip_force_and_kick``.
+It is the case for the cooling, therefore you will need to add it in 
+``engine_skip_force_and_kick()``.
+Additionally, the tasks will be marked as 'to be skipped' once they've been
+executed during a time step, and then reactivated during the next time step if
+they need to be executed again. This way, all the created tasks can be kept and
+don't need to be recreated every time step. In order to be unskipped however, 
+you need to add the unskipping manually to ``engine_do_unskip_mapper()`` in 
+``engine_unskip.c``.
+
+
+
+Finally, you also need to initialize your new variables and pointers in 
+``space_rebuild_recycle_mapper`` in ``space_recycle.c``.
+
+
+
+
+
 
 Implementing your Task
 ----------------------
 
-The last part is situated in ``runner.c``.
+The last part is situated in ``runner_main.c``.
 
 You will need to implement a function ``runner_do_cooling``
 (do not forget to time it)::
@@ -255,8 +314,21 @@ in the switch::
     break;
 
 
+Adding your task to the analysis tools
+--------------------------------------
+
+To produce the task graphs, the analysis scripts need to know about
+the new task. You will need to edit the python file that contains the
+hardcoded data of swift: ``tools/task_plots/swift_hardcoded_data.py``. 
+You will need to add the name of the new task to the lists in there. 
+*The order of this list needs to be the same as the enum type in the 
+task.h file!*
+
+
 Finalizing your Task
 --------------------
 
-Now that you have done the easiest part, you can start debugging by implementing a test and/or an example.
-Before creating your merge request with your new task, do not forget the most funny part that consists in writing a nice and beautiful documentation ;)
+Now that you have done the easiest part, you can start debugging by implementing
+a test and/or an example. Before creating your merge request with your new task,
+do not forget the most funny part that consists in writing a nice and beautiful
+documentation ;)

doc/RTD/source/Task/adding_your_own_neighbour_loop.rst

+.. Neighbour Loop Task
+   Mladen Ivkovic Sep 2020
+
+.. _task_adding_your_own_neighbour_loop:
+.. highlight:: c
+
+Adding a Particle Interaction/Neighbour Loop Task
+=================================================
+
+There are quite a few subtle and not so subtle differences when adding tasks that include
+particle neighbour loops and/or particle interactions compared to "independent" tasks, where
+no particles interact with each other, but work is done on them individually.
+
+Particle interactions are handled on a cell basis. If only particles of one cell
+interact with each other, the task is referred to as a  ``self`` task. When particles of
+two different cells interact with each other, we call the task type ``pair``. For a 
+particle neighbour loop, one typically requires both ``self`` and ``pair`` tasks.
+Furthermore, sometimes interaction tasks which have too much of a workload can get split 
+into smaller units of work with tasks of the type ``sub_self`` and ``sub_pair``.
+
+For the next paragraphs, let's assume that we want to implement the task ``new_iact``.
+
+
+
+
+Adding it to the Task List
+--------------------------
+First you will need to add it to the task list situated in ``task.h`` and ``task.c``.
+Here we aren't adding a task **type**, as the type will always be a ``pair`` or a 
+``self``, but **subtypes**.
+
+In ``task.h``, you need to provide an additional entry to the enum ``task_subtypes`` 
+(e.g. ``task_subtype_new_iact``). The last entry ``task_subtype_count`` should always 
+stay at the end as it is a counter of the number of elements.
+For example::
+
+    enum task_subtypes {
+      task_subtype_none = 0,        task_subtype_density,           task_subtype_gradient,
+      task_subtype_force,           task_subtype_limiter,           task_subtype_grav,
+      task_subtype_external_grav,   task_subtype_tend_part,         task_subtype_tend_gpart,
+      task_subtype_tend_spart,      task_subtype_tend_sink,         task_subtype_tend_bpart,
+      task_subtype_xv,              task_subtype_rho,               task_subtype_part_swallow,
+      task_subtype_bpart_merger,    task_subtype_gpart,             task_subtype_multipole,
+      task_subtype_spart_density,   task_subtype_stars_density,     task_subtype_stars_feedback,
+      task_subtype_new_iact,        task_subtype_count
+    } __attribute__((packed));
+
+
+In ``task.c``, you will find an array containing the name of each task and need to add your own (e.g. ``new_iact``).
+Be careful with the order that should be the same than in the previous list.
+For example::
+
+    const char *subtaskID_names[task_subtype_count] = {
+            "none",         "density",      "gradient",         "force", 
+            "limiter",      "grav",         "external_grav",    "tend_part",
+            "tend_gpart",   "tend_spart",   "tend_sink",        "tend_bpart",
+            "xv",           "rho",          "part_swallow",     "bpart_merger",
+            "gpart",        "multipole",    "spart",            "stars_density",
+            "stars_feedback","sf_count",    "bpart_rho",        "sink",
+            "new_iact"
+   };
+
+
+
+Adding it to the Cells
+----------------------
+
+Each cell contains a list to its tasks and therefore you need to provide a link for it.
+
+In ``cell_<particle_type>.h``, add a pointer to a task in the structure. For 
+example, cooling couples to the hydro particles, so we'll be adding our task
+to ``cell_hydro.h``.
+
+We won't be adding just one task though, but an entire (linked) list of them, since we're
+going to need a ``self`` type task and multiple ``pair`` type tasks to have a complete
+neighbour loop. So instead of pointing to a single task, we store a struct ``link`` in
+the cell struct.  For example::
+
+  /**
+   * @brief Hydro-related cell variables.
+   */
+  struct cell_hydro {
+
+    /*! Pointer to the #part data. */
+    struct part *parts;
+
+    /*! Pointer to the #xpart data. */
+    struct xpart *xparts;
+
+    /* Lot of stuff */
+
+    /*! Task for sorting the stars again after a SF event */
+    struct task *stars_resort;
+
+    /*! My new interaction task */
+    struct link *new_iact;
+
+    /* Lot of stuff after */
+
+  };
+
+
+Adding a new Timer
+------------------
+
+As SWIFT is HPC oriented, any new task need to be optimized.
+It cannot be done without timing the function.
+
+In ``timers.h``, you will find an enum that contains all the tasks.
+You will need to add yours inside it.
+For example::
+
+  enum {
+    timer_none = 0,
+    timer_prepare,
+    timer_init,
+    timer_drift_part,
+    timer_drift_gpart,
+    timer_kick1,
+    timer_kick2,
+    timer_timestep,
+    timer_endforce,
+    timer_dosort,
+    timer_doself_density,
+    timer_doself_gradient,
+    timer_doself_force,
+    timer_dopair_density,
+    timer_dopair_gradient,
+    timer_dopair_force,
+    timer_dosub_self_density,
+    timer_dosub_self_gradient,
+    timer_dosub_self_force,
+    timer_dosub_pair_density,
+    timer_dosub_pair_gradient,
+    timer_dosub_pair_force,
+    timer_doself_subset,
+    timer_dopair_subset,
+    timer_dopair_subset_naive,
+    timer_dosub_subset,
+    timer_do_ghost,
+    timer_do_extra_ghost,
+    timer_dorecv_part,
+    timer_do_cooling,
+    timer_gettask,
+    timer_qget,
+    timer_qsteal,
+    timer_locktree,
+    timer_runners,
+    timer_step,
+    timer_cooling,
+    timer_new_iact,
+    timer_count,
+  };
+
+As for ``task.h``,
+you will need to give a name to your timer in ``timers.c``::
+
+  const char* timers_names[timer_count] = {
+    "none",
+    "prepare",
+    "init",
+    "drift_part",
+    "kick1",
+    "kick2",
+    "timestep",
+    "endforce",
+    "dosort",
+    "doself_density",
+    "doself_gradient",
+    "doself_force",
+    "dopair_density",
+    "dopair_gradient",
+    "dopair_force",
+    "dosub_self_density",
+    "dosub_self_gradient",
+    "dosub_self_force",
+    "dosub_pair_density",
+    "dosub_pair_gradient",
+    "dosub_pair_force",
+    "doself_subset",
+    "dopair_subset",
+    "dopair_subset_naive",
+    "dosub_subset",
+    "do_ghost",
+    "do_extra_ghost",
+    "dorecv_part",
+    "gettask",
+    "qget",
+    "qsteal",
+    "locktree",
+    "runners",
+    "step",
+    "cooling",
+    "new_iact",
+  };
+
+
+You can now easily time
+your functions by using::
+
+  TIMER_TIC;
+  /* Your complicated functions */
+  if (timer) TIMER_TOC(timer_new_iact);
+
+
+Adding your Task to the System
+------------------------------
+
+Now the tricky part happens.
+SWIFT is able to deal automatically with the conflicts between tasks, but unfortunately 
+cannot understand the dependencies.
+
+To implement your new task in the task system, you will need to modify a few functions 
+in ``engine_maketasks.c``.
+
+First, you will need to add mainly three functions: ``scheduler_addtask``, ``engine_addlink``,
+``scheduler_addunlocks`` in the ``engine_make_extra_hydroloop_tasks_mapper`` functions 
+(depending on the type of task you implement, you will need to write it to a different 
+function). The (hydro) particle interaction tasks are first created only for the density 
+loop, and then replicated in ``engine_make_extra_hydroloop_tasks_mapper`` for everything
+else.
+
+In ``engine_make_extra_hydroloop_tasks_mapper``, we add the task through the following 
+call::
+
+
+    struct task *t_new_iact = NULL;
+
+    /* ... lots of stuff ... */
+
+    /* Self-interaction? */
+    else if (t_type == task_type_self && t_subtype == task_subtype_density) {
+
+      /* ... lots of stuff ... */
+
+      t_new_iact = scheduler_addtask(sched, task_type_self, task_subtype_new_iact, 
+                                     flags, 0, ci, NULL);
+
+      /* Link the tasks to the cells */
+      engine_addlink(e, &ci->new_iact, t_new_iact);
+
+      /* Create the task dependencies */
+      scheduler_addunlock(sched, ci->task_that_unlocks_this_one, t_new_iact);
+      scheduler_addunlock(sched, t_new_iact, ci->task_that_will_be_unlocked_by_this_one);
+    }
+
+    /* Otherwise, pair interaction? */
+    else if (t_type == task_type_pair && t_subtype == task_subtype_density) {
+
+      /* ... lots of stuff ... */
+
+      t_new_iact = 
+        scheduler_addtask(sched, task_type_pair, task_subtype_new_iact, 
+                          flags, 0, ci, cj);
+      engine_addlink(e, &ci->new_iact, t_new_iact);
+      engine_addlink(e, &cj->new_iact, t_new_iact);
+
+      /* ... lots of stuff ... */
+
+      if (ci->nodeID == nodeID) {
+
+        /* ... lots of stuff ... */
+
+        scheduler_addunlock(sched, ci->task_that_unlocks_this_one, t_new_iact);
+        scheduler_addunlock(sched, t_new_iact, ci->task_that_will_be_unlocked_by_this_one);
+      }
+
+      if (cj->nodeID == nodeID) {
+
+        if (ci->hydro.super != cj->hydro.super) {
+
+          /* ... lots of stuff ... */
+
+          scheduler_addunlock(sched, cj->task_that_unlocks_this_one, t_new_iact);
+          scheduler_addunlock(sched, t_new_iact, cj->task_that_will_be_unlocked_by_this_one);
+
+        }
+      }
+    }
+
+    /* Otherwise, sub-self interaction? */
+    else if (t_type == task_type_sub_self &&
+             t_subtype == task_subtype_density) {
+
+        /* You need to do the same as for task_type_self above */
+    }
+
+    /* Otherwise, sub-pair interaction? */
+    else if (t_type == task_type_sub_pair &&
+             t_subtype == task_subtype_density) {
+
+        /* You need to do the same as for task_type_pair above */
+    }
+
+
+
+
+
+The next step is to activate your task in the relevant section of
+ ``cell_unskip.c`` (things are split by the type of particles the
+ tasks run on)::
+
+
+  /* Single-cell task? */
+  if (t_type == task_type_self || t_type == task_type_sub_self) {
+
+    /* ... lots of stuff ...  */
+
+    else if (t_subtype == task_subtype_new_iact) {
+      scheduler_activate(s, t);
+    }
+  }
+
+  /* Pair? */
+  else if (t_type == task_type_pair || t_type == task_type_sub_pair) {
+
+    /* ... lots of stuff ...  */
+
+    else if (t_subtype == task_subtype_new_iact) {
+      scheduler_activate(s, t);
+    }
+  }
+
+
+Then you will need to update the estimate for the number of tasks in 
+``engine_estimate_nr_tasks`` in ``engine.c`` by modifying ``n1`` or ``n2``.
+``n1`` is the expected maximal number of tasks per top-level/super cell. ``n2``
+``n2`` is the expected maximum number of tasks for all other cells, independent
+of the depth of the tree. Most likely ``n2`` won't need updating, and you will
+only need to update ``n1``. As to how to update ``n1``, you just need to count
+the number of tasks that you will be adding, e.g. 1 self + (3^3-1)/2 = 13 pair 
+tasks + 1 ghost, etc... All these numbers can be overwritten at run time by 
+the user anyway in the parameter file (``Scheduler: tasks_per_cell``).
+
+Then give the task an estimate of the computational cost that it will have in 
+``scheduler_reweight`` in  ``scheduler.c``::
+
+      case task_type_self:
+        if (t->subtype == task_subtype_grav) {
+          cost = 1.f * (wscale * gcount_i) * gcount_i;
+        /* ... lots of stuff ... */
+        else if (t->subtype == task_subtype_new_iact)
+          cost = 1.f * wscale * scount_i * count_i;
+        else
+          error("Untreated sub-type for selfs: %s",
+                subtaskID_names[t->subtype]);
+        break;
+
+Similarly, you'll need to update ``case task_type_sub_self``, ``task_type_pair``, 
+and ``task_type_sub_pair`` as well.
+
+This activates your tasks once they've been created.
+
+
+Initially, the engine will need to skip the task that updates the particles.
+If this is the case for your task, you will need to add it in 
+``engine_skip_force_and_kick``.
+Additionally, the tasks will be marked as 'to be skipped' once they've been
+executed during a time step, and then reactivated during the next time step if
+they need to be executed again. This way, all the created tasks can be kept and
+don't need to be recreated every time step. In order to be unskipped however, 
+you need to add the unskipping manually to ``engine_do_unskip_mapper()`` in 
+``engine_unskip.c``.
+
+
+Finally, you also need to initialize your new variables and pointers in 
+``space_rebuild_recycle_mapper`` in ``space_recycle.c``. Additionally, you need 
+to initialize the ``link`` structs in ``cell_clean_links`` in ``cell.c``.
+
+
+
+
+Implementing your Task
+----------------------
+
+The last part is situated in ``runner_main.c``, where the actual functions executed
+by the task are called inside the function in ``runner_main`` in the switch::
+
+    /* Different types of tasks... */
+    switch (t->type) {
+      case task_type_self:
+        if (t->subtype == task_subtype_density)
+          runner_doself1_branch_density(r, ci);
+        /* ... lots of stuff ... */
+        else if (t->subtype == task_subtype_new_iact)
+          runner_doself_branch_new_iact(r, ci, 1);
+        else
+          error("Unknown/invalid task subtype (%s).",
+                subtaskID_names[t->subtype]);
+        break;
+        
+      case task_type_pair:
+        /* ... lots of stuff ... */
+        else if (t->subtype == task_subtype_new_iact)
+          runner_dopair_branch_new_iact(r, ci, cj, 1);
+        else
+          error("Unknown/invalid task subtype (%s/%s).",
+                taskID_names[t->type], subtaskID_names[t->subtype]);
+        break;
+
+      case task_type_sub_self:
+        /* ... lots of stuff ... */
+        else if (t->subtype == task_subtype_new_iact)
+          runner_dosub_self_new_iact(r, ci, 1);
+        else
+          error("Unknown/invalid task subtype (%s/%s).",
+                taskID_names[t->type], subtaskID_names[t->subtype]);
+        break;
+
+      case task_type_sub_pair:
+        /* ... lots of stuff ... */
+        else if (t->subtype == task_subtype_new_iact)
+          runner_dosub_pair_new_iact(r, ci, cj, 1);
+        else
+          error("Unknown/invalid task subtype (%s/%s).",
+                taskID_names[t->type], subtaskID_names[t->subtype]);
+        break;
+
+
+
+The functions ``runner_doself1_branch_density``, ``runner_dopair_branch_new_iact``,
+``runner_dosub_self_new_iact``,  and ``runner_dosub_pair_new_iact`` still need to be
+implemented by you. If you only plan on doing this type of particle interaction once
+per time step, you can get away with directly implementing these functions and call 
+it a day. But if you intend to use the same kind of particle loop more than once, as 
+it's done in e.g. the hydro density and force loops, it's better to construct the
+functions using macros.
+For example, you could have a file ``runner_doiact_my_stuff.h``::
+
+    /* File runner_doiact_my_stuff.h */
+
+    #define PASTE(x, y) x##_##y
+
+    #define _DOSELF1_BRANCH_NEW(f) PASTE(runner_doself_branch, f)
+    #define DOSELF1_BRANCH_NEW _DOSELF1_BRANCH_NEW(FUNCTION)
+
+    #define _DOPAIR1_BRANCH_NEW(f) PASTE(runner_dopair_branch, f)
+    #define DOPAIR1_BRANCH_NEW _DOPAIR1_BRANCH_NEW(FUNCTION)
+
+    #define _DOSUB_PAIR1_NEW(f) PASTE(runner_dosub_pair, f)
+    #define DOSUB_PAIR1_NEW _DOSUB_PAIR1_NEW(FUNCTION)
+
+    #define _DOSUB_SELF1_NEW(f) PASTE(runner_dosub_self, f)
+    #define DOSUB_SELF1_NEW _DOSUB_SELF1_NEW(FUNCTION)
+
+    #define _IACT_NEW(f) PASTE(runner_iact, f)
+    #define IACT_NEW _IACT_NEW(FUNCTION)
+
+    void DOSELF1_BRANCH_NEW(struct runner *r, struct cell *c, int timer);
+    void DOPAIR1_BRANCH_NEW(struct runner *r, struct cell *ci, struct cell *cj, 
+                           int timer);
+
+    void DOSUB_SELF1_NEW(struct runner *r, struct cell *ci, int timer);
+    void DOSUB_PAIR1_NEW(struct runner *r, struct cell *ci, struct cell *cj,
+                           int timer);
+
+
+
+And a second file, ``runner_doiact_function_my_stuff.h``, where you define those
+functions which have been declared using the macros, e.g. ::
+
+
+    #include "runner_doiact_my_stuff.h"
+
+    void DOSELF1_BRANCH_NEW(struct runner *r, struct cell *c, int timer) {
+      /* do your stuff, call IACT_NEW(...) at some point...*/
+    }
+
+    void DOPAIR1_BRANCH_NEW(struct runner *r, struct cell *ci, struct cell *cj, int timer) {
+      /* do your stuff, call IACT_NEW(...) at some point...*/
+    }
+
+    void DOSUB_SELF1_NEW(struct runner *r, struct cell *c, int timer) {
+      /* do your stuff, call IACT_NEW(...) at some point...*/
+    }
+
+    void DOSUB_PAIR1_NEW(struct runner *r, struct cell *ci, struct cell *cj, int timer) {
+      /* do your stuff, call IACT_NEW(...) at some point...*/
+    }
+
+
+Then we also need a ``runner_doiact_my_suff.c`` file where the functions declared in
+``runner_doiact_my_suff.h`` are defined by including them with ``FUNCTION`` defined::
+
+
+    #include <config.h>
+    /* other includes too... */
+
+    /* Import the new interaction loop functions. */
+    #define FUNCTION new_iact
+    #include "runner_doiact_functions_my_stuff.h"
+    #undef FUNCTION
+
+
+
+
+Finally, we include them in ``runner_main.c`` as follows::
+
+    /* ... lots of includes and stuff ... */
+
+    /* Import new interaction loop functions. */
+    #define FUNCTION new_iact
+    #include "runner_doiact_my_suff.h"
+    #undef FUNCTION
+
+    /**
+     * @brief The #runner main thread routine.
+     *
+     * @param data A pointer to this thread's data.
+     */
+    void *runner_main(void *data) {
+        /* ... */
+    }
+
+
+The functions ``runner_doself_branch_density``, ``runner_dopair_branch_new_iact``,
+``runner_dosub_self_new_iact``,  and ``runner_dosub_pair_new_iact`` will be properly
+found and linked this way. All that's left for you to do is to write the function
+into which ``IACT_NEW`` will expand, in the above case it would be ``runner_iact_new_iact``.
+
+If you intend on using existing neighbour loops with a new different name, or plan
+on including your newly written header files multiple times with slight changes that can
+be dealt with a preprocessing macro, you can define a new macro like this::
+
+
+    /* ... lots of includes and stuff ... */
+
+    /* Import new interaction loop functions. */
+    #define FUNCTION new_iact
+    #define FUNCTION_TASK_LOOP TASK_LOOP_NEW_IACT
+    #include "runner_doiact_my_suff.h"
+    #undef FUNCTION
+    #undef FUNCTION_TASK_LOOP
+
+    /**
+     * @brief The #runner main thread routine.
+     *
+     * @param data A pointer to this thread's data.
+     */
+    void *runner_main(void *data) {
+        /* ... */
+
+    }
+
+However, in that case you will also need to give your ``TASK_LOOP_NEW_IACT`` a unique
+definition in ``src/runner.h``, e.g.::
+
+
+    /* Unique identifier of loop types */
+    #define TASK_LOOP_DENSITY 0
+    #define TASK_LOOP_GRADIENT 1
+    #define TASK_LOOP_FORCE 2
+    #define TASK_LOOP_LIMITER 3
+    #define TASK_LOOP_FEEDBACK 4
+    #define TASK_LOOP_SWALLOW 5
+    #define TASK_LOOP_SINK_FORMATION 6
+    #define TASK_LOOP_NEW_IACT 7
+
+
+
+
+Finalizing your Task
+--------------------
+
+Now that you have done the easiest part, you can start debugging by implementing a 
+test and/or an example. Before creating your merge request with your new task, do 
+not forget the most funny part that consists in writing a nice and beautiful 
+documentation ;)
+
+
+Things to Keep in Mind
+----------------------
+
+- If you are inserting a new neighbour loop in between existing loops, or want to
+  insert more than one neighbour loop, usually a new ghost task in between them is
+  also needed.
+
+- Neighbour loops may also require MPI communication tasks.

doc/RTD/source/Task/current_dependencies.rst

+.. Current task dependencies
+   Loic Hausammann, 2020
+
+
+.. _current_dependencies:
+
+
+
+Current Task Dependencies
+=========================
+
+In order to compute the physics in the correct order, SWIFT uses dependencies in between the tasks.
+In :ref:`Analysis_Tools`, we describe our tool to generate a graph of the dependencies but,
+unfortunately, the graphs tend to be too large.
+Therefore in this section we show some smaller graphs created by hand (thus they are not necessarily reflecting an actual run depending on the physics simulated).
+The task in the form of ellipses are computational tasks while the diamond are communications.
+The colors are picked depending on the type of physics computed (blue for hydro, yellow for stars, red for gravity and black for the rest).
+The gray tasks are implicit tasks that do not compute anything but are useful to simplify the task dependencies.
+The first graph shows the full graph (without AGN and sink particles) but with some tasks collapsed into a single meta-task
+(hydrodynamics, gravity and stellar feedback):
+
+.. figure:: reduced.png
+    :width: 400px
+    :align: center
+    :figclass: align-center
+    :alt: Task dependencies for a run with gravity, hydro, cooling, star formation and stellar feedback.
+
+    This figure shows the task dependencies for a run with gravity, hydro, cooling, star formation and stellar feedback.
+    The tasks for the limiter, hydrodynamics, stellar feedback and gravity are collapsed into a single meta-task.
+    The other graphs are showing the details of the collapsed tasks except for the limiter that is done in the same way than the density loop.
+    The first tasks to be executed are at the top (without any incoming links) and then follow the order of the links
+    until the last tasks without any outgoing links.
+    This was done with SWIFT v0.9.0.
+
+As the hydrodynamics are described in :ref:`hydro`, we are only showing the gravity and stellar feedback here:
+
+.. figure:: grav.png
+    :width: 400px
+    :align: center
+    :figclass: align-center
+    :alt: Task dependencies for the gravity.
+
+    This figure shows the task dependencies for the gravity.
+    The rectangle represents the computation of the forces through either the direct computation or the multipole expansion.
+    This was done with SWIFT v0.9.0.
+
+.. figure:: stars.png
+    :width: 400px
+    :align: center
+    :figclass: align-center
+    :alt: Task dependencies for the stellar feedback.
+
+    This figure shows the task dependencies for the stellar feedback.
+    The first rectangle groups the tasks that compute the smoothing length of the stars and
+    the second one the tasks that deposit the energy into the surrounding gas.
+    This was done with SWIFT v0.9.0.
+
+
+.. figure:: sink.png
+    :width: 400px
+    :align: center
+    :figclass: align-center
+    :alt: Task dependencies for the sink scheme.
+
+    This figure shows the task dependencies for the sink scheme.
+    The first rectangle groups the tasks that determine if sink particles will swallow other
+    sink particles or gas particles.
+    In the second one, the gas particles tagged as "to be swallowed" are effectively swallowed.
+    In the third one, the sink particles tagged as "to be swallowed" are effectively swallowed.
+    This was done with SWIFT v0.9.0.
+
+For documentation on the radiative transfer tasking system, please refer to its 
+:ref:`own page <rt_task_system>`.

doc/RTD/source/Task/grav.dot

+digraph task_dep {
+	 # Header
+	 compound=true;
+	 ratio=1.41;
+	 node[nodesep=0.15, fontsize=30, penwidth=3.];
+	 edge[fontsize=0, penwidth=0.5];
+	 ranksep=0.8;
+
+	 # Special tasks
+	 self_grav[color=red3];
+	 pair_grav[color=red3];
+	 init_grav[color=red3];
+	 init_grav_out[style=filled,fillcolor=grey90,color=red3];
+	 drift_gpart[color=red3];
+	 drift_gpart_out[style=filled,fillcolor=grey90,color=red3];
+	 kick2[color=black];
+	 send_gpart[shape=diamond,style=filled,fillcolor=azure,color=red3];
+	 recv_gpart[shape=diamond,style=filled,fillcolor=azure,color=red3];
+	 grav_long_range[color=red3];
+	 grav_mm[color=red3];
+	 grav_down_in[style=filled,fillcolor=grey90,color=red3];
+	 grav_down[color=red3];
+	 grav_end_force[color=red3];
+	 recv_tend_gpart[shape=diamond,style=filled,fillcolor=azure,color=red3];
+
+	 subgraph clusterGravity {
+		 label="";
+		 bgcolor="grey99";
+		 grav_long_range;
+		 grav_mm;
+		 pair_grav;
+		 self_grav;
+	 };
+
+
+	 # Dependencies
+	 self_grav->grav_down_in[fontcolor=red3,color=red3]
+	 pair_grav->grav_down_in[fontcolor=red3,color=red3]
+	 pair_grav->recv_tend_gpart[fontcolor=red3,color=red3]
+	 init_grav->grav_long_range[fontcolor=red3,color=red3]
+	 init_grav->init_grav_out[fontcolor=red3,color=red3]
+	 init_grav_out->self_grav[fontcolor=red3,color=red3]
+	 init_grav_out->pair_grav[fontcolor=red3,color=red3]
+	 init_grav_out->init_grav_out[fontcolor=red3,color=red3]
+	 init_grav_out->grav_mm[fontcolor=red3,color=red3]
+	 drift_gpart->drift_gpart_out[fontcolor=red3,color=red3]
+	 drift_gpart->send_gpart[fontcolor=red3,color=red3]
+	 drift_gpart_out->self_grav[fontcolor=red3,color=red3]
+	 drift_gpart_out->pair_grav[fontcolor=red3,color=red3]
+	 drift_gpart_out->drift_gpart_out[fontcolor=red3,color=red3]
+	 send_gpart->grav_down[fontcolor=red3,color=red3]
+	 recv_gpart->pair_grav[fontcolor=red3,color=red3]
+	 grav_long_range->grav_down[fontcolor=red3,color=red3]
+	 grav_mm->grav_down_in[fontcolor=red3,color=red3]
+	 grav_down_in->grav_down[fontcolor=red3,color=red3]
+	 grav_down_in->grav_down_in[fontcolor=red3,color=red3]
+	 grav_down->grav_end_force[fontcolor=red3,color=red3]
+	 grav_end_force->kick2[fontcolor=red3,color=red3]
+}
\ No newline at end of file

doc/RTD/source/Task/index.rst

@@ -17,4 +17,8 @@ Everything is described in :ref:`Analysis_Tools`.
    :maxdepth: 2
    :caption: Contents:
 
+   current_dependencies
    adding_your_own
+   adding_your_own_neighbour_loop
+   adding_to_dependency_plotting_tool
+   locks

doc/RTD/source/Task/locks.rst

+.. locks
+   Mladen Ivkovic Feb 2021
+
+.. _task_locks:
+.. highlight:: c
+
+
+Task Locks
+=============
+
+Swift doesn't only deal with dependencies, but with data conflicts as
+well.  We can't have two independent tasks working on the same data
+concurrently, even if the tasks have no dependencies between them.
+For this reason, each task locks the data it works on when it begins,
+and unlocks the data again once it's finished. Data here refers to the
+particle (gas, gravity, stars,...) content of a cell as well as of the
+hierarchy of cells in the tree at levels closer to the root than the
+cell itself.  By default, it is assumed that the task is doing work on
+hydro particles. If your task requires other particle types, you will
+need to specify that in ``src/task.c``. Suppose you have implemented a
+task with type ``task_type_new`` that requires both stars and hydro
+particles: ::
+
+
+    /**
+     * @brief Try to lock the cells associated with this task.
+     *
+     * @param t the #task.
+     */
+    int task_lock(struct task *t) {
+
+      const enum task_types type = t->type;
+      const enum task_subtypes subtype = t->subtype;
+      struct cell *ci = t->ci, *cj = t->cj;
+
+      switch (type) {
+        /* lots of stuff */
+
+        case task_type_new:
+
+          /* is the data locked already? */
+          if (ci->hydro.hold || ci->stars.hold) return 0;
+
+          /* if it's not locked already, lock first particle type (hydro)
+           * if something locked it in the meantime, exit with failure. */
+          if (cell_locktree(ci) != 0) return 0;
+
+          /* if it's not locked already, lock first particle type (stars)
+           * if something locked it in the meantime, first unlock what you
+           * locked, and then exit with failure. */
+          if (cell_slocktree(ci) != 0) {
+            cell_unlocktree(ci);
+            return 0;
+          }
+
+      /* lots of other stuff */
+      }
+
+      /* lots of other stuff */
+
+    }
+
+
+
+Similarly, don't forget to write your unlocking routines in ``task_unlock()`` !

doc/RTD/source/Task/reduced.dot

+digraph task_dep {
+	 # Header
+	 compound=true;
+	 node[nodesep=0.1, fontsize=20, penwidth=3.];
+	 edge[fontsize=0, penwidth=0.5];
+	 ranksep=0.8;
+
+	 # Special tasks
+	 hydro[color=blue3,shape=folder];
+	 limiter[color=black,shape=folder];
+	 grav[color=red3,shape=folder];
+	 stars[color=darkorange1,shape=folder];
+	 drift_part[color=blue3];
+	 drift_spart[color=darkorange1];
+	 drift_gpart[color=red3];
+	 kick2[color=black];
+	 timestep[color=black];
+	 timestep_limiter[color=black];
+	 timestep_sync[color=black];
+	 recv_limiter[shape=diamond,style=filled,fillcolor=azure,color=black];
+	 recv_tend_part[shape=diamond,style=filled,fillcolor=azure,color=blue3];
+	 recv_sf_count[shape=diamond,style=filled,fillcolor=azure,color=darkorange1];
+	 cooling[color=blue3];
+	 cooling_in[style=filled,fillcolor=grey90,color=blue3];
+	 cooling_out[style=filled,fillcolor=grey90,color=blue3];
+	 star_formation[color=blue3];
+	 kick1[color=black];
+	 recv_tend_gpart[shape=diamond,style=filled,fillcolor=azure,color=red3];
+	 recv_tend_spart[shape=diamond,style=filled,fillcolor=azure,color=darkorange1];
+	 send_tend_part[shape=diamond,style=filled,fillcolor=azure,color=blue3,rank=min];
+	 send_limiter[shape=diamond,style=filled,fillcolor=azure,color=black,rank=min];
+	 send_tend_spart[shape=diamond,style=filled,fillcolor=azure,color=darkorange1,rank=min];
+	 send_tend_gpart[shape=diamond,style=filled,fillcolor=azure,color=red3,rank=min];
+	 send_sf_count[shape=diamond,style=filled,fillcolor=azure,color=darkorange1];
+
+	 # Dependencies
+	 hydro->stars[fontcolor=blue3,color=blue3]
+	 limiter->kick1[fontcolor=black,color=black]
+	 limiter->timestep_limiter[fontcolor=black,color=black]
+	 stars->timestep_sync[fontcolor=darkorange1,color=darkorange1]
+	 hydro->recv_tend_part[fontcolor=blue3,color=blue3]
+	 grav->recv_tend_gpart[fontcolor=red3,color=red3]
+	 stars->recv_tend_spart[fontcolor=darkorange1,color=darkorange1]
+	 stars->recv_tend_spart[fontcolor=darkorange1,color=darkorange1]
+	 drift_part->hydro[fontcolor=blue3,color=blue3]
+	 drift_part->stars[fontcolor=blue3,color=blue3]
+	 drift_part->limiter[fontcolor=blue3,color=blue3]
+	 drift_spart->kick2[fontcolor=darkorange1,color=darkorange1]
+	 drift_spart->stars[fontcolor=darkorange1,color=darkorange1]
+	 drift_gpart->grav[fontcolor=red3,color=red3]
+	 hydro->cooling_in[fontcolor=blue3,color=blue3]
+	 kick2->timestep[fontcolor=black,color=black]
+	 kick2->stars[fontcolor=black,color=black]
+	 kick2->star_formation[fontcolor=black,color=black]
+	 timestep->kick1[fontcolor=black,color=black]
+	 timestep->timestep_limiter[fontcolor=black,color=black]
+	 timestep->timestep_sync[fontcolor=black,color=black]
+	 timestep->limiter[fontcolor=black,color=black]
+	 timestep->send_tend_part[fontcolor=black,color=black]
+	 timestep->send_limiter[fontcolor=black,color=black]
+	 timestep->send_tend_spart[fontcolor=black,color=black]
+	 timestep->send_tend_gpart[fontcolor=black,color=black]
+	 timestep_limiter->kick1[fontcolor=black,color=black]
+	 timestep_limiter->timestep_sync[fontcolor=black,color=black]
+	 timestep_sync->kick1[fontcolor=black,color=black]
+	 recv_limiter->limiter[fontcolor=black,color=black]
+	 recv_tend_part->limiter[fontcolor=blue3,color=blue3]
+	  recv_sf_count->stars[fontcolor=darkorange1,color=darkorange1]
+	 grav->kick2[fontcolor=red3,color=red3]
+	 cooling->cooling_out[fontcolor=blue3,color=blue3]
+	 cooling_in->cooling[fontcolor=blue3,color=blue3]
+	 cooling_out->kick2[fontcolor=blue3,color=blue3]
+	 star_formation->timestep[fontcolor=blue3,color=blue3]
+	 star_formation->stars[fontcolor=blue3,color=blue3]
+	 star_formation->send_sf_count[fontcolor=blue3,color=blue3]
+	 stars->timestep[fontcolor=darkorange1,color=darkorange1]
+
+   # style
+   timestep_limiter->send_tend_part[style=invis];
+   timestep_limiter->send_tend_gpart[style=invis];
+   timestep_limiter->send_tend_spart[style=invis];
+   timestep_limiter->send_limiter[style=invis];
+}
\ No newline at end of file

doc/RTD/source/Task/sink.dot

+digraph task_dep {
+	 # Header
+	 compound=true;
+	 ratio=0.66;
+	 node[nodesep=0.15, fontsize=18, penwidth=3.];
+	 edge[fontsize=12, penwidth=0.5];
+	 ranksep=0.8;
+
+	 # Special tasks
+	 self_sink_do_sink_swallow[color=forestgreen];
+	 self_sink_swallow[color=forestgreen];
+	 self_sink_do_gas_swallow[color=forestgreen];
+	 pair_sink_do_sink_swallow[color=forestgreen];
+	 pair_sink_swallow[color=forestgreen];
+	 pair_sink_do_gas_swallow[color=forestgreen];
+	 sub_self_sink_do_sink_swallow[color=forestgreen];
+	 sub_self_sink_swallow[color=forestgreen];
+	 sub_self_sink_do_gas_swallow[color=forestgreen];
+	 sub_pair_sink_do_sink_swallow[color=forestgreen];
+	 sub_pair_sink_swallow[color=forestgreen];
+	 sub_pair_sink_do_gas_swallow[color=forestgreen];
+	 drift_sink[color=lightseagreen];
+	 star_formation_sink[color=lightseagreen];
+	 sink_in[style=filled,fillcolor=grey90,color=lightseagreen];
+	 sink_ghost1[style=filled,fillcolor=grey90,color=lightseagreen];
+	 sink_ghost2[style=filled,fillcolor=grey90,color=lightseagreen];
+	 sink_out[style=filled,fillcolor=grey90,color=lightseagreen];
+	 sink_formation[color=lightseagreen];
+
+	 # Clusters
+	 subgraph clusterSinkAccretion {
+		 label="";
+		 bgcolor="grey99";
+		 pair_sink_do_gas_swallow;
+		 self_sink_do_gas_swallow;
+		 sub_pair_sink_do_gas_swallow;
+		 sub_self_sink_do_gas_swallow;
+	 };
+
+	 subgraph clusterSinkFormation {
+		 label="";
+		 bgcolor="grey99";
+		 pair_sink_swallow;
+		 self_sink_swallow;
+		 sub_pair_sink_swallow;
+		 sub_self_sink_swallow;
+	 };
+
+	 subgraph clusterSinkMerger {
+		 label="";
+		 bgcolor="grey99";
+		 pair_sink_do_sink_swallow;
+		 self_sink_do_sink_swallow;
+		 sub_pair_sink_do_sink_swallow;
+		 sub_self_sink_do_sink_swallow;
+	 };
+
+
+	 # Dependencies
+	 self_sink_do_sink_swallow->sink_out[fontcolor=forestgreen]
+	 self_sink_swallow->sink_ghost1[fontcolor=forestgreen]
+	 self_sink_do_gas_swallow->sink_ghost2[fontcolor=forestgreen]
+	 pair_sink_do_sink_swallow->sink_out[fontcolor=forestgreen]
+	 pair_sink_swallow->sink_ghost1[fontcolor=forestgreen]
+	 pair_sink_do_gas_swallow->sink_ghost2[fontcolor=forestgreen]
+	 sub_self_sink_do_sink_swallow->sink_out[fontcolor=forestgreen]
+	 sub_self_sink_swallow->sink_ghost1[fontcolor=forestgreen]
+	 sub_self_sink_do_gas_swallow->sink_ghost2[fontcolor=forestgreen]
+	 sub_pair_sink_do_sink_swallow->sink_out[fontcolor=forestgreen]
+	 sub_pair_sink_swallow->sink_ghost1[fontcolor=forestgreen]
+	 sub_pair_sink_do_gas_swallow->sink_ghost2[fontcolor=forestgreen]
+	 drift_part->sink_formation[fontcolor=blue3]
+	 drift_sink->kick2[fontcolor=lightseagreen]
+	 drift_sink->sink_formation[fontcolor=lightseagreen]
+	 kick2->sink_in[fontcolor=black]
+	 kick2->star_formation_sink[fontcolor=black]
+	 star_formation_sink->timestep[fontcolor=lightseagreen]
+	 sink_in->sink_formation[fontcolor=lightseagreen]
+	 sink_ghost1->self_sink_do_gas_swallow[fontcolor=lightseagreen]
+	 sink_ghost1->sub_self_sink_do_gas_swallow[fontcolor=lightseagreen]
+	 sink_ghost1->sub_pair_sink_do_gas_swallow[fontcolor=lightseagreen]
+	 sink_ghost1->pair_sink_do_gas_swallow[fontcolor=lightseagreen]
+	 sink_ghost2->self_sink_do_sink_swallow[fontcolor=lightseagreen]
+	 sink_ghost2->sub_self_sink_do_sink_swallow[fontcolor=lightseagreen]
+	 sink_ghost2->sub_pair_sink_do_sink_swallow[fontcolor=lightseagreen]
+	 sink_ghost2->pair_sink_do_sink_swallow[fontcolor=lightseagreen]
+	 sink_out->timestep[fontcolor=lightseagreen]
+	 sink_out->star_formation_sink[fontcolor=lightseagreen]
+	 sink_formation->self_sink_swallow[fontcolor=lightseagreen]
+	 sink_formation->sub_self_sink_swallow[fontcolor=lightseagreen]
+	 sink_formation->sub_pair_sink_swallow[fontcolor=lightseagreen]
+	 sink_formation->pair_sink_swallow[fontcolor=lightseagreen]
+}
\ No newline at end of file

doc/RTD/source/Task/stars.dot

+digraph task_dep {
+	 # Header
+	 compound=true;
+	 ratio=1.41;
+	 node[nodesep=0.15, fontsize=30, penwidth=3.];
+	 edge[fontsize=0, penwidth=0.5];
+	 ranksep=0.8;
+
+	 # Special tasks
+	 sort[color=blue3];
+	 self_stars_density[color=darkorange1];
+	 self_stars_feedback[color=darkorange1];
+	 pair_stars_density[color=darkorange1];
+	 pair_stars_feedback[color=darkorange1];
+	 sub_self_stars_density[color=darkorange1];
+	 sub_self_stars_feedback[color=darkorange1];
+	 sub_pair_stars_density[color=darkorange1];
+	 sub_pair_stars_feedback[color=darkorange1];
+	 drift_part[color=blue3];
+	 drift_spart[color=darkorange1];
+	 kick2[color=black];
+	 timestep[color=black];
+	 timestep_sync[color=black];
+	 send_spart[shape=diamond,style=filled,fillcolor=azure,color=darkorange1];
+	 recv_rho[shape=diamond,style=filled,fillcolor=azure,color=blue3];
+	 recv_spart[shape=diamond,style=filled,fillcolor=azure,color=darkorange1];
+	 recv_sf_count[shape=diamond,style=filled,fillcolor=azure,color=darkorange1];
+	 star_formation[color=blue3];
+	 stars_in[style=filled,fillcolor=grey90,color=darkorange1];
+	 stars_out[style=filled,fillcolor=grey90,color=darkorange1];
+	 stars_ghost[color=darkorange1];
+	 stars_sort[color=darkorange1];
+	 stars_resort[color=darkorange1];
+	 recv_tend_spart[shape=diamond,style=filled,fillcolor=azure,color=darkorange1];
+	 send_sf_count[shape=diamond,style=filled,fillcolor=azure,color=darkorange1];
+
+	 subgraph clusterStarsDensity {
+		 label="";
+		 bgcolor="grey99";
+		 pair_stars_density;
+		 self_stars_density;
+		 sub_pair_stars_density;
+		 sub_self_stars_density;
+	 };
+
+	 subgraph clusterStarsFeedback {
+		 label="";
+		 bgcolor="grey99";
+		 pair_stars_feedback;
+		 self_stars_feedback;
+		 sub_pair_stars_feedback;
+		 sub_self_stars_feedback;
+	 };
+
+
+	 # Dependencies
+	 sort->pair_stars_density[fontcolor=blue3,color=blue3]
+	 sort->sub_self_stars_density[fontcolor=blue3,color=blue3]
+	 sort->sub_pair_stars_density[fontcolor=blue3,color=blue3]
+	 self_stars_density->stars_ghost[fontcolor=darkorange1,color=darkorange1]
+	 self_stars_feedback->stars_out[fontcolor=darkorange1,color=darkorange1]
+	 self_stars_feedback->timestep_sync[fontcolor=darkorange1,color=darkorange1]
+	 pair_stars_density->stars_ghost[fontcolor=darkorange1,color=darkorange1]
+	 pair_stars_density->recv_spart[fontcolor=darkorange1,color=darkorange1]
+	 pair_stars_feedback->stars_out[fontcolor=darkorange1,color=darkorange1]
+	 pair_stars_feedback->timestep_sync[fontcolor=darkorange1,color=darkorange1]
+	 pair_stars_feedback->recv_tend_spart[fontcolor=darkorange1,color=darkorange1]
+	 sub_self_stars_density->stars_ghost[fontcolor=darkorange1,color=darkorange1]
+	 sub_self_stars_feedback->stars_out[fontcolor=darkorange1,color=darkorange1]
+	 sub_self_stars_feedback->timestep_sync[fontcolor=darkorange1,color=darkorange1]
+	 sub_pair_stars_density->stars_ghost[fontcolor=darkorange1,color=darkorange1]
+	 sub_pair_stars_density->recv_spart[fontcolor=darkorange1,color=darkorange1]
+	 sub_pair_stars_feedback->stars_out[fontcolor=darkorange1,color=darkorange1]
+	 sub_pair_stars_feedback->timestep_sync[fontcolor=darkorange1,color=darkorange1]
+	 sub_pair_stars_feedback->recv_tend_spart[fontcolor=darkorange1,color=darkorange1]
+	 drift_part->self_stars_density[fontcolor=blue3,color=blue3]
+	 drift_part->pair_stars_density[fontcolor=blue3,color=blue3]
+	 drift_part->sub_self_stars_density[fontcolor=blue3,color=blue3]
+	 drift_part->sub_pair_stars_density[fontcolor=blue3,color=blue3]
+	 drift_spart->kick2[fontcolor=darkorange1,color=darkorange1]
+	 drift_spart->self_stars_density[fontcolor=darkorange1,color=darkorange1]
+	 drift_spart->pair_stars_density[fontcolor=darkorange1,color=darkorange1]
+	 drift_spart->stars_sort[fontcolor=darkorange1,color=darkorange1]
+	 drift_spart->send_spart[fontcolor=darkorange1,color=darkorange1]
+	 drift_spart->sub_self_stars_density[fontcolor=darkorange1,color=darkorange1]
+	 drift_spart->sub_pair_stars_density[fontcolor=darkorange1,color=darkorange1]
+	 kick2->stars_in[fontcolor=black,color=black]
+	 kick2->star_formation[fontcolor=black,color=black]
+	 send_spart->stars_out[fontcolor=darkorange1,color=darkorange1]
+	 recv_rho->pair_stars_density[fontcolor=blue3,color=blue3]
+	 recv_spart->stars_sort[fontcolor=darkorange1,color=darkorange1]
+	 recv_spart->pair_stars_feedback[fontcolor=darkorange1,color=darkorange1]
+	 recv_spart->sub_pair_stars_feedback[fontcolor=darkorange1,color=darkorange1]
+	 recv_sf_count->recv_spart[fontcolor=darkorange1,color=darkorange1]
+	 star_formation->stars_resort[fontcolor=blue3,color=blue3]
+	 star_formation->send_sf_count[fontcolor=blue3,color=blue3]
+	 stars_in->self_stars_density[fontcolor=darkorange1,color=darkorange1]
+	 stars_in->pair_stars_density[fontcolor=darkorange1,color=darkorange1]
+	 stars_in->sub_self_stars_density[fontcolor=darkorange1,color=darkorange1]
+	 stars_in->sub_pair_stars_density[fontcolor=darkorange1,color=darkorange1]
+	 stars_out->timestep[fontcolor=darkorange1,color=darkorange1]
+	 stars_ghost->self_stars_feedback[fontcolor=darkorange1,color=darkorange1]
+	 stars_ghost->pair_stars_feedback[fontcolor=darkorange1,color=darkorange1]
+	 stars_ghost->send_spart[fontcolor=darkorange1,color=darkorange1]
+	 stars_ghost->sub_self_stars_feedback[fontcolor=darkorange1,color=darkorange1]
+	 stars_ghost->sub_pair_stars_feedback[fontcolor=darkorange1,color=darkorange1]
+	 stars_sort->pair_stars_density[fontcolor=darkorange1,color=darkorange1]
+	 stars_sort->pair_stars_feedback[fontcolor=darkorange1,color=darkorange1]
+	 stars_sort->sub_self_stars_density[fontcolor=darkorange1,color=darkorange1]
+	 stars_sort->sub_pair_stars_density[fontcolor=darkorange1,color=darkorange1]
+	 stars_sort->sub_pair_stars_feedback[fontcolor=darkorange1,color=darkorange1]
+	 stars_resort->stars_in[fontcolor=darkorange1,color=darkorange1]
+}
\ No newline at end of file