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

+.. Sink particles in GEAR model
+   Darwin Roduit, 14 July 2024
+
+.. sink_GEAR_model:
+
+Snapshots ouputs
+----------------
+
+Here, we provide a summary of the quantities written in the snapshots, in addition to positions, velocities, masses, smoothing lengths and particle IDs.
+
+
+Sink particles
+~~~~~~~~~~~~~~
+
++---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
+| Name                                  | Description                                 | Units                  | Comments                                          |
++=======================================+=============================================+========================+===================================================+
+| ``NumberOfSinkSwallows``              | | Number of merger events with other sinks  | [-]                    |                                                   |
++---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
+| ``NumberOfGasSwallows``               | | Number of gas particles accreted          | [-]                    |                                                   |
++---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
+| ``TargetMass``                        | | Target mass required to spawn the next    | [U_M]                  | | You can use it to determine if the target mass  |
+|                                       | | star particle                             |                        | | is so huge that the sink's mass cannot spawn    |
+|                                       |                                             |                        | | such a star. Such rare behaviour may bias the   |
+|                                       |                                             |                        | | IMF towards high masses.                        |
++---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
+| ``Nstars``                            | | Number of stars created by this sink      | [-]                    |                                                   |
++---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
+| ``SwallowedAngularMomentum``          | | Total angular momentum of accreted        | [U_M U_L^2 U_T^{-1}]   |                                                   |
+|                                       | | material                                  |                        |                                                   |
++---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
+| ``MetalMassFractions``                | | Mass fraction of each tracked metal       | [-]                    | | *Only in GEAR chemistry module.*                |
+|                                       | | element                                   |                        | | Array of length ``N`` (number of elements),     |
+|                                       |                                             |                        | | set at compile time via                         |
+|                                       |                                             |                        | | ``--with-chemistry=GEAR_N``.                    |
++---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
+| ``BirthScaleFactors``                 | | Scale factor at the time of sink creation | [-]                    | | Only used in *cosmological* runs.               |
++---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
+| ``BirthTimes``                        | | Time when the sink was created            | [U_T]                  | | Only used in *non-cosmological* runs.           |
++---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
+
+
+
+Stars
+~~~~~
+
++---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
+| Name                                  | Description                                 | Units                  | Comments                                          |
++=======================================+=============================================+========================+===================================================+
+| ``BirthScaleFactors``                 | | Scale-factors when the stars were born    | [-]                    | | Only used in cosmological runs.                 |
++---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
+| ``BirthTimes``                        | | Time when the stars were born             | [U_T]                  | | Only used in non-cosmological runs.             |
++---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
+| ``BirthMasses``                       | | Masses of the stars at birth time         | [U_M]                  | | SF and sinks modules                            |
++---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
+| ``ProgenitorIDs``                     | | ID of the progenitor sinks or gas         | [-]                    | | SF and sinks modules                            |
+|                                       | | particles                                 |                        |                                                   |
++---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
+| ``BirthDensities``                    | | Gas density at star formation             | [U_M U_L^{-3}]         | | *Only in SF module*                             |
++---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
+| ``BirthTemperatures``                 | | Gas temperature at star formation         | [K]                    | | *Only in SF module*                             |
++---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
+| ``Potentials``                        | | Gravitational potential of the star       | [U_L^2 U_T^{-2}]       |                                                   |
++---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
+| ``StellarParticleType``               | | Type of the stellar particle:             | [-]                    | | 0: (discrete) single star                       |
+|                                       | |                                           |                        | | 1: continuous IMF part star                     |
+|                                       | |                                           |                        | | 2: single population star                       |
+|                                       | |                                           |                        | | The last type corresponds to legacy IMF stars.  |
++---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
+| ``MetalMassFractions``                | | Mass fraction of each metal element       | [-]                    | | *Only in GEAR chemistry module*.                |
+|                                       | |                                           |                        | | Array of length ``N`` (number of elements),     |
+|                                       | |                                           |                        | | set at compile time by                          |
+|                                       | |                                           |                        | | ``--with-chemistry=GEAR_N``.                    |
++---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
+
+Gas particles
+~~~~~~~~~~~~~
+
+Since hydro scheme writes its own set of outputs, we only provide the outputs that ``GEAR`` writes for gas particles. 
+
++------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
+| Name                                     | Description                                                 | Units                  | Comments                                           |
++==========================================+=============================================================+========================+====================================================+
+| ``SmoothedMetalMassFractions``           | | Mass fraction of each metal element                       | [-]                    | | *Only in GEAR chemistry module.*                 |
+|                                          | | smoothed over the SPH kernel                              |                        | | Array of length ``N``, set at compile time by    |
+|                                          |                                                             |                        | | ``--with-chemistry=GEAR_N``                      |
++------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
+| ``MetalMassFractions``                   | | Raw (non-smoothed) mass fraction of                       | [-]                    | | *Only in GEAR chemistry module.*                 |
+|                                          | | each metal element                                        |                        | | Same layout as above.                            |
++------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
+| ``HI``                                   | | Mass fraction of neutral H (:math:`\mathrm{H}`)           | [-]                    | | *Only if* ``GRACKLE_1 to 3``                     |
++------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
+| ``HII``                                  | | Mass fraction of ionized H (:math:`\mathrm{H}^+`)         | [-]                    | | *Only if* ``GRACKLE_1 to 3``                     |
++------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
+| ``HeI``                                  | | Mass fraction of neutral He (:math:`\mathrm{He}`)         | [-]                    | | *Only if* ``GRACKLE_1 to 3``                     |
++------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
+| ``HeII``                                 | | Mass fraction of singly ionized He (:math:`\mathrm{He}^+`)| [-]                    | | *Only if* ``GRACKLE_1 to 3``                     |
++------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
+| ``HeIII``                                | | Mass fraction of doubly ionized He                        | [-]                    | | *Only if* ``GRACKLE_1 to 3``                     |
+|                                          | | (:math:`\mathrm{He}^{++}`)                                |                        |                                                    |
++------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
+| ``e``                                    | | Free electron mass fraction (:math:`\mathrm{e}^-`)        | [-]                    | | *Only if* ``GRACKLE_1 to 3``                     |
++------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
+| ``HM``                                   | | Mass fraction of :math:`\mathrm{H}^-`                     | [-]                    | | *Only if* ``GRACKLE_2 or 3``                     |
++------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
+| ``H2I``                                  | | Mass fraction of neutral :math:`\mathrm{H}_2`             | [-]                    | | *Only if* ``GRACKLE_2 or 3``                     |
++------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
+| ``H2II``                                 | | Mass fraction of ionized :math:`\mathrm{H}_2^+`           | [-]                    | | *Only if* ``GRACKLE_2 or 3``                     |
++------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
+| ``DI``                                   | | Mass fraction of neutral D (:math:`\mathrm{D}`)           | [-]                    | | *Only if* ``GRACKLE_3``                          |
++------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
+| ``DII``                                  | | Mass fraction of ionized D (:math:`\mathrm{D}^+`)         | [-]                    | | *Only if* ``GRACKLE_3``                          |
++------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
+| ``HDI``                                  | | Mass fraction of :math:`\mathrm{HD}`                      | [-]                    | | *Only if* ``GRACKLE_3``                          |
++------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+

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

