holodeck.discrete.evolution

Module for binary evolution from the time of formation/galaxy-merger until BH coalescence.

#!NOTE: much of this documentation needs to be updated to reflect that much of the material in this #! file was moved to holodeck.hardening.

In holodeck, initial binary populations are typically defined near the time of galaxy-galaxy merger, when two MBHs come together at roughly kiloparsec scales. Environmental ‘hardening’ mechanisms are required to dissipate orbital energy and angular momentum, allowing the binary separation to shrink (‘harden’). Typically dynamical friction (DF) is most important at large scales (\(\sim \mathrm{kpc}\)). Near where the pair of MBHs become a gravitationally-bound binary, the DF approximations break down, and individual stellar scattering events (from stars in the ‘loss cone’ of parameter space) need to be considered. In the presence of significant gas (i.e. from accretion), a circumbinary accretion disk (CBD) may form, and gravitational circumbinary disk torques from the gas distribution (typically spiral waves) may become important. Finally, at the smaller separations, GW emission takes over. The classic work describing the overall process of binary evolution in its different stages is [BBR1980]. [Kelley2017a] goes into significant detail on implementations and properties of each component of binary evolution also. Triple MBH interactions, and perturbations from other massive objects (e.g. molecular clouds) can also be important.

The holodeck.evolution submodule provides tools for modeling the binary evolution from the time of binary ‘formation’ (i.e. galaxy merger) until coalescence. Models for binary evolutionary can range tremendously in complexity. In the simplest models, binaries are assumed to coalesce instantaneously (in that the age of the universe is the same at formation and coalescence), and are also assumed to evolve purely due to GW emission (in that the time spent in any range of orbital frequencies can be calculated from the GW hardening timescale). Note that these two assumptions are contradictory, but commonly employed in the literature. The ideal, self-consistent approach, is to model binary evolution using self-consistently derived environments (e.g. gas and stellar mass distributions), and applying the same time-evolution prescription to both the redshift evolution of each binary, and also the GW calculation. Note that GWs can only be calculated based on some sort of model for binary evolution. The model may be extremely simple, in which case it is sometimes glanced over.

The core component of the evolution module is the Evolution class. This class combines a population of initial binary parameters (i.e. from the holodeck.population.Pop_Illustris class), along with a specific binary hardening model (_Hardening subclass), and performs the numerical integration of each binary over time - from large separations to coalescence. The Evolution class also acts to store the binary evolution histories (‘trajectories’ or ‘tracks’), which are then used to calculate GW signals (e.g. the GWB). To facilitate GW and similar calculations, Evolution also provides interpolation functionality along binary trajectories.

To-Do

  • General

    • evolution modifiers should act at each step, instead of after all steps? This would be a way to implement a changing accretion rate, for example; or set a max/min hardening rate.

    • re-implement “magic” hardening models that coalesce in zero change-of-redshift or fixed amounts of time.

  • Evolution

    • _sample_universe() : sample in comoving-volume instead of redshift

References

class holodeck.discrete.evolution.Evolution(pop, hard, nsteps: int = 100, mods=None, debug: bool = False, acc=None)[source]

Bases: object

Base class to evolve discrete binary populations forward in time.

NOTE: This class is primary built to be used with holodeck.population.Pop_Illustris.

The Evolution class is instantiated with a holodeck.population._Population_Discrete instance, and a particular binary hardening model (subclass of _Hardening). It then numerically integrates each binary from their initial separation to coalescence, determining how much time that process takes, and thus the rate of redshift/time evolution. NOTE: at initialization, a regular range of binary separations are chosen for each binary being evolved, and the integration calculates the time it takes for each step to complete. This is unlike most types of dynamical integration in which there is a prescribed time-step, and the amount of distance (etc) traveled over that time is then calculated.

Initialization: all attributes are set to empty arrays of the appropriate size. NOTE: the 0th step is not initialized at this time, it happens in Evolution.evolve().

Evolution: binary evolution is performed by running the Evolution.evolve() function. This function first calls Evolution._init_step_zero(), which sets the 0th step values, and then iterates over each subsequent step, calling Evolution._take_next_step(). Once all steps are taken (integration is completed), then Evolution._finalize() is called, at which points any stored modifiers (utils._Modifier subclasses, in the Evolution._mods attribute) are applied.

NOTE: whenever frequencies are used (rest-frame or observer-frame), they refer to orbital frequencies, not GW frequencies. For circular binaries, GW-freq = 2 * orb-freq.

Additional Notes

acc: Instance of accretion class. This supplies the method by which total accretion

