doc/RTD/source/Planetary/index.rst

 .. Planetary Simulations
-   Jacob Kegerreis, 3rd October 2020
+   Jacob Kegerreis, 14th July 2022
 
 .. _planetary:
-   
+
 Planetary Simulations
 =====================
 
 SWIFT is also designed for running planetary simulations
-with a current focus on giant impacts, as introduced in 
-`Kegerreis et al. (2019) <https://doi.org/10.1093/mnras/stz1606>`_, MNRAS 487:4.
+such as of giant impacts, as introduced in
+`Kegerreis et al. (2019)  <https://doi.org/10.1093/mnras/stz1606>`_,
+and of any non-planetary simulations with more complicated equations of state 
+and/or multiple materials, etc.
 
-Active development is ongoing of more features and examples for planetary 
-simulations, so please let us know if you are interested in using SWIFT 
-or have any implementation requests. 
+Active development is ongoing of more features, supporting tools, and examples,
+so please let us know if you are interested in using SWIFT
+or have any implementation requests.
 
-You can find an example simulation in ``swiftsim/examples/Planetary/``
-under ``EarthImpact/``.
-The tabulated equation of state files can be downloaded using 
+You can find an example of creating initial conditions and running an impact
+simulation in ``swiftsim/examples/Planetary/`` under ``DemoImpactInitCond/``
+and ``DemoImpact/``. Check out their `README.md` files for more info.
+The tabulated equation of state files can be downloaded there using
 ``EoSTables/get_eos_tables.sh``.
 
-Planetary simulations are currently intended to be run with 
-SWIFT configured to use the planetary hydrodynamics scheme and equations of state:
-``--with-hydro=planetary`` and ``--with-equation-of-state=planetary``.
-These allow for multiple materials to be used,
-chosen from the several available equations of state.
+To run planetary or similar simulations with a traditional SPH formulation,
+SWIFT should be configured to use the planetary hydrodynamics scheme and
+equations of state: ``--with-hydro=planetary`` and
+``--with-equation-of-state=planetary``. These allow for multiple materials to be
+used, chosen from the several available equations of state. 
+
+Planetary simulations can also be carried out using the advanced 
+:ref:`remix_sph` scheme, which improves the treatment of mixing and material 
+interfaces (Sandnes et al. 2025). 
+To configure within REMIX, use ``--with-hydro=remix`` and
+``--with-equation-of-state=planetary``. Note that REMIX simulations require
+initial particle densities in the initial conditions. We also recommend
+configuring with ``--with-kernel=wendland-C2`` and with ``resolution_eta: 1.487``
+in the parameter file for improved hydrodynamic behaviour and since this is the
+configuration used for the validation simulations of Sandnes et al. (2025).
 
 .. toctree::
-   :maxdepth: 2
+   :maxdepth: 1
    :caption: More information:
-   
+
    Hydro Scheme <hydro_scheme>
    Equations of State <equations_of_state>
    Initial Conditions <initial_conditions>
+   Removed Particles <removed_particles>

doc/RTD/source/Planetary/initial_conditions.rst

 .. Planetary Initial Conditions
-   Jacob Kegerreis, 13th March 2020
+   Jacob Kegerreis, 14th July 2022
 
 .. _planetary_initial_conditions:
-   
+
 Initial Conditions
 ==================
 
-Tools for the creation of initial conditions are available via 
-the 
-`WoMa <https://github.com/srbonilla/WoMa>`_ and 
-`SEAGen <https://github.com/jkeger/seagen>`_ open-source python packages, 
+Tools for the creation of initial conditions are available via the
+`WoMa <https://github.com/srbonilla/WoMa>`_ and
+`SEAGen <https://github.com/jkeger/seagen>`_ open-source python packages,
 including: creating spherical or spinning planetary (or similar) profiles;
 placing particles to match arbitrary profiles with precise SPH densities;
 and setting the initial target and impactor positions and velocities,
-as presented in 
+as presented in
 `Kegerreis et al. (2019)  <https://doi.org/10.1093/mnras/stz1606>`_ and
-Ruiz-Bonilla et al. (2020).
+`Ruiz-Bonilla et al. (2020)  <https://doi.org/10.1093/mnras/staa3385>`_ .
 
-They are available with documentation and examples at 
+They are available with documentation and examples at
 https://github.com/srbonilla/WoMa and https://github.com/jkeger/seagen,
 or can be installed directly with ``pip``