+.. Sink particles in GEAR model
+   Darwin Roduit, 17 April 2024
+
+.. sink_GEAR_model:
+
+.. _sink_GEAR_model:
+
+Sink particles and star formation in GEAR model
+===============================================
+
+.. toctree::
+   :caption: Table of Contents
+
+   introduction
+   theory
+   sink_timesteps
+   params

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

+.. Sink particles in GEAR model
+   Darwin Roduit, 15 March 2024
+
+.. sink_GEAR_model:
+
+Introduction
+------------
+
+GEAR sink particles provide an alternative model to star formation. Instead of stochastically transforming gas particles into stars as is done in GEAR star formation scheme under some conditions, we transform a gas into a sink particle. The main property of the sink particle is its accretion radius. When gas particles within this accretion radius are eligible to be swallowed by the sink, we remove them and transfer their mass, momentum, angular momentum, chemistry properties, etc to the sink particle.
+
+With the sink particles, the IMF splits into two parts: the continuous part and the discrete part. Those parts will correspond to two kinds of stars. Particles in the discrete part of the IMF represent individual stars. It means that discrete IMF-sampled stars have different masses. Particles in the continuous part represent a population of stars, all with the same mass.
+
+The sink particle will randomly choose a target mass, accrete gas until it reaches this target mass and finally spawn a star. Then, the sink chooses a new target mass and repeats the same procedure. When stars are spawned, they are given a new position and velocity as well as chemical properties.
+
+In ``theory.rst`` we outline all of the theory which is implemented as part of the model. This includes how sink particles are formed, how the gas is accreted, how sinks are merged and how stars are spawned. In ``params.rst`` we list and discuss all parameters used by the model. Below we outline how to configure and run the model.
+
+Compiling and running the model
+-------------------------------
+
+You can configure the model with ``--with-sink=GEAR`` in combination with other configure options of the GEAR model. The model will then be used when the ``--sinks`` flag is among the runtime options. In particular, the sink particles require ``--feedback`` runtime option.
+
+Notice that you also need to compile with ``--with-star-formation=GEAR``. The star formation module is required to collect and write star formation data. However, you do not need to run Swift with the option ``--star-formation``.
+
+Then, you do not need to do anything special. Sink particles will be created during your runs. If you want, you can have sink particles in your ICs. At the moment, sink particles do not have any special fields to be set.
+
+A full list of all relevant parameters of the model is in :ref:`sink_GEAR_parameters`. We also briefly describe the most important parameters which need to be set to run the model, as well as how to run it in different configurations.
+
+.. warning::
+   Currently, MPI is not implemented for the sink particles. If you try to run with MPI, you will encounter an error. We thus recommend configuring Swift with ``--disable-mpi`` to avoid any surprises.

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