rates are divided into individual accretion rates for each BH. By default, accretion rates are calculated at every step as a fraction of the Eddington limit. If acc contains a path to an accretion rate file which already stores total accretion rates at every timestep, then we omit the step where we calculate mdot_total as a fraction of the Eddington limit. This gives the flexibility to include accretion rates motivated by e.g. Illustris or other cosmological simulations.

_EVO_PARS = ['mass', 'sepa', 'eccen', 'scafa', 'tlook', 'dadt', 'dedt']
_LIN_INTERP_PARS = ['eccen', 'scafa', 'tlook', 'dadt', 'dedt']
_SELF_CONSISTENT = None
_STORE_FROM_POP = ['_sample_volume']
_pop

initial binary population instance

_debug

debug flag for performing extra diagnostics and output

_nsteps

number of integration steps for each binary

_mods

modifiers to be applied after evolution is completed

scafa

scale-factor of the universe, set to 1.0 after z=0

tlook

negative after redshift zero

Type:

lookback time [sec], NOTE

sepa

semi-major axis (separation) [cm]

mass

mass of BHs [g], 0-primary, 1-secondary

mdot

accretion rate onto each component of binary [g/s]

dadt

hardening rate in separation [cm/s]

eccen

eccentricity [], None if not being evolved

dedt

eccen evolution rate [1/s], None if not evolved

evolve(progress=False)[source]

Evolve binary population from initial separation until coalescence in discrete steps.

Each binary has a fixed number of ‘steps’ from initial separation until coalescence. The role of the evolve() method is to determine the amount of time each step takes, based on the ‘hardening rate’ (in separation and possible eccentricity i.e. \(da/dt\) and \(de/dt\)). The hardening rate is calculated from the stored _Hardening instances in the iterable Evolution._hard attribute.

When Evolution.evolve() is called, the 0th step is initialized separately, using the Evolution._init_step_zero() method, and then each step is integrated by calling the Evolution._take_next_step() method. Once all steps are completed, the Evolution._finalize() method is called, where any stored modifiers are applied.

Parameters:

progress (bool,) – Show progress-bar using tqdm package.

modify(mods=None)[source]

Apply and modifiers after evolution has been completed.

at(xpar, targets, params=None, coal=False, lin_interp=None)[source]

Interpolate evolution to the given target locations in either separation or frequency.

The variable of interpolation is specified with the xpar argument. The evolutionary tracks are interpolated in that variable, to the new target locations given by targets. We use ‘x’ to refer to the independent variable, and ‘y’ to refer to the dependent variable that is being interpolated. Which values are interpolated are specified with the params argument.

The behavior of this function is broken into three sub-functions, that are only used here:

Parameters:

  • xpar (str, in ['fobs', 'sepa']) – String specifying the variable of interpolation.

  • targets (array_like,) –

    Locations to interpolate to:

    • if xpar == sepa : binary separation, units of [cm],

    • if xpar == fobs : binary orbital freq, observer-frame, units of [1/sec],

  • params (None or (list of str)) – Names of the parameters that should be interpolated. If None, defaults to Evolution._EVO_PARS attribute.

  • coal (bool,) – Only store evolution values for binaries coalescing before redshift zero. Interpolated values for other binaries are set to np.nan.

  • lin_interp (None or bool,) –

    Interpolate parameters in linear space.

    • True : all parameters interpolated in lin-lin space.

    • False: all parameters interpolated in log-log space.

    • None : parameters are interpolated in log-log space, unless they’re included in the Evolution._LIN_INTERP_PARS attribute.

Returns:

vals – Dictionary of arrays for each interpolated parameter. The returned shape is (N, T), where T is the number of target locations to interpolate to, and N is the total number of binaries. Each data array is filled with np.nan values if the targets are outside of its evolution track. If coal=True, then binaries that do not coalesce before redshift zero also have their data array values fillwed with np.nan.

Return type:

dict,

Notes

  • Out of bounds values are set to np.nan.

  • Interpolation is 1st-order in log-log space in general, but properties which are in the _LIN_INTERP_PARS array are interpolated at 1st-order in lin-lin space. Parameters which can be negative should be interpolated in linear space. Passing a boolean for the lin_interp parameter will override the behavior (see Parameters above).

_at__inputs(xpar, targets, params, lin_interp)[source]

Parse/sanitize the inputs of the Evolution.at() method.

Parameters:

  • xpar (str, in ['fobs', 'sepa']) – String specifying the variable of interpolation.

  • targets (array_like,) – Locations to interpolate to. One of: * if xpar == sepa : binary separation, units of [cm], * if xpar == fobs : binary orbital-frequency, observer-frame, units of [1/sec].

  • params (None or list[str]) – Names of parameters that should be interpolated. If None, defaults to Evolution._EVO_PARS attribute.

