"""Synthesis: one table and one figure collecting, for every handle, the upper limit on the enclosed fraction, the population it refers to, the mass range, the volume, the enclosure temperatures for which the handle applies, its status and its main systematic. Revised after the referee report: * every row states WHICH fraction it bounds (companions at ~1 au, local stellar mass, published planet-host lenses, stars in high-latitude cones); they are not one quantity; * rows are tagged limit / ceiling from counts / indicative / forecast / not applicable, and the figure separates limits (left panel) from everything else (right panel); * the optical-depth 'ceiling' is withdrawn; * the published infrared limits (Suazo et al. 2022) are shown as what they are: limits on PARTIAL spheres around stars with an optical counterpart, not applicable to complete enclosures. """ 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 DATA, FIGS def load(name): with open(os.path.join(DATA, name)) as f: return json.load(f) def main(): ir = load("ir_detectability_summary.json") bn = load("binaries_summary.json") ws = load("wise_thermal_summary.json") oo = load("oort_summary.json") ml = load("microlensing_summary.json") wd = load("w4dark_summary.json") w4 = pd.read_csv(os.path.join(DATA, "w4dark_ceiling_by_T.csv")) lim = {r["M2_bin"]: r for r in bn["limits_baseline"]} tmin = ws["baseline_T_min_percentiles_5_50_95"] rows = [] for b in ("0.1-0.5", "0.5-0.8", "0.8-1.4", "0.1-0.4"): r = lim[b] if b == "0.5-0.8": sysm = "white-dwarf degeneracy (limit = white dwarfs + enclosures); selection model (factor ~2)" elif b == "0.1-0.5": sysm = "selection model (factor ~2); upper edge contains the low-mass wing of the white-dwarf peak" elif b == "0.1-0.4": sysm = "selection model (factor ~2); sub-bin of 0.1-0.5 chosen to be free of white dwarfs" else: sysm = "selection model (factor ~2)" rows.append({ "handle": f"Gaia DR3 astrometric binaries, companions {b} Msun" + (" (sub-bin)" if b == "0.1-0.4" else ""), "status": "limit", "f95": r["f95_conservative"], "f95_alt": r["f95_nonI"], "alt_meaning": "with detection threshold fitted to the class-I distance/magnitude distributions", "range_over_variants": f"{r['f95_nonI_min']:.3f}-{r['f95_nonI_max']:.3f}", "N_obs": int(r["N_obs_nonI"]), "N_expected_if_all_enclosed": r["N_pred_published_cuts_only"], "population": "stellar companions (less massive than the visible star) at P = 100-1000 d of 0.2-1.5 Msun " "main-sequence primaries", "mass_range_Msun": b, "volume": "d < 300 pc", "enclosure_T_K": f"any emission law; a thermal enclosure must be hotter than {tmin[0]:.0f}-{tmin[2]:.0f} K " "(Roche lobe) and is then excluded by AllWISE", "main_systematic": sysm}) rows.append({ "handle": "Same candidates, thermally radiating enclosure (AllWISE)", "status": "test", "f95": np.nan, "f95_alt": np.nan, "alt_meaning": "", "range_over_variants": "", "N_obs": ws["baseline_n_not_excluded_with_wise"], "N_expected_if_all_enclosed": ws["baseline_n_with_allwise"], "population": f"{ws['baseline_n_thermal_excluded']} of {ws['baseline_n_with_allwise']} candidates with AllWISE " f"excluded (lowest {ws['baseline_min_nsig_thermal']:.1f} sigma, evaluated at the {ws['T_cap_K']:.0f} K cap)", "mass_range_Msun": "0.1-1.4", "volume": "d < 300 pc", "enclosure_T_K": f"{tmin[0]:.0f}-{ws['T_cap_K']:.0f}, blackbody", "main_systematic": "assumes isotropic blackbody re-emission of the full luminosity; hotter shells not tested"}) b = oo["baseline"] rows.append({ "handle": "Local mass budget (Oort limit): dynamical disc - visible baryons", "status": "limit", "f95": b["f95_enclosed_fraction_of_stars"], "f95_alt": oo["S3_f95_range"][1], "alt_meaning": "if there were no dark-matter halo within 1.1 kpc", "range_over_variants": f"{oo['S1_f95_range'][0]:.2f}-{oo['S1_f95_range'][1]:.2f}", "N_obs": np.nan, "N_expected_if_all_enclosed": np.nan, "population": f"local stellar mass (Sigma_enc < {b['UL95_positive_prior']:.0f} Msun/pc2 against 25.4 visible)", "mass_range_Msun": "all (mass-weighted)", "volume": "column |z| < 1.1 kpc at the Sun", "enclosure_T_K": "any", "main_systematic": "pre-Gaia dynamical determinations only; baryon inventory (+-3.4 Msun/pc2); disc-halo split; " "disequilibrium"}) for cname, short in (("Sun-like 0.8-1.2 Msun", "Sun-like stars"), ("M dwarfs 0.2-0.6 Msun", "M dwarfs")): c = wd["ceiling_f95"][cname] c2 = wd["ceiling_f95_thin_disc_only"][cname] for T in ("50", "100", "200"): rows.append({ "handle": f"AllWISE source counts: 22-um-bright, 3.4-um-dark sources ({short}, {T} K)", "status": "ceiling from counts", "f95": c[T], "f95_alt": c2[T], "alt_meaning": "thin disc only, scale height 250 pc", "range_over_variants": "", "N_obs": wd["N_obs_by_T"][T], "N_expected_if_all_enclosed": (wd["N_stars_sunlike_by_T"] if short.startswith("Sun") else wd["N_stars_mdwarf_by_T"])[T], "population": f"{cname} stars in {wd['area_deg2']:.0f} deg2 at |b| > 31 deg whose enclosure would exceed " f"W4 = {wd['selection']['W4_lim_mag']} mag", "mass_range_Msun": cname.split()[-2], "volume": "out to the distance where W4 = 7.0 (see table)", "enclosure_T_K": f"{T} (isotropic single-temperature blackbody)", "main_systematic": "blackbody assumption; Galactic disc model; assumed completeness 0.8; counts are not vetted"}) ll = ml["lens_light"] rows.append({ "handle": "Microlensing: published lenses with measured mass and high-resolution imaging", "status": "indicative", "f95": ll["f95_dark_conservative"], "f95_alt": ll["f95_dark_nominal"], "alt_meaning": "nominal: ambiguous events dropped (2 dark of 28)", "range_over_variants": f"{ll['f95_dark_nominal']:.2f}-{ll['f95_dark_most_conservative']:.2f}", "N_obs": ll["n_dark_conservative"], "N_expected_if_all_enclosed": ll["n_total_conservative"], "population": "published follow-ups of (mostly planet-host) lenses; not a representative sample of stars", "mass_range_Msun": f"{ll['luminous_mass_range'][0]}-{ll['luminous_mass_range'][1]}", "volume": f"lenses at {ll['luminous_DL_range_kpc'][0]:.1f}-{ll['luminous_DL_range_kpc'][1]:.1f} kpc towards the bulge", "enclosure_T_K": "about 100-450 (R_s < 0.6 r_E and dark at 2 um)", "main_systematic": "hand-compiled sample, probably biased towards luminous lenses"}) fc = ml["occultation_forecast"] rows.append({ "handle": "Occultations of bulge stars by cold enclosures", "status": "forecast", "f95": fc["f95_if_zero_found_eff10pct_smith_norm"], "f95_alt": fc["f95_if_zero_found_eff10pct_rate_norm"], "alt_meaning": "normalised to the measured microlensing rate instead of Smith et al. 2021", "range_over_variants": "", "N_obs": np.nan, "N_expected_if_all_enclosed": fc["reference_N_per_unit_f_smith_norm"], "population": "stars along bulge sightlines; reach of a null search with an ASSUMED 10% efficiency", "mass_range_Msun": "all", "volume": "kpc-scale sightlines to the bulge", "enclosure_T_K": "< about 50-150 (R_s > r_E)", "main_systematic": "no efficiency-calibrated search exists; survey totals are order-of-magnitude"}) rows.append({ "handle": "Published infrared limit, Suazo et al. 2022 (for comparison)", "status": "not applicable to complete enclosures", "f95": 2e-5, "f95_alt": 8e-4, "alt_meaning": "within 5 kpc instead of 100 pc", "range_over_variants": "", "N_obs": np.nan, "N_expected_if_all_enclosed": np.nan, "population": "stars WITH a Gaia counterpart (G <= 21) hosting a partial sphere covering 90% (limits quoted for 1-90%)", "mass_range_Msun": "-", "volume": "100 pc / 5 kpc", "enclosure_T_K": "~300 (100-1000)", "main_systematic": "a complete enclosure has no Gaia source and is outside the sample by construction"}) tab = pd.DataFrame(rows) tab.to_csv(os.path.join(DATA, "synthesis_limits.csv"), index=False, float_format="%.4g") pd.set_option("display.width", 250) print(tab[["handle", "status", "f95", "f95_alt", "mass_range_Msun", "volume"]].to_string(index=False)) # ---- figure --------------------------------------------------------------------------------- fig, axes = plt.subplots(1, 2, figsize=(14.5, 7.8), sharey=True) T_lo = ir["Tmax_hidden_beyond_100pc_L1_WISE_bandavg"] T_hi = ir["Tmax_hidden_beyond_100pc_L1_WISE"] for ax in axes: ax.axvspan(10, T_hi, color="0.9", zorder=0) ax.set_xscale("log") ax.set_yscale("log") ax.set_xlim(10, 1500) ax.set_ylim(1e-6, 2.0) ax.set_xlabel("temperature of the enclosure if it radiates as a blackbody [K]") ax.grid(alpha=0.25, which="both") axes[0].text(10.6, 1.6e-6, f"1 L$_\\odot$ sphere below the AllWISE\n22 um flux limit beyond 100 pc\n(T < {T_lo:.0f}-{T_hi:.0f} K)", fontsize=8) def bar(ax, x0, x1, y, color, label, ls="-", lw=3, yerr=None): ax.plot([x0, x1], [y, y], color=color, lw=lw, ls=ls, label=label, solid_capstyle="butt") if yerr is not None: xm = np.sqrt(x0 * x1) ax.plot([xm, xm], yerr, color=color, lw=1) ax = axes[0] bar(ax, 10, 1500, b["f95_enclosed_fraction_of_stars"], "#2ca02c", f"Oort limit: fraction of LOCAL STELLAR MASS < {b['f95_enclosed_fraction_of_stars']:.2f} (any emission)", yerr=(oo["S1_f95_range"][0], oo["S3_f95_range"][1])) cols = {"0.1-0.5": "#1f77b4", "0.5-0.8": "#17becf", "0.8-1.4": "#9467bd", "0.1-0.4": "#7f7f7f"} for k in ("0.5-0.8", "0.1-0.5", "0.8-1.4", "0.1-0.4"): r = lim[k] bar(ax, tmin[0], 1500, r["f95_conservative"], cols[k], f"Gaia binaries: fraction of ~1 au COMPANIONS of {k} M$_\\odot$ < {r['f95_conservative']:.3f}" + (" (= white dwarfs)" if k == "0.5-0.8" else "") + (" (sub-bin)" if k == "0.1-0.4" else ""), yerr=(r["f95_nonI"], r["f95_conservative"]), lw=2.5 if k != "0.1-0.4" else 1.5) Tg = w4.T_K.values base = "Juric+2008 thin+thick" alt = "thin disc only, h = 250 pc" for cname, c, lab in (("Sun-like 0.8-1.2 Msun", "#d62728", "Sun-like stars"), ("M dwarfs 0.2-0.6 Msun", "#ff7f0e", "M dwarfs")): y = w4[f"f95 | {cname} | {base}"].values y2 = w4[f"f95 | {cname} | {alt}"].values ok = y < 1.5 ax.plot(Tg[ok], y[ok], "o-", color=c, ms=3.5, lw=1.8, label=f"AllWISE source-count ceiling: fraction of {lab.upper()} (blackbody enclosures only)") ax.fill_between(Tg[ok], y[ok], np.minimum(y2[ok], 2.0), color=c, alpha=0.2, lw=0) ax.set_ylabel("95% upper limit on the enclosed fraction\n(the population differs by handle; see legend)") ax.set_title("Limits and count-based ceilings", fontsize=10) ax.legend(fontsize=7.4, loc="upper center", bbox_to_anchor=(0.5, -0.1), ncol=1, framealpha=0.95) ax = axes[1] ax.plot([100, 450], [ll["f95_dark_conservative"]] * 2, color="#8c564b", lw=2, ls="--", marker="o", mfc="none", label=f"INDICATIVE: dark fraction of published planet-host microlenses < {ll['f95_dark_conservative']:.2f} " f"({ll['n_dark_conservative']} of {ll['n_total_conservative']}; nominal {ll['f95_dark_nominal']:.2f})") lo_, hi_ = sorted([fc["f95_if_zero_found_eff10pct_rate_norm"], fc["f95_if_zero_found_eff10pct_smith_norm"]]) ax.fill_between([10, 100], [lo_, lo_], [hi_, hi_], facecolor="none", edgecolor="#d62728", hatch="//", lw=1, label="FORECAST (not a limit): reach of a null occultation search in OGLE-IV\n" "if the detection efficiency were 10% (two normalisations)") ax.plot([100, 1000], [2e-5, 2e-5], color="0.35", lw=1.5, ls="-.", label="NOT APPLICABLE to complete enclosures: Suazo+2022, partial spheres (90% covering)\n" "around Gaia stars, < 2e-5 within 100 pc") ax.plot([100, 1000], [8e-4, 8e-4], color="0.35", lw=1.0, ls=":", label="same, < 8e-4 within 5 kpc") ax.set_title("Indicative values, forecasts and published limits on a different population", fontsize=10) ax.legend(fontsize=7.4, loc="upper center", bbox_to_anchor=(0.5, -0.1), ncol=1, framealpha=0.95) fig.suptitle("Upper limits on fully enclosed stars: what each handle bounds, and for which enclosure temperatures", fontsize=11) fig.tight_layout() fig.savefig(os.path.join(FIGS, "fig5_synthesis.png"), dpi=150) plt.close(fig) if __name__ == "__main__": main()