+.. Sink particles in GEAR model
+   Darwin Roduit, 15 March 2024
+
+.. sink_GEAR_model:
+
+.. _sink_GEAR_parameters:
+
+Model parameters
+----------------
+
+The parameters of the GEAR sink model are grouped into the ``GEARSink`` section of the parameter file. 
+
+The first three parameters are:
+
+* Whether we are using a fixed cut-off radius for gas and sink accretion: ``use_fixed_cut_off_radius``,
+* The sink cut-off radius: ``cut_off_radius``,
+* The sink inner accretion radius fraction in terms of the cut-off radius: ``f_acc``.
+
+The ``use_fixed_cut_off_radius`` is mandatory and should be set to 0 or 1. If set to 1, the GEAR model will use a fixed cutoff radius equal to the value of ``cut_off_radius``. If not, the cutoff radius is allowed to vary according to the local gas density. In the code, the cutoff radius is always equal to the sink's smoothing length multiplied by a constant factor ``kernel_gamma`` - setting a fixed cutoff will fix the smoothing length at the appropriate value for all sinks.
+
+The ``cut_off_radius`` parameter is optional, and is ignored if ``use_fixed_cut_off_radius`` is 0. If a fixed cut-off is used, every sink's smoothing length will be permanently set to this number divided by ``kernel_gamma`` on formation.
+
+The ``f_acc`` parameter is also optional. Its default value is :math:`0.1`. Its value must respect :math:`0 \leq f_\text{acc} \leq 1` . It describes the inner radius :math:`f_{\text{acc}} \cdot r_{\text{cut-off}}` in which gas particles are swallowed without any further checks, as explained below.
+
+The next three mandatory parameters are:
+
+* the gas temperature threshold to form a sink when :math:`\rho_\text{threshold} < \rho_\text{gas} < \rho_\text{maximal}` :``temperature_threshold_K``,
+* the minimal gas threshold density required to form a sink:``density_threshold_Hpcm3``,
+* the maximal density at which the temperature check is not performed:``maximal_density_threshold_Hpcm3`` (Default: ``FLT_MAX``).
+
+These three parameters govern the first two criteria of the sink formation scheme. If these criteria are not passed, sink particles are not created. If they are passed, the code performs further checks to form sink particles. Some of those criteria checks can be disabled, as explained below.
+
+The next set of parameters deals with the sampling of the IMF and the representation of star particles:
+
+* minimal mass of stars represented by discrete particles: ``minimal_discrete_mass_Msun``,
+* mass of the stellar particle representing the continuous part of the IMF: ``stellar_particle_mass_Msun``,
+* minimal mass of the first stars represented by discrete particles: ``minimal_discrete_mass_first_stars_Msun``,
+* mass of the stellar particle (first stars) representing the continuous part of the IMF:: ``stellar_particle_mass_first_stars_Msun``.
+
+With sink particles, star particles can represent either a single star or a population of stars in the low mass part of the IMF (continuous IMF sampling). The stars in the continuous part of the IMF are put together in a particle of mass ``stellar_particle_mass_Msun`` or ``stellar_particle_mass_first_stars_Msun``, while individual stars in the discrete part have their mass sampled from the IMF. The limit between the continuous and discrete sampling of the IMF is controlled by  ``minimal_discrete_mass_Msun`` and ``minimal_discrete_mass_first_stars_Msun``.
+
+The next set of parameters controls the sink formation scheme. More details are provided in the GEAR documentation. Here is a brief overview:
+
+* whether or not the gas must be contracting: ``sink_formation_contracting_gas_criterion`` (default: 1), 
+* whether or not the gas smoothing length must be small enough: ``sink_formation_smoothing_length_criterion`` (default: 1),
+* whether or not the gas must be in a Jeans unstable state: ``sink_formation_jeans_instability_criterion`` (default: 1),
+* whether or not the gas must be in a bound state: ``sink_formation_bound_state_criterion`` (default: 1),
+* whether or not a new sink can be formed in a region where its ``cut_off_radius`` and the one of an existing sink overlap: ``sink_formation_overlapping_sink_criterion`` (default: 1).
+
+Those criteria are checked if the density and temperature criteria are successfully passed. They control the behaviour of the sink formation scheme. By default, they are all activated and set to ``1``.
+
+The next parameter is ``disable_sink_formation`` (default: 0). It controls whether sinks are formed or not in the simulation. The main purpose is when sinks are put in initial conditions and sinks are not wanted to be added during the run. This parameter is set to ``0`` by default, i.e. sink formation is *enabled*.
+
+The next set of parameters deals with the sink time-steps:
+
+* Courant-Friedrich-Levy constant for the CFL-like time-step constraint:``CFL_condition``,
+* age (in Myr) at which a sink is considered dead (no accretion) and without time-step limitations, except for 2-body encounters involving another young/old sink and gravity: ``timestep_age_threshold_unlimited_Myr`` (default: FLT_MAX),
+* age (in Myr) at which sinks switch from young to old for time-stepping purposes:``timestep_age_threshold``  (default: FLT_MAX),
+* maximal time-step length of young sinks (in Myr): ``max_timestep_young_Myr``  (default: FLT_MAX),
+* maximal time-step length of old sinks (in Myr): ``max_timestep_old_Myr`` (default: FLT_MAX),
+* number of times the IMF mass can be swallowed in a single time-step: ``n_IMF`` (default: FLT_MAX).
+* tolerance parameter for SF timestep constraint: ``tolerance_SF_timestep`` (default: 0.5)
+
+The last parameter is ``sink_minimal_mass_Msun``. This parameter is mainly intended for low-resolution simulations with :math:`m_\text{gas} > 100 \; M_\odot`. It prevents :math:`m_\text{sink} \ll m_\text{gas}` simulations when sinks spawn stars, which can lead to gravity run away.
+
+The full section is:
+
+.. code:: YAML
+
+   GEARSink:
+     use_fixed_cut_off_radius: 1                 # Are we using a fixed cutoff radius? If we are, in GEAR the cutoff radius is fixed at the value specified below, and the sink smoothing length is fixed at this value divided by kernel_gamma. If not, the cutoff radius varies with the sink smoothing length as r_cut = h*kernel_gamma.
+     cut_off_radius: 1e-3                        # Cut off radius of all the sinks in internal units. Ignored if use_fixed_cut_off_radius is 0. 
+     f_acc: 0.1                                  # (Optional) Fraction of the cut_off_radius that determines if a gas particle should be swallowed wihtout additional check. It has to respect 0 <= f_acc <= 1. (Default: 0.1)
+     temperature_threshold_K:        100         # Max temperature (in K) for forming a sink when density_threshold_Hpcm3 <= density <= maximal_density_threshold_Hpcm3.
+     density_threshold_Hpcm3: 1e3                # Minimum gas density (in Hydrogen atoms/cm3) required to form a sink particle.
+     maximal_density_threshold_Hpcm3: 1e5        # If the gas density exceeds this value (in Hydrogen atoms/cm3), a sink forms regardless of temperature if all other criteria are passed. (Default: FLT_MAX)
+     sink_minimal_mass_Msun:     0.              # (Optional) Sink minimal mass in Msun. This parameter prevents m_sink << m_gas in low resolution simulations. (Default: 0.0)
+     stellar_particle_mass_Msun:  20             # Mass of the stellar particle representing the low mass stars (continuous IMF sampling) (in solar mass)
+     minimal_discrete_mass_Msun: 8               # Minimal mass of stars represented by discrete particles (in solar mass)
+     stellar_particle_mass_first_stars_Msun: 20      # Mass of the stellar particle representing the low mass stars (continuous IMF sampling) (in solar mass). First stars
+     minimal_discrete_mass_first_stars_Msun: 8       # Minimal mass of stars represented by discrete particles (in solar mass). First stars
+     star_spawning_sigma_factor: 0.2                 # Factor to rescale the velocity dispersion of the stars when they are spawned. (Default: 0.2)
+     sink_formation_contracting_gas_criterion: 1     # (Optional) Activate the contracting gas check for sink formation. (Default: 1)
+     sink_formation_smoothing_length_criterion: 1    # (Optional) Activate the smoothing length check for sink formation. (Default: 1)
+     sink_formation_jeans_instability_criterion: 1   # (Optional) Activate the two Jeans instability checks for sink formation. (Default: 1)
+     sink_formation_bound_state_criterion: 1         # (Optional) Activate the bound state check for sink formation. (Default: 1)
+     sink_formation_overlapping_sink_criterion: 1    # (Optional) Activate the overlapping sink check for sink formation. (Default: 1)
+     disable_sink_formation: 0                       # (Optional) Disable sink formation. (Default: 0)
+     CFL_condition:                        0.5       # Courant-Friedrich-Levy condition for time integration.
+     timestep_age_threshold_unlimited_Myr: 100.      # (Optional) Age above which sinks no longer have time-step restrictions, except for 2-body encounters involving another young/old sink and gravity (in Mega-years). (Default: FLT_MAX)
+     timestep_age_threshold_Myr:           25.       # (Optional) Age at which sink switch from young to old for time-stepping purposes (in Mega-years). (Default: FLT_MAX)
+     max_timestep_young_Myr:               1.0       # (Optional) Maximal time-step length of young sinks (in Mega-years). (Default: FLT_MAX)
+     max_timestep_old_Myr:                 5.0       # (Optional) Maximal time-step length of old sinks (in Mega-years). (Default: FLT_MAX)
+     n_IMF:                                 2        # (Optional) Number of times the IMF mass can be swallowed in a single timestep. (Default: FLTM_MAX)
+     tolerance_SF_timestep:                 0.5      # (Optional) Tolerance parameter for SF timestep constraint. (Default: 0.5)
+
+.. warning::
+   Some parameter choices can greatly impact the outcome of your simulations. Think twice when choosing them.
+
+Sink accretion radius
+~~~~~~~~~~~~~~~~~~~~~
+
+The most critical parameter is ``cut_off_radius``. As explained in the theory, to form a sink, the gas smoothing kernel edge :math:`\gamma_k h` (:math:`\gamma_k` is a kernel dependent constant) must be smaller than ``cut_off_radius`` (if this criterion is enabled). Therefore, the cut-off radius strongly depends on the resolution of your simulations. Moreover, if you use a minimal gas smoothing length `h`, and plan to use sink particles, consider whether the cut-off radius will meet the smoothing length criterion. If `h` never meets the aforementioned criterion, you will never form sinks and thus never have stars.
+
+On the contrary, if you set a too high cut-off radius, then sinks will accrete a lot of gas particles and spawn a lot of stars in the same cell, which the code might not like and crash with the error:
+
+``runner_others.c:runner_do_star_formation_sink():274: Too many stars in the cell tree leaf! The sorting task will not be able to perform its duties. Possible solutions: (1) The code need to be run with different star formation parameters to reduce the number of star particles created. OR (2) The size of the sorting stack must be increased in runner_sort.c.``
+
+This problem can be mitigated by choosing a higher value of ``stellar_particle_mass_Msun`` and ``stellar_particle_mass_first_stars_Msun``, or higher values of ``minimal_discrete_mass_Msun`` and ``minimal_discrete_mass_first_stars_Msun``. Of course, this comes at the price of having fewer individual stars. Finally, all parameters will depend on your needs.
+
+*If you do not want to change your parameters*, you can increase the ``sort_stack_size`` variable at the beginning ``runner_sort.c``. The default value is 10 in powers of 2 (so the stack size is 1024 particles). Increase it to the desired value. Be careful to not overestimate this.
+
+Note that if a cutoff radius is not specified, and the radius is instead left to vary with the local gas density, the smoothing length criterion is always satisfied.
+
+Guide to choose the the accretion radius or the density threshold
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+We provide some advice to help you set up the sink accretion radius or the threshold density appropriately.
+
+First, you must choose either the sink accretion radius or the threshold density. Choosing the density might be easier based on your previous work or if you have an expected star formation density. Once you fix the density or the accretion radius, you can use the following formula to estimate the remaining parameter. In the code, the gas smoothing length is determined with:
+
+.. math::
+   h = \eta \left( \frac{X_{\text{H}} m_B}{m_{\text{H}} n_{\text{H}}} \right)^{1/3} \, ,
+
+where :math:`\eta` is a constant related to the number of neighbours in the kernel, :math:`X_{\text{H}}` is the hydrogen mass fraction, :math:`m_B` the gas particle's mass, :math:`m_{\text{H}}` the hydrogen particle mass and :math:`n_{\text{H}}` the hydrogen number density.
+
+Let us provide an example. In GEAR, we do not model physical processes below the parsec scale. Hence, let us take :math:`h \sim 1` pc. In zoom-in simulations we have :math:`m_B \simeq 95 \; M_{\odot}`. The remaining parameters are :math:`\eta = 1.2348` and :math:`X_{\text{H}} = 0.76`. So, after inverting the formula, we find :math:`n_H \simeq 5500 \text{ hydrogen atoms/cm}^3`. In practice, we use :math:`n_H = 1000 \text{ hydrogen atoms/cm}^3`, close to the estimation, and an accretion radius :math:`r_{\text{acc}} = 10` pc. These values are slightly different for safety reasons, but they are consistent.
+
+Remember that this was a way, among others, to determine good accretion radius and threshold density. It can help you with your first runs with sink particles.
+
+Comment on star formation efficiency
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Notice that this model does not have parameters to control the star formation rate of the sink. The SFR is self-regulated by the gas/sink accretion and other feedback mechanisms. Supernovae tend to create bubbles of lower density at the site of star formation, removing the gas and preventing further gas accretion. However, the sink might run into this stack size problem by the time the first supernovae explode. Other pre-stellar feedback mechanisms could do the job earlier, though they are not implemented in GEAR.
+
+.. note:: 
+   We provide a piece of general advice: do some calibration on low-resolution simulations. This will help to see what works and what does not work. Keep in mind that you might want to put a higher ``stellar_particle_mass_X_Msun`` at the beginning to avoid spawning too many stars. For the high-resolution simulations, you then can lower the particle's mass.

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/Task/adding_your_own.rst

