holodeck.sams.sam

Semi Analytic Modeling (SAM) submodule.

The core element of the SAM module is the Semi_Analytic_Model class. This class requires four components as arguments:

1. Galaxy Stellar Mass Function (GSMF): gives the comoving number-density of galaxies as a function of stellar mass. This is implemented as subclasses of the _Galaxy_Stellar_Mass_Function base class.

2. Galaxy Pair Fraction (GPF): gives the fraction of galaxies that are in a ‘pair’ with a given mass ratio (and typically a function of redshift and primary-galaxy mass). Implemented as subclasses of the _Galaxy_Pair_Fraction subclass.

3. Galaxy Merger Time (GMT): gives the characteristic time duration for galaxy ‘mergers’ to occur. Implemented as subclasses of the _Galaxy_Merger_Time subclass.

4. M_bh - M_bulge Relation (mmbulge): gives MBH properties for a given galaxy stellar-bulge mass. Implemented as subcalsses of the holodeck.host_relations._MMBulge_Relation subclass.

The Semi_Analytic_Model class defines a grid in parameter space of total MBH mass (\(M=M_1 + M_2\)), MBH mass ratio (\(q \\equiv M_1/M_2\)), redshift (\(z\)), and at times binary separation (semi-major axis \(a\)) or binary rest-frame orbital-frequency (\(f_r\)). Over this grid, the distribution of comoving number-density of MBH binaries in the Universe is calculated. Methods are also provided that interface with the kalepy package to draw ‘samples’ (discretized binaries) from the distribution, and to calculate GW signatures.

The step of going from a number-density of binaries in \((M, q, z)\) space, to also the distribution in \(a\) or \(f\) is subtle, as it requires modeling the binary evolution (i.e. hardening rate).

To-Do (sam.py)

• Allow SAM class to take M-sigma in addition to M-Mbulge.

References

• [Sesana2008] Sesana, Vecchio, Colacino 2008.

• [Chen2019] Chen, Sesana, Conselice 2019.

holodeck.sams.sam.REDZ_SAMPLE_VOLUME = True

get redshifts by sampling uniformly in 3D spatial volume, and converting

holodeck.sams.sam.GSMF_USES_MTOT = False

the mass used in the GSMF is interpretted as M=m1+m2, otherwise use primary m1

holodeck.sams.sam.GPF_USES_MTOT = False

the mass used in the GPF is interpretted as M=m1+m2, otherwise use primary m1

holodeck.sams.sam.GMT_USES_MTOT = False

the mass used in the GMT is interpretted as M=m1+m2, otherwise use primary m1

class holodeck.sams.sam.Semi_Analytic_Model(mtot=(1.988409870698051e+37, 1.988409870698051e+45, 91), mrat=(0.001, 1.0, 81), redz=(0.001, 10.0, 101), shape=None, log=None, gsmf=<class 'holodeck.sams.components.GSMF_Schechter'>, gpf=None, gmt=None, gmr=None, mmbulge=<class 'holodeck.host_relations.MMBulge_KH2013'>, **kwargs)

Bases: object

Semi-Analytic Model (SAM) of MBH Binary populations.

This class produces simulated MBH binary populations using idealized (semi-)analytic functions starting from galaxy populations, to massive black holes, to merger rates. Using SAMs, MBH binary populations are calculated over a fixed, rectilinear grid of 3 or 4 dimensional parameter-space. The starting parameter space is total mass, mass ratio, and redshift; but often the distribution of binaries are desired at particular orbital frequencies or separations, which adds a dimension. Some parameters are calculated at grid edges (e.g. binary number densities), while others are calculated at grid centers (e.g. number of binaries in a universe). Ultimately, what the SAM calculates in the number (or number-density) of binaries at each point in this 3-4 dimensional parameter space.

Conceptually, three components are required to build SAMs.

1. Galaxies and stellar-masses (i.e. how many galaxies there are as a function of stellar mass and redshift). This component is provided by the Galaxy Stellar-Mass Function (GSMF), which are implemented as subclasses of the holodeck.sams.components._Galaxy_Stellar_Mass_Function base-class.

2. Galaxy merger rates (GMRs; i.e. how often galaxies merge as a function of stellar mass, mass ratio, and redshift.) This component can be provided in two ways:

   • Subclasses of holodeck.sams.components._Galaxy_Merger_Rate, which provide merger rates directly.

   • Both a galaxy pair fraction (GPF; subclasses of holodeck.sams.components._Galaxy_Pair_Fraction) which give the fraction of galaxies in the process of merger, and a galaxy merger time (GMT; subclasses of _Galaxy_Merger_Time) which gives the duration of time that galaxies spend in the merger process.

