Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
  • 840-unit-test-testtimeline-fails
  • 875-wendland-c6-missing-neighbour-contributions
  • 887-code-does-not-compile-with-parmetis-installed-locally-but-without-metis
  • CubeTest
  • FS_Del
  • GEARRT_Iliev1
  • GEARRT_Iliev3
  • GEARRT_Iliev4
  • GEARRT_Iliev5
  • GEARRT_Iliev5-fixed-nr-subcycles
  • GEARRT_Iliev7
  • GEARRT_Iliev_static
  • GEARRT_Ivanova
  • GEARRT_fixed_nr_subcycles
  • GEARRT_injection_tests_Iliev0
  • GPU_swift
  • GrackleCoolingUpdates2
  • Lambda-T-table
  • MAGMA2
  • MAGMA2_matthieu
  • MHD_FS
  • MHD_FS_TESTs
  • MHD_FS_VP_AdvectGauge
  • MHD_Orestis
  • MHD_canvas
  • MHD_canvas_RF_128
  • MHD_canvas_RF_growth_rate
  • MHD_canvas_RobertsFlow
  • MHD_canvas_SPH_errors
  • MHD_canvas_matthieu
  • MHD_canvas_nickishch
  • MHD_canvas_nickishch_Lorentz_force_test
  • MHD_canvas_nickishch_track_everything
  • MHD_canvas_sid
  • OAK/CPAW_updates
  • OAK/LoopAdvectionTest
  • OAK/adaptive_divv
  • OAK/kinetic_dedner
  • REMIX_cosmo
  • RT_dualc
  • RT_recombination_radiation
  • RT_test_mladen
  • SIDM
  • SIDM_wKDSDK
  • SNdust
  • SPHM1RT_CosmologicalStromgrenSphere
  • SPHM1RT_bincheck
  • SPHM1RT_smoothedRT
  • TangoSIDM
  • TestPropagation3D
  • Test_fixedhProb
  • activate_fewer_comms
  • active_h_max_optimization
  • adaptive_softening_Lieuwe
  • add_2p5D
  • add_black_holes_checks
  • adding_sidm_to_master
  • agn_crksph
  • agn_crksph_subtask_speedup
  • amd-optimization
  • arm_vec
  • automatic_tasks
  • better_ray_RNG
  • black_holes_accreted_angular_momenta_from_gas
  • burkert-potential
  • c11
  • c11_atomics_copy
  • cancel_all_sorts
  • cell_exchange_improvements
  • cell_types
  • cherry-pick-cd1c39e0
  • comm_tasks_are_special
  • conduction_velocities
  • cpp-fixes
  • cuda_test
  • darwin/adaptive_softening
  • darwin/gear_chemistry_fluxes
  • darwin/gear_mechanical_feedback
  • darwin/gear_preSN_feedback
  • darwin/gear_radiation
  • darwin/simulations
  • darwin/sink_formation_proba
  • darwin/sink_mpi
  • darwin/sink_mpi_physics
  • dead-time-stats
  • derijcke_cooling
  • dev_cms
  • do-not-activate-empty-star-pairs
  • domain_zoom_nometis
  • drift_flag_debug_check
  • driven_turbulence
  • driven_turbulence_forcings
  • engineering
  • eos_updates
  • evrard_disc
  • expand_fof_2022
  • explict_bkg_cdim
  • fewer_gpart_comms
  • fewer_star_comms
  • fewer_timestep_comms_no_empty_pairs
  • v0.0
  • v0.1
  • v0.1.0-pre
  • v0.2.0
  • v0.3.0
  • v0.4.0
  • v0.5.0
  • v0.6.0
  • v0.7.0
  • v0.8.0
  • v0.8.1
  • v0.8.2
  • v0.8.3
  • v0.8.4
  • v0.8.5
  • v0.9.0
  • v1.0.0
  • v2025.01
  • v2025.04
119 results

Target

Select target project
  • dc-oman1/swiftsim
  • swift/swiftsim
  • pdraper/swiftsim
  • tkchan/swiftsim
  • dc-turn5/swiftsim
5 results
Select Git revision
  • 840-unit-test-testtimeline-fails
  • 875-wendland-c6-missing-neighbour-contributions
  • CubeTest
  • FS_Del
  • GEARRT_Iliev1
  • GEARRT_Iliev3
  • GEARRT_Iliev4
  • GEARRT_Iliev5
  • GEARRT_Iliev5-fixed-nr-subcycles
  • GEARRT_Iliev7
  • GEARRT_Iliev_static
  • GEARRT_Ivanova
  • GEARRT_Stan_project_cosmology
  • GEARRT_cosmo_IonMassFraction_example
  • GEARRT_cosmo_redshifting
  • GEARRT_cosmo_subcycling_Stan
  • GEARRT_cosmo_thermochem
  • GEARRT_fixed_nr_subcycles
  • GEARRT_injection_tests_Iliev0
  • GPU_swift
  • GrackleCoolingUpdates2
  • Lambda-T-table
  • MAGMA2
  • MAGMA2_matthieu
  • MHD_FS
  • MHD_FS_TESTs
  • MHD_FS_VP_AdvectGauge
  • MHD_Orestis
  • MHD_canvas
  • MHD_canvas_nickishch
  • MHD_canvas_sid
  • OAK/CPAW_updates
  • OAK/LoopAdvectionTest
  • OAK/kinetic_dedner
  • RT_dualc
  • RT_recombination_radiation
  • RT_test_mladen
  • SIDM
  • SIDM_wKDSDK
  • SNdust
  • SPHM1RT_CosmologicalStromgrenSphere
  • SPHM1RT_bincheck
  • SPHM1RT_smoothedRT
  • TangoSIDM
  • TestPropagation3D
  • Test_fixedhProb
  • active_h_max_optimization
  • adaptive_softening_Lieuwe
  • add_black_holes_checks
  • adding_sidm_to_master
  • agn_crksph
  • agn_crksph_subtask_speedup
  • amd-optimization
  • arm_vec
  • automatic_tasks
  • better_ray_RNG
  • black_holes_accreted_angular_momenta_from_gas
  • burkert-potential
  • c11
  • c11_atomics_copy
  • cell_types
  • cherry-pick-cd1c39e0
  • comm_tasks_are_special
  • conduction_velocities
  • cuda_test
  • dead-time-stats
  • derijcke_cooling
  • dev_cms
  • do-not-activate-empty-star-pairs
  • domain_zoom_nometis
  • drift_flag_debug_check
  • driven_turbulence
  • engineering
  • eos_updates
  • evrard_disc
  • expand_fof
  • expand_fof_2022
  • explict_bkg_cdim
  • fewer_timestep_comms_no_empty_pairs
  • fix-velcut
  • fix_max_toplevel_cell_rounding
  • fixed-bh-accretion
  • fixed_hSIDM
  • flux-counter
  • for_isabel
  • foreign_gpart
  • format_sh_eagle_stars
  • fstasys/Clean_Blast_now_with_VP
  • fstasys/Clean_Fast_Rotor_now_with_VP
  • fstasys/MHD_cosmo_run
  • fstasys/RobertsFlow_plain_forcing
  • fstasys/VP_CosmoRun.GalileanInvariance
  • fstasys/cleanout_gauge_effects_in_VP
  • gear_sink_imf_sampling
  • gear_sink_imf_sampling_merged
  • gear_sink_regulated_accretion
  • genetic_partitioning2
  • gizmo
  • gizmo-new-timestep
  • gizmo-timestep
  • v0.0
  • v0.1
  • v0.1.0-pre
  • v0.2.0
  • v0.3.0
  • v0.4.0
  • v0.5.0
  • v0.6.0
  • v0.7.0
  • v0.8.0
  • v0.8.1
  • v0.8.2
  • v0.8.3
  • v0.8.4
  • v0.8.5
  • v0.9.0
  • v1.0.0
