Stellar-Remnant Population Synthesis

Scrollable formula. Use Left and Right arrow keys to reveal hidden content.

Scrollable table. Use Left and Right arrow keys to reveal additional columns.

Scrollable figure. Use Left and Right arrow keys to inspect hidden content.

Population synthesis is a forward-modeling discipline: draw stellar and orbital initial conditions, evolve them through uncertain physical branches, and retain the weighted remnants and binaries that satisfy a declared selection. Its central scientific product is not one histogram but a conditional distribution together with the assumptions that generated it. Agreement with an observed mass distribution can be informative, yet it is never unique when winds, mass transfer, common envelopes, collapse, kicks, dynamics, cosmic history, and detector selection can compensate for one another.

Helpful background. Black-Hole Geometry and Mechanics supplies the mass and spin quantities assigned to each remnant; Relativistic Matter and Compact Objects supplies stellar and binary ingredients; Black-Hole Astrophysics supplies survey selection and population interpretation; and Black-Hole Dynamics and Collapse supplies the merger, encounter, and collapse stages represented by transition rules in a fast model.

A forward map with declared conditioning

Let $\boldsymbol\theta_{\rm ZAMS}$ collect initial masses, mass ratio, orbital period, eccentricity, metallicity, spins, and environment. Let $\boldsymbol\Lambda$ denote hyperparameters such as an initial-mass-function slope, wind normalization, common-envelope efficiency, and kick model. A compact representation is

p(\boldsymbol\theta_{\rm rem}\mid\boldsymbol\Lambda) =\int d\boldsymbol\theta_{\rm ZAMS}\, p(\boldsymbol\theta_{\rm ZAMS}\mid\boldsymbol\Lambda) K(\boldsymbol\theta_{\rm rem}\mid \boldsymbol\theta_{\rm ZAMS},\boldsymbol\Lambda),

where the kernel $K$ includes deterministic evolution, stochastic draws, branch probabilities, and delay times. The normalization needs care. If $p(\boldsymbol\theta_{\rm rem}\mid\boldsymbol\Lambda)$ is a probability distribution over all outcomes, the outcome space must include systems that do not make the desired remnant or merger. If those failures are omitted, the same integral instead gives an unnormalized successful-outcome density whose integral is the formation probability. Normalizing it to one removes that probability and retains only the conditional shape. The output cannot be interpreted without both the input measure and the evolution kernel.

For a cosmic merger rate, another layer weights metallicity-dependent formation efficiency by the star-formation and metallicity history and maps formation redshift to merger redshift through a delay-time distribution. A detector-facing distribution then multiplies by a selection probability before normalization. Intrinsic remnant production, mergers per formed stellar mass, and detected events per observing time are therefore different denominators.

For example, let $p_{\rm int}(\theta)$ be a normalized intrinsic distribution and $P_{\rm det}(\theta)$ the probability that a source with parameters $\theta$ enters the sample. Then

p_{\rm det}(\theta) =\frac{P_{\rm det}(\theta)p_{\rm int}(\theta)}{\alpha}, \qquad \alpha=\int P_{\rm det}(\theta)p_{\rm int}(\theta)\,d\theta.

The factor $\alpha$ is the detected fraction in this specified source ensemble. Equal intrinsic numbers in two classes with detection probabilities $0.8$ and $0.2$ produce an expected detected ratio of $4:1$ . Recovering a birth population therefore requires more than interpreting the peaks of the observed histogram. For a cosmological event rate, source-time rates must additionally be converted to observer time using $dt_{\rm obs}=(1+z)dt_{\rm src}$ .

Physical branches and adversarial checks

A single-star branch maps winds, core structure, pair instability, explosion, and fallback to a remnant. An isolated-binary branch additionally treats stable and unstable mass transfer, tides, common-envelope evolution, supernova mass loss, natal kicks, and gravitational-radiation delay. Dense-environment models add encounters, exchanges, ejections, stellar collisions, and hierarchical mergers. Belczynski and collaborators give a detailed example of these coupled prescriptions and their internal comparisons in Belczynski et al. 2008, §§ 2–7 and § 8.

Every branch needs local checks. Conservative mass transfer should conserve the modeled mass and orbital angular momentum; symmetric instantaneous collapse should recover the half-mass disruption criterion; zero wind should not remove angular momentum; and a refinement or sample-size study should stabilize weighted observables. Global checks compare more than the target distribution: binary fractions, compact-object rates, metallicity trends, delay times, and known limiting cases can expose compensating errors.

Rapid codes trade detailed stellar structure for coverage. That trade is legitimate when interpolation domains, branch boundaries, random seeds, sampling weights, code version, and input files are published. Dominik and collaborators vary common-envelope, wind, kick, and explosion assumptions to show how model families change compact-binary populations Dominik et al. 2012, §§ 2–3. Such ensembles are more informative than a single “fiducial” curve with no alternatives.

The following calculation is a teaching reconstruction, not a claimed rerun of StarTrack, COMPAS, or another research code. It asks for the fraction of progenitors above a toy black-hole threshold $20M_\odot$ within $8\leq M/M_\odot\leq100$ , under $p(M)\propto M^{-2.35}$ . The composite trapezoidal algorithm is compared with the independent exact solution of the power-law integral. The displayed output was generated with Python 3.9.6 and NumPy==2.0.2. Save the script as population_fraction_check.py.

import numpy as np

ALPHA = 2.35
LOW, CUT, HIGH = 8.0, 20.0, 100.0

def power_integral(a, b):
    exponent = 1.0 - ALPHA
    return (b**exponent - a**exponent) / exponent