3. MBH-Host relationships which determine MBH properties for a given host galaxy. Currently these relationships are only fully implemented as Mbh-MBulge (MMBulge) relationships, which are subclasses of holodeck.host_relations._MMBulge_Relation.

__init__(mtot=(1.988409870698051e+37, 1.988409870698051e+45, 91), mrat=(0.001, 1.0, 81), redz=(0.001, 10.0, 101), shape=None, log=None, gsmf=<class 'holodeck.sams.components.GSMF_Schechter'>, gpf=None, gmt=None, gmr=None, mmbulge=<class 'holodeck.host_relations.MMBulge_KH2013'>, **kwargs)

Construct a new Semi_Analytic_Model instance.

Parameters:

• mtot ((3,) tuple) – Specification for the domain of the grid in total-mass. Three arguments must be included, the 0th giving the lower-limit [grams], the 1th giving the upper-limit [grams], and the 2th giving the number of bin-edges (i.e. the number-of-bins plus one).

• mrat ((3,) tuple) – Specification for the domain of the grid in mass-ratio. Three arguments must be included, the 0th giving the lower-limit, the 1th giving the upper-limit, and the 2th giving the number of bin-edges (i.e. the number-of-bins plus one).

• redz ((3,) tuple) – Specification for the domain of the grid in redshift. Three arguments must be included, the 0th giving the lower-limit, the 1th giving the upper-limit, and the 2th giving the number of bin-edges (i.e. the number-of-bins plus one).

• shape (int or (3,) tuple) –

  The shape of the grid in total-mass, mass-ratio, and redshift. This argument specifies the number of grid-edges in each dimension, and overrides the shape arguments of mtot, mrat, and redz. * If a single int is given, then this is the shape applied to all dimensions. * If a (3,) iterable of values is given, then each value specifies the size of the grid

  in the corresponding dimension. None values can be provided which indicate to use the default sizes (provided by the mtot, mrat, and redz arguments.) For example, shape=(12, None, 14) would produce 12 grid edges in total mass, the default number of grid edges in mass ratio, and 14 grid edges in redshit.

• gsmf (None or _Galaxy_Stellar_Mass_Function subclass instance) –

• gpf (None or _Galaxy_Pair_Fraction subclass instance) –

• gmt (None or _Galaxy_Merger_Time subclass instance) –

• gmr (None or _Galaxy_Merger_Rate subclass instance) –

• mmbulge (None or _MMBulge_Relation subclass instance) –

_gsmf

Galaxy Stellar-Mass Function (_Galaxy_Stellar_Mass_Function instance)

_mmbulge

Mbh-Mbulge relation (host_relations._MMBulge_Relation instance)

_gpf

Galaxy Pair Fraction (_Galaxy_Pair_Fraction instance)

_gmt

Galaxy Merger Time (_Galaxy_Merger_Time instance)

_gmr

Galaxy Merger Rate (_Galaxy_Merger_Rate instance)

_density

Binary comoving number-density

_shape

Shape of the parameter-space domain (mtot, mrat, redz)

_gmt_time

GMT timescale of galaxy mergers [sec]

_redz_prime

redshift following galaxy merger process

property edges

[mtot, mrat, redz])

Type:

The grid edges defining the domain (list of

property shape

Shape of the parameter space domain (number of edges in each dimension), (3,) tuple

mass_stellar()

Calculate stellar masses for each MBH based on the M-MBulge relation.

Returns:

masses – Galaxy total stellar masses for all MBH. [0, :] is primary, [1, :] is secondary [grams].

Return type:

(2, N) ndarray of scalar,

property static_binary_density

The number-density of binaries at the edges of the grid in mtot, mrat, redz.

The ‘number density’ is a density both in terms of volume (i.e. number of binaries per unit comoving-volume, \(n = dN/dV_c\)), and in terms of binary parameters (e.g. binaries per unit of log10 mass, \(d n /d \log_{10} M\)). Specifically, the values returned are:

\[d^3 n / [d \log_{10} M d q d z]\]

For each Semi_Analytic_Model instance, this value is calculated once and cached.

Returns:

density – Number density of binaries, per unit redshift, mass-ratio, and log10 of mass. The values are in units of inverse cubic, comoving-Mpc.

Return type:

(M, Q, Z) ndarray [\(cMpc^{-3}\)]

Notes

• This function effectively calculates Eq.21 & 5 of [Chen2019]; or equivalently, Eq. 6 of [Sesana2008].

• Bins which ‘merge’ after redshift zero are set to zero density (using the self._gmt instance).

dynamic_binary_number_at_fobs(hard, fobs_orb, **kwargs)

_dynamic_binary_number_at_fobs_consistent(hard, fobs_orb, steps=200, details=False)

Get correct redshifts for full binary-number calculation.

Slower but more correct than old dynamic_binary_number. Same as new cython implementation sam_cyutils.dynamic_binary_number_at_fobs, which is more than 10x faster. LZK 2023-05-11

# BUG doesn’t work for Fixed_Time_2PL

_dynamic_binary_number_at_fobs_inconsistent(hard, fobs_orb)

_dynamic_binary_number_at_sepa_consistent(hard, target_sepa, steps=200, details=False)

Get correct redshifts for full binary-number calculation.

Slower but more correct than old dynamic_binary_number. Same as new cython implementation sam_cyutils.dynamic_binary_number_at_fobs, which is more than 10x faster. LZK 2023-05-11

gwb_new(fobs_gw_edges, hard=<holodeck.hardening.Hard_GW object>, realize=100)

Calculate GWB using new cython implementation, 10x faster!

gwb_old(fobs_gw_edges, hard=<class 'holodeck.hardening.Hard_GW'>, realize=100)

Calculate GWB using new dynamic_binary_number_at_fobs method, better, but slower.

gwb_ideal(fobs_gw, sum=True, redz_prime=None)

Calculate the idealized, continuous GWB amplitude.

Calculation follows [Phinney2001] (Eq.5) or equivalently [Enoki+Nagashima-2007] (Eq.3.6). This calculation assumes a smooth, continuous population of binaries that are purely GW driven. * There are no finite-number effects. * There are no environmental or non-GW driven evolution effects. * There is no coalescence of binaries cutting them off at high-frequencies.

gwb(fobs_gw_edges, hard=<holodeck.hardening.Hard_GW object>, realize=100, loudest=1, params=False)

Calculate the (smooth/semi-analytic) GWB and CWs at the given observed GW-frequencies.

Parameters:

• fobs_gw_edges ((F,) array_like of scalar,) – Observer-frame GW-frequencies. [1/sec] These are the frequency bin edges, which are integrated across to get the number of binaries in each frequency bin.

• hard (holodeck.evolution._Hardening class or instance) – Hardening mechanism to apply over the range of fobs_gw.

• realize (int) – Specification of how many discrete realizations to construct. Realizations approximate the finite-source effects of a realistic population.

• loudest (int) – Number of loudest single sources to distinguish from the background.

• params (Boolean) – Whether or not to return astrophysical parameters of the binaries.

Returns:

• hc_ss ((F, R, L) NDarray of scalars) – The characteristic strain of the L loudest single sources at each frequency.

• hc_bg ((F, R) NDarray of scalars) – Characteristic strain of the GWB.

• sspar ((4, F, R, L) NDarray of scalars) – Astrophysical parametes (total mass, mass ratio, initial redshift, final redshift) of each loud single sources, for each frequency and realization. Returned only if params = True.

• bgpar ((7, F, R) NDarray of scalars) – Average effective binary astrophysical parameters (total mass, mass ratio, initial redshift, final redshift, final comoving distance, final separation, final angular separation) for background sources at each frequency and realization, Returned only if params = True.

rate_chirps(hard=None, integrate=True)

Find the event rate of binary coalescences (‘chirps’).

Get the number of coalescence events per unit time, in units of [1/sec].

Parameters:

• hard (None or _Hardening subclass instance) –

• integrate (bool) –

Returns:

• redz_final ((M, Q, Z)) – Redshift of binary coalescence. Binaries stalling before z=0, have values set to -1.0.

• rate (ndarray) – Rate of coalescence events in each bin, in units of [1/sec]. The shape and meaning depends on the value of the integrate flag:

  • if integrate == True, then the returned values is dN/dt, with shape (M-1, Q-1, Z-1)

  • if integrate == False, then the returned values is dN/[dlog10M dq dz dt], with shape (M, Q, Z)

_integrate_event_rate(rate)

_ndens_gal(mass_gal, mrat_gal, redz)

_ndens_mbh(mass_gal, mrat_gal, redz)

_integrated_binary_density(ndens=None, sum=True)

dynamic_binary_number(*args, **kwargs)

new_gwb(*args, **kwargs)

holodeck.sams.sam.sample_sam_with_hardening(sam, hard, fobs_orb=None, sepa=None, sample_threshold=10.0, cut_below_mass=None, limit_merger_time=None, **sample_kwargs)

Discretize Semi-Analytic Model into sampled binaries assuming the given binary hardening rate.

Parameters:

• sam (Semi_Analytic_Model) – Instance of an initialized semi-analytic model.

• hard (holodeck.evolution._Hardening) – Binary hardening model for calculating binary hardening rates (dadt or dfdt).

• fobs_orb (ArrayLike) – Observer-frame orbital-frequencies. Units of [1/sec]. NOTE: Either fobs_orb or sepa must be provided, and not both.

• sepa (ArrayLike) – Binary orbital separation. Units of [cm]. NOTE: Either fobs_orb or sepa must be provided, and not both.

Returns:

• vals ((4, S) ndarray of scalar) – Parameters of sampled binaries. Four parameters are: * mtot : total mass of binary (m1+m2) in [grams] * mrat : mass ratio of binary (m2/m1 <= 1) * redz : redshift of binary * fobs_orb / sepa : observer-frame orbital-frequency [1/s] or binary separation [cm]

• weights ((S,) ndarray of scalar) – Weights of each sample point.

• edges ((4,) of list of scalars) – Edges of parameter-space grid for each of above parameters (mtot, mrat, redz, fobs_orb) The lengths of each list will be [(M,), (Q,), (Z,), (F,)]

• dnum ((M, Q, Z, F) ndarray of scalar) – Number-density of binaries over grid specified by edges.

holodeck.sams.sam.evolve_eccen_uniform_single(sam, eccen_init, sepa_init, nsteps)

Evolve binary eccentricity from an initial value along a range of separations.

Parameters:

• sam (holodeck.sam.Semi_Analytic_Model instance) – The input semi-analytic model. All this does is provide the range of total-masses to determine the minimum ISCO radius, which then determines the smallest separations to evolve until.

• eccen_init (float,) – Initial eccentricity of binaries at the given initial separation sepa_init. Must be between [0.0, 1.0).

• sepa_init (float,) – Initial binary separation at which evolution begins. Units of [cm].

• nsteps (int,) – Number of (log-spaced) steps in separation between the initial separation sepa_init, and the final separation which is determined as the minimum ISCO radius based on the smallest total-mass of binaries in the sam instance.

Returns:

• sepa ((E,) ndarray of float) – The separations at which the eccentricity evolution is defined over. This is the independent variable of the evolution. The shape E is the value of the nsteps parameter.

• eccen ((E,)) – The eccentricity of the binaries at each location in separation given by sepa. The shape E is the value of the nsteps parameter.

holodeck.sams.sam.add_scatter_to_masses(mtot, mrat, dens, scatter, refine=4, log=None)

Add the given scatter to masses m1 and m2, for the given distribution of binaries.

The procedure is as follows (see dev-notebooks/sam-ndens-scatter.ipynb):

• 1. The density is first interpolated to a uniform, regular grid in (m1, m2) space. A 2nd order interpolant is used first. A 0th-order interpolant is used to fill-in bad values. In-between, a 1st-order interpolant is used if linear_interp_backup is True.

• 2. The density distribution is convolved with a smoothing function along each axis (m1, m2) to account for scatter.

• 3. The new density distribution is interpolated back to the original (mtot, mrat) grid.

Parameters:

• mtot ((M,) ndarray) – Total masses in grams.

• mrat ((Q,) ndarray) – Mass ratios.

• dens ((M, Q) ndarray) – Density of binaries over the given mtot and mrat domain.

• scatter (float) – Amount of scatter in the M-MBulge relationship, in dex (i.e. over log10 of masses).

• refine (int,) – The increased density of grid-points used in the intermediate (m1, m2) domain, in step (1).

• linear_interp_backup (bool,) – Whether a linear interpolant is used to fill-in bad values after the 2nd order interpolant. This generally doesn’t seem to fix any values.

• logspace_interp (bool,) – Whether interpolation should be performed in the log-space of masses. NOTE: strongly recommended.

Returns:

m1m2_dens – Binary density with scatter introduced.

Return type:

(M, Q) ndarray,