Running an isolated elliptical galaxy with all the available COLIBRE physics

The default IC for the Colibre example is the Mh13_c9_M5_3kpc.hdf5 IC file. 

#####################################################################

How are the ICs setup:
For the dark matter and the stars we use a Hernquist profile, in the case of 
the stellar Hernquist profile, this is generated using MakeNewDisk from 
Springel+ (2005), the Hernquist stellar bulge is setup assuming that the 
we have a dark matter Hernquist profile. In our runs we use a static dark 
matter profile to save in computational cost. When we setup the stars + 
dark matter we setup the CGM around our elliptical galaxy using hydrostatic
equilibrium using the hydrostatic relation (following Stern+ (2019)):

d ln P / d ln r = - \gamma v_c^2 / c_s^2. (1)

We assume a few things when calculating the pressure profile, firstly we 
assume that the gas has a negligible effect on the mass of the system. 
This means that Eq. (1) is scale invariant in terms of the density of the 
gas. It is not invariant under a temperature change because the sound 
speed depends on the temperature of the gas. In our case we also assume 
that the ratio between the circular velocity and the sound speed is equal 
to unity at large radii, e.g.:

v_c^2 = c_s^2, 

-> T = \mu * m_p / (k_B * \gamma) * G M_encl (r) / r. (2)

Using the assumption of Eq. (2), it is easy to solve the differential 
equation of Eq. (1):

P \propto r^{-\gamma}. (3)

Than the model only requires an additional boundary density to get the 
correct normalization for the pressure and density, for this we use the 
gas fraction at R500 to get the correct enclosed mass at this radius. 
We determine R500 when no gas is yet added to the IC and check that R500 
changes less than a few percentage, typical lower than 0.5% when we include 
gas. The model we use for the correct gas fraction within R500 is the best
fit following Debackere+ (2019).

In general the assumption of always having that the temperature is given by
the rotational velocity is incorrect especially in the center of the halo,
because the galaxy is here and feedback processes will make this assumption
invalid and impractical because of a few drawbacks:
1. We could get zero initial temperature because the rotational velocity is
   approaching zero. 
2. We will get a cooling catastrophy because the temperature at the center 
   is very small, and therefore get a big burst of initial star formation.

To prevent the problem of having a minimum temperature it is possible to 
specify in the code a minimum temperature in the center. This temperature 
is than added on top of the temperature based on the circular velocity,
in general we use a s-curve to add the additional temperature to the 
circular velocity:

T_additional = Tmin / (1 + exp((r-a*rmin)/rmin)), a = 2.5 (4)

In general for the 10^13 halo Tmin should be around 1e6 and rmin=15 kpc. 

In the case of adding this additional temperature it is no longer possible
to find a simple analytical solution for the pressure therefore we solve 
the differential equation numerical, the solution reduces at large radii to
Eq. (3). 

The parameters of the central temperature can be changed in the file as 
variable:
Tmin_model = 1e6    # K (central minimum temperature)
Tscale_length = 15. # kpc (length the core temperature decreases)

#####################################################################

Angular momentum distribution in the gas:
For the angular momentum distribution in the galaxy we follow Bullock+ (2001)
Following Bullock+ (2001) we first find the normalization of the angular 
momentum which is given by:

j0 = sqrt(2) * V200 * R200 * \lambda / (- \mu * \ln(1 - 1/\mu) - 1. ) (5)

In which \lambda is the spin parameter, \mu is a property of the angular 
momentum profile which is on average \mu = 1.25. Using this we assume that 
the specific angular momentum is lower closer to the center because otherwise
we do get velocities that increase by 1/r when we get closer to the center.
We use than:

M(<r) = M200 * \mu * j(r) / (j0 + j(r)) (6)

Rewriting to find j(r):

j = j0 * M(<r)/M200 / ( \mu - M(<r)/M200 ). (7)

This than gives us the specific angular momentum inside R200, because the work
of Bullock+ (2001) does not go farther than R200. For radii larger than R200,
we assume that the specific angular momentum remains constant, implying the 
velocity will decrease as 1/r.

Now that we know the angular momentum for each particle it is possible to 
calculate the velocity (j = v x r):

v = j/r. (8)

#####################################################################

Running with different gas fractions:
In general the gas fraction that is used is the mean gas fraction for the used
M500 mass, however it is possible to run with a certain amount of standard 
deviations higher or lower gas mass if desired. This can be done by changing 
the parameter amount_of_sigma_used in the Debackere model section of the code.
This basically allows the user to use 1 or 2 sigma lower or higher gas fraction
as is based on the best fit. For this it is assumed that the uncertainty in the 
gas fraction is given by 0.02. 

Note: If for a low mass halo -1 or -2 sigma is used and the best fit is already
lower than this factor times the uncertainty this will produce a negative mass,
in this case the behaviour of the code is unclear. For example: Mh = 10^13 and 
amount_of_sigma_used = -2 will cause a problem.

#####################################################################

Boundary of IC gas:

Small radius:
For the gas profile we use a minimum radius of 1 kpc, the main reason for this 
is to prevent the inner profile to get higher densities near the center, in the 
case that we get high densities in the center we often start the simulation with 
a very large AGN heating event that quenches the galaxy for a long time, because
of this we start the run with a gap in the density profile such that no gas is 
present within 1 kpc, this will allow the simulation to start more gentle in 
terms of the AGN feedback. 

Large radius:
The current version of the IC generator only places gas within the minimum of 
1 Mpc or 3 * R200, this is to prevent the simulation from containing too much 
gas particles, because things that happen at very large radii are not very 
important for the inner revolution for at least the first few Gyr.

#####################################################################

Running the IC generator code:
The IC generator script requires that you give an IC file that contains the 
stellar bulge. Using this the IC generator can be called as:

python3 calculate_elliptical_CGM.py input_HDF5_file.hdf5 output_HDF5_file.hdf5 mass_scheme

For the default example we recommend to use:
python3 calculate_elliptical_CGM.py Mh13_c9_M5_3kpc.hdf5 Mh13_c9_M5_3kpc_CGM.hdf5 default

The last argument is the mass scheme, available schemes are default, linear1,
linear5, linear10 and quadratic.

#####################################################################

Available mass schemes:
In order to do large haloes we need to be able to run the code with larger mass 
at larger radii from the center, this is implemented in the IC generator code,
if we do not want a variable mass, the default scheme will be used, that uses the
same gas mass as the stellar mass in the input hdf5 file. If we want to run the 
code with a variable mass at large radii (r > 100 kpc) we can use 4 different 
schemes that have a mass dependence as:

Mres = Mmin *(  a*(r/100 kpc)^n - b ), b = a-1 (9)

The available schemes have:
linear1 : a = 1, n = 1
linear5 : a = 5, n = 1 
linear10 : a = 10, n = 1
quadratic: a = 1, n=2

Amount of memory saved when using the mass scheme (snapshots):
Fid : 5.2 Gb
linear1: 1.3 Gb
linear10: 0.465 Gb
quadratic: 0.622 Gb

Runtime impact for first 500 Myr (a bit uncertain):
Fid:     ~11 hours
linear1: ~6 hours
others:  ~5 hours

#####################################################################

Other parameters that can be modified:
mBH = 10^8.5          # Black hole mass
lambda_angular = 0.05 # the spin parameter of the gas at R200
mu_angular = 1.25     # the Bullock+ (2001) free parameter of the angular 
                      # momentum profile

#####################################################################

Chosen default parameters galaxy:
For the default elliptical galaxy we have chosen the most average properties
of an elliptical galaxy at a total mass of 10^13 Msun. We have chosen to use
a stellar mass of 10^11 Msun, a stellar effective size of 3 kpc, and no 
angular momentum in the CGM. The other properties of the galaxy are fixed by
the above assumptions (e.g. gas distribution and total mass is not a free 
parameter under our assumptions)

The default ICs:
Name IC file:             eps:        Mstar:
Mh13_c9_M5_3kpc.hdf5      0.2 (kpc)   10^5 Msun
Mh13_c9_M6_3kpc.hdf5      0.4 (kpc)   10^6 Msun
Mh13_c9_M4_3kpc.hdf5      0.1 (kpc)   10^4 Msun

The ICs with a more compact bulge:
Name IC file:             eps:        Mstar:
Mh13_c9_M5_2kpc.hdf5      0.2 (kpc)   10^5 Msun
Mh13_c9_M6_2kpc.hdf5      0.4 (kpc)   10^6 Msun
Mh13_c9_M4_2kpc.hdf5      0.1 (kpc)   10^4 Msun

The ICs with a more diffuse bulge:
Name IC file:             eps:        Mstar:
Mh13_c9_M5_4.5kpc.hdf5      0.2 (kpc)   10^5 Msun
Mh13_c9_M6_4.5kpc.hdf5      0.4 (kpc)   10^6 Msun
Mh13_c9_M4_4.5kpc.hdf5      0.1 (kpc)   10^4 Msun

#####################################################################

Compiling the code and black hole behaviour:
In general the code can be run with a mobile black hole that can move around
and do reposition if wanted, however sometimes it is desired to place the 
black hole exactly at the center, this is possible by setting the black hole 
as a fixed boundary particle, e.g. a particle that remains fixed under all 
conditions, in which it cannot receive acceleration or a velocity due to 
conservation laws. In the generated IC files the black hole by construction 
has the lowest ID, such that we can easily turn the black hole into a fixed
boundary particle.

Recommended config option mobile BH:
./configure --with-subgrid=COLIBRE --with-hydro=sphenix --with-kernel-wendland-C2 --with-ext-potential=hernquist --with-tbbmalloc --enable-ipo

Recommended config option fixed black hole:
./configure --with-subgrid=COLIBRE --with-hydro=sphenix --with-kernel-wendland-C2 --with-ext-potential=hernquist --with-tbbmalloc --enable-ipo --enable-fixed-boundary-particles=2