exact = power_integral(CUT, HIGH) / power_integral(LOW, HIGH)
print("exact_fraction", f"{exact:.10f}")
print("N  numerical_fraction  absolute_error")
for n in (256, 512, 1024):
    all_mass = np.linspace(LOW, HIGH, n + 1)
    bh_mass = np.linspace(CUT, HIGH, n + 1)
    denominator = np.trapezoid(all_mass**(-ALPHA), all_mass)
    numerator = np.trapezoid(bh_mass**(-ALPHA), bh_mass)
    estimate = numerator / denominator
    print(n, f"{estimate:.10f}", f"{abs(estimate-exact):.3e}")

The program must produce the exact benchmark and three rows containing interval count, numerical fraction, and absolute error. The exact save-and-run recipe is:

python3 -m venv .venv-pop
. .venv-pop/bin/activate
python3 -m pip install numpy==2.0.2
python3 population_fraction_check.py

Grid intervals $N$	Numerical fraction	Absolute error
256	0.2658697542	0.0001273
512	0.2659652043	0.00003185
1024	0.2659890916	0.000007965

The exact value is $0.2659970562$ , and doubling $N$ reduces the error by about four, consistent with second-order trapezoidal convergence on each smooth interval. This verifies integration and normalization only. The principal numerical failure mode in a real calculation is missing a narrow or discontinuous evolutionary branch; important physical sensitivity enters through the invented $20M_\odot$ threshold itself. A converged toy integral is not a validated remnant model.

Results and limitations through 2026-08-18

Forward modeling and the need to propagate selection are durable. The efficiencies of common envelopes, mass transfer, natal kicks, pair-instability edges, cluster processing, and cosmic metallicity evolution remain unresolved. Competing formulations can fit overlapping observables in different parameter regimes and yield different source populations.

Position or evidence	Result	Scope and limitation
Detailed binary-evolution suites	Demonstrate viable isolated routes and quantify conditional efficiencies	Sensitive to common-envelope survival, transfer stability, winds, and kicks
Rapid Monte Carlo or weighted grids	Cover broad initial conditions and cosmic histories efficiently	Interpolation and discrete branches can hide structural-model error
Dense-cluster calculations	Produce exchanges, ejections, collisions, and hierarchical mergers	Depend on cluster initial conditions, stellar evolution, and relativistic encounter treatment
Detected compact-binary populations	Constrain mixtures of masses, spins, redshifts, and rates	Conclusions depend on waveform model, inference prior, selection function, and channel mixture

Stevenson and collaborators illustrate how a rapid binary code can connect a modeled intrinsic population to gravitational-wave observations while exposing its assumptions Stevenson et al. 2017, Results, Fig. 1, and Methods. This treatment should be reviewed after public code comparisons isolate a major branch discrepancy, after selection-corrected observations decisively distinguish field and dynamical mixtures, or after revised stellar, cluster, or cosmic-metallicity calculations materially change the predicted mass, spin, redshift, or delay-time distributions.

Interpreting an output population

A credible result names the code and commit, input distributions and weights, random seed or deterministic grid, interpolation rules, physical switches, cosmology, metallicity history, normalization, and selection function. It reports convergence in both sampling and grid resolution and supplies machine-readable output. A posterior constraint on $\boldsymbol\Lambda$ also requires priors and a likelihood; forward histograms alone are not inference.

This page does not adjudicate common-envelope physics or perform a gravitational-wave hierarchical analysis. It defines the calculation that must precede those claims and the checks that keep conditional predictions honest.

Common pitfalls

Normalizing away the rate. Two models can have similar normalized mass shapes but efficiencies that differ by orders of magnitude. Shape and absolute rate should be tested separately.

Counting simulated systems equally after biased sampling. If the simulation draws initial conditions from $q(\theta)$ but the target population has density $p(\theta)$ , the weight is proportional to $p(\theta)/q(\theta)$ , not merely $1/q(\theta)$ . The sampling distribution must be nonzero wherever the target contributes.

Exercise

A cosmic mixture forms $70\%$ of its stellar mass in a high-metallicity component with merger yield $10^{-3}$ per solar mass of formed stars and $30\%$ in a low-metallicity component with yield $10^{-2}$ per solar mass. Treat these as toy yields integrated over all delay times. Find the total yield and the fraction of mergers contributed by the low-metallicity component.

Solution

The mass-weighted yield is

\eta=\left[0.70(10^{-3})+0.30(10^{-2})\right]M_\odot^{-1} =3.7\times10^{-3}M_\odot^{-1}.

The low-metallicity contribution is $3.0\times10^{-3}M_\odot^{-1}$ , so its share of mergers is $3.0/3.7=0.811$ , or about $81\%$ . A minority of star formation dominates because its yield is ten times larger. A merger rate at one epoch also requires the formation history convolved with delay times; detection selection can change the observed share again.

References

Belczynski, K., Kalogera, V., Rasio, F. A., Taam, R. E., Zezas, A., Bulik, T., Maccarone, T. J., and Ivanova, N. (2008). “Compact Object Modeling with the StarTrack Population Synthesis Code.” The Astrophysical Journal Supplement Series 174, 223–260. DOI.
Dominik, M., Belczynski, K., Fryer, C., Holz, D. E., Berti, E., Bulik, T., Mandel, I., and O’Shaughnessy, R. (2012). “Double Compact Objects. I. The Significance of the Common Envelope on Merger Rates.” The Astrophysical Journal 759, 52. DOI.
Stevenson, S., Vigna-Gómez, A., Mandel, I., Barrett, J. W., Neijssel, C. J., Perkins, D., and de Mink, S. E. (2017). “Formation of the First Three Gravitational-wave Observations through Isolated Binary Evolution.” Nature Communications 8, 14906. DOI.