diff --git a/doc/RTD/source/SubgridModels/GEAR/chemistry.rst b/doc/RTD/source/SubgridModels/GEAR/chemistry.rst
new file mode 100644
index 0000000000000000000000000000000000000000..aefa7a3608b9253cb2dae38ac77cf35fdad43bb9
--- /dev/null
+++ b/doc/RTD/source/SubgridModels/GEAR/chemistry.rst
@@ -0,0 +1,27 @@
+.. GEAR sub-grid model chemistry
+ Darwin Roduit, 30th March 2025
+
+.. gear_chemistry:
+
+.. _gear_chemistry:
+
+Chemistry
+=========
+
+Feedback mechanisms such as supernova feedback transfer metals to the nearby gas particles. There is no mass exchange (mass advection) in SPH methods, as in grid-based or moving-mesh hydrodynamics solvers. As a result, the metals received by the gas particles are locked in these particles. Grid or moving-mesh codes have mass advection and often implement passive scalar advection of metals to model metal mixing. GEAR implements different methods to model metal mixing.
+
+.. _gear_smoothed_metallicity:
+
+Smoothed metallicity
+--------------------
+
+The smoothed metallicity scheme consists in using the SPH to smooth the metallicity of each particle over the neighbors. It is worth to point the fact that we are *not exchanging* any metals but only smoothing it. The parameter ``GEARChemistry:initial_metallicity`` set the (non smoothed) initial mass fraction of each element for all the particles and ``GEARChemistry:scale_initial_metallicity`` use the feedback table to scale the initial metallicity of each element according the Sun's composition. If ``GEARChemistry:initial_metallicity`` is negative, then the metallicities are read from the initial conditions.
+
+For this chemistry scheme the parameters are:
+
+.. code:: YAML
+
+ GEARChemistry:
+ initial_metallicity: 1 # Initial metallicity of the gas (mass fraction)
+ scale_initial_metallicity: 1 # Should we scale the initial metallicity with the solar one?
+
diff --git a/doc/RTD/source/SubgridModels/GEAR/feedback.rst b/doc/RTD/source/SubgridModels/GEAR/feedback.rst
new file mode 100644
index 0000000000000000000000000000000000000000..72a8996d29e51575e367609868e5ced72ac2a48d
--- /dev/null
+++ b/doc/RTD/source/SubgridModels/GEAR/feedback.rst
@@ -0,0 +1,144 @@
+.. GEAR sub-grid model feedback and stellar evolution
+ Loic Hausammann, 17th April 2020
+ Darwin Roduit, 30th March 2025
+
+.. _gear_stellar_evolution_and_feedback:
+
+Stellar evolution and feedback
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The feedback is composed of a few different models:
+ - The initial mass function (IMF) defines the quantity of each type of stars,
+ - The lifetime of a star defines when a star will explode (or simply die),
+ - The supernovae of type II (SNII) defines the rates and yields,
+ - The supernovae of type Ia (SNIa) defines the rates and yields,
+ - The energy injection that defines how to inject the energy / metals into the particles.
+
+Most of the parameters are defined inside a table (``GEARFeedback:yields_table``) but can be override with some parameters in the YAML file.
+I will not describe theses parameters more than providing them at the end of this section.
+Two different models exist for the supernovae (``GEARFeedback:discrete_yields``).
+In the continuous mode, we integrate the quantities over the IMF and then explodes a floating point number of stars (can be below 1 in some cases).
+In the discrete mode, we avoid the problem of floating points by rounding the number of supernovae (using a floor and randomly adding a supernovae depending on the fractional part) and then compute the properties for a single star at a time.
+
+Initial mass function
+^^^^^^^^^^^^^^^^^^^^^
+
+GEAR is using the IMF model from `Kroupa (2001) <https://ui.adsabs.harvard.edu/abs/2001MNRAS.322..231K/abstract>`_.
+We have a difference of 1 in the exponent due to the usage of IMF in mass and not in number.
+We also restrict the mass of the stars to be inside :math:`[0.05, 50] M_\odot`.
+Here is the default model used, but it can be easily adapted through the initial mass function parameters:
+
+.. math::
+ \xi(m) \propto m^{-\alpha_i}\, \textrm{where}\,
+ \begin{cases}
+ \alpha_0 = 0.3,\, & 0.01 \leq m / M_\odot < 0.08, \\
+ \alpha_1 = 1.3,\, & 0.08 \leq m / M_\odot < 0.50, \\
+ \alpha_2 = 2.3,\, & 0.50 \leq m / M_\odot < 1.00, \\
+ \alpha_3 = 2.3,\, & 1.00 \leq m / M_\odot,
+ \end{cases}
+
+
+Lifetime
+^^^^^^^^
+
+The lifetime of a star in GEAR depends only on two parameters: first its mass and then its metallicity.
+
+.. math::
+ \log(\tau(m)) = a(Z) \log^2(m) + b(Z) \log(m) + c(Z) \\ \\
+ a(Z) = -40.110 Z^2 + 5.509 Z + 0.7824 \\
+ b(Z) = 141.929 Z^2 - 15.889 Z - 3.2557 \\
+ c(Z) = -261.365 Z^2 + 17.073 Z + 9.8661
+
+where :math:`\tau` is the lifetime in years, :math:`m` is the mass of the star (in solar mass) and Z the metallicity of the star.
+The parameters previously given are the default ones, they can be modified in the parameters file.
+
+Supernovae II
+^^^^^^^^^^^^^
+
+The supernovae rate is simply given by the number of stars massive enough that end their life at the required time.
+
+.. math::
+ \dot{N}_\textrm{SNII}(t) = \int_{M_l}^{M_u} \delta(t - \tau(m)) \frac{\phi(m)}{m} \mathrm{d}m
+
+where :math:`M_l` and :math:`M_u` are the lower and upper mass limits for a star exploding in SNII, :math:`\delta` is the Dirac function and :math:`\phi` is the initial mass function (in mass).
+
+The yields for SNII cannot be written in an analytical form, they depend on a few different tables that are based on the work of `Kobayashi et al. (2000) <https://ui.adsabs.harvard.edu/abs/2000ApJ...539...26K/abstract>`_ and `Tsujimoto et al. (1995) <https://ui.adsabs.harvard.edu/abs/1995MNRAS.277..945T/abstract>`_.
+
+Supernovae Ia
+^^^^^^^^^^^^^
+
+The supernovae Ia are a bit more complicated as they involve two different stars.
+
+.. math::
+ \dot{N}_\textrm{SNIa}(t) = \left( \int_{M_{p,l}}^{M_{p,u}} \frac{\phi(m)}{m} \mathrm{d}m \right) \sum_i b_i \int_{M_{d,l,i}}^{M_{d,u,i}}
+ \delta(t-\tau(m)) \frac{\phi_d(m)}{m}\mathrm{d}m
+
+.. math::
+ \phi_d(m) \propto m^{-0.35}
+
+where :math:`M_{p,l}` and :math:`M_{p,u}` are the mass limits for a progenitor of a white dwarf, :math:`b_i` is the probability to have a companion and
+:math:`M_{d,l,i}` and :math:`M_{d,u,i}` are the mass limits for each type of companion.
+The first parenthesis represents the number of white dwarfs and the second one the probability to form a binary.
+
++------------------+--------------------+-------------------+------------------+
+| Companion | :math:`M_{d,l,i}` | :math:`M_{d,u,i}` | :math:`b_i` |
++==================+====================+===================+==================+
+| Red giant | 0.9 | 1.5 | 0.02 |
++------------------+--------------------+-------------------+------------------+
+| Main sequence | 1.8 | 2.5 | 0.05 |
++------------------+--------------------+-------------------+------------------+
+
+The yields are based on the same papers than the SNII.
+
+Supernova energy injection
+^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+When a star goes into a supernova (type II and Ia), the stellar evolution determines how much energy (``GEARFeedback:supernovae_energy_erg`` TO BE CHECKED), mass and metals are released during the explosion. The energy can be ditributed as internal/thermal energy or as momentum. Thus, we need to distribute internal energy, momentum, mass and metals to the 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. We describe the two schemes in the :ref:`gear_sn_feedback_models` page.
+
+
+Generating a new table
+^^^^^^^^^^^^^^^^^^^^^^
+
+The feedback table is an HDF5 file with the following structure:
+
+.. graphviz:: feedback_table.dot
+
+where the solid (dashed) squares represent a group (a dataset) with the name of the object underlined and the attributes written below. Everything is in solar mass or without units (e.g. mass fraction or unitless constant).
+In ``Data``, the attribute ``elts`` is an array of string with the element names (the last should be ``Metals``, it corresponds to the sum of all the elements), ``MeanWDMass`` is the mass of the white dwarfs
+and ``SolarMassAbundances`` is an array of float containing the mass fraction of the different element in the sun.
+In ``IMF``, ``n + 1`` is the number of part in the IMF, ``as`` are the exponent (``n+1`` elements), ``ms`` are the mass limits between each part (``n`` elements) and
+``Mmin`` (``Mmax``) is the minimal (maximal) mass of a star.
+In ``LifeTimes``, the coefficient are given in the form of a single table (``coeff_z`` with a 3x3 shape).
+In ``SNIa``, ``a`` is the exponent of the distribution of binaries, ``bb1`` and ``bb2`` are the coefficient :math:`b_i` and the other attributes follow the same names than in the SNIa formulas.
+The ``Metals`` group from the ``SNIa`` contains the name of each elements (``elts``) and the metal mass fraction ejected by each supernovae (``data``) in the same order. They must contain the same elements than in ``Data``.
+Finally for the ``SNII``, the mass limits are given by ``Mmin`` and ``Mmax``. For the yields, the datasets required are ``Ej`` (mass fraction ejected [processed]), ``Ejnp`` (mass fraction ejected [non processed]) and one dataset for each element present in ``elts``. The datasets should all have the same size, be uniformly sampled in log and contains the attributes ``min`` (mass in log for the first element) and ``step`` (difference of mass in log between two elements).
+
+.. code:: YAML
+
+ GEARFeedback:
+ supernovae_energy_erg: 0.1e51 # Energy released by a single supernovae.
+ yields_table: chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5 # Table containing the yields.
+ discrete_yields: 0 # Should we use discrete yields or the IMF integrated one?
+ GEARInitialMassFunction:
+ number_function_part: 4 # Number of different part in the IMF
+ exponents: [0.7, -0.8, -1.7, -1.3] # Exponents of each part of the IMF
+ mass_limits_msun: [0.05, 0.08, 0.5, 1, 50] # Limits in mass between each part of the IMF
+ GEARLifetime:
+ quadratic: [-40.1107, 5.50992, 0.782432] # Quadratic terms in the fit
+ linear: [141.93, -15.8895, -3.25578] # Linear terms in the fit
+ constant: [-261.366, 17.0735, 9.86606] # Constant terms in the fit
+ GEARSupernovaeIa:
+ exponent: -0.35 # Exponent for the distribution of companions
+ min_mass_white_dwarf_progenitor: 3 # Minimal mass of a progenitor of white dwarf
+ max_mass_white_dwarf_progenitor: 8 # Maximal mass of a progenitor of white dwarf
+ max_mass_red_giant: 1.5 # Maximal mass for a red giant
+ min_mass_red_giant: 0.9 # Minimal mass for a red giant
+ coef_red_giant: 0.02 # Coefficient for the distribution of red giants companions
+ max_mass_main_sequence: 2.6 # Maximal mass for a main sequence star
+ min_mass_main_sequence: 1.8 # Minimal mass for a main sequence star
+ coef_main_sequence: 0.05 # Coefficient for the distribution of main sequence companions
+ white_dwarf_mass: 1.38 # Mass of a white dwarf
+ GEARSupernovaeII:
+ interpolation_size: 200 # Number of elements for the interpolation of the data
diff --git a/doc/RTD/source/SubgridModels/GEAR/gear_model.rst b/doc/RTD/source/SubgridModels/GEAR/gear_model.rst
new file mode 100644
index 0000000000000000000000000000000000000000..8714763a895ce78d32a938d83e0ea62b546d4ed5
--- /dev/null
+++ b/doc/RTD/source/SubgridModels/GEAR/gear_model.rst
@@ -0,0 +1,204 @@
+.. GEAR sub-grid model
+ Loic Hausammann, 17th April 2020
+ Darwin Roduit, 30th March 2025
+
+.. _gear_pressure_floor:
+
+Pressure Floor
+~~~~~~~~~~~~~~
+
+In order to avoid the artificial collapse of unresolved clumps, a minimum in pressure is applied to the particles.
+This additional pressure can be seen as the pressure due to unresolved hydrodynamics turbulence and is given by:
+
+.. math::
+ P_\textrm{Jeans} = \frac{\rho}{\gamma} \frac{4}{\pi} G h^2 \rho N_\textrm{Jeans}^{2/3}
+
+where :math:`\rho` is the density, :math:`\gamma` the adiabatic index, :math:`G` is the gravitational constant,
+:math:`h` the kernel support and :math:`N_\textrm{Jeans}` (``GEARPressureFloor:jeans_factor`` in the parameter file) is the number of particle required in order to resolve a clump.
+
+
+This must be directly implemented into the hydro schemes, therefore only a subset of schemes (Gadget-2, SPHENIX and Pressure-Energy) have the floor available.
+In order to implement it, you need equation 12 in `Hopkins 2013 <https://arxiv.org/abs/1206.5006>`_:
+
+.. math::
+ m_i \frac{\mathrm{d}v_i}{\mathrm{d}t} = - \sum_j x_i x_j \left[ \frac{P_i}{y_i^2} f_{ij} \nabla_i W_{ij}(h_i) + \frac{P_j}{y_j^2} f_{ji} \nabla_j W_{ji}(h_j) \right]
+
+and simply replace the :math:`P_i, P_j` by the pressure with the floor (when the pressure is below the floor).
+Here the :math:`x, y` are simple weights that should never have the pressure floor included even if they are related to the pressure (e.g. pressure-entropy).
+
+.. code:: YAML
+
+ GEARPressureFloor:
+ jeans_factor: 10. # Number of particles required to suppose a resolved clump and avoid the pressure floor.
+
+
+.. _gear_star_formation:
+
+
+Star formation
+~~~~~~~~~~~~~~
+
+The star formation is done in two steps: first we check if a particle is in the star forming regime and then we use a stochastic approach to transform the gas particles into stars.
+
+A particle is in the star forming regime if:
+ - The velocity divergence is negative (:math:`\nabla\cdot v < 0`),
+ - The temperature is lower than a threshold (:math:`T < T_t` where :math:`T_t` is defined with ``GEARStarFormation:maximal_temperature_K``),
+ - The gas density is higher than a threshold (:math:`\rho > \rho_t` where :math:`\rho_t` is defined with ``GEARStarFormation:density_threshold_Hpcm3``)
+ - The particle reaches the pressure floor (:math:`\rho > \frac{\pi}{4 G N_\textrm{Jeans}^{2/3} h^2}\frac{\gamma k_B T}{\mu m_p}` where :math:`N_\textrm{Jeans}` is defined in the pressure floor).
+
+If ``GEARStarFormation:star_formation_mode`` is set to ``agora``, the condition on the pressure floor is ignored. Its default value is ``default``.
+
+A star will be able to form if a randomly drawn number is below :math:`\frac{m_g}{m_\star}\left(1 - \exp\left(-c_\star \Delta t / t_\textrm{ff}\right)\right)` where :math:`t_\textrm{ff}` is the free fall time, :math:`\Delta t` is the time step of the particle and :math:`c_\star` is the star formation coefficient (``GEARStarFormation:star_formation_efficiency``), :math:`m_g` the mass of the gas particle and :math:`m_\star` the mass of the possible future star. The mass of the star is computed from the average gas mass in the initial conditions divided by the number of possible stars formed per gas particle (``GEARStarFormation:n_stars_per_particle``). When we cannot have enough mass to form a second star (defined with the fraction of mass ``GEARStarFormation:min_mass_frac``), we fully convert the gas particle into a stellar particle. Once the star is formed, we move it a bit in a random direction and fraction of the smoothing length in order to avoid any division by 0.
+
+Currently, only the following hydro schemes are compatible: SPHENIX, Gadget2, minimal SPH, Gasoline-2 and Pressure-Energy.
+Implementing the other hydro schemes is not complicated but requires some careful thinking about the cosmological terms in the definition of the velocity divergence (comoving vs non comoving coordinates and if the Hubble flow is included or not).
+
+.. code:: YAML
+
+ GEARStarFormation:
+ star_formation_efficiency: 0.01 # star formation efficiency (c_*)
+ maximal_temperature_K: 3e4 # Upper limit to the temperature of a star forming particle
+ density_threshold_Hpcm3: 10 # Density threshold (Hydrogen atoms/cm^3) for star formation
+ n_stars_per_particle: 4 # Number of stars that an hydro particle can generate
+ min_mass_frac: 0.5 # Minimal mass for a stellar particle as a fraction of the average mass for the stellar particles.
+
+Initial Conditions
+++++++++++++++++++
+
+Note that if in the initial conditions, the time of formation of a stellar particle is given (``BirthTime``)
+and set to a negative value, the stellar particle will provide no feedback.
+A similar behavior will be obtained if the parameter ``Stars:overwrite_birth_time`` is set to 1 and
+``Stars:birth_time`` to -1.
+
+
+.. _gear_grackle_cooling:
+
+Cooling: Grackle
+~~~~~~~~~~~~~~~~
+
+Grackle is a chemistry and cooling library presented in `B. Smith et al. 2017 <https://ui.adsabs.harvard.edu/abs/2017MNRAS.466.2217S>`_
+(do not forget to cite if used). Four different modes are available:
+equilibrium, 6 species network (H, H\\( ^+ \\), e\\( ^- \\), He, He\\( ^+ \\)
+and He\\( ^{++} \\)), 9 species network (adds H\\(^-\\), H\\(_2\\) and
+H\\(_2^+\\)) and 12 species (adds D, D\\(^+\\) and HD). Following the same
+order, the swift cooling options are ``grackle_0``, ``grackle_1``, ``grackle_2``
+and ``grackle_3`` (the numbers correspond to the value of
+``primordial_chemistry`` in Grackle). It also includes some metal cooling (on/off with ``GrackleCooling:with_metal_cooling``), self-shielding
+methods and UV background (on/off with ``GrackleCooling:with_UV_background``). In order to use the Grackle cooling, you will need
+to provide a HDF5 table computed by Cloudy (``GrackleCooling:cloudy_table``).
+
+Configuring and compiling SWIFT with Grackle
+++++++++++++++++++++++++++++++++++++++++++++
+
+In order to compile SWIFT with Grackle, you need to provide the options ``with-chemistry=GEAR`` and ``with-grackle=$GRACKLE_ROOT``
+where ``$GRACKLE_ROOT`` is the root of the install directory (not the ``lib``).
+
+.. warning::
+ (State 2023) Grackle is experiencing current development, and the API is subject
+ to changes in the future. For convenience, a frozen version is hosted as a fork
+ on github here: https://github.com/mladenivkovic/grackle-swift .
+ The version available there will be tried and tested and ensured to work with
+ SWIFT.
+
+ Additionally, that repository hosts files necessary to install that specific
+ version of grackle with spack.
+
+To compile it, run
+the following commands from the root directory of Grackle:
+``./configure; cd src/clib``.
+Update the variables ``LOCAL_HDF5_INSTALL`` and ``MACH_INSTALL_PREFIX`` in
+the file ``src/clib/Make.mach.linux-gnu``.
+Finish with ``make machine-linux-gnu; make && make install``.
+Note that we require the 64 bit float version of Grackle, which should be the default setting.
+(The precision can be set while compiling grackle with ``make precision-64``).
+If you encounter any problem, you can look at the `Grackle documentation <https://grackle.readthedocs.io/en/latest/>`_
+
+You can now provide the path given for ``MACH_INSTALL_PREFIX`` to ``with-grackle``.
+
+Parameters
+++++++++++
+
+When starting a simulation without providing the different element fractions in the non equilibrium mode, the code supposes an equilibrium and computes them automatically.
+The code uses an iterative method in order to find the correct initial composition and this method can be tuned with two parameters. ``GrackleCooling:max_steps`` defines the maximal number of steps to reach the convergence and ``GrackleCooling:convergence_limit`` defines the tolerance in the relative error.
+
+In the parameters file, a few different parameters are available.
+
+- ``GrackleCooling:redshift`` defines the redshift to use for the UV background (for cosmological simulation, it must be set to -1 in order to use the simulation's redshift).
+
+- ``GrackleCooling:provide_*_heating_rates`` can enable the computation of user provided heating rates (such as with the radiative transfer) in either volumetric or specific units.
+
+- Feedback can be made more efficient by turning off the cooling during a few Myr (``GrackleCooling:thermal_time_myr``) for the particles touched by a supernovae.
+
+- The self shielding method is defined by ``GrackleCooling:self_shielding_method`` where 0 means no self shielding, > 0 means a method defined in Grackle (see Grackle documentation for more information) and -1 means GEAR's self shielding that simply turn off the UV background when reaching a given density (``GrackleCooling:self_shielding_threshold_atom_per_cm3``).
+
+- The initial elemental abundances can be specified with ``initial_nX_to_nY_ratio``, e.g. if you want to specify an initial HII to H abundance. A negative value ignores the parameter. The complete list can be found below.
+
+A maximal (physical) density must be set with the ``GrackleCooling:maximal_density_Hpcm3 parameter``. The density passed to Grackle is *the minimum of this density and the gas particle (physical) density*. A negative value (:math:`< 0`) deactivates the maximal density, i.e. there is no maximal density limit.
+The purpose of this parameter is the following. The Cloudy tables provided by Grackle are limited in density (typically to :math:`10^4 \; \mathrm{hydrogen \; atoms/cm}^3`). In high-resolution simulations, particles can have densities higher than :math:`10^4 \; \mathrm{hydrogen \; atoms/cm}^3`. This maximal density ensures that we pass a density within the interpolation ranges of the table, should the density exceed it.
+It can be a solution to some of the following errors (with a translation of what the values mean):
+
+.. code:: text
+
+ inside if statement solve rate cool: 0 0
+ MULTI_COOL iter > 10000 at j,k = 1 1
+ FATAL error (2) in MULTI_COOL
+ dt = 1.092E-08 ttmin = 7.493E-12
+ 2.8E-19 // sub-cycling timestep
+ 7.5E-12 // time elapsed (in the sub-cycle)
+ 2.2E+25 // derivative of the internal energy
+ T
+
+.. note::
+ This problem is particularly relevant with metal cooling enabled. Another solution is to modify the tables. But one is not exempted from exceeding the table maximal density value, since Grackle does not check if the particle density is larger than the table maximal density.
+
+
+Complete parameter list
+-----------------------
+
+Here is the complete section in the parameter file:
+
+.. code:: YAML
+
+ GrackleCooling:
+ cloudy_table: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
+ with_UV_background: 1 # Enable or not the UV background
+ redshift: 0 # Redshift to use (-1 means time based redshift)
+ with_metal_cooling: 1 # Enable or not the metal cooling
+ provide_volumetric_heating_rates: 0 # (optional) User provide volumetric heating rates
+ provide_specific_heating_rates: 0 # (optional) User provide specific heating rates
+ max_steps: 10000 # (optional) Max number of step when computing the initial composition
+ convergence_limit: 1e-2 # (optional) Convergence threshold (relative) for initial composition
+ thermal_time_myr: 5 # (optional) Time (in Myr) for adiabatic cooling after a feedback event.
+ self_shielding_method: -1 # (optional) Grackle (1->3 for Grackle's ones, 0 for none and -1 for GEAR)
+ self_shielding_threshold_atom_per_cm3: 0.007 # Required only with GEAR's self shielding. Density threshold of the self shielding
+ maximal_density_Hpcm3: 1e4 # 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.
+ HydrogenFractionByMass : 1. # Hydrogen fraction by mass (default is 0.76)
+
+ use_radiative_transfer : 1 # Arrays of ionization and heating rates are provided
+ RT_heating_rate_cgs : 0 # heating rate in units of / nHI_cgs
+ RT_HI_ionization_rate_cgs : 0 # HI ionization rate in cgs [1/s]
+ RT_HeI_ionization_rate_cgs : 0 # HeI ionization rate in cgs [1/s]
+ RT_HeII_ionization_rate_cgs: 0 # HeII ionization rate in cgs [1/s]
+ RT_H2_dissociation_rate_cgs: 0 # H2 dissociation rate in cgs [1/s]
+
+ volumetric_heating_rates_cgs: 0 # Volumetric heating rate in cgs [erg/s/cm3]
+ specific_heating_rates_cgs: 0 # Specific heating rate in cgs [erg/s/g]
+ H2_three_body_rate : 1 # Specific the H2 formation three body rate (0->5,see Grackle documentation)
+ H2_cie_cooling : 0 # Enable/disable H2 collision-induced emission cooling from Ripamonti & Abel (2004)
+ H2_on_dust: 0 # Flag to enable H2 formation on dust grains
+ local_dust_to_gas_ratio : -1 # The ratio of total dust mass to gas mass in the local Universe (-1 to use the Grackle default value).
+ cmb_temperature_floor : 1 # Enable/disable an effective CMB temperature floor
+
+ initial_nHII_to_nH_ratio: -1 # initial nHII to nH ratio (number density ratio). Value is ignored if set to -1.
+ initial_nHeI_to_nH_ratio: -1 # initial nHeI to nH ratio (number density ratio). Value is ignored if set to -1.
+ initial_nHeII_to_nH_ratio: -1 # initial nHeII to nH ratio (number density ratio). Value is ignored if set to -1.
+ initial_nHeIII_to_nH_ratio: -1 # initial nHeIII to nH ratio (number density ratio). Value is ignored if set to -1.
+ initial_nDI_to_nH_ratio: -1 # initial nDI to nH ratio (number density ratio). Value is ignored if set to -1.
+ initial_nDII_to_nH_ratio: -1 # initial nDII to nH ratio (number density ratio). Value is ignored if set to -1.
+ initial_nHM_to_nH_ratio: -1 # initial nHM to nH ratio (number density ratio). Value is ignored if set to -1.
+ initial_nH2I_to_nH_ratio: -1 # initial nH2I to nH ratio (number density ratio). Value is ignored if set to -1.
+ initial_nH2II_to_nH_ratio: -1 # initial nH2II to nH ratio (number density ratio). Value is ignored if set to -1.
+ initial_nHDI_to_nH_ratio: -1 # initial nHDI to nH ratio (number density ratio). Value is ignored if set to -1.
+
+.. note::
+ A simple example running SWIFT with Grackle can be find in ``examples/Cooling/CoolingBox``. A more advanced example combining heating and cooling (with heating and ionization sources) is given in ``examples/Cooling/CoolingHeatingBox``. ``examples/Cooling/CoolingWithPrimordialElements/`` runs a uniform cosmological box with imposed abundances and let them evolve down to redshift 0.
diff --git a/doc/RTD/source/SubgridModels/GEAR/index.rst b/doc/RTD/source/SubgridModels/GEAR/index.rst
index 30b79689a28527f975a739925ae7835f90e398c4..26c8b2474e893221a882fcca41306f57cbf9adf7 100644
--- a/doc/RTD/source/SubgridModels/GEAR/index.rst
+++ b/doc/RTD/source/SubgridModels/GEAR/index.rst
@@ -6,341 +6,16 @@ GEAR model
===========
GEAR's model are mainly described in `Revaz \& Jablonka <https://ui.adsabs.harvard.edu/abs/2018A%26A...616A..96R/abstract>`_.
-This model can be selected with the configuration option ``--with-subgrid=GEAR`` and run with the option ``--gear``. A few examples exist and can be found in ``examples/GEAR``.
+This model can be selected with the configuration option ``--with-subgrid=GEAR`` and run with the option ``--gear``. A few examples exist and can be found in ``examples/GEAR``, ``examples/SinkParticles`` or ``IsolatedGalaxy/IsolatedGalaxy_multi_component``.
-.. _gear_pressure_floor:
-Pressure Floor
-~~~~~~~~~~~~~~
+.. toctree::
+ :caption: Table of Contents
-In order to avoid the artificial collapse of unresolved clumps, a minimum in pressure is applied to the particles.
-This additional pressure can be seen as the pressure due to unresolved hydrodynamics turbulence and is given by:
+ gear_model
+ chemistry
+ feedback
+ supernova_feedback
+ sinks/index
+ output
-.. math::
- P_\textrm{Jeans} = \frac{\rho}{\gamma} \frac{4}{\pi} G h^2 \rho N_\textrm{Jeans}^{2/3}
-
-where :math:`\rho` is the density, :math:`\gamma` the adiabatic index, :math:`G` is the gravitational constant,
-:math:`h` the kernel support and :math:`N_\textrm{Jeans}` (``GEARPressureFloor:jeans_factor`` in the parameter file) is the number of particle required in order to resolve a clump.
-
-
-This must be directly implemented into the hydro schemes, therefore only a subset of schemes (Gadget-2, SPHENIX and Pressure-Energy) have the floor available.
-In order to implement it, you need equation 12 in `Hopkins 2013 <https://arxiv.org/abs/1206.5006>`_:
-
-.. math::
- m_i \frac{\mathrm{d}v_i}{\mathrm{d}t} = - \sum_j x_i x_j \left[ \frac{P_i}{y_i^2} f_{ij} \nabla_i W_{ij}(h_i) + \frac{P_j}{y_j^2} f_{ji} \nabla_j W_{ji}(h_j) \right]
-
-and simply replace the :math:`P_i, P_j` by the pressure with the floor (when the pressure is below the floor).
-Here the :math:`x, y` are simple weights that should never have the pressure floor included even if they are related to the pressure (e.g. pressure-entropy).
-
-.. code:: YAML
-
- GEARPressureFloor:
- jeans_factor: 10. # Number of particles required to suppose a resolved clump and avoid the pressure floor.
-
-
-
-.. _gear_grackle_cooling:
-
-Cooling: Grackle
-~~~~~~~~~~~~~~~~
-
-Grackle is a chemistry and cooling library presented in `B. Smith et al. 2017 <https://ui.adsabs.harvard.edu/abs/2017MNRAS.466.2217S>`_
-(do not forget to cite if used). Four different modes are available:
-equilibrium, 6 species network (H, H\\( ^+ \\), e\\( ^- \\), He, He\\( ^+ \\)
-and He\\( ^{++} \\)), 9 species network (adds H\\(^-\\), H\\(_2\\) and
-H\\(_2^+\\)) and 12 species (adds D, D\\(^+\\) and HD). Following the same
-order, the swift cooling options are ``grackle_0``, ``grackle_1``, ``grackle_2``
-and ``grackle_3`` (the numbers correspond to the value of
-``primordial_chemistry`` in Grackle). It also includes some metal cooling (on/off with ``GrackleCooling:with_metal_cooling``), self-shielding
-methods and UV background (on/off with ``GrackleCooling:with_UV_background``). In order to use the Grackle cooling, you will need
-to provide a HDF5 table computed by Cloudy (``GrackleCooling:cloudy_table``).
-
-
-When starting a simulation without providing the different element fractions in the non equilibrium mode, the code supposes an equilibrium and computes them automatically.
-The code uses an iterative method in order to find the correct initial composition and this method can be tuned with two parameters. ``GrackleCooling:max_steps`` defines the maximal number of steps to reach the convergence and ``GrackleCooling:convergence_limit`` defines the tolerance in the relative error.
-
-In order to compile SWIFT with Grackle, you need to provide the options ``with-chemistry=GEAR`` and ``with-grackle=$GRACKLE_ROOT``
-where ``$GRACKLE_ROOT`` is the root of the install directory (not the ``lib``).
-
-.. warning::
- (State 2023) Grackle is experiencing current development, and the API is subject
- to changes in the future. For convenience, a frozen version is hosted as a fork
- on github here: https://github.com/mladenivkovic/grackle-swift .
- The version available there will be tried and tested and ensured to work with
- SWIFT.
-
- Additionally, that repository hosts files necessary to install that specific
- version of grackle with spack.
-
-
-To compile it, run
-the following commands from the root directory of Grackle:
-``./configure; cd src/clib``.
-Update the variables ``LOCAL_HDF5_INSTALL`` and ``MACH_INSTALL_PREFIX`` in
-the file ``src/clib/Make.mach.linux-gnu``.
-Finish with ``make machine-linux-gnu; make && make install``.
-Note that we require the 64 bit float version of Grackle, which should be the default setting.
-(The precision can be set while compiling grackle with ``make precision-64``).
-If you encounter any problem, you can look at the `Grackle documentation <https://grackle.readthedocs.io/en/latest/>`_
-
-You can now provide the path given for ``MACH_INSTALL_PREFIX`` to ``with-grackle``.
-
-In the parameters file, a few different parameters are available.
-``GrackleCooling:redshift`` defines the redshift to use for the UV background (for cosmological simulation, it must be set to -1 in order to use the simulation's redshift) and ``GrackleCooling:provide_*_heating_rates`` can enable the computation of user provided heating rates (such as with the radiative transfer) in either volumetric or specific units.
-
-For the feedback, it can be made more efficient by turning off the cooling during a few Myr (``GrackleCooling:thermal_time_myr``) for the particles touched by a supernovae.
-
-The self shielding method is defined by ``GrackleCooling:self_shielding_method`` where 0 means no self shielding, > 0 means a method defined in Grackle (see Grackle documentation for more information) and -1 means GEAR's self shielding that simply turn off the UV background when reaching a given density (``GrackleCooling:self_shielding_threshold_atom_per_cm3``).
-
-A maximal (physical) density must be set with the ``GrackleCooling:maximal_density_Hpcm3 parameter``. The density passed to Grackle is *the minimum of this density and the gas particle (physical) density*. A negative value (:math:`< 0`) deactivates the maximal density, i.e. there is no maximal density limit.
-The purpose of this parameter is the following. The Cloudy tables provided by Grackle are limited in density (typically to :math:`10^4 \; \mathrm{hydrogen \; atoms/cm}^3`). In high-resolution simulations, particles can have densities higher than :math:`10^4 \; \mathrm{hydrogen \; atoms/cm}^3`. This maximal density ensures that we pass a density within the interpolation ranges of the table, should the density exceed it.
-It can be a solution to some of the following errors (with a translation of what the values mean):
-
-.. code:: text
-
- inside if statement solve rate cool: 0 0
- MULTI_COOL iter > 10000 at j,k = 1 1
- FATAL error (2) in MULTI_COOL
- dt = 1.092E-08 ttmin = 7.493E-12
- 2.8E-19 // sub-cycling timestep
- 7.5E-12 // time elapsed (in the sub-cycle)
- 2.2E+25 // derivative of the internal energy
- T
-
-.. note::
- This problem is particularly relevant with metal cooling enabled. Another solution is to modify the tables. But one is not exempted from exceeding the table maximal density value, since Grackle does not check if the particle density is larger than the table maximal density.
-
-Here is the complete section in the parameter file:
-
-.. code:: YAML
-
- GrackleCooling:
- cloudy_table: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
- with_UV_background: 1 # Enable or not the UV background
- redshift: 0 # Redshift to use (-1 means time based redshift)
- with_metal_cooling: 1 # Enable or not the metal cooling
- provide_volumetric_heating_rates: 0 # (optional) User provide volumetric heating rates
- provide_specific_heating_rates: 0 # (optional) User provide specific heating rates
- max_steps: 10000 # (optional) Max number of step when computing the initial composition
- convergence_limit: 1e-2 # (optional) Convergence threshold (relative) for initial composition
- thermal_time_myr: 5 # (optional) Time (in Myr) for adiabatic cooling after a feedback event.
- self_shielding_method: -1 # (optional) Grackle (1->3 for Grackle's ones, 0 for none and -1 for GEAR)
- self_shielding_threshold_atom_per_cm3: 0.007 # Required only with GEAR's self shielding. Density threshold of the self shielding
- maximal_density_Hpcm3: 1e4 # Maximal density (in 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.
- HydrogenFractionByMass : 1. # Hydrogen fraction by mass (default is 0.76)
- use_radiative_transfer : 1 # Arrays of ionization and heating rates are provided
- RT_heating_rate_cgs : 0 # heating rate in units of / nHI_cgs
- RT_HI_ionization_rate_cgs : 0 # HI ionization rate in cgs [1/s]
- RT_HeI_ionization_rate_cgs : 0 # HeI ionization rate in cgs [1/s]
- RT_HeII_ionization_rate_cgs: 0 # HeII ionization rate in cgs [1/s]
- RT_H2_dissociation_rate_cgs: 0 # H2 dissociation rate in cgs [1/s]
- volumetric_heating_rates_cgs: 0 # Volumetric heating rate in cgs [erg/s/cm3]
- specific_heating_rates_cgs: 0 # Specific heating rate in cgs [erg/s/g]
- H2_three_body_rate : 1 # Specific the H2 formation three body rate (0->5,see Grackle documentation)
- H2_cie_cooling : 0 # Enable/disable H2 collision-induced emission cooling from Ripamonti & Abel (2004)
- cmb_temperature_floor : 1 # Enable/disable an effective CMB temperature floor
-
-.. note::
- A simple example running SWIFT with Grackle can be find in ``examples/Cooling/CoolingBox``. A more advanced example combining heating and cooling (with heating and ionization sources) is given in ``examples/Cooling/CoolingHeatingBox``.
-
-
-.. _gear_star_formation:
-
-Star formation
-~~~~~~~~~~~~~~
-
-The star formation is done in two steps: first we check if a particle is in the star forming regime and then we use a stochastic approach to transform the gas particles into stars.
-
-A particle is in the star forming regime if:
- - The velocity divergence is negative (:math:`\nabla\cdot v < 0`),
- - The temperature is lower than a threshold (:math:`T < T_t` where :math:`T_t` is defined with ``GEARStarFormation:maximal_temperature_K``),
- - The gas density is higher than a threshold (:math:`\rho > \rho_t` where :math:`\rho_t` is defined with ``GEARStarFormation:density_threshold_Hpcm3``)
- - The particle reaches the pressure floor (:math:`\rho > \frac{\pi}{4 G N_\textrm{Jeans}^{2/3} h^2}\frac{\gamma k_B T}{\mu m_p}` where :math:`N_\textrm{Jeans}` is defined in the pressure floor).
-
-If ``GEARStarFormation:star_formation_mode`` is set to ``agora``, the condition on the pressure floor is ignored. Its default value is ``default``.
-
-A star will be able to form if a randomly drawn number is below :math:`\frac{m_g}{m_\star}\left(1 - \exp\left(-c_\star \Delta t / t_\textrm{ff}\right)\right)` where :math:`t_\textrm{ff}` is the free fall time, :math:`\Delta t` is the time step of the particle and :math:`c_\star` is the star formation coefficient (``GEARStarFormation:star_formation_efficiency``), :math:`m_g` the mass of the gas particle and :math:`m_\star` the mass of the possible future star. The mass of the star is computed from the average gas mass in the initial conditions divided by the number of possible stars formed per gas particle (``GEARStarFormation:n_stars_per_particle``). When we cannot have enough mass to form a second star (defined with the fraction of mass ``GEARStarFormation:min_mass_frac``), we fully convert the gas particle into a stellar particle. Once the star is formed, we move it a bit in a random direction and fraction of the smoothing length in order to avoid any division by 0.
-
-Currently, only the following hydro schemes are compatible: SPHENIX and Gadget2.
-Implementing the other hydro schemes is not complicated but requires some careful thinking about the cosmological terms in the definition of the velocity divergence (comoving vs non comoving coordinates and if the Hubble flow is included or not).
-
-.. code:: YAML
-
- GEARStarFormation:
- star_formation_efficiency: 0.01 # star formation efficiency (c_*)
- maximal_temperature_K: 3e4 # Upper limit to the temperature of a star forming particle
- density_threshold_Hpcm3: 10 # Density threshold (Hydrogen atoms/cm^3) for star formation
- n_stars_per_particle: 4 # Number of stars that an hydro particle can generate
- min_mass_frac: 0.5 # Minimal mass for a stellar particle as a fraction of the average mass for the stellar particles.
-
-Sink particles
-~~~~~~~~~~~~~~
-
-GEAR now implements sink particles for star formation. Instead of stochastically transforming gas particles into stars as is done in the star formation scheme above when some criteria are met, 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.
-
-This was a brief overview of the model. More details can be found in :ref:`sink_GEAR_model` pages. Please refer to these for configuration, compilations and details of the model.
-
-
-Chemistry
-~~~~~~~~~
-
-In the chemistry, we are using the smoothed metallicity scheme that consists in using the SPH to smooth the metallicity of each particle over the neighbors. It is worth to point the fact that we are not exchanging any metals but only smoothing it. The parameter ``GEARChemistry:initial_metallicity`` set the (non smoothed) initial mass fraction of each element for all the particles and ``GEARChemistry:scale_initial_metallicity`` use the feedback table to scale the initial metallicity of each element according the Sun's composition. If ``GEARChemistry:initial_metallicity`` is negative, then the metallicities are read from the initial conditions.
-
-.. code:: YAML
-
- GEARChemistry:
- initial_metallicity: 1 # Initial metallicity of the gas (mass fraction)
- scale_initial_metallicity: 1 # Should we scale the initial metallicity with the solar one?
-
-
-.. _gear_feedback:
-
-Feedback
-~~~~~~~~
-
-The feedback is composed of a few different models:
- - The initial mass function (IMF) defines the quantity of each type of stars,
- - The lifetime of a star defines when a star will explode (or simply die),
- - The supernovae of type II (SNII) defines the rates and yields,
- - The supernovae of type Ia (SNIa) defines the rates and yields,
- - The energy injection that defines how to inject the energy / metals into the particles.
-
-Most of the parameters are defined inside a table (``GEARFeedback:yields_table``) but can be override with some parameters in the YAML file.
-I will not describe theses parameters more than providing them at the end of this section.
-Two different models exist for the supernovae (``GEARFeedback:discrete_yields``).
-In the continuous mode, we integrate the quantities over the IMF and then explodes a floating point number of stars (can be below 1 in some cases).
-In the discrete mode, we avoid the problem of floating points by rounding the number of supernovae (using a floor and randomly adding a supernovae depending on the fractional part) and then compute the properties for a single star at a time.
-
-Initial mass function
-^^^^^^^^^^^^^^^^^^^^^
-
-GEAR is using the IMF model from `Kroupa (2001) <https://ui.adsabs.harvard.edu/abs/2001MNRAS.322..231K/abstract>`_.
-We have a difference of 1 in the exponent due to the usage of IMF in mass and not in number.
-We also restrict the mass of the stars to be inside :math:`[0.05, 50] M_\odot`.
-Here is the default model used, but it can be easily adapted through the initial mass function parameters:
-
-.. math::
- \xi(m) \propto m^{-\alpha_i}\, \textrm{where}\,
- \begin{cases}
- \alpha_0 = 0.3,\, & 0.01 \leq m / M_\odot < 0.08, \\
- \alpha_1 = 1.3,\, & 0.08 \leq m / M_\odot < 0.50, \\
- \alpha_2 = 2.3,\, & 0.50 \leq m / M_\odot < 1.00, \\
- \alpha_3 = 2.3,\, & 1.00 \leq m / M_\odot,
- \end{cases}
-
-
-Lifetime
-^^^^^^^^
-
-The lifetime of a star in GEAR depends only on two parameters: first its mass and then its metallicity.
-
-.. math::
- \log(\tau(m)) = a(Z) \log^2(m) + b(Z) \log(m) + c(Z) \\ \\
- a(Z) = -40.110 Z^2 + 5.509 Z + 0.7824 \\
- b(Z) = 141.929 Z^2 - 15.889 Z - 3.2557 \\
- c(Z) = -261.365 Z^2 + 17.073 Z + 9.8661
-
-where :math:`\tau` is the lifetime in years, :math:`m` is the mass of the star (in solar mass) and Z the metallicity of the star.
-The parameters previously given are the default ones, they can be modified in the parameters file.
-
-Supernovae II
-^^^^^^^^^^^^^
-
-The supernovae rate is simply given by the number of stars massive enough that end their life at the required time.
-
-.. math::
- \dot{N}_\textrm{SNII}(t) = \int_{M_l}^{M_u} \delta(t - \tau(m)) \frac{\phi(m)}{m} \mathrm{d}m
-
-where :math:`M_l` and :math:`M_u` are the lower and upper mass limits for a star exploding in SNII, :math:`\delta` is the Dirac function and :math:`\phi` is the initial mass function (in mass).
-
-The yields for SNII cannot be written in an analytical form, they depend on a few different tables that are based on the work of `Kobayashi et al. (2000) <https://ui.adsabs.harvard.edu/abs/2000ApJ...539...26K/abstract>`_ and `Tsujimoto et al. (1995) <https://ui.adsabs.harvard.edu/abs/1995MNRAS.277..945T/abstract>`_.
-
-Supernovae Ia
-^^^^^^^^^^^^^
-
-The supernovae Ia are a bit more complicated as they involve two different stars.
-
-.. math::
- \dot{N}_\textrm{SNIa}(t) = \left( \int_{M_{p,l}}^{M_{p,u}} \frac{\phi(m)}{m} \mathrm{d}m \right) \sum_i b_i \int_{M_{d,l,i}}^{M_{d,u,i}}
- \delta(t-\tau(m)) \frac{\phi_d(m)}{m}\mathrm{d}m
-
-.. math::
- \phi_d(m) \propto m^{-0.35}
-
-where :math:`M_{p,l}` and :math:`M_{p,u}` are the mass limits for a progenitor of a white dwarf, :math:`b_i` is the probability to have a companion and
-:math:`M_{d,l,i}` and :math:`M_{d,u,i}` are the mass limits for each type of companion.
-The first parenthesis represents the number of white dwarfs and the second one the probability to form a binary.
-
-+------------------+--------------------+-------------------+------------------+
-| Companion | :math:`M_{d,l,i}` | :math:`M_{d,u,i}` | :math:`b_i` |
-+==================+====================+===================+==================+
-| Red giant | 0.9 | 1.5 | 0.02 |
-+------------------+--------------------+-------------------+------------------+
-| Main sequence | 1.8 | 2.5 | 0.05 |
-+------------------+--------------------+-------------------+------------------+
-
-The yields are based on the same papers than the SNII.
-
-Energy injection
-^^^^^^^^^^^^^^^^
-
-All the supernovae (type II and Ia) inject the same amount of energy into the surrounding gas (``GEARFeedback:supernovae_energy_erg``) and distribute it according to the hydro kernel.
-The same is done with the metals and the mass.
-
-
-Generating a new table
-^^^^^^^^^^^^^^^^^^^^^^
-
-The feedback table is an HDF5 file with the following structure:
-
-.. graphviz:: feedback_table.dot
-
-where the solid (dashed) squares represent a group (a dataset) with the name of the object underlined and the attributes written below. Everything is in solar mass or without units (e.g. mass fraction or unitless constant).
-In ``Data``, the attribute ``elts`` is an array of string with the element names (the last should be ``Metals``, it corresponds to the sum of all the elements), ``MeanWDMass`` is the mass of the white dwarfs
-and ``SolarMassAbundances`` is an array of float containing the mass fraction of the different element in the sun.
-In ``IMF``, ``n + 1`` is the number of part in the IMF, ``as`` are the exponent (``n+1`` elements), ``ms`` are the mass limits between each part (``n`` elements) and
-``Mmin`` (``Mmax``) is the minimal (maximal) mass of a star.
-In ``LifeTimes``, the coefficient are given in the form of a single table (``coeff_z`` with a 3x3 shape).
-In ``SNIa``, ``a`` is the exponent of the distribution of binaries, ``bb1`` and ``bb2`` are the coefficient :math:`b_i` and the other attributes follow the same names than in the SNIa formulas.
-The ``Metals`` group from the ``SNIa`` contains the name of each elements (``elts``) and the metal mass fraction ejected by each supernovae (``data``) in the same order. They must contain the same elements than in ``Data``.
-Finally for the ``SNII``, the mass limits are given by ``Mmin`` and ``Mmax``. For the yields, the datasets required are ``Ej`` (mass fraction ejected [processed]), ``Ejnp`` (mass fraction ejected [non processed]) and one dataset for each element present in ``elts``. The datasets should all have the same size, be uniformly sampled in log and contains the attributes ``min`` (mass in log for the first element) and ``step`` (difference of mass in log between two elements).
-
-.. code:: YAML
-
- GEARFeedback:
- supernovae_energy_erg: 0.1e51 # Energy released by a single supernovae.
- yields_table: chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5 # Table containing the yields.
- discrete_yields: 0 # Should we use discrete yields or the IMF integrated one?
- GEARInitialMassFunction:
- number_function_part: 4 # Number of different part in the IMF
- exponents: [0.7, -0.8, -1.7, -1.3] # Exponents of each part of the IMF
- mass_limits_msun: [0.05, 0.08, 0.5, 1, 50] # Limits in mass between each part of the IMF
- GEARLifetime:
- quadratic: [-40.1107, 5.50992, 0.782432] # Quadratic terms in the fit
- linear: [141.93, -15.8895, -3.25578] # Linear terms in the fit
- constant: [-261.366, 17.0735, 9.86606] # Constant terms in the fit
- GEARSupernovaeIa:
- exponent: -0.35 # Exponent for the distribution of companions
- min_mass_white_dwarf_progenitor: 3 # Minimal mass of a progenitor of white dwarf
- max_mass_white_dwarf_progenitor: 8 # Maximal mass of a progenitor of white dwarf
- max_mass_red_giant: 1.5 # Maximal mass for a red giant
- min_mass_red_giant: 0.9 # Minimal mass for a red giant
- coef_red_giant: 0.02 # Coefficient for the distribution of red giants companions
- max_mass_main_sequence: 2.6 # Maximal mass for a main sequence star
- min_mass_main_sequence: 1.8 # Minimal mass for a main sequence star
- coef_main_sequence: 0.05 # Coefficient for the distribution of main sequence companions
- white_dwarf_mass: 1.38 # Mass of a white dwarf
- GEARSupernovaeII:
- interpolation_size: 200 # Number of elements for the interpolation of the data
-
-Initial Conditions
-~~~~~~~~~~~~~~~~~~
-
-Note that if in the initial conditions, the time of formation of a stellar particle is given (``BirthTime``)
-and set to a negative value, the stellar particle will provide no feedback.
-A similar behavior will be obtained if the parameter ``overwrite_birth_time`` is set to 1 and
-``birth_time`` to -1.
diff --git a/doc/RTD/source/SubgridModels/GEAR/output.rst b/doc/RTD/source/SubgridModels/GEAR/output.rst
new file mode 100644
index 0000000000000000000000000000000000000000..e534bb52ac7f3e12299195f55688e477f762e9e1
--- /dev/null
+++ b/doc/RTD/source/SubgridModels/GEAR/output.rst
@@ -0,0 +1,115 @@
+.. 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`` |
++------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
diff --git a/doc/RTD/source/SubgridModels/GEAR/sinks/index.rst b/doc/RTD/source/SubgridModels/GEAR/sinks/index.rst
index 699dfa53dcdb7991c914bcc4c5c91ff2382758ea..aa4313dda9ee78a0c9b22d57c98aeac7dca2aa24 100644
--- a/doc/RTD/source/SubgridModels/GEAR/sinks/index.rst
+++ b/doc/RTD/source/SubgridModels/GEAR/sinks/index.rst
@@ -15,4 +15,3 @@ Sink particles and star formation in GEAR model
theory
sink_timesteps
params
- output
diff --git a/doc/RTD/source/SubgridModels/GEAR/sinks/output.rst b/doc/RTD/source/SubgridModels/GEAR/sinks/output.rst
deleted file mode 100644
index fda4c6334cbef571c55cbb4d562830a9eef7b9c6..0000000000000000000000000000000000000000
--- a/doc/RTD/source/SubgridModels/GEAR/sinks/output.rst
+++ /dev/null
@@ -1,73 +0,0 @@
-.. 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 and particle IDs.
-
-Sink particles
-~~~~~~~~~~~~~~
-
-+---------------------------------------+-------------------------------------+-----------+---------------------------------------------------+
-| Name | Description | Units | Comments |
-+=======================================+=====================================+===========+===================================================+
-| ``NumberOfSinkSwallows`` | | Number of sink merger events | [-] | |
-| | | | | |
-+---------------------------------------+-------------------------------------+-----------+---------------------------------------------------+
-| ``NumberOfGasSwallows`` | | Number of gas swallowed | [-] | |
-| | | | | |
-+---------------------------------------+-------------------------------------+-----------+---------------------------------------------------+
-| ``TargetMass`` | | Sink target mass to spawn the | [U_M] | | You can use it to determine if the target mass |
-| | | next 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 the | [-] | |
-| | | the sink particles | | |
-| | | | | |
-+---------------------------------------+-------------------------------------+-----------+---------------------------------------------------+
-| ``SwallowedAngularMomentum`` | | Total angular momentum swallowed | [U_M U_L | |
-| | | by the sink particles | ^2 U_T | |
-| | | | ^-1] | |
-+---------------------------------------+-------------------------------------+-----------+---------------------------------------------------+
-| ``MetalMassFractions`` | | Mass fraction of each metal | [-] | | Array of length ``N`` for each particles. The |
-| | | element | | | number of elements ``N`` is determined at |
-| | | | | | compile time by ``--with-chemistry=GEAR_N``. |
-+---------------------------------------+-------------------------------------+-----------+---------------------------------------------------+
-| ``BirthScaleFactors`` | | Scale-factors when the sinks were | [-] | Only used in cosmological runs. |
-| | | born | | |
-| | | | | |
-+---------------------------------------+-------------------------------------+-----------+---------------------------------------------------+
-| ``BirthTimes`` | | Time when the sinks were | [U_T] | Only used in non-cosmological runs. |
-| | | born | | |
-+---------------------------------------+-------------------------------------+-----------+---------------------------------------------------+
-
-
-
-Stars
-~~~~~
-
-+---------------------------------------+-------------------------------------+-----------+---------------------------------------------------+
-| Name | Description | Units | Comments |
-+=======================================+=====================================+===========+===================================================+
-| ``BirthScaleFactors`` | | Scale-factors when the stars were | [-] | Only used in cosmological runs. |
-| | | born | | |
-| | | | | |
-+---------------------------------------+-------------------------------------+-----------+---------------------------------------------------+
-| ``BirthTimes`` | | Time when the stars were | [U_T] | Only used in non-cosmological runs. |
-| | | born | | |
-+---------------------------------------+-------------------------------------+-----------+---------------------------------------------------+
-| ``BirthMasses`` | | Masses of the stars at brith time | [U_M] | |
-| | | | | |
-+---------------------------------------+-------------------------------------+-----------+---------------------------------------------------+
-| ``ProgenitorIDs`` | | ID of the progenitor sinks | [-] | |
-| | | | | |
-+---------------------------------------+-------------------------------------+-----------+---------------------------------------------------+
-| ``StellarParticleType`` | | Type of the stellar particle: | [-] | | The last type correspond to star particles in |
-| | | 0: (discrete) single star | | | the previous model, i.e. representing the full |
-| | | 1: continuous IMF part star | | | IMF |
-| | | 2: single population star | | | |
-+---------------------------------------+-------------------------------------+-----------+---------------------------------------------------+
diff --git a/doc/RTD/source/SubgridModels/GEAR/sinks/theory.rst b/doc/RTD/source/SubgridModels/GEAR/sinks/theory.rst
index cbb8337a5bf3afc665f10feb76f2f9d5a26e9dd3..5991e330666414e33c2dfb4aff698cfbc95c9c76 100644
--- a/doc/RTD/source/SubgridModels/GEAR/sinks/theory.rst
+++ b/doc/RTD/source/SubgridModels/GEAR/sinks/theory.rst
@@ -239,7 +239,7 @@ 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_feedback` section. Here, we only provide the changes from the previous model.
+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.
diff --git a/doc/RTD/source/SubgridModels/GEAR/supernova_feedback.rst b/doc/RTD/source/SubgridModels/GEAR/supernova_feedback.rst
new file mode 100644
index 0000000000000000000000000000000000000000..9a83d47c21da93165c3d3fee5edb0fe8cb7c90c6
--- /dev/null
+++ b/doc/RTD/source/SubgridModels/GEAR/supernova_feedback.rst
@@ -0,0 +1,40 @@
+.. 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.
diff --git a/doc/RTD/source/SubgridModels/index.rst b/doc/RTD/source/SubgridModels/index.rst
index 69daf0eb2c0435b488b3ce9cb6b7c971d615787f..017614e5fcc3adc84fcabd48bc4ce0945aa35878 100644
--- a/doc/RTD/source/SubgridModels/index.rst
+++ b/doc/RTD/source/SubgridModels/index.rst
@@ -16,7 +16,6 @@ be use as an empty canvas to be copied to create additional models.
EAGLE/index
QuickLymanAlpha/index
GEAR/index
- GEAR/sinks/index
Basic/index
AGNSpinJets/index
AGORA/index