117 results
Show changes
Showing
with 1060 additions and 661 deletions
......@@ -34,3 +34,8 @@ In order to add a new scheme, you will need to:
``nobase_noinst_HEADERS``, add your new header files.
6. Update the documentation. Add your equations/documentation to ``doc/RTD``.
.. toctree::
:caption: Table of Contents
sink_adding_new_scheme
.. Adding new schemes
Darwin Roduit, 16 Ocotber 2024
.. _new_option_sink:
How to add your sink scheme
-------------------------------
Here, we provide comprehensive information to guide you in adding your sink scheme into Swift. To better understand how to add new schemes within Swift, read the general information provided on :ref:`new_option` page.
The default sink scheme is empty and gives you an idea of the minimum required fields and functions for the code to compile. The GEAR sink module has the base functions plus some extra ones for its operations. It can provide you with a working example. However, it can only work with the GEAR feedback module since it relies on IMF properties that are only located there.
As discussed in the GEAR sink :ref:`sink_GEAR_model_summary`, the physics relies on the following tasks: sink formation, gas and sink particle flagging, gas swallowing, sink swallowing and star formation. You do not need to care about the tasks, only the core functions within the sink module. However, you may need to locate where the code calls these functions. The file ``src/runner_others.c`` contains the ``runner_do_star_formation_sink()`` and ``runner_do_sink_formation()``. These functions are responsible for generating stars out of sinks and sinks from gas particles. The other general task-related functions are in ``src/runner_sinks.c``.
The following presents the most essential functions you need to implement. This will give you an idea of the workload.
Sink formation
~~~~~~~~~~~~~~
Before forming a sink, the code loops over all gas and sink particles to gather data about its neighbourhood. This is performed in ``sink_prepare_part_sink_formation_gas_criteria()`` and ``sink_prepare_part_sink_formation_sink_criteria()`` functions. For instance, in GEAR, we compute the total potential energy, thermal energy, etc.
Then, to decide if we can turn a gas particle into a sink particle, the function ``sink_is_forming()`` is called. Before forming a sink particle, there is a call to ``sink_should_convert_to_sink()``. This function determines whether the gas particle must transform into a sink. Both functions return either 0 or 1.
Gas-sink density interactions
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The first interaction task to be run for the sinks is the density task. This task updates the smoothing length for the sink particle, unless a fixed cutoff radius is being used (coming soon). It can also calculate the contributions made by neigbouring gas particles to the density, sound speed, velocity etc. at the location of the sink. Code for these interactions should be added to ``sink_iact.h/runner_iact_nonsym_sinks_gas_density()``.
Once the contributions of all neigbouring gas particles have been calculated, the density calculation is completed by the sink density ghost task. You can set what this task does with the functions ``sink_end_density()`` and ``sink_prepare_swallow()`` in ``sink.h``.
The ``sink_end_density()`` function completes the calculation of the smoothing length (coming soon), and this is where you can finish density-based calculations by e.g. dividing mass-weighted contributions to the velocity field by the total density in the kernel. For examples of this, see the equivalent task for the black hole particles.
The ``sink_prepare_swallow()`` task is where you can calculate density-based quantities that you might need to use in swallowing interactions later. For example, a Bondi-Hoyle accretion prescription should calculate an accretion rate and target mass to be accreted here.
Gas and sink flagging: finding whom to eat
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Before accreting the gas/sink particles, the sink needs to look for eligible particles. The gas swallow interactions are performed within ``runner_iact_nonsym_sinks_gas_swallow()`` and the sink swallow in ``runner_iact_nonsym_sinks_sink_swallow()``.
Gas and sink swallowing: updating the sink properties
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
When the sink swallows gas particles, it updates its internal properties based on the gas particles' properties. The ``sink_swallow_part()`` function takes care of this.
Similarly, when the sink swallows sink particles, it updates its properties from the to-be-swallowed sink particles. The ``sink_swallow_sink()`` performs the update.
There is no more than that. The code properly removes the swallowed particles.
Star formation
~~~~~~~~~~~~~~
The most important function is ``sink_spawn_star()``. It controls whether the code should continue creating stars and must return 0 or 1.
The ``sink_copy_properties_to_star()`` does what it says. This function is also responsible for adequately initialising the stars' properties. In GEAR, we give new positions and velocities within this function.
The following three functions allow you to update your sink particles (e.g. their masses) before, during and after the star formation loop: ``sink_update_sink_properties_before_star_formation()``, ``sink_update_sink_properties_during_star_formation()`` and ``sink_update_sink_properties_after_star_formation()``.
These functions are located in ``sink/Default/sink.h``
......@@ -15,11 +15,12 @@ instead of, the lossless gzip compression filter.
**These compression filters are lossy, meaning that they modify the
data written to disk**
*The filters will reduce the accuracy of the data stored. No check is
made inside SWIFT to verify that the applied filters make sense. Poor
choices can lead to all the values of a given array reduced to 0, Inf,
or to have lost too much accuracy to be useful. The onus is entirely
on the user to choose wisely how they want to compress their data.*
.. warning::
The filters will reduce the accuracy of the data stored. No check is
made inside SWIFT to verify that the applied filters make sense. Poor
choices can lead to all the values of a given array reduced to 0, Inf,
or to have lost too much accuracy to be useful. The onus is entirely
on the user to choose wisely how they want to compress their data.
The filters are not applied when using parallel-hdf5.
......@@ -27,6 +28,14 @@ The name of any filter applied is carried by each individual field in
the snapshot using the meta-data attribute ``Lossy compression
filter``.
.. warning::
Starting with HDF5 version 1.14.4, filters which compress the data
by more than 2x are flagged as problematic (see their
`doc <https://docs.hdfgroup.org/hdf5/v1_14/group___f_a_p_l.html#gafa8e677af3200e155e9208522f8e05c0>`_
). SWIFT can nevertheless write files with them by setting the
appropriate file-level flags. However, some tools (such as
``h5py``) may *not* be able to read these fields.
The available filters are listed below.
N-bit filters for long long integers
......@@ -171,6 +180,8 @@ Same for a ``double`` (64 bits) output:
+=================+==============+==============+=============+===================================================+===================+
| No filter | 52 | 11 | 15.9 digits | :math:`[2.2\times 10^{-308}, 1.8\times 10^{308}]` | --- |
+-----------------+--------------+--------------+-------------+---------------------------------------------------+-------------------+
| ``DMantissa21`` | 21 | 11 | 6.62 digits | :math:`[2.2\times 10^{-308}, 1.8\times 10^{308}]` | 1.93x |
+-----------------+--------------+--------------+-------------+---------------------------------------------------+-------------------+
| ``DMantissa13`` | 13 | 11 | 4.21 digits | :math:`[2.2\times 10^{-308}, 1.8\times 10^{308}]` | 2.56x |
+-----------------+--------------+--------------+-------------+---------------------------------------------------+-------------------+
| ``DMantissa9`` | 9 | 11 | 3.01 digits | :math:`[2.2\times 10^{-308}, 1.8\times 10^{308}]` | 3.05x |
......
......@@ -74,6 +74,9 @@ CGS. Entries in the file look like:
SmoothingLengths_Gas: on # Co-moving smoothing lengths (FWHM of the kernel) of the particles : a U_L [ cm ]
...
This can also be used to set the outputs produced by the
:ref:`fof_stand_alone_label`.
For cosmological simulations, users can optionally add the ``--cosmology`` flag
to generate the field names appropriate for such a run.
......
......@@ -466,7 +466,10 @@ can be either drawn randomly by setting the parameter ``generate_random_ids``
newly generated IDs do not clash with any other pre-existing particle. If this
option is set to :math:`0` (the default setting) then the new IDs are created in
increasing order from the maximal pre-existing value in the simulation, hence
preventing any clash.
preventing any clash. Finally, if the option
``particle_splitting_log_extra_splits`` is set, the code will log all the splits
that go beyond the maximal allowed (typically 64) in a file so that the split tree
for these particles can still be reconstructed.
The final set of parameters in this section determine the initial and minimum
temperatures of the particles.
......@@ -503,7 +506,7 @@ The full section to start a typical cosmological run would be:
minimal_temperature: 100 # U_T
H_mass_fraction: 0.755
H_ionization_temperature: 1e4 # U_T
particle_splitting: 1
particle_splitting: 1
particle_splitting_mass_threshold: 5e-3 # U_M
.. _Parameters_Stars:
......@@ -514,7 +517,7 @@ Stars
The ``Stars`` section is used to set parameters that describe the Stars
calculations when doing feedback or enrichment. Note that if stars only act
gravitationally (i.e. SWIFT is run *without* ``--feedback``) no parameters
in this section are used.
in this section are used.
The first four parameters are related to the neighbour search:
......@@ -576,6 +579,26 @@ specified, SWIFT will start and use the birth times specified in the
ICs. If no values are given in the ICs, the stars' birth times will be
zeroed, which can cause issues depending on the type of run performed.
.. _Parameters_Sinks:
Sinks
-----
Currently, there are two models for the sink particles, the Default model and the GEAR one. Their parameters are described below. To choose a model, configure the code with ``--with-sink=<model>``, where ``<model>`` can be ``none`` or ``GEAR``. To run with sink particles, add the option ``--sinks``.
Below you will find the description of the ``none`` which is the default model. For ``GEAR`` model, please refer to :ref:`sink_GEAR_model`.
By default, the code is configured with ``--with-sink=none``. Then, the ``DefaultSink`` section is used to set parameters that describe the sinks in this model. The unique parameter is the sink accretion radius (also called cut-off radius): ``cut_off_radius``.
Note that this model does not create sink particles or accrete gas.
The full section is:
.. code:: YAML
DefaultSink:
cut_off_radius: 1e-3 # Cut off radius of the sink particles (in internal units). This parameter should be adapted with the resolution..
.. _Parameters_time_integration:
Time Integration
......@@ -623,7 +646,7 @@ the start and end times or scale factors from the parameter file.
``max_dt_RMS_factor`` (default: ``0.25``)
* Whether or not only the gas particle masses should be considered for
the baryon component of the calculation: ``dt_RMS_use_gas_only`` (default: ``0``)
These values rarely need altering. The second parameter is only
meaningful if a subgrid model produces star (or other) particles with
masses substantially smaller than the gas masses. See the theory
......@@ -648,7 +671,7 @@ Whilst for a cosmological run, one would need:
dt_min: 1e-10
max_dt_RMS_factor: 0.25 # Default optional value
dt_RMS_use_gas_only: 0 # Default optional value
.. _Parameters_ICs:
Initial Conditions
......@@ -693,7 +716,7 @@ Finally, SWIFT also offers these options:
* Whether to replicate the box along each axis: ``replicate`` (default: ``1``).
* Whether to re-map the IDs to the range ``[0, N]`` and hence discard
the original IDs from the IC file: ``remap_ids`` (default: ``0``).
The shift is expressed in internal units and will be written to the header of
the snapshots. The option to replicate the box is especially useful for
weak-scaling tests. When set to an integer >1, the box size is multiplied by
......@@ -858,24 +881,42 @@ that can be used as if it was a non-distributed snapshot. In this case, the
HDF5 library itself can figure out which file is needed when manipulating the
snapshot.
On Lustre filesystems [#f4]_ it is important to properly stripe files to achieve
a good writing speed. If the parameter ``lustre_OST_count`` is set to the number
of OSTs present on the system, then SWIFT will set the `stripe count` of each
distributed file to `1` and set each file's `stripe index` to the MPI rank
generating it modulo the OST count [#f5]_. If the parameter is not set then the
files will be created with the default system policy (or whatever was set for
the directory where the files are written). This parameter has no effect on
non-Lustre file systems and no effect if distributed snapshots are not used.
* The number of Lustre OSTs to distribute the single-striped distributed
snapshot files over: ``lustre_OST_count`` (default: ``0``)
On Lustre filesystems [#f4]_ it is important to properly stripe files to
achieve a good writing and reading speed. If the parameter
``lustre_OST_checks`` is set and the lustre API is available SWIFT will
determine the number of OSTs available and rank these by free space, it will
then set the `stripe count` of each file to `1` and choose an OST
`offset` so each rank writes to a different OST, unless there are more ranks
than OSTs in which case the assignment wraps. In this way OSTs should be
filled evenly and written to using an optimal access pattern.
If the parameter is not set then the files will be created with the default
system policy (or whatever was set for the directory where the files are
written). This parameter has no effect on non-Lustre file systems.
Other parameters are also provided to handle the cases when individual OSTs do
not have sufficient free space to write a file: ``lustre_OST_free`` and
when OSTs are closed for administrative reasons: ``lustre_OST_test``, in which
case they cannot be written. This is important as the `offset` assignment in
this case is not used by lustre which picks the next writable OST, so in our
scheme such OSTs will be used more times than we intended.
* Use the lustre API to assign a stripe and offset to the distributed snapshot
files:
``lustre_OST_checks`` (default: ``0``)
* Do not use OSTs that do not have a certain amount of free space in MiB.
Zero disables and -1 activates a guess based on the size of the process:
``lustre_OST_free`` (default: ``0``)
* Check OSTs can be written to and remove those from consideration:
``lustre_OST_test`` (default: ``0``)
Users can optionally ask to randomly sub-sample the particles in the snapshots.
This is specified for each particle type individually:
* Whether to switch on sub-sampling: ``subsample``
* Whether to switch on sub-sampling: ``subsample_fraction``
* Whether to switch on sub-sampling: ``subsample``
* Whether to switch on sub-sampling: ``subsample_fraction``
These are arrays of 7 elements defaulting to seven 0s if left unspecified. Each
entry corresponds to the particle type used in the initial conditions and
......@@ -885,7 +926,7 @@ indicating the fraction of particles to keep in the outputs. Note that the
selection of particles is selected randomly for each individual
snapshot. Particles can hence not be traced back from output to output when this
is switched on.
Users can optionally specify the level of compression used by the HDF5 library
using the parameter:
......@@ -1103,7 +1144,7 @@ output field will be enabled. If this is 0 lossy compression is not applied.
* Whether to use lossless compression in the particle outputs: ``particles_gzip_level``
If this is non-zero the HDF5 deflate filter will be applied to lightcone particle output with
the compression level set to the specified value.
the compression level set to the specified value.
* HEALPix map resolution: ``nside``
......@@ -1144,7 +1185,7 @@ filter names. Set the filter name to ``on`` to disable compression.
* Whether to use lossless compression in the HEALPix map outputs: ``maps_gzip_level``
If this is non-zero the HDF5 deflate filter will be applied to the lightcone map output with
the compression level set to the specified value.
the compression level set to the specified value.
The following shows a full set of light cone parameters for the case where we're making two
light cones which only differ in the location of the observer:
......@@ -1158,7 +1199,7 @@ light cones which only differ in the location of the observer:
buffer_chunk_size: 100000
max_particles_buffered: 1000000
hdf5_chunk_size: 10000
# Redshift ranges for particle types
z_range_for_Gas: [0.0, 0.05]
z_range_for_DM: [0.0, 0.05]
......@@ -1166,7 +1207,7 @@ light cones which only differ in the location of the observer:
z_range_for_Stars: [0.0, 0.05]
z_range_for_BH: [0.0, 0.05]
z_range_for_Neutrino: [0.0, 0.05]
# Healpix map parameters
nside: 512
radius_file: ./shell_radii.txt
......@@ -1191,7 +1232,7 @@ light cones which only differ in the location of the observer:
enabled: 1
basename: lightcone1
observer_position: [74.2, 10.80, 53.59]
An example of the radius file::
......@@ -1217,11 +1258,11 @@ Equation of State (EoS)
The ``EoS`` section contains options for the equations of state.
Multiple EoS can be used for :ref:`planetary`,
see :ref:`planetary_eos` for more information.
see :ref:`planetary_eos` for more information.
To enable one or multiple EoS, the corresponding ``planetary_use_*:``
flag(s) must be set to ``1`` in the parameter file for a simulation,
along with the path to any table files, which are set by the
along with the path to any table files, which are set by the
``planetary_*_table_file:`` parameters.
For the (non-planetary) isothermal EoS, the ``isothermal_internal_energy:``
......@@ -1280,6 +1321,7 @@ The options are:
* The number of grid foldings to use: ``num_folds``.
* The factor by which to fold at each iteration: ``fold_factor`` (default: 4)
* The order of the window function: ``window_order`` (default: 3)
* Whether or not to correct the placement of the centre of the k-bins for small k values: ``shift_centre_small_k_bins`` (default: 1)
The window order sets the way the particle properties get assigned to the mesh.
Order 1 corresponds to the nearest-grid-point (NGP), order 2 to cloud-in-cell
......@@ -1319,6 +1361,20 @@ are mutually exclusive. The particles selected in each half are different in
each output. Note that neutrino PS can only be computed when neutrinos are
simulated using particles.
SWIFT uses bins of integer :math:`k`, with bins :math:`[0.5,1.5]`, :math:`[1.5,2.5]` etc. The
representative :math:`k` values used to be assigned to the bin centres (so k=1, 2, etc), which are
then transformed to physical :math:`k` by a factor :math:`kL/(2*pi)`. For the first few bins, only
few modes contribute to each bin. It is then advantageous to move the "centre" of the bin to the
actual location correponding to the mean of the contributing modes. The :math:`k` label of the bin
is thus shifted by a small amount. The way to calculate these shifts is to consider a 3D cube of
:math:`(kx,ky,kz)` cells and check which cells fall inside a spherical shell with boundaries
:math:`(i+0.5,i+1.5)`, then calculate the average :math:`k=\sqrt{kx^2+ky^2+kz^2}`. So for
:math:`i=0` there cells :math:`k=1` and 12 cells :math:`k=\sqrt(2)`, so the weighted k becomes
:math:`(6 * 1 + 12 * \sqrt(2)) / 18 = 1.2761424`. Note that only the first 7 (22) bins require a
correction larger than 1 (0.1) percent. We apply a correction to the first 128 terms. This
correction is activated when ``shift_centre_small_k_bins`` is switched on (that's the default
behaviour).
An example of a valid power-spectrum section of the parameter file looks like:
.. code:: YAML
......@@ -1392,17 +1448,35 @@ the MPI-rank. SWIFT writes one file per MPI rank. If the ``save`` option has
been activated, the previous set of restart files will be named
``basename_000000.rst.prev``.
On Lustre filesystems [#f4]_ it is important to properly stripe files to achieve
a good writing and reading speed. If the parameter ``lustre_OST_count`` is set
to the number of OSTs present on the system, then SWIFT will set the `stripe
count` of each restart file to `1` and set each file's `stripe index` to the MPI
rank generating it modulo the OST count [#f5]_. If the parameter is not set then
the files will be created with the default system policy (or whatever was set
for the directory where the files are written). This parameter has no effect on
non-Lustre file systems.
On Lustre filesystems [#f4]_ it is important to properly stripe files to
achieve a good writing and reading speed. If the parameter
``lustre_OST_checks`` is set and the lustre API is available SWIFT will
determine the number of OSTs available and rank these by free space, it will
then set the `stripe count` of each restart file to `1` and choose an OST
`offset` so each rank writes to a different OST, unless there are more ranks
than OSTs in which case the assignment wraps. In this way OSTs should be
filled evenly and written to using an optimal access pattern.
* The number of Lustre OSTs to distribute the single-striped restart files over:
``lustre_OST_count`` (default: ``0``)
If the parameter is not set then the files will be created with the default
system policy (or whatever was set for the directory where the files are
written). This parameter has no effect on non-Lustre file systems.
Other parameters are also provided to handle the cases when individual OSTs do
not have sufficient free space to write a restart file: ``lustre_OST_free`` and
when OSTs are closed for administrative reasons: ``lustre_OST_test``, in which
case they cannot be written. This is important as the `offset` assignment in
this case is not used by lustre which picks the next writable OST, so in our
scheme such OSTs will be used more times than we intended.
* Use the lustre API to assign a stripe and offset to restart files:
``lustre_OST_checks`` (default: ``0``)
* Do not use OSTs that do not have a certain amount of free space in MiB.
Zero disables and -1 activates a guess based on the size of the process:
``lustre_OST_free`` (default: ``0``)
* Check OSTs can be written to and remove those from consideration:
``lustre_OST_test`` (default: ``0``)
SWIFT can also be stopped by creating an empty file called ``stop`` in the
directory where the restart files are written (i.e. the directory speicified by
......@@ -1810,11 +1884,11 @@ necessary and one would use:
invoke_stf: 1 # We want VELOCIraptor to be called when snapshots are dumped.
# ...
# Rest of the snapshots properties
StructureFinding:
config_file_name: my_stf_configuration_file.cfg # See the VELOCIraptor manual for the content of this file.
basename: ./haloes/ # Write the catalogs in this sub-directory
If one additionally want to call VELOCIraptor at times not linked with
snapshots, the additional parameters need to be supplied.
......@@ -1894,7 +1968,7 @@ Fermi-Dirac momenta will be generated if ``generate_ics`` is used. The
``use_delta_f``. Finally, a random seed for the Fermi-Dirac momenta can
be set with ``neutrino_seed``.
For mode details on the neutrino implementation, refer to :ref:`Neutrinos`.
For mode details on the neutrino implementation, refer to :ref:`Neutrinos`.
A complete specification of the model looks like
.. code:: YAML
......@@ -1906,7 +1980,7 @@ A complete specification of the model looks like
------------------------
.. [#f1] The thorough reader (or overly keen SWIFT tester) would find that the speed of light is :math:`c=1.8026\times10^{12}\,\rm{fur}\,\rm{ftn}^{-1}`, Newton's constant becomes :math:`G_N=4.896735\times10^{-4}~\rm{fur}^3\,\rm{fir}^{-1}\,\rm{ftn}^{-2}` and Planck's constant turns into :math:`h=4.851453\times 10^{-34}~\rm{fur}^2\,\rm{fir}\,\rm{ftn}^{-1}`.
......
......@@ -13,6 +13,8 @@ Every SPH particle then requires and carries the additional ``MaterialID`` flag
from the initial conditions file. This flag indicates the particle's material
and which EoS it should use.
If you have another EoS that you would like us to add, then just let us know!
It is important to check that the EoS you use are appropriate
for the conditions in the simulation that you run.
Please follow the original sources of these EoS for more information and
......@@ -24,15 +26,18 @@ So far, we have implemented several Tillotson, ANEOS, SESAME,
and Hubbard \& MacFarlane (1980) materials, with more on the way.
Custom materials in SESAME-style tables can also be provided.
The material's ID is set by a somewhat arbitrary base type ID
(multiplied by 100) plus an individual value:
(multiplied by 100) plus an individual value, matching our code for making
planetary initial conditions, `WoMa <https://github.com/srbonilla/WoMa>`_:
+ Ideal gas: ``0``
+ Default (\\(\\gamma\\) set using ``--with-adiabatic-index``, default 5/3): ``0``
+ Default (Set :math:`\gamma` using ``--with-adiabatic-index``, default 5/3): ``0``
+ Tillotson (Melosh, 2007): ``1``
+ Iron: ``100``
+ Granite: ``101``
+ Water: ``102``
+ Basalt: ``103``
+ Ice: ``104``
+ Custom user-provided parameters: ``190``, ``191``, ..., ``199``
+ Hubbard \& MacFarlane (1980): ``2``
+ Hydrogen-helium atmosphere: ``200``
+ Ice H20-CH4-NH3 mix: ``201``
......@@ -42,6 +47,10 @@ The material's ID is set by a somewhat arbitrary base type ID
+ Basalt (7530): ``301``
+ Water (7154): ``302``
+ Senft \& Stewart (2008) water: ``303``
+ AQUA, Haldemann, J. et al. (2020) water: ``304``
+ Chabrier, G. et al. (2019) Hydrogen: ``305``
+ Chabrier, G. et al. (2019) Helium: ``306``
+ Chabrier & Debras (2021) H/He mixture Y=0.245 (Jupiter): ``307``
+ ANEOS (in SESAME-style tables): ``4``
+ Forsterite (Stewart et al. 2019): ``400``
+ Iron (Stewart, zenodo.org/record/3866507): ``401``
......@@ -53,7 +62,7 @@ The data files for the tabulated EoS can be downloaded using
the ``examples/Planetary/EoSTables/get_eos_tables.sh`` script.
To enable one or multiple EoS, the corresponding ``planetary_use_*:``
flag(s) must be set to ``1`` in the parameter file for a simulation,
flag(s) must be set from ``0`` to ``1`` in the parameter file for a simulation,
along with the path to any table files, which are set by the
``planetary_*_table_file:`` parameters,
as detailed in :ref:`Parameters_eos` and ``examples/parameter_example.yml``.
......@@ -61,18 +70,18 @@ as detailed in :ref:`Parameters_eos` and ``examples/parameter_example.yml``.
Unlike the EoS for an ideal or isothermal gas, these more complicated materials
do not always include transformations between the internal energy,
temperature, and entropy. At the moment, we have implemented
\\(P(\\rho, u)\\) and \\(c_s(\\rho, u)\\) (and more in some cases),
:math:`P(\rho, u)` and :math:`c_s(\rho, u)` (and more in some cases),
which is sufficient for the :ref:`planetary_sph` hydro scheme,
but some materials may thus currently be incompatible with
e.g. entropy-based schemes.
The Tillotson sound speed was derived using
\\(c_s^2 = \\left. ( \\partial P / \\partial \\rho ) \\right|_S \\)
:math:`c_s^2 = \left. ( \partial P / \partial \rho ) \right|_S`
as described in
`Kegerreis et al. (2019) <https://doi.org/10.1093/mnras/stz1606>`_.
Note that there is a typo in the sign of
\\(du = T dS - P dV = T dS + (P / \\rho^2) d\\rho \\) in the appendix;
the correct version was used in the actual derivation.
:math:`du = T dS - P dV = T dS + (P / \rho^2) d\rho` in the appendix,
but the correct version was used in the actual derivation.
The ideal gas uses the same equations detailed in :ref:`equation_of_state`.
......
......@@ -16,5 +16,18 @@ and which EoS it should use.
The Balsara viscosity switch is used by default, but can be disabled by
compiling SWIFT with ``make CFLAGS=-DPLANETARY_SPH_NO_BALSARA``.
Note: the boundary-improvement methods presented in Ruiz-Bonilla+2022 are
in the process of being improved for release with better documentation etc soon.
Note: to access the boundary-improvement method presented in Ruiz-Bonilla+2022,
use the ``planetary_imbalance_RB22`` git branch and compile with
``--with-hydro=planetary-gdf``. However, we instead recommend using the REMIX
SPH scheme, as it has effectively replaced this method.
.. _planetary_remix_hydro:
REMIX SPH
======================
REMIX is an SPH scheme designed to alleviate effects that typically suppress
mixing and instability growth at density discontinuities in SPH simulations
(Sandnes et al. 2025), and also includes the multiple EoS options as the base
Planetary scheme. For more information on what is included in the REMIX
scheme and how to configure SWIFT to use REMIX, see: :ref:`remix_sph`.
......@@ -7,24 +7,36 @@ Planetary Simulations
=====================
SWIFT is also designed for running planetary simulations
such as of giant impacts, as introduced in Kegerreis+2019,
and any other types of simulations with more complicated equations of state
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,
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/``, as well as some hydro tests for comparison with other
schemes. 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: 1
......
......@@ -249,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
......@@ -267,7 +273,9 @@ 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 [#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 |
......@@ -324,6 +332,12 @@ 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
---------------------------
......@@ -357,12 +371,16 @@ 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.
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.
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
-----------------------------------------
......
......@@ -10,22 +10,41 @@ In order to launch jets, we introduce a jet reservoir that functions identically
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:`\vec{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`.
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 then add the kick vector to the particle's current velocity. We do this in a way that conserves energy, so that the magnitude of the final velocity is computed from
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::
\frac{1}{2}m_\mathrm{gas}\vec{v}_\mathrm{fin}^2=\frac{1}{2}m_\mathrm{gas}\vec{v}_\mathrm{0}^2 + \frac{1}{2}m_\mathrm{gas}\vec{v}_\mathrm{kick}^2,
(\vec{v}_0+\Delta\vec{v})^2 = \vert \vec{v}_0\vert ^2 + v_\mathrm{j}^2,
while its direction is computed from conservation of momentum:
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::
m_\mathrm{gas}\vec{v}_\mathrm{fin}=m_\mathrm{gas}\vec{v}_\mathrm{0} + m_\mathrm{gas}\vec{v}_\mathrm{kick}.
\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
---------------------
Black holes will generally have time steps based on their gravitational interactions, but also based on their current accretion rate and the expected time interval until which the thermal feedback reservoir (representing radiative feedback) will have grown enough to heat 1 particle. We introduce a similar time step, but based on when the jet reservoir grows enough to kick :math:`N_\mathrm{j}` particles. We then take the minimum of those two.
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},
On top of that, we add a time step that makes sure the BH spin doesn't get evolved too much over one time step. Its magnitude is actually not a problem, since the growth of the spin magnitude is always tied with the growth of mass. However, the direction is more problematic (especially in the thin disk case, where alignment can occur very quickly, with little mass or spin growth, due to large warp radii). For this reason, we make sure that the amount of warp angular momentum interacting with the BH over the next time-step, :math:`\Delta_J=(J_\mathrm{warp}/M_\mathrm{warp})\Delta M`, is a small fraction of the current BH angular momentum :math:`J_\mathrm{BH}` (e.g. :math:`0.01`).
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.
......@@ -49,37 +49,37 @@ Below we give an example of parameter choices applicable for e.g. a 50 Mpc box.
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 thin and thick disk, as dot(m) = 0.2 * alpha^2.
mdot_crit_ADAF: 0.03 # 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.
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_cs_ratio: 10. # This sets the jet velocity to v_jet_cs_ratio times the sound speed of the hot gas of the parent halo the black hole is in. This is used if 'AGN_jet_velocity_model' is 'BlackHoleMass'.
v_jet_BH_mass_scaling_reference_mass_Msun: 3.4e3 # 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.65 # 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_mass_loading: 400. # The constant mass loading to use if 'AGN_jet_velocity_model' is 'MassLoading'.
v_jet_xi: 1. # The numerical multiplier by which the jet velocity formula is scaled, if 'AGN_jet_velocity_model' is 'Local'. If a value of 1 is used, the formulas are used as derived, i.e. they are not rescaled.
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: 1e5 # The maximal jet velocity. This is used if 'AGN_jet_velocity_model' is 'BlackHoleMass', 'MassLoading' or 'Local'.
temperature_hot_gas_min_K: 1e4 # The minimum temperature (in Kelvin) of hot gas to include in the calculation of the hot gas sound speed around the BH. This is only used if 'AGN_jet_velocity_model' is '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'
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]. If used, jets will be launched exclusively along the z axis. Should be set to 1 only for tests.
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: 0.1 # The accretion efficiency (suppression factor of the accretion rate) to use in the thick disc (ADAF), to represent the effects of subgrid winds that take away most of the mass flowing through the accretion disc.
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.
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.
include_ADIOS_suppression: 0 # Whether to suppress the accretion rate in the fully thick disc regime [1] (Eddington rate below 0.2alpha^2) by the amount expected to be taken away by isotropic kinetic disk winds.
turn_off_secondary_feedback: 1 # If set to 1, there will be only radiative (thermal) feedback in the thin disk mode, and only jets in the thick disk mode.
jet_h_r_slope: 1. # The slope of the dependence of jet efficiency on aspect ratio of the subgrid accretion disk, H/R. Default value is 1, and another reasonable value is 0 (same jet efficiency for all disks). Reality could be anything in between. This parameter is only used if turn_off_secondary_feedback is set to 0.
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: 0 # 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.
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), 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.
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.
If you set ``use_var_v_jet: 0``, you will need to change ``v_jet``, the kicking velocity of particles, depending on what system you are simulating. You should typically choose values at least 10 times larger than the sound speed of the hot gas in your most massive haloes (e.g. 1500 km/s for a MW-type galaxy and 10 000 km/s for a :math:`10^{14}` :math:`\mathrm{M}_\odot` halo). If, on the other hand, you set ``use_var_v_jet: 1``, the launching velocities will vary on their own depending on the typical sound speed (virial velocity) of the hot gas in the haloes. You then need to set ``v_jet_cs_ratio`` to values :math:`\gg1` (10-30 works well) in order to have significant shocking.
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).
......@@ -188,6 +188,39 @@ def L_adv(x, alpha):
) / 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)
......@@ -196,6 +229,7 @@ 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
......@@ -204,17 +238,17 @@ 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.028 for x in a], color="blue", alpha=0.2)
plt.fill_between(a, [0.028 for x in a], mdotcrit1, color="red", alpha=0.2)
plt.fill_between(a, mdotcrit1, [375 for x in a], color="orange", alpha=0.2)
plt.ylabel("$\dot{m}$", fontsize=24, usetex=True)
plt.xlabel("$a$", fontsize=24, usetex=True)
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.001, 100])
plt.text(-0.22, 0.004, "Thick disc", fontsize=20)
plt.text(-0.2, 0.33, "Thin disc", fontsize=20)
plt.text(-0.2, 18, "Slim disc", fontsize=20)
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",
......@@ -260,6 +294,7 @@ plt.tick_params(
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 = (
......@@ -271,14 +306,14 @@ jet_factor = (
* horizon_ang_vel ** 2
* (1.0 + 1.38 * horizon_ang_vel ** 2 - 9.2 * horizon_ang_vel ** 4)
)
Z1_j = np.array(
Z_1 = np.array(
[
1 + (1 - x ** 2) ** 0.333 * ((1 + abs(x)) ** 0.333 + (1 - abs(x)) ** 0.333)
for x in a
]
)
Z2_j = np.array(np.sqrt(3 * a ** 2 + Z1_j ** 2))
r_iso = 3 + Z2_j - np.sign(np.array(a)) * np.sqrt((3 - Z1_j) * (3 + Z1_j + 2 * Z2_j))
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)
......@@ -286,6 +321,7 @@ 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
......@@ -317,109 +353,82 @@ 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)
fig = plt.figure(figsize=(18, 4))
fig.subplots_adjust(wspace=0, hspace=0, top=1, bottom=0)
gs = gridspec.GridSpec(1, 3, width_ratios=[1, 1, 1])
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,
eps_ADAF2,
100 * 0.005 * (1 + 3 * (phi / 50) ** 2 * (horizon_ang_vel / 0.2) ** 2),
linewidth=2,
label="$\epsilon_\mathrm{rad}(\dot{m}<0.0044)$",
color="red",
label=r"$\epsilon_\mathrm{wind,thick}$",
color="blue",
)
plt.plot(
a,
eps_ADAF1,
100 * 0.1 * (1 - np.sqrt(1 - 2 / (3 * r_iso))),
linewidth=2,
label="$\epsilon_\mathrm{rad}(\dot{m}=0.028)$",
label=r"$\epsilon_\mathrm{f}\epsilon_\mathrm{rad,NT}$ $\mathrm{(thin}$ $\mathrm{disc})$",
color="red",
linestyle="--",
)
plt.plot(a, 0.97 * Jet_ADAF, linewidth=2, label="$\epsilon_\mathrm{jet}$", color="blue")
plt.fill_between(a, eps_ADAF1, eps_ADAF2, color="red", alpha=0.2)
plt.ylabel("$\epsilon_\mathrm{feedback}$", fontsize=24, usetex=True)
plt.xlabel("$a$", fontsize=24, usetex=True)
plt.tick_params(axis="y", right=True, direction="in")
plt.legend(loc="upper left", prop={"size": 13})
plt.xticks([-1, -0.5, 0, 0.5, 1], [-1, -0.5, 0, 0.5, 1])
plt.yticks(
[0.0001, 0.001, 0.01, 0.1, 1, 10, 100],
["", 10 ** (-3), 10 ** (-2), 10 ** (-1), 10 ** (-0), 10 ** (1), 10 ** (2)],
)
plt.minorticks_on()
plt.axis([-1, 1, 0.0001, 10])
plt.yscale("log")
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.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.title("Thick disc", fontsize=16)
plt.subplot(gs[1])
plt.plot(a, 0.97 * eps_TD, linewidth=2, label="$\epsilon_\mathrm{rad}$", color="red")
plt.plot(
a,
Jet_TD2,
linewidth=2,
label="$\epsilon_\mathrm{jet}(\dot{m}=0.028,M_\mathrm{BH}=10^9 \mathrm{M}_\odot)$",
color="blue",
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,
Jet_TD1,
linewidth=2,
label="$\epsilon_\mathrm{jet}(\dot{m}=1,M_\mathrm{BH}=10^6 \mathrm{M}_\odot)$",
color="blue",
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.fill_between(a, Jet_TD1, Jet_TD2, color="blue", alpha=0.2)
plt.xlabel("$a$", fontsize=24, usetex=True)
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")
plt.yscale("log")
plt.legend(loc="upper left", prop={"size": 13})
plt.xticks([-1, -0.5, 0, 0.5, 1], ["", -0.5, 0, 0.5, 1])
plt.axis([-1, 1, 0.0001, 10])
plt.yticks([0.001, 0.01, 0.1, 1, 10], ["", "", "", "", ""])
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",
......@@ -460,29 +469,58 @@ plt.tick_params(
which="minor",
labelsize=16,
)
plt.title("Thin disc", fontsize=16)
plt.title("Wind efficiency", fontsize=16)
plt.subplot(gs[2])
plt.subplot(gs[1])
plt.plot(
a, eps_SD1, linewidth=2, label="$\epsilon_\mathrm{rad}(\dot{m}=1)$", color="red"
a, 100 * Jet_ADAF, linewidth=2, label=r"$\epsilon_\mathrm{jet,thick}$", color="blue"
)
plt.plot(
a,
eps_SD2,
linewidth=2,
label="$\epsilon_\mathrm{rad}(\dot{m}=50)$",
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, Jet_SD, linewidth=2, label="$\epsilon_\mathrm{jet}$", color="blue")
plt.fill_between(a, eps_SD1, eps_SD2, color="red", alpha=0.2)
plt.xlabel("$a$", fontsize=24, usetex=True)
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")
plt.yscale("log")
plt.legend(loc="upper left", prop={"size": 13})
plt.xticks([-1, -0.5, 0, 0.5, 1], ["", -0.5, 0, 0.5, 1])
plt.axis([-1, 1, 0.0001, 10])
plt.yticks([0.001, 0.01, 0.1, 1, 10], ["", "", "", "", ""])
pylab.legend(loc="upper left", prop={"size": 15})
plt.axis([-1, 1, 0, 200])
plt.minorticks_on()
plt.tick_params(
axis="x",
......@@ -524,7 +562,7 @@ plt.tick_params(
which="minor",
labelsize=16,
)
plt.title("Slim disc", fontsize=16)
plt.title("Jet efficiency", fontsize=16)
plt.savefig("efficiencies.png", bbox_inches="tight")
......@@ -532,7 +570,7 @@ L_isco1 = [2 / 3 * 1 / np.sqrt(3) * (1 + 2 * np.sqrt(3 * r_isco(x) - 2)) for x i
plt.style.use("classic")
fig = plt.figure(figsize=(8, 6), linewidth=4)
plt.plot(a, L_isco1, linewidth=2, label="$\ell_\mathrm{ISCO}$", color="red")
plt.plot(a, L_isco1, linewidth=2, label=r"$\ell_\mathrm{ISCO}$", color="red")
plt.plot(
a,
0.45 * np.array(L_isco1),
......@@ -562,8 +600,8 @@ plt.plot(
label=r"$\ell_\mathrm{adv},\alpha=0.3$",
color="purple",
)
plt.ylabel("$\ell_\mathrm{in}$", fontsize=24, usetex=True)
plt.xlabel("$a$", fontsize=24, usetex=True)
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()
......@@ -670,125 +708,75 @@ da_SD_Benson = (
)
fig = plt.figure(figsize=(18, 4))
fig.subplots_adjust(wspace=0, hspace=0, top=1, bottom=0)
gs = gridspec.GridSpec(1, 3, width_ratios=[1, 1, 1])
fig = plt.figure(figsize=(7, 5))
plt.style.use("classic")
plt.subplot(gs[0])
plt.plot(a, da_TD_acc_only, linewidth=2, label="Accretion only", color="black")
plt.plot(a, da_TD_Benson, linewidth=1.5, label="Jet spindown included", color="blue")
plt.plot(a, [0 for x in a], linewidth=1.5, color="black", linestyle="--")
plt.ylabel("$\mathrm{d}a/\mathrm{d}\ln M_\mathrm{BH,0}$", fontsize=24, usetex=True)
plt.xlabel("$a$", fontsize=24, usetex=True)
plt.tick_params(axis="y", right=True, direction="in")
plt.legend(loc="lower left", prop={"size": 15})
plt.minorticks_on()
plt.axis([-1, 1, -4, 7])
plt.tick_params(
axis="x",
direction="in",
bottom=True,
top=True,
length=8,
width=1.2,
which="major",
labelsize=16,
z1 = np.array(
[
1 + (1 - x ** 2) ** 0.333 * ((1 + abs(x)) ** 0.333 + (1 - abs(x)) ** 0.333)
for x in a
]
)
plt.tick_params(
axis="y",
direction="in",
left=True,
right=True,
length=8,
width=1.2,
which="major",
labelsize=16,
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.tick_params(
axis="x",
direction="in",
bottom=True,
top=True,
length=4,
width=0.9,
which="minor",
labelsize=16,
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.tick_params(
axis="y",
direction="in",
left=True,
right=True,
length=4,
width=0.9,
which="minor",
labelsize=16,
plt.plot(
a,
s(0.1, a),
linewidth=2,
label=r"Thin disk, $f_\mathrm{Edd}=0.1$",
linestyle="-.",
color="red",
)
plt.title("Thin disc", fontsize=16)
plt.subplot(gs[1])
plt.plot(a, da_ADAF_acc_only, linewidth=2, label="Accretion only", color="black")
plt.plot(a, da_ADAF_Benson, linewidth=1.5, label="Jet spindown included", color="blue")
plt.plot(a, [0 for x in a], linewidth=1.5, color="black", linestyle="--")
plt.xlabel("$a$", fontsize=24, usetex=True)
plt.tick_params(axis="y", right=True, direction="in")
plt.xticks([-1.0, -0.5, 0.0, 0.5, 1.0], ["", -0.5, 0.0, 0.5, 1.0])
plt.yticks([-8, -6, -4, -2, 0, 2, 4, 6, 8], ["", "", "", "", "", "", "", "", ""])
plt.minorticks_on()
plt.axis([-1, 1, -4, 7])
plt.tick_params(
axis="x",
direction="in",
bottom=True,
top=True,
length=8,
width=1.2,
which="major",
labelsize=16,
plt.plot(
a,
s(1, a),
linewidth=2,
label=r"Thin disk, $f_\mathrm{Edd}=1$",
linestyle="--",
color="red",
)
plt.tick_params(
axis="y",
direction="in",
left=True,
right=True,
length=8,
width=1.2,
which="major",
labelsize=16,
plt.plot(
a,
s(10, a),
linewidth=2,
label=r"Slim disk, $f_\mathrm{Edd}=10$",
linestyle="--",
color="orange",
)
plt.tick_params(
axis="x",
direction="in",
bottom=True,
top=True,
length=4,
width=0.9,
which="minor",
labelsize=16,
plt.plot(
a,
s(100, a),
linewidth=2,
label=r"Slim disk, $f_\mathrm{Edd}=100$",
linestyle="-",
color="orange",
)
plt.tick_params(
axis="y",
direction="in",
left=True,
right=True,
length=4,
width=0.9,
which="minor",
labelsize=16,
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.title("Thick disc", fontsize=16)
plt.subplot(gs[2])
plt.plot(a, da_SD_acc_only, linewidth=2, label="Accretion only", color="black")
plt.plot(a, da_SD_Benson, linewidth=1.5, label="Jet spindown included", color="blue")
plt.plot(a, [0 for x in a], linewidth=1.5, color="black", linestyle="--")
plt.xlabel("$a$", fontsize=24, usetex=True)
plt.tick_params(axis="y", right=True, direction="in")
plt.xticks([-1.0, -0.5, 0.0, 0.5, 1.0], ["", -0.5, 0.0, 0.5, 1.0])
plt.yticks([-8, -6, -4, -2, 0, 2, 4, 6, 8], ["", "", "", "", "", "", "", "", ""])
pylab.legend(loc="lower left", prop={"size": 13})
plt.minorticks_on()
plt.axis([-1, 1, -4, 7])
plt.axis([-1, 1, -10, 10])
plt.tick_params(
axis="x",
direction="in",
......@@ -829,6 +817,5 @@ plt.tick_params(
which="minor",
labelsize=16,
)
plt.title("Slim disc", fontsize=16)
plt.savefig("spinup.png", bbox_inches="tight")
......@@ -11,7 +11,7 @@ Here we provide a comprehensive summary of the model. In order to more easily vi
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. 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.
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.
......@@ -24,42 +24,63 @@ Transitions from one accretion state to another
:figclass: align-center
:alt: Accretion regimes
The type of accretion disk surrounding the BHs depending on their accretion rates and spins. The transition between the thick and thin disk is calculated assuming the viscosity parameter :math:`\alpha=0.2`, while the transition from thin to slim disk is assumed to occur when the latter is :math:`F=0.5` times as radiatively efficienct as the former.
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:`\dot{m}`. We assume that the subgrid accretion disk is thick for :math:`\dot{m}<0.03`, based on observations (`Russell et al. 2013 <https://ui.adsabs.harvard.edu/abs/2013MNRAS.432..530R/abstract>`_). This also allows us to constrain the value of one of the main parameters in any accretion model: the viscosity parameter :math:`\alpha` (which is related to the kinematic viscosity :math:`\nu`and sound speed :math:`c_\mathrm{s}` through :math:`\nu=\alpha c_\mathrm{s}H`, with :math:`c_\mathrm{s}` the sound speed and :math:`H` the disk half-thickness). Numerical calculations suggest that thick disks are present for :math:`\dot{m}<0.4\alpha^2` (`Yuan & Narayan 2014 <https://ui.adsabs.harvard.edu/abs/2014ARA%26A..52..529Y/abstract>`_), and this agrees with observations if :math:`\alpha=0.25-0.3`. These values agree very well with more direct observational estimates, which suggest :math:`\alpha=0.2-0.3` (`Martin et al. 2019 <https://ui.adsabs.harvard.edu/abs/2014ARA%26A..52..529Y/abstract>`_).
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`.
The transition from the thin to the slim disk should occur around :math:`\dot{m}\approx 1`. However, the exact physics of this transition is not well understood. There is likely some spin dependence of the critical accretion rate, due to different radiative physics depending on spin. One of the main properties of slim disks is that they are less radiatively efficient than thin disks (`Sadowski et al. 2014 <https://ui.adsabs.harvard.edu/abs/2014MNRAS.439..503S/abstract>`_). We thus assume that the transition occurs when the radiative efficiency of a slim disk, :math:`\epsilon_\mathrm{r,SD}`, falls below some fraction of the radiative efficiency of a thin disk, :math:`\epsilon_\mathrm{r,TD}`. We quantify this as :math:`\epsilon_\mathrm{r,SD}<F\epsilon_\mathrm{r,SD}`, with :math:`F\approx 0.5` a free parameter. We give the expressions for both of the efficiencies below.
Accretion efficiencies
-----------------------------------------------
Jet 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
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::
\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{j}=\frac{\kappa}{4\pi}\bigg(\frac{H/R}{0.3}\bigg)^\eta \phi_\mathrm{BH}^2\Omega_\mathrm{BH}^2\big(1+1.38\Omega_\mathrm{BH}^2-9.2\Omega_\mathrm{BH}^4\big),
\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
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::
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::
\phi_\mathrm{BH}(a)=-20.2a^3-14.9a^2+34a+52.6.
R_\mathrm{tr} = R_1 \bigg(\frac{0.01}{f_\mathrm{Edd}}\bigg)^2,
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.
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:
The dependence of the jet efficiency on the type of accretion disk comes through the factor that depends on the aspect ratio :math:`H/R`, since accretion disks differ in this quantity. Theoretical, self-similar models of thick disks suggest :math:`H/R=0.5` (`Narayan & Yi 1995b <https://ui.adsabs.harvard.edu/abs/1995ApJ...452..710N/abstract>`_), but we instead take :math:`H/R=0.3`, more in line with simulations. For slim disks, which have received less attention in simulations, we assume the value :math:`H/R=1/(2\sqrt{5})\approx 0.2` (based on the self-similar model by `Wang & Zhou. 1999 <https://ui.adsabs.harvard.edu/abs/1999ApJ...516..420W/abstract>`_).
.. 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
----------------
Thin disks are, not surprisingly, much thinner. The value of :math:`H/R` in this regime is not a constant, but rather depends on the BH mass and accretion rate, slightly on radius and also on the viscosity parameter :math:`\alpha`. Thin disks have three different regions in the `Shakura & Sunyaev (1973) <https://ui.adsabs.harvard.edu/abs/1973A%26A....24..337S/abstract>`_ model. For simplicity, we model the whole disk as being represented with only one region. In region a), the innermost one, radiation dominates over gas pressure. It is typically very small or doesn't exist at all, so we disregard it as a possibility. In regions b) and c), gas pressure dominates over radiation pressure. In b), electrons dominate in the opacity, while in c), free-free absorption dominates. We leave both regions as a possibility, and leave the choice as a free parameter in the model (not likely to lead to large differences in galaxy/BH evolution). The expressions for the aspect ratio in these regions are
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::
\bigg(\frac{H}{R}\bigg)_\mathrm{TD,b} = 1.25\times10^{-3} \alpha^{-1/10}\dot{m}^{1/5}\bigg(\frac{M_\mathrm{BH}}{10^8\hspace{0.5mm}\mathrm{M}_\odot}\bigg)^{-1/10}\bigg(\frac{R}{2R_\mathrm{G}}\bigg)^{1/20}
\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),
in region b) and
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::
\bigg(\frac{H}{R}\bigg)_\mathrm{TD,c} = 1.15\times10^{-3} \alpha^{-1/10}\dot{m}^{3/20}\bigg(\frac{M_\mathrm{BH}}{10^8\hspace{0.5mm}\mathrm{M}_\odot}\bigg)^{-1/10}\bigg(\frac{R}{2R_\mathrm{G}}\bigg)^{1/8}
\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.
in region c).
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
The jet efficiency also depends on the slope :math:`\eta`. Classical jet theory (`Meier 2001 <https://ui.adsabs.harvard.edu/abs/2001Sci...291...84M/abstract>`_) suggests that jet powers depend on the aspect ratio linearly, so :math:`\eta=1`. This is also in line with some simulations finding a reduction in jet efficiencies with the aspect ratio (e.g. `Tchekhovskoy et al. 2014 <https://ui.adsabs.harvard.edu/abs/2014MNRAS.437.2744T/abstract>`_). In this scenario, jets launched from thin disks are of order :math:`\approx100` times less powerful than those launched from thick disks. On the other hand, some simulations of thin disks have found jet efficiencies similar to thick disk ones (e.g. `Liska et al. 2019 <https://ui.adsabs.harvard.edu/abs/2019MNRAS.487..550L/abstract>`_), which is supported by observations of blazars (`Ghisellini et al. 2014 <https://ui.adsabs.harvard.edu/abs/2014Natur.515..376G/abstract>`_). In this picture, thin jets are approximately as efficient as thick disk ones, which can in our case be implemented as :math:`\eta=0`. The reality is likely to be somwhere in between. Note that the choice of :math:`\eta` likely has a strong impact on the evolution of galaxies and BHs; our default choice is the classical picture in which :math:`\eta=1`.
.. 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
......@@ -69,8 +90,8 @@ The jet efficiency also depends on the slope :math:`\eta`. Classical jet theory
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 efficiencies
----------------------
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
......@@ -79,7 +100,7 @@ In the EAGLE and COLIBRE models, all subgrid accretion disks are effectively thi
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:`\dot{m}<\dot{m}_\mathrm{crit,visc}`, viscous heating dominates the heating of electrons, whereas for :math:`\dot{m}_\mathrm{crit,visc}<\dot{m}<\dot{m}_\mathrm{crit,ADAF}`, it is dominated by ion-electron heating. Here, :math:`\dot{m}_\mathrm{crit,visc}` is the transitional value between the two thick disc (ADAF) regimes, and :math:`\dot{m}_\mathrm{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
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),
......@@ -87,30 +108,32 @@ In the thick disk regime, radiative efficiencies are lower by a factor :math:`\a
while in the ion-heating regime it is given by
.. math::
\epsilon_\mathrm{r,ADAF}=0.2\epsilon_\mathrm{r,TD}\bigg(\frac{\dot{m}}{\alpha^2}\bigg)\bigg(\frac{\beta}{0.5}\bigg)\bigg(\frac{6}{r_\mathrm{ISCO}}\bigg).
\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::
\dot{m}_\mathrm{crit,visc}=0.0002\bigg(\frac{\delta_\mathrm{ADAF}}{0.0005}\bigg)\bigg(\frac{1-\beta}{\beta}\bigg)\alpha^2.
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}{\dot{m}}A(a)\bigg( \frac{0.985}{1.6/\dot{m}+B(a)}+\frac{0.015}{1.6/\dot{m}+C(a)}\bigg),
\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:`\dot{m}\gtrapprox1`, but depending on spin), the radiative efficiency drops. These two formulas are used to decide when a disk transitions from thin to slim.
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.
Evolution of the black hole spin magnitude
------------------------------------------
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>`_:
.. figure:: spec_ang_mom.png
:width: 600px
:align: center
:figclass: align-center
:alt: Angular momenta
.. 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:
Dimensionless pecific angular momentum of the thin disk at the innermost stable circular orbit (ISCO, solid red line), compared with the specific angular momentum at the inner radius (the event horizon) for advection-dominated flows (the thick and slim disk) for a few values of the viscosity parameter :math:`\alpha`. The dashed red line shows that the latter can be approximated as :math:`45` per cent of the former.
.. 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.
......@@ -119,16 +142,34 @@ In the absence of jet spindown, the evolution of angular momentum is given simpl
.. 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).
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.
For the thin disc, the inner radius :math:`R_\mathrm{in}` can be taken to be the radius of the ISCO, since orbits should quickly degrade within it. We thus take :math:`\ell_\mathrm{in}` as the specific angular momentum at the ISCO for the thin disc (the expression for which is given in e.g. Appendix B of `Fiacconi et al. 2018 <https://ui.adsabs.harvard.edu/abs/2018MNRAS.477.3807F/abstract>`_). For the thin disk, the spinup function (the equation shown above) is always positive, meaning that the BH will always be spun up to positive values. This means that the BH will be spun down if spin is negative, or spun up to an equilibrium value of :math:`a_\mathrm{eq}=1` if spin is positive. For advection-dominated disks (the thick and slim disk), we assume that :math:`\ell_\mathrm{in}` is :math:`45` per cent of the ISCO value, based on numerical GR calculations by `Popham & Gammie (1998) <https://ui.adsabs.harvard.edu/abs/1998ApJ...504..419P/abstract>`_. We base this assumption on fits of the `Popham & Gammie (1998) <https://ui.adsabs.harvard.edu/abs/1998ApJ...504..419P/abstract>`_ results done by `Benson & Babul (2009) <https://ui.adsabs.harvard.edu/abs/2009MNRAS.397.1302B/abstract>`_. We independently compared these fits to the ISCO values and found :math:`\ell_\mathrm{in}\approx0.45\ell_\mathrm{ISCO}` with no more than :math:`10` per cent error for all values of spin and relevant values of :math:`\alpha=0.1-0.3`.
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>`_:
For the thick and slim disk, these lower specific angular momenta lead to a non-zero equilibrium spin value :math:`a_\mathrm{eq}<1`. If :math:`a>a_\mathrm{eq}`, the BH will be spun down due to frame-dragging and viscosity; the frame-dragging rotationally accelerates any accreting gas (on account of the BH angular momentum), while viscosity carries away some of that angular momentum. Including jets into the model leads to further spindown. The jet spindown term (to be added to the spinup equation above) can be derived as
.. 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}\ln M_\mathrm{BH,0}}\bigg)_\mathrm{j}=-\epsilon_\mathrm{j}(a)\frac{\sqrt{1-a^2}}{a}\bigg[\Big(\sqrt{1-a^2}+1 \Big)^2+a^2 \bigg]
\bigg(\frac{\mathrm{d}a}{\mathrm{d}M_\mathrm{BH,0}/M_\mathrm{BH}}\bigg)_\mathrm{thin/slim}=s_\mathrm{HD} - s_\mathrm{EM},
(see `Benson & Babul 2009 <https://ui.adsabs.harvard.edu/abs/2009MNRAS.397.1302B/abstract>`_ for a derivation, which we have independently verified). Including jet spindown leads to even lower equilibrium spin values; e.g. for the thick disk this is only :math:`a_\mathrm{eq}\approx0.25`.
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
......@@ -136,7 +177,7 @@ For the thick and slim disk, these lower specific angular momenta lead to a non-
:figclass: align-center
:alt: Spinup/spindown function
Spinup/spindown function (the rate of black hole spin evolution) as a function of spin for all three accretion disk types. Black lines show evolution with only accretion included, while blue lines show the total including jet spindown. These plots show that the thin disk is always spun up to to :math:`a_\mathrm{eq}=1`, even with jets (due to low jet efficiencies). The advection-dominated disks (thick and slim disk) are spun up to positive equilibrium values :math:`a_\mathrm{eq}<1`, or spun down to such an equilibrium value if :math:`a>a_\mathrm{eq}`. This is due to extraction of rotational energy from the BH by frame dragging and transport of the angular momentum to large distances through viscous forces. Including jet spindown pushes these equilibrium spins to even smaller values.
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
------------------------------------------
......@@ -149,7 +190,7 @@ In all cases, Lense-Thirring torques are effective only within some radius :math
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.
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
----------------------------------------------------
......@@ -172,34 +213,34 @@ As mentioned already, Lense-Thirring torques have different effects depending on
(`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} \dot{m}^{3 / 5}\left(\frac{M_{\mathrm{BH}}}{10^{8} M_{\odot}}\right)^{1 / 8}\left(\frac{R}{R_{\mathrm{S}}}\right)^{-3 / 5},
\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} \dot{m}^{7/10}\left(\frac{M_{\mathrm{BH}}}{10^{8} M_{\odot}}\right)^{1 / 20}\left(\frac{R}{R_{\mathrm{S}}}\right)^{-3 / 4}.
\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} \dot{m}^{-1 / 4}\left(\frac{M_{\mathrm{BH}}}{10^{8} M_{\odot}}\right)^{1 / 8}
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}\dot{m}^{-6/35}\bigg(\frac{M_\mathrm{BH}}{10^8\hspace{0.5mm}\mathrm{M}_\odot} \bigg)^{4/35},
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} \dot{m}^{-18/51}\left(\frac{M_{\mathrm{BH}}}{10^{8} M_{\odot}}\right)^{-49/51}
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} \dot{m}^{-22/45}\left(\frac{M_{\mathrm{BH}}}{10^{8} M_{\odot}}\right)^{-52/45}
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.
......@@ -212,7 +253,10 @@ The exact behaviour of the thick and slim disk (which we will collectively call
.. 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 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
......@@ -234,9 +278,9 @@ where :math:`v_\mathrm{r}=-\alpha v_0 v_\mathrm{K}` is the radial velocity. Here
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. (2007) <https://journals.aps.org/prd/abstract/10.1103/PhysRevD.78.044002>`_, 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. The only flaw in these formulas is that they do not include the effects of gravitational radiation. However, the effects of this radiation is confined to a :math:`\approx10\%` level, and only if either of the spin vectors is aligned or counteraligned with the direction of the orbital angular momentum (if it is not, the fits are even more accurate).
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. (2007) <https://journals.aps.org/prd/abstract/10.1103/PhysRevD.78.044002>`_ can be calculated as
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),
......@@ -248,4 +292,14 @@ where :math:`q=M_2/M_1` is the mass ratio (such that :math:`M_2<M_1`), :math:`\m
\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.129`, :math:`s_5 = -0.384`, :math:`t_0 = -2.686`, :math:`t_2 = -3.454`, :math:`t_3 = 2.353`. 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}}`.
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}.
......@@ -36,7 +36,7 @@ Star formation
The AGORA model uses the :ref:`GEAR model <gear_star_formation>` scheme, however with the
``GEARStarFormation:star_formation_mode`` parameter set to ``agora``. Instead of requiring the gas
density to reach the pressure floor, we simply require it to be denser than a density
threshold defined by ``GEARStarFormation:density_threshold``.
threshold defined by ``GEARStarFormation:density_threshold_Hpcm3``.
Recommended parameters for the AGORA model should be:
......@@ -47,10 +47,10 @@ Recommended parameters for the AGORA model should be:
GEARStarFormation:
star_formation_mode: agora
star_formation_efficiency: 0.01
maximal_temperature: 1e10
maximal_temperature_K: 1e10
n_stars_per_particle: 1
min_mass_frac: 0.5
density_threshold: 1.67e-23
density_threshold_Hpcm3: 10
......@@ -127,9 +127,6 @@ Recommended parameters for the AGORA model should be:
.. _agora_pressure_floor:
Pressure Floor
......@@ -138,4 +135,12 @@ Pressure Floor
The AGORA model uses precisely the same pressure floor than the :ref:`GEAR model <gear_pressure_floor>`.
.. _agora_initial_conditions:
Initial Conditions
~~~~~~~~~~~~~~~~~~
Note that if in the initial conditions, the time of formation of a stellar particle is given (``BirthTime``)
and set to a negative value, the stellar particle will provide no feedback.
A similar behavior will be obtained if the parameter ``overwrite_birth_time`` is set to 1 and
``birth_time`` to -1.
......@@ -5,6 +5,25 @@
Basic model (others)
====================
Sinks: Simple Bondi-Hoyle accretion
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The ``Basic`` sink model provides a foundation on which new sink implementations could be built. It includes a prescription for Bondi-Hoyle gas accretion, and a method for sink-sink mergers that is a slightly simplified version of the implementation used in GEAR.
No other physics is implemented for this model. Sinks cannot form - to use this model, sink particles must already be present in the initial conditions. They also cannot spawn stars from the gas they accrete.
Bondi-Hoyle accretion can be done in one of two ways:
* Gas particles within the sink's kernel are stochastically swallowed entirely, with a probability set by the Bondi-Hoyle rate. Specifically, the probability is set by the current difference between the sink's subgrid mass (determined by the accretion rate) and its dynamical mass (which tracks the number of particles/sinks actually swallowed). This mode is equivalent to the EAGLE black hole accretion model.
* Gas particles within the sink's kernel are "nibbled" down to some minimal mass, which can be specified by the user. This method is equivalent to the black hole accretion model of Bahe et al. 2022.
This model has only two parameters that must be specified in your parameter ``yml`` file:
* ``BasicSink:use_nibbling``: determines whether accretion is done by "nibbling" or by swallowing outright.
* ``BasicSink:min_gas_mass_for_nibbling_Msun``: if using "nibbling", the minimum mass to which gas particles can be nibbled. A good default is half the original particle mass.
For an even more bare-bones starting point, the ``Default`` sink model contains no physics at all, and is a totally blank canvas on which to build your sink model.
Cooling: Analytic models
~~~~~~~~~~~~~~~~~~~~~~~~
......
......@@ -537,6 +537,15 @@ Note that the star formation rates are expressed in internal units and not in
solar masses per year as is the case in many other codes. This choice ensures
consistency between all the fields written to the snapshots.
Finally, the star formation model can also create more than one star particle
per gas particle in a star formation event. The number of particles generated is
controlled by a runtime parameter. All the stars generated share the same
properties. They are slightly displaced from each others using a random vector
of magnitude :math:`0.1h` added to the position of the gas particle. If only
one star is formed (as is the default), no displacement is added. If N stars are
formed per gas particle, each star is born with a mass of 1/N of the gas
particle's mass.
For a normal EAGLE run, that section of the parameter file reads:
.. code:: YAML
......@@ -559,7 +568,8 @@ For a normal EAGLE run, that section of the parameter file reads:
threshold_max_density_H_p_cm3: 10.0 # Maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
min_over_density: 57.7 # Over-density above which star-formation is allowed.
EOS_entropy_margin_dex: 0.5 # (Optional) Logarithm base 10 of the maximal entropy above the EOS at which stars can form.
num_of_stars_per_gas_particle: 1 # (Optional) The number star particles to form per gas particle converted to stars. (Defaults to 1. Must be > 0)
Alternatively, the code can also use a simple Schmidt law for the SF rate
:math:`\dot{m}_* = \epsilon_{ff} \times m_g \times \frac{3 \pi} {32 G
......
.. GEAR sub-grid model chemistry
Darwin Roduit, 30th March 2025
.. gear_chemistry:
.. _gear_chemistry:
Chemistry
=========
Feedback mechanisms such as supernova feedback transfer metals to the nearby gas particles. There is no mass exchange (mass advection) in SPH methods, as in grid-based or moving-mesh hydrodynamics solvers. As a result, the metals received by the gas particles are locked in these particles. Grid or moving-mesh codes have mass advection and often implement passive scalar advection of metals to model metal mixing. GEAR implements different methods to model metal mixing.
.. _gear_smoothed_metallicity:
Smoothed metallicity
--------------------
The smoothed metallicity scheme consists in using the SPH to smooth the metallicity of each particle over the neighbors. It is worth to point the fact that we are *not exchanging* any metals but only smoothing it. The parameter ``GEARChemistry:initial_metallicity`` set the (non smoothed) initial mass fraction of each element for all the particles and ``GEARChemistry:scale_initial_metallicity`` use the feedback table to scale the initial metallicity of each element according the Sun's composition. If ``GEARChemistry:initial_metallicity`` is negative, then the metallicities are read from the initial conditions.
For this chemistry scheme the parameters are:
.. code:: YAML
GEARChemistry:
initial_metallicity: 1 # Initial metallicity of the gas (mass fraction)
scale_initial_metallicity: 1 # Should we scale the initial metallicity with the solar one?
.. GEAR sub-grid model feedback and stellar evolution
Loic Hausammann, 17th April 2020
Darwin Roduit, 30th March 2025
.. _gear_stellar_evolution_and_feedback:
Stellar evolution and feedback
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The feedback is composed of a few different models:
- The initial mass function (IMF) defines the quantity of each type of stars,
- The lifetime of a star defines when a star will explode (or simply die),
- The supernovae of type II (SNII) defines the rates and yields,
- The supernovae of type Ia (SNIa) defines the rates and yields,
- The energy injection that defines how to inject the energy / metals into the particles.
Most of the parameters are defined inside a table (``GEARFeedback:yields_table``) but can be override with some parameters in the YAML file.
I will not describe theses parameters more than providing them at the end of this section.
Two different models exist for the supernovae (``GEARFeedback:discrete_yields``).
In the continuous mode, we integrate the quantities over the IMF and then explodes a floating point number of stars (can be below 1 in some cases).
In the discrete mode, we avoid the problem of floating points by rounding the number of supernovae (using a floor and randomly adding a supernovae depending on the fractional part) and then compute the properties for a single star at a time.
Initial mass function
^^^^^^^^^^^^^^^^^^^^^
GEAR is using the IMF model from `Kroupa (2001) <https://ui.adsabs.harvard.edu/abs/2001MNRAS.322..231K/abstract>`_.
We have a difference of 1 in the exponent due to the usage of IMF in mass and not in number.
We also restrict the mass of the stars to be inside :math:`[0.05, 50] M_\odot`.
Here is the default model used, but it can be easily adapted through the initial mass function parameters:
.. math::
\xi(m) \propto m^{-\alpha_i}\, \textrm{where}\,
\begin{cases}
\alpha_0 = 0.3,\, & 0.01 \leq m / M_\odot < 0.08, \\
\alpha_1 = 1.3,\, & 0.08 \leq m / M_\odot < 0.50, \\
\alpha_2 = 2.3,\, & 0.50 \leq m / M_\odot < 1.00, \\
\alpha_3 = 2.3,\, & 1.00 \leq m / M_\odot,
\end{cases}
Lifetime
^^^^^^^^
The lifetime of a star in GEAR depends only on two parameters: first its mass and then its metallicity.
.. math::
\log(\tau(m)) = a(Z) \log^2(m) + b(Z) \log(m) + c(Z) \\ \\
a(Z) = -40.110 Z^2 + 5.509 Z + 0.7824 \\
b(Z) = 141.929 Z^2 - 15.889 Z - 3.2557 \\
c(Z) = -261.365 Z^2 + 17.073 Z + 9.8661
where :math:`\tau` is the lifetime in years, :math:`m` is the mass of the star (in solar mass) and Z the metallicity of the star.
The parameters previously given are the default ones, they can be modified in the parameters file.
Supernovae II
^^^^^^^^^^^^^
The supernovae rate is simply given by the number of stars massive enough that end their life at the required time.
.. math::
\dot{N}_\textrm{SNII}(t) = \int_{M_l}^{M_u} \delta(t - \tau(m)) \frac{\phi(m)}{m} \mathrm{d}m
where :math:`M_l` and :math:`M_u` are the lower and upper mass limits for a star exploding in SNII, :math:`\delta` is the Dirac function and :math:`\phi` is the initial mass function (in mass).
The yields for SNII cannot be written in an analytical form, they depend on a few different tables that are based on the work of `Kobayashi et al. (2000) <https://ui.adsabs.harvard.edu/abs/2000ApJ...539...26K/abstract>`_ and `Tsujimoto et al. (1995) <https://ui.adsabs.harvard.edu/abs/1995MNRAS.277..945T/abstract>`_.
Supernovae Ia
^^^^^^^^^^^^^
The supernovae Ia are a bit more complicated as they involve two different stars.
.. math::
\dot{N}_\textrm{SNIa}(t) = \left( \int_{M_{p,l}}^{M_{p,u}} \frac{\phi(m)}{m} \mathrm{d}m \right) \sum_i b_i \int_{M_{d,l,i}}^{M_{d,u,i}}
\delta(t-\tau(m)) \frac{\phi_d(m)}{m}\mathrm{d}m
.. math::
\phi_d(m) \propto m^{-0.35}
where :math:`M_{p,l}` and :math:`M_{p,u}` are the mass limits for a progenitor of a white dwarf, :math:`b_i` is the probability to have a companion and
:math:`M_{d,l,i}` and :math:`M_{d,u,i}` are the mass limits for each type of companion.
The first parenthesis represents the number of white dwarfs and the second one the probability to form a binary.
+------------------+--------------------+-------------------+------------------+
| Companion | :math:`M_{d,l,i}` | :math:`M_{d,u,i}` | :math:`b_i` |
+==================+====================+===================+==================+
| Red giant | 0.9 | 1.5 | 0.02 |
+------------------+--------------------+-------------------+------------------+
| Main sequence | 1.8 | 2.5 | 0.05 |
+------------------+--------------------+-------------------+------------------+
The yields are based on the same papers than the SNII.
Supernova energy injection
^^^^^^^^^^^^^^^^^^^^^^^^^^
When a star goes into a supernova (type II and Ia), the stellar evolution determines how much energy (``GEARFeedback:supernovae_energy_erg`` TO BE CHECKED), mass and metals are released during the explosion. The energy can be ditributed as internal/thermal energy or as momentum. Thus, we need to distribute internal energy, momentum, mass and metals to the gas particles. We will group all these in the “fluxes” term.
We have two models for the distribution of these fluxes and the subgrid modelling of the supernovae : GEAR model and GEAR mechanical model. We describe the two schemes in the :ref:`gear_sn_feedback_models` page.
Generating a new table
^^^^^^^^^^^^^^^^^^^^^^
The feedback table is an HDF5 file with the following structure:
.. graphviz:: feedback_table.dot
where the solid (dashed) squares represent a group (a dataset) with the name of the object underlined and the attributes written below. Everything is in solar mass or without units (e.g. mass fraction or unitless constant).
In ``Data``, the attribute ``elts`` is an array of string with the element names (the last should be ``Metals``, it corresponds to the sum of all the elements), ``MeanWDMass`` is the mass of the white dwarfs
and ``SolarMassAbundances`` is an array of float containing the mass fraction of the different element in the sun.
In ``IMF``, ``n + 1`` is the number of part in the IMF, ``as`` are the exponent (``n+1`` elements), ``ms`` are the mass limits between each part (``n`` elements) and
``Mmin`` (``Mmax``) is the minimal (maximal) mass of a star.
In ``LifeTimes``, the coefficient are given in the form of a single table (``coeff_z`` with a 3x3 shape).
In ``SNIa``, ``a`` is the exponent of the distribution of binaries, ``bb1`` and ``bb2`` are the coefficient :math:`b_i` and the other attributes follow the same names than in the SNIa formulas.
The ``Metals`` group from the ``SNIa`` contains the name of each elements (``elts``) and the metal mass fraction ejected by each supernovae (``data``) in the same order. They must contain the same elements than in ``Data``.
Finally for the ``SNII``, the mass limits are given by ``Mmin`` and ``Mmax``. For the yields, the datasets required are ``Ej`` (mass fraction ejected [processed]), ``Ejnp`` (mass fraction ejected [non processed]) and one dataset for each element present in ``elts``. The datasets should all have the same size, be uniformly sampled in log and contains the attributes ``min`` (mass in log for the first element) and ``step`` (difference of mass in log between two elements).
.. code:: YAML
GEARFeedback:
supernovae_energy_erg: 0.1e51 # Energy released by a single supernovae.
yields_table: chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5 # Table containing the yields.
discrete_yields: 0 # Should we use discrete yields or the IMF integrated one?
GEARInitialMassFunction:
number_function_part: 4 # Number of different part in the IMF
exponents: [0.7, -0.8, -1.7, -1.3] # Exponents of each part of the IMF
mass_limits_msun: [0.05, 0.08, 0.5, 1, 50] # Limits in mass between each part of the IMF
GEARLifetime:
quadratic: [-40.1107, 5.50992, 0.782432] # Quadratic terms in the fit
linear: [141.93, -15.8895, -3.25578] # Linear terms in the fit
constant: [-261.366, 17.0735, 9.86606] # Constant terms in the fit
GEARSupernovaeIa:
exponent: -0.35 # Exponent for the distribution of companions
min_mass_white_dwarf_progenitor: 3 # Minimal mass of a progenitor of white dwarf
max_mass_white_dwarf_progenitor: 8 # Maximal mass of a progenitor of white dwarf
max_mass_red_giant: 1.5 # Maximal mass for a red giant
min_mass_red_giant: 0.9 # Minimal mass for a red giant
coef_red_giant: 0.02 # Coefficient for the distribution of red giants companions
max_mass_main_sequence: 2.6 # Maximal mass for a main sequence star
min_mass_main_sequence: 1.8 # Minimal mass for a main sequence star
coef_main_sequence: 0.05 # Coefficient for the distribution of main sequence companions
white_dwarf_mass: 1.38 # Mass of a white dwarf
GEARSupernovaeII:
interpolation_size: 200 # Number of elements for the interpolation of the data
.. GEAR sub-grid model
Loic Hausammann, 17th April 2020
Darwin Roduit, 30th March 2025
.. _gear_pressure_floor:
Pressure Floor
~~~~~~~~~~~~~~
In order to avoid the artificial collapse of unresolved clumps, a minimum in pressure is applied to the particles.
This additional pressure can be seen as the pressure due to unresolved hydrodynamics turbulence and is given by:
.. math::
P_\textrm{Jeans} = \frac{\rho}{\gamma} \frac{4}{\pi} G h^2 \rho N_\textrm{Jeans}^{2/3}
where :math:`\rho` is the density, :math:`\gamma` the adiabatic index, :math:`G` is the gravitational constant,
:math:`h` the kernel support and :math:`N_\textrm{Jeans}` (``GEARPressureFloor:jeans_factor`` in the parameter file) is the number of particle required in order to resolve a clump.
This must be directly implemented into the hydro schemes, therefore only a subset of schemes (Gadget-2, SPHENIX and Pressure-Energy) have the floor available.
In order to implement it, you need equation 12 in `Hopkins 2013 <https://arxiv.org/abs/1206.5006>`_:
.. math::
m_i \frac{\mathrm{d}v_i}{\mathrm{d}t} = - \sum_j x_i x_j \left[ \frac{P_i}{y_i^2} f_{ij} \nabla_i W_{ij}(h_i) + \frac{P_j}{y_j^2} f_{ji} \nabla_j W_{ji}(h_j) \right]
and simply replace the :math:`P_i, P_j` by the pressure with the floor (when the pressure is below the floor).
Here the :math:`x, y` are simple weights that should never have the pressure floor included even if they are related to the pressure (e.g. pressure-entropy).
.. code:: YAML
GEARPressureFloor:
jeans_factor: 10. # Number of particles required to suppose a resolved clump and avoid the pressure floor.
.. _gear_star_formation:
Star formation
~~~~~~~~~~~~~~
The star formation is done in two steps: first we check if a particle is in the star forming regime and then we use a stochastic approach to transform the gas particles into stars.
A particle is in the star forming regime if:
- The velocity divergence is negative (:math:`\nabla\cdot v < 0`),
- The temperature is lower than a threshold (:math:`T < T_t` where :math:`T_t` is defined with ``GEARStarFormation:maximal_temperature_K``),
- The gas density is higher than a threshold (:math:`\rho > \rho_t` where :math:`\rho_t` is defined with ``GEARStarFormation:density_threshold_Hpcm3``)
- The particle reaches the pressure floor (:math:`\rho > \frac{\pi}{4 G N_\textrm{Jeans}^{2/3} h^2}\frac{\gamma k_B T}{\mu m_p}` where :math:`N_\textrm{Jeans}` is defined in the pressure floor).
If ``GEARStarFormation:star_formation_mode`` is set to ``agora``, the condition on the pressure floor is ignored. Its default value is ``default``.
A star will be able to form if a randomly drawn number is below :math:`\frac{m_g}{m_\star}\left(1 - \exp\left(-c_\star \Delta t / t_\textrm{ff}\right)\right)` where :math:`t_\textrm{ff}` is the free fall time, :math:`\Delta t` is the time step of the particle and :math:`c_\star` is the star formation coefficient (``GEARStarFormation:star_formation_efficiency``), :math:`m_g` the mass of the gas particle and :math:`m_\star` the mass of the possible future star. The mass of the star is computed from the average gas mass in the initial conditions divided by the number of possible stars formed per gas particle (``GEARStarFormation:n_stars_per_particle``). When we cannot have enough mass to form a second star (defined with the fraction of mass ``GEARStarFormation:min_mass_frac``), we fully convert the gas particle into a stellar particle. Once the star is formed, we move it a bit in a random direction and fraction of the smoothing length in order to avoid any division by 0.
Currently, only the following hydro schemes are compatible: SPHENIX, Gadget2, minimal SPH, Gasoline-2 and Pressure-Energy.
Implementing the other hydro schemes is not complicated but requires some careful thinking about the cosmological terms in the definition of the velocity divergence (comoving vs non comoving coordinates and if the Hubble flow is included or not).
.. code:: YAML
GEARStarFormation:
star_formation_efficiency: 0.01 # star formation efficiency (c_*)
maximal_temperature_K: 3e4 # Upper limit to the temperature of a star forming particle
density_threshold_Hpcm3: 10 # Density threshold (Hydrogen atoms/cm^3) for star formation
n_stars_per_particle: 4 # Number of stars that an hydro particle can generate
min_mass_frac: 0.5 # Minimal mass for a stellar particle as a fraction of the average mass for the stellar particles.
Initial Conditions
++++++++++++++++++
Note that if in the initial conditions, the time of formation of a stellar particle is given (``BirthTime``)
and set to a negative value, the stellar particle will provide no feedback.
A similar behavior will be obtained if the parameter ``Stars:overwrite_birth_time`` is set to 1 and
``Stars:birth_time`` to -1.
.. _gear_grackle_cooling:
Cooling: Grackle
~~~~~~~~~~~~~~~~
Grackle is a chemistry and cooling library presented in `B. Smith et al. 2017 <https://ui.adsabs.harvard.edu/abs/2017MNRAS.466.2217S>`_
(do not forget to cite if used). Four different modes are available:
equilibrium, 6 species network (H, H\\( ^+ \\), e\\( ^- \\), He, He\\( ^+ \\)
and He\\( ^{++} \\)), 9 species network (adds H\\(^-\\), H\\(_2\\) and
H\\(_2^+\\)) and 12 species (adds D, D\\(^+\\) and HD). Following the same
order, the swift cooling options are ``grackle_0``, ``grackle_1``, ``grackle_2``
and ``grackle_3`` (the numbers correspond to the value of
``primordial_chemistry`` in Grackle). It also includes some metal cooling (on/off with ``GrackleCooling:with_metal_cooling``), self-shielding
methods and UV background (on/off with ``GrackleCooling:with_UV_background``). In order to use the Grackle cooling, you will need
to provide a HDF5 table computed by Cloudy (``GrackleCooling:cloudy_table``).
Configuring and compiling SWIFT with Grackle
++++++++++++++++++++++++++++++++++++++++++++
In order to compile SWIFT with Grackle, you need to provide the options ``with-chemistry=GEAR`` and ``with-grackle=$GRACKLE_ROOT``
where ``$GRACKLE_ROOT`` is the root of the install directory (not the ``lib``).
.. warning::
(State 2023) Grackle is experiencing current development, and the API is subject
to changes in the future. For convenience, a frozen version is hosted as a fork
on github here: https://github.com/mladenivkovic/grackle-swift .
The version available there will be tried and tested and ensured to work with
SWIFT.
Additionally, that repository hosts files necessary to install that specific
version of grackle with spack.
To compile it, run
the following commands from the root directory of Grackle:
``./configure; cd src/clib``.
Update the variables ``LOCAL_HDF5_INSTALL`` and ``MACH_INSTALL_PREFIX`` in
the file ``src/clib/Make.mach.linux-gnu``.
Finish with ``make machine-linux-gnu; make && make install``.
Note that we require the 64 bit float version of Grackle, which should be the default setting.
(The precision can be set while compiling grackle with ``make precision-64``).
If you encounter any problem, you can look at the `Grackle documentation <https://grackle.readthedocs.io/en/latest/>`_
You can now provide the path given for ``MACH_INSTALL_PREFIX`` to ``with-grackle``.
Parameters
++++++++++
When starting a simulation without providing the different element fractions in the non equilibrium mode, the code supposes an equilibrium and computes them automatically.
The code uses an iterative method in order to find the correct initial composition and this method can be tuned with two parameters. ``GrackleCooling:max_steps`` defines the maximal number of steps to reach the convergence and ``GrackleCooling:convergence_limit`` defines the tolerance in the relative error.
In the parameters file, a few different parameters are available.
- ``GrackleCooling:redshift`` defines the redshift to use for the UV background (for cosmological simulation, it must be set to -1 in order to use the simulation's redshift).
- ``GrackleCooling:provide_*_heating_rates`` can enable the computation of user provided heating rates (such as with the radiative transfer) in either volumetric or specific units.
- Feedback can be made more efficient by turning off the cooling during a few Myr (``GrackleCooling:thermal_time_myr``) for the particles touched by a supernovae.
- The self shielding method is defined by ``GrackleCooling:self_shielding_method`` where 0 means no self shielding, > 0 means a method defined in Grackle (see Grackle documentation for more information) and -1 means GEAR's self shielding that simply turn off the UV background when reaching a given density (``GrackleCooling:self_shielding_threshold_atom_per_cm3``).
- The initial elemental abundances can be specified with ``initial_nX_to_nY_ratio``, e.g. if you want to specify an initial HII to H abundance. A negative value ignores the parameter. The complete list can be found below.
A maximal (physical) density must be set with the ``GrackleCooling:maximal_density_Hpcm3 parameter``. The density passed to Grackle is *the minimum of this density and the gas particle (physical) density*. A negative value (:math:`< 0`) deactivates the maximal density, i.e. there is no maximal density limit.
The purpose of this parameter is the following. The Cloudy tables provided by Grackle are limited in density (typically to :math:`10^4 \; \mathrm{hydrogen \; atoms/cm}^3`). In high-resolution simulations, particles can have densities higher than :math:`10^4 \; \mathrm{hydrogen \; atoms/cm}^3`. This maximal density ensures that we pass a density within the interpolation ranges of the table, should the density exceed it.
It can be a solution to some of the following errors (with a translation of what the values mean):
.. code:: text
inside if statement solve rate cool: 0 0
MULTI_COOL iter > 10000 at j,k = 1 1
FATAL error (2) in MULTI_COOL
dt = 1.092E-08 ttmin = 7.493E-12
2.8E-19 // sub-cycling timestep
7.5E-12 // time elapsed (in the sub-cycle)
2.2E+25 // derivative of the internal energy
T
.. note::
This problem is particularly relevant with metal cooling enabled. Another solution is to modify the tables. But one is not exempted from exceeding the table maximal density value, since Grackle does not check if the particle density is larger than the table maximal density.
Complete parameter list
-----------------------
Here is the complete section in the parameter file:
.. code:: YAML
GrackleCooling:
cloudy_table: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
with_UV_background: 1 # Enable or not the UV background
redshift: 0 # Redshift to use (-1 means time based redshift)
with_metal_cooling: 1 # Enable or not the metal cooling
provide_volumetric_heating_rates: 0 # (optional) User provide volumetric heating rates
provide_specific_heating_rates: 0 # (optional) User provide specific heating rates
max_steps: 10000 # (optional) Max number of step when computing the initial composition
convergence_limit: 1e-2 # (optional) Convergence threshold (relative) for initial composition
thermal_time_myr: 5 # (optional) Time (in Myr) for adiabatic cooling after a feedback event.
self_shielding_method: -1 # (optional) Grackle (1->3 for Grackle's ones, 0 for none and -1 for GEAR)
self_shielding_threshold_atom_per_cm3: 0.007 # Required only with GEAR's self shielding. Density threshold of the self shielding
maximal_density_Hpcm3: 1e4 # Maximal density (in hydrogen atoms/cm^3) for cooling. Higher densities are floored to this value to ensure grackle works properly when interpolating beyond the cloudy_table maximal density. A value < 0 deactivates this parameter.
HydrogenFractionByMass : 1. # Hydrogen fraction by mass (default is 0.76)
use_radiative_transfer : 1 # Arrays of ionization and heating rates are provided
RT_heating_rate_cgs : 0 # heating rate in units of / nHI_cgs
RT_HI_ionization_rate_cgs : 0 # HI ionization rate in cgs [1/s]
RT_HeI_ionization_rate_cgs : 0 # HeI ionization rate in cgs [1/s]
RT_HeII_ionization_rate_cgs: 0 # HeII ionization rate in cgs [1/s]
RT_H2_dissociation_rate_cgs: 0 # H2 dissociation rate in cgs [1/s]
volumetric_heating_rates_cgs: 0 # Volumetric heating rate in cgs [erg/s/cm3]
specific_heating_rates_cgs: 0 # Specific heating rate in cgs [erg/s/g]
H2_three_body_rate : 1 # Specific the H2 formation three body rate (0->5,see Grackle documentation)
H2_cie_cooling : 0 # Enable/disable H2 collision-induced emission cooling from Ripamonti & Abel (2004)
H2_on_dust: 0 # Flag to enable H2 formation on dust grains
local_dust_to_gas_ratio : -1 # The ratio of total dust mass to gas mass in the local Universe (-1 to use the Grackle default value).
cmb_temperature_floor : 1 # Enable/disable an effective CMB temperature floor
initial_nHII_to_nH_ratio: -1 # initial nHII to nH ratio (number density ratio). Value is ignored if set to -1.
initial_nHeI_to_nH_ratio: -1 # initial nHeI to nH ratio (number density ratio). Value is ignored if set to -1.
initial_nHeII_to_nH_ratio: -1 # initial nHeII to nH ratio (number density ratio). Value is ignored if set to -1.
initial_nHeIII_to_nH_ratio: -1 # initial nHeIII to nH ratio (number density ratio). Value is ignored if set to -1.
initial_nDI_to_nH_ratio: -1 # initial nDI to nH ratio (number density ratio). Value is ignored if set to -1.
initial_nDII_to_nH_ratio: -1 # initial nDII to nH ratio (number density ratio). Value is ignored if set to -1.
initial_nHM_to_nH_ratio: -1 # initial nHM to nH ratio (number density ratio). Value is ignored if set to -1.
initial_nH2I_to_nH_ratio: -1 # initial nH2I to nH ratio (number density ratio). Value is ignored if set to -1.
initial_nH2II_to_nH_ratio: -1 # initial nH2II to nH ratio (number density ratio). Value is ignored if set to -1.
initial_nHDI_to_nH_ratio: -1 # initial nHDI to nH ratio (number density ratio). Value is ignored if set to -1.
.. note::
A simple example running SWIFT with Grackle can be find in ``examples/Cooling/CoolingBox``. A more advanced example combining heating and cooling (with heating and ionization sources) is given in ``examples/Cooling/CoolingHeatingBox``. ``examples/Cooling/CoolingWithPrimordialElements/`` runs a uniform cosmological box with imposed abundances and let them evolve down to redshift 0.
......@@ -6,293 +6,16 @@ GEAR model
===========
GEAR's model are mainly described in `Revaz \& Jablonka <https://ui.adsabs.harvard.edu/abs/2018A%26A...616A..96R/abstract>`_.
This model can be selected with the configuration option ``--with-subgrid=GEAR`` and run with the option ``--gear``. A few examples exist and can be found in ``examples/GEAR``.
This model can be selected with the configuration option ``--with-subgrid=GEAR`` and run with the option ``--gear``. A few examples exist and can be found in ``examples/GEAR``, ``examples/SinkParticles`` or ``IsolatedGalaxy/IsolatedGalaxy_multi_component``.
.. _gear_pressure_floor:
Pressure Floor
~~~~~~~~~~~~~~
.. toctree::
:caption: Table of Contents
In order to avoid the artificial collapse of unresolved clumps, a minimum in pressure is applied to the particles.
This additional pressure can be seen as the pressure due to unresolved hydrodynamics turbulence and is given by:
gear_model
chemistry
feedback
supernova_feedback
sinks/index
output
.. math::
P_\textrm{Jeans} = \frac{\rho}{\gamma} \frac{4}{\pi} G h^2 \rho N_\textrm{Jeans}^{2/3}
where :math:`\rho` is the density, :math:`\gamma` the adiabatic index, :math:`G` is the gravitational constant,
:math:`h` the kernel support and :math:`N_\textrm{Jeans}` (``GEARPressureFloor:jeans_factor`` in the parameter file) is the number of particle required in order to resolve a clump.
This must be directly implemented into the hydro schemes, therefore only a subset of schemes (Gadget-2, SPHENIX and Pressure-Energy) have the floor available.
In order to implement it, you need equation 12 in `Hopkins 2013 <https://arxiv.org/abs/1206.5006>`_:
.. math::
m_i \frac{\mathrm{d}v_i}{\mathrm{d}t} = - \sum_j x_i x_j \left[ \frac{P_i}{y_i^2} f_{ij} \nabla_i W_{ij}(h_i) + \frac{P_j}{y_j^2} f_{ji} \nabla_j W_{ji}(h_j) \right]
and simply replace the :math:`P_i, P_j` by the pressure with the floor (when the pressure is below the floor).
Here the :math:`x, y` are simple weights that should never have the pressure floor included even if they are related to the pressure (e.g. pressure-entropy).
.. code:: YAML
GEARPressureFloor:
jeans_factor: 10. # Number of particles required to suppose a resolved clump and avoid the pressure floor.
.. _gear_grackle_cooling:
Cooling: Grackle
~~~~~~~~~~~~~~~~
Grackle is a chemistry and cooling library presented in `B. Smith et al. 2017 <https://ui.adsabs.harvard.edu/abs/2017MNRAS.466.2217S>`_
(do not forget to cite if used). Four different modes are available:
equilibrium, 6 species network (H, H\\( ^+ \\), e\\( ^- \\), He, He\\( ^+ \\)
and He\\( ^{++} \\)), 9 species network (adds H\\(^-\\), H\\(_2\\) and
H\\(_2^+\\)) and 12 species (adds D, D\\(^+\\) and HD). Following the same
order, the swift cooling options are ``grackle_0``, ``grackle_1``, ``grackle_2``
and ``grackle_3`` (the numbers correspond to the value of
``primordial_chemistry`` in Grackle). It also includes some metal cooling (on/off with ``GrackleCooling:with_metal_cooling``), self-shielding
methods and UV background (on/off with ``GrackleCooling:with_UV_background``). In order to use the Grackle cooling, you will need
to provide a HDF5 table computed by Cloudy (``GrackleCooling:cloudy_table``).
When starting a simulation without providing the different element fractions in the non equilibrium mode, the code supposes an equilibrium and computes them automatically.
The code uses an iterative method in order to find the correct initial composition and this method can be tuned with two parameters. ``GrackleCooling:max_steps`` defines the maximal number of steps to reach the convergence and ``GrackleCooling:convergence_limit`` defines the tolerance in the relative error.
In order to compile SWIFT with Grackle, you need to provide the options ``with-chemistry=GEAR`` and ``with-grackle=$GRACKLE_ROOT``
where ``$GRACKLE_ROOT`` is the root of the install directory (not the ``lib``).
.. warning::
The actual Grackle version fully supported by SWIFT is 3.2.1. It can be downloaded from
`the official Grackle git repository <https://github.com/grackle-project/grackle/archive/refs/tags/grackle-3.2.1.tar.gz>`_.
However, this version still had a bug when using threadsafe functions. Alternately, it is possible to get a fixed version
using `the following fork frozen for compatibility with SWIFT <https://github.com/mladenivkovic/grackle-swift>`_.
To compile it, run
the following commands from the root directory of Grackle:
``./configure; cd src/clib``.
Update the variables ``LOCAL_HDF5_INSTALL`` and ``MACH_INSTALL_PREFIX`` in
the file ``src/clib/Make.mach.linux-gnu``.
Finish with ``make machine-linux-gnu; make && make install``.
Note that we require the 64 bit float version of Grackle, which should be the default setting.
(The precision can be set while compiling grackle with ``make precision-64``).
If you encounter any problem, you can look at the `Grackle documentation <https://grackle.readthedocs.io/en/latest/>`_
You can now provide the path given for ``MACH_INSTALL_PREFIX`` to ``with-grackle``.
In the parameters file, a few different parameters are available.
``GrackleCooling:redshift`` defines the redshift to use for the UV background (for cosmological simulation, it must be set to -1 in order to use the simulation's redshift) and ``GrackleCooling:provide_*_heating_rates`` can enable the computation of user provided heating rates (such as with the radiative transfer) in either volumetric or specific units.
For the feedback, it can be made more efficient by turning off the cooling during a few Myr (``GrackleCooling:thermal_time_myr``) for the particles touched by a supernovae.
The self shielding method is defined by ``GrackleCooling:self_shielding_method`` where 0 means no self shielding, > 0 means a method defined in Grackle (see Grackle documentation for more information) and -1 means GEAR's self shielding that simply turn off the UV background when reaching a given density (``GrackleCooling:self_shielding_threshold_atom_per_cm3``).
.. code:: YAML
GrackleCooling:
cloudy_table: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
with_UV_background: 1 # Enable or not the UV background
redshift: 0 # Redshift to use (-1 means time based redshift)
with_metal_cooling: 1 # Enable or not the metal cooling
provide_volumetric_heating_rates: 0 # (optional) User provide volumetric heating rates
provide_specific_heating_rates: 0 # (optional) User provide specific heating rates
max_steps: 10000 # (optional) Max number of step when computing the initial composition
convergence_limit: 1e-2 # (optional) Convergence threshold (relative) for initial composition
thermal_time_myr: 5 # (optional) Time (in Myr) for adiabatic cooling after a feedback event.
self_shielding_method: -1 # (optional) Grackle (1->3 for Grackle's ones, 0 for none and -1 for GEAR)
self_shielding_threshold_atom_per_cm3: 0.007 # Required only with GEAR's self shielding. Density threshold of the self shielding
HydrogenFractionByMass : 1. # Hydrogen fraction by mass (default is 0.76)
use_radiative_transfer : 1 # Arrays of ionization and heating rates are provided
RT_heating_rate_cgs : 0 # heating rate in units of / nHI_cgs
RT_HI_ionization_rate_cgs : 0 # HI ionization rate in cgs [1/s]
RT_HeI_ionization_rate_cgs : 0 # HeI ionization rate in cgs [1/s]
RT_HeII_ionization_rate_cgs: 0 # HeII ionization rate in cgs [1/s]
RT_H2_dissociation_rate_cgs: 0 # H2 dissociation rate in cgs [1/s]
volumetric_heating_rates_cgs: 0 # Volumetric heating rate in cgs [erg/s/cm3]
specific_heating_rates_cgs: 0 # Specific heating rate in cgs [erg/s/g]
H2_three_body_rate : 1 # Specific the H2 formation three body rate (0->5,see Grackle documentation)
H2_cie_cooling : 0 # Enable/disable H2 collision-induced emission cooling from Ripamonti & Abel (2004)
cmb_temperature_floor : 1 # Enable/disable an effective CMB temperature floor
.. note::
A simple example running SWIFT with Grackle can be find in ``examples/Cooling/CoolingBox``. A more advanced example combining heating and cooling (with heating and ionization sources) is given in ``examples/Cooling/CoolingHeatingBox``.
.. _gear_star_formation:
Star formation
~~~~~~~~~~~~~~
The star formation is done in two steps: first we check if a particle is in the star forming regime and then we use a stochastic approach to transform the gas particles into stars.
A particle is in the star forming regime if:
- The velocity divergence is negative (:math:`\nabla\cdot v < 0`),
- The temperature is lower than a threshold (:math:`T < T_t` where :math:`T_t` is defined with ``GEARStarFormation:maximal_temperature``),
- The gas density is higher than a threshold (:math:`\rho > \rho_t` where :math:`\rho_t` is defined with ``GEARStarFormation:density_threshold``)
- The particle reaches the pressure floor (:math:`\rho > \frac{\pi}{4 G N_\textrm{Jeans}^{2/3} h^2}\frac{\gamma k_B T}{\mu m_p}` where :math:`N_\textrm{Jeans}` is defined in the pressure floor).
If ``GEARStarFormation:star_formation_mode`` is set to ``agora``, the condition on the pressure floor is ignored. Its default value is ``default``.
A star will be able to form if a randomly drawn number is below :math:`\frac{m_g}{m_\star}\left(1 - \exp\left(-c_\star \Delta t / t_\textrm{ff}\right)\right)` where :math:`t_\textrm{ff}` is the free fall time, :math:`\Delta t` is the time step of the particle and :math:`c_\star` is the star formation coefficient (``GEARStarFormation:star_formation_efficiency``), :math:`m_g` the mass of the gas particle and :math:`m_\star` the mass of the possible future star. The mass of the star is computed from the average gas mass in the initial conditions divided by the number of possible stars formed per gas particle (``GEARStarFormation:n_stars_per_particle``). When we cannot have enough mass to form a second star (defined with the fraction of mass ``GEARStarFormation:min_mass_frac``), we fully convert the gas particle into a stellar particle. Once the star is formed, we move it a bit in a random direction and fraction of the smoothing length in order to avoid any division by 0.
Currently, only the following hydro schemes are compatible: SPHENIX and Gadget2.
Implementing the other hydro schemes is not complicated but requires some careful thinking about the cosmological terms in the definition of the velocity divergence (comoving vs non comoving coordinates and if the Hubble flow is included or not).
.. code:: YAML
GEARStarFormation:
star_formation_efficiency: 0.01 # star formation efficiency (c_*)
maximal_temperature: 3e4 # Upper limit to the temperature of a star forming particle
n_stars_per_particle: 4 # Number of stars that an hydro particle can generate
min_mass_frac: 0.5 # Minimal mass for a stellar particle as a fraction of the average mass for the stellar particles.
Chemistry
~~~~~~~~~
In the chemistry, we are using the smoothed metallicity scheme that consists in using the SPH to smooth the metallicity of each particle over the neighbors. It is worth to point the fact that we are not exchanging any metals but only smoothing it. The parameter ``GEARChemistry:initial_metallicity`` set the (non smoothed) initial mass fraction of each element for all the particles and ``GEARChemistry:scale_initial_metallicity`` use the feedback table to scale the initial metallicity of each element according the Sun's composition.
.. code:: YAML
GEARChemistry:
initial_metallicity: 1 # Initial metallicity of the gas (mass fraction)
scale_initial_metallicity: 1 # Should we scale the initial metallicity with the solar one?
Feedback
~~~~~~~~
The feedback is composed of a few different models:
- The initial mass function (IMF) defines the quantity of each type of stars,
- The lifetime of a star defines when a star will explode (or simply die),
- The supernovae of type II (SNII) defines the rates and yields,
- The supernovae of type Ia (SNIa) defines the rates and yields,
- The energy injection that defines how to inject the energy / metals into the particles.
Most of the parameters are defined inside a table (``GEARFeedback:yields_table``) but can be override with some parameters in the YAML file.
I will not describe theses parameters more than providing them at the end of this section.
Two different models exist for the supernovae (``GEARFeedback:discrete_yields``).
In the continuous mode, we integrate the quantities over the IMF and then explodes a floating point number of stars (can be below 1 in some cases).
In the discrete mode, we avoid the problem of floating points by rounding the number of supernovae (using a floor and randomly adding a supernovae depending on the fractional part) and then compute the properties for a single star at a time.
Initial mass function
^^^^^^^^^^^^^^^^^^^^^
GEAR is using the IMF model from `Kroupa (2001) <https://ui.adsabs.harvard.edu/abs/2001MNRAS.322..231K/abstract>`_.
We have a difference of 1 in the exponent due to the usage of IMF in mass and not in number.
We also restrict the mass of the stars to be inside :math:`[0.05, 50] M_\odot`.
Here is the default model used, but it can be easily adapted through the initial mass function parameters:
.. math::
\xi(m) \propto m^{-\alpha_i}\, \textrm{where}\,
\begin{cases}
\alpha_0 = 0.3,\, & 0.01 \leq m / M_\odot < 0.08, \\
\alpha_1 = 1.3,\, & 0.08 \leq m / M_\odot < 0.50, \\
\alpha_2 = 2.3,\, & 0.50 \leq m / M_\odot < 1.00, \\
\alpha_3 = 2.3,\, & 1.00 \leq m / M_\odot,
\end{cases}
Lifetime
^^^^^^^^
The lifetime of a star in GEAR depends only on two parameters: first its mass and then its metallicity.
.. math::
\log(\tau(m)) = a(Z) \log^2(m) + b(Z) \log(m) + c(Z) \\ \\
a(Z) = -40.110 Z^2 + 5.509 Z + 0.7824 \\
b(Z) = 141.929 Z^2 - 15.889 Z - 3.2557 \\
c(Z) = -261.365 Z^2 + 17.073 Z + 9.8661
where :math:`\tau` is the lifetime in years, :math:`m` is the mass of the star (in solar mass) and Z the metallicity of the star.
The parameters previously given are the default ones, they can be modified in the parameters file.
Supernovae II
^^^^^^^^^^^^^
The supernovae rate is simply given by the number of stars massive enough that end their life at the required time.
.. math::
\dot{N}_\textrm{SNII}(t) = \int_{M_l}^{M_u} \delta(t - \tau(m)) \frac{\phi(m)}{m} \mathrm{d}m
where :math:`M_l` and :math:`M_u` are the lower and upper mass limits for a star exploding in SNII, :math:`\delta` is the Dirac function and :math:`\phi` is the initial mass function (in mass).
The yields for SNII cannot be written in an analytical form, they depend on a few different tables that are based on the work of `Kobayashi et al. (2000) <https://ui.adsabs.harvard.edu/abs/2000ApJ...539...26K/abstract>`_ and `Tsujimoto et al. (1995) <https://ui.adsabs.harvard.edu/abs/1995MNRAS.277..945T/abstract>`_.
Supernovae Ia
^^^^^^^^^^^^^
The supernovae Ia are a bit more complicated as they involve two different stars.
.. math::
\dot{N}_\textrm{SNIa}(t) = \left( \int_{M_{p,l}}^{M_{p,u}} \frac{\phi(m)}{m} \mathrm{d}m \right) \sum_i b_i \int_{M_{d,l,i}}^{M_{d,u,i}}
\delta(t-\tau(m)) \frac{\phi_d(m)}{m}\mathrm{d}m
.. math::
\phi_d(m) \propto m^{-0.35}
where :math:`M_{p,l}` and :math:`M_{p,u}` are the mass limits for a progenitor of a white dwarf, :math:`b_i` is the probability to have a companion and
:math:`M_{d,l,i}` and :math:`M_{d,u,i}` are the mass limits for each type of companion.
The first parenthesis represents the number of white dwarfs and the second one the probability to form a binary.
+------------------+--------------------+-------------------+------------------+
| Companion | :math:`M_{d,l,i}` | :math:`M_{d,u,i}` | :math:`b_i` |
+==================+====================+===================+==================+
| Red giant | 0.9 | 1.5 | 0.02 |
+------------------+--------------------+-------------------+------------------+
| Main sequence | 1.8 | 2.5 | 0.05 |
+------------------+--------------------+-------------------+------------------+
The yields are based on the same papers than the SNII.
Energy injection
^^^^^^^^^^^^^^^^
All the supernovae (type II and Ia) inject the same amount of energy into the surrounding gas (``GEARFeedback:supernovae_energy_erg``) and distribute it according to the hydro kernel.
The same is done with the metals and the mass.
Generating a new table
^^^^^^^^^^^^^^^^^^^^^^
The feedback table is an HDF5 file with the following structure:
.. graphviz:: feedback_table.dot
where the solid (dashed) squares represent a group (a dataset) with the name of the object underlined and the attributes written below. Everything is in solar mass or without units (e.g. mass fraction or unitless constant).
In ``Data``, the attribute ``elts`` is an array of string with the element names (the last should be ``Metals``, it corresponds to the sum of all the elements), ``MeanWDMass`` is the mass of the white dwarfs
and ``SolarMassAbundances`` is an array of float containing the mass fraction of the different element in the sun.
In ``IMF``, ``n + 1`` is the number of part in the IMF, ``as`` are the exponent (``n+1`` elements), ``ms`` are the mass limits between each part (``n`` elements) and
``Mmin`` (``Mmax``) is the minimal (maximal) mass of a star.
In ``LifeTimes``, the coefficient are given in the form of a single table (``coeff_z`` with a 3x3 shape).
In ``SNIa``, ``a`` is the exponent of the distribution of binaries, ``bb1`` and ``bb2`` are the coefficient :math:`b_i` and the other attributes follow the same names than in the SNIa formulas.
The ``Metals`` group from the ``SNIa`` contains the name of each elements (``elts``) and the metal mass fraction ejected by each supernovae (``data``) in the same order. They must contain the same elements than in ``Data``.
Finally for the ``SNII``, the mass limits are given by ``Mmin`` and ``Mmax``. For the yields, the datasets required are ``Ej`` (mass fraction ejected [processed]), ``Ejnp`` (mass fraction ejected [non processed]) and one dataset for each element present in ``elts``. The datasets should all have the same size, be uniformly sampled in log and contains the attributes ``min`` (mass in log for the first element) and ``step`` (difference of mass in log between two elements).
.. code:: YAML
GEARFeedback:
supernovae_energy_erg: 0.1e51 # Energy released by a single supernovae.
yields_table: chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5 # Table containing the yields.
discrete_yields: 0 # Should we use discrete yields or the IMF integrated one?
GEARInitialMassFunction:
number_function_part: 4 # Number of different part in the IMF
exponents: [0.7, -0.8, -1.7, -1.3] # Exponents of each part of the IMF
mass_limits_msun: [0.05, 0.08, 0.5, 1, 50] # Limits in mass between each part of the IMF
GEARLifetime:
quadratic: [-40.1107, 5.50992, 0.782432] # Quadratic terms in the fit
linear: [141.93, -15.8895, -3.25578] # Linear terms in the fit
constant: [-261.366, 17.0735, 9.86606] # Constant terms in the fit
GEARSupernovaeIa:
exponent: -0.35 # Exponent for the distribution of companions
min_mass_white_dwarf_progenitor: 3 # Minimal mass of a progenitor of white dwarf
max_mass_white_dwarf_progenitor: 8 # Maximal mass of a progenitor of white dwarf
max_mass_red_giant: 1.5 # Maximal mass for a red giant
min_mass_red_giant: 0.9 # Minimal mass for a red giant
coef_red_giant: 0.02 # Coefficient for the distribution of red giants companions
max_mass_main_sequence: 2.6 # Maximal mass for a main sequence star
min_mass_main_sequence: 1.8 # Minimal mass for a main sequence star
coef_main_sequence: 0.05 # Coefficient for the distribution of main sequence companions
white_dwarf_mass: 1.38 # Mass of a white dwarf
GEARSupernovaeII:
interpolation_size: 200 # Number of elements for the interpolation of the data