Returns:

  • xnew (np.ndarray) – (T,) Log10 of the target locations to interpolate to, i.e. log10(targets).

  • xold (np.ndarray) – (N, M) Log10 of the x-values at which to evaluate the target interpolation points. Either log10(sepa) or log10(freq_orb_obs). NOTE: these values will be returned in increasing order. If they have been reversed, then rev=True.

  • params (list[str]) – Names of parameters that should be interpolated.

  • rev (bool) – Whether or not the xold array has been reversed.

  • squeeze (bool) – Whether or not the targets were a single scalar value (i.e. T=1, as opposed to an iterable). If targets were a scalar, then the data returned from Evolution.at() will be shaped as (N,) instead of (N,T); since in this case, T=1.

_at__index_frac(xnew, xold)[source]

Find indices bounding target locations, and the fractional distance to go between them.

Parameters:

  • xnew (np.ndarray) – Target locations to interplate to. Shape (T,).

  • xold (np.ndarray) – Values of the x-coordinate between which to interpolate. Shape (N, M). These are the x-values of either sepa or fobs from the evolutionary tracks of each binary.

Returns:

  • cut_idx (np.ndarray) – For each binary, the step-number indices between which to interpolate, for each target interpolation point. shape (2, N, T); where the 0th dimension, the 0th value is the low/before index, and the 1th value is the high/after index. i.e. for binary ‘i’ and target ‘j’, we need to interpolate between indices given by [0, i, j] and [1, i, j].

  • interp_frac (np.ndarray) – The fractional distance between the low value and the high value, to interpolate to. Shape (2, N, M). For binary ‘i’ and target ‘j’, interp_frac[i, j] is how the fraction of the way, from index cut_idx[0, i, j] to cut_idx[1, i, j] to interpolate to, in the data array.

  • valid (np.ndarray) – Array of boolean values, giving whether or not the target interpolation points are within the bounds of each binary’s evolution track. Shape (N, T).

_at__interpolate_array(yold, cut_idx, interp_frac, lin_interp_flag)[source]

Interpolate a parameter to a fraction between integration steps.

Parameters:

  • yold (np.ndarray) – The data to be interpolated. This is the raw evolution data, for each binary and each step. Shaped either as (N, M) or (N, M, 2) if parameter is mass.

  • cut_idx (np.ndarray) – For each binary, the step-number indices between which to interpolate, for each target interpolation point. shape (2, N, T); where the 0th dimension, the 0th value is the low/before index, and the 1th value is the high/after index. i.e. for binary ‘i’ and target ‘j’, we need to interpolate between indices given by [0, i, j] and [1, i, j].

  • interp_frac (np.ndarray) – The fractional distance between the low value and the high value, to interpolate to. Shape (2, N, M). For binary ‘i’ and target ‘j’, interp_frac[i, j] is how the fraction of the way, from index cut_idx[0, i, j] to cut_idx[1, i, j] to interpolate to, in the data array.

  • lin_interp_flag (bool,) – Whether data should be interpolated in lin-lin space (True), or log-log space (False).

Returns:

ynew – The input data interpolated to the new target locations. Shape is (N, T) or (N, T, 2) for N-binaries, T-target points. A third dimension is present if the input data was 3D.

Return type:

np.ndarray

sample_universe(fobs_orb_edges, down_sample=None)[source]

Construct a full universe of binaries based on resampling this population.

Parameters:

  • fobs (array_like,) – Observer-frame orbital-frequencies at which to sample population. Units of [1/sec].

  • down_sample (None or float,) – Factor by which to downsample the resulting population. For example, 10.0 will produce 10x fewer output binaries.

Returns:

  • names ((4,) list of str,) – Names of the returned binary parameters (i.e. each array in samples and vals).

  • samples ((4, S) ndarray,) – Sampled binary data. For each binary samples S, 4 parameters are returned: [‘mtot’, ‘mrat’, ‘redz’, ‘fobs’] (these are listed in the names returned value.) NOTE: fobs is observer-frame orbital-frequencies. These values correspond to all of the binaries in an observer’s Universe (i.e. light-cone), within the given frequency bins. The number of samples S is roughly the sum of the weights — but the specific number is drawn from a Poisson distribution around the sum of the weights.

  • vals ((4,) list of (V,) ndarrays or float) – Binary parameters (log10 of parameters specified in the names return values) at each frequency bin. Binaries not reaching the target frequency bins before redshift zero, or before coalescing, are not returned. Thus the number of values V may be less than F*N for F frequency bins and N binaries.

  • weights ((V,) ndarray of float) – The weight of each binary-frequency sample. i.e. number of observer-universe binaries corresponding to this binary in the simulation, at the target frequency.

  • To-Do

  • —–

  • * This should sample in volume instead of redz, see how it’s done in sam module.