-(https://pypi.org/project/woma/, https://pypi.org/project/seagen/).
\ No newline at end of file
+(https://pypi.org/project/woma/, https://pypi.org/project/seagen/).
+
+
+Settling initial conditions with fixed entropies
+------------------------------------------------
+
+If the particles' equations of state include specific entropies,
+and the initial conditions file includes specific entropies for each particle
+(in ``PartType0/Entropies``),
+then configuring SWIFT with ``--enable-planetary-fixed-entropy``
+will override the internal energy of each particle each step such that its
+specific entropy remains constant.
+
+This should be used with caution, but may be a convenient way to maintain an
+entropy profile while initial conditions settle to equilibrium with their
+slightly different SPH densities.

doc/RTD/source/Planetary/removed_particles.rst

+.. Planetary Removed Particles
+   Jacob Kegerreis, 27th November 2020
+
+.. _planetary_removed_particles:
+   
+Removed Particles
+=================
+
+If a particle leaves a non-periodic simulation box then it is removed. 
+Currently, the particle's information is simply printed to stdout in a csv 
+format (see ``hydro_remove_part()`` in ``src/hydro/Planetary/hydro.h``).
+This could be extracted for analysis e.g. using grep:
+``grep '## Removed' -A 1 --no-group-separator output.txt > removed.txt``, 
+then read e.g. with python:
+
+.. code-block:: python
+
+    import numpy as np
+    id, x, y, z, vx, vy, vz, mass, u, P, rho, h, mat_id, time = np.loadtxt(
+        "removed.txt", delimiter=",", unpack=True
+    )
+    pos = np.transpose([x, y, z])
+    vel = np.transpose([vx, vy, vz])
+    id = id.astype(int)
+    mat_id = mat_id.astype(int)
\ No newline at end of file

doc/RTD/source/RadiativeTransfer/GEAR_RT.rst

+.. GEAR Radiative Transfer
+    Mladen Ivkovic 05.2021
+
+.. _rt_GEAR:
+   
+GEAR RT
+-------
+
+.. warning::
+    The radiative transfer schemes are still in development and are not useable
+    at this moment. This page is currently a placeholder to document new
+    features and requirements as the code grows.
+
+
+Compiling for GEAR RT
+~~~~~~~~~~~~~~~~~~~~~
+
+-   To compile swift to be able to run with GEAR RT, you need to configure with
+    ``--with-rt=GEAR_N`` where ``N`` is the integer number of photon groups that 
+    you intend to use in your simulation.
+
+-   You need to choose a Riemann solver for the RT equations. You can choose
+    between the ``GLF`` and ``HLL`` solver. For the time being, I recommend 
+    sticking to the ``GLF`` solver as the ``HLL`` solver is more expensive,
+    but seemingly offers no advantage, although this remains to be confirmed
+    in further testing.
+
+-   GEAR RT is only compatible with the Meshless Finite Volume scheme. You'll
+    need to compile using ``--with-hydro=gizmo-mfv``, which will also require
+    you to select a hydro Riemann solver, e.g ``--with-riemann-solver=hllc``.
+
+-   The thermochemistry requires the `grackle <https://github.com/grackle-project/grackle>`_ 
+    library version above 3.2.1. [#f4]_ Grackle is a chemistry and cooling library presented in 
+    `B. Smith et al. 2017 <https://ui.adsabs.harvard.edu/abs/2017MNRAS.466.2217S>`_.
+    Please note that the current implementation is not (yet) as
+    advanced as the :ref:`GEAR subgrid model grackle cooling <gear_grackle_cooling>`, 
+    and the parameters listed as available there are not applicable for the 
+    grackle cooling in combination with GEAR RT. You can however follow the Grackle 
+    installation instructions documented there.
+
+.. 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
+    GEAR-RT. 
+
+    Additionally, that repository hosts files necessary to install that specific 
+    version of grackle with spack.
+
+
+
+
+Compulsory Runtime Parameters
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+You need to provide the following runtime parameters in the yaml file (which will 
+be further explained below):
+
+.. code:: yaml
+
+   GEARRT:
+       photon_groups_Hz: [3.288e15, 5.945e15, 13.157e15]  # Photon frequency group bin edges in Hz
+       stellar_spectrum_type: 0                           # Which radiation spectrum to use. 
+                                                          #   0: constant. 
+                                                          #   1: blackbody spectrum.
+       stellar_luminosity_model: const                    # Which luminosity model to use.
+       const_stellar_luminosities_LSol: [1., 1., 1.]      # stellar emission rates for each photon 
+                                                          #   frequency bin in units of solar luminosity
+                                                          #   for the 'const' luminosity model
+       f_reduce_c: 1e-3                                   # reduce the speed of light by this factor
+       CFL_condition: 0.9                                 # CFL condition for time integration
+       hydrogen_mass_fraction:  0.76                      # total hydrogen (H + H+) mass fraction in the 
+                                                          # metal-free portion of the gas
+
+
+   TimeIntegration:
+       max_nr_rt_subcycles: 128         # maximal number of RT subcycles per hydro step
+
+
+The ``photon_groups_Hz`` need to be ``N`` frequency edges (floats) to separate 
+the spectrum into ``N`` groups, where ``N`` is the same number you configured
+with using ``--with_rt=GEAR_N``. The edges are **lower** edges of the bins, and
+need to be sorted in increasing order. The final upper edge is defined in a 
+different manner, and depends on the stellar spectrum type you assume (see below
+for more details).
+
+To specify the radiation emitted by stars, there are two main parameters:
+``stellar_luminosity_model`` defines which model to use to obtain star 
+luminosities, while ``stellar_spectrum_type`` determines the spectrum of the
+radiation.
+At the moment, the only way to define star emission rates is to use constant
+stellar luminosities by setting ``stellar_luminosity_model: const``. [#f3]_
+The constant star emission rates need to be provided in the parameter file and
+to be defined for each photon frequency group individually using the 
+``const_stellar_luminosities_LSol`` parameter. The luminosities are expected to
+be in units of solar luminosities. Each star particle will then emit the given 
+luminosities, independent of their other properties, e.g. the stellar age, 
+metallicity, redshift, etc.
+
+When solving the thermochemistry, we need to assume some form of stellar
+spectrum so we may integrate over frequency bins to obtain average interaction
+rates. The parameter ``stellar_spectrum_type`` is hence required, and allows you
+to select between:
+
+- constant spectrum (``stellar_spectrum_type: 0``)
+    - Assume same energy density for any frequency.
+    - This choice additionally requires you to provide a maximal frequency for
+      the spectrum after which it'll be cut off via the 
+      ``stellar_spectrum_const_max_frequency_Hz`` parameter
+
+- blackbody spectrum (``stellar_spectrum_type: 1``)
+    - Assume the spectrum is a blackbody spectrum
+    - In this case, you need to provide also temperature of the blackbody via the 
+      ``stellar_spectrum_blackbody_temperature_K`` parameter.
+    - The assumed maximal considered frequency :math:`\nu_{max}` for this spectrum 
+      is equal to 10 times :math:`\nu_{peak}`, the frequency at which the blackbody 
+      spectrum has its maximum, i.e.
+
+.. math::
+
+     \nu_{peak} = 2.82144 \times k_{B} \times T / h_{Planck}
+
+     \nu_{max} = 10 \times \nu_{peak}
+
+
+.. warning::
+   The ``stellar_spectrum_type`` parameter also determines the averaged photon 
+   interaction cross sections, as they are being computed by integrating a 
+   parametrization of the cross section multiplied by the assumed spectrum. See
+   e.g. equations 9 - 11 in `Rosdahl et al. 2013. 
+   <https://ui.adsabs.harvard.edu/abs/2013MNRAS.436.2188R/abstract>`_
+
+Finally, you will also need to provide an upper threshold for the number of 
+RT-subcycles w.r.t. a single hydro step via ``TimeIntegration:max_nr_rt_subcycles``.
+For more details, refer to :ref:`the subcycling documentation <rt_subcycling>`.
+
+
+
+
+
+
+Optional Runtime Parameters
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+There are several optional parameters which can be set in the ``GEARRT`` block in the
+yaml parameter file:
+
+- ``f_limit_cooling_time``: (Default: ``0.6``) The cooling time computed by grackle 
+  will be  multipled by this factor when estimating the particle's RT time step 
+  size. If set to ``0.0``, the computation of cooling time is turned off.
+- ``skip_thermochemistry``: (Default: ``0``) If set to ``1``, the entire 
+  thermochemistry part of radiative transfer will be skipped. This is intended 
+  only for debugging and testing the radiation transport, as it breaks the 
+  purpose of RT.
+- ``stars_max_timestep``: (Default: ``-1.``) Restrict the maximal timestep of stars 
+  to this value (in internal units). Set to negative to turn off.
+- ``grackle_verbose``: (Default: ``0``) set grackle to verbose.
+- ``case_B_recombination``: (Default: ``1``) If ``1``, use case B recombination rates. 
+  Otherwise, case A recombination rates are used.
+- ``max_tchem_recursion``: (Default: ``0``) If set to positive nonzero value, the
+  thermochemistry is computed using a "10% rule" *in addition to the ones grackle
+  already uses*. If during a thermochemistry step the internal energy :math:``u`` of
+  the gas changes by more than 10%, i.e. :math:`|u_{new}/u_{old} - 1| > 0.1`, then
+  the thermochemistry step is repeated twice using half the time step size. This
+  parameter sets the maximal recursion depth of halving the time step size and 
+  repeating the entire thermochemistry step.
+
+There are some further optional parameters related to setting up initial ion mass
+fractions which are detailed in the 
+:ref:`corresponding section of this documentation <rt_GEAR_set_ion_mass_fractions>`.
+
+
+
+
+
+
+
+
+
+
+Choice of Internal Units
+~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The choice of internal units requires a bit of special attention. Part of the 
+reason is that the exponents of the gas and radiation variables can quickly 
+change by several dozens and cause overflows and other errors. Furthermore, the 
+grackle library may have some other troubles with the units, e.g. when trying to
+find a converging solution. [#f2]_
+
+For this reason, I **strongly encourage** you to run the Internal Units check for 
+GEAR-RT which you can find in the 
+`swiftsim-rt-tools <https://github.com/SWIFTSIM/swiftsim-rt-tools/GEARRTUnitCheck>`_ 
+repository under ``/GEARRTUnitsCheck``. The test should take no more than a 
+minute to run, and requires only two yaml parameter files: the yaml parameter 
+file that you intend to run your simulation with, and one that a provided script 
+can extract automatically from the initial conditions hdf5 file. This test can 
+save you a lot of headaches down the line.
+
+
+
+
+
+Initial Conditions
+~~~~~~~~~~~~~~~~~~
+
+Setting Up Initial Conditions for RT
+````````````````````````````````````
+
+Optionally, you may want to provide initial conditions for the radiation field
+and/or the mass fraction of the ionizing species.
+To do so, you need to add the following datasets to the ``/PartType0`` particle
+group:
+
+.. code:: 
+
+   PhotonEnergiesGroup1
+   PhotonEnergiesGroup2 
+   .
+   .
+   .
+   PhotonEnergiesGroupN
+   PhotonFluxesGroup1
+   PhotonFluxesGroup2
+   .
+   .
+   .
+   PhotonFluxesGroupN
+   MassFractionHI
+   MassFractionHII
+   MassFractionHeI
+   MassFractionHeII
+   MassFractionHeIII
+
+
+-   The ``PhotonEnergies*`` datasets need to have dimension ``nparts``, while the
+    ``PhotonFluxesGroup*`` datasets need to have dimension ``(nparts, 3)``, where
+    ``nparts`` is the number of hydro particles. 
+-   Note that the GEAR-RT scheme expects the ``PhotonEnergies*`` to be total 
+    energies, not energy densities. 
+-   If you are writing initial conditions where the fields have units [#f1]_, then 
+    ``PhotonEnergies*`` are expected to have units of energy 
+    :math:`[M L^2 T^{-2}]`), while the ``PhotonFluxes*`` fields should be in units 
+    of energy times velocity (i.e. energy per unit time per unit area times volume, 
+    :math:`[M L^3 T^{-3}]`).
+-   The ``MassFraction*`` datasets need to have dimension ``nparts`` as well, and
+    are all unitless.
+
+
+Example using Python and ``swiftsimio``
+````````````````````````````````````````
+
+If you are using `swiftsimio <https://github.com/SWIFTSIM/swiftsimio>`_ to write
+the initial condition files, then the easiest way of adding the RT initial
+conditions is to first use the swiftsimio routines to write a file, then open it
+up again and write the additional RT fields again using ``h5py`` routines.
+
+Here is an example:
+
+.. code:: python
+
+    from swiftsimio import Writer
+    import unyt
+    import numpy as np
+    import h5py
+
+    # define unit system to use.
+    unitsystem = unyt.unit_systems.cgs_unit_system
+
+    # number of photon groups
+    nPhotonGroups = 4
+
+    # filename of ICs to be generated
+    outputfilename = "my_rt_ICs.hdf5"
+
+    # open a swiftsimio.Writer object
+    w = Writer(...)
+
+    # do your IC setup for gas, gravity etc now
+    # ... 
+
+    # write the IC file without doing anything RT related.
+    w.write(outputfilename)
+
+    # Now open file back up again and add RT data.
+    F = h5py.File(outputfilename, "r+")
+    header = F["Header"]
+    nparts = header.attrs["NumPart_ThisFile"][0]
+    parts = F["/PartType0"]
+
+    # Create initial photon energies and fluxes. You can leave them unitless, 
+    # the units have already been written down with w.write(). In this case, 
+    # it's in cgs.
+    for grp in range(nPhotonGroups):
+        dsetname = "PhotonEnergiesGroup{0:d}".format(grp + 1)
+        energydata = np.ones((nparts), dtype=np.float32) * some_value_you_want
+        parts.create_dataset(dsetname, data=energydata)
+
+        dsetname = "PhotonFluxesGroup{0:d}".format(grp + 1)
+        fluxdata = np.zeros((nparts, 3), dtype=np.float32) * some_value_you_want
+        parts.create_dataset(dsetname, data=fluxdata)
+
+    # Create initial ionization species mass fractions.     
+    HIdata = np.ones((nparts), dtype=np.float32) * 0.4
+    parts.create_dataset("MassFractionHI", data=HIdata)
+    HIIdata = np.ones((nparts), dtype=np.float32) * 0.1
+    parts.create_dataset("MassFractionHII", data=HIIdata)
+    HeIdata = np.ones((nparts), dtype=np.float32) * 0.3
+    parts.create_dataset("MassFractionHeI", data=HeIdata)
+    HeIIdata = np.ones((nparts), dtype=np.float32) * 0.15
+    parts.create_dataset("MassFractionHeII", data=HeIIdata)
+    HeIIIdata = np.ones((nparts), dtype=np.float32) * 0.05
+    parts.create_dataset("MassFractionHeIII", data=HeIIIdata)
+
+    # close up, and we're done!
+    F.close()
+
+
+
+.. _rt_GEAR_set_ion_mass_fractions:
+
+Generate Ionization Mass Fractions Using SWIFT
+``````````````````````````````````````````````
+
+.. warning:: Using SWIFT to generate initial ionization mass fractions will
+   overwrite the mass fractions that have been read in from the initial 
+   conditions.
+
+Optionally, you can use SWIFT to generate the initial mass fractions of the
+ionizing species. To set the initial mass fractions of all particles to the same
+value, use the following parameters in the yaml parameter file:
+
+.. code:: yaml
+
+    set_initial_ionization_mass_fractions: 1    # (Optional) manually overwrite initial mass fractions 
+                                                # (using the values you set below)
+    mass_fraction_HI: 0.76                      # set initial HI mass fractions to this value
+    mass_fraction_HII: 0.                       # set initial HII mass fractions to this value
+    mass_fraction_HeI: 0.24                     # set initial HeI mass fractions to this value
+    mass_fraction_HeII: 0.                      # set initial HeII mass fractions to this value
+    mass_fraction_HeIII: 0.                     # set initial HeIII mass fractions to this value
+
+Alternatively, you can make SWIFT compute the initial ionization mass fractions
+for you assuming ionization equilibrium, following `Katz, et al. 1996 
+<ui.adsabs.harvard.edu/abs/1996ApJS..105...19K>`_ by setting
+
+.. code:: yaml
+
+    set_equilibrium_initial_ionization_mass_fractions: 1    # (Optional) set the initial ionization fractions 
+                                                            # depending on gas temperature assuming ionization 
+                                                            # equilibrium.
+    hydrogen_mass_fraction:  0.76                           # total hydrogen (H + H+) mass fraction in the 
+                                                            # metal-free portion of the gas
+
+The ``hydrogen_mass_fraction`` (which is a compulsory argument in any case) will
+determine the hydrogen and helium mass fractions, while SWIFT will determine the
+equilibrium ionizations.
+
+.. warning:: If you have somewhat sophisticated initial conditions (e.g. proper galaxies 
+   etc) it is strongly recommended to set up the initial mass fractions to equilibrium
+   values. Otherwise, the initial states can be too far off for GRACKLE's thermochemistry
+   to handle it well, leading to all sorts of troubles, including crashes.
+
+
+
+
+Accessing Output Data
+~~~~~~~~~~~~~~~~~~~~~~
+
+We recommend using `swiftsimio <https://github.com/SWIFTSIM/swiftsimio>`_ to 
+access the RT related snapshot data. The compatibility is being maintained.
+Here's an example how to access some specific quantities that you might find
+useful:
+
+
+.. code:: python
+
+    #!/usr/bin/env python3
+
+    import swiftsimio
+    import unyt
+
+    data = swiftsimio.load("output_0001.hdf5")
+    meta = data.metadata
+
+
+
+    # Accessing RT Related Metadata
+    # ---------------------------------
+
+    # get scheme name: "GEAR M1closure"
+    scheme = str(meta.subgrid_scheme["RT Scheme"].decode("utf-8"))
+
+    # number of photon groups used
+    ngroups = int(meta.subgrid_scheme["PhotonGroupNumber"])
+
+    # get the reduced speed of light that was used. Will have unyts.
+    reduced_speed_of_light = meta.reduced_lightspeed
+
+
+
+
+    # Accessing Photon Data
+    # ------------------------
+
+    # accessing a photon group directly
+    # NOTE: group names start with 1
+    group_1_photon_energies = data.gas.photon_energies.group1
+    group_1_photon_fluxes_x = data.gas.photon_fluxes.Group1X
+    group_1_photon_fluxes_y = data.gas.photon_fluxes.Group1Y
+    group_1_photon_fluxes_z = data.gas.photon_fluxes.Group1Z
+
+    # want to stack all fluxes into 1 array?
+    group1fluxes = swiftsimio.cosmo_array(
+        unyt.uvstack(
+            (group_1_photon_fluxes_x, group_1_photon_fluxes_y, group_1_photon_fluxes_z)
+        ),
+        group_1_photon_fluxes_x.units,
+    ).T
+    # group1fluxes.shape = (npart, 3)
+
+
+    # Load all photon energies in a list
+    photon_energies = [
+        getattr(data.gas.photon_energies, "group" + str(g + 1)) for g in range(ngroups)
+    ]
+
+
+
+    # Accessing Ion Mass Fractions
+    # -------------------------------
+    fHI = data.gas.ion_mass_fractions.HI
+    fHII = data.gas.ion_mass_fractions.HII
+    fHeI = data.gas.ion_mass_fractions.HeI
+    fHeII = data.gas.ion_mass_fractions.HeII
+    fHeIII = data.gas.ion_mass_fractions.HeIII
+
+
+
+
+.. rubric:: Footnotes
+
+.. [#f1] To avoid possible confusions, here are some notes and equations
+   regarding this choice of units.
+
+   One of the RT equations solved by the GEAR RT is the zeroth moment of the
+   equation of radiative transfer for each photon frequency group :math:`i` :
+
+   :math:`\frac{\partial E_i}{\partial t} + \nabla \cdot \mathbf{F}_i = 0`
+
+   where
+
+   - :math:`E_i` : photon energy density; with :math:`[E_i] = erg / cm^3 = M L^{-1} T^{-2}`
+   - :math:`F_i` : radiation flux (energy per unit time per unit surface); with :math:`[F_i] = erg / cm^2 / s = M T^{-3}` 
+
+   and we neglect possible source and sink terms in this footnote.
+
+   These dimensions are also used internally when solving the equations.
+   For the initial conditions however, we require these quantities multiplied by
+   the particle volume. The reason for this choice is so that the photon
+   energies for each particle can be set by the users exactly, while the
+   particle volume computation can be left to SWIFT to worry about internally.
+   The addition of the particle volume term for the radiation flux was made so
+   that the initial conditions are compatible with the SPHM1RT conventions, and
+   both methods can run on the exact same ICs.
+
+
+.. [#f2] For example, choosing cgs units as the internal units may lead to
+   trouble with grackle. (Trouble like a gas at 10^6K without any heating
+   sources heating up instead of cooling down.) The library is set up to work 
+   with units geared towards cosmology. According to Britton Smith (private comm), 
+   a decent rule of thumb is density_units ~ proton mass in g, time_units ~ 1 Myr 
+   to 1 Gyr in s, length_units ~ 1 kpc to 1 Mpc in cm. This should keep you in a 
+   relatively safe range.
+   This is the state of things at 08.2022, with grackle being at version 3.2 (commit
+   ``a089c837b8649c97b53ed3c51c84b1decf5073d8``)
+    
+.. [#f3] Technically there is also the model used for "Test 4" from the 
+   `I. Iliev et al. 2006 <https://ui.adsabs.harvard.edu/abs/2006MNRAS.369.1625I>`_ 
+   paper, but that is very specialized and shouldn't have much use in real 
+   applications.
+
+.. [#f4] Grackle version 3.2.1 still contained a bug related to the use of their
+   "threadsafe functions" that could lead to disastrous outcomes. That was fixed
+   in commit `a59489f`. So versions after 3.2.1 should work as expected.
+   To be safe, we recommend you use the forked grackle repository specifically
+   intended to freeze a stable version for use with SWIFT. You can find that fork
+   on github: https://github.com/mladenivkovic/grackle-swift .

doc/RTD/source/RadiativeTransfer/RT_general.rst

+.. Radiative Transfer Scheme Requirements
+    Mladen Ivkovic 05.2021
+
+.. _rt_general:
+
+
+.. warning::
+    The radiative transfer schemes are still in development and are not useable
+    at this moment. This page is currently a placeholder to document new
+    features and requirements as the code grows.
+
+
+
+
+Radiative Transfer in SWIFT
+---------------------------
+
+We currently support two methods for radiative transfer in SWIFT:
+:ref:`GEAR-RT <rt_GEAR>` and :ref:`SPHM1RT <rt_SPHM1>`. They are both
+moment-based methods. For documentation on each of these methods, please
+refer to their respective pages.
+
+This page contains some general documentation on running SWIFT with RT, which
+is valid irrespective of the exact method used.
+
+
+
+
+
+Requirements
+~~~~~~~~~~~~
+
+
+To be able to run with radiative transfer, you'll need to run swift with the
+following flags:
+
+.. code::
+
+    swift --radiation --hydro --feedback --stars --self-gravity [and/or] --external-gravity
+
+
+
+Some notes on these runtime flags:
+
+
+- The radiation data is coupled to the gas data, so you can't run without
+  ``--hydro``. (Also, the whole point of these schemes is the interaction between
+  gas and radiation, so why would you want to?)
+
+- Currently the only source of radiation are stars, so you need ``--stars``.
+  If you want to run without radiative sources, still run the code with
+  ``--stars``, even if you don't have any stars in your initial conditions.
+
+- Running with ``--stars`` requires some form of gravity, be it self-gravity or
+  external gravity. Since we need stars, we inherit this requirement. If you want
+  no gravity, run with ``--external-gravity`` and set the external potential to
+  zero.
+
+- We need ``--feedback`` in order to have meaningful smoothing lengths for
+  stars. However, you don't need any specific feedback model; It'll work even if
+  you configured ``--with-feedback=none``.
+
+
+
+
+
+.. _rt_subcycling:
+
+RT Sub-Cycling
+~~~~~~~~~~~~~~
+
+SWIFT allows to sub-cycle the solution of radiative transfer steps (both
+photon propagation and thermochemistry) with respect to the hydrodynamics
+time steps. Basically you can tell SWIFT to run up to X radiative transfer
+steps during a single hydrodynamics step for all particles in the simulation.
+The aim is to not waste time doing unnecessary hydrodynamics updates, which
+typically allow for much higher time steps compared to radiation due to the
+propagation speed of the respective advected quantity.
+
+You will need to provide an upper limit on how many RT sub-cycles per hydro
+step you want to allow. That is governed by the
+
+.. code:: yaml
+
+   TimeIntegration:
+       max_nr_rt_subcycles: 128         # maximal number of RT sub-cycles per hydro step
+
+parameter, which is mandatory for any RT runs. To turn off sub-cycling and
+couple the radiative transfer and the hydrodynamics time steps one-to-one,
+set this parameter to either 0 or 1.
+
+Due to the discretization of individual particle time steps in time bins
+with a factor of 2 difference in time step size from a lower to a higher
+time bin, the ``max_nr_rt_subcycles`` parameter itself is required to be
+a power of 2 as well.
+
+Note that this parameter will set an upper limit to the number of sub-cycles
+per hydro step. If the ratio of hydro-to-RT time step is greater than what
+``max_nr_rt_subcycles`` allows for, then the hydro time step will be reduced
+to fit the maximal threshold. If it is smaller, the particle will simply do
+fewer sub-cycles.
+
+
+
+
+
+
+
+Reading the output of time-step information
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+When running with :ref:`sub-cycling enabled <rt_subcycling>`, additional time
+step information is written.
+
+Firstly, whenever you run with a SWIFT compiled with any RT method, it will also
+create a file called ``rtsubcycles.txt``, which  contains analogous data to the
+``timesteps.txt`` file. However, if you run without sub-cycling (i.e. with
+``TimeIntegration:max_nr_rt_subcycles`` set to ``0`` or ``1``), **then the file
+will remain empty** (save for the header information). That's because its contents
+would be redundant - all the data will be identical to what will be written in
+``timesteps.txt``.
+
+
+
+Secondly, when running with RT and sub-cycling enabled, the output to ``STDOUT``
+contains additional information. More concretely, it changes from something like this:
+
+
+.. code-block::
+
+    [00003.0] engine_dump_snapshot: Dumping snapshot at t=0.000000e+00
+    [00003.1] engine_print_stats: Saving statistics at t=0.000000e+00
+    #   Step           Time Scale-factor     Redshift      Time-step Time-bins      Updates    g-Updates    s-Updates sink-Updates    b-Updates  Wall-clock time [ms]  Props    Dead time [ms]
+           0   0.000000e+00    1.0000000    0.0000000   0.000000e+00    1   56        18000        31000        13000            0            0              2610.609    281           101.971
+           1   1.220703e-05    1.0000000    0.0000000   1.220703e-05   43   43           11           11            0            0            0                61.686      1             1.324
+           2   2.441406e-05    1.0000000    0.0000000   1.220703e-05   43   44        12685        12685            0            0            0              1043.433      0            35.461
+           3   3.662109e-05    1.0000000    0.0000000   1.220703e-05   43   43           11           11            0            0            0                51.340      1             1.628
+           4   4.882813e-05    1.0000000    0.0000000   1.220703e-05   43   45        18000        18000            0            0            0              1342.531      0            36.831
+           5   6.103516e-05    1.0000000    0.0000000   1.220703e-05   43   43           11           11            0            0            0                48.412      1             1.325
+           6   7.324219e-05    1.0000000    0.0000000   1.220703e-05   43   44        12685        12685            0            0            0              1037.307      0            34.718
+           7   8.544922e-05    1.0000000    0.0000000   1.220703e-05   43   43           11           11            0            0            0                47.791      1             1.362
+           8   9.765625e-05    1.0000000    0.0000000   1.220703e-05   43   46        18000        18004            4            0            0              1410.851      0            35.005
+           9   1.098633e-04    1.0000000    0.0000000   1.220703e-05   43   43           11           11            0            0            0                48.322      1             1.327
+          10   1.220703e-04    1.0000000    0.0000000   1.220703e-05   43   44        12685        12685            0            0            0              1109.944      0            33.691
+
+
+To something like this:
+
+
+.. code-block::
+
+     [rt-sc] 0    0.000000e+00    1.000000    0.000000  1.220703e-05    1   56        18000            -            -            -            -                     -      -                 -
+     [rt-sc] 1    1.220703e-05    1.000000    0.000000  1.220703e-05   43   43           11            -            -            -            -                 0.683      -             0.170
+     [rt-sc] 2    2.441406e-05    1.000000    0.000000  1.220703e-05   43   44        12685            -            -            -            -               100.543      -             5.070
+     [rt-sc] 3    3.662109e-05    1.000000    0.000000  1.220703e-05   43   43           11            -            -            -            -                 0.762      -             0.208
+     [rt-sc] 4    4.882813e-05    1.000000    0.000000  1.220703e-05   43   45        18000            -            -            -            -               124.011      -             6.396
+     [rt-sc] 5    6.103516e-05    1.000000    0.000000  1.220703e-05   43   43           11            -            -            -            -                 0.831      -             0.254
+     [rt-sc] 6    7.324219e-05    1.000000    0.000000  1.220703e-05   43   44        12685            -            -            -            -               107.172      -             5.572
+     [rt-sc] 7    8.544922e-05    1.000000    0.000000  1.220703e-05   43   43           11            -            -            -            -                 0.759      -             0.227
+    [00003.0] engine_dump_snapshot: Dumping snapshot at t=0.000000e+00
+    [00003.1] engine_print_stats: Saving statistics at t=0.000000e+00
+    #   Step           Time Scale-factor     Redshift      Time-step Time-bins      Updates    g-Updates    s-Updates sink-Updates    b-Updates  Wall-clock time [ms]  Props    Dead time [ms]
+           0   0.000000e+00    1.0000000    0.0000000   0.000000e+00    1   56        18000        31000        13000            0            0              2941.254    281           120.261
+     [rt-sc] 0    9.765625e-05    1.000000    0.000000  1.220703e-05   43   46        18000            -            -            -            -                     -      -                 -
+     [rt-sc] 1    1.098633e-04    1.000000    0.000000  1.220703e-05   43   43           11            -            -            -            -                 0.990      -             0.417
+     [rt-sc] 2    1.220703e-04    1.000000    0.000000  1.220703e-05   43   44        12685            -            -            -            -               104.155      -             5.744
+     [rt-sc] 3    1.342773e-04    1.000000    0.000000  1.220703e-05   43   43           11            -            -            -            -                 0.765      -             0.176
+     [rt-sc] 4    1.464844e-04    1.000000    0.000000  1.220703e-05   43   45        18000            -            -            -            -               125.237      -             5.605
+     [rt-sc] 5    1.586914e-04    1.000000    0.000000  1.220703e-05   43   43           11            -            -            -            -                 0.856      -             0.282
+     [rt-sc] 6    1.708984e-04    1.000000    0.000000  1.220703e-05   43   44        12685            -            -            -            -               112.171      -             5.251
+     [rt-sc] 7    1.831055e-04    1.000000    0.000000  1.220703e-05   43   43           11            -            -            -            -                 0.861      -             0.241
+           1   9.765625e-05    1.0000000    0.0000000   9.765625e-05   46   46            4            8            4            0            0               546.225      1            24.648
+     [rt-sc] 0    1.953125e-04    1.000000    0.000000  1.220703e-05   43   47        18000            -            -            -            -                     -      -                 -
+     [rt-sc] 1    2.075195e-04    1.000000    0.000000  1.220703e-05   43   43           11            -            -            -            -                 0.842      -             0.212
+     [rt-sc] 2    2.197266e-04    1.000000    0.000000  1.220703e-05   43   44        12685            -            -            -            -               126.674      -             6.295
+     [rt-sc] 3    2.319336e-04    1.000000    0.000000  1.220703e-05   43   43           11            -            -            -            -                 0.797      -             0.289
+     [rt-sc] 4    2.441406e-04    1.000000    0.000000  1.220703e-05   43   45        18000            -            -            -            -               142.086      -             5.511
+     [rt-sc] 5    2.563477e-04    1.000000    0.000000  1.220703e-05   43   43           11            -            -            -            -                 0.919      -             0.196
+     [rt-sc] 6    2.685547e-04    1.000000    0.000000  1.220703e-05   43   44        12685            -            -            -            -               131.550      -             5.896
+     [rt-sc] 7    2.807617e-04    1.000000    0.000000  1.220703e-05   43   43           11            -            -            -            -                 0.809      -             0.186
+           2   1.953125e-04    1.0000000    0.0000000   9.765625e-05   46   47           27           43           16            0            0               558.226      0            27.711
+     [rt-sc] 0    2.929688e-04    1.000000    0.000000  1.220703e-05   43   46        18000            -            -            -            -                     -      -                 -
+     [rt-sc] 1    3.051758e-04    1.000000    0.000000  1.220703e-05   43   43           11            -            -            -            -                 0.738      -             0.207
+     [rt-sc] 2    3.173828e-04    1.000000    0.000000  1.220703e-05   43   44        12685            -            -            -            -               122.572      -             5.170
+     [rt-sc] 3    3.295898e-04    1.000000    0.000000  1.220703e-05   43   43           11            -            -            -            -                 1.063      -             0.345
+     [rt-sc] 4    3.417969e-04    1.000000    0.000000  1.220703e-05   43   45        18000            -            -            -            -               147.110      -             5.409
+     [rt-sc] 5    3.540039e-04    1.000000    0.000000  1.220703e-05   43   43           11            -            -            -            -                 1.091      -             0.350
+     [rt-sc] 6    3.662109e-04    1.000000    0.000000  1.220703e-05   43   44        12685            -            -            -            -               134.273      -             6.561
+     [rt-sc] 7    3.784180e-04    1.000000    0.000000  1.220703e-05   43   43           11            -            -            -            -                 0.825      -             0.298
+           3   2.929688e-04    1.0000000    0.0000000   9.765625e-05   46   46            4            8            4            0            0               557.164      0            24.760
+
+
+Here's what's going on here:
+
+- All lines beginning with the prefix ``[rt-sc]`` are time step data of RT
+  sub-cycling steps, i.e. of sub-cycles.
+- The sub-cycling index follows the prefix ``[rt-sc]``. For example, a line
+  beginning with ``[rt-sc] 2`` is the sub-cycle with index ``2`` of some time step.
+- The "sub-cycle" with index ``0`` is the one performed alongside all other tasks
+  (e.g. hydro, gravity, stellar feedback, etc.) All other sub-cycle indices
+  indicate actual sub-cycles, i.e. actual intermediate steps where only radiative
+  transfer is being solved (or put differently: where only RT tasks are being launched).
+- The sub-cycling lines are written to ``STDOUT`` *before* the line of the full
+  time step data. More precisely, in the above example, time step ``1`` with all its
+  RT sub-cycles is written to ``STDOUT`` as this block:
+
+.. code-block::
+
+    #   Step           Time Scale-factor     Redshift      Time-step Time-bins      Updates    g-Updates    s-Updates sink-Updates    b-Updates  Wall-clock time [ms]  Props    Dead time [ms]
+
+    [ ... some lines omitted ... ]
+
+     [rt-sc] 0    9.765625e-05    1.000000    0.000000  1.220703e-05   43   46        18000            -            -            -            -                     -      -                 -
+     [rt-sc] 1    1.098633e-04    1.000000    0.000000  1.220703e-05   43   43           11            -            -            -            -                 0.990      -             0.417
+     [rt-sc] 2    1.220703e-04    1.000000    0.000000  1.220703e-05   43   44        12685            -            -            -            -               104.155      -             5.744
+     [rt-sc] 3    1.342773e-04    1.000000    0.000000  1.220703e-05   43   43           11            -            -            -            -                 0.765      -             0.176
+     [rt-sc] 4    1.464844e-04    1.000000    0.000000  1.220703e-05   43   45        18000            -            -            -            -               125.237      -             5.605
+     [rt-sc] 5    1.586914e-04    1.000000    0.000000  1.220703e-05   43   43           11            -            -            -            -                 0.856      -             0.282
+     [rt-sc] 6    1.708984e-04    1.000000    0.000000  1.220703e-05   43   44        12685            -            -            -            -               112.171      -             5.251
+     [rt-sc] 7    1.831055e-04    1.000000    0.000000  1.220703e-05   43   43           11            -            -            -            -                 0.861      -             0.241
+           1   9.765625e-05    1.0000000    0.0000000   9.765625e-05   46   46            4            8            4            0            0               546.225      1            24.648
+
+- Let's have a closer look at the written data:
+
+  - In step ``1``, 4 hydro particles, 8 gravity particles, and 4 star particles
+    were updated. (You can see that in the last line.)
+  - The integration of the full step was performed over a time step with
+    size ``9.765625e-05``. **That is valid for all physics except radiative
+    transfer.** The RT was integrated 8 times with a time step size of ``1.220703e-05``.
+    You can see this in the sub-cycling output lines.
+  - Each RT sub-cycling line also tells you the minimal and maximal time bin
+    size that was worked on, as well as how many hydro particles underwent RT
+    updates.
+  - RT sub-cycles only ever update radiative transfer. There will never be any
+    gravity, star, sink, or black hole particle updates in it.
+  - Since the sub-cycle with index ``0`` is performed alongside all other physics
+    during the main step, the isolated wall-clock time and dead time fields are
+    not available for it.
+  - The wall-clock time and dead time fields in the full step line (the one starting
+    *without* ``[rt-sc]``) include the data of the sub-cycles for this step as well.
+
+
+
+
+

doc/RTD/source/RadiativeTransfer/RT_notes_for_developers.rst

+.. RT developer notes
+    Mladen Ivkovic 07.2022
+
+.. _rt_dev:
+   
+Notes for Developers
+========================
+
+
+.. _rt_workflow:
+
+The Radiative Transfer Workflow
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+This section is intended to give a superficial overview on the general workflow 
+that the RT tasks are going through. It should be sufficient to give you a general
+idea of how things work, and allow you to plan and implement your own RT scheme
+into SWIFT.
+
+For a more complete documentation on the tasking system used for RT, please refer
+to the :ref:`subsequent section <rt_task_system>`.
+
+.. figure:: RTWorkflow.png
+    :width: 400px
+    :align: center
+    :figclass: align-center
+    :alt: Workflow for the Radiative Transfer in SWIFT.
+
+    This figure shows the general RT workflow in SWIFT.
+    The green nodes are the ones related to the RT.
+
+
+- There are some additions to star tasks: 
+
+  1)    During the star density loop, we gather additional neighbour information
+        that are required for the star-gas interaction.
+
+  2)    In the star ghost tasks, stars compute how much radiation they will 
+        inject in the surrounding gas. The radiation injection scheme follows 
+        the philosophy of the stellar feedback, and only stars that are 
+        currently active will inject radiation into neighbouring gas particles, 
+        regardless of whether the gas particle is currently active (as opposed 
+        to active gas particles "pulling" radiation from possibly inactive 
+        stars). So the each star needs to compute how much radiation it will 
+        emit during its current time step.
+
+  3)    During the stellar feedback interaction loop, the stars inject radiation
+        onto neighbouring gas particles.
+
+- ``Ghost1`` tasks operate on individual particles, and are intended to finish up
+  any leftover work from the injection.
+
+- ``Gradient`` tasks are particle-particle interaction tasks, intended for 
+  particles to gather data from its own neighbours, e.g. so we can estimate the 
+  current local gradients.  This is an interaction of "type 1", meaning that any 
+  particle will only interact with neighbours which are within its own compact 
+  support radius.
+
+- ``Ghost2`` tasks operate on individual particles, and are intended to finish 
+  up any leftover work from the "gradients".
+
+- ``Transport`` tasks are particle-particle interaction tasks, intended to 
+  propagate the radiation. This is an interaction of "type 2", meaning that any 
+  particle will interact with any neighbours within either particles' compact 
+  support radius.
+
+- ``thermochemistry`` tasks operate on individual particles, and are intended
+  to solve the thermochemistry equations.
+
+
+
+
+
+
+
+.. _rt_task_system:
+
+Current Task System
+~~~~~~~~~~~~~~~~~~~~
+
+Some RT tasks featured in the full task graphs below, like the 
+``rt_advance_cell_time``, ``rt_collect_times``, and ``rt_sorts``, have not been 
+mentioned in the previous section. They are necessary for internal machinations 
+of the RT subcycling scheme, and do not affect the RT scheme itself. If you are
+implementing a new RT scheme into SWIFT, you should not need to touch those
+tasks. For more documentation on them, please refer to the :ref:`subsequent
+section <rt_subcycling_documentation>`.
+
+
+.. figure:: RTTaskDependenciesFull-simplified.png
+    :width: 400px
+    :align: center
+    :figclass: align-center
+    :alt: Task dependencies for the sink scheme.
+
+    This figure shows the task dependencies for the radiative transfer scheme.
+    Some tasks with little or no relevance to the RT scheme have been simplified
+    and condensed for clarity.
+    This was done with SWIFT v0.9.0.
+
+
+.. figure:: full_dependency_graph_RT.png
+    :width: 400px
+    :align: center
+    :figclass: align-center
+    :alt: Task dependencies for the sink scheme.
+
+    This figure shows the full task dependencies for the radiative transfer scheme
+    with self-gravity.
+    This was done with SWIFT v0.9.0.
+
+
+
+
+
+
+.. _rt_subcycling_documentation:
+
+Notes on Subcycling
+~~~~~~~~~~~~~~~~~~~~~
+
+Note: This section is directed towards developers and maintainers, not
+necessarily towards users.
+
+How it works
+`````````````````
+
+A subcycle is basically a SWIFT time step where only radiative transfer is being
+run.
+
+After a normal SWIFT time step (i.e. after a call to ``engine_launch()`` and the
+global collection and communication) is complete, the starting time of the
+following global time step is known. We also collect the current minimal RT time 
+step size, which allows us to determine how many sub-cycles we need to complete
+before the next normal SWIFT time step is launched. Particles are not drifted 
+during a subcycle, and the propagation velocity (aka the speed of light) is 
+taken to be constant, so the number of subcycles is fixed at the end of a normal 
+step. For each subcycle, we then unskip the RT tasks, and make a new call to
+``engine_launch()``.
+
+For the time integration to work correctly, the time integration variables of
+particles like the time-bins are kept independently from the hydro ones. The same
+goes for the respective quantities of cells, like the next integer end time of
+the cell, or the minimal RT time step size in the cell. Furthermore, the global
+time variables that are stored in the engine (e.g. current integer time, current
+max active bin...) have a copy that is being kept up-to-date outside of normal
+SWIFT steps in the same manner the non-RT variables are being updated each
+normal step. The RT subcycling scheme never touches or changes any global time
+integration related variable.
+
+Since the time stepping variables of particles and cells are taken to be
+constant during subcycles (because there are no drifts, constant speed of
+light), the ``timestep`` tasks are not being run during a sub-cycle. This
+effectively means that the particle time bins can only be changed in a normal
+step when the particle is also hydro-active. Furthermore, there are no MPI 
+communications after the tasks have finished executing to update any global 
+times etc. for the same reason. There are some functionalities of the
+``timestep`` and the ``collect`` tasks which are still necessary though:
+
+- The ``timestep`` task also updates the cell's next integer end time after it
+  has been determined during the task. During a subcycle, the next end time is
+  simply the current time plus the minimal time step size of the cell, but we
+  need a task to actually update the cell at the end of each subcycle. The
+  ``rt_advance_cell_time`` task does exactly that, and in that sense does the
+  ``timestep`` task's job during subcycles.
+
+- The ``collect`` task propagates sub-cell data like the minimal end time or the
+  RT time step size from the super level to the top level. This functionality is
+  replaced with the ``rt_collect_times`` tasks during subcycles. Note that the
+  ``rt_collect_times`` tasks aren't being activated during normal steps, as the
+  ``collect`` tasks already do the job just fine.
+
+Something special about the ``rt_advance_cell_time`` tasks is that they are
+also created and run on foreign cells. During a subcycle, the ``tend`` tasks
+don't run and don't update the cell time variables from the original cell, so
+during the subsequent unskipping, the data will be wrong, leading to all sorts
+of trouble. We can do that on foreign cells during sub-cycles because all the
+cell's time step sizes stay fixed between two regular SWIFT steps, and hence
+the number of sub-cycles all the sub-cycles' end times are predictable.
+
+
+
+
+RT Sorts
+````````````````
+
+The sorting of particles required for pair-type interaction tasks requires some
+special attention. The issues arise because a subcycle step of a cell can
+coincide with the main step of another cell. To illustrate, suppose we have two
+cells, ``A`` and ``B``. Let cell ``A`` have a hydro time step of size 4, and 
+cell ``B`` a hydro time step of size 8. Let both cells do 2 RT subcycles per
+hydro step each. In the graph below, an ``X`` represents when a cell will be
+updated:
+
+.. code::
+
+   Cell A
+     Hydro active:    X               X               X               X               X       
+     RT active:       X       X       X       X       X       X       X       X       X       X
+
+   Cell B
+     Hydro active:    X                               X                               X
+     RT active:       X               X               X               X               X       
+    
+    
+    ------------------|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
+    t                 0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18
+
+
+Note that e.g. at ``t`` = 4, cell ``B`` is only RT active, while cell ``A`` is
+also hydro active. It being hydro active means that we will have a SWIFT main
+step at that time.
+
+Now suppose cell cells ``A`` and ``B`` are neighbouring cells that undergo hydro
+interactions, but are on different MPI ranks. For the hydro interactions in the
+normal SWIFT step at ``t`` = 4, cell ``B`` will be sent over to the rank of 
+cell ``A``. Once it was received, it will be sorted, because after being 
+received, the ``recv_xv`` task resets the arrived cell's sort tasks, and a hydro 
+sort is activated. This is the default hydrodynamics workflow.
+
+
+
+Complications however arise when several conditions coincide:
+
+- a foreign cell has been drifted on its "home" domain due to some reason other 
+  than hydro, e.g. for gravity or for stellar feedback. 
+- the foreign cell and all of its neighbours are not hydro active at the current
+  main step, so no ``recv_xv`` task has been run on the local rank, no ``sort``
+  task has been run, and the sorting flags have not been reset for the cell.
+- the foreign cell is involved in an RT interaction that coincides with the 
+  current main step, i.e. the cell or one of its neighbours is RT active during
+  a main step. (E.g. like cell ``A`` at ``t`` = 2 in the graph above.)
+
+If these conditions are met, the cell will undergo an interaction while
+unsorted. Obviously that's a no-no.
+
+
+To illustrate this problem, consider the following scenario. Say we have cells ``A``,
+``B``, and ``C``. Cell ``A`` is "home" on MPI rank 0, while cell ``B`` is on MPI rank 1. 
+In the current step in this scenario, cell ``A`` and ``B`` are both active for RT, and 
+interact with each other. ``C`` has an active star, which inject energy into particles of
+``B``.  Cell ``A`` and ``C`` do *not* interact with each other:
+
+
+.. code::
+
+   rank 0       | rank 1  
+                |
+      RT interaction        Star Feedback  
+   A <----------|---> B <------------------ C
+                |
+                |
+
+
+Since ``A`` and ``B`` interact, but are separated between different MPI ranks, both ``A``
+and ``B`` will have local copies of each other available on their local ranks, respectively.
+These cells are referred to as "foreign cells". Let's denote them with an additional ``f``:
+
+.. code::
+
+    rank 0                      | rank 1  
+                                |
+             RT interaction     |        RT interaction      Star Feedback  
+           A <-------------> Bf |   Af  <-------------> B <------------------ C
+                             ^  |   ^
+    (this is a foreign cell) |  |   | (this is a foreign cell)
+
+
+Now let's first have a look at the default workflow when hydrodynamics is involved. 
+Assume that all cells ``A``, ``B``, and ``C`` are hydro *and* RT active for this step.
+Once again cells ``A`` and ``B`` interact, and ``B`` and ``C`` interact, but ``A`` 
+and ``C`` *don't* interact. Cell ``C`` contains an active star particle.
+
+**Without** MPI, the order of operations for each cell would look like this: (operations 
+where two cells interact with each other are marked with an arrow)
+
+
+.. code::
+
+    A                   B                  C
+    ----------------------------------------------------
+    Drift               Drift               Drift
+    Sort                Sort                Sort
+    Density Loop <----> Density Loop <----> Density Loop
+    Ghost               Ghost               Ghost
+    Force Loop <------> Force Loop <------> Force Loop
+    End Force           End Force           End Force
+    Kick2               Kick2               Kick2
+
+                                            Star Drift
+                                            Star Sort
+                        (inactive) <------> Star Density
+                                            Star Ghost
+                        (inactive) <------> Star Feedback
+
+    RT Ghost1           RT Ghost1           RT Ghost1
+    RT Gradient <-----> RT Gradient <-----> RT Gradient
+    RT Ghost2           RT Ghost2           RT Ghost2
+    RT Transport <----> RT Transport <----> RT Transport
+    RT Tchem            RT Tchem            RT Tchem
+    Timestep            Timestep            Timestep
+    Kick1               Kick1               Kick1
+
+
+Now  **with** MPI communications, cells ``A`` and ``B`` need to send over the 
+up-to-date data to their foreign counterparts ``Af`` and ``Bf``, respectively, 
+*before* each interaction type task (the ones with arrows in the sketch above). 
+The order of operations should look like this:
+(The foreign cell ``Af`` on rank 1 is omitted for clarity, but follows the same
+principle as ``Bf`` on rank 0)
+
+.. code::
+
+    rank 0                                  |  rank 1
+    A                   Bf                  |  B                   C
+    ----------------------------------------|----------------------------------------
+    Drift                                   |  Drift               Drift
+                        Recv XV  <------------ Send XV
+    Sort                Sort                |  Sort                Sort
+    Density Loop <----> Density Loop        |  Density Loop <----> Density Loop
+    Ghost                                   |  Ghost               Ghost
+                        Recv Density <-------- Send Density
+    Force Loop <------> Force Loop          |  Force Loop <------> Force Loop
+    End Force                               |  End Force           End Force
+    Kick2                                   |  Kick2               Kick2
+                                            |
+                                            |                      Star Drift
+                                            |                      Star Sort
+                                            |  (inactive) <------> Star Density
+                                            |                      Star Ghost
+                                            |  (inactive) <------> Star Feedback
+                                            |
+    RT Ghost1                               |  RT Ghost1           RT Ghost1
+                        Recv RT Gradient <---- Send RT Gradient
+    RT Gradient <-----> RT Gradient         |  RT Gradient <-----> RT Gradient
+    RT Ghost2                               |  RT Ghost2           RT Ghost2
+                        Recv RT Transport <--- Send RT Transport
+    RT Transport <----> RT Transport        |  RT Transport <----> RT Transport
+    RT Tchem                                |  RT Tchem            RT Tchem
+    Timestep                                |  Timestep            Timestep
+                        Recv tend <----------- Send tend
+    Kick1                                   |  Kick1               Kick1
+
+
+Finally, let's look at the scenario which causes problems with the sort. This 
+is the case, as described above, when (a) cells ``A`` and ``B`` are RT active during
+a main step (like in the sketch above), (b) aren't hydro active during a main step 
+(unlike what is drawn above), (c) one of these cells is foreign (in this case, ``Bf``),
+while the "home" cell (cell ``B``) get drifted during a main step for some reason
+other than hydrodynamics, e.g. because a star interaction with cell ``C`` requested it.
+
+In this case, the workflow looks like this:
+
+.. code::
+
+    rank 0                                  |  rank 1
+    A                   Bf                  |  B                   C
+    ----------------------------------------|----------------------------------------
+                                            |  Drift               Drift
+                                            |  Sort                Sort
+                                            |                      Kick2
+                                            |
+                                            |                      Star Drift
+                                            |                      Star Sort
+                                            |  (inactive) <------> Star Density
+                                            |                      Star Ghost
+                                            |  (inactive) <------> Star Feedback
+                                            |
+    RT Ghost1                               |  RT Ghost1           RT Ghost1
+                        Recv RT Gradient <---- Send RT Gradient
+    RT Gradient <-----> RT Gradient         |  RT Gradient <-----> RT Gradient
+    RT Ghost2                               |  RT Ghost2           RT Ghost2
+                        Recv RT Transport <--- Send RT Transport
+    RT Transport <----> RT Transport        |  RT Transport <----> RT Transport
+    RT Tchem                                |  RT Tchem            RT Tchem
+    Timestep                                |  Timestep            Timestep
+                        Recv tend <----------- Send tend
+    Kick1                                   |  Kick1               Kick1
+
+The issue is that with the missing hydro communication tasks, the first communication
+between cell ``B`` and its foreign counterpart ``Bf`` is the ``recv_rt_gradient`` task.
+Recall that during a communication, we always send over all particle data of a cell.
+This includes all the particle positions, which may have been updated during a drift.
+However, the sorting information is not stored in particles, but in the cell itself.
+For this reason, a ``sort`` task is *always* run directly after a cell finishes the 
+``recv_xv`` task, which, until the sub-cycling was added, was always the first task any 
+foreign cell would run, with a subsequent ``sort``.
+
+The RT sub-cycling now however allows the ``recv_xv`` task to not run at all
+during a main step, since it's not always necessary, as shown in the example above.
+All the required data for the RT interactions can be sent over with
+``send/recv_rt_gradient`` tasks. An unintended consequence however is that in a
+scenario as sketched above, at the time of the ``A <---> Bf`` RT Gradient
+interaction, cell ``Bf`` will not be sorted. That's a problem.
+
+To solve this issue, a new task has been added, named ``rt_sorts``. It is only
+required for foreign cells, like cell ``Bf``, and only during normal/main steps 
+(as we don't drift during subcycles, there won't be any reason to re-sort.) On
+local cells, each time a drift is activated for an interaction type task, the 
+sort task is also activated. So there is no need for ``rt_sorts`` tasks on local
+cells.
+An additional advantage to adding a new sort task like the ``rt_sorts`` is that
+it allows us to sidestep possible deadlocks. Suppose that as an alternative to
+the ``rt_sort`` tasks we instead use the regular hydro ``sort`` task. The default
+hydro ``sort`` task is set up to run before the other hydro tasks, and in 
+particular before the ``kick2`` task. However the RT and star tasks are executed 
+*after* the ``kick2``. This means that there are scenarios where a cell with a 
+foreign counterpart like cells ``B`` and ``Bf`` can deadlock when ``Bf`` is waiting
+for the ``recv_rt_gradient`` to arrive so it may sort the data, while ``B`` is
+waiting for ``Bf`` to finish the sorting and proceed past the ``kick2`` stage so
+it can run the ``send_rt_gradient`` data which would allow ``Bf`` to run the
+sorts.
+
+
+The ``rt_sorts`` tasks are executed after the first RT related ``recv``, in this 
+case the ``recv_rt_gradient``. The order of operations should now look like this:
+
+
+.. code::
+
+    rank 0                                  |  rank 1
+    A                   Bf                  |  B                   C
+    ----------------------------------------|----------------------------------------
+                                            |  Drift               Drift
+                                            |  Sort                Sort
+                                            |                      Kick2
+                                            |
+                                            |                      Star Drift
+                                            |                      Star Sort
+                                            |  (inactive) <------> Star Density
+                                            |                      Star Ghost
+                                            |  (inactive) <------> Star Feedback
+                                            |
+    RT Ghost1                               |  RT Ghost1           RT Ghost1
+                        Recv RT Gradient <---- Send RT Gradient
+                        rt_sort             |
+    RT Gradient <-----> RT Gradient         |  RT Gradient <-----> RT Gradient
+    RT Ghost2                               |  RT Ghost2           RT Ghost2
+                        Recv RT Transport <--- Send RT Transport
+    RT Transport <----> RT Transport        |  RT Transport <----> RT Transport
+    RT Tchem                                |  RT Tchem            RT Tchem
+    Timestep                                |  Timestep            Timestep
+                        Recv tend <----------- Send tend
+    Kick1                                   |  Kick1               Kick1
+
+
+
+In order to minimize unnecessary work, three new cell flags concerning the RT
+sorts have been added:
+ 
+- ``cell_flag_do_rt_sub_sort``: tracks whether we need an RT sub sort, which is 
+  equivalent to the ``cell_flag_do_sub_sort`` flag for hydro. We can't use the 
+  hydro flag though because the hydro flag is also used to early-exit walking up 
+  the cell hierarchy when activating hydro subcell sorts. In particular, this
+  condition in ``cell_unskip.c:cell_activate_hydro_sorts_up():``
+
+.. code:: 
+
+   void cell_activate_hydro_sorts_up(struct cell *c, struct scheduler *s) {
+       /* omitted lines */
+
+       for (struct cell *parent = c->parent;
+           parent != null && !cell_get_flag(parent, cell_flag_do_hydro_sub_sort); parent = parent->parent) {
+       /* !! this is the problem ---^ */
+   
+       /* omitted lines */
+     }
+   }
+
+The sort activation for RT and for hydro can run concurrently. So there is no
+guarantee that when the hydro sort activation sets the
+``cell_flag_do_hydro_sub_sort`` flag, the RT sorting tasks will be activated
+correctly, which occurs at the top of the cell hierarchy tree walk.
+So we need an independent flag for RT here to not abort the tree walk early
+and in error.
+
+- ``cell_flag_do_rt_sort``: tracks whether the call to the 
+  ``runner_do_hydro_sort()`` function was requested by an RT sort. (Both the (hydro) 
+  ``sort`` and the ``rt_sort`` tasks call the same function.) This flag is used to 
+  discriminate during the actual sorting in the ``runner_do_hydro_sort()``
+  function if the internal check whether the cell is drifted to the current time
+  may be disregarded. When an RT subcycle coincides with a main step, the particles 
+  won't necessarily be drifted to the current time as there is no need to drift them 
+  for RT only. This is intended behaviour, so we allow ``runner_do_hydro_sort()`` to
+  skip this drift check in the case where the sorting was requested for RT purposes.
+
+- ``cell_flag_skip_rt_sort``: Tracks whether a regular hydro sort has been
+  activated for this cell. If it has, then there is no need to run an RT sort as
+  well, and we skip it.
+

doc/RTD/source/RadiativeTransfer/SPHM1_RT.rst

+.. SPHM1RT Radiative Transfer
+    Tsang Keung Chan 01.2022
+
+.. _rt_SPHM1:
+   
+SPHM1 RT
+--------
+
+SPHM1RT is the first two-moment radiative transfer on smoothed particle hydrodynamics (`Chan et al. 2021
+<https://ui.adsabs.harvard.edu/abs/2021MNRAS.505.5784C/abstract>`_). It solves the radiation energy and flux equations with a modified Eddington tensor closure. It is adaptive, efficient, and easy to parallelize.
+
+.. warning::
+    The radiative transfer schemes are still in development and are not useable
+    at this moment. This page is currently a placeholder to document new
+    features and requirements as the code grows.
+
+
+Compiling for SPHM1-RT
+~~~~~~~~~~~~~~~~~~~~~~
+
+-   To compile swift to be able to run with SPHM1-RT, you need to configure with
+    ``--with-rt=SPHM1RT_N`` where ``N`` is the integer number of photon groups that 
+    you intend to use in your simulation. The number of photon groups for SPHM1RT has
+    to be four for now (since it is hard-coded in thermo-chemistry solver).
+
+-   SPHM1-RT is compatible with any SPH scheme. You'll
+    need to compile using ``--with-hydro=sphenix`` or other SPH schemes, e.g. we have tested gadget2, minimal, and sphenix.
+
+-   SPHM1-RT solves non-equilibrium with the `SUNDIALS <https://computing.llnl.gov/projects/sundials>` library, 
+    which is SUite of Nonlinear and DIfferential/ALgebraic Equation Solvers. The SUNDIALS version has to be  5 . 
+    You'll need to compile using ``--with-sundials=$SUNDIALS_ROOT``    
+    SUNDIALS_ROOT is the root directory that contains the lib and include directories, e.g. on cosma:
+    SUNDIALS_ROOT=/cosma/local/sundials/5.1.0/
+    (IMPORTANT: Need SUNDIALS version  = 5). 
+    The instructions of installing Sundials can be found, e.g., 
+    `here <https://sundials.readthedocs.io/en/latest/Install_link.html>` or in this `website 
+    <https://richings.bitbucket.io/chimes/user_guide/GettingStarted/sundials.html>`.
+
+
+
+Runtime Parameters
+~~~~~~~~~~~~~~~~~~
+
+You need to provide the following runtime parameters in the yaml file:
+
+.. code:: yaml
+
+    SPHM1RT:
+        cred: 2.99792458e10                                 # value of reduced speed of light for the RT solver in code unit
+        CFL_condition: 0.1                                  # CFL condition for RT, independent of hydro 
+        chi:  [0, 0, 0]                                     # (Optional) initial opacity in code unit for all gas particles
+        photon_groups_Hz: [3.288e15, 5.945e15, 13.157e15]   # Photon frequency group bin edges in Hz.
+        use_const_emission_rates: 1                         # (Optional) use constant emission rates for stars as defined with star_emission_rates_erg_LSol parameter
+        star_emission_rates: [1e-32, 1e-32, 1e-32]          # (Optional) constant star emission rates (internal unit: energy/time) for each photon frequency group to use if use_constant_emission_rates is set.
+        stellar_spectrum_type: 0                            # Which radiation spectrum to use. 0: constant from 0 until some max frequency set by stellar_spectrum_const_max_frequency_Hz. 1: blackbody spectrum.
+        stellar_spectrum_const_max_frequency_Hz: 1.e17      # (Conditional) if stellar_spectrum_type=0, use this maximal frequency for the constant spectrum. 
+        stars_max_timestep: -1.                             # (Optional) restrict the maximal timestep of stars to this value (in internal units). Set to negative to turn off.
+        reinject:               1                           # (Optional) gather energy around injection radius and re-inject the energy
+
+
+The ``photon_groups_Hz`` need to be ``N - 1`` frequency edges (floats) to separate 
+the spectrum into ``N`` groups. The outer limits of zero and infinity are 
+assumed.
+
+At the moment, the only way to define star emission rates is to use constant
+star emission rates that need to be provided in the parameter file. The star 
+emission rates need to be defined for each photon frequency group individually.
+The first entry of the array is for the photon group with frequency 
+``[0, <first entry of photon_groups_Hz>)``. Each star particle will then emit
+the given energies, independent of their other properties.
+
+Furthermore, even though the parameter ``use_const_emission_rates`` is 
+intended to be optional in the future, **for now it needs to be set to 1.**, and
+it requires you to manually set the stellar emission rates via the
+``star_emission_rates`` parameter.
+
+When solving the thermochemistry, we need to assume some form of stellar
+spectrum so we may integrate over frequency bins to obtain average interaction
+rates. The parameter ``stellar_spectrum_type`` is hence required, and allows you
+to select between:
+
+- constant spectrum (``stellar_spectrum_type: 0``)
+    - This choice additionally requires you to provide a maximal frequency for
+      the spectrum after which it'll be cut off via the 
+      ``stellar_spectrum_const_max_frequency_Hz`` parameter
+
+- blackbody spectrum (``stellar_spectrum_type: 1``)
+    - In this case, you need to provide also temperature of the blackbody via the 
+      ``stellar_spectrum_blackbody_temperature_K`` parameter.
+
+If ``reinject`` is set to 1, then the star will correct the radiation energy and 
+flux within the injection radius and re-inject radiation. Combined with larger 
+smoothing length, this can increase the isotropy of the radiation and the central 
+radiation distribution. 
+
+
+
+Thermo-chemistry parameters for RT
+``````````````````````````````````
+.. code:: yaml
+
+    relativeTolerance:          1e-3                    # (Optional) Relative tolerance for SPHM1RT thermo-chemistry intergration
+    absoluteTolerance:          1e-10                   # (Optional) Absolute tolerance for SPHM1RT thermo-chemistry integration
+    explicitTolerance:          0.1                     # (Optional) Tolerance below which we use the explicit solution in SPHM1RT thermo-chemistry
+    ionizing_photon_energy_erg: [3.0208e-11, 5.61973e-11, 1.05154e-10]  # (Optional) ionizing photon energy in erg averaged over frequency bins
+    skip_thermochemistry: 0                             # (Optional) skip the thermochemistry. This is intended only for debugging and testing the radiation transport, as it breaks the purpose of RT.
+    coolingon:              1                           # (Optional) switch for cooling (and photoheating), but photo-ionization will be ongoing even if coolingon==0 
+    useparams:              1                           # (Optional) switch to use thermo-chemistry parameters from the parameter file
+    sigma_cross:            [8.13e-18, 1e-32, 1e-32] # (Conditional) (if useparams=1) The cross section of ionizing photons for hydrogen (cm^2)
+    alphaB:                 2.59e-13                    # (Conditional) (if useparams=1) The case B recombination coefficient for hydrogen (cgs)
+    beta:                   3.1e-16                   # (Conditional) (if useparams=1) The collisional ionization coefficient for hydrogen (cgs)
+
+
+ionizing_photon_energy_erg is the photon energy averaged over frequency within a frequency bin, given the radiation spectrum. In the default case, 
+the first value corresponds to the bin from HI ionizing frequency to HeI ionizing frequency.
+The second value is from HeI ionizing frequency to HeII ionizing frequency.
+The third value is above HeII ionizing frequency.
+The default values are calculated with T=1e5 K blackbody spectrum.
+We currently assume the spectral shape is unchanged and universal, i.e. T=1e5 K blackbody spectrum everywhere. 
+ionizing_photon_energy_erg is used to convert photon energy density to photon number density, photo-heating, and photo-ionization.
+
+
+
+sigma_cross is also cross-section averaged within a frequency bin. 
+
+Currently, SPHM1RT uses CVODE in SUNDIALS to solve non-equilibrium hydrogen and helium thermochemistry in three frequency bins,
+from HI-HeII, HeII-HeIII and HeIII-inf. The precise coefficients will be published in Chan et al. in prep.,
+but they can be found in src/rt_cooling_rates.h
+
+Note that the first parameter in the thermo-chemistry array 
+corresponds to the second parameter in injection array. For example, if
+star_emission_rates: [0.0, 1.0, 0.0, 0.0], 
+the star emits in the HI-HeII frequency and interacts with the first bin (8.13e-18):
+sigma_cross:            [8.13e-18, 0.0, 0.0]
+
+relativeTolerance, absoluteTolerance, and explicitTolerance are tolerances used in the CVODE calculation. 
+These tolerances can be relaxed to increase the calculation speed, which could sacrifice accuracy.
+
+We can also turn off thermochemistry or cooling for testing purpose by skip_thermochemistry and coolingon.
+For testing purpose, we can also overwrite the thermo-chemistry parameters by setting useparams to 1
+Currently, useparams==1 only works for pure hydrogen gas.
+
+
+
+
+Initial Conditions
+~~~~~~~~~~~~~~~~~~
+
+
+Setting Up Initial Conditions for RT
+````````````````````````````````````
+
+Optionally, you may want to provide initial conditions for the radiation field
+and/or the mass fraction of the ionizing species.
+To do so, you need to add the following datasets to the ``/PartType0`` particle
+group:
+
+.. code:: 
+
+   PhotonEnergiesGroup1
+   PhotonEnergiesGroup2 
+   .
+   .
+   .
+   PhotonEnergiesGroupN
+   PhotonFluxesGroup1
+   PhotonFluxesGroup2
+   .
+   .
+   .
+   PhotonFluxesGroupN
+
+
+The ``PhotonEnergies*`` datasets need to have dimension ``nparts``, while the
+``PhotonFluxesGroup*`` datasets need to have dimension ``(nparts, 3)``, where
+``nparts`` is the number of hydro particles. If you are writing initial
+conditions where the fields have units, then ``PhotonEnergies*`` are expected to
+have units of energy :math:`[M L^2 T^{-2}]`), while the ``PhotonFluxes*`` fields
+should be in units of energy times speed, :math:`[M L^3
+T^{-3}]`).
+
+
+Example using Python and ``swiftsimio``
+````````````````````````````````````````
+
+If you are using `swiftsimio <https://github.com/SWIFTSIM/swiftsimio>`_ to write
+the initial condition files, then the easiest way of adding the RT initial
+conditions is to first use the swiftsimio routines to write a file, then open it
+up again and write the additional RT fields again using ``h5py`` routines.
+
+Here is an example:
+
+.. code:: python
+
+    from swiftsimio import Writer
+    import unyt
+    import numpy as np
+    import h5py
+
+    # define unit system to use.
+    unitsystem = unyt.unit_systems.cgs_unit_system
+
+    # number of photon groups
+    nPhotonGroups = 4
+
+    # filename of ICs to be generated
+    outputfilename = "my_rt_ICs.hdf5"
+
+    # open a swiftsimio.Writer object
+    w = Writer(...)
+
+    # do your IC setup for gas, gravity etc now
+    # ... 
+
+    # write the IC file without doing anything RT related.
+    w.write(outputfilename)
+
+    # Now open file back up again and add RT data.
+    F = h5py.File(outputfilename, "r+")
+    header = F["Header"]
+    nparts = header.attrs["NumPart_ThisFile"][0]
+    parts = F["/PartType0"]
+
+    # Create initial photon energies and fluxes. You can leave them unitless, 
+    # the units have already been written down with w.write(). In this case, 
+    # it's in cgs.
+    for grp in range(nPhotonGroups):
+        dsetname = "PhotonEnergiesGroup{0:d}".format(grp + 1)
+        energydata = np.ones((nparts), dtype=np.float32) * some_value_you_want
+        parts.create_dataset(dsetname, data=energydata)
+
+        dsetname = "PhotonFluxesGroup{0:d}".format(grp + 1)
+        fluxdata = np.zeros((nparts, 3), dtype=np.float32) * some_value_you_want
+        parts.create_dataset(dsetname, data=fluxdata)
+
+    # Create initial element mass fractions.
+    # Can be overwritten in parameter file if init_mass_fraction_metal is not -1.f (or set)
+    # the element order: [Hydrogen, Helium]
+    mfracH = np.ones(numPart)
+    mfracHe = np.ones(numPart) * 0.0
+    EMFdata = np.stack((mfracH, mfracHe), axis=-1)
+    parts.create_dataset("RtElementMassFractions", data=EMFdata)
+
+    # Create initial species abundances.
+    # abundance is in n_X/n_H unit.
+    # Can be overwritten in parameter file if useabundances = 1
+    # the abundance order: [e, HI, HII, HeI, HeII, HeIII]
+    Ae = np.ones(numPart) * 0.0   
+    AHI = np.ones(numPart) * 1.0  
+    AHII = np.ones(numPart) * 0.0 
+    AHeI = np.ones(numPart) * 0.0 
+    AHeII = np.ones(numPart) * 0.0 
+    AHeIII = np.ones(numPart) * 0.0 
+    SAdata = np.stack((Ae, AHI, AHII, AHeI, AHeII, AHeIII), axis=-1)    
+    parts.create_dataset("RtSpeciesAbundances", data=SAdata)
+
+    # close up, and we're done!
+    F.close()
+
+
+
+Generate Ionization Mass Fractions Using SWIFT
+``````````````````````````````````````````````
+
+.. warning:: Using SWIFT to generate initial ionization mass fractions will
+   overwrite the mass fractions that have been read in from the initial 
+   conditions.
+
+Optionally, you can use SWIFT to generate the initial mass fractions of the
+elements. To set the initial mass fractions of all particles to the same
+value, use the following parameters in the yaml parameter file:
+
+.. code:: yaml
+
+  init_mass_fraction_metal:     0.                    # (Optional) Inital mass fraction of particle mass in *all* metals (if it is set or not equal to -1.F, the initial fraction will be over-written.)
+  init_mass_fraction_Hydrogen:  1.0                   # (Conditional) (if init_mass_fraction_metal != -1.0f) Inital mass fraction of particle mass in Hydrogen
+  init_mass_fraction_Helium:    0.0                   # (Conditional) (if init_mass_fraction_metal != -1.0f) Inital mass fraction of particle mass in Helium
+
+To set the species abundances of all particles to the same
+value, use the following parameters in the yaml parameter file:
+
+.. code:: yaml
+
+  useabundances:              1                       # (Optional) use the species abundances below, instead of reading from initial condition
+  init_species_abundance_e:        1e-5               # (Conditional) (if useabundances==1) free electron abundances (in unit hydrogen number density:nH)
+  init_species_abundance_HI:       0.99999            # (Conditional) (if useabundances==1) HI abundances (in unit hydrogen number density:nH)
+  init_species_abundance_HII:      1e-5               # (Conditional) (if useabundances==1) HII abundances (in unit hydrogen number density:nH)
+  init_species_abundance_HeI:      0.0                # (Conditional) (if useabundances==1) HeI abundances (in unit hydrogen number density:nH)
+  init_species_abundance_HeII:     0.0                # (Conditional) (if useabundances==1) HeII abundances (in unit hydrogen number density:nH)
+  init_species_abundance_HeIII:    0.0                # (Conditional) (if useabundances==1) HeIII abundances (in unit hydrogen number density:nH)
+
+
+Accessing Output Data
+~~~~~~~~~~~~~~~~~~~~~~
+
+We recommend using `swiftsimio <https://github.com/SWIFTSIM/swiftsimio>`_ to 
+access the RT related snapshot data. The compatibility is being maintained.
+Here's an example how to access some specific quantities that you might find
+useful:
+
+
+.. code:: python
+
+    #!/usr/bin/env python3
+
+    import swiftsimio
+    import unyt
+
+    data = swiftsimio.load("output_0001.hdf5")
+    meta = data.metadata
+
+
+
+    # Accessing RT Related Metadata
+    # ---------------------------------
+
+    # get scheme name: "SPH M1closure"
+    scheme = str(meta.subgrid_scheme["RT Scheme"].decode("utf-8"))
+
+    # number of photon groups used
+    ngroups = int(meta.subgrid_scheme["PhotonGroupNumber"])
+
+    # get the reduced speed of light that was used. Will have unyts.
+    reduced_speed_of_light = meta.reduced_lightspeed
+
+
+
+
+    # Accessing Photon Data
+    # ------------------------
+
+    # accessing a photon group directly
+    # NOTE: group names start with 1
+    group_1_photon_energies = data.gas.photon_energies.group1
+    group_1_photon_fluxes_x = data.gas.photon_fluxes.Group1X
+    group_1_photon_fluxes_y = data.gas.photon_fluxes.Group1Y
+    group_1_photon_fluxes_z = data.gas.photon_fluxes.Group1Z
+
+    # want to stack all fluxes into 1 array?
+    group1fluxes = swiftsimio.cosmo_array(
+        unyt.uvstack(
+            (group_1_photon_fluxes_x, group_1_photon_fluxes_y, group_1_photon_fluxes_z)
+        ),
+        group_1_photon_fluxes_x.units,
+    ).T
+    # group1fluxes.shape = (npart, 3)
+
+
+    # Load all photon energies in a list
+    photon_energies = [
+        getattr(data.gas.photon_energies, "group" + str(g + 1)) for g in range(ngroups)
+    ]
+
+
+    # Accessing Element mass fraction
+    fH = data.gas.rt_element_mass_fractions.hydrogen
+    fHe = data.gas.rt_element_mass_fractions.helium
+
+    # Accessing Species Abundances 
+    # abundance is in n_X/n_H unit.
+    # -------------------------------
+    Ae = data.gas.rt_species_abundances.e
+    AHI = data.gas.rt_species_abundances.HI
+    AHII = data.gas.rt_species_abundances.HII
+    AHeI = data.gas.rt_species_abundances.HeI
+    AHeII = data.gas.rt_species_abundances.HeII
+    AHeIII = data.gas.rt_species_abundances.HeIII

doc/RTD/source/RadiativeTransfer/additional_tools.rst

+.. GEAR Radiative Transfer
+    Mladen Ivkovic 08.2022
+
+.. _rt_additional_tools:
+
+
+Additional Tools
+-------------------
+
+A plethora of additional RT related tools is hosted in the 
+`the swiftsim-rt-tools <https://github.com/SWIFTSIM/swiftsim-rt-tools>`_
+repository. For example, it contains tools to convert injected photon number 
+rates into luminosities, compute integrated interaction cross sections for 
+various spectra, some rudimentary thermochemistry python functions, etc.
+

doc/RTD/source/RadiativeTransfer/full_dependency_graph_RT.dot

+digraph task_dep {
+	 # Header
+	 label="Task dependencies for SWIFT v0.9.0-1189-g3424fa24";
+	 compound=true;
+	 ratio=0.66;
+	 node[nodesep=0.15, fontsize=18, penwidth=3.];
+	 edge[fontsize=12, penwidth=0.5];
+	 ranksep=0.8;
+
+	 # Special tasks
+	 sort[color=blue3];
+	 self_density[color=blue3];
+	 self_gradient[color=blue3];
+	 self_force[color=blue3];
+	 self_grav[color=red3];
+	 self_stars_density[color=darkorange1];
+	 self_stars_feedback[color=darkorange1];
+	 self_rt_gradient[color=springgreen];
+	 self_rt_transport[color=springgreen];
+	 pair_density[color=blue3];
+	 pair_gradient[color=blue3];
+	 pair_force[color=blue3];
+	 pair_grav[color=red3];
+	 pair_stars_density[color=darkorange1];
+	 pair_stars_feedback[color=darkorange1];
+	 pair_rt_gradient[color=springgreen];
+	 pair_rt_transport[color=springgreen];
+	 sub_self_density[color=blue3];
+	 sub_self_gradient[color=blue3];
+	 sub_self_force[color=blue3];
+	 sub_self_stars_density[color=darkorange1];
+	 sub_self_stars_feedback[color=darkorange1];
+	 sub_self_rt_gradient[color=springgreen];
+	 sub_self_rt_transport[color=springgreen];
+	 sub_pair_density[color=blue3];
+	 sub_pair_gradient[color=blue3];
+	 sub_pair_force[color=blue3];
+	 sub_pair_stars_density[color=darkorange1];
+	 sub_pair_stars_feedback[color=darkorange1];
+	 sub_pair_rt_gradient[color=springgreen];
+	 sub_pair_rt_transport[color=springgreen];
+	 init_grav[color=red3];
+	 init_grav_out[style=filled,fillcolor=grey90,color=red3];
+	 ghost_in[style=filled,fillcolor=grey90,color=blue3];
+	 ghost[color=blue3];
+	 ghost_out[style=filled,fillcolor=grey90,color=blue3];
+	 extra_ghost[color=blue3];
+	 drift_part[color=blue3];
+	 drift_spart[color=darkorange1];
+	 drift_gpart[color=red3];
+	 drift_gpart_out[style=filled,fillcolor=grey90,color=red3];
+	 hydro_end_force[color=blue3];
+	 kick2[color=black];
+	 timestep[color=black];
+	 collect[color=black];
+	 send_gradient[shape=diamond,style=filled,fillcolor=azure,color=blue3];
+	 send_xv[shape=diamond,style=filled,fillcolor=azure,color=blue3];
+	 send_rho[shape=diamond,style=filled,fillcolor=azure,color=blue3];
+	 send_gpart[shape=diamond,style=filled,fillcolor=azure,color=red3];
+	 send_spart_density[shape=diamond,style=filled,fillcolor=azure,color=darkorange1];
+	 send_rt_gradient[shape=diamond,style=filled,fillcolor=azure,color=springgreen];
+	 send_rt_transport[shape=diamond,style=filled,fillcolor=azure,color=springgreen];
+	 recv_gradient[shape=diamond,style=filled,fillcolor=azure,color=blue3];
+	 recv_xv[shape=diamond,style=filled,fillcolor=azure,color=blue3];
+	 recv_rho[shape=diamond,style=filled,fillcolor=azure,color=blue3];
+	 recv_gpart[shape=diamond,style=filled,fillcolor=azure,color=red3];
+	 recv_spart_density[shape=diamond,style=filled,fillcolor=azure,color=darkorange1];
+	 recv_rt_gradient[shape=diamond,style=filled,fillcolor=azure,color=springgreen];
+	 recv_rt_transport[shape=diamond,style=filled,fillcolor=azure,color=springgreen];
+	 grav_long_range[color=red3];
+	 grav_mm[color=red3];
+	 grav_down_in[style=filled,fillcolor=grey90,color=red3];
+	 grav_down[color=red3];
+	 grav_end_force[color=red3];
+	 stars_in[style=filled,fillcolor=grey90,color=darkorange1];
+	 stars_out[style=filled,fillcolor=grey90,color=darkorange1];
+	 stars_density_ghost[color=darkorange1];
+	 stars_sort[color=darkorange1];
+	 rt_in[style=filled,fillcolor=grey90,color=springgreen];
+	 rt_out[style=filled,fillcolor=grey90,color=springgreen];
+	 rt_ghost1[color=springgreen];
+	 rt_ghost2[color=springgreen];
+	 rt_transport_out[style=filled,fillcolor=grey90,color=springgreen];
+	 rt_tchem[color=springgreen];
+	 rt_advance_cell_time[color=springgreen];
+	 rt_sorts[color=springgreen];
+	 recv_tend[shape=diamond,style=filled,fillcolor=azure,color=black];
+	 kick1[color=black];
+	 send_tend[shape=diamond,style=filled,fillcolor=azure,color=black];
+	 rt_collect_times[color=springgreen];
+
+	 # Clusters
+	 subgraph clusterDensity {
+		 label="";
+		 bgcolor="grey99";
+		 pair_density;
+		 self_density;
+		 sub_pair_density;
+		 sub_self_density;
+	 };
+
+	 subgraph clusterForce {
+		 label="";
+		 bgcolor="grey99";
+		 pair_force;
+		 self_force;
+		 sub_pair_force;
+		 sub_self_force;
+	 };
+
+	 subgraph clusterGradient {
+		 label="";
+		 bgcolor="grey99";
+		 pair_gradient;
+		 self_gradient;
+		 sub_pair_gradient;
+		 sub_self_gradient;
+	 };
+
+	 subgraph clusterGravity {
+		 label="";
+		 bgcolor="grey99";
+		 grav_long_range;
+		 grav_mm;
+		 pair_grav;
+		 self_grav;
+	 };
+
+	 subgraph clusterRTgradient {
+		 label="";
+		 bgcolor="grey99";
+		 pair_rt_gradient;
+		 self_rt_gradient;
+		 sub_pair_rt_gradient;
+		 sub_self_rt_gradient;
+	 };
+
+	 subgraph clusterRTtransport {
+		 label="";
+		 bgcolor="grey99";
+		 pair_rt_transport;
+		 self_rt_transport;
+		 sub_pair_rt_transport;
+		 sub_self_rt_transport;
+	 };
+
+	 subgraph clusterStarsDensity {
+		 label="";
+		 bgcolor="grey99";
+		 pair_stars_density;
+		 self_stars_density;
+		 sub_pair_stars_density;
+		 sub_self_stars_density;
+	 };
+
+	 subgraph clusterStarsFeedback {
+		 label="";
+		 bgcolor="grey99";
+		 pair_stars_feedback;
+		 self_stars_feedback;
+		 sub_pair_stars_feedback;
+		 sub_self_stars_feedback;
+	 };
+
+
+	 # Dependencies
+	 sort->pair_density[color=blue3,fontcolor=blue3]
+	 sort->pair_force[color=blue3,fontcolor=blue3]
+	 sort->pair_rt_transport[color=blue3,fontcolor=blue3]
+	 sort->pair_stars_density[color=blue3,fontcolor=blue3]
+	 sort->pair_rt_gradient[color=blue3,fontcolor=blue3]
+	 sort->sub_pair_density[color=blue3,fontcolor=blue3]
+	 sort->sub_pair_force[color=blue3,fontcolor=blue3]
+	 sort->sub_pair_rt_transport[color=blue3,fontcolor=blue3]
+	 sort->sub_pair_stars_density[color=blue3,fontcolor=blue3]
+	 sort->sub_pair_rt_gradient[color=blue3,fontcolor=blue3]
+	 sort->rt_sorts[color=blue3,fontcolor=blue3]
+	 sort->recv_rt_gradient[color=blue3,fontcolor=blue3]
+	 sort->recv_rho[color=blue3,fontcolor=blue3]
+	 sort->recv_gradient[color=blue3,fontcolor=blue3]
+	 sort->sub_self_density[color=blue3,fontcolor=blue3]
+	 sort->sub_self_stars_density[color=blue3,fontcolor=blue3]
+	 sort->sub_self_rt_gradient[color=blue3,fontcolor=blue3]
+	 self_density->ghost_in[color=blue3,fontcolor=blue3]
+	 self_gradient->extra_ghost[color=blue3,fontcolor=blue3]
+	 self_force->hydro_end_force[color=blue3,fontcolor=blue3]
+	 self_grav->grav_down_in[color=red3,fontcolor=red3]
+	 self_stars_density->stars_density_ghost[color=darkorange1,fontcolor=darkorange1]
+	 self_stars_feedback->stars_out[color=darkorange1,fontcolor=darkorange1]
+	 self_rt_gradient->rt_ghost2[color=springgreen,fontcolor=springgreen]
+	 self_rt_transport->rt_transport_out[color=springgreen,fontcolor=springgreen]
+	 pair_density->ghost_in[color=blue3,fontcolor=blue3]
+	 pair_density->recv_rho[color=blue3,fontcolor=blue3]
+	 pair_gradient->extra_ghost[color=blue3,fontcolor=blue3]
+	 pair_gradient->recv_gradient[color=blue3,fontcolor=blue3]
+	 pair_force->hydro_end_force[color=blue3,fontcolor=blue3]
+	 pair_force->recv_tend[color=blue3,fontcolor=blue3]
+	 pair_force->recv_rt_gradient[color=blue3,fontcolor=blue3]
+	 pair_grav->grav_down_in[color=red3,fontcolor=red3]
+	 pair_grav->recv_tend[color=red3,fontcolor=red3]
+	 pair_stars_density->stars_density_ghost[color=darkorange1,fontcolor=darkorange1]
+	 pair_stars_density->recv_spart_density[color=darkorange1,fontcolor=darkorange1]
+	 pair_stars_feedback->stars_out[color=darkorange1,fontcolor=darkorange1]
+	 pair_stars_feedback->recv_tend[color=darkorange1,fontcolor=darkorange1]
+	 pair_rt_gradient->rt_ghost2[color=springgreen,fontcolor=springgreen]
+	 pair_rt_gradient->recv_rt_transport[color=springgreen,fontcolor=springgreen]
+	 pair_rt_transport->rt_transport_out[color=springgreen,fontcolor=springgreen]
+	 pair_rt_transport->recv_tend[color=springgreen,fontcolor=springgreen]
+	 pair_rt_transport->rt_advance_cell_time[color=springgreen,fontcolor=springgreen]
+	 sub_self_density->ghost_in[color=blue3,fontcolor=blue3]
+	 sub_self_gradient->extra_ghost[color=blue3,fontcolor=blue3]
+	 sub_self_force->hydro_end_force[color=blue3,fontcolor=blue3]
+	 sub_self_stars_density->stars_density_ghost[color=darkorange1,fontcolor=darkorange1]
+	 sub_self_stars_feedback->stars_out[color=darkorange1,fontcolor=darkorange1]
+	 sub_self_rt_gradient->rt_ghost2[color=springgreen,fontcolor=springgreen]
+	 sub_self_rt_transport->rt_transport_out[color=springgreen,fontcolor=springgreen]
+	 sub_pair_density->ghost_in[color=blue3,fontcolor=blue3]
+	 sub_pair_density->recv_rho[color=blue3,fontcolor=blue3]
+	 sub_pair_gradient->extra_ghost[color=blue3,fontcolor=blue3]
+	 sub_pair_gradient->recv_gradient[color=blue3,fontcolor=blue3]
+	 sub_pair_force->hydro_end_force[color=blue3,fontcolor=blue3]
+	 sub_pair_force->recv_tend[color=blue3,fontcolor=blue3]
+	 sub_pair_force->recv_rt_gradient[color=blue3,fontcolor=blue3]
+	 sub_pair_stars_density->stars_density_ghost[color=darkorange1,fontcolor=darkorange1]
+	 sub_pair_stars_density->recv_spart_density[color=darkorange1,fontcolor=darkorange1]
+	 sub_pair_stars_feedback->stars_out[color=darkorange1,fontcolor=darkorange1]
+	 sub_pair_stars_feedback->recv_tend[color=darkorange1,fontcolor=darkorange1]
+	 sub_pair_rt_gradient->rt_ghost2[color=springgreen,fontcolor=springgreen]
+	 sub_pair_rt_gradient->recv_rt_transport[color=springgreen,fontcolor=springgreen]
+	 sub_pair_rt_transport->rt_transport_out[color=springgreen,fontcolor=springgreen]
+	 sub_pair_rt_transport->recv_tend[color=springgreen,fontcolor=springgreen]
+	 sub_pair_rt_transport->rt_advance_cell_time[color=springgreen,fontcolor=springgreen]
+	 init_grav->grav_long_range[color=red3,fontcolor=red3]
+	 init_grav->init_grav_out[color=red3,fontcolor=red3]
+	 init_grav_out->pair_grav[color=red3,fontcolor=red3]
+	 init_grav_out->self_grav[color=red3,fontcolor=red3]
+	 init_grav_out->init_grav_out[color=red3,fontcolor=red3]
+	 init_grav_out->grav_mm[color=red3,fontcolor=red3]
+	 ghost_in->ghost[color=blue3,fontcolor=blue3]
+	 ghost->ghost_out[color=blue3,fontcolor=blue3]
+	 ghost_out->pair_gradient[color=blue3,fontcolor=blue3]
+	 ghost_out->sub_self_gradient[color=blue3,fontcolor=blue3]
+	 ghost_out->sub_pair_gradient[color=blue3,fontcolor=blue3]
+	 ghost_out->send_rho[color=blue3,fontcolor=blue3]
+	 ghost_out->self_gradient[color=blue3,fontcolor=blue3]
+	 extra_ghost->pair_force[color=blue3,fontcolor=blue3]
+	 extra_ghost->sub_self_force[color=blue3,fontcolor=blue3]
+	 extra_ghost->sub_pair_force[color=blue3,fontcolor=blue3]
+	 extra_ghost->send_gradient[color=blue3,fontcolor=blue3]
+	 extra_ghost->self_force[color=blue3,fontcolor=blue3]
+	 drift_part->pair_density[color=blue3,fontcolor=blue3]
+	 drift_part->pair_stars_density[color=blue3,fontcolor=blue3]
+	 drift_part->pair_rt_gradient[color=blue3,fontcolor=blue3]
+	 drift_part->sub_self_density[color=blue3,fontcolor=blue3]
+	 drift_part->sub_self_stars_density[color=blue3,fontcolor=blue3]
+	 drift_part->sub_self_rt_gradient[color=blue3,fontcolor=blue3]
+	 drift_part->sub_pair_density[color=blue3,fontcolor=blue3]
+	 drift_part->sub_pair_stars_density[color=blue3,fontcolor=blue3]
+	 drift_part->sub_pair_rt_gradient[color=blue3,fontcolor=blue3]
+	 drift_part->sort[color=blue3,fontcolor=blue3]
+	 drift_part->send_rho[color=blue3,fontcolor=blue3]
+	 drift_part->send_xv[color=blue3,fontcolor=blue3]
+	 drift_part->send_rt_gradient[color=blue3,fontcolor=blue3]
+	 drift_part->send_rt_transport[color=blue3,fontcolor=blue3]
+	 drift_part->self_density[color=blue3,fontcolor=blue3]
+	 drift_part->self_stars_density[color=blue3,fontcolor=blue3]
+	 drift_part->self_rt_gradient[color=blue3,fontcolor=blue3]
+	 drift_spart->kick2[color=darkorange1,fontcolor=darkorange1]
+	 drift_spart->pair_stars_density[color=darkorange1,fontcolor=darkorange1]
+	 drift_spart->sub_self_stars_density[color=darkorange1,fontcolor=darkorange1]
+	 drift_spart->sub_pair_stars_density[color=darkorange1,fontcolor=darkorange1]
+	 drift_spart->stars_sort[color=darkorange1,fontcolor=darkorange1]
+	 drift_spart->send_spart_density[color=darkorange1,fontcolor=darkorange1]
+	 drift_spart->self_stars_density[color=darkorange1,fontcolor=darkorange1]
+	 drift_gpart->drift_gpart_out[color=red3,fontcolor=red3]
+	 drift_gpart->send_gpart[color=red3,fontcolor=red3]
+	 drift_gpart_out->pair_grav[color=red3,fontcolor=red3]
+	 drift_gpart_out->self_grav[color=red3,fontcolor=red3]
+	 drift_gpart_out->drift_gpart_out[color=red3,fontcolor=red3]
+	 hydro_end_force->kick2[color=blue3,fontcolor=blue3]
+	 kick2->timestep[color=black,fontcolor=black]
+	 kick2->stars_in[color=black,fontcolor=black]
+	 kick2->rt_in[color=black,fontcolor=black]
+	 timestep->kick1[color=black,fontcolor=black]
+	 timestep->collect[color=black,fontcolor=black]
+	 collect->send_tend[color=black,fontcolor=black]
+	 send_gradient->hydro_end_force[color=blue3,fontcolor=blue3]
+	 send_xv->send_rho[color=blue3,fontcolor=blue3]
+	 send_xv->ghost_in[color=blue3,fontcolor=blue3]
+	 send_xv->send_rt_gradient[color=blue3,fontcolor=blue3]
+	 send_rho->send_gradient[color=blue3,fontcolor=blue3]
+	 send_rho->extra_ghost[color=blue3,fontcolor=blue3]
+	 send_gpart->grav_down[color=red3,fontcolor=red3]
+	 send_spart_density->stars_out[color=darkorange1,fontcolor=darkorange1]
+	 send_rt_gradient->send_rt_transport[color=springgreen,fontcolor=springgreen]
+	 send_rt_gradient->rt_ghost2[color=springgreen,fontcolor=springgreen]
+	 send_rt_transport->rt_transport_out[color=springgreen,fontcolor=springgreen]
+	 recv_gradient->recv_rt_gradient[color=blue3,fontcolor=blue3]
+	 recv_gradient->pair_force[color=blue3,fontcolor=blue3]
+	 recv_gradient->sub_pair_force[color=blue3,fontcolor=blue3]
+	 recv_xv->recv_rho[color=blue3,fontcolor=blue3]
+	 recv_xv->recv_gradient[color=blue3,fontcolor=blue3]
+	 recv_xv->recv_rt_gradient[color=blue3,fontcolor=blue3]
+	 recv_xv->sort[color=blue3,fontcolor=blue3]
+	 recv_xv->pair_density[color=blue3,fontcolor=blue3]
+	 recv_xv->sub_pair_density[color=blue3,fontcolor=blue3]
+	 recv_rho->recv_gradient[color=blue3,fontcolor=blue3]
+	 recv_rho->recv_rt_gradient[color=blue3,fontcolor=blue3]
+	 recv_rho->pair_gradient[color=blue3,fontcolor=blue3]
+	 recv_rho->pair_stars_density[color=blue3,fontcolor=blue3]
+	 recv_rho->sub_pair_gradient[color=blue3,fontcolor=blue3]
+	 recv_rho->sub_pair_stars_density[color=blue3,fontcolor=blue3]
+	 recv_gpart->pair_grav[color=red3,fontcolor=red3]
+	 recv_spart_density->stars_sort[color=darkorange1,fontcolor=darkorange1]
+	 recv_spart_density->pair_stars_feedback[color=darkorange1,fontcolor=darkorange1]
+	 recv_spart_density->sub_pair_stars_feedback[color=darkorange1,fontcolor=darkorange1]
+	 recv_rt_gradient->rt_sorts[color=springgreen,fontcolor=springgreen]
+	 recv_rt_gradient->recv_rt_transport[color=springgreen,fontcolor=springgreen]
+	 recv_rt_gradient->rt_advance_cell_time[color=springgreen,fontcolor=springgreen]
+	 recv_rt_gradient->pair_rt_gradient[color=springgreen,fontcolor=springgreen]
+	 recv_rt_gradient->sub_pair_rt_gradient[color=springgreen,fontcolor=springgreen]
+	 recv_rt_transport->rt_advance_cell_time[color=springgreen,fontcolor=springgreen]
+	 recv_rt_transport->pair_rt_transport[color=springgreen,fontcolor=springgreen]
+	 recv_rt_transport->sub_pair_rt_transport[color=springgreen,fontcolor=springgreen]
+	 grav_long_range->grav_down[color=red3,fontcolor=red3]
+	 grav_mm->grav_down_in[color=red3,fontcolor=red3]
+	 grav_down_in->grav_down[color=red3,fontcolor=red3]
+	 grav_down_in->grav_down_in[color=red3,fontcolor=red3]
+	 grav_down->grav_end_force[color=red3,fontcolor=red3]
+	 grav_end_force->kick2[color=red3,fontcolor=red3]
+	 stars_in->pair_stars_density[color=darkorange1,fontcolor=darkorange1]
+	 stars_in->sub_self_stars_density[color=darkorange1,fontcolor=darkorange1]
+	 stars_in->sub_pair_stars_density[color=darkorange1,fontcolor=darkorange1]
+	 stars_in->self_stars_density[color=darkorange1,fontcolor=darkorange1]
+	 stars_out->timestep[color=darkorange1,fontcolor=darkorange1]
+	 stars_out->rt_in[color=darkorange1,fontcolor=darkorange1]
+	 stars_density_ghost->pair_stars_feedback[color=darkorange1,fontcolor=darkorange1]
+	 stars_density_ghost->sub_self_stars_feedback[color=darkorange1,fontcolor=darkorange1]
+	 stars_density_ghost->sub_pair_stars_feedback[color=darkorange1,fontcolor=darkorange1]
+	 stars_density_ghost->send_spart_density[color=darkorange1,fontcolor=darkorange1]
+	 stars_density_ghost->self_stars_feedback[color=darkorange1,fontcolor=darkorange1]
+	 stars_sort->pair_stars_feedback[color=darkorange1,fontcolor=darkorange1]
+	 stars_sort->sub_pair_stars_feedback[color=darkorange1,fontcolor=darkorange1]
+	 stars_sort->pair_stars_density[color=darkorange1,fontcolor=darkorange1]
+	 stars_sort->sub_self_stars_density[color=darkorange1,fontcolor=darkorange1]
+	 stars_sort->sub_pair_stars_density[color=darkorange1,fontcolor=darkorange1]
+	 rt_in->rt_ghost1[color=springgreen,fontcolor=springgreen]
+	 rt_out->timestep[color=springgreen,fontcolor=springgreen]
+	 rt_out->collect[color=springgreen,fontcolor=springgreen]
+	 rt_ghost1->pair_rt_gradient[color=springgreen,fontcolor=springgreen]
+	 rt_ghost1->sub_self_rt_gradient[color=springgreen,fontcolor=springgreen]
+	 rt_ghost1->sub_pair_rt_gradient[color=springgreen,fontcolor=springgreen]
+	 rt_ghost1->send_rt_gradient[color=springgreen,fontcolor=springgreen]
+	 rt_ghost1->self_rt_gradient[color=springgreen,fontcolor=springgreen]
+	 rt_ghost2->pair_rt_transport[color=springgreen,fontcolor=springgreen]
+	 rt_ghost2->sub_self_rt_transport[color=springgreen,fontcolor=springgreen]
+	 rt_ghost2->sub_pair_rt_transport[color=springgreen,fontcolor=springgreen]
+	 rt_ghost2->send_rt_transport[color=springgreen,fontcolor=springgreen]
+	 rt_ghost2->self_rt_transport[color=springgreen,fontcolor=springgreen]
+	 rt_transport_out->rt_tchem[color=springgreen,fontcolor=springgreen]
+	 rt_tchem->rt_advance_cell_time[color=springgreen,fontcolor=springgreen]
+	 rt_advance_cell_time->rt_collect_times[color=springgreen,fontcolor=springgreen]
+	 rt_advance_cell_time->rt_out[color=springgreen,fontcolor=springgreen]
+	 rt_advance_cell_time->send_tend[color=springgreen,fontcolor=springgreen]
+	 rt_advance_cell_time->recv_tend[color=springgreen,fontcolor=springgreen]
+	 rt_sorts->recv_rt_transport[color=springgreen,fontcolor=springgreen]
+	 rt_sorts->pair_rt_gradient[color=springgreen,fontcolor=springgreen]
+	 rt_sorts->sub_pair_rt_gradient[color=springgreen,fontcolor=springgreen]
+}
\ No newline at end of file

doc/RTD/source/RadiativeTransfer/index.rst

+.. Radiative Transfer Schemes
+    Mladen Ivkovic 05.2021
+
+.. _rt:
+   
+Radiative Transfer Schemes
+==========================
+
+This section of the documentation includes information on the radiative transfer
+schemes. 
+
+
+.. warning::
+    The radiative transfer schemes are still in development and are not useable
+    at this moment. This page is currently a placeholder to document new
+    features and requirements as the code grows.
+
+
+
+.. toctree::
+   :maxdepth: 2
+   :caption: Contents:
+
+   RT_general
+   GEAR_RT
+   SPHM1_RT
+   RT_notes_for_developers
+   additional_tools

doc/RTD/source/RadiativeTransfer/make_graphs.py

+import os
+
+# Yes... this is overkill... but the main sphinx conf script looks for python
+# scripts to execute... So we piggy back on this.
+os.system("dot -Tpng full_dependency_graph_RT.dot -o full_dependency_graph_RT.png")

doc/RTD/source/Snapshots/index.rst

@@ -24,24 +24,85 @@ that is always set to the string ``SWIFT``, which can be used to identify
 SWIFT-generated snapshots and hence make use of all the extensions to the file
 format described below.
 
-The most important quantity of the header is the array ``NumPart_ThisFile``
-which contains the number of particles of each type in this snapshot. This is an
-array of 6 numbers; one for each of the 5 supported types and a dummy "type 3"
-field only used for compatibility reasons but always containing a zero.
+The most important quantity of the header is the array ``NumPart_Total`` which
+contains the number of particles of each type in this snapshot. This is an array
+of 6 numbers; one for each of the supported types. Following the Gadget-2
+convention, if that number is larger than 2^31, SWIFT will use the
+``NumPart_HighWord`` field to store the high-word bits of the total number of
+particles. The field ``NumPart_ThisFile`` contains the number of particles in
+this sub-snapshot file when the user asked for distributed snapshots (see
+:ref:`Parameters_snapshots`); otherwise it contains the same information as
+``NumPart_Total``. Note however, that there is no high-word for this field. We
+store it as a 64-bits integer [#f1]_. The field ``NumFilesPerSnapshot`` specifies the
+number of sub-snapshot files (always 1 unless a distributed snapshot was asked
+for) and ``ThisFile`` the id of that specific file (always 0 unless a distributed
+snapshot was asked for). 
+
+The field ``TotalNumberOfParticles`` gives the total number of particles of each type
+as a 64 bit integer. This allows the total number of particles to be read directly
+with no calculation required even if there are 2^31 or more particles. This field is
+equal to ``NumPart_ThisFile`` if the snapshot is not distributed over multiple files.
+
+The field ``InitialMassTable`` contains the *mean* initial mass of each of the
+particle types present in the initial conditions. This can be used as estimator
+of the mass resolution of the run. The masses are expressed in internal units.
+
+The field ``OutputType`` contains information about the kind of output this
+snapshot is. The possible values are:
+
++---------------------+-----------------------------------------------------+
+| OutputType          | Definition                                          |
++=====================+=====================================================+
+| ``FullVolume``      | Regular vanilla snapshot                            |
++---------------------+-----------------------------------------------------+
+| ``SubSampled``      | Snapshot where some particle types were sub-sampled |
++---------------------+-----------------------------------------------------+
+| ``LineOfSight``     | Line-of-sight snapshot                              |
++---------------------+-----------------------------------------------------+
+| ``FOF``             | Friends-Of-Friends Halo Catalogue                   |
++---------------------+-----------------------------------------------------+
+
 
 The ``RunName`` field contains the name of the simulation that was specified as
 the ``run_name`` in the :ref:`Parameters_meta_data` section of the YAML
 parameter file.
 
+The ``System`` field contains the name of the machine where the MPI rank 0 was
+placed. This name is whatever UNIX's ``gethostname()`` function returns on that
+system. Similarly, the ``SnapshotDate`` field contains the date and time when
+the file was written.
+
 The ``TimeBase_dloga`` field contains the change in logarithm of the
 scale-factor corresponding to a time-step of length 1 on the integer
 time-line. This is the smallest time-step size that the code can use. This field
 is zero in non-cosmological runs. Similarly, the field ``TimeBase_dt`` contains
 the smallest time-step size (in internal units) that the code can take. This
 would be the increase in time a particle in the time-bin one would have. Note
-that in cosmological runs this quantity evolves with redhsift as the (logarithm
+that in cosmological runs this quantity evolves with redshift as the (logarithm
 of the) scale-factor is used on the integer time-line.
 
+The field ``SelectOutput`` will contain the name of the
+:ref:`Output_selection_label` used for this specific output and will take the value
+``Default`` if no such selection (or the default one) was used.
+
+If a sub-sampling of the particle fields was used, then the header additionally
+contains a field describing the fraction of the particles of each type that were
+written to the snapshot. Note, however, that when sub-sampling the fields 
+``NumPart_Total``, ``NumPart_HighWord``, and ``NumPart_ThisFile`` contain the number
+of particles actually written (i.e. after sub-sampling), not the total number of
+particles in the run.
+
+The field ``CanHaveTypes`` contains information about whether a given particle
+type is to be expected in snapshots of the run. For instance, a simulation with
+star formation switched on, the code may not have formed a star yet but might in
+future snapshots. This allows reading tools to distinguish fields they will
+never expect to find in a given simulation from fields that may be present in
+other outputs.
+
+Finally, the shift that may have been applied to all the particles upon reading
+the ICs (See :ref:`Parameters_ICs`) is added to the header in the field
+``Shift``. This is expressed in internal units.
+
 Meta-data about the code and run
 --------------------------------
 
@@ -151,10 +212,10 @@ Structure of the particle arrays
 
 There are several groups that contain 'auxiliary' information, such as
 ``Header``.  Particle data is placed in separate groups depending of the type of
-the particles. The type use the naming convention of Gadget-2 (with
-the OWLS and EAGLE extensions). A more intuitive naming convention is
-given in the form of aliases within the file. The aliases are shown in
-the third column of the table.
+the particles. There are currently 6 particle types available. The type use the
+naming convention of Gadget-2 (with the OWLS and EAGLE extensions). A more
+intuitive naming convention is given in the form of aliases within the file. The
+aliases are shown in the third column of the table.
 
 +---------------------+------------------------+-----------------------------+----------------------------------------+
 | HDF5 Group Name     | Physical Particle Type | HDF5 alias                  | In code ``enum part_type``             |
@@ -165,14 +226,22 @@ the third column of the table.
 +---------------------+------------------------+-----------------------------+----------------------------------------+
 | ``/PartType2/``     | Background Dark Matter | ``/DMBackgroundParticles/`` | ``swift_type_dark_matter_background``  |
 +---------------------+------------------------+-----------------------------+----------------------------------------+
+| ``/PartType3/``     | Sinks                  | ``/SinkParticles/``         | ``swift_type_sink``                    |
++---------------------+------------------------+-----------------------------+----------------------------------------+
 | ``/PartType4/``     | Stars                  | ``/StarsParticles/``        | ``swift_type_star``                    |
 +---------------------+------------------------+-----------------------------+----------------------------------------+
 | ``/PartType5/``     | Black Holes            | ``/BHParticles/``           | ``swift_type_black_hole``              |
 +---------------------+------------------------+-----------------------------+----------------------------------------+
+| ``/PartType6/``     | Neutrino Dark Matter   | ``/NeutrinoParticles/``     | ``swift_type_neutrino``                |
++---------------------+------------------------+-----------------------------+----------------------------------------+
 
 The last column in the table gives the ``enum`` value from ``part_type.h``
 corresponding to a given entry in the files.
 
+For completeness, the list of particle type names is stored in the snapshot
+header in the array ``/Header/PartTypeNames``. The number of types (aka. the
+length of this array) is stored as the attribute ``/Header/NumPartTypes``.
+
 Each group contains a series of arrays corresponding to each field of the
 particles stored in the snapshots. The exact list of fields depends on what
 compile time options were used and what module was activated. A full list can be
@@ -180,6 +249,12 @@ obtained by running SWIFT with the ``-o`` runtime option (See
 :ref:`Output_selection_label` for details). Each field contains a short
 description attribute giving a brief summary of what the quantity represents.
 
+Note that the HDF5 names of some fields differ from the GADGET-2 format for
+initial condition files (see :ref:`Initial_Conditions_label`) that mixes
+singular and plural names, which in snapshot files are all plural by default 
+(e.g. ``InternalEnergies`` in snapshots versus ``InternalEnergy`` in initial
+conditions).
+
 All the individual arrays created by SWIFT have had the Fletcher 32 check-sum
 filter applied by the HDF5 library when writing them. This means that any
 eventual data corruption on the disks will be detected and reported by the
@@ -196,9 +271,11 @@ Each particle field contains meta-data about the units and how to
 convert it to CGS in physical or co-moving frames. The meta-data is in
 part designed for users to directly read and in part for machine
 reading of the information. Each field contains the exponent of the
-scale-factor, reduced Hubble constant [#f1]_ and each of the 5 base units
+scale-factor, reduced Hubble constant [#f2]_ and each of the 5 base units
 that is required to convert the field values to physical CGS
-units. These fields are:
+units. The base assumption is that all fields are written in the
+co-moving frame (see below for exceptions).
+These fields are:
 
 +----------------------+---------------------------------------+
 | Meta-data field name | Description                           |
@@ -254,6 +331,56 @@ case of the densities and assuming the usual system of units
 
 In the case of a non-cosmological simulation, these two expressions
 are identical since :math:`a=1`.
+
+In some special cases, the fields cannot be meaningfully expressed as
+co-moving quantities. In these exceptional circumstances, we set the
+value of the attribute ``Value stored as physical`` to ``1``. And we
+additionally set the attribute ``Property can be converted to
+comoving`` to ``0``.
+
+Particle splitting metadata
+---------------------------
+
+When particle splitting is turned on (see :ref:`Parameters_basics`; by using
+``particle_splitting=1`` in the parameter file) some particles in the output
+may have been created from the 'splitting' of a single, over-massive, particle.
+
+There are three fields, associated with all gas, star, and black hole particles,
+that can be used to understand if, and how, these particles were split.
+
+These three fields are:
+
++ ``ProgenitorIDs``, the IDs of the gas particles in the initial conditions
+  that is the direct progenitor of this particle.
++ ``SplitCounts``, the number of times this gas particle has been split; or,
+  if a star or black hole, how many times the gas particle that became this
+  star (or black hole seed) was split before becoming so.
++ ``SplitTrees``, a binary tree (encoded as a 64 bit integer) showing how this
+  particle was split. Each item in the tree shows whether this particle retained
+  its original ID (encoded as 0) or was given a new ID (encoded as 1) in the
+  splitting event. This data is enough to completely reconstruct the splitting 
+  history of the particles.
+
+For example, if a particle has been split 5 times (``SplitCounts=5`` for this
+particle), and has a binary tree of "10010", it retained its original ID in
+the first event, was given a new one in the second event, for the next two
+events it retained its new ID (obtained in the second event), and finally was
+given a new ID in the final event. Throughout this process, the value of
+``ProgenitorIDs`` remained the same. Through this system, we can ensure that
+the combination of ``ProgenitorID`` and this binary tree corresponds to a
+fully traceable, unique, identifier for every particle in the simulation volume.
+
+Note that we can only track 64 splitting events for a given particle, and after
+this the binary tree is meaningless. In practice, however, such a high number of
+splitting events is extremely unlikely to occur. The logging of extra splits can
+optionally be activated. When particles reach 64 splits, their tree information
+is reset but the status prior to the reset is stored in a log file allowing for
+the reconstruction of the full history even in the cases where the maximum is
+reached.
+
+An example is provided in ``examples/SubgridTests/ParticleSplitting``, with a
+figure showing how one particle is split (eventually) into 16 descendants that
+makes use of this metadata.
    
 Quick access to particles via hash-tables
 -----------------------------------------
@@ -286,9 +413,9 @@ expressed in the unit system used for the snapshots (see above) and are hence
 consistent with the particle positions themselves. 
 
 Once the cell(s) containing the region of interest has been located,
-users can use the ``/Cells/Offsets/PartTypeN/Files``,
-``/Cells/Offsets/PartTypeN/Counts`` and
-``/Cells/Offsets/PartTypeN/OffsetsInFile`` to retrieve the location of
+users can use the ``/Cells/Files/PartTypeN/``,
+``/Cells/Counts/PartTypeN/`` and
+``/Cells/OffsetsInFile/PartTypeN/`` to retrieve the location of
 the particles of type ``N`` in the ``/PartTypeN`` arrays.  These
 contain information about which file contains the particles of a given
 cell. It also gives the offset from the start of the ``/PartTypeN``
@@ -305,6 +432,19 @@ In the case of a single-file snapshot, the ``Files`` array is just an array of
 zeroes since all the particles will be in the 0-th file. Note also that in the
 case of a multi-files snapshot, a cell is always contained in a single file.
 
+As noted above, particles can (slightly) drift out of their cells. This can be
+problematic in cases where one wants to find precisely all the particles in a
+given region. To help with this, the meta-data also contains a "cell bounding
+box". The arrays ``/Cells/MinPositions/PartTypeN`` and
+``/Cells/MaxPositions/PartTypeN`` contain the minimal (maximal) x,y,z
+coordinates of all the particles of this type in the cells. Note that these
+coordinates can be outside of the cell itself. When using periodic boundary
+conditions, no box-wrapping is applied.
+
+If a snapshot used a sub-sampled output, then the counts and offsets are
+adjusted accordingly and correspond to the actual content of the file
+(i.e. after the sub-sampling was applied).
+
 As an example, if one is interested in retriving all the densities of the gas
 particles in the cell around the position `[1, 1, 1]` in a single-file
 snapstshot one could use a piece of code similar to:
@@ -354,9 +494,47 @@ from the disk.
 Note that this is all automated in the ``swiftsimio`` python library
 and we highly encourage its use.
 
-.. [#f1] Note that all quantities in SWIFT are always "h-free" in the
-	 sense that they are expressed in units withouy any h
-	 terms. This implies that the ``h-scale exponent`` field value
-	 is always 0. SWIFT nevertheless includes this field to be
-	 comprehensive and to prevent confusion with other software
-         packages that express their quantities with h-full units.
+Meta-file for distributed snapshots
+-----------------------------------
+
+If distributed snapshots are chosen for an MPI parallel run (see
+:ref:`Parameters_snapshots`), N snapshot files are produced, where N is the
+number of MPI ranks. When HDF5 1.10.0 or higher is available, an
+additional meta-snapshot is produced that uses HDF5's virtual dataset
+feature to present these N files as if they were a single, regular
+snapshot file.
+
+The meta-snapshot contains all the meta-data (including the top level
+cell hash-tables) contained in a regular snapshot, but does not store
+any actual particle data. Instead, the particle datasets contain virtual
+links to the corresponding particle data in the distributed snapshot
+files. Since this is a feature of the HDF5 library itself, this is
+entirely transparent to modules like ``h5py`` that try to read the data.
+A user only needs to access the meta-snapshot, and the HDF5 library
+takes care of the rest.
+
+The virtual links in the meta-snapshot only work if the HDF5 library
+knows the location of the distributed snapshots. These are stored within
+the meta-snapshot as relative paths. When SWIFT produces a distributed
+snapshot, all files are placed within the same directory. This means
+that the meta-snapshot can only be safely read if the other N files are
+also present in the same directory.
+
+The header of a meta-snapshot looks exactly like the header of a normal,
+non-distributed snapshot (i.e. ``NumFilesPerSnapshot`` is 1). However,
+the attribute ``Virtual`` is set to 1 to distinguish it from a normal
+snapshot file.
+
+.. [#f1] In the rare case where an output
+	 selection (see :ref:`Output_selection_label`) disabling a given particle type in
+	 its entirety was used, the corresponding entry in ``NumPart_ThisFile`` will be 0
+	 whilst the ``NumPart_Total`` field will still contain the number of
+	 particles present in the run.
+
+
+.. [#f2] Note that all quantities in SWIFT are always "h-free" in the sense that
+	 they are expressed in units withouy any h terms. This implies that the
+	 ``h-scale exponent`` field value is always 0. SWIFT nevertheless
+	 includes this field to be comprehensive and to prevent confusion with
+	 other software packages that express their quantities with h-full
+	 units.

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

+.. AGN spin and jet model
+   Filip Husko, 1 April 2022
+
+.. AGN_spin_jet:
+
+AGN jets and black hole spin in hydrodynamical simulations
+==========================================================
+
+Model summary
+-------------
+
+
+The main feature of this model, in terms of the effects on galaxy populations, is the addition of an AGN jet 
+mode of feedback. In order to launch realistic jets, black hole spin is tracked and evolved for all BH 
+particles in the simulation.
+
+Jet powers, in addition to depending on spin, also depend on which accretion state the black hole is in. We 
+include three accretion states: the thick, thin and slim disk. The thick disk appears at low accretion rates, 
+has very strong jets and is inefficient at spinning up the black hole. The thin disk, appearing at 
+intermediate accretion rates, has weak jets, strong radiation and efficiently spins up the black hole. The 
+slim disk, corresponding to super-Eddington accretion, has features of both, and has both strong radiation and 
+jets. Slim disks can be turned off in the model, but the thick and thin disks are intimitely tied to their 
+feedback modes (jets and radiation, respectively).
+
+In ``theory.rst`` we outline all of the theory which is implemented as part of the model. This includes when 
+the black holes transition from one state to another, the strength of feedback in each state, how spin is 
+evolved in terms of magnitude and direction, etc. In ``numerics.rst`` we discuss how jet launching is 
+implemented, and additional black hole time steps introduced into the code. In ``params.rst`` we list and 
+discuss all parameters used by the model. In ``output.rst`` we list additional arrays output for the BHs and 
+tracers. Below we outline how to configure and run the model.
+
+Compiling and running the model
+-------------------------------
+
+The model can be run with either the EAGLE or COLIBRE models. You can configure the model with ``--with-black-holes=SPIN_JET`` in combination with other configure options, or you can configure the full EAGLE or COLIBRE models with the new spin/jet physics as ``--with-subgrid=SPIN_JET_EAGLE`` and ``--with-subgrid=SPIN_JET_COLIBRE``, respectively. The model will then run as long as ``--black-holes`` is among the runtime options.
+
+For cosmological simulations you do not need to do anything special, but for isolated runs (or any runs with black holes in the initial conditions), the ICs must include two new fields for all black holes: a scalar field representing black hole spins called ``Spins`` and a vector field representing the directions of the spin called ``AngularMomentumDirections``. The former should be between 0 and 1, while the latter should be normalized to 1.
+
+A full list of all relevant parameters of the model is in ``params.rst``. 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.
+
+.. toctree::
+
+  theory
+  numerics
+  params
+  output
+  variable_heating_temperatures

doc/RTD/source/SubgridModels/AGNSpinJets/numerics.rst

+.. AGN spin and jet model
+   Filip Husko, 1 April 2022
+
+.. AGN_spin_jet:
+
+Jet launching algorithm
+-----------------------
+
+In order to launch jets, we introduce a jet reservoir that functions identically to the thermal reservoir used in EAGLE and COLIBRE. When the jet reservoir exceeds the value :math:`N_\mathrm{j}\overline{m}_\mathrm{ngb}v_\mathrm{j}^2/2`, kicks are given out to :math:`N_\mathrm{j}` particles. Here, :math:`N_\mathrm{j}` is the target number of particles to kick per each kicking event (typically equal to 2, but should always be a multiple of 2 to ensure approximately symmetrical jets). :math:`\overline{m}_\mathrm{ngb}v_\mathrm{j}^2/2` is the average kinetic energy of one kicking event, with :math:`\overline{m}_\mathrm{ngb}` the mean neighbour mass of gas particles in the BH smoothing kernel and :math:`v_\mathrm{j}` the target kicking velocity.
+
+These kicks are handed out to particles in a symmetric way with respect to the spin vector of the BH. :math:`N_\mathrm{j}/2` particles are kicked from the 'upper' hemisphere relative to the spin vector, and the other half from the lower hemisphere. The particles to be kicked can be any in the smoothing kernel. We include four different choices: the particles kicked are: 1) the closest to the BH, 2) the farthest from the BH, 3) the ones of minimal density and 4) the ones closest to the spin axis, in terms of angular distance. Note that these sortings are done for each hemisphere seperately. 
+
+The particles chosen are always given velocities based on the same algorithm, regardless of their positions in the kernel. We perform the actual kicks in the following way. Velocity kicks are chosen at random from a cone around the current spin vector with a (half-)opening angle of :math:`\theta_\mathrm{j}`. In particular, we first choose the kick vector around the z-axis as :math:`\hat{v}_\mathrm{kick}=(\sin\theta\cos\phi,\hspace{0.3mm}\sin\theta\sin\phi,\hspace{0.3mm}\cos \theta)`. Here, :math:`\cos\theta` is chosen uniformly from the interval :math:`[\cos\theta_\mathrm{j},1]`, and :math:`\sin\theta=\sqrt{1-\cos\theta^2}`. :math:`\phi` is chosen uniformly from :math:`[0,2\pi]`. This random vector, now representing a random kick within a cone around the z-axis, is rotated into the frame of the spin vector so that the cone is pointing in the right direction. For particles being kicked from the 'negative' side of the BH hemisphere, the final kick vector is simply multiplied by :math:`-1`.
+
+We increase the particle's velocity in the chosen kick direction by an amount :math:`\vert \Delta \vec{v} \vert` such that its energy increases by :math:`(1/2)m_\mathrm{gas}v_\mathrm{jet}^2`. For this reason, :math:`\vert \Delta \vec{v} \vert< v_\mathrm{jet}` generally holds. We calculate :math:`\vert \Delta \vec{v} \vert` from the equation
+
+.. math::
+    (\vec{v}_0+\Delta\vec{v})^2 = \vert \vec{v}_0\vert ^2 + v_\mathrm{j}^2,
+    
+which follows from conservation of energy. This vector equation can be solved to yield the necessary magnitude of the velocity increase that should be applied (in the chosen kick direction)
+
+.. math::
+    \vert \Delta\vec{v}\vert = \sqrt{v_\mathrm{j}^2 + v_\mathrm{0,j}^2} - v_\mathrm{0,j},
+    
+where :math:`v_\mathrm{0,j} = \sin \theta\vert \vec{v}_0\vert` is the magnitude of the initial velocity projected onto the chosen kick direction, with :math:`\sin \theta` the angle between the direction of the initial velocity and the chosen kick direction.
+
+Black hole time steps
+---------------------
+
+Given the changes to the BH model in the form of AGN jets and BH spin evolution, a few additional time-step criteria need to be implemented. The minimum of these time-steps is taken to actually evolve the BH, alongside the other time-steps already used for the BH in the code. We introduce a jet-related time-step that is given by:
+
+.. math::
+    \Delta t_\mathrm{jet}=\frac{\Delta E_\mathrm{jet}}{P_\mathrm{jet}}.
+
+This time-step ensures that the BH is woken up by the time it needs to 'hand out' a pair of kicks. In the above equation, :math:`P_\mathrm{jet}` is the current, instantenous jet power, while :math:`\Delta E_\mathrm{jet}=2\times m_\mathrm{ngb}v_\mathrm{jet}^2` is the energy to be handed out to a pair of particles, with :math:`m_\mathrm{ngb}` the average gas particle mass in the BH's kernel, and :math:`v_\mathrm{jet}` the target jet velocity.
+
+We also introduce two time-steps related to the angular momentum of the BH. The first of these ensures that the magnitude of spin does not change too much over a single time-step, and it is given by
+
+.. math::
+    \Delta t_\mathrm{a}=0.1\frac{\vert a\vert M_\mathrm{BH}}{s \dot{M}_\mathrm{BH,0}},
+
+where :math:`s` is the spinup/spindown function. The numerical factor :math:`0.1` quantifies how finely we want to evolve spin; it ensures that the value of spin changes no more than :math:`10` per cent (relative to the current value) over the next time-step.
+
+We also introduce a time-step related to the redirection of the spin vector. Since the spin vector may be redirected very quickly relative to its magnitude (due to LT torques), this criterion is separate to the one mentioned above. This time-step is given
+
+.. math::
+    \Delta t_\mathrm{a}=0.1\frac{M_\mathrm{warp}J_\mathrm{BH}}{\dot{M}_\mathrm{BH,0}J_\mathrm{warp}\sin\theta},
+
+where :math:`\theta` is the angle between the current BH spin vector and the angular momentum of gas in the accretion disc on large scales. The numerical prefactor is again present to ensure a fine enough evolution of the spin vector direction. In particular, in the case that the spin vector and the gas angular momentum are perpendicular (:math:`\sin\theta=1`), this criterion will lead to a change of no more than :math:`\approx5\degree` in the spin vector direction per time-step.
+

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

+.. AGN spin and jet model
+   Filip Husko, 1 April 2022
+
+.. AGN_spin_jet:
+
+Black holes
+-----------
+
+
++---------------------------------------+-------------------------------------+-----------+-----------------------------+
+| Name                                  | Description                         | Units     | Comments                    |
++=======================================+=====================================+===========+=============================+
+| ``AngularMomentumDirections``         | | The direction of the angular      | [-]       | | Array of length           |
+|                                       | | momentum (spin) vector of the BH  |           | | 3 for each particle       |
++---------------------------------------+-------------------------------------+-----------+-----------------------------+
+| ``AccretionDiscAspectRatios``         | | Aspect ratio, H/R, of the subgrid | [-]       |                             |
+|                                       | | accretion disk surrounding each   |           |                             |
+|                                       | | black hole                        |           |                             |
++---------------------------------------+-------------------------------------+-----------+-----------------------------+
+| ``AccretionModes``                    | | Type of subgrid accretion disk    | [-]       |                             |
+|                                       | | surrounding the black holes       |           |                             |
+|                                       | | 0 - Thick disk, 1 - Thin disk,    |           |                             |
+|                                       | | 2 - Slim disk                     |           |                             |
++---------------------------------------+-------------------------------------+-----------+-----------------------------+
+| ``CosAccretionDiskAngle``             | | Cosine of the angle between the   | [-]       |                             |
+|                                       | | spin vector and the gas angular   |           |                             |
+|                                       | | momentum vector around the BH     |           |                             |
++---------------------------------------+-------------------------------------+-----------+-----------------------------+
+| ``InjectedJetEnergies``               | | Total jet energy injected into    | [U_M U_L  |                             |
+|                                       | | surroundings of this BH           | ^2 U_t^-2]|                             |
++---------------------------------------+-------------------------------------+-----------+-----------------------------+
+| ``JetEfficiencies``                   | | The efficiency of jet launching,  | [-]       |                             |
+|                                       | | i.e. the jet power divided by the |           |                             |
+|                                       | | accretion rate times c*c          |           |                             |
++---------------------------------------+-------------------------------------+-----------+-----------------------------+
+| ``JetReservoirs``                     | | The remaining jet energy left to  | [U_M U_L  |                             |
+|                                       | | be launched from the BH           | ^2 U_t^-2]|                             |
++---------------------------------------+-------------------------------------+-----------+-----------------------------+
+| ``JetTimeSteps``                      | | Jet-limited time steps of the BHs | [U_t]     |                             |
++---------------------------------------+-------------------------------------+-----------+-----------------------------+
+| ``LastAGNJetScaleFactors``            | | Last AGN jet scale factors when   | [-]       |                             |
+|                                       | | the BH did jet feedback, if       |           |                             |
+|                                       | | cosmology is turned on            |           |                             |
++---------------------------------------+-------------------------------------+-----------+-----------------------------+
+| ``LastAGNJetTimes``                   | | Last AGN jet times when the BH    | [U_t]     |                             |
+|                                       | | did jet feedback, if cosmology is |           |                             |
+|                                       | | turned off                        |           |                             |
++---------------------------------------+-------------------------------------+-----------+-----------------------------+
+| ``NumberOfAGNJetEvents``              | | Total number of times this BH did | [-]       |                             |
+|                                       | | jet feedback                      |           |                             |
++---------------------------------------+-------------------------------------+-----------+-----------------------------+
+| ``NumberOfJetParticlesLaunched``      | | Total number of times this BH     | [-]       |                             |
+|                                       | | launched particles as part of jet |           |                             |
+|                                       | | feedback                          |           |                             |
++---------------------------------------+-------------------------------------+-----------+-----------------------------+
+| ``Spins``                             | | Dimensionless spin parameters of  | [-]       |                             |
+|                                       | | the black holes. Negative values  |           |                             |
+|                                       | | indicate retrograde accretion     |           |                             |
++---------------------------------------+-------------------------------------+-----------+-----------------------------+
+| ``RadiativeEfficiencies``             | | The efficiency of radiative       | [-]       |                             |
+|                                       | | feedback, i.e. the radiative      |           |                             |
+|                                       | | power divided by the accretion    |           |                             |
+|                                       | | rate times c*c                    |           |                             |
++---------------------------------------+-------------------------------------+-----------+-----------------------------+
+
+
+Tracers (gas and stars)
+-----------------------
+
++---------------------------------------+-------------------------------------+-----------+-----------------------------+
+| Name                                  | Description                         | Units     | Comments                    |
++=======================================+=====================================+===========+=============================+
+| ``EnergiesReceivedFromJetFeedback``   | | The total energy this particle    | [U_M U_L  |                             |
+|                                       | | has received by being kicked as   | ^2 U_t^-2]|                             |
+|                                       | | part of jet feedback              |           |                             |
++---------------------------------------+-------------------------------------+-----------+-----------------------------+
+| ``LastAGNJetFeedbackScaleFactors``    | | Scale factor when the particle    | [-]       |                             |
+|                                       | | was last kicked as part of jet    |           |                             |
+|                                       | | feedback, if with cosmology       |           |                             |
++---------------------------------------+-------------------------------------+-----------+-----------------------------+
+| ``LastAGNJetFeedbackTimes``           | | Times when the particle was last  | [U_t]     |                             |
+|                                       | | last kicked as part of jet        |           |                             |
+|                                       | | feedback, if without cosmology    |           |                             |
++---------------------------------------+-------------------------------------+-----------+-----------------------------+
+| ``KickedByJetFeedback``               | | How many times this particle has  | [-]       |                             |
+|                                       | | been kicked as                    |           |                             |
+|                                       | | part of jet feedback              |           |                             |
++---------------------------------------+-------------------------------------+-----------+-----------------------------+

doc/RTD/source/SubgridModels/AGNSpinJets/params.rst

+.. AGN spin and jet model
+   Filip Husko, 1 April 2022
+
+.. AGN_spin_jet:
+
+Model parameters
+----------------
+
+Below we give an example of parameter choices applicable for e.g. a 50 Mpc box. The new parameters are from ``include_jets`` and below. Their descriptions are given next to the parameters.
+
+.. code:: YAML
+
+    SPINJETAGN:
+        subgrid_seed_mass_Msun:             1e5        # Black hole subgrid mass at creation time in solar masses.
+        use_subgrid_mass_from_ics:          1          # (Optional) Use subgrid masses specified in ICs [1, default], or initialise them to particle masses [0]?
+        with_subgrid_mass_check:            1          # (Optional) Verify that initial black hole subgrid masses are positive [1, default]. Only used if use_subgrid_mass_from_ics is 1.
+        use_multi_phase_bondi:              0          # Compute Bondi rates per neighbour particle [1] or for the smoothed ambient gas around the black hole [0]?
+        use_subgrid_gas_properties:         1          # Use subgrid density [1] or dynamical density [0] to calculate BH accretion rates?
+        use_krumholz:                       1          # Use Krumholz et al. (2006) [1] or standard Bondi-Hoyle-Lyttleton formula [0] for black hole accretion rates? Only used if multi_phase_bondi is 0.
+        with_krumholz_vorticity:            0          # Include the vorticity term in Krumholz et al. formula? Only used if use_multi_phase_bondi is 0.
+        with_angmom_limiter:                0          # Are we applying the Rosas-Guevara (2015) viscous time-scale reduction term?
+        viscous_alpha:                      1e6        # Normalisation constant of the viscous time-scale in the accretion reduction term. Only used if with_angmom_limiter is 1.
+        with_boost_factor:                  0          # Are we using the model from Booth, Schaye (2009)?
+        boost_alpha:                        1.         # Lowest value for the accretion effeciency for the Booth, Schaye 2009 accretion model.
+        boost_beta:                         2.         # Slope of the power law for the Booth, Schaye 2009 model, set beta to zero for constant alpha models.
+        boost_n_h_star_cm3:                 0.1        # Normalization of the power law for the Booth Schaye 2009 model in cgs (cm^-3).
+        eddington_fraction_for_recording:   0.1        # Record the last time BHs reached an Eddington ratio above this threshold.
+        use_nibbling:                       1          # Continuously transfer small amounts of mass from all gas neighbours to a black hole [1] or stochastically swallow whole gas particles [0]? 
+        min_gas_mass_for_nibbling_Msun:     9e5        # Minimum mass for a gas particle to be nibbled from [M_Sun]. Only used if use_nibbling is 1.
+        coupling_efficiency:                0.15       # Fraction of the radiated energy that couples to the gas in feedback events.
+        AGN_delta_T_K:                      1e8        # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
+        AGN_num_ngb_to_heat:                1.         # Target number of gas neighbours to heat in an AGN feedback event.
+        with_potential_correction:          1          # Subtract BH's own contribution to the potential of neighbours when determining repositioning targets.
+        max_reposition_mass_Msun:           2e8        # Maximal BH mass considered for BH repositioning in solar masses.
+        max_reposition_distance_ratio:      3.0        # Maximal distance a BH can be repositioned, in units of the softening length.
+        with_reposition_velocity_threshold: 0          # Should we only reposition to particles that move slowly w.r.t. the black hole?
+        max_reposition_velocity_ratio:      0.25       # Maximal velocity offset of a particle to reposition a BH to, in units of the ambient sound speed of the BH. Only meaningful if with_reposition_velocity_ratio is 1.
+        min_reposition_velocity_threshold_km_p_s: -1.0 # Minimal value of the velocity threshold for repositioning [km/s], set to < 0 for no effect. Only meaningful if with_reposition_velocity_ratio is 1.
+        set_reposition_speed:               0          # Should we reposition black holes with (at most) a prescribed speed towards the potential minimum?
+        reposition_coefficient_upsilon:     0.001      # Repositioning speed normalisation [km/s/M_sun]. Only meaningful if set_reposition_speed is 1.
+        reposition_exponent_xi:             1.0        # (Optional) Scaling of repositioning velocity with BH subgrid mass (default: 1.0, linear). Only meaningful if set_reposition_speed is 1.
+        threshold_major_merger:             0.333      # Mass ratio threshold to consider a BH merger as 'major'
+        threshold_minor_merger:             0.1        # Mass ratio threshold to consider a BH merger as 'minor'
+        merger_threshold_type:              DynamicalEscapeVelocity   # Type of velocity threshold for BH mergers (CircularVelocity as in EAGLE, EscapeVelocity, or DynamicalEscapeVelocity)
+        merger_max_distance_ratio:          3.0        # Maximal distance over which two BHs can merge, in units of the softening length.
+        AGN_use_deterministic_feedback:     1          # Deterministic (1) or stochastic (0) AGN feedback model
+        AGN_feedback_model:                 Isotropic  # AGN feedback model (Isotropic or MinimumDistance)
+        minimum_timestep_yr:                1000.0     # Minimum time-step of black-hole particles
+        max_eddington_fraction:             1.         # Maximal allowed accretion rate in units of the Eddington rate.
+        include_jets:                       1          # Global switch whether to include jet feedback [1] or not [0].
+        turn_off_radiative_feedback:        0          # Global switch whether to turn off radiative (thermal) feedback [1] or not [0]. This should only be used if 'include_jets' is set to 1, since we want feedback in some form or another.
+        alpha_acc:                          0.2        # Viscous alpha of the subgrid accretion disks. Likely to be within the 0.1-0.3 range. The main effect is that it sets the transition accretion rate between the thin and thick disk, as dot(m) = 0.2 * alpha^2.
+        mdot_crit_ADAF:                     0.01       # The transition normalized accretion rate (Eddington ratio) at which the disc goes from thick (low accretion rates) to thin (high accretion rates). The feedback also changes from kinetic jets to thermal isotropic, respectively.
+        seed_spin:                          0.01       # The (randomly-directed) black hole spin assigned to BHs when they are seeded. Should be strictly between 0 and 1.
+        AGN_jet_velocity_model:             Constant   # How AGN jet velocities are calculated. If 'Constant', a single value is used. If 'BlackHoleMass', then an empirical relation between halo mass and black hole mass is used to calculate jet velocities. 'HaloMass' is currently not supported. 
+        v_jet_km_p_s:                       3160.      # Jet velocity to use if 'AGN_jet_velocity_model' is 'Constant'. Units are km/s.
+        v_jet_BH_mass_scaling_reference_mass_Msun: 1e9 # The reference mass used in the relation between halo mass and BH mass used to calculate jet velocities. Only used if 'AGN_jet_velocity_model' is 'BlackHoleMass'.
+        v_jet_BH_mass_scaling_slope:        0.5        # The slope of the relation between halo mass and BH mass used to calculate jet velocities. Only used if 'AGN_jet_velocity_model' is 'BlackHoleMass'.
+        v_jet_min_km_p_s:                   500        # The minimal jet velocity. This is used if 'AGN_jet_velocity_model' is 'BlackHoleMass', 'MassLoading' or 'Local'.
+        v_jet_max_km_p_s:                   1e4        # The maximal jet velocity. This is used if 'AGN_jet_velocity_model' is 'BlackHoleMass', 'MassLoading' or 'Local'.
+        opening_angle_in_degrees:           7.5        # The half-opening angle of the jet in degrees. Should use values < 15 unless for tests.
+        N_jet:                              2          # Target number of particles to kick as part of a single jet feedback event. Should be a multiple of 2 to ensure approximate momentum conservation (we always kick particles in pairs, one from each 'side' of the BH, relative to the spin vector).
+        AGN_jet_feedback_model:             MinimumDistance # Which particles to kick from the black hole smoothing kernels. Should be 'SpinAxis', 'MinimumDistance', 'MaximumDistance' or 'MinimumDensity'
+        eps_f_jet:                          1.         # Coupling efficiency for jet feedback. No reason to expect this to be less than 1.
+        fix_jet_efficiency:                 0          # Global switch whether to fix jet efficiency to a particular value [1], or use a spin-dependant formula [0].
+        jet_efficiency:                     0.1        # The constant jet efficiency used if 'fix_jet_efficiency' is set to 1.
+        fix_jet_direction:                  0          # Global switch whether to fix the jet direction to be along the z-axis, instead of along the spin vector.
+        accretion_efficiency_mode:          Variable   # How the accretion efficiencies are calculated for the thick accretion disc. If 'Constant', the value of 'accretion_efficiency_thick' will be used. If 'Variable', the accretion efficiency will scale with Eddington ratio.
+        accretion_efficiency_thick:         0.01       # The accretion efficiency (suppression factor of the accretion rate) to use in the thick disc (ADAF), to represent the effects of subgrid ADIOS winds that take away most of the mass flowing through the accretion disc.
+        accretion_efficiency_slim:          1          # The constant accretion efficiency to use in the slim disc, at super-Eddington rates.
+        ADIOS_s:                            0.5        # The exponent of the scaling between accretion efficiency and transition radius of the accretion disc, used if 'accretion_efficiency_mode' is 'Variable'.
+        ADIOS_R_in:                         1e4        # The normalisation (the value) of the transition radius of the accretion disc at the critical Eddington ratio (0.01), used if 'accretion_efficiency_mode' is 'Variable'.
+        fix_radiative_efficiency:           0          # Global switch whether to fix the radiative efficiency to a particular value [1], or use a spin-dependant formula [0]. 
+        radiative_efficiency:               0.1        # The constant jet efficiency used if 'fix_radiative_efficiency' is set to 1. Otherwise, this value is used to define the Eddington accretion rate.
+        TD_region:                          B          # How to treat the subgrid accretion disk if it is thin, according to the Shakura & Sunyaev (1973) model. If set to B, region b will be used. If set to C, region c will be used.
+        include_GRMHD_spindown:             1          # Whether to include high jet spindown rates from GRMHD simulations [1], or use an analytical formula that assumes extraction of energy from the rotational mass/energy of the BH.
+        delta_ADAF:                         0.2        # Electron heating parameter, which controls the strength of radiative feedback in thick disks. Should be between 0.1 and 0.5. This parameter is only used if turn_off_secondary_feedback is set to 0.
+        include_slim_disk:                  1          # Global switch whether to include super-Eddington accretion, modeled as the slim disk. If set to 0, disks will be considered thin even at very large accretion rates.
+        use_jets_in_thin_disc:              1          # Whether to use jets alongside radiation in the thin disc at moderate Eddington ratios.
+        use_ADIOS_winds:                    1          # Whether to include ADIOS winds in the thick disc as thermal isotropic feedback (same channel as thin disc quasar feedback, but with a different efficiency). 
+        slim_disc_wind_factor:              1          # The relative efficiency of slim disc winds at super-Eddington rates. If '1', full winds will be used, while '0' will lead to no winds. Any value in between those can also be used. The wind is implemented in the thermal isotropic feedback channel.
+
+Most of these parameters should work well generally, and should not be changed except for tests. We will discuss only some of the more important ones. You can choose whether to have only the thick and thin disk (low and high BH accretion rates, respectively, separated by a value of ``mdot_crit_ADAF``), or you can also include the slim disk at super-Eddington rates with ``include_slim_disk``. You can control what type of feedback you (do not) want with ``include_jets`` and ``turn_off_radiative_feedback``. If you choose to turn off jets, everything will be modeled as a thin disk (regardless of accretion rate), since jets go hand-in-hand with the thick and the slim disk. Similarly, if you turn off radiation, everything will be treated as a thick disk.
+
+Turning on ``use_jets_in_thin_disc`` or ``use_ADIOS_winds`` will cause jets to also be used in the thin disk and winds (thermal isotropic feedback) in the thick disk. Similarly, ``use_ADIOS_winds`` will lead to winds in the slim disk, with the value of the parameter being used to rescale the formula that is implemented (a value of 1 leading to maximal winds).

doc/RTD/source/SubgridModels/AGNSpinJets/plots.py

+import os
+
+if (
+    os.path.exists("efficiencies.png")
+    and os.path.exists("modes.png")
+    and os.path.exists("spec_ang_mom.png")
+    and os.path.exists("spinup.png")
+):
+    # no need to rerun script
+    exit()
+
+import numpy as np
+from scipy.optimize import fsolve
+
+
+def Z1(x):
+    return 1 + (1 - x ** 2) ** 0.3333 * (
+        (1 + np.absolute(x)) ** 0.3333 + (1 - np.absolute(x)) ** 0.3333
+    )
+
+
+def Z2(x):
+    return np.sqrt(3 * x ** 2 + Z1(x) ** 2)
+
+
+def r_hor(x):
+    return 1 + np.sqrt(1 - x ** 2)
+
+
+def r_isco(x):
+    return 3 + Z2(x) - np.sign(x) * np.sqrt((3 - Z1(x)) * (3 + Z1(x) + 2 * Z2(x)))
+
+
+def eps_NT(x):
+    return 1 - np.sqrt(1 - 2 / 3 * 1 / r_isco(x))
+
+
+def eff_sd(m, a, limit):
+    return limit - (
+        1
+        / 10
+        * 1
+        / m
+        * (
+            0.985 / ((4.627 - 4.445 * a) ** -0.5524 + 1.6 * 1 / m)
+            + 0.015 / ((827.3 - 718.1 * a) ** -0.706 + 1.6 * 1 / m)
+        )
+        * (0.9663 - 0.9292 * a) ** -0.5693
+    ) / (1 - np.sqrt(1 - 2 / (3 * r_isco(a))))
+
+
+def find_root(a, limit):
+    return fsolve(eff_sd, x0=1, args=(a, limit))[0]
+
+
+def A(x):
+    return (0.9663 - 0.9292 * x) ** -0.5693
+
+
+def B(x):
+    return (4.627 - 4.445 * x) ** -0.5524
+
+
+def C(x):
+    return (827.3 - 718.1 * x) ** -0.706
+
+
+def m_dot_crit1(x, limit):
+    C1 = C(x) / B(x)
+    eps_1 = 16 * limit / (0.985 * A(x)) * eps_NT(x)
+    N1 = 0.015 / 0.985
+    res1 = (
+        1.6
+        / B(x)
+        * 1
+        / (2 * C1 * eps_1)
+        * (
+            np.sqrt(
+                (C1 * (1 - eps_1) + N1 - eps_1) ** 2 + 4 * eps_1 * C1 * (N1 - eps_1 + 1)
+            )
+            + C1 * (1 - eps_1)
+            + N1
+            - eps_1
+        )
+    )
+    return res1
+
+
+def m_dot_crit2(x, limit):
+    C1 = C(x) / B(x)
+    eps_1 = 16 * limit / (0.985 * A(x)) * eps_NT(x)
+    N1 = 0.015 / 0.985
+    res1 = (
+        1.6
+        / B(x)
+        * 1
+        / (2 * C1 * eps_1)
+        * (
+            -1
+            * np.sqrt(
+                (C1 * (1 - eps_1) + N1 - eps_1) ** 2 + 4 * eps_1 * C1 * (N1 - eps_1 + 1)
+            )
+            + C1 * (eps_1 - 1)
+            + N1
+            - eps_1
+        )
+    )
+    return res1
+
+
+def beta(alfa):
+    return 1 / (1 + 2 * alfa)
+
+
+def gamma(alpha):
+    return (8 - 3 * beta(alpha)) / (6 - 3 * beta(alpha))
+
+
+def t1(alpha):
+    return -0.2703 * gamma(alpha) + 1.3603
+
+
+def t2(alpha):
+    return -0.94 + 4.475 * (gamma(alpha) - 1.444) - 5.1402 * (gamma(alpha) - 1.444) ** 2
+
+
+def t3(alpha):
+    return -0.84 * np.log10(alpha) - 0.919 + 0.643 * np.exp(-0.209 / alpha)
+
+
+def t4(x):
+    return (0.6365 * r_isco(x) - 0.4828) * (1 + 11.9 * np.exp(-0.838 * r_isco(x) ** 4))
+
+
+def t5(x):
+    return 1.444 * np.exp(-1.01 * r_isco(x) ** 0.86) + 0.1
+
+
+def T(x, alpha):
+    return (
+        0.31
+        * (1 + (t4(x) / r_hor(x)) ** 0.9) ** (t2(alpha) + t3(alpha))
+        * 1
+        / (r_hor(x) - t5(x)) ** (t1(alpha))
+    )
+
+
+def eta(x, alpha):
+    return 1 + gamma(alpha) / (gamma(alpha) - 1) * T(x, alpha)
+
+
+def f1(x):
+    return 0.0871 * r_isco(x) - 0.1082
+
+
+def f2(alpha):
+    return 0.5 - 7.798 * (gamma(alpha) - 1.333) ** 1.26
+
+
+def f3(x):
+    return 0.153 * (r_isco(x) - 0.6) ** 0.3 + 0.105
+
+
+def f4(x, alpha):
+    return (
+        f3(x)
+        * (0.9 * gamma(alpha) - 0.2996)
+        * (1.202 - 0.08 * (np.log10(alpha) + 2.5) ** 2.6)
+    )
+
+
+def f5(alpha):
+    return -1.8 * gamma(alpha) + 4.299 - 0.018 + 0.018 * (np.log10(alpha) + 2) ** 3.571
+
+
+def f6(x, alpha):
+    return (
+        f4(x, alpha)
+        * (((0.14 * np.log10(r_hor(x) ** f5(alpha)) + 0.23) / f4(x, alpha)) ** 10 + 1)
+        ** 0.1
+    )
+
+
+def L_adv(x, alpha):
+    return (
+        f2(alpha)
+        + (f1(x) + 10 ** f6(x, alpha)) * (1.15 - 0.03 * (np.log10(alpha) + 3) ** 2.37)
+    ) / eta(x, alpha)
+
+
+def jet_eff(f_Edd, a):
+    horizon_ang_vel = abs(a) / (2.0 * (1.0 + np.sqrt(1 - a * a)))
+    phi = -20.2 * a ** 3 - 14.9 * a ** 2 + 34.0 * a + 52.6
+    phi = phi * (f_Edd / 1.88) ** 1.29 / (1 + (f_Edd / 1.88) ** 1.29)
+    return (
+        0.05
+        / (4.0 * np.pi)
+        * phi ** 2
+        * horizon_ang_vel ** 2
+        * (1.0 + 1.38 * horizon_ang_vel ** 2 - 9.2 * horizon_ang_vel ** 4)
+    )
+
+
+def s_HD(f_Edd, a):
+    xi = f_Edd * 0.017
+    s_min = 0.86 - 1.94 * a
+    L_ISCO = 0.385 * (1.0 + 2.0 * np.sqrt(3.0 * r_isco(a) - 2.0))
+    s_thin = L_ISCO - 2.0 * a * (1.0 - eps_NT(a))
+    return (s_thin + s_min * xi) / (1 + xi)
+
+
+def s(f_Edd, a):
+    horizon_ang_vel = abs(a) / (2.0 * (1.0 + np.sqrt(1 - a * a)))
+    k_EM = 0.23 * np.ones(np.size(a))
+    k_EM[a > 0] = np.minimum(0.1 + 0.5 * a[a > 0], 0.35 * np.ones(np.size(a[a > 0])))
+
+    s_EM = (
+        -1 * a / abs(a) * jet_eff(f_Edd, a) * (1.0 / (k_EM * horizon_ang_vel) - 2.0 * a)
+    )
+
+    return s_HD(f_Edd, a) + s_EM
+
+
+a = np.arange(-1, 1, 0.0001)
+mdotcrit1 = m_dot_crit1(a, 0.5)
+mdotcrit2 = m_dot_crit2(a, 0.5)
+m_a_90 = [find_root(x, 0.9) for x in a]
+m_a_75 = [find_root(x, 0.75) for x in a]
+m_a_50 = [find_root(x, 0.5) for x in a]
+
+import matplotlib
+import pylab
+
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import matplotlib.gridspec as gridspec
+
+fig = plt.figure(figsize=(8, 6))
+
+plt.style.use("classic")
+plt.fill_between(a, [0.0001 for x in a], [0.01 for x in a], color="blue", alpha=0.2)
+plt.fill_between(a, [0.01 for x in a], [1 for x in a], color="red", alpha=0.2)
+plt.fill_between(a, [1 for x in a], [375 for x in a], color="orange", alpha=0.2)
+plt.ylabel(r"$f_\mathrm{Edd}$", fontsize=24, usetex=True)
+plt.xlabel(r"$a$", fontsize=24, usetex=True)
+plt.tick_params(axis="y", right=True, direction="in")
+plt.yscale("log")
+plt.axis([-1, 1, 0.0001, 100])
+plt.text(-0.22, 0.0008, "Thick disk", fontsize=20)
+plt.text(-0.2, 0.08, "Thin disk", fontsize=20)
+plt.text(-0.2, 8, "Slim disk", fontsize=20)
+plt.minorticks_on()
+plt.tick_params(
+    axis="x",
+    direction="in",
+    bottom=True,
+    top=True,
+    length=8,
+    width=1.2,
+    which="major",
+    labelsize=16,
+)
+plt.tick_params(
+    axis="y",
+    direction="in",
+    left=True,
+    right=True,
+    length=8,
+    width=1.2,
+    which="major",
+    labelsize=16,
+)
+plt.tick_params(
+    axis="x",
+    direction="in",
+    bottom=True,
+    top=True,
+    length=4,
+    width=0.9,
+    which="minor",
+    labelsize=16,
+)
+plt.tick_params(
+    axis="y",
+    direction="in",
+    left=True,
+    right=True,
+    length=4,
+    width=0.9,
+    which="minor",
+    labelsize=16,
+)
+
+plt.savefig("modes.png", bbox_inches="tight")
+plt.close()
+
+a = np.arange(-1, 1, 0.0001)
+phi = -20.2 * a ** 3 - 14.9 * a ** 2 + 34.0 * a + 52.6
+horizon_ang_vel = a / (2 * (1 + np.sqrt(1 - a ** 2)))
+jet_factor = (
+    0.05
+    / (4.0 * np.pi)
+    * 1
+    / 0.3
+    * phi ** 2
+    * horizon_ang_vel ** 2
+    * (1.0 + 1.38 * horizon_ang_vel ** 2 - 9.2 * horizon_ang_vel ** 4)
+)
+Z_1 = np.array(
+    [
+        1 + (1 - x ** 2) ** 0.333 * ((1 + abs(x)) ** 0.333 + (1 - abs(x)) ** 0.333)
+        for x in a
+    ]
+)
+Z_2 = np.array(np.sqrt(3 * a ** 2 + Z_1 ** 2))
+r_iso = 3 + Z_2 - np.sign(np.array(a)) * np.sqrt((3 - Z_1) * (3 + Z_1 + 2 * Z_2))
+eps_TD = 1 - np.sqrt(1 - 2 / (3 * r_iso))
+eps_ADAF1 = 0.144 * (6 / r_iso) * eps_TD * min(1, 0.028 / 0.0044)
+eps_ADAF2 = 0.144 * (6 / r_iso) * eps_TD * min(1, 0.001 / 0.0044)
+Jet_ADAF = jet_factor * 0.3
+Jet_SD = 0.22 * jet_factor
+Jet_TD1 = 10 ** -3 * 0.1 ** (-0.1) * 100 ** 0.2 * 10 ** (2 * 0.1) * jet_factor
+Jet_TD2 = 10 ** -3 * 0.1 ** (-0.1) * 10 ** (-1 * 0.1) * jet_factor
+
+eps_SD1 = (
+    1
+    / 10
+    * 1
+    / 1
+    * (
+        0.985 / ((4.627 - 4.445 * a) ** -0.5524 + 1.6 * 1 / 1)
+        + 0.015 / ((827.3 - 718.1 * a) ** -0.706 + 1.6 * 1 / 1)
+    )
+    * (0.9663 - 0.9292 * a) ** -0.5693
+)
+eps_SD2 = (
+    1
+    / 10
+    * 1
+    / 50
+    * (
+        0.985 / ((4.627 - 4.445 * a) ** -0.5524 + 1.6 * 1 / 50)
+        + 0.015 / ((827.3 - 718.1 * a) ** -0.706 + 1.6 * 1 / 50)
+    )
+    * (0.9663 - 0.9292 * a) ** -0.5693
+)
+
+mdot_bh_ADAF1 = (1 - Jet_ADAF / 4.447) * (1 - eps_ADAF1 - Jet_ADAF)
+mdot_bh_ADAF2 = (1 - Jet_ADAF / 4.447) * (1 - eps_ADAF2 - Jet_ADAF)
+mdot_bh_ADAF3 = (1 - Jet_ADAF / 11.56) * (1 - eps_ADAF1 - Jet_ADAF)
+mdot_bh_ADAF4 = (1 - Jet_ADAF / 11.56) * (1 - eps_ADAF2 - Jet_ADAF)
+mdot_bh_TD1 = (1 - Jet_TD1 / 4.447) * (1 - eps_TD - Jet_TD1)
+mdot_bh_TD2 = (1 - Jet_TD2 / 4.447) * (1 - eps_TD - Jet_TD2)
+
+
+def omega(spin):
+    return spin / (2 * (1 + np.sqrt(1 - spin ** 2)))
+
+
+fig = plt.figure(figsize=(13, 4))
+fig.subplots_adjust(top=1, bottom=0, wspace=0.25)
+gs = gridspec.GridSpec(1, 2, width_ratios=[1, 1])
+plt.style.use("classic")
+
+plt.subplot(gs[0])
+plt.plot(
+    a,
+    100 * 0.005 * (1 + 3 * (phi / 50) ** 2 * (horizon_ang_vel / 0.2) ** 2),
+    linewidth=2,
+    label=r"$\epsilon_\mathrm{wind,thick}$",
+    color="blue",
+)
+plt.plot(
+    a,
+    100 * 0.1 * (1 - np.sqrt(1 - 2 / (3 * r_iso))),
+    linewidth=2,
+    label=r"$\epsilon_\mathrm{f}\epsilon_\mathrm{rad,NT}$ $\mathrm{(thin}$ $\mathrm{disc})$",
+    color="red",
+)
+plt.plot(
+    a,
+    100
+    * 0.0635
+    * (1 + ((1 / 1.88) ** 1.29 / (1 + (1 / 1.88) ** 1.29) * phi / 50) ** 2)
+    * np.maximum((1 - 8 * omega(a) ** 2 + 1 * omega(a)), np.zeros(np.size(a))),
+    linestyle=":",
+    linewidth=1.5,
+    label=r"$\epsilon_\mathrm{wind,slim},$ $f_\mathrm{Edd}=1$",
+    color="orange",
+)
+plt.plot(
+    a,
+    100
+    * 0.0635
+    * (1 + ((10 / 1.88) ** 1.29 / (1 + (10 / 1.88) ** 1.29) * phi / 50) ** 2)
+    * np.maximum((1 - 8 * omega(a) ** 2 + 1 * omega(a)), np.zeros(np.size(a))),
+    linestyle="-.",
+    linewidth=1.5,
+    label=r"$\epsilon_\mathrm{wind,slim},$ $f_\mathrm{Edd}=10$",
+    color="orange",
+)
+plt.plot(
+    a,
+    100
+    * 0.0635
+    * (1 + ((100 / 1.88) ** 1.29 / (1 + (100 / 1.88) ** 1.29) * phi / 50) ** 2)
+    * np.maximum((1 - 8 * omega(a) ** 2 + 1 * omega(a)), np.zeros(np.size(a))),
+    linestyle="--",
+    linewidth=1.5,
+    label=r"$\epsilon_\mathrm{wind,slim},$ $f_\mathrm{Edd}=100$",
+    color="orange",
+)
+plt.plot(
+    a,
+    100
+    * 0.0635
+    * (1 + ((1000 / 1.88) ** 1.29 / (1 + (1000 / 1.88) ** 1.29) * phi / 50) ** 2)
+    * np.maximum((1 - 8 * omega(a) ** 2 + 1 * omega(a)), np.zeros(np.size(a))),
+    linestyle="-",
+    linewidth=1.5,
+    label=r"$\epsilon_\mathrm{wind,slim},$ $f_\mathrm{Edd}=1000$",
+    color="orange",
+)
+
+plt.fill_between(a, eps_ADAF1, eps_ADAF2, color="red", alpha=0.2)
+plt.ylabel(r"$\epsilon_\mathrm{wind}$ $[\%]$", fontsize=24, usetex=True)
+plt.xlabel(r"$a$", fontsize=24, usetex=True)
+plt.tick_params(axis="y", right=True, direction="in")
+pylab.legend(loc="upper left", prop={"size": 12}, ncol=2)
+plt.minorticks_on()
+plt.axis([-1, 1, 0, 25])
+plt.tick_params(
+    axis="x",
+    direction="in",
+    bottom=True,
+    top=True,
+    length=8,
+    width=1.2,
+    which="major",
+    labelsize=16,
+)
+plt.tick_params(
+    axis="y",
+    direction="in",
+    left=True,
+    right=True,
+    length=8,
+    width=1.2,
+    which="major",
+    labelsize=16,
+)
+plt.tick_params(
+    axis="x",
+    direction="in",
+    bottom=True,
+    top=True,
+    length=4,
+    width=0.9,
+    which="minor",
+    labelsize=16,
+)
+plt.tick_params(
+    axis="y",
+    direction="in",
+    left=True,
+    right=True,
+    length=4,
+    width=0.9,
+    which="minor",
+    labelsize=16,
+)
+plt.title("Wind efficiency", fontsize=16)
+
+plt.subplot(gs[1])
+plt.plot(
+    a, 100 * Jet_ADAF, linewidth=2, label=r"$\epsilon_\mathrm{jet,thick}$", color="blue"
+)
+plt.plot(
+    a,
+    100 * Jet_ADAF * ((0.01 / 1.88) ** 1.29 / (1 + (0.01 / 1.88) ** 1.29)) ** 2,
+    linewidth=1.5,
+    linestyle=":",
+    label=r"$\epsilon_\mathrm{jet,thin},$ $f_\mathrm{Edd}=0.01$",
+    color="red",
+)
+plt.plot(
+    a,
+    100 * Jet_ADAF * ((0.1 / 1.88) ** 1.29 / (1 + (0.1 / 1.88) ** 1.29)) ** 2,
+    linewidth=1.5,
+    linestyle="-.",
+    label=r"$\epsilon_\mathrm{jet,thin},$ $f_\mathrm{Edd}=0.1$",
+    color="red",
+)
+plt.plot(
+    a,
+    100 * Jet_ADAF * ((1 / 1.88) ** 1.29 / (1 + (1 / 1.88) ** 1.29)) ** 2,
+    linewidth=1.5,
+    linestyle="--",
+    label=r"$\epsilon_\mathrm{jet,thin},$ $f_\mathrm{Edd}=1$",
+    color="red",
+)
+plt.plot(
+    a,
+    100 * Jet_ADAF * ((10 / 1.88) ** 1.29 / (1 + (10 / 1.88) ** 1.29)) ** 2,
+    linewidth=1.5,
+    linestyle="--",
+    label=r"$\epsilon_\mathrm{jet,slim},$ $f_\mathrm{Edd}=10$",
+    color="orange",
+)
+plt.plot(
+    a,
+    100 * Jet_ADAF * ((100 / 1.88) ** 1.29 / (1 + (100 / 1.88) ** 1.29)) ** 2 - 2,
+    linewidth=1.5,
+    linestyle="-",
+    label=r"$\epsilon_\mathrm{jet,slim},$ $f_\mathrm{Edd}=100$",
+    color="orange",
+)
+
+plt.ylabel(r"$\epsilon_\mathrm{jet}$ $[\%]$", fontsize=24, usetex=True)
+plt.xlabel(r"$a$", fontsize=24, usetex=True)
+plt.tick_params(axis="y", right=True, direction="in")
+pylab.legend(loc="upper left", prop={"size": 15})
+plt.axis([-1, 1, 0, 200])
+plt.minorticks_on()
+plt.tick_params(
+    axis="x",
+    direction="in",
+    bottom=True,
+    top=True,
+    length=8,
+    width=1.2,
+    which="major",
+    labelsize=16,
+)
+plt.tick_params(
+    axis="y",
+    direction="in",
+    left=True,
+    right=True,
+    length=8,
+    width=1.2,
+    which="major",
+    labelsize=16,
+)
+plt.tick_params(
+    axis="x",
+    direction="in",
+    bottom=True,
+    top=True,
+    length=4,
+    width=0.9,
+    which="minor",
+    labelsize=16,
+)
+plt.tick_params(
+    axis="y",
+    direction="in",
+    left=True,
+    right=True,
+    length=4,
+    width=0.9,
+    which="minor",
+    labelsize=16,
+)
+plt.title("Jet efficiency", fontsize=16)
+
+plt.savefig("efficiencies.png", bbox_inches="tight")
+
+L_isco1 = [2 / 3 * 1 / np.sqrt(3) * (1 + 2 * np.sqrt(3 * r_isco(x) - 2)) for x in a]
+
+plt.style.use("classic")
+fig = plt.figure(figsize=(8, 6), linewidth=4)
+plt.plot(a, L_isco1, linewidth=2, label=r"$\ell_\mathrm{ISCO}$", color="red")
+plt.plot(
+    a,
+    0.45 * np.array(L_isco1),
+    linewidth=3,
+    linestyle="--",
+    label=r"$0.45\ell_\mathrm{ISCO}$",
+    color="red",
+)
+plt.plot(
+    a,
+    [L_adv(x, 0.1) for x in a],
+    linewidth=2,
+    label=r"$\ell_\mathrm{adv},\alpha=0.1$",
+    color="green",
+)
+plt.plot(
+    a,
+    [L_adv(x, 0.2) for x in a],
+    linewidth=2,
+    label=r"$\ell_\mathrm{adv},\alpha=0.2$",
+    color="teal",
+)
+plt.plot(
+    a,
+    [L_adv(x, 0.3) for x in a],
+    linewidth=2,
+    label=r"$\ell_\mathrm{adv},\alpha=0.3$",
+    color="purple",
+)
+plt.ylabel(r"$\ell_\mathrm{in}$", fontsize=24, usetex=True)
+plt.xlabel(r"$a$", fontsize=24, usetex=True)
+plt.tick_params(axis="y", right=True, direction="in")
+plt.legend(loc="upper right", prop={"size": 14})
+plt.minorticks_on()
+plt.axis([-1, 1, 0, 5])
+plt.tick_params(
+    axis="x",
+    direction="in",
+    bottom=True,
+    top=True,
+    length=8,
+    width=1.2,
+    which="major",
+    labelsize=16,
+)
+plt.tick_params(
+    axis="y",
+    direction="in",
+    left=True,
+    right=True,
+    length=8,
+    width=1.2,
+    which="major",
+    labelsize=16,
+)
+plt.tick_params(
+    axis="x",
+    direction="in",
+    bottom=True,
+    top=True,
+    length=4,
+    width=0.9,
+    which="minor",
+    labelsize=16,
+)
+plt.tick_params(
+    axis="y",
+    direction="in",
+    left=True,
+    right=True,
+    length=4,
+    width=0.9,
+    which="minor",
+    labelsize=16,
+)
+
+plt.savefig("spec_ang_mom.png", bbox_inches="tight")
+plt.close()
+
+z1 = np.array(
+    [
+        1 + (1 - x ** 2) ** 0.333 * ((1 + abs(x)) ** 0.333 + (1 - abs(x)) ** 0.333)
+        for x in a
+    ]
+)
+z2 = np.array(np.sqrt(3 * a ** 2 + z1 ** 2))
+r_iso = 3 + z2 - np.sign(np.array(a)) * np.sqrt((3 - z1) * (3 + z1 + 2 * z2))
+
+phi = -20.2 * a ** 3 - 14.9 * a ** 2 + 34.0 * a + 52.6
+horizon_ang_vel = a / (2 * (1 + np.sqrt(1 - a ** 2)))
+jet_factor = (
+    0.05
+    / (4.0 * np.pi)
+    * 1
+    / 0.3
+    * phi ** 2
+    * horizon_ang_vel ** 2
+    * (1.0 + 1.38 * horizon_ang_vel ** 2 - 9.2 * horizon_ang_vel ** 4)
+)
+
+da_TD_acc_only = 2 / 3 * 1 / np.sqrt(3) * (
+    1 + 2 * np.sqrt(3 * r_iso - 2)
+) - 2 * a * np.sqrt(1 - 2 / (3 * r_iso))
+da_TD_Benson = (
+    2 / 3 * 1 / np.sqrt(3) * (1 + 2 * np.sqrt(3 * r_iso - 2))
+    - 2 * a * np.sqrt(1 - 2 / (3 * r_iso))
+    - (1.25 * 10 ** -3 * 0.1 ** (-0.1) * 100 ** 0.2 * 10 ** (2 * 0.1) * jet_factor)
+    * 2
+    / a
+    * (np.sqrt(1 - a ** 2))
+    * (1 + np.sqrt(1 - a ** 2))
+)
+da_ADAF_acc_only = L_adv(a, 0.1) - 2 * a
+da_ADAF_Benson = (
+    L_adv(a, 0.1)
+    - 2 * a
+    - (jet_factor * 0.3) * 2 / a * (np.sqrt(1 - a ** 2)) * (1 + np.sqrt(1 - a ** 2))
+)
+eps_SD = (
+    1
+    / 10
+    * 1
+    / 10
+    * (
+        0.985 / ((4.627 - 4.445 * a) ** -0.5524 + 1.6 * 1 / 10)
+        + 0.015 / ((827.3 - 718.1 * a) ** -0.706 + 1.6 * 1 / 10)
+    )
+    * (0.9663 - 0.9292 * a) ** -0.5693
+)
+da_SD_acc_only = L_adv(a, 0.1) - 2 * a * (1 - eps_SD)
+da_SD_Benson = (
+    L_adv(a, 0.1)
+    - 2 * a * (1 - eps_SD)
+    - (jet_factor * 0.22) * 2 / a * (np.sqrt(1 - a ** 2)) * (1 + np.sqrt(1 - a ** 2))
+)
+
+
+fig = plt.figure(figsize=(7, 5))
+plt.style.use("classic")
+
+z1 = np.array(
+    [
+        1 + (1 - x ** 2) ** 0.333 * ((1 + abs(x)) ** 0.333 + (1 - abs(x)) ** 0.333)
+        for x in a
+    ]
+)
+z2 = np.array(np.sqrt(3 * a ** 2 + z1 ** 2))
+r_iso = 3 + z2 - np.sign(np.array(a)) * np.sqrt((3 - z1) * (3 + z1 + 2 * z2))
+da_TD_acc_only = 2 / 3 * 1 / np.sqrt(3) * (
+    1 + 2 * np.sqrt(3 * r_iso - 2)
+) - 2 * a * np.sqrt(1 - 2 / (3 * r_iso))
+da_ADAF_Narayan = (
+    0.45 - 12.53 * a - 7.8 * a ** 2 + 9.44 * a ** 3 + 5.71 * a ** 4 - 4.03 * a ** 5
+)
+
+plt.plot(a, da_ADAF_Narayan, linewidth=2, label="Thick disk", color="blue")
+plt.plot(
+    a,
+    s(0.01, a),
+    linewidth=2,
+    label=r"Thin disk, $f_\mathrm{Edd}=0.01$",
+    linestyle=":",
+    color="red",
+)
+plt.plot(
+    a,
+    s(0.1, a),
+    linewidth=2,
+    label=r"Thin disk, $f_\mathrm{Edd}=0.1$",
+    linestyle="-.",
+    color="red",
+)
+plt.plot(
+    a,
+    s(1, a),
+    linewidth=2,
+    label=r"Thin disk, $f_\mathrm{Edd}=1$",
+    linestyle="--",
+    color="red",
+)
+plt.plot(
+    a,
+    s(10, a),
+    linewidth=2,
+    label=r"Slim disk, $f_\mathrm{Edd}=10$",
+    linestyle="--",
+    color="orange",
+)
+plt.plot(
+    a,
+    s(100, a),
+    linewidth=2,
+    label=r"Slim disk, $f_\mathrm{Edd}=100$",
+    linestyle="-",
+    color="orange",
+)
+plt.plot(a, [0 for x in a], linewidth=1.0, color="black", linestyle="--")
+plt.plot([-0.0001, 0.0001], [-200, 200], linewidth=1.0, color="black", linestyle="--")
+plt.ylabel(
+    r"$\mathrm{d}a/(\mathrm{d} M_\mathrm{BH,0}/M_\mathrm{BH})$", fontsize=24, usetex=True
+)
+plt.xlabel("$a$", fontsize=24, usetex=True)
+plt.tick_params(axis="y", right=True, direction="in")
+pylab.legend(loc="lower left", prop={"size": 13})
+plt.minorticks_on()
+plt.axis([-1, 1, -10, 10])
+plt.tick_params(
+    axis="x",
+    direction="in",
+    bottom=True,
+    top=True,
+    length=8,
+    width=1.2,
+    which="major",
+    labelsize=16,
+)
+plt.tick_params(
+    axis="y",
+    direction="in",
+    left=True,
+    right=True,
+    length=8,
+    width=1.2,
+    which="major",
+    labelsize=16,
+)
+plt.tick_params(
+    axis="x",
+    direction="in",
+    bottom=True,
+    top=True,
+    length=4,
+    width=0.9,
+    which="minor",
+    labelsize=16,
+)
+plt.tick_params(
+    axis="y",
+    direction="in",
+    left=True,
+    right=True,
+    length=4,
+    width=0.9,
+    which="minor",
+    labelsize=16,
+)
+
+plt.savefig("spinup.png", bbox_inches="tight")

doc/RTD/source/SubgridModels/AGNSpinJets/theory.rst

+.. AGN spin and jet model
+   Filip Husko, 1 April 2022
+
+.. AGN_spin_jet:
+
+
+Model summary
+-------------
+
+Here we provide a comprehensive summary of the model. In order to more easily visualize the model, it is recommended to run the python script in order to generate plots that show various aspects of the model.
+
+Any model for realistic AGN jets must include black hole spin since jet powers depend steeply on spin, and because including spin provides a well-defined direction for the jets to be launched in. The spin (angular momentum) of BHs is best represented through the dimensionlesss spin parameter :math:`a=J_\mathrm{BH}c/M_\mathrm{BH}^2 G`, where :math:`J_\mathrm{BH}` is its angular momentum. For theoretical reasons, its magnitude cannot grow above 1. It can be positive, representing prograde accretion, or negative, representing retrograde accretion.
+
+Jet powers, in addition to depending on spin, also depend on which accretion state the black hole is in. We refer to these states by the shape of the accretion disk that surrounds the BH. We include three accretion states: the thick (or advection-dominated accretion flow; ADAF), thin (standard) and slim disk (super-Eddington accretion). Our main reference points for these disks are the following papers: `Shakura & Sunyaev (1973) <https://ui.adsabs.harvard.edu/abs/1973A%26A....24..337S/abstract>`_, `Narayan & Yi (1994) <https://ui.adsabs.harvard.edu/abs/1994ApJ...428L..13N/abstract>`_ and `Wang & Zhou. (1999) <https://ui.adsabs.harvard.edu/abs/1999ApJ...516..420W/abstract>`_, respectively.
+
+The thick disk appears at low accretion rates, has very strong jets and is inefficient at spinning up the black hole. The thin disk, appearing at intermediate accretion rates, typically has weak jets, strong radiation and efficiently spins up the black hole. The slim disk, corresponding to super-Eddington accretion, has features of both: in terms of geometry, orbits and angular momentum, it is similar to the thick disk. It is optically thin, leading to strong radiation. However, it also has strong jets. We assume that each subgrid accretion disk launches jets and radiates at the same time, regardless of the type it is. However, we use expressions for the jet and radiative efficiencies that depend on the type of the disk, and which are physically motivated.
+
+Transitions from one accretion state to another
+-----------------------------------------------
+
+.. figure:: modes.png
+    :width: 600px
+    :align: center
+    :figclass: align-center
+    :alt: Accretion regimes
+
+    The type of accretion disk surrounding the BHs depending on their accretion rates and spins.
+
+The state of the subgrid accretion disk depends mostly on the Eddington fraction, i.e. the (dimensionless) accretion rate of the BH in units of the Eddington accretion rate, which we denote as :math:`f_\mathrm{Edd}`. We assume that the subgrid accretion disk is thick for :math:`f_\mathrm{Edd}<f_\mathrm{Edd,crit,thick}`, where :math:`f_\mathrm{Edd,crit,thick}\approx0.01-0.03` (based on observations; `Russell et al. 2013 <https://ui.adsabs.harvard.edu/abs/2013MNRAS.432..530R/abstract>`_) is a free parameter in the model. The accretion disc is assumed to be slim for :math:`f_\mathrm{Edd}>1`.
+
+Accretion efficiencies
+-----------------------------------------------
+
+Our model requires the usage of accretion efficiencies, minimally in the thick disk regime. These accretion efficiencies arise due to winds that take away most of the accreting mass as it falls towards the BH. We supress the large-scale accretion rate (e.g. the Bondi rate), :math:`\dot{M}_\mathrm{BH,0}:` such that the net accretion rate is equal to
+
+.. math::
+    \dot{M}_\mathrm{BH} = (1 - \epsilon_\mathrm{rad} - \epsilon_\mathrm{wind} - \epsilon_\mathrm{jet})\epsilon_\mathrm{acc}\dot{M}_\mathrm{BH,0},
+
+where the terms in parenthesis are feedback efficiencies (discussed below), defined as :math:`\epsilon_i=P_i/\dot{M}_\mathrm{BH}c^2`, while :math:`\epsilon_\mathrm{acc}` is the accretion efficiency.
+
+In the thick and slim disk, we allow options where the accretion efficiencies are free parameters in the model, to be tuned somehow (in practice, for the thick disc efficiency, to the local AGN bolometric luminosity function). We also allow a more complex scaling for the thick disk, based on GRMHD simulations (e.g. `Cho et al. 2024 <https://arxiv.org/abs/2405.13887>`_). These simulations show that the accretion efficiency in thick disks can be calculated as
+
+.. math::
+    \epsilon_\mathrm{acc} = \bigg(\frac{R_0}{R_\mathrm{thick}}\bigg)^s,
+    
+where :math:`R_0\approx5-10` (in units of :math:`R_\mathrm{G}`), :math:`R_\mathrm{thick}` is the radius of the hot accretion flow and :math:`s=0.5` in recent such simulations. The radius :math:`R_\mathrm{thick}` is the same as the Bondi radius if the accretion rate is very low. However, at high accretion rates (but still below :math:`f_\mathrm{Edd}\approx0.01`), the accretion disk may be thin at large distances and thick at smaller ones, with the transition occuring at some radius :math:`R_\mathrm{tr}`. In this case, the winds (and mass loading) operate only at smaller radii. Given these considerations, we write
+
+.. math::
+    R_\mathrm{thick} = \min(R_\mathrm{B},R_\mathrm{tr}),
+    
+and we choose to parametrize :math:`R_\mathrm{tr}` based on original calculations by `Narayan & Yi (1995) <https://ui.adsabs.harvard.edu/abs/1995ApJ...452..710N/abstract>`_, who found the transition radius as the radius where half of the energy gets radiated away, and half advected inwards. Their formula can be written in the form
+
+.. math::
+    R_\mathrm{tr} = R_1 \bigg(\frac{0.01}{f_\mathrm{Edd}}\bigg)^2,
+    
+where :math:`R_1` is some normalisation radius that depends strongly on the assumed value of the accretion disk viscosity parameter :math:`\alpha`. We instead leave :math:`R_1` as a free parameter in our formulation. Note that at moderate Eddington ratios, where the Bondi radius is not the limiting radius (i.e. where a thin disk component exists outside the thick disk), we may write the accretion efficiency as:
+
+.. math::
+    \epsilon_\mathrm{acc} = \sqrt{\frac{R_0}{R_1}}\bigg(\frac{f_\mathrm{Edd}}{0.01}\bigg),
+    
+where we have assumed :math:`s=0.5`.
+
+Jet efficiencies
+----------------
+
+The jet efficiency is related to the jet power through :math:`\epsilon_\mathrm{j}=P_\mathrm{j}/\dot{M}_\mathrm{BH,0}c^2`, where :math:`\dot{M}_\mathrm{BH,0}` is the accretion rate measured in the simulation, e.g. the Bondi rate). We use the formula for the jet efficiency based on general-relativistic, magneto-hydrodynamical (GRMHD) simulations by `Tchekhovskoy et al. (2010) <https://ui.adsabs.harvard.edu/abs/2010ApJ...711...50T/abstract>`_:
+
+.. math::
+    \epsilon_\mathrm{j}=\frac{\kappa}{4\pi} \phi_\mathrm{BH}^2\Omega_\mathrm{BH}^2\big(1+1.38\Omega_\mathrm{BH}^2-9.2\Omega_\mathrm{BH}^4\big),
+
+where :math:`\kappa\approx0.05` is a numerical factor which depends on the initial geometry of the magnetic field, :math:`\phi_\mathrm{BH}` is the dimensionless magnetic flux threading the horizon (see original paper for precise definition), and :math:`\Omega_\mathrm{BH}=a/2r_\mathrm{H}` is the (dimensionless) angular velocity of the black hole event horizon. Here, :math:`r_\mathrm{H}=1+\sqrt{1-a^2}` is the radius of the horizon in units of the gravitational radius :math:`R_\mathrm{G}=M_\mathrm{BH}G/c^2`. The formula above, for the jet efficiency, agrees very well with the results from higher-resolution simulations performed by `Narayan et al. (2021) <https://ui.adsabs.harvard.edu/abs/2010ApJ...711...50T/abstract>`_, who provide the following fit for the magnetic flux as a function of spin:
+
+.. math::
+    \phi_\mathrm{BH,MAD}(a)=-20.2a^3-14.9a^2+34a+52.6.
+    
+The `Tchekhovskoy et al. (2010) <https://ui.adsabs.harvard.edu/abs/2010ApJ...711...50T/abstract>`_ jet efficiency depends very steeply on spin (:math:`\epsilon_\mathrm{j}\propto a^2` for small spin and :math:`\epsilon_\mathrm{j}\propto a^6` near :math:`a=1`). It can reach values above 100 per cent for large spins, and is also different (weaker) for negative spins.
+
+The dependence of the jet efficiency on the type of accretion disk is encoded in the fact that thick disks are thought to be in a magnetically-arred state (so-called MAD. see `Narayan et al. 2003 <https://ui.adsabs.harvard.edu/abs/2003PASJ...55L..69N/abstract>`_), while thin disks are likely not, because they do not feature strong advection. The slim disk, on the other hand, is thought to be similar to the thick disk in terms of advection, and thus probably in terms of jet powers. Recent simulations by `Ricarte et al. (2023) <https://ui.adsabs.harvard.edu/abs/2023ApJ...954L..22R/abstract>`_ have found an increase of :math:`\phi_\mathrm{BH}` in the thin and slim disk regime as the Eddington ratio increases, and they parametrise this increase as
+
+.. math::
+    \phi_\mathrm{BH,thin,slim} = \frac{(f_\mathrm{Edd}/1.88)^{1.29}}{1+(f_\mathrm{Edd}/1.88)^{1.29}}\phi_\mathrm{BH,MAD}.
+
+The magnetic flux eventually saturates (at very high :math:`f_\mathrm{Edd}`) at the same value as that reached in the thick disc; :math:`\phi_\mathrm{BH,MAD}`.
+
+.. figure:: efficiencies.png
+    :width: 1200px
+    :align: center
+    :figclass: align-center
+    :alt: Efficiencies
+
+    Feedback efficiencies (jet - blue, radiation - red) for all three accretion disk types. Shaded regions represent likely ranges of efficiencies (where the efficiencies depend on mass and/or accretion rate). The thin disk jet efficiencies were computed assuming the slope of the efficiency vs. aspect ratio relation is :math:`\eta=1`, and the aspect ratios were computed for region b) of the Shakura & Sunyaev solution. Radiative efficiencies in the thick disk were computed assuming the electron heating parameter :math:`\delta=0.2`.
+
+Radiative/wind efficiencies
+---------------------------
+
+In the EAGLE and COLIBRE models, all subgrid accretion disks are effectively thin, and the BH is always assumed to be in this regime. In our model, the radiative efficiency (defined in an analagous way to the jet efficiency, but using the luminosity) is no longer fixed at a value of order :math:`10` per cent. Instead, we use spin-dependant formulas that vary with the type of disk. In the thin disk, the radiative efficiency :math:`\epsilon_\mathrm{r,TD}` is related to the binding energy at the innermost stable circular orbit (ISCO) and is given by
+
+.. math::
+    \epsilon_\mathrm{r,TD}(a) = 1-e_\mathrm{ISCO}(a)=1-\sqrt{1-\frac{2}{3r_\mathrm{ISCO}(a)}}.
+    
+Here, :math:`r_\mathrm{ISCO}` is the radius of the ISCO in gravitational radii (see e.g. appendix B of `Fiacconi et al. 2018 <https://ui.adsabs.harvard.edu/abs/2018MNRAS.477.3807F/abstract>`_ for an expression giving the spin dependence). The radiative efficiency of the thin disk grows slowly from its minimum value of :math:`\approx4` per cent for :math:`a=-1` to :math:`\approx5.5` per cent for :math:`a=0`. For positive spins it grows more steeply; it is :math:`10` per cent by :math:`a=0.65`. Beyond that the dependence steepens even further, with values of :math:`20`, :math:`30` and :math:`40` per cent reached at :math:`a=0.95`, :math:`a=0.997` and :math:`a=1`, respectively.
+
+In the thick disk regime, radiative efficiencies are lower by a factor :math:`\approx100` than jet efficiencies. The formulas we use are based on results by `Mahadevan (1997) <https://ui.adsabs.harvard.edu/abs/1997ApJ...477..585M/abstract>`_, who studied cooling processes of electrons (which dominate in the radiation) in the context of the original thick disc solution. They found two different regimes: for :math:`f_\mathrm{Edd}<f_\mathrm{Edd,crit,visc}`, viscous heating dominates the heating of electrons, whereas for :math:`f_\mathrm{Edd,crit,visc}<f_\mathrm{Edd}<f_\mathrm{Edd,crit,ADAF}`, it is dominated by ion-electron heating. Here, :math:`f_\mathrm{Edd,crit,visc}` is the transitional value between the two thick disc (ADAF) regimes, and :math:`f_\mathrm{Edd,crit,ADAF}=0.4\alpha^2` is the transitional accretion rate which separates thin and thick discs. The radiative efficiency in the viscous heating regime is given by
+
+.. math::
+    \epsilon_\mathrm{r,ADAF}=0.0002\epsilon_\mathrm{r,TD}\bigg(\frac{\delta_\mathrm{ADAF}}{0.0005}\bigg)\bigg(\frac{1-\beta}{0.5}\bigg)\bigg(\frac{6}{r_\mathrm{ISCO}}\bigg),
+
+while in the ion-heating regime it is given by
+
+.. math::
+    \epsilon_\mathrm{r,ADAF}=0.2\epsilon_\mathrm{r,TD}\bigg(\frac{f_\mathrm{Edd}}{\alpha^2}\bigg)\bigg(\frac{\beta}{0.5}\bigg)\bigg(\frac{6}{r_\mathrm{ISCO}}\bigg).
+    
+Here, :math:`\beta` is the ratio of gas pressure and total pressure (which includes the magnetic pressure). `Yuan & Narayan (2014) <https://ui.adsabs.harvard.edu/abs/2014ARA%26A..52..529Y/abstract>`_ define a somewhat different parameter, :math:`\beta_\mathrm{ADAF}`, as the ratio of gas pressure and magnetic pressure. The two parameters are related by :math:`\beta=\beta_\mathrm{ADAF}/(1+\beta_\mathrm{ADAF})`. :math:`\beta_\mathrm{ADAF}` is not an independent parameter; many simulations have found that :math:`\alpha\beta_\mathrm{ADAF}\approx0.5` (e.g. `Begelman et al. 2021 <https://ui.adsabs.harvard.edu/abs/2022MNRAS.511.2040B/abstract>`_, see also `Yuan & Narayan 2014 <https://ui.adsabs.harvard.edu/abs/2014ARA%26A..52..529Y/abstract>`_ for a review), which we adopt. :math:`\delta_\mathrm{ADAF}` represents the fraction of viscous energy transferred to the electrons, and is constrained in theoretical studies between 0.1 and 0.5 (`Yuan & Narayan 2014 <https://ui.adsabs.harvard.edu/abs/2014ARA%26A..52..529Y/abstract>`_, `Sharma et al. 2007 <https://ui.adsabs.harvard.edu/abs/2007ApJ...667..714S/abstract>`_). Observations imply a value close to 0.2 (`Yuan et al. 2003 <https://ui.adsabs.harvard.edu/abs/2003ApJ...598..301Y/abstract>`_, `Liu & Wu 2013 <https://ui.adsabs.harvard.edu/abs/2013ApJ...764...17L/abstract>`_). The critical accretion rate between the two thick disc regimes can be found by ensuring that both formulas presented above yield the same radiative efficiency (at that accretion rate). This gives an accretion rate equal to
+
+.. math::
+    f_\mathrm{Edd,crit,visc}=0.0002\bigg(\frac{\delta_\mathrm{ADAF}}{0.0005}\bigg)\bigg(\frac{1-\beta}{\beta}\bigg)\alpha^2.
+    
+For slim disks we take the radiative efficiency based on GRMHD simulations of super-Eddington accretion (for various BH spins) performed by `Sadowski et al. (2014) <https://ui.adsabs.harvard.edu/abs/2014MNRAS.439..503S/abstract>`_. `Madau et al. (2014) <https://ui.adsabs.harvard.edu/abs/2014ApJ...784L..38M/abstract>`_ found the following fitting function which represents the `Sadowski et al. (2014) <https://ui.adsabs.harvard.edu/abs/2014MNRAS.439..503S/abstract>`_ results:
+
+.. math::
+    \epsilon_\mathrm{r,SD}=\frac{0.1}{f_\mathrm{Edd}}A(a)\bigg( \frac{0.985}{1.6/f_\mathrm{Edd}+B(a)}+\frac{0.015}{1.6/f_\mathrm{Edd}+C(a)}\bigg),
+    
+where the three spin-dependant functions are given by :math:`A(a)=(0.9663-0.9292a)^{-0.5639}`, :math:`B(a)=(4.627-4.445a)^{-0.5524}` and :math:`C(a)=(827.3-718.1a)^{-0.7060}`. The radiative efficiency of slim disks, based on this formula, matches the thin disk radiative efficiency (given at the beginning of the section) at low accretion rates. At high accretion rates (:math:`f_\mathrm{Edd}\gtrapprox1`, but depending on spin), the radiative efficiency drops.
+
+The thin disc radiative efficiency is used to source feedback in the simulations. In the thin disk regime, a fraction :math:`\epsilon_\mathrm{f}\approx0.1` of all of the radiation released by black holes couples to the gas in the form of thermal energy. In the thick and slim disk, we do not use radiation to source feedback. We do, however, assume that winds launched from the accretion disk are present in these two states. In the thick disk, winds are thought to be launched on account of a combination of gas pressure and MHD effects. We use the formula from `Sadowski et al. (2013) <https://ui.adsabs.harvard.edu/abs/2013MNRAS.436.3856S/abstract>`_:
+
+.. math::
+    \epsilon_\mathrm{wind,thick} = 0.005\bigg[1+0.3\bigg(\frac{\phi_\mathrm{BH,MAD}}{50}\bigg)\bigg(\frac{\Omega_\mathrm{H}}{0.2}\bigg) \bigg].
+	
+For the slim disk, we again use results from `Ricarte et al. (2023) <https://ui.adsabs.harvard.edu/abs/2023ApJ...954L..22R/abstract>`_, as we did for the jet efficiency. We use their total MHD efficiency and subtract from that the analytical jet efficiency as given by the formula we use as a function of spin and magnetic flux. We then found a simple fitting function to the remaining efficiency, representing the wind:
+
+.. math::
+    \epsilon_\mathrm{wind,slim} = 0.065\bigg[1+\bigg(\frac{\phi_\mathrm{BH,thin,slim}}{50}\bigg)^2\bigg] \big(1+\Omega_\mathrm{H}-8\Omega_\mathrm{H}^2\big).
+
+Evolution of the black hole spin magnitude
+------------------------------------------
+
+The BH spin (or angular momentum) is, naturally, a vector. However, due to Lense-thirring torques (we discuss these in more detail below), the accretion disk is always either aligned or counteraligned with the rotational axis of the black hole. This means that almost all relevant quantities, such as the efficiencies discussed above, can be expressed as depending only on the magnitude of spin, but also allowing for a negative sign to account for counteraligned disks (retrograde accretion). This is also true for the evolution of the magnitude of spin.
+
+In the absence of jet spindown, the evolution of angular momentum is given simply by :math:`\dot{J}_\mathrm{BH}=L_\mathrm{in}\dot{M}_\mathrm{BH}`, where :math:`L_\mathrm{in}` is the specific angular momentum at the inner radius of the accretion disk. This can be transformed into an equation for spin evolution, yielding
+
+.. math::
+    \frac{\mathrm{d}a}{\mathrm{d}\ln M_\mathrm{BH,0}}=\ell_\mathrm{in}-2a e_\mathrm{in},
+    
+where :math:`\ell_\mathrm{in}` is the specific angular momentum in units where :math:`G` and :math:`c` are equal to unity, and :math:`\mathrm{d}\ln M_\mathrm{BH,0}=\mathrm{d}M_\mathrm{BH,0}/M_\mathrm{BH}` is the logarithmic change in mass, not including losses due to radiation (`Fanidakis et al. 2011 <https://ui.adsabs.harvard.edu/abs/2011MNRAS.410...53F/abstract>`_). The specific binding energy can be related to the radiative efficiency through :math:`e_\mathrm{in}=1-\epsilon_\mathrm{r}` for all three accretion states (for the thick disc, the radiative efficiency is negligible for this application). All of the above quantities are evaluated at some inner radius beyond which gas orbits are unstable.
+
+To be consistent with what we assumed for feedback efficiencies, we take results for the spinup/spindown function directly from GRMHD simulations. For the thick disc, we use the formula from `Narayan et al. (2021) <https://ui.adsabs.harvard.edu/abs/2010ApJ...711...50T/abstract>`_:
+
+.. math::
+    \bigg(\frac{\mathrm{d}a}{\mathrm{d}M_\mathrm{BH,0}/M_\mathrm{BH}}\bigg)_\mathrm{thick}=0.45 - 12.53a - 7.8a^2 +9.44a^3 + 5.71a^4 -4.03a^5.
+  
+For the slim and thin disc, we use results from `Ricarte et al. (2023) <https://ui.adsabs.harvard.edu/abs/2023ApJ...954L..22R/abstract>`_, who find a fitting formula that smoothly interpolates between the thin disc regime without significant jet feedback (for :math:`f_\mathrm{Edd}` not close to super-Eddington values), and that where jet feedback essentially matches the thick disc (and so jet spindown should also be similar). Their formula takes the form
+
+.. math::
+    \bigg(\frac{\mathrm{d}a}{\mathrm{d}M_\mathrm{BH,0}/M_\mathrm{BH}}\bigg)_\mathrm{thin/slim}=s_\mathrm{HD} - s_\mathrm{EM},
+    
+where the first term is a pure hydrodynamical term, while the second is an electromagnetic term. The first term is given by
+
+.. math::
+    s_\mathrm{HD}=\frac{s_\mathrm{thin}+s_\mathrm{min}\xi}{1+\xi},
+
+where :math:`\xi=0.017f_\mathrm{Edd}`, :math:`s_\mathrm{min}=0.86-1.94a` and :math:`s_\mathrm{thin}=\ell_\mathrm{ISCO}-2a e_\mathrm{ISCO}` is the spinup/spindown function of the 'pure' thin disc (with no outflows and outside the MAD regime), in which :math:`\ell_\mathrm{ISCO}` and :math:`e_\mathrm{ISCO}` are the (dimensionless) specific angular momentum and binding energy, respectively, at the ISCO. The EM term is given by
+
+.. math::
+    s_\mathrm{EM}=\mathrm{sgn}(a)\epsilon_\mathrm{EM}\bigg(\frac{1}{k\Omega_\mathrm{H}}-2a\bigg),
+    
+where :math:`\epsilon_\mathrm{EM}` is the total (jet+wind) EM efficiency, and :math:`k` is given by 
+
+.. math::
+    k=\min(0.35,0.1+0.5a)
+    
+for positive spins :math:`a>0` and by :math:`k=0.23` for negative spins :math:`a<0`.
+
+.. figure:: spinup.png
+    :width: 1200px
+    :align: center
+    :figclass: align-center
+    :alt: Spinup/spindown function
+
+    Spinup/spindown function (the dimensionless rate of black hole spin evolution) as a function of spin for all three accretion disk types. For the thin and slim disk, we show several curves for different choices of the Eddington ratio.
+
+Evolution of the black hole spin direction
+------------------------------------------
+
+In the previous section we claimed that the evolution of the magnitude of spin depends only on whether accretion is prograde or retrograde. The two remaining questions are: 1) what about its direction, and 2) how to decide whether accretion is prograde or retrograde. We will now address the first of these.
+
+Lense-Thirring torques (`Lense & Thirring 1918 <https://ui.adsabs.harvard.edu/abs/1918PhyZ...19..156L/abstract>`_) arise from additional GR forces that operate near spinning BHs, related to the frame-dragging of spacetime. In isolation, they cause the precession of a parcel of gas as it orbits around the BH. For accretion disks, their effects depend on the type of disk (see `Nixon & King 2016 <https://ui.adsabs.harvard.edu/abs/2016LNP...905...45N/abstract>`_ for a review). Lense-Thirring torques do not have a component in the direction of the BH spin vector, which is why they do not play a role in the evolution of the magnitude of spin.
+
+In all cases, Lense-Thirring torques are effective only within some radius :math:`R_\mathrm{warp}`, which marks the boundary between the outer disk and an inner region, within which the BH can 'communicate' through these torques with the disk. Within this radius, the disk is on average aligned or counteraligned with the BH, whereas outside it, it is aligned with some large-scale angular momentum direction (which we can measure in the simulation) - hence the name warp radius. Given some surface density, one can also define the warp mass :math:`M_\mathrm{warp}` and the warp angular momentum :math:`J_\mathrm{warp}` as the total mass and angular momentum within :math:`R_\mathrm{warp}`, respectively. We will discuss how all of these warp-related quantities are calculated in each of the accretion disks further below, but for now we focus on how these warped disks feature in our model.
+
+In terms of the evolution of the spin direction, the main assumption of our model is as follows (see `King et al. 2005 <https://ui.adsabs.harvard.edu/abs/2005MNRAS.363...49K/abstract>`_ for the original argument, and additional discussions in e.g. `Fanidakis et al. 2011 <https://ui.adsabs.harvard.edu/abs/2011MNRAS.410...53F/abstract>`_, `Fiacconi et al. 2018 <https://ui.adsabs.harvard.edu/abs/2018MNRAS.477.3807F/abstract>`_ and `Griffin et al. 2019a <https://ui.adsabs.harvard.edu/abs/2019MNRAS.487..198G/abstract>`_). All matter that flows through an accretion disk is aligned or counteraligned with the BH spin vector in the accretion process. Due to conservation of angular momentum, the spin vector itself also has to adjust to keep the total angular momentum conserved. In the process of consuming one warp mass :math:`M_\mathrm{warp}`, the direction of the BH spin vector is aligned to match the direction of the total angular momentum of the system comprising the BH and the disk out to the warp radius. The direction of the BH spin vector can then be determined from :math:`\vec{J}_\mathrm{warp}=\vec{J}_\mathrm{BH}+J_\mathrm{warp}\hat{J}_\mathrm{d}`, where :math:`\vec{J}_\mathrm{BH}` is the old BH angular momentum vector, and :math:`\hat{J}_\mathrm{d}` is the direction of the large-scale accretion disk (which we assume matches the direction of the angular momentum of the gas in the BH smoothing kernel).
+
+In practice, the BH will consume parcels of mass that differ from :math:`M_\mathrm{warp}`. We assume that any such parcel of mass :math:`\Delta M` (e.g. the mass to be consumed within a single time step) can be split up onto :math:`n=\Delta M / M_\mathrm{warp}` individual increments of accretion, so the total angular momentum of the system within that time step is :math:`\vec{J}_\mathrm{warp}=\vec{J}_\mathrm{BH}+n J_\mathrm{warp}\hat{J}_\mathrm{d}`, i.e. :math:`n` warp angular momenta are consumed, with an angular momentum of :math:`\Delta \vec{J}=n J_\mathrm{warp}\hat{J}_\mathrm{d}=(J_\mathrm{warp}/M_\mathrm{warp})\Delta M`. This can also be viewed as the BH consuming material with a specific angular momentum of :math:`L_\mathrm{warp}=J_\mathrm{warp}/M_\mathrm{warp}`. Note that this picture is only valid if the BH spin vector does not change much during this process (in both magnitude and direction), which can be ensured with wisely chosen time steps.
+
+Deciding whether accretion is prograde or retrograde
+----------------------------------------------------
+
+We now discuss how to decide whether the sign of spin is positive or negative. In the process of communicating with the inner disk through Lense-Thirring torques, the disk either aligns or counteraligns with the BH spin vector. The condition for which of the two occurs can be derived by assuming that the magnitude of the spin does not change during this alignment (`King et al. 2005 <https://ui.adsabs.harvard.edu/abs/2005MNRAS.363...49K/abstract>`_). Accretion is retrograde if
+
+.. math::
+    \cos \theta<-\frac{J_{\mathrm{warp}}}{2 J_{\mathrm{BH}}},
+    
+where :math:`\cos \theta=\hat{J}_\mathrm{BH}\cdot\hat{J}_\mathrm{d}` is the angle between the initial spin vector and the large-scale angular momentum of the disk. If this condition is not fulfilled, accretion is assumed to be prograde. Note that retrograde accretion is only possible if the angle between the spin vector and the large-scale accretion disk is larger than :math:`90^\circ`, and if the warp angular momentum is comparable to the BH one.
+
+Structure of the warped and precessing accretion disk
+-----------------------------------------------------
+
+As mentioned already, Lense-Thirring torques have different effects depending on the type of accretion disk. In particular, their effects depend on the ratio of the viscosity parameter :math:`\alpha` and the aspect ratio :math:`H/R`. For thin discs (:math:`\alpha\gg H/R`), the disk is exactly warped as in the manner described in the preceeding two sections (`Bardeen & Peterson 1975 <https://ui.adsabs.harvard.edu/abs/1975ApJ...195L..65B/abstract>`_). The radius :math:`R_\mathrm{warp}` which separates the inner and outer accretion disc can be calculated by equating the Lense-Thirring precession time-scale (:math:`t_\mathrm{p}=2\pi/\Omega_\mathrm{p}`, with :math:`\Omega_\mathrm{p}=2GJ_\mathrm{BH}/c^2R^3` the precession rate) and the vertical warp propagation time-scale (:math:`t_\mathrm{warp}=R^2/\nu_2`, with :math:`\nu_2` the kinematic viscosity in the vertical direction) (e.g. `Martin et al. 2007 <https://ui.adsabs.harvard.edu/abs/2007MNRAS.381.1617M/abstract>`_). The vertical kinematic viscosity :math:`\nu_2` can be related to the horizontal one, :math:`\nu_1`, by :math:`\nu_2=\xi\nu_1`, with :math:`\xi` a numerical parameter given by
+
+.. math::
+    \xi=\frac{2}{\alpha^2}\frac{1+7\alpha^2}{4+\alpha^2}
+
+(`Ogilvie 1999 <https://ui.adsabs.harvard.edu/abs/1999MNRAS.304..557O/abstract>`_, see also `Lodato et al. 2010 <https://ui.adsabs.harvard.edu/abs/2010MNRAS.405.1212L/abstract>`_ for a detailed discussion). We use the relation :math:`\dot{M}=3\pi\nu_1 \Sigma` to calculate :math:`\nu_1`, and therefore :math:`\nu_2`. The warp radius will depend on which region of the thin disc we assume, with each having its own expression for :math:`\Sigma`. In region b) of the `Shakura & Sunyaev (1973) <https://ui.adsabs.harvard.edu/abs/1973A%26A....24..337S/abstract>`_ thin disk, the surface density can be expressed as
+
+.. math::
+    \Sigma_\mathrm{TD,b}=6.84 \times 10^{5} \mathrm{~g} \mathrm{~cm}^{-2} \alpha^{-4 / 5} f_\mathrm{Edd}^{3 / 5}\left(\frac{M_{\mathrm{BH}}}{10^{8} M_{\odot}}\right)^{1 / 8}\left(\frac{R}{R_{\mathrm{S}}}\right)^{-3 / 5},
+    
+while in region c) we have
+
+.. math::
+    \Sigma_\mathrm{TD,c}=3.41 \times 10^{4} \mathrm{~g} \mathrm{~cm}^{-2} \alpha^{-4 / 5} f_\mathrm{Edd}^{7/10}\left(\frac{M_{\mathrm{BH}}}{10^{8} M_{\odot}}\right)^{1 / 20}\left(\frac{R}{R_{\mathrm{S}}}\right)^{-3 / 4}.
+    
+These relations lead to the following expressions for :math:`R_\mathrm{warp}`:
+
+.. math::
+    R_{\text {warp,TD,b}}=3410 R_{S} a^{5 / 8} \xi^{-5/8}\alpha^{-1 / 2} f_\mathrm{Edd}^{-1 / 4}\left(\frac{M_{\mathrm{BH}}}{10^{8} M_{\odot}}\right)^{1 / 8}
+    
+(in region b) and
+
+.. math::
+    R_\mathrm{warp,TD,c}=2629R_\mathrm{S}a^{4/7}\xi^{-4/7}\alpha^{-16/35}f_\mathrm{Edd}^{-6/35}\bigg(\frac{M_\mathrm{BH}}{10^8\hspace{0.5mm}\mathrm{M}_\odot}  \bigg)^{4/35},
+    
+(in region c), with :math:`R_\mathrm{S}=2R_\mathrm{G}` the Schwarzschild radius. These warp radii are generally of order :math:`\approx1000R_\mathrm{G}`, which can lead to fairly quick alignment of the thin disk with the large-scale angular momentum direction (quicker than any significant evolution in mass or spin magnitude, illustrating why the inclusion of the effects of Lense-Thirring torques is important).
+
+In the context of thin disks, there is a futher complication. The self-gravity of the disk may become important at large radii (see `Lodato 2007 <https://www.sif.it/riviste/sif/ncr/econtents/2007/030/07/article/0>`_ for a review). The disk will fragment in the region where the Toomre parameter is :math:`Q(R)>1`. We thus assume that the disk extends out to where :math:`Q(R_\mathrm{sg})=1`. The self-gravity radius :math:`R_\mathrm{sg}` can be calculated from this condition and the definition of the Toomre parameter :math:`Q=\Omega c_{\mathrm{s}} /(\pi G \Sigma)`, yielding
+
+.. math::
+    R_{\text {sg,TD,b}}=6460 R_{S} \alpha^{28/51} f_\mathrm{Edd}^{-18/51}\left(\frac{M_{\mathrm{BH}}}{10^{8} M_{\odot}}\right)^{-49/51}
+    
+in region b) and
+
+.. math::
+    R_\mathrm{sg,TD,c}=2456 R_{S} \alpha^{28/45} f_\mathrm{Edd}^{-22/45}\left(\frac{M_{\mathrm{BH}}}{10^{8} M_{\odot}}\right)^{-52/45}
+    
+in region c). In all our calculations involving :math:`R_\mathrm{warp}` (for deciding the sign of spin and evolving the direction of angular momentum, as described in the preceeding sections), we always take the minimum of :math:`R_\mathrm{warp}` and :math:`R_\mathrm{sg}`. This is because if :math:`R_\mathrm{sg}<R_\mathrm{warp}`, the entire disk of extent :math:`R_\mathrm{sg}` will be warped.
+
+The thick disk does not experience the Bardeen-Peterson effect, i.e. it is never truly aligned nor counter-aligned in its inner regions. Instead, the disk precesses out to several dozen :math:`R_\mathrm{G}`, as seen in simulations (e.g. `Fragile et al. 2007 <https://ui.adsabs.harvard.edu/abs/2007ApJ...668..417F/abstract>`_), and likely observations through quasi-periodic oscillations (QPO; e.g. `Ingram et al. 2012 <https://ui.adsabs.harvard.edu/abs/2012MNRAS.419.2369I/abstract>`_). The slim disk has received much less attention in both simulations and observations (it is both harder to simulate and observe), but its similarity to the thick disk in its geometric aspects likely means that it precesses in a similar manner.
+
+The exact behaviour of the thick and slim disk (which we will collectively call the advection-dominated disks) again depends on the ratio of :math:`\alpha` and :math:`H/R`. Unfortunately, the advection-dominated disks both satisfy :math:`\alpha\approx H/R`, and in this regime, the effects of Lense-Thirring torques are not well understood from a theoretical perspective. However, if :math:`\alpha\ll H/R` (the so-called bending-wave regime), Lense-Thirring torques are known to cause precession of the entire inner disk as a solid body, as seen in observations and simulations. For simplicity, we will thus assume this to be the case for advection-dominated disks.
+
+`Lubow et al. (2002) <https://ui.adsabs.harvard.edu/abs/2002MNRAS.337..706L/abstract>`_ studied the bending-wave regime. In the inner regions, the disk precesses around the spin axis, while in the outer regions, it is aligned with the large-scale angular momentum of the disk. Based on their results the transition radius between the precessing and non-precessing regions of the disk given by
+
+.. math::
+    R_\mathrm{warp,adv}=R_\mathrm{G}\bigg(\frac{384a}{25(H/R)^2}\bigg)^{2/5}.
+    
+In our model, we assume that the inner regions of the disks are on
+average aligned or counteraligned with the spin vector (one can think
+of this as averaging over the precession, which has periods of
+:math:`\approx` days, over long enough time scales). For simplicity, we  also refer to the radii within which this is true as the warp radii. For both of the advection-dominated disks, these radii are only of order several :math:`R_\mathrm{G}`. Note that similar values are found if one assumes that the Bardeen-Peterson effect operates in these disks. While there are some uncertainties in the assumptions we have made, we point out that using any of these values is much more physically motivated than using thin disk equations (the warp radii of order thousands of :math:`R_\mathrm{G}`), which is what is often done (e.g. `Griffin et al. 2019a <https://ui.adsabs.harvard.edu/abs/2019MNRAS.487..198G/abstract>`_, `Dubois et al. 2012 <https://ui.adsabs.harvard.edu/abs/2014MNRAS.440.1590D/abstract>`_).
+
+In order to determine the sign of spin and evolve the angular momentum direction, expressions for the warp mass :math:`M_\mathrm{warp}` and warp angular momentum :math:`J_\mathrm{warp}` are also needed. We calculate this using surface integrals as
+
+.. math::
+    M_\mathrm{warp}(R_\mathrm{warp})=2\pi\int_0^{R_\mathrm{warp}}\Sigma(R)R\mathrm{d}R,
+    
+and
+    
+.. math::
+    J_\mathrm{warp}(R_\mathrm{warp})=2\pi\int_0^{R_\mathrm{warp}}L(R)\Sigma(R)R\mathrm{d}R,
+    
+respectively. Here, :math:`L(R)` is the specific angular momentum. In the case of the thin disk, we assume Keplerian orbits, i.e. :math:`L(R)=\sqrt{M_\mathrm{BH}G R}`. For the advection-dominated disks, we assume that they are smaller by a numerical factor :math:`\Omega_0`, which is given in the self-similar solutions for the thick (`Narayan & Yi 1995b <https://ui.adsabs.harvard.edu/abs/1995ApJ...452..710N/abstract>`_) and slim disk (`Wang & Zhou 1999 <https://ui.adsabs.harvard.edu/abs/1999ApJ...516..420W/abstract>`_), seperately. The surface densities in both of these accretion disks are given by the same formula in the self-similar solutions, which is
+
+.. math::
+    \Sigma_\mathrm{adv}=\frac{\dot{M}}{2\pi R\vert v_\mathrm{r} \vert},
+
+where :math:`v_\mathrm{r}=-\alpha v_0 v_\mathrm{K}` is the radial velocity. Here, :math:`v_\mathrm{K}=\sqrt{M_\mathrm{BH}G/R}` is the Keplerian velocity, and :math:`v_0` is another numerical coefficient that differs between the two solutions. In the thick disk, the numerical coefficients are given by :math:`v_0=3/(5+2\varepsilon)` and :math:`\Omega_0=\sqrt{2\varepsilon/(5+2\varepsilon)}`, where :math:`\varepsilon=(5/3-\gamma)/(\gamma-1)`. The adiabatic index depends on how magnetized the disk is. In particular, it depends on the gas-to-total pressure ratio as :math:`\gamma = (8-3\beta)/(6-3\beta)`, and :math:`\beta` itself depends on :math:`\alpha` (see discussion above on radiative efficiency in the thin disk). :math:`v_0` varies weakly with :math:`\alpha`; for :math:`\alpha=0.05`, it is :math:`0.56`, whereas for :math:`\alpha=0.3`, it evaluates to 0.5. :math:`\Omega_0` depends on :math:`\alpha` somewhat more strongly; we obtain :math:`0.27` and :math:`0.41` for the same values of :math:`\alpha`. The latter value agrees well with the ratio of actual to Keplerian (ISCO) orbital velocity at the event horizon, which is :math:`0.45`. For the slim disc, :math:`v_0=\Omega_0=1/\sqrt{\gamma}`, with :math:`\gamma=5`.
+
+Black hole mergers
+------------------
+
+In the process of merging, BHs interact in a very complicated manner. Their final spin is not trivial to predict, and it can depend on a very large parameter space (including the mass ratio of the black holes and the relative orientation and magnitude of the spins). Orbital angular momentum plays a role in the merger as well. We use the fitting function found by `Rezzolla et al. (2009) <https://ui.adsabs.harvard.edu/abs/2009CQGra..26i4023R/abstract>`_, whose results have been found to be very accurate in newer and more sophisticated studies that sweep the huge parameter space of possible merger configurations. These formulas are also applicable to cosmological simulations, since they cover the scenario of inspiral from very large distances.
+
+The final spin, according to `Rezzolla et al. (2009) <https://ui.adsabs.harvard.edu/abs/2009CQGra..26i4023R/abstract>`_ can be calculated as
+
+.. math::
+    \mathbf{a}_\mathrm{fin} = \frac{1}{(1+q)^2}(\mathbf{a}_1+\mathbf{a}_2q^2+\mathbf{l}q),
+
+where :math:`q=M_2/M_1` is the mass ratio (such that :math:`M_2<M_1`), :math:`\mathbf{a}_1` and :math:`\mathbf{a}_2` are the spin vectors, and :math:`\mathbf{l}` is a vector whose direction is the same as that of the orbital angular momentum :math:`\mathbf{L}` (in the centre-of-mass frame), while its magnitude is given by
+
+.. math::
+    |\mathbf{l}|=\frac{s_{4}}{\left(1+q^{2}\right)^{2}}\left(\left|\mathbf{a}_{1}\right|^{2}+|\mathbf{a}|_{1}^{2} q^{4}+2\left|\mathbf{a}_{1} \| \mathbf{a}_{2}\right| q^{2} \cos \phi\right)+ \\
+    \left(\frac{s_{5} \mu+t_{0}+2}{1+q^{2}}\right)\left(\left|\mathbf{a}_{1}\right| \cos \theta+\left|\mathbf{a}_{2}\right| q^{2} \cos \xi\right)+ \\
+    2 \sqrt{3}+t_{2} \mu+t_{3} \mu^{2}.
+
+Here, :math:`\mu=q/(1+q)^2` is the symmetric mass ratio, and :math:`s_4 = -0.1229`, :math:`s_5 = -0.4537`, :math:`t_0 = -2.8904`, :math:`t_2 = -3.5171`, :math:`t_3 = 2.5763`. The three cosines depend on the angles between the different vectors which play a role in the merger: :math:`\cos \phi=\hat{\mathbf{a}_{1}} \cdot \hat{\mathbf{a}_{\mathbf{2}}}`, :math:`\cos \theta=\hat{\mathbf{a}_{1}} \cdot \hat{\mathbf{l}}`, :math:`\cos \xi=\hat{\mathbf{a}_{2}} \cdot \hat{\mathbf{l}}`.
+
+Given the information available within the model, we could in principle calculate the recoil velocity of the remnant, as well as the total mass fraction lost to gravitational waves. We do not implement the former at this stage since we cannot reliably track the movement of black holes in their host galaxies. However, we do implement the latter. We use results from the same series of numerical relativity simulations as above (`Barausse et al. 2012 <https://ui.adsabs.harvard.edu/abs/2012ApJ...758...63B/abstract>`_) and write the final mass of the remnant as:
+
+.. math::
+    M_\mathrm{BH,fin} = (M_\mathrm{BH,1}+M_\mathrm{BH,2})\Big\{1 - [1 - e_\mathrm{ISCO}(\tilde{a})]\mu  - 4\mu^2[4p_0+16p_1\tilde{a}(\tilde{a}+1)+e_\mathrm{ISCO}(\tilde{a})-1]\Big\},
+
+where :math:`p_0=0.04827`, :math:`p_1=0.01707` and :math:`e_\mathrm{ISCO}(\tilde{a})` is the dimensionless specific binding energy at the innermost stable circular orbit calculated using an effective spin variable defined as 
+
+.. math::
+    \tilde{a} = \frac{|\mathbf{a_1}|\cos\theta+|\mathbf{a_2}|\cos\xi}{(1+q)^2}.