Add nuclear star cluster component

[liempi]

Add nuclear star cluster model and BH seed formation scenario within NSCs due to stellar collisions.

[abensonca]

Adding comments from Matias here for reference:

1. First, in the objects.nodes.components.black_hole.standard.F90, the Node_Component_Black_Hole_Standard_Scale_Set subroutine has 1 M⊙ for the mass scale, it is because otherwise it leads to negative masses. I think it would be nice to introduce a new relative mass scale with the stellar mass of the NSC, maybe when the NSC already formed a seed.
2. The BH seed formed under the model is created in the nodes.operators.physics.stellarCollisions.NSCs.F90 file. I tried to do move this node operator to the blackHoleSeedClass and it is not possible because as there are not BHs in all the galaxies, it ends in an error in routines using the potential in their calculations (as the mass of the BH is 0.0).
3. In the objects.nodes.components.black_hole.simple.F90, I'm not sure if to introduce a modification. In the rate task, in principle we should remove the gas accreted also for the hot halo as in the standard file but it is not implemented. If it makes sense, I can do it.

and not related to my projects, I think there is a mistake in the cooling.cooling_radius.beta_profile.F90 file. In line 381, the variable is called densityZero which is supposed to be the density evaluated at radius 0.0, but it is using the coordinates of the outerRadius of the Halo (previously used to compute the density at the outer radius).

[abensonca]

and not related to my projects, I think there is a mistake in the cooling.cooling_radius.beta_profile.F90 file. In line 381, the variable is called densityZero which is supposed to be the density evaluated at radius 0.0, but it is using the coordinates of the outerRadius of the Halo (previously used to compute the density at the outer radius).

Yes - definitely a mistake. I opened issue #781 for this and will push a fix.

[abensonca]

1. First, in the `objects.nodes.components.black_hole.standard.F90`, the `Node_Component_Black_Hole_Standard_Scale_Set` subroutine has 1 M⊙ for the mass scale, it is because otherwise it leads to negative masses. I think it would be nice to introduce a new relative mass scale with the stellar mass of the NSC, maybe when the NSC already formed a seed.

This seems reasonable. I'd suggest making it some fraction of the NSC stellar mass.

[abensonca]

2. The BH seed formed under the model is created in the `nodes.operators.physics.stellarCollisions.NSCs.F90` file. I tried to do move this node operator to the `blackHoleSeedClass` and it is not possible because as there are not BHs in all the galaxies, it ends in an error in routines using the potential in their calculations (as the mass of the BH is 0.0).

I caught this problem in the PR refactoring black hole physics, and it is now implemented correctly there. We'll need to be sure to move the equivalent NSC behavior to this new class (after I merge in that PR).

[abensonca]

2. The BH seed formed under the model is created in the `nodes.operators.physics.stellarCollisions.NSCs.F90` file. I tried to do move this node operator to the `blackHoleSeedClass` and it is not possible because as there are not BHs in all the galaxies, it ends in an error in routines using the potential in their calculations (as the mass of the BH is 0.0).

It would be better to move this to the blackHoleSeedClass if possible. If the only issue is the black hole potential, then we should change the massDistributionBlackHole class so that it correctly handles cases of zero mass black holes (just returning a zero potential in such cases). I'm actually surprised this happens as in the Node_Component_Black_Hole_Standard_Mass_Distribution function it checks that the black hole mass is $>0$ - and if it isn't returns a null massDistribution - in which case this problem should never occur. But maybe that's not working for some reason.

[abensonca]

Hi Matias: I finished reviewing the PR. You've done an impressive amount of work on this, and the code looks good and clean overall. I also merged in the recent changes from the master branch which move some of the black hole physics out of the nodeComponentBlackHole classes into nodeOperators and associated functions. This didn't require too much modification on your branch, so I suspect it will work successfully.

One important change. In the new NSC component, I changed:

   <isDefault>true</isDefault>

to

   <isDefault>false</isDefault>

This will cause Galacticus to default to a null NSC component (i.e. no NSC), which means that the default behavior will be unchanged compared to earlier versions (without this, some models would break as they wouldn't have the correct parameters set for NSCs, and all models would have their results changed because they would now include NSC physics when previously the didn't).

I've left numerous review comments on the PR in GitHub. Some of these are quite minor (e.g. naming a constant instead of using a value directly, adding a citation, an unused property that could be removed, etc.). A few are pointing out possible bugs (e.g. wrong calculation of angular momentum, failure of mass conservation, etc.). Then a few are more structural - suggesting moving some behavior into separate classes, or making some property a meta-property. I'm happy to discuss these if anything isn't clear, but a couple of notes:

• A "meta-property" is some property associated with a nodeComponent, but which isn't defined in that nodeComponent directly. Instead, it's attached to the nodeComponent by some other class (often a nodeOperator) which needs to store some data in the nodeComponent in order to do its calculation. Given the modular nature of Galacticus, this is useful as it means we don't have to carry around every property that could ever be needed by every module - we only create and use them if a module is actively used in a given model. This save memory and CPU time.

• When deciding what functionality to put into which classes, I think about it this way. Most behavior in Galacticus can be written as a differential equation. For example, star formation in the NSC could be written as:

$$ \mathrm{d}M_\star/\mathrm{d}t = +\phi(\cdots) ; \mathrm{d}M_\mathrm{g}/\mathrm{d}t = -\phi(\cdots) $$

where $M_\star$ is the stellar mass of the NSC, $M_\mathrm{g}$ is the gas mass of the NSC, and $\phi(\cdots)$ is the rate of star formation (dependent on some set of properties of the galaxy indicated by " $\cdots$ ").

There are three pieces to this:

1. The properties involved, $M_\star$ and $M_\mathrm{g}$;
2. The equations defining how this process (star formation) moves mass (or other quantities) between properties;
3. The rate of the process, $\phi(\cdots)$.

In general:

1. The properties themselves are defined in a nodeComponent;
2. The equations are defined in a nodeOperator;
3. The rate of the process is defined in its own class (e.g. starFormationRateNuclearStarClusters in this case).

Currently, for the nodeOperatorstellarCollisionsNSC and nodeOperatorgasMassRateNSC classes, the rates are also being implemented within the nodeOperator. It would be better to move these into their own classes (so, more like how nodeOperatorStarFormationNuclearStarClusters works by utilizing the starFormationRateNuclearStarClusters class).

This will make it easier to use different physical models for the rates of these processes.

[abensonca]

I only see this property set - it is never used. I assume the intent here is that it is output as a quantity of interest? If so, it might be better to make it a meta-properties (so it is only present in the output if requested). The same applies to the massCritical property below.

[abensonca]

Define 1.0d0 as a parameter here, e.g.:

double precision, parameter :: surfaceDensityGasNormalization=1.0d0 ! M☉/pc²

to make it more clear what we're doing - otherwise dividing by 1 seems pointless!

[abensonca]

This function is essentially identical to the one in source/star_formation.rate_surface_density.disks.Krumholz2009.F90. To avoid repeating code, we could move this function to a new module (in a new file), called Star_Formation_Rate_Krumholz2009_Utilities or something similar - and then have both of these classes use the function from that new module.

[abensonca]

These calculations might also be good to put into a separate utility module (see below) since they are repeated in the disks Krumholz2009 star formation rate class.

[abensonca]

I moved metallicityRelativeToSolar, chi, and s out of the class and made them local to this function, as they didn't need to be in the class directly. (In the disk Krumholz2009 class these values are memorized for re-use on subsequent calls, so it makes sense to define them in the class there).