_sample_universe__at_values_weights(fobs_orb_edges)[source]

Interpolate binary histories to target frequency bins, obtaining parameters and weights.

The weights correspond to the number of binaries in an observer’s Universe (light-cone) corresponding to each simulated binary sample.

Parameters:

fobs_orb_edges ((F+1,) arraylike) – Edges of target frequency bins to sample population. These are observer-frame orbital frequencies. Binaries are interpolated to frequency bin centers, calculated from the midpoints of the provided bin edges.

Returns:

  • names ((4,) list of str,) – Names of the returned binary parameters (i.e. each array in vals).

  • vals ((4,) list of (V,) ndarrays or float) – Binary parameters (log10 of parameters specified in the names return values) at each frequency bin. Binaries not reaching the target frequency bins before redshift zero, or before coalescing, are not returned. Thus the number of values V may be less than F*N for F frequency bins and N binaries.

  • weights ((V,) ndarray of float) – The weight of each binary-frequency sample. i.e. number of observer-universe binaries corresponding to this binary in the simulation, at the target frequency.

_sample_universe__resample(fobs_orb_edges, vals, weights, down_sample)[source]
_init_step_zero()[source]

Set the initial conditions of the binaries at the 0th step.

Transfers attributes from the stored holodeck.population._Population_Discrete instance to the 0th index of the evolution arrays. The attributes are [sepa, scafa, mass, and optionally eccen]. The hardening model is also used to calculate the 0th hardening rates dadt and dedt. The initial lookback time, tlook is also set.

_take_next_step(step)[source]

Integrate the binary population forward (to smaller separations) by one step.

For an integration step s, we are moving from index s-1 to index s. These correspond to the ‘left’ and ‘right’ edges of the step. The evolutionary trajectory values have already been calculated on the left edges (during either the previous time step, or the initial time step). Each subsequent integration step then proceeds as follows:

  1. The hardening rate is calculated at the right edge of the step.

  2. The time it takes to move from the left to right edge is calculated using a trapezoid rule in log-log space.

  3. The right edge evolution values are stored and updated.

Parameters:

step (int) – The destination integration step number, i.e. step=1 means integrate from 0 to 1.

_hardening_rate(step, store_debug=True)[source]

Calculate the net hardening rate for the given integration step.

The hardening rates (_Hardening subclasses) stored in the Evolution._hard attribute are called in sequence, their _Hardening.dadt_dedt() methods are called, and the \(da/dt\) and \(de/dt\) hardening rates are added together. NOTE: the da/dt and de/dt values are added together to get the net rate, this is an approximation.

Parameters:

step (int) – Current step number (the destination of the current step, i.e. step=1 is for integrating from 0 to 1.)

Returns:

  • dadt (np.ndarray) – The hardening rate in separation, \(da/dt\), in units of [cm/s]. The shape is (N,) where N is the number of binaries.

  • dedt (np.ndarray or None) – If eccentricity is not being evolved, this is None. If eccentricity is being evolved, this is the hardening rate in eccentricity, \(de/dt\), in units of [1/s]. In this case, the shape is (N,) where N is the number of binaries.

_check()[source]

Perform basic diagnostics on parameter validity after evolution.

_finalize()[source]

Perform any actions after completing all of the integration steps.

_update_derived()[source]

Update any derived quantities after modifiers are applied.

property shape

The number of binaries and number of steps (N, S).

property size

The number of binaries

property steps

The number of evolution steps

property coal

Indices of binaries that coalesce before redshift zero.

property tage

Age of the universe [sec] for each binary-step.

Derived from Evolution.tlook.

Returns:

ta – (B, S). Age of the universe.

Return type:

np.ndarray,

property mtmr

Total-mass and mass-ratio.

Returns:

  • mt (np.ndarray) – Total mass (\(M = m_1 + m_2\)) in [gram].

  • mr (np.ndarray) – Mass ratio (\(q = m_2/m_1 \leq 1.0\)).

property freq_orb_rest

Rest-frame orbital frequency. [1/s]

property freq_orb_obs

Observer-frame orbital frequency. [1/s]

_check_evolved()[source]

Raise an error if this instance has not yet been evolved.