"""Handle 1: dark companions in Gaia DR3 astrometric binaries within 300 pc. Steps A. Build the observed sample: NSS astrometric orbits with main-sequence primaries, d < 300 pc, G < 17.5, 100 < P < 1000 d. Compute the photocentre semi-major axis a0, the astrometric mass-ratio function A (AMRF), the class limits and the companion mass under the hypothesis that the companion emits no optical light. B. Identify 'optically dark' candidates (A above the single-main-sequence limit by >3 sigma), cross-match with Shahaf+2023 (class probabilities, Table 2 compact-object candidates) and Shahaf+2024 (red colour excess). C. Denominator: how many companions of each mass WOULD have been catalogued as dark candidates if every stellar companion were enclosed. The catalogue's own class-I (luminous) binaries calibrate the unmodelled part of the Gaia selection function in cells of (period, primary mass); the model only supplies the RELATIVE detectability of dark versus luminous companions and the companion mass-ratio distribution. D. 95% upper limits on the enclosed fraction of companions per companion-mass bin. E. Rough expectation for white-dwarf companions from progenitor counting. F. (revision) A white-dwarf-free sub-bin, 0.1-0.4 Msun, and a breakdown of the low-mass candidates: how many lie above the limit for a close main-sequence PAIR (class III), i.e. cannot be hierarchical triples, and how many are within 1 sigma of the white-dwarf peak. The model counts enclosed companions only for mass ratios 0.1 <= q = M2/M1 <= 1 (companions less massive than the visible star), whereas the observed dark candidates are counted for any q; this makes the limits conservative. """ import json import os import sys import numpy as np import pandas as pd from scipy import stats from scipy.special import ndtr import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) from common import RAW, DATA, FIGS from mslib import (MainSequence, amrf_dark, q_from_amrf_dark, amrf_luminous, amrf_limits, a0_from_thiele_innes) DMAX = 300.0 GMAX = 17.5 PMIN, PMAX = 100.0, 1000.0 M1_EDGES = np.array([0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.5]) P_EDGES = np.logspace(2, 3, 5) M2_BINS = [(0.1, 0.5), (0.5, 0.8), (0.8, 1.4), (0.1, 0.4)] # the 4th is a sub-bin of the 1st (revision) N_MAIN_BINS = 3 NSIG = 3.0 # significance by which A must exceed the single-MS limit DMS_CUT = -0.2 # 'single-light' cut: no more than 0.2 mag above the dwarf locus M2_GRID = np.arange(0.105, 1.45, 0.01) SIG_M1_REL = 0.10 # adopted relative uncertainty of photometric primary masses RED_THRESH = 0.64 # Shahaf+2024 (erratum) threshold on Pr(red) for 'red colour excess' SIGLOG = 0.20 # log-normal scatter of sigma_a0 and sigma_parallax at fixed G (measured below) KTHR = 1.0 # effective multiplier on the published a0/sigma_a0 threshold (fitted below) def poisson_ul(n, cl=0.95): """Classical one-sided Poisson upper limit on the mean for n observed events.""" return stats.chi2.ppf(cl, 2 * (n + 1)) / 2.0 def ms_offset(ms, bprp, MG): """Magnitude offset from the mean dwarf locus at the star's colour (negative = brighter).""" # locus: M_G as a function of BP-RP from the Mamajek table order = np.argsort(ms._bprp) return MG - np.interp(bprp, ms._bprp[order], ms._c_MG[order]) def load_observed(ms): d = pd.read_csv(os.path.join(RAW, "nss_astrometric.csv")) n_all = len(d) d = d.sort_values("significance", ascending=False).drop_duplicates("source_id") d["a0"] = a0_from_thiele_innes(d.a_thiele_innes.values, d.b_thiele_innes.values, d.f_thiele_innes.values, d.g_thiele_innes.values) d["sig_a0"] = d.a0 / d.significance d["dist"] = 1000.0 / d.plx_nss d["MG"] = d.phot_g_mean_mag + 5 * np.log10(d.plx_nss) - 10 d["dMS"] = ms_offset(ms, d.bp_rp.values, d.MG.values) d["M1"] = ms.mass_of_MG(d.MG.values) cuts = [("all astrometric NSS solutions (rows)", n_all), ("unique sources", len(d))] d = d[d.nss_solution_type.isin(["Orbital", "AstroSpectroSB1"])] cuts.append(("Orbital or AstroSpectroSB1", len(d))) d = d[(d.dist < DMAX)] cuts.append((f"distance 1/parallax < {DMAX:.0f} pc", len(d))) d = d[(d.phot_g_mean_mag < GMAX) & d.bp_rp.notna()] cuts.append((f"G < {GMAX} and BP-RP available", len(d))) d = d[(d.period > PMIN) & (d.period < PMAX)] cuts.append((f"{PMIN:.0f} < P < {PMAX:.0f} d", len(d))) d = d[(d.dMS > -1.3) & (d.dMS < 0.7)] cuts.append(("within -1.3/+0.7 mag of the dwarf locus (main-sequence primary)", len(d))) d = d[(d.M1 >= M1_EDGES[0]) & (d.M1 < M1_EDGES[-1])] cuts.append((f"photometric primary mass {M1_EDGES[0]}-{M1_EDGES[-1]} Msun", len(d))) d = d.copy() # AMRF, its uncertainty, class limits, dark-companion mass Pyr = d.period / 365.25 d["A"] = (d.a0 / d.plx_nss) * d.M1 ** (-1 / 3) * Pyr ** (-2 / 3) rel = np.sqrt((1 / d.significance) ** 2 + (d.plx_nss_err / d.plx_nss) ** 2 + (2 / 3 * d.period_error / d.period) ** 2 + (SIG_M1_REL / 3) ** 2) d["sig_A"] = d.A * rel grid = np.linspace(0.15, 1.6, 146) a_ms, a_tr = amrf_limits(ms, grid) d["A_ms"] = np.interp(d.M1, grid, a_ms) d["A_tr"] = np.interp(d.M1, grid, a_tr) d["M2_dark"] = q_from_amrf_dark(d.A.values) * d.M1 d["M2_dark_lo"] = q_from_amrf_dark(np.maximum(d.A - d.sig_A, 1e-3).values) * d.M1 * (1 - SIG_M1_REL * 2 / 3) d["M2_dark_hi"] = q_from_amrf_dark((d.A + d.sig_A).values) * d.M1 * (1 + SIG_M1_REL * 2 / 3) d["nsigI"] = (d.A - d.A_ms) / d.sig_A d["nsigIII"] = (d.A - d.A_tr) / d.sig_A d["nonI"] = d.nsigI > NSIG d["classIII"] = d.nsigIII > NSIG d["classI"] = d.A <= d.A_ms return d, cuts, (grid, a_ms, a_tr) def attach_external(d): bm = pd.read_csv(os.path.join(RAW, "binary_masses.csv")) bm = bm[bm.combination_method.isin(["Orbital+M1", "AstroSpectroSB1+M1"])].drop_duplicates("source_id") d = d.merge(bm[["source_id", "m1", "m1_lower", "m1_upper"]], on="source_id", how="left") s1 = pd.read_csv(os.path.join(RAW, "shahaf2023_table1.csv")) s1 = s1.rename(columns={"GaiaDR3": "source_id", "A": "A_shahaf", "M1": "M1_shahaf", "M2min": "M2_shahaf"}) d = d.merge(s1[["source_id", "A_shahaf", "M1_shahaf", "M2_shahaf", "PII", "PIII"]], on="source_id", how="left") s2t = pd.read_csv(os.path.join(RAW, "shahaf2023_table2.csv")) s2t = s2t.rename(columns={"GaiaDR3": "source_id", "Class": "shahaf23_class", "n_Class": "shahaf23_note", "M2min": "shahaf23_M2min"}) s2t["shahaf23_note"] = s2t.shahaf23_note.fillna("").astype(str).str.strip() d = d.merge(s2t[["source_id", "shahaf23_class", "shahaf23_note", "shahaf23_M2min"]].drop_duplicates("source_id"), on="source_id", how="left") from astropy.table import Table import warnings with warnings.catch_warnings(): warnings.simplefilter("ignore") t = Table.read(os.path.join(RAW, "shahaf2024_tableA1.fits")) s2 = pd.DataFrame({"source_id": np.array(t["source_id"]).astype(np.int64), "red_prob": np.array(t["red_prob"], float), "color_excess": np.array(t["color_excess"], float)}) d = d.merge(s2, on="source_id", how="left") return d def load_parent(ms): p = pd.read_csv(os.path.join(RAW, "parent_cells_300pc.csv")) d1, d2 = p.dbin * 25.0, (p.dbin + 1) * 25.0 p["dist"] = ((d1 ** 3 + d2 ** 3) / 2) ** (1 / 3) p["G"] = (p.gbin + 0.5) * 0.5 p["MG"] = (p.mgbin + 0.5) * 0.25 p["bprp"] = (p.cbin + 0.5) * 0.1 p["dMS"] = ms_offset(ms, p.bprp.values, p.MG.values) p["M1"] = ms.mass_of_MG(p.MG.values) n_all = p.n.sum() p = p[(p.dMS > -1.3) & (p.dMS < 0.7) & (p.M1 >= M1_EDGES[0]) & (p.M1 < M1_EDGES[-1]) & (p.dist < DMAX) & (p.G < GMAX)] # collapse colour: cells of (dist, G, MG) g = p.groupby(["dbin", "gbin", "mgbin"], as_index=False).agg(n=("n", "sum"), dist=("dist", "first"), G=("G", "first"), MG=("MG", "first"), M1=("M1", "first")) return g, int(n_all), int(p.n.sum()) def sigma_models(allnss): """Median sigma_a0 and sigma_parallax of published orbital solutions versus G, measured on solutions well above the acceptance threshold (so the threshold does not bias them).""" thr = np.maximum(5, 158 / np.sqrt(allnss.period)) hi = allnss[allnss.significance > 2 * thr] edges = np.arange(5.0, 18.01, 1.0) cen, sa, sp, sc = [], [], [], [] for a, b in zip(edges[:-1], edges[1:]): s = hi[(hi.phot_g_mean_mag >= a) & (hi.phot_g_mean_mag < b)] if len(s) < 30: continue cen.append(0.5 * (a + b)) sa.append(np.median(s.sig_a0)) sp.append(np.median(s.plx_nss_err)) sc.append(0.5 * (np.log10(s.sig_a0.quantile(0.84)) - np.log10(s.sig_a0.quantile(0.16)))) return np.array(cen), np.array(sa), np.array(sp), float(np.median(sc)) def qdist(name, q): """Un-normalised companion mass-ratio distributions p(q) on 0.1 <= q <= 1.""" if name == "flat": return np.ones_like(q) if name == "MDS17": # Moe & Di Stefano 2017-like for solar-type, P ~ 100-1000 d p = np.where(q < 0.3, (q / 0.3) ** 0.3, (q / 0.3) ** -0.5) return p * np.where(q > 0.95, 1.6, 1.0) # ~10% excess twins if name == "rising": return q raise ValueError(name) def detect_prob(a0, plx, P, G, smod): """Probability that an orbit with photocentre semi-major axis a0 [mas] passes the two published DR3 significance cuts, given the empirical sigma(G) and its scatter.""" cen, sa, sp, _ = smod lsa = np.interp(G, cen, np.log10(sa)) lsp = np.interp(G, cen, np.log10(sp)) smin = np.maximum(5.0, 158.0 / np.sqrt(P)) with np.errstate(divide="ignore"): p1 = ndtr((np.log10(np.maximum(a0, 1e-12) / (KTHR * smin)) - lsa) / SIGLOG) p2 = ndtr((np.log10(plx * P / 20000.0) - lsp) / SIGLOG) return p1 * p2 def forward_model(ms, parent, smod, lim, nsig=NSIG, nq=91, nP=3): """For every (P bin, M1 bin), summed over the parent stars: lum[qname] : detectability of luminous single-MS companions, averaged over p(q) dark[k] : detectability, as a non-class-I system, of a dark companion of mass M2_GRID[k] dark3[k] : same for class III encl[qname][k] : dark[k] weighted by the number of main-sequence companions per unit M2 implied by p(q) (q = M2/M1 <= 1), i.e. the yield if all were enclosed """ grid, a_ms, a_tr = lim q = np.linspace(0.1, 1.0, nq) M1 = parent.M1.values plx = 1000.0 / parent.dist.values G = parent.G.values w = parent.n.values.astype(float) Ams = np.interp(M1, grid, a_ms) Atr = np.interp(M1, grid, a_tr) S = np.array([ms.flux_ratio_G(q * m, m) for m in M1]) # (cells, nq) A_lum = amrf_luminous(q[None, :], S) cen, sa, sp, _ = smod sig_a0 = 10 ** np.interp(G, cen, np.log10(sa)) sig_px = 10 ** np.interp(G, cen, np.log10(sp)) qd = M2_GRID[None, :] / M1[:, None] # (cells, nm2) A_drk = amrf_dark(qd) pqs = {} for name in ("flat", "MDS17", "rising"): pq = qdist(name, q) pqs[name] = pq / np.trapezoid(pq, q) out = {} for ip in range(len(P_EDGES) - 1): Ps = np.logspace(np.log10(P_EDGES[ip]), np.log10(P_EDGES[ip + 1]), 2 * nP + 1)[1::2] for im in range(len(M1_EDGES) - 1): sel = (M1 >= M1_EDGES[im]) & (M1 < M1_EDGES[im + 1]) if not sel.any(): continue lum = {name: 0.0 for name in pqs} lum_cell = np.zeros(sel.sum()) dark = np.zeros(len(M2_GRID)) dark3 = np.zeros(len(M2_GRID)) encl = {name: np.zeros(len(M2_GRID)) for name in pqs} for P in Ps: # luminous companions a_au = ((M1[sel, None] * (1 + q[None, :])) * (P / 365.25) ** 2) ** (1 / 3) a0_l = plx[sel, None] * a_au / (1 + q[None, :]) ** (1 / 3) * A_lum[sel] pl = detect_prob(a0_l, plx[sel, None], P, G[sel, None], smod) for name, pq in pqs.items(): lum[name] += np.trapezoid((w[sel, None] * pq[None, :] * pl).sum(0), q) / len(Ps) lum_cell += w[sel] * np.trapezoid(pqs["flat"][None, :] * pl, q, axis=1) / len(Ps) # dark companions on the M2 grid (q may exceed 1) a_au = ((M1[sel, None] + M2_GRID[None, :]) * (P / 365.25) ** 2) ** (1 / 3) a0_d = plx[sel, None] * a_au * qd[sel] / (1 + qd[sel]) pd_ = detect_prob(a0_d, plx[sel, None], P, G[sel, None], smod) rel = np.sqrt((sig_a0[sel, None] / np.maximum(a0_d, 1e-9)) ** 2 + (sig_px[sel] / plx[sel])[:, None] ** 2 + (SIG_M1_REL / 3) ** 2) nonI = A_drk[sel] * (1 - nsig * rel) > Ams[sel, None] cIII = A_drk[sel] * (1 - nsig * rel) > Atr[sel, None] dark += (w[sel, None] * pd_ * nonI).sum(0) / len(Ps) dark3 += (w[sel, None] * pd_ * cIII).sum(0) / len(Ps) inq = (qd[sel] >= 0.1) & (qd[sel] <= 1.0) for name in pqs: dens = np.interp(qd[sel], q, pqs[name]) / M1[sel, None] * inq # companions per unit M2 encl[name] += (w[sel, None] * pd_ * nonI * dens).sum(0) / len(Ps) encl.setdefault(name + "_III", np.zeros(len(M2_GRID))) encl[name + "_III"] += (w[sel, None] * pd_ * cIII * dens).sum(0) / len(Ps) out[(ip, im)] = {"lum": lum, "dark": dark, "dark3": dark3, "encl": encl, "nstars": w[sel].sum(), "lum_cell": lum_cell, "sel": np.where(sel)[0]} return out def binfrac(m): """Binary-fraction scaling with primary mass (as assumed by El-Badry et al. 2024).""" m = np.asarray(m, float) lo = 0.40 * np.log10(np.maximum(m, 0.0801) / 0.08) hi = 0.40 + (0.57 - 0.40) * (np.minimum(m, 3.0) - 1.2) / (3.0 - 1.2) return np.where(m < 0.8, lo, np.where(m < 1.2, 0.40, hi)) def wd_expectation(m1, tau_gal=10.0): """White-dwarf companions per main-sequence companion (same period range today) for a visible star of mass m1, and their mass distribution on M2_GRID, from progenitor counting: Kroupa IMF, constant star-formation rate over tau_gal Gyr, main-sequence lifetime 10 Gyr (M/Msun)^-2.5, flat mass-ratio distribution, binary fraction as in binfrac(), the Cummings et al. (2018) initial-final mass relation, and NO change of the orbital period when the original primary becomes a white dwarf (the dominant uncertainty: at these separations the pair must have interacted).""" mi = np.linspace(max(m1, 0.85), 8.0, 6000) # original primary mass xi = mi ** -2.3 life = 10.0 * mi ** -2.5 frac_dead = np.clip(1 - life / tau_gal, 0, 1) qorig = m1 / mi okq = (qorig >= 0.1) & (qorig <= 1.0) wgt = xi * binfrac(mi) * frac_dead * okq / (0.9 * mi) # per unit m1 alive1 = 1 - float(np.clip(1 - (10.0 * m1 ** -2.5) / tau_gal, 0, 1)) den = m1 ** -2.3 * float(binfrac(m1)) * max(alive1, 0.05) mwd = ifmr(mi) dens = np.zeros(len(M2_GRID)) dm = M2_GRID[1] - M2_GRID[0] idx = np.clip(np.round((mwd - M2_GRID[0]) / dm).astype(int), 0, len(M2_GRID) - 1) np.add.at(dens, idx, wgt * (mi[1] - mi[0])) return float(np.trapezoid(wgt, mi) / den), dens / den / dm # ratio, density per unit M2 def ifmr(mi): """Initial-final mass relation (Cummings et al. 2018, PARSEC-based piecewise linear).""" mi = np.asarray(mi, float) return np.where(mi < 2.85, 0.080 * mi + 0.489, np.where(mi < 3.60, 0.187 * mi + 0.184, 0.107 * mi + 0.471)) def count_obs(obs, nsig, dms_cut, which="nonI"): col = "nsigI" if which == "nonI" else "nsigIII" sel = obs[col] > nsig if dms_cut is not None: sel &= obs.dMS > dms_cut out = [] for b in M2_BINS: out.append(int((sel & (obs.M2_dark >= b[0]) & (obs.M2_dark < b[1])).sum())) return out, sel def main(): global KTHR, SIGLOG ms = MainSequence() obs, cuts, lim = load_observed(ms) obs = attach_external(obs) grid, a_ms, a_tr = lim print("[census] cuts:") for c in cuts: print(" %-70s %7d" % c) pd.DataFrame(cuts, columns=["cut", "n_remaining"]).to_csv(os.path.join(DATA, "cuts_gaia_nss.csv"), index=False) # ---- cross-checks against published work ------------------------------------------------- chk = {} has = obs.A_shahaf.notna() Aown = (obs.a0 / obs.plx_nss) * obs.M1_shahaf ** (-1 / 3) * (obs.period / 365.25) ** (-2 / 3) chk["n_in_shahaf2023_table1"] = int(has.sum()) chk["median_ratio_my_A_to_shahaf_A_same_M1"] = float(np.nanmedian(Aown[has] / obs.A_shahaf[has])) chk["median_ratio_my_M1_to_gaia_m1"] = float(np.nanmedian(obs.M1 / obs.m1)) chk["scatter_my_M1_to_gaia_m1_1684"] = [float(x) for x in np.nanpercentile((obs.M1 / obs.m1).dropna(), [16, 84])] chk["sum_shahaf_PII_plus_PIII_in_sample"] = float((obs.PII.fillna(0) + obs.PIII.fillna(0)).sum()) chk["n_my_nonI_3sigma"] = int(obs.nonI.sum()) chk["n_my_classIII_3sigma"] = int(obs.classIII.sum()) both = has & obs.nonI chk["my_nonI_in_shahaf_sample"] = int(both.sum()) chk["my_nonI_mean_shahaf_P_nonI"] = float((obs.PII + obs.PIII)[both].mean()) s24 = obs.red_prob.notna() chk["n_in_shahaf2024_tableA1"] = int(s24.sum()) chk["shahaf2024_in_sample_recovered_as_nonI_3sigma"] = int((s24 & obs.nonI).sum()) chk["shahaf2024_no_red_excess_in_sample"] = int((obs.red_prob < RED_THRESH).sum()) chk["median_ecc_classI"] = float(obs[obs.classI].eccentricity.median()) chk["median_ecc_S24_no_red_excess"] = float(obs[obs.red_prob < RED_THRESH].eccentricity.median()) chk["median_ecc_S24_red_excess"] = float(obs[obs.red_prob >= RED_THRESH].eccentricity.median()) print("[census] checks:", json.dumps(chk, indent=1)) # ---- parent population -------------------------------------------------------------------- parent, n_parent_all, n_parent_ms = load_parent(ms) pr = pd.read_csv(os.path.join(RAW, "parent_random_1in40.csv")) pr["MG"] = pr.phot_g_mean_mag + 5 * np.log10(pr.parallax) - 10 pr["dMS"] = ms_offset(ms, pr.bp_rp.values, pr.MG.values) pr["M1"] = ms.mass_of_MG(pr.MG.values) pr = pr[(pr.dMS > -1.3) & (pr.dMS < 0.7) & (pr.M1 >= M1_EDGES[0]) & (pr.M1 < M1_EDGES[-1])] print(f"[census] parent stars within {DMAX:.0f} pc: {n_parent_all}; main-sequence 0.2-1.5 Msun: " f"{n_parent_ms}; cells {len(parent)}; random 1-in-40 subsample on the main sequence: {len(pr)}") def retention(cut): if cut is None: return np.ones(len(M1_EDGES) - 1) r = [] for im in range(len(M1_EDGES) - 1): s_ = pr[(pr.M1 >= M1_EDGES[im]) & (pr.M1 < M1_EDGES[im + 1])] r.append(float((s_.dMS > cut).mean())) return np.array(r) # ---- empirical sigma model ------------------------------------------------------------------ allnss = pd.read_csv(os.path.join(RAW, "nss_astrometric.csv"), usecols=["source_id", "period", "significance", "plx_nss", "plx_nss_err", "phot_g_mean_mag", "a_thiele_innes", "b_thiele_innes", "f_thiele_innes", "g_thiele_innes", "nss_solution_type"]) allnss = allnss[allnss.nss_solution_type.isin(["Orbital", "AstroSpectroSB1"])].copy() allnss["a0"] = a0_from_thiele_innes(allnss.a_thiele_innes.values, allnss.b_thiele_innes.values, allnss.f_thiele_innes.values, allnss.g_thiele_innes.values) allnss["sig_a0"] = allnss.a0 / allnss.significance smod = sigma_models(allnss) pd.DataFrame({"G": smod[0], "median_sigma_a0_mas": smod[1], "median_sigma_parallax_mas": smod[2]}).to_csv( os.path.join(DATA, "sigma_model.csv"), index=False, float_format="%.4g") print("[census] measured scatter of log sigma_a0 at fixed G (dex):", round(smod[3], 3), "adopted", SIGLOG) obsI = obs[obs.classI] dm2 = M2_GRID[1] - M2_GRID[0] # ---- fit the effective detection threshold --------------------------------------------------- # The two published cuts are necessary, not sufficient: the RUWE > 1.4 pre-selection and the # acceleration-solution cascade remove more low-significance orbits (El-Badry et al. 2024; # Lam et al. 2025). We absorb this into an effective multiplier KTHR on the a0/sigma_a0 # threshold and a width SIGLOG, fitted so that the model (calibrated per period and # primary-mass cell only) reproduces the DISTANCE and MAGNITUDE distributions of class-I orbits. dobs = np.histogram(obsI.dist, bins=np.arange(0, 301, 25))[0] gobs = np.histogram(obsI.phot_g_mean_mag, bins=np.arange(3, 19, 1))[0] def validate(k, wdt, nP=2): global KTHR, SIGLOG KTHR, SIGLOG = k, wdt fm_ = forward_model(ms, parent, smod, lim, nsig=NSIG, nP=nP) dpred, gpred = np.zeros(12), np.zeros(15) for (ip, im), v in fm_.items(): selo = ((obsI.period >= P_EDGES[ip]) & (obsI.period < P_EDGES[ip + 1]) & (obsI.M1 >= M1_EDGES[im]) & (obsI.M1 < M1_EDGES[im + 1])) contrib = selo.sum() / v["lum"]["flat"] * v["lum_cell"] cells = parent.iloc[v["sel"]] np.add.at(dpred, cells.dbin.values.astype(int), contrib) np.add.at(gpred, np.floor(cells.G.values).astype(int) - 3, contrib) chi = (np.sum((dobs - dpred) ** 2 / np.maximum(dpred, 1)) + np.sum((gobs - gpred) ** 2 / np.maximum(gpred, 1))) return chi, dpred, gpred fit_rows = [] for k in (1.0, 2.0, 2.5, 3.0, 3.5, 4.0, 5.0, 6.0): for wdt in (0.15, 0.2, 0.25, 0.3, 0.4): fit_rows.append({"KTHR": k, "SIGLOG": wdt, "chi2": validate(k, wdt)[0]}) fit = pd.DataFrame(fit_rows).sort_values("chi2") fit.to_csv(os.path.join(DATA, "threshold_fit.csv"), index=False, float_format="%.5g") K_BEST, W_BEST = float(fit.iloc[0].KTHR), float(fit.iloc[0].SIGLOG) chi_nominal = float(fit[(fit.KTHR == 1.0) & (fit.SIGLOG == 0.2)].chi2.iloc[0]) print(f"[census] effective threshold fit: KTHR={K_BEST}, SIGLOG={W_BEST}, chi2={fit.iloc[0].chi2:.0f} " f"(nominal published cuts only: chi2={chi_nominal:.0f})") _, dpred_nom, gpred_nom = validate(1.0, 0.2, nP=3) _, dpred_fit, gpred_fit = validate(K_BEST, W_BEST, nP=3) pd.DataFrame({"d_lo": np.arange(0, 300, 25), "N_obs_classI": dobs, "N_model_published_cuts_only": dpred_nom, "N_model_fitted_threshold": dpred_fit}).to_csv( os.path.join(DATA, "validation_distance.csv"), index=False, float_format="%.5g") pd.DataFrame({"G_lo": np.arange(3, 18), "N_obs_classI": gobs, "N_model_published_cuts_only": gpred_nom, "N_model_fitted_threshold": gpred_fit}).to_csv( os.path.join(DATA, "validation_magnitude.csv"), index=False, float_format="%.5g") KTHR, SIGLOG = K_BEST, W_BEST def run_variant(nsig, dms_cut, qname, fm): """Predicted number of dark detections per M2 bin if every MS companion were enclosed, and the observed number of dark candidates, for one analysis variant.""" ret = retention(dms_cut) pred = np.zeros(len(M2_BINS)) pred3 = np.zeros(len(M2_BINS)) calibs = {} for (ip, im), v in fm.items(): selo = ((obsI.period >= P_EDGES[ip]) & (obsI.period < P_EDGES[ip + 1]) & (obsI.M1 >= M1_EDGES[im]) & (obsI.M1 < M1_EDGES[im + 1])) calib = selo.sum() / v["lum"][qname] if v["lum"][qname] > 0 else 0.0 calibs[(ip, im)] = calib for i, b in enumerate(M2_BINS): inb = (M2_GRID >= b[0]) & (M2_GRID < b[1]) pred[i] += calib * ret[im] * v["encl"][qname][inb].sum() * dm2 pred3[i] += calib * ret[im] * v["encl"][qname + "_III"][inb].sum() * dm2 n_obs, _ = count_obs(obs, nsig, dms_cut, "nonI") n_obs3, _ = count_obs(obs, nsig, dms_cut, "III") return pred, pred3, n_obs, n_obs3, calibs rows = [] fms = {} for kthr in (1.0, K_BEST, 6.0): KTHR = kthr for nsig in (2.0, 3.0, 5.0): fm_k = forward_model(ms, parent, smod, lim, nsig=nsig) if kthr == K_BEST: fms[nsig] = fm_k for dms_cut in (None, -0.3, -0.2, -0.1): for qname in ("flat", "MDS17", "rising"): pred, pred3, n_obs, n_obs3, _ = run_variant(nsig, dms_cut, qname, fm_k) for i, b in enumerate(M2_BINS): rows.append({"KTHR": kthr, "nsig": nsig, "dMS_cut": "none" if dms_cut is None else dms_cut, "qdist": qname, "M2_bin": f"{b[0]}-{b[1]}", "N_obs_nonI": n_obs[i], "N_pred_nonI_if_all_enclosed": pred[i], "f95_nonI": poisson_ul(n_obs[i]) / pred[i], "N_obs_classIII": n_obs3[i], "N_pred_classIII_if_all_enclosed": pred3[i], "f95_classIII": poisson_ul(n_obs3[i]) / pred3[i] if pred3[i] > 0 else np.nan}) KTHR = K_BEST var = pd.DataFrame(rows) var.to_csv(os.path.join(DATA, "limits_binaries_variants.csv"), index=False, float_format="%.5g") base = var[(var.KTHR == K_BEST) & (var.nsig == NSIG) & (var.dMS_cut == DMS_CUT) & (var.qdist == "flat")].copy() rng = var[(var.dMS_cut != "none") & (var.KTHR <= K_BEST)].groupby("M2_bin").agg(f95_nonI_min=("f95_nonI", "min"), f95_nonI_max=("f95_nonI", "max"), f95_III_min=("f95_classIII", "min"), f95_III_max=("f95_classIII", "max")) base = base.merge(rng, on="M2_bin") nocut = var[(var.KTHR == K_BEST) & (var.nsig == NSIG) & (var.dMS_cut == "none") & (var.qdist == "flat")][["M2_bin", "N_obs_nonI", "N_pred_nonI_if_all_enclosed", "f95_nonI"]] nocut.columns = ["M2_bin", "N_obs_nonI_nocut", "N_pred_nocut", "f95_nonI_nocut"] base = base.merge(nocut, on="M2_bin") # the same analysis with the published acceptance cuts only (KTHR = 1): weaker limits, quoted as # the conservative headline because it does not rely on the fitted threshold pub = var[(var.KTHR == 1.0) & (var.nsig == NSIG) & (var.dMS_cut == DMS_CUT) & (var.qdist == "flat")][ ["M2_bin", "N_pred_nonI_if_all_enclosed", "f95_nonI", "f95_classIII"]] pub.columns = ["M2_bin", "N_pred_published_cuts_only", "f95_nonI_published_cuts_only", "f95_classIII_published_cuts_only"] base = base.merge(pub, on="M2_bin") base["f95_conservative"] = np.maximum(base.f95_nonI, base.f95_nonI_published_cuts_only) # x = N_dark/N_lum scaled is f/(1-f); the fraction itself is x/(1+x) (quoting x is conservative) base["f95_conservative_as_fraction"] = base.f95_conservative / (1 + base.f95_conservative) print(base.T.to_string()) # ---- baseline candidates, calibration table and white-dwarf expectation ---------------------- fm = fms[NSIG] pred, pred3, n_obs, n_obs3, calibs = run_variant(NSIG, DMS_CUT, "flat", fm) _, selc = count_obs(obs, NSIG, DMS_CUT, "nonI") cand = obs[selc].copy() cand["shahaf24"] = np.where(cand.red_prob.isna(), "not in S24", np.where(cand.red_prob >= RED_THRESH, "red excess", "no red excess")) keep = ["source_id", "nss_solution_type", "ra", "dec", "dist", "phot_g_mean_mag", "bp_rp", "MG", "dMS", "M1", "m1", "period", "eccentricity", "a0", "significance", "A", "sig_A", "A_ms", "A_tr", "nsigI", "nsigIII", "M2_dark", "M2_dark_lo", "M2_dark_hi", "classIII", "PII", "PIII", "red_prob", "shahaf24", "shahaf23_class", "shahaf23_note", "shahaf23_M2min"] cand[keep].to_csv(os.path.join(DATA, "dark_companion_candidates.csv"), index=False, float_format="%.6g") comp_rows = [] for i, b in enumerate(M2_BINS): inb = (cand.M2_dark >= b[0]) & (cand.M2_dark < b[1]) comp_rows.append({"M2_bin": f"{b[0]}-{b[1]}", "N_dark_candidates": int(inb.sum()), "of_which_S24_no_red_excess": int((inb & (cand.shahaf24 == "no red excess")).sum()), "of_which_S24_red_excess": int((inb & (cand.shahaf24 == "red excess")).sum()), "of_which_not_in_S24": int((inb & (cand.shahaf24 == "not in S24")).sum()), "of_which_classIII": int((inb & (cand.nsigIII > NSIG)).sum()), "median_ecc": float(cand[inb].eccentricity.median()), "median_M1": float(cand[inb].M1.median())}) comp = pd.DataFrame(comp_rows) print(comp.to_string(index=False)) # ---- (revision) low-mass candidates and the objects above 1.4 Msun ---------------------------- low_rows = [] for lo, hi in ((0.1, 0.3), (0.3, 0.4), (0.1, 0.4), (0.4, 0.5), (0.1, 0.5), (0.5, 0.8), (0.8, 1.4), (1.4, 99.0)): s_ = cand[(cand.M2_dark >= lo) & (cand.M2_dark < hi)] low_rows.append({ "M2_range": f"{lo}-{hi}", "N": int(len(s_)), "N_classIII_above_close_MS_pair_limit": int((s_.nsigIII > NSIG).sum()), "N_below_close_MS_pair_limit": int((s_.nsigIII <= NSIG).sum()), "N_consistent_with_close_MS_pair_at_1sigma": int((s_.nsigIII <= 1.0).sum()), "N_M2_within_1sigma_of_0.5": int(((s_.M2_dark_hi >= 0.5) & (s_.M2_dark_lo <= 0.5)).sum()) if hi <= 0.5 else np.nan, "N_in_S24": int((s_.shahaf24 != "not in S24").sum()), "N_S24_red_excess": int((s_.shahaf24 == "red excess").sum()), "N_S24_no_red_excess": int((s_.shahaf24 == "no red excess").sum()), "N_q_gt_1": int((s_.M2_dark > s_.M1).sum()), "median_ecc": float(s_.eccentricity.median()) if len(s_) else np.nan, "median_M1": float(s_.M1.median()) if len(s_) else np.nan, "N_in_Shahaf23_table2": int(s_.shahaf23_class.notna().sum())}) low = pd.DataFrame(low_rows) low.to_csv(os.path.join(DATA, "candidate_breakdown.csv"), index=False, float_format="%.4g") print(low.to_string(index=False)) heavy = cand[cand.M2_dark >= 1.4][["source_id", "dist", "M1", "M2_dark", "M2_dark_lo", "M2_dark_hi", "period", "eccentricity", "nsigIII", "shahaf24", "shahaf23_class", "shahaf23_note", "shahaf23_M2min"]] heavy.to_csv(os.path.join(DATA, "candidates_above_1p4.csv"), index=False, float_format="%.6g") print("[census] candidates above 1.4 Msun:") print(heavy.to_string(index=False)) t2 = cand[cand.shahaf23_class.notna()] t2_by = pd.crosstab(pd.cut(t2.M2_dark, [0.1, 0.5, 0.8, 1.4, 99]), t2.shahaf23_class + " " + t2.shahaf23_note) print("[census] candidates in Shahaf+2023 Table 2, by mass bin and class:") print(t2_by.to_string()) ret = retention(DMS_CUT) cal_rows, wd_pred = [], np.zeros(len(M2_BINS)) wd_by_m1 = [] for im in range(len(M1_EDGES) - 1): m1c = 0.5 * (M1_EDGES[im] + M1_EDGES[im + 1]) rwd, dens = wd_expectation(m1c) tot = np.zeros(len(M2_BINS)) for ip in range(len(P_EDGES) - 1): v = fm[(ip, im)] selo = ((obsI.period >= P_EDGES[ip]) & (obsI.period < P_EDGES[ip + 1]) & (obsI.M1 >= M1_EDGES[im]) & (obsI.M1 < M1_EDGES[im + 1])) cal_rows.append({"P_lo": P_EDGES[ip], "P_hi": P_EDGES[ip + 1], "M1_lo": M1_EDGES[im], "M1_hi": M1_EDGES[im + 1], "parent_stars": v["nstars"], "N_obs_classI": int(selo.sum()), "mean_detectability_luminous_flatq": v["lum"]["flat"] / v["nstars"], "companion_freq_x_unmodelled_eff": calibs[(ip, im)], "single_light_retention": ret[im]}) for i, b in enumerate(M2_BINS): inb = (M2_GRID >= b[0]) & (M2_GRID < b[1]) tot[i] += calibs[(ip, im)] * ret[im] * (v["dark"] * dens)[inb].sum() * dm2 wd_pred += tot wd_by_m1.append({"M1_lo": M1_EDGES[im], "M1_hi": M1_EDGES[im + 1], "WD_per_MS_companion": rwd, **{f"N_WD_exp_{b[0]}-{b[1]}": tot[i] for i, b in enumerate(M2_BINS)}}) cal = pd.DataFrame(cal_rows) cal.to_csv(os.path.join(DATA, "selection_calibration.csv"), index=False, float_format="%.5g") wd = pd.DataFrame(wd_by_m1) wd.to_csv(os.path.join(DATA, "wd_expectation.csv"), index=False, float_format="%.4g") print(wd.to_string(index=False)) print("[census] WD expected per M2 bin (no orbital evolution):", wd_pred) byP = cal.groupby("P_lo").apply(lambda x: pd.Series({ "N_obs_classI": x.N_obs_classI.sum(), "freq_x_eff": (x.companion_freq_x_unmodelled_eff * x.parent_stars).sum() / x.parent_stars.sum(), "mean_detectability": (x.mean_detectability_luminous_flatq * x.parent_stars).sum() / x.parent_stars.sum()}), include_groups=False) print(byP) base["N_WD_expected_no_orbital_evolution"] = wd_pred base["N_obs_minus_WD_expected"] = base.N_obs_nonI - wd_pred base.to_csv(os.path.join(DATA, "limits_binaries.csv"), index=False, float_format="%.5g") summary = {"cuts": cuts, "checks": chk, "n_parent_all": n_parent_all, "n_parent_ms": n_parent_ms, "n_obs_sample": int(len(obs)), "n_obs_classI": int(obs.classI.sum()), "baseline": {"nsig": NSIG, "dMS_cut": DMS_CUT, "qdist": "flat", "KTHR": K_BEST, "SIGLOG": W_BEST}, "threshold_fit_chi2_best": float(fit.iloc[0].chi2), "threshold_fit_chi2_published_cuts_only": chi_nominal, "single_light_retention_by_M1_bin": [float(x) for x in ret], "limits_baseline": base.to_dict(orient="records"), "candidate_composition": comp.to_dict(orient="records"), "candidate_breakdown": low.to_dict(orient="records"), "candidates_above_1p4": heavy.to_dict(orient="records"), "n_candidates_in_shahaf23_table2": int(len(t2)), "selection_model_factor_published_over_fitted": [float(x) for x in (base.f95_nonI_published_cuts_only / base.f95_nonI)], "wd_expected_by_bin": [float(x) for x in wd_pred], "freq_x_eff_by_period_bin": byP.reset_index().to_dict(orient="records"), "sigma_scatter_dex_measured": smod[3]} with open(os.path.join(DATA, "binaries_summary.json"), "w") as f: json.dump(summary, f, indent=1, default=float) obs.to_pickle(os.path.join(RAW, "obs_sample.pkl")) # ---- figure 2: AMRF plane and dark-companion mass spectrum ---------------------------------- fig, ax = plt.subplots(1, 2, figsize=(12.5, 4.8)) ax[0].scatter(obs.M1, obs.A, s=1, c="0.78", rasterized=True, label=f"all {len(obs)} orbits (d < 300 pc)") oth = obs[(obs.nsigI > NSIG) & ~selc] ax[0].scatter(oth.M1, oth.A, s=3, c="#ff7f0e", rasterized=True, label=f"above single-star limit, over-luminous ({len(oth)})") ax[0].scatter(cand.M1, cand.A, s=3, c="#1f77b4", rasterized=True, label=f"dark-companion candidates ({len(cand)})") ax[0].plot(grid, a_ms, "k--", lw=1.2, label="maximum for one main-sequence companion") ax[0].plot(grid, a_tr, "k-", lw=1.2, label="maximum for a close main-sequence pair") for m2 in (0.3, 0.6, 1.0, 1.4): ax[0].plot(grid, amrf_dark(m2 / grid), color="#2ca02c", lw=0.7) ax[0].text(1.505, amrf_dark(m2 / 1.5), f"{m2}", color="#2ca02c", fontsize=8, va="center") ax[0].set_xlim(0.2, 1.6) ax[0].set_ylim(0, 1.3) ax[0].set_xlabel("primary mass M$_1$ [M$_\\odot$]") ax[0].set_ylabel("astrometric mass-ratio function A") ax[0].legend(fontsize=7, loc="upper right") ax[0].set_title("Green: locus of a dark companion of the labelled mass [M$_\\odot$]", fontsize=9) bins = np.arange(0.1, 1.6, 0.05) ax[1].hist(cand.M2_dark, bins=bins, color="#1f77b4", label="observed dark-companion candidates") scale_enc = 0.1 encl_tot = np.zeros(len(M2_GRID)) wd_tot = np.zeros(len(M2_GRID)) for (ip, im), v in fm.items(): encl_tot += calibs[(ip, im)] * ret[im] * v["encl"]["flat"] m1c = 0.5 * (M1_EDGES[im] + M1_EDGES[im + 1]) wd_tot += calibs[(ip, im)] * ret[im] * v["dark"] * wd_expectation(m1c)[1] ax[1].plot(M2_GRID, encl_tot * 0.05 * scale_enc, color="#d62728", lw=1.5, label="expected if 10% of main-sequence companions were enclosed") kern = np.exp(-0.5 * (np.arange(-15, 16) * dm2 / 0.05) ** 2) wd_s = np.convolve(wd_tot, kern / kern.sum(), mode="same") ax[1].plot(M2_GRID, wd_s * 0.05, color="k", lw=1.2, ls="--", label="white dwarfs expected with no orbital evolution") for b in M2_BINS[:N_MAIN_BINS]: ax[1].axvline(b[0], color="k", ls=":", lw=0.6) ax[1].axvline(M2_BINS[N_MAIN_BINS - 1][1], color="k", ls=":", lw=0.6) ax[1].axvline(0.4, color="0.5", ls="--", lw=0.6) c3 = cand[cand.nsigIII > NSIG] ax[1].hist(c3.M2_dark, bins=bins, color="#17becf", label="of which above the limit for a close main-sequence pair (class III)") ax[1].set_xlabel("companion mass if optically dark [M$_\\odot$]") ax[1].set_ylabel("systems per 0.05 M$_\\odot$") ax[1].set_xlim(0.1, 1.5) ax[1].set_ylim(0, 1250) ax[1].legend(fontsize=7.5, loc="upper right") ax[1].set_title("Mass spectrum of the dark-companion candidates", fontsize=9) fig.tight_layout() fig.savefig(os.path.join(FIGS, "fig2_amrf_dark_companions.png"), dpi=150) plt.close(fig) if __name__ == "__main__": main()