@@ -217,7 +217,8 @@ and the second kick cannot be done before the cooling::
 
 
 The next step is to activate your task
-in ``engine_marktasks_mapper`` in ``engine_marktasks.c``::
+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);

doc/RTD/source/Task/adding_your_own_neighbour_loop.rst

@@ -296,8 +296,9 @@ call::
 
 
 
-The next step is to activate your task
-in ``engine_marktasks_mapper`` in ``engine_marktasks.c``::
+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? */

doc/RTD/source/conf.py

@@ -23,9 +23,9 @@ copyright = "2014-2023, SWIFT Collaboration"
 author = "SWIFT Team"
 
 # The short X.Y version
-version = "1.0"
+version = "2025.04"
 # The full version, including alpha/beta/rc tags
-release = "1.0.0"
+release = "2025.04"
 
 # -- Find additional scripts to run as part of the documentation build -------
 import glob
@@ -147,7 +147,7 @@ latex_documents = [
     (
         master_doc,
         "SWIFT-user-manual.tex",
-        "SWIFT user \& developer documentation",
+        "SWIFT user & developer documentation",
         "SWIFT Collaboration",
         "manual",
     )

doc/onboardingGuide/source/dependencies.rst

@@ -8,9 +8,9 @@ To compile SWIFT, you will need the following libraries:
 HDF5
 ~~~~
 
-Version 1.8.x or higher is required. Input and output files are stored as HDF5
+Version 1.10.x or higher is required. Input and output files are stored as HDF5
 and are compatible with the GADGET-2 specification. A parallel-HDF5 build
-and HDF5 > 1.10.x is strongly recommended when running over MPI.
+and HDF5 >= 1.12.x is recommended when running over MPI.
 
 MPI
 ~~~
@@ -33,7 +33,15 @@ GSL
 ~~~
 The GSL 2.x is required for cosmological integration.
 
+In most cases the configuration script will be able to detect the libraries
+installed on the system. If that is not the case, the script can be pointed
+towards the libraries' location using the following parameters
 
+.. code-block:: bash
+
+  ./configure --with-gsl=<PATH-TO-GSL>
+
+and similar for the other libaries.
 
 Optional Dependencies
 =====================
@@ -48,10 +56,6 @@ TCmalloc/Jemalloc/TBBmalloc
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 TCmalloc/Jemalloc/TBBmalloc are used for faster memory allocations when available.
 
-DOXYGEN
-~~~~~~~
-You can build documentation for SWIFT with DOXYGEN.
-
 Python
 ~~~~~~
 To run the examples, you will need python 3 and some of the standard scientific libraries (numpy, matplotlib).

doc/onboardingGuide/source/initial_setup.rst

@@ -13,13 +13,13 @@ We use autotools for setup. To get a basic running version of the code (the exec
 MacOS Specific Oddities
 ~~~~~~~~~~~~~~~~~~~~~~~
 
-To build on MacOS you will need to disable compiler warnings due to an
+To build on MacOS you will need to enable compiler warnings due to an
 incomplete implementation of pthread barriers. DOXYGEN also has some issues on
 MacOS, so it is best to leave it out. To configure:
 
 .. code-block:: bash
 
-  ./configure --disable-compiler-warnings \ 
+  ./configure --enable-compiler-warnings \
       --disable-doxygen-doc
 
 When using the ``clang`` compiler, the hand-written vectorized routines

examples/AGORA/AgoraDisk/agora_disk.yml

@@ -79,15 +79,16 @@ GrackleCooling:
   max_steps: 1000
   convergence_limit: 1e-2
   thermal_time_myr: 0
+  maximal_density_Hpcm3: -1   # Maximal density (in hydrogen atoms/cm^3) for cooling. Higher densities are floored to this value to ensure grackle works properly when interpolating beyond the cloudy_table maximal density. A value < 0 deactivates this parameter.
 
 
 GEARStarFormation:
   star_formation_mode: agora        # default or agora
   star_formation_efficiency: 0.01   # star formation efficiency (c_*)
-  maximal_temperature:  1e10        # Upper limit to the temperature of a star forming particle
+  maximal_temperature_K:     1e10   # Upper limit to the temperature of a star forming particle
+  density_threshold_Hpcm3:   10     # Density threshold in Hydrogen atoms/cm3
   n_stars_per_particle: 1
   min_mass_frac: 0.5
-  density_threshold:   1.67e-23     # Density threashold in g/cm3
 
 
 GEARPressureFloor:

examples/AGORA/AgoraDisk/getGrackleCoolingTable.sh

 #!/bin/bash
 
-wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/CoolingTables/CloudyData_UVB=HM2012.h5
+wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/CoolingTables/CloudyData_UVB=HM2012.h5

examples/AGORA/AgoraDisk/getIC.sh

@@ -6,6 +6,5 @@ if [ "$#" -ne 1 ]; then
     exit
 fi
 
-wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/AgoraDisk/snap.$1.hdf5
+wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/AgoraDisk/snap.$1.hdf5
 mv snap.$1.hdf5 agora_disk.hdf5
-

examples/ChemistryTests/MetalAdvectionTestEAGLE/getGlass.sh

 #!/bin/bash
-wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassPlane_64.hdf5
+wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassPlane_64.hdf5