"""Thermal test of the dark-companion candidates with AllWISE. An enclosure around the companion of an astrometric binary must fit inside the companion's Roche lobe at periastron, so it cannot be cold: R_s <= R_L -> T >= T_min = T_1au sqrt(sqrt(L)/R_L). If it re-radiates the companion's luminosity L(M2) isotropically as a blackbody at any T_min <= T <= 1500 K, the system must show an excess in G-W1, G-W2, G-W3 or G-W4 relative to the single-star locus. For each candidate we find the temperature that hides best and ask whether even that is excluded at > 3 sigma in at least one band. Scope (revision): the test covers blackbody enclosures from T_min up to T_CAP = 1500 K. For most candidates the best-hiding temperature IS the cap, so the quoted significances are significances at 1500 K. A hotter shell is smaller than about 0.07 au sqrt(L/Lsun), radiates mainly at 1-3 um and is in effect a cool star: it is tested by the optical/near-infrared colours (the Shahaf et al. 2024 red-excess test), not here. """ import json import os import sys import numpy as np import pandas as pd 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, H_PLANCK, C_LIGHT, K_B, SIGMA_SB, L_SUN, AU, PC) from mslib import MainSequence WISE = {"w1": (3.3526, 309.540), "w2": (4.6028, 171.787), "w3": (11.5608, 31.674), "w4": (22.0883, 8.363)} T1AU = (L_SUN / (4 * np.pi * SIGMA_SB * AU ** 2)) ** 0.25 T_CAP = 1500.0 def bnu(lam_um, T): nu = C_LIGHT / (lam_um * 1e-6) return 2 * H_PLANCK * nu ** 3 / C_LIGHT ** 2 / np.expm1(H_PLANCK * nu / (K_B * T)) def enclosure_flux_jy(lam_um, T, L, d_pc): R = np.sqrt(L * L_SUN / (4 * np.pi * SIGMA_SB * T ** 4)) return np.pi * bnu(lam_um, T) * (R / (d_pc * PC)) ** 2 / 1e-26 def roche_lobe(q): """Eggleton (1983) Roche-lobe radius of the star with mass ratio q = M_this/M_other, in units of a.""" q3 = q ** (1 / 3) return 0.49 * q3 ** 2 / (0.6 * q3 ** 2 + np.log1p(q3)) def running_locus(x, y, grid, half=0.15, minn=30): med, sig = np.full(len(grid), np.nan), np.full(len(grid), np.nan) for i, g in enumerate(grid): s = y[(x > g - half) & (x < g + half)] s = s[np.isfinite(s)] if len(s) >= minn: med[i] = np.median(s) sig[i] = 1.4826 * np.median(np.abs(s - med[i])) ok = np.isfinite(med) return grid[ok], med[ok], sig[ok] def main(): ms = MainSequence() obs = pd.read_pickle(os.path.join(RAW, "obs_sample.pkl")) w = pd.read_csv(os.path.join(RAW, "nss_allwise.csv")).sort_values("angular_distance").drop_duplicates("source_id") o = obs.merge(w, on="source_id", how="left") o["has_wise"] = o.w1mpro.notna() bands = ["w1", "w2", "w3", "w4"] for i, b in enumerate(bands): q = o.ph_qual.fillna("XXXX").str[i] o[b + "_ok"] = q.isin(["A", "B"]) & o[b + "mpro"].notna() & o[b + "mpro_error"].notna() o["g_" + b] = o.phot_g_mean_mag - o[b + "mpro"] # empirical single-light locus of G - W_k versus BP-RP, from all orbits in the sample grid = np.arange(0.4, 3.61, 0.05) locus = {} for b in bands: s = o[o[b + "_ok"]] gx, med, sig = running_locus(s.bp_rp.values, s["g_" + b].values, grid) locus[b] = (gx, med, sig) o["exc_" + b] = o["g_" + b] - np.interp(o.bp_rp, gx, med) o["sig_" + b] = np.sqrt(o[b + "mpro_error"].fillna(0.3) ** 2 + np.interp(o.bp_rp, gx, sig) ** 2) cand = o[o.nonI].copy() # all orbits above the single-star AMRF limit at 3 sigma cand["single_light"] = cand.dMS > -0.2 # the baseline dark-companion candidates of 03_dark_census # enclosure geometry a_au = ((cand.M1 + cand.M2_dark) * (cand.period / 365.25) ** 2) ** (1 / 3) cand["a_au"] = a_au cand["R_max_au"] = a_au * (1 - cand.eccentricity) * roche_lobe(cand.M2_dark / cand.M1) cand["L2"] = 10 ** ms.logL_of_mass(np.clip(cand.M2_dark, 0.075, 3.0)) cand["T_min"] = T1AU * np.sqrt(np.sqrt(cand.L2) / cand.R_max_au) Tgrid = np.logspace(np.log10(100), np.log10(T_CAP), 60) nsig_best = np.full(len(cand), np.nan) # the smallest (over T) of the largest (over bands) discrepancy T_best = np.full(len(cand), np.nan) pred_w3_at_tmin = np.full(len(cand), np.nan) for j, (_, r) in enumerate(cand.iterrows()): if not r.has_wise: continue worst = np.inf for T in Tgrid[Tgrid >= r.T_min]: best_band = -np.inf for b in bands: if not r[b + "_ok"] and b != "w4": continue lam, f0 = WISE[b] fe = enclosure_flux_jy(lam, T, r.L2, r.dist) gx, med, _ = locus[b] # photospheric flux predicted from G and the locus colour mphot = r.phot_g_mean_mag - np.interp(r.bp_rp, gx, med) fp = f0 * 10 ** (-0.4 * mphot) pred_exc = 2.5 * np.log10(1 + fe / fp) if r[b + "_ok"]: nsig = (pred_exc - r["exc_" + b]) / r["sig_" + b] else: # W4 not detected: the enclosure alone would have to be fainter than the # AllWISE 5-sigma limit (5 mJy); count as 5 sigma x (flux / limit) nsig = 5.0 * fe / 5.0e-3 if fe > 5.0e-3 else -np.inf best_band = max(best_band, nsig) if best_band < worst: worst, tb = best_band, T nsig_best[j] = worst T_best[j] = tb lam, f0 = WISE["w3"] gx, med, _ = locus["w3"] fp = f0 * 10 ** (-0.4 * (r.phot_g_mean_mag - np.interp(r.bp_rp, gx, med))) pred_w3_at_tmin[j] = 2.5 * np.log10(1 + enclosure_flux_jy(lam, max(r.T_min, 100.0), r.L2, r.dist) / fp) cand["nsig_thermal_min"] = nsig_best cand["T_hides_best"] = T_best cand["pred_W3_excess_at_Tmin"] = pred_w3_at_tmin cand["thermal_excluded"] = cand.nsig_thermal_min > 3 keep = ["source_id", "dist", "phot_g_mean_mag", "bp_rp", "M1", "M2_dark", "period", "eccentricity", "a_au", "R_max_au", "L2", "T_min", "w1mpro", "w2mpro", "w3mpro", "w4mpro", "ph_qual", "exc_w1", "exc_w2", "exc_w3", "exc_w4", "sig_w3", "pred_W3_excess_at_Tmin", "nsig_thermal_min", "T_hides_best", "thermal_excluded", "classIII", "red_prob", "dMS", "single_light"] cand[keep].to_csv(os.path.join(DATA, "wise_thermal_test.csv"), index=False, float_format="%.5g") bl = cand[cand.single_light] blw = bl[bl.has_wise] summ = { "T_cap_K": T_CAP, "baseline_T_min_percentiles_5_50_95": [float(x) for x in np.nanpercentile(bl.T_min, [5, 50, 95])], "baseline_R_max_au_percentiles_5_50_95": [float(x) for x in np.nanpercentile(bl.R_max_au, [5, 50, 95])], "baseline_pred_W3_excess_mag_percentiles_5_50_95": [float(x) for x in np.nanpercentile( blw.pred_W3_excess_at_Tmin.dropna(), [5, 50, 95])], "baseline_obs_W3_excess_mag_percentiles_5_50_95": [float(x) for x in np.nanpercentile(bl.exc_w3[bl.w3_ok], [5, 50, 95])], "baseline_n_best_hiding_T_at_cap": int((blw.T_hides_best >= T_CAP * 0.999).sum()), "baseline_frac_best_hiding_T_at_cap": float((blw.T_hides_best >= T_CAP * 0.999).mean()), "baseline_min_nsig_thermal": float(blw.nsig_thermal_min.min()), "baseline_median_nsig_thermal": float(blw.nsig_thermal_min.median()), "baseline_n_excluded_gt5sigma": int((blw.nsig_thermal_min > 5).sum()), "n_candidates": int(len(cand)), "n_with_allwise": int(cand.has_wise.sum()), "n_thermal_excluded": int(cand.thermal_excluded.sum()), "n_not_excluded_with_wise": int((cand.has_wise & ~cand.thermal_excluded).sum()), "baseline_n_candidates": int(cand.single_light.sum()), "baseline_n_with_allwise": int((cand.single_light & cand.has_wise).sum()), "baseline_n_thermal_excluded": int((cand.single_light & cand.thermal_excluded).sum()), "baseline_n_not_excluded_with_wise": int((cand.single_light & cand.has_wise & ~cand.thermal_excluded).sum()), "baseline_n_without_allwise": int((cand.single_light & ~cand.has_wise).sum()), "baseline_min_nsig_thermal_percentiles_5_50": [float(x) for x in np.nanpercentile( cand.nsig_thermal_min[cand.single_light & cand.has_wise], [5, 50])], "T_min_percentiles_5_50_95": [float(x) for x in np.nanpercentile(cand.T_min, [5, 50, 95])], "R_max_au_percentiles_5_50_95": [float(x) for x in np.nanpercentile(cand.R_max_au, [5, 50, 95])], "pred_W3_excess_mag_percentiles_5_50_95": [float(x) for x in np.nanpercentile(cand.pred_W3_excess_at_Tmin.dropna(), [5, 50, 95])], "obs_W3_excess_mag_percentiles_5_50_95_candidates": [float(x) for x in np.nanpercentile(cand.exc_w3[cand.w3_ok], [5, 50, 95])], "obs_W3_excess_mag_percentiles_5_50_95_all": [float(x) for x in np.nanpercentile(o.exc_w3[o.w3_ok], [5, 50, 95])], "n_candidates_W3_excess_gt_1mag": int(((cand.exc_w3 > 1.0) & cand.w3_ok).sum()), "n_all_W3_excess_gt_1mag": int(((o.exc_w3 > 1.0) & o.w3_ok).sum()), "n_all_with_w3": int(o.w3_ok.sum()), } with open(os.path.join(DATA, "wise_thermal_summary.json"), "w") as f: json.dump(summ, f, indent=1) print("[wise]", json.dumps(summ, indent=1)) ne = cand[cand.has_wise & ~cand.thermal_excluded] if len(ne): print(ne[keep].to_string()) fig, ax = plt.subplots(1, 2, figsize=(11.5, 4.4)) s = o[o.w3_ok] ax[0].scatter(s.bp_rp, s.g_w3, s=1, c="0.8", rasterized=True, label="all orbits with W3") c = cand[cand.w3_ok] ax[0].scatter(c.bp_rp, c.g_w3, s=4, c="#1f77b4", label="dark-companion candidates (observed)") pred = c.g_w3 - c.exc_w3 + c.pred_W3_excess_at_Tmin ax[0].scatter(c.bp_rp, pred, s=4, c="#d62728", label="same systems if the companion were an\nenclosed star radiating at T$_{min}$") ax[0].set_xlabel("BP - RP") ax[0].set_ylabel("G - W3") ax[0].set_xlim(0.4, 3.2) ax[0].legend(fontsize=8) ax[0].set_title("12 um test of the dark companions", fontsize=10) ax[1].hist(cand[cand.single_light].T_min.clip(0, 1500), bins=40, color="#8c564b") ax[1].set_xlabel("minimum enclosure temperature allowed by the Roche lobe [K]") ax[1].set_ylabel("baseline dark-companion candidates") ax[1].set_title("Enclosures inside ~1 au binaries cannot be cold", fontsize=10) fig.tight_layout() fig.savefig(os.path.join(FIGS, "fig3_wise_thermal_test.png"), dpi=150) plt.close(fig) if __name__ == "__main__": main()