initialize_PBHs_from_galaxy_potential

eridanus_pbh.scripts.WithPBHs.PBH_in_galaxy.initialize_PBHs_from_galaxy_potential(total_potential: (<class 'galpy.potential.Potential.Potential'>, <class 'galpy.potential.amuse.galpy_profile'>), mass_fraction: float, imf_func: (unit<MSun>, <class 'function'>) = quantity<1000 MSun>, imf_args: (unit<MSun>, <class 'list'>) = quantity<10 MSun>, imf_kwargs: dict = {}, radius_cutoff: unit<kpc> = quantity<100 kpc>, distr_args: list = [], distr_kwargs: dict = {}, position: unit<kpc> = quantity<[0, 0, 0] kpc>, velocity: unit<kms> = quantity<[0, 0, 0] kms>, obj_radius: unit<AU> = quantity<0 AU>, evln_func: (<class 'bool'>, <class 'function'>) = False, evln_kwargs: dict = {}, gravity_func: (<class 'bool'>, <class 'function'>) = <class 'amuse.community.bhtree.interface.BHTree'>, gravity_args: list = [], gravity_kwargs: dict = {}, smoothing_length: unit<parsec> = quantity<1.0 parsec>, opening_angle: float = 0.6, number_of_workers: int = 8, use_self_gravity=False, timestep: unit<Myr> = quantity<1.0 Myr>, random=True, channel_attrs=None, logger=<utilipy.utils.logging._LogPrint.LogPrint object at 0x120dd6590>, verbose=None)[source]

Create PBHs in a galaxy.

Parameters
total_potential: :class:`~galpy.potential.amuse.galpy_profile`

The original potential

mass_fraction: float
imf_func: amuse quantity or Callable

initial mass function if a quantity, then treated as average mass per pbh

imf_args: list

arguments into imf_func if imf_func is a quantity, then scale in random.normal

imf_kwargs: dict

key-word arguments into imf_func not used if imf_func is a quantity

radius_cutoff: amuse distance quantity

cutoff radius for calculating enclosed mass

distr_args: list, optional

the arguments for distr_func (default [])

distr_kwargs: dictionary, optional

the kwargs for distr_func (default dict())

Rvirial: Quantity, optional

the virial radius of the system (default 10 pc)

position: VectorQuantity, optional

the position of the system in GC coordinates (default [0,0,0] pc)

velocity: VectorQuantity, optional

the velocity of the system in GC coordinates (default [0, 0, 0] kms)

obj_radius: Quantity, optional

the radius of the individual objects (default 0 AU)

evln_func: function or False, optional

object evolution function (default False) signature of function should be

func(metallicity, **evln_kwargs)
metallicity: float, optional

the metallicity parameter for the evolution ex SeBa

imf_kwargs: dictionary, optional

the kwargs for evln_func

gravity_func: function or False, optional

gravity code (default False) signature of function should be

func(converter, *gravity_args,
     number_of_workers=number_of_workers,
     **gravity_kwargs)
gravity_args: list, optional

the arguments for gravity_func

gravity_kwargs: dictionary, optional

the kwargs for gravity_func

smoothing_length: distance quantity, optional

the smoothing length used in scaling and the gravity (default 0 pc)

opening_angle: float, optional

the gravity code opening angle (default 0.6)

number_of_workers: int, optional

number of gravity workers (default 1)

timestep: time quantity, optional

the timestep for evolving the gravity and stellar evolution (default 1 Myr)

random: True or random number generator, optional

(default True) ex np.random.default_rng(seed=0) will default to random seed

channel_attrs: list or None, optional

(default None)

Returns
system: datamodel.System

a dataclass object with parameters

  • particles

  • evolution

  • gravity

  • converter

  • channel_attrs

will try to automatically make channels to/from all things the amuse particles / evolution / gravity classes are proxied in a datamodel.Container that adds .name, .channel_to/from

adj_potential: amuse potential

adjusted smooth potential since the PBHs are part of the total_potential

Notes

For calculating the velocity distribution

  1. Assuming the PBH contribute smoothly to the potential, not as individual masses

  2. Assuming the system is virialized, then 2 KE = -PE.

  3. \(v^2 \sim \sigma^2\)

  4. the velocities are a perfect Gaussian (related to the above).

Putting these together gives \(\sigma(r)\) at any r. For each PBH, drawing 3 velocities for the sigmas.

Can use the total, not adjusted, potential.

Method: \(v^2 = - \rm{Potential}\) and \(v^2 = \sqrt{3} \sigma^2\) This gets the magnitude. For the direction draw a Gaussian for x, y, z, normalize to unit magnitude, and multiply by \(\sigma\)