"""Part 3: expected numbers of total occultations, 95% upper limits on the enclosed fraction f from the OGLE-II null result, the conditional MACHO (Drake & Cook 2003) limit, forecasts for every other existing data set, and the comparison with waste-heat and gravity limits. Expected events in OGLE-II for enclosed fraction f: N = f * sum_fields sum_magbins N_stars(field, bin) * eps_blend(bin) * eps_fg(field) * eps_incl(bin) * norm * sum_pairs w_rate * R_eff * T_eff(duration; bin) with (w_rate, R_eff, duration) from the Galactic model (galmodel.Sightline), T_eff from the end-to-end injection-recovery of 04_ogle2_search.py (including the requirement that the star be present in the reference image), eps_incl the probability that the star is in the catalogue of candidate variables with its centroid on a reference-image star (04b_catalogue_inclusion.py), and the remaining factors defined in CASES below. Outputs: data/blend_efficiency.csv, data/ogle2_expected_events.csv, data/limits_f_vs_R.csv, data/limits_f_vs_T.csv, data/robustness.csv, data/forecast_events.csv, data/limits_summary.json, figures/fig4_limits.png """ import json import os import numpy as np import pandas as pd import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from common import (DATA, FIGS, RAW, poisson_ul, radius_au, H_PLANCK, C_LIGHT, K_B, AU, PC) import galmodel as g import popmodel as pm R_GRID = [1, 2, 3, 5, 7, 10, 20, 50, 100, 200, 500, 1000] T_GRID = [20, 30, 40, 50, 75, 100, 150, 200, 300] POPS = {"all MS stars 0.08-1 Msun": (pm.M_MS_LO, pm.M_MS_HI), "Sun-like 0.7-1 Msun": (0.7, 1.0)} N_SAMP, N_MASS = 20000, 16 MAG_EDGES = [0, 15, 16, 17, 17.5, 18] BINS = [f"{MAG_EDGES[k]}_{MAG_EDGES[k + 1]}" for k in range(len(MAG_EDGES) - 1)] OGLE2_AREA_ARCSEC2 = 49 * 2048 * 8192 * 0.417 ** 2 # 49 drift-scan fields, 0.417 arcsec pixels N_SEASONS_OGLE2 = 4 # Efficiency cases. Each factor is described in the report; "conservative" gives the weaker limit. # incl probability that an occulted star is in the catalogue of candidate variables and passes the # centroid cut C2. "eb": the fraction of independently catalogued deep (> 0.75 mag) eclipsing # binaries of the same magnitude that the catalogue contains with ID_RAD <= 1.5 px (measured, # 04b_catalogue_inclusion.py). "min": the smaller of (eb minus two binomial sigma) and the # first-version estimate, 174/214 (recovery of known microlensing events, Wozniak et al. 2001) # times the fraction of all catalogue objects of that magnitude with ID_RAD <= 1.5 px. # teff which effective exposure: star in the reference image at I_ref < 18.5 or < 19.5. # ramp_bright apply the measured loss for a 4-day ingress and egress to the stars with I < 15 (giants). CASES = { "conservative": dict( incl="min", teff="ref185", ramp_bright=True, norm="p16", # 16th percentile of observed/model OGLE-IV microlensing rate over 81 fields fg="blue", # every star bluer than the clump by > 0.4 mag is a foreground star with zero rate blend=dict(r_blend=1.3, i_to=20.0)), # everything within one seeing FWHM (1.3 arcsec) is merged; LF steep to I = 20 "optimistic": dict( incl="eb", teff="ref195", ramp_bright=False, norm="weighted", # event-weighted observed/model ratio fg="none", blend=dict(r_blend=0.65, i_to=19.0)), # everything within half a FWHM is merged } EPS_CAT_MICROLENSING = 174 / 214 X_MIN_BASE, X_MIN_VARIANT = 0.6, 0.5 # minimum share of the catalogued flux the occulted star must supply ID_FRAC = {} # per magnitude bin: fraction of catalogue objects passing the centroid cut C2 (filled in main) INCL = {} # case -> per magnitude bin inclusion probability (filled in main) # ------------------------------------------------------------------------------------------------ def blend_mc(counts, r_blend, i_to, x_min, n=300000, seed=5, i_sky=20.0): """Number of real stars that supply at least x_min of the flux of their catalogue object, divided by the number of catalogue objects, per magnitude bin of the object. Monte Carlo of random superpositions: every star within r_blend [arcsec] of a target star is merged with it into one catalogue object (stars farther away are separate objects). Stars fainter than i_sky are below the detection limit; their mean surface brightness is absorbed in the local sky, so only its fluctuation is added. Luminosity function: OGLE-II map counts to I = 18.5, continued with the slope measured between I = 17 and 18.5 down to the turn-off magnitude i_to, then with slope 0.3 dex/mag to I = 24.5.""" rng = np.random.default_rng(seed) tot = counts.sum(numeric_only=True) dens = lambda col, width: tot[col] / OGLE2_AREA_ARCSEC2 / width # per arcsec^2 per mag edges = np.arange(12.0, 24.5001, 0.25) mid = 0.5 * (edges[1:] + edges[:-1]) lf = np.zeros_like(mid) for lo, hi, col in ((12, 15, "n_I_0_15"), (15, 16, "n_I_15_16"), (16, 17, "n_I_16_17"), (17, 17.5, "n_I_17_17.5"), (17.5, 18, "n_I_17.5_18"), (18, 18.5, "n_I_18_18.5")): lf[(mid > lo) & (mid < hi)] = dens(col, hi - lo) s1 = np.log10(dens("n_I_18_18.5", 0.5) / dens("n_I_17_17.5", 0.5)) / 1.0 d185 = dens("n_I_18_18.5", 0.5) f1 = mid > 18.5 lf[f1] = d185 * 10 ** (s1 * (np.minimum(mid[f1], i_to) - 18.25)) * 10 ** (0.3 * np.clip(mid[f1] - i_to, 0, None)) flux = 10 ** (-0.4 * (mid - 18.0)) area = np.pi * r_blend ** 2 lam = lf * 0.25 * area # mean number of neighbours per bin inside the blend radius psel = mid < 19.75 # targets: every real star brighter than I = 19.75 ip = rng.choice(np.where(psel)[0], n, p=lf[psel] / lf[psel].sum()) fp = flux[ip] cnt = rng.poisson(lam, (n, len(mid))) fnb = cnt @ flux brightest_nb = np.where(cnt.any(axis=1), flux[np.argmax(cnt > 0, axis=1)], 0.0) faint = mid > i_sky total = np.clip(fp + fnb - float(lam[faint] @ flux[faint]), 1e-6, None) itot = 18.0 - 2.5 * np.log10(total) frac = fp / total is_brightest = fp >= brightest_nb out = {} for k in range(len(MAG_EDGES) - 1): inb = (itot >= max(MAG_EDGES[k], 12.0)) & (itot < MAG_EDGES[k + 1]) den = (is_brightest & inb).sum() out[BINS[k]] = float(((frac >= x_min) & is_brightest & inb).sum() / den) if den > 50 else 1.0 inb = itot < 18 out["all_I_lt_18"] = float(((frac >= x_min) & is_brightest & inb).sum() / (is_brightest & inb).sum()) out["lf_slope_17_18.5"] = float(s1) out["stars_per_arcsec2_to_I24.5"] = float(lf.sum() * 0.25) out["mean_neighbours_brighter_than_I20_in_blend_area"] = float(lam[~faint].sum()) return out # ------------------------------------------------------------------------------------------------ def load_teff(): """Effective exposures from the injection runs: key (low_frac, blend, ramp) -> condition -> dict(logd, all, bins).""" e = pd.read_csv(os.path.join(DATA, "ogle2_efficiency.csv")) out = {} for (lf, bl, rp), s in e.groupby(["low_frac", "blend_residual", "ramp_d"]): s = s.sort_values("duration_d") out[(lf, bl, rp)] = {} for cond in ("any", "ref185", "ref195", "ingress"): tag = "" if cond == "any" else f"_{cond}" out[(lf, bl, rp)][cond] = dict(logd=np.log10(s.duration_d.values), all=s[f"T_eff_yr{tag}"].values, bins=[s[f"T_eff_yr{tag}_I_{b}"].values for b in BINS]) return out def with_ramp_for_bright(te, teff, cond): """Copy of te in which the I < 15 bin (giants, whose ingress can last days) is scaled by the measured ratio of effective exposures for a 4-day and an instantaneous ingress (paired runs).""" num, den = teff[(0.5, 0.0, 4.0)][cond]["all"], teff[(0.5, -1.0, 0.0)][cond]["all"] ratio = np.where(den > 0, np.clip(num / np.maximum(den, 1e-9), 0.0, 1.0), 0.0) out = dict(logd=te["logd"], all=te["all"], bins=[b.copy() for b in te["bins"]]) out["bins"][0] = out["bins"][0] * ratio return out def fold(sl, R, M, tes, lensing=True): """Rate per star-yr (f = 1) and, for each effective-exposure table in tes, rate x T_eff per magnitude bin and for the bin-averaged T_eff.""" t, w = sl.durations(R, M, lensing) lt = np.log10(np.maximum(t, 1e-3)) ks = [] for te in tes: ks.append(np.array([np.sum(w * np.interp(lt, te["logd"], b, left=0.0)) for b in te["bins"]] + [np.sum(w * np.interp(lt, te["logd"], te["all"], left=0.0))])) return float(w.sum()), ks def field_kernels(fields, tes, Rs=R_GRID, Ts=T_GRID, n=N_SAMP, n_m=N_MASS, seed0=1000, lensing=True, lens_mass=0.3, beta=1.0, src_pop="model"): """For every field: rate and (rate x T_eff) per magnitude bin for each fixed radius and for each (population, temperature), for every effective-exposure table in the dict tes. Returned as dict name -> dict key -> array (n_fields, 1 + n_bins + 1).""" names = list(tes) tel = [tes[k] for k in names] keys = [("R", R) for R in Rs] + [(p, T) for p in POPS for T in Ts] out = {nm: {k: np.zeros((len(fields), len(BINS) + 2)) for k in keys} for nm in names} for j, f in enumerate(fields.itertuples()): sl = g.Sightline(f.l, f.b, n=n, seed=seed0 + j, beta=beta, src_pop=src_pop) for R in Rs: r, ks = fold(sl, R, lens_mass, tel, lensing) for nm, k in zip(names, ks): out[nm][("R", R)][j] = np.concatenate([[r], k]) for pname, (lo, hi) in POPS.items(): for T in Ts: m, w, Rm = pm.population_radii(T, n=n_m, lo=lo, hi=hi) acc = [np.zeros(len(BINS) + 2) for _ in names] for mi, wi, Ri in zip(m, w, Rm): r, ks = fold(sl, Ri, mi, tel, lensing) for a, k in zip(acc, ks): a += wi * np.concatenate([[r], k]) for nm, a in zip(names, acc): out[nm][(pname, T)][j] = a return out def calibration(): c = json.load(open(os.path.join(DATA, "rate_sanity_checks.json")))["microlensing_calibration"] return {"weighted": c["ratio_gamma_weighted"], "p16": c["ratio_gamma_16_84"][0], "median": c["ratio_gamma_median"]} def expected_ogle2(fields, kern, case, blend_eff, cal, bins=BINS, incl=None): """Expected events for f = 1 in OGLE-II for one efficiency case; dict key -> N.""" c = CASES[case] incl = incl or INCL[case] if c["fg"] == "blue": okv = fields.frac_I18_with_V > 0.85 blue = np.where(okv, fields.frac_I18_bluer_than_clump, fields.frac_I18_bluer_than_clump[okv].median()) fg = 1.0 - blue else: fg = np.ones(len(fields)) out = {} for key, arr in kern.items(): n = 0.0 for b in bins: kb = BINS.index(b) n += np.sum(fields[f"n_I_{b}"].values * blend_eff[b] * incl[b] * fg * arr[:, 1 + kb]) out[key] = float(n * cal[c["norm"]]) return out def bnu_jy(lam_um, T, L, d_pc): nu = C_LIGHT / (lam_um * 1e-6) b = 2 * H_PLANCK * nu ** 3 / C_LIGHT ** 2 / np.expm1(H_PLANCK * nu / (K_B * T)) return float(np.pi * b * (radius_au(L, T) * AU / (d_pc * PC)) ** 2 / 1e-26) # ------------------------------------------------------------------------------------------------ def main(): fields = pd.read_csv(os.path.join(DATA, "ogle2_fields.csv")) summ = json.load(open(os.path.join(DATA, "ogle2_search_summary.json"))) ebs = json.load(open(os.path.join(DATA, "eb_inclusion_summary.json"))) ebt = pd.read_csv(os.path.join(DATA, "eb_inclusion.csv")) ebt = ebt[ebt.depth_bin.str.contains("deeper than 0.75")].set_index("I_bin") counts = fields[[c for c in fields.columns if c.startswith("n_I_")]] ID_FRAC.update(summ["frac_catalogue_objects_with_id_rad_le_cut_by_bin"]) eb_frac = {b: float(ebt.loc[b, "frac_in_catalogue_and_pass_C2"]) for b in BINS} eb_lo = {b: float(ebt.loc[b, "frac_in_catalogue_and_pass_C2"] - 2 * ebt.loc[b, "binomial_err_C2"]) for b in BINS} first_version = {b: EPS_CAT_MICROLENSING * ID_FRAC[b] for b in BINS} INCL["optimistic"] = eb_frac INCL["conservative"] = {b: min(eb_lo[b], first_version[b]) for b in BINS} teff = load_teff() cal = calibration() base, var = teff[(0.5, 0.0, 0.0)], teff[(0.7, 0.5, 0.0)] # "ingress": conservative exposure restricted to events whose ingress falls inside the survey # (star in the reference image at I_ref < 18.5, 4-day ramp loss for I < 15, star visible on the first frame) tes = {"any": base["any"], "ingress": with_ramp_for_bright(base["ingress"], teff, "ingress")} tes_var = {} for case, c in CASES.items(): te = base[c["teff"]] tes[case] = with_ramp_for_bright(te, teff, c["teff"]) if c["ramp_bright"] else te tes_var[case] = var[c["teff"]] ul0 = poisson_ul(0) # ---- blending efficiency brow, blend = [], {} for case, c in CASES.items(): for tag, xm in (("baseline", X_MIN_BASE), ("variant", X_MIN_VARIANT)): b = blend_mc(counts, c["blend"]["r_blend"], c["blend"]["i_to"], xm) blend[(case, tag)] = b brow.append({"case": case, "search": tag, "x_min": xm, "blend_radius_arcsec": c["blend"]["r_blend"], "turnoff_I": c["blend"]["i_to"], **b}) pd.DataFrame(brow).to_csv(os.path.join(DATA, "blend_efficiency.csv"), index=False, float_format="%.4g") print(pd.DataFrame(brow).to_string(index=False, float_format=lambda x: f"{x:.3g}")) # ---- expected events and limits, OGLE-II kern = field_kernels(fields, tes) kern_var = field_kernels(fields, tes_var, Rs=[3, 10, 100], Ts=[30, 50, 100], seed0=1000) exp = {case: expected_ogle2(fields, kern[case], case, blend[(case, "baseline")], cal) for case in CASES} # alternative exposure definition: the conservative case restricted to events whose ingress falls inside the survey exp_ing = expected_ogle2(fields, kern["ingress"], "conservative", blend[("conservative", "baseline")], cal) rows = [] nstar = fields.n_I_lt_18.sum() kany = kern["any"] for key in kany: raw_rate = float(np.sum(fields.n_I_lt_18.values * kany[key][:, 0]) / nstar) # mean rate per star-yr, model raw_exposure = float(np.sum(fields.n_I_lt_18.values * kany[key][:, -1])) # stars x rate x T_eff(cuts only), no efficiencies rows.append({"kind": "radius" if key[0] == "R" else key[0], "R_au_or_T_K": key[1], "mean_rate_per_star_yr_f1_model": raw_rate, "mean_T_eff_yr": raw_exposure / (raw_rate * nstar) if raw_rate > 0 else np.nan, "N_f1_no_efficiency_factors": raw_exposure, **{f"N_f1_{case}": exp[case][key] for case in CASES}, **{f"f95_{case}": min(1.0, ul0 / exp[case][key]) if exp[case][key] > 0 else 1.0 for case in CASES}, "N_f1_conservative_ingress_only": exp_ing[key], "f95_conservative_ingress_only": min(1.0, ul0 / exp_ing[key]) if exp_ing[key] > 0 else 1.0}) ev = pd.DataFrame(rows) ev.to_csv(os.path.join(DATA, "ogle2_expected_events.csv"), index=False, float_format="%.4g") limR = ev[ev.kind == "radius"].rename(columns={"R_au_or_T_K": "R_au"}) limT = ev[ev.kind != "radius"].rename(columns={"R_au_or_T_K": "T_K", "kind": "population"}) pd.set_option("display.width", 250) print(ev.to_string(index=False, float_format=lambda x: f"{x:.3g}")) # ---- Drake & Cook (2003) conditional limit: their own detected-event number for the KBS disc # clouds (1e5 toward the bulge) scaled by the model ratio of enclosure to cloud event rates kbs = json.load(open(os.path.join(DATA, "rate_sanity_checks.json")))["kbs_drake_cook"] ratios = np.array([v["ratio_enclosure_R7_to_kbs"] for v in kbs.values()]) occ = pd.read_csv(os.path.join(DATA, "occultation_rates.csv")) ob = occ[occ.sightline == "bulge l=1 b=-3"].set_index("R_au").rate_per_star_yr_f1 n_obs_dc = 83 ul_dc = poisson_ul(n_obs_dc) dc = [] for R in (3, 5, 10, 20, 50): scale = ob[R] / np.interp(7.0, ob.index.values, ob.values) nlo, nhi = 1e5 * ratios.min() * scale, 1e5 * ratios.max() * scale dc.append({"R_au": R, "N_f1_low": nlo, "N_f1_high": nhi, "f95_if_all_83_are_candidates_weak": min(1, ul_dc / nlo), "f95_if_all_83_are_candidates_strong": min(1, ul_dc / nhi)}) n7 = 1e5 * ratios dc7 = {"R_au": 7, "N_f1_range": [float(n7.min()), float(n7.max())], "n_obs_upper": n_obs_dc, "poisson_ul_83": float(ul_dc), "f95_range_all_83_counted": [float(ul_dc / n7.max()), float(min(1, ul_dc / n7.min()))], "status": "order of magnitude only: rests on the cloud-model normalisation of Drake & Cook (implied efficiencies " "0.2-1.5 in rate_sanity_checks.json); their light curves were not inspected", "named_candidate": "MACHO 120.21786.958: drop > 3.5 mag, passed their dust-reddening test, i.e. chromatic"} dc = pd.DataFrame(dc) limR = limR.merge(dc, on="R_au", how="left") limR.to_csv(os.path.join(DATA, "limits_f_vs_R.csv"), index=False, float_format="%.4g") limT.to_csv(os.path.join(DATA, "limits_f_vs_T.csv"), index=False, float_format="%.4g") # ---- robustness rob = [] def add(name, key, n_cons, n_opt, n_obs=0, note=""): ul = poisson_ul(n_obs) rob.append({"check": name, "benchmark": f"{key[0]} {key[1]}", "n_obs": n_obs, "N_f1_conservative": n_cons, "N_f1_optimistic": n_opt, "f95_conservative": min(1, ul / n_cons), "f95_optimistic": min(1, ul / n_opt), "note": note}) bench = [("R", 10), ("R", 100), ("R", 1000), ("all MS stars 0.08-1 Msun", 50), ("Sun-like 0.7-1 Msun", 30)] for key in bench: add("baseline", key, exp["conservative"][key], exp["optimistic"][key]) bright = BINS[:3] e17 = {case: expected_ogle2(fields, kern[case], case, blend[(case, "baseline")], cal, bins=bright) for case in CASES} for key in bench: add("stars with I < 17 only", key, e17["conservative"][key], e17["optimistic"][key], note=f"{int(summ['n_stars_I_lt_17'])} stars; blending matters least") ev_var = {case: expected_ogle2(fields, kern_var[case], case, blend[(case, "variant")], cal) for case in CASES} nvar = int(summ["variant_n_pass_all"]) for key in [k for k in bench if k in kern_var["conservative"]]: add("looser variant (drop >= 30%), all survivors counted as real", key, ev_var["conservative"][key], ev_var["optimistic"][key], n_obs=nvar, note=f"{nvar} leftovers taken as events although every one has a mundane explanation") # effect of the exposure definition: same efficiency factors, other effective-exposure tables for nm, lab, note in (("any", "no reference-image condition (as in the first version)", "every event recovered by the cuts is counted"), ("ingress", "only events whose ingress falls inside the survey", "events already in progress on the first frame are not counted (referee's alternative); reference-image condition I_ref < 18.5 and 4-day ramp for I < 15 applied in both columns")): alt = {case: expected_ogle2(fields, kern[nm], case, blend[(case, "baseline")], cal) for case in CASES} for key in bench: add("exposure: " + lab, key, alt["conservative"][key], alt["optimistic"][key], note=note) # first-version inclusion factors (microlensing recovery x centroid fraction of all catalogue objects) fv_c = expected_ogle2(fields, kern["conservative"], "conservative", blend[("conservative", "baseline")], cal, incl=first_version) fv_o = expected_ogle2(fields, kern["optimistic"], "optimistic", blend[("optimistic", "baseline")], cal, incl={b: (208 / 214) * ID_FRAC[b] for b in BINS}) for key in bench: add("inclusion factors of the first version (inferred from microlensing)", key, fv_c[key], fv_o[key], note="0.81 or 0.97 times the centroid-cut fraction of all catalogue objects") # model variations evaluated on a subset of fields (every 4th) and applied as ratios sub = fields.iloc[::4].reset_index(drop=True) tsub = {"conservative": tes["conservative"]} ref = field_kernels(sub, tsub, Rs=[3, 10, 100], Ts=[50], seed0=2000)["conservative"] variants = {"light bending ignored": dict(lensing=False), "lens mass 1.0 Msun instead of 0.3": dict(lens_mass=1.0), "sources weighted by D^2 (beta = 0: giants seen through the whole bulge)": dict(beta=0.0), "bar sources only (no disc sources)": dict(src_pop="bar"), "independent Monte Carlo seed": dict(seed0=9000)} wsub = sub.n_I_lt_18.values for name, kw in variants.items(): kw2 = dict(seed0=2000) kw2.update(kw) alt = field_kernels(sub, tsub, Rs=[3, 10, 100], Ts=[50], **kw2)["conservative"] for key in (("R", 3), ("R", 10), ("R", 100), ("all MS stars 0.08-1 Msun", 50)): ratio = float(np.sum(wsub * alt[key][:, -1]) / np.sum(wsub * ref[key][:, -1])) rob.append({"check": "model: " + name, "benchmark": f"{key[0]} {key[1]}", "n_obs": 0, "N_f1_conservative": exp["conservative"].get(key, np.nan) * ratio if key in exp["conservative"] else np.nan, "N_f1_optimistic": np.nan, "f95_conservative": np.nan, "f95_optimistic": np.nan, "note": f"expected events change by a factor {ratio:.2f}"}) # finite ingress for every source, and residual blend light: ratios of effective exposure (paired injection runs) refrun = teff[(0.5, -1.0, 0.0)]["ref185"] for keyrun, lab in (((0.5, 0.0, 4.0), "4-day ingress and egress for every source"), ((0.5, 0.0, 10.0), "10-day ingress and egress for every source"), ((0.5, 0.3, 0.0), "30% of the catalogued flux remains during the event"), ((0.5, 0.4, 0.0), "40% of the catalogued flux remains during the event")): alt_te = teff[keyrun]["ref185"] kk = field_kernels(sub, {"ref": refrun, "alt": alt_te}, Rs=[3, 10, 100], Ts=[50], seed0=2000) for key in (("R", 3), ("R", 10), ("R", 100), ("all MS stars 0.08-1 Msun", 50)): ratio = float(np.sum(wsub * kk["alt"][key][:, -1]) / np.sum(wsub * kk["ref"][key][:, -1])) rob.append({"check": "injection: " + lab, "benchmark": f"{key[0]} {key[1]}", "n_obs": 0, "N_f1_conservative": exp["conservative"][key] * ratio, "N_f1_optimistic": np.nan, "f95_conservative": np.nan, "f95_optimistic": np.nan, "note": f"expected events change by a factor {ratio:.2f}"}) rob = pd.DataFrame(rob) rob.to_csv(os.path.join(DATA, "robustness.csv"), index=False, float_format="%.4g") print(rob.to_string(index=False, float_format=lambda x: f"{x:.3g}")) # ---- forecasts for other existing data sets ----------------------------------------------------- popr = pd.read_csv(os.path.join(DATA, "occultation_rates_population.csv")) def pop_ratio(sight, pop, T): # rate(population at T) / rate(R = 10 au) on a reference sightline a = popr[(popr.sightline == sight) & (popr.population == pop) & (popr.T_K == T)].rate_per_star_yr_f1.iloc[0] return a / occ[(occ.sightline == sight) & (occ.R_au == 10)].rate_per_star_yr_f1.iloc[0] def rate(sight, R=10): return float(occ[(occ.sightline == sight) & (occ.R_au == R)].rate_per_star_yr_f1.iloc[0]) fld = "field star at 1 kpc, l=90 b=30" fr50, fr30s = pop_ratio(fld, "all MS stars 0.08-1 Msun", 50), pop_ratio(fld, "Sun-like 0.7-1 Msun", 30) br50 = pop_ratio("bulge l=1 b=-3", "all MS stars 0.08-1 Msun", 50) br30s = pop_ratio("bulge l=1 b=-3", "Sun-like 0.7-1 Msun", 30) fc = [] def add_fc(name, stars, years, n10, n50, n30s, records, searched, basis): fc.append({"data_set": name, "stars": stars, "effective_years": years, "N_f1_R10au": n10, "N_f1_allMS_50K": n50, "N_f1_sunlike_30K": n30s, "pipeline_records_vanished_star": records, "searched_for_total_occultations": searched, "basis": basis}) k10, k50, k30s = ("R", 10), ("all MS stars 0.08-1 Msun", 50), ("Sun-like 0.7-1 Msun", 30) gm = lambda key: float(np.sqrt(exp["conservative"][key] * exp["optimistic"][key])) add_fc("OGLE-II bulge DIA 1997-2000, I<18", float(nstar), float(ev[(ev.kind == "radius") & (ev.R_au_or_T_K == 10)].mean_T_eff_yr.iloc[0]), gm(k10), gm(k50), gm(k30s), "yes", "yes (this work): 0 found", "this work; geometric mean of conservative and optimistic efficiencies") # OGLE-IV: per-field star counts (I<18) from Mroz et al. 2019 table 5, same seasonal duty cycle as OGLE-II t5 = pd.read_csv(os.path.join(RAW, "mroz2019_table5.csv")) f4 = t5.groupby("Field").agg(l=("GLON", "mean"), b=("GLAT", "mean"), n18=("N18", "sum"), n21=("N21", "sum")).reset_index() k4 = field_kernels(f4, {"any": tes["any"]}, Rs=[10, 100], Ts=[30, 50], n=12000, n_m=12, seed0=4000)["any"] eff_mid = float(np.sqrt(exp["conservative"][k10] * exp["optimistic"][k10]) / ev[(ev.kind == "radius") & (ev.R_au_or_T_K == 10)].N_f1_no_efficiency_factors.iloc[0]) eff_cons = float(exp["conservative"][k10] / ev[(ev.kind == "radius") & (ev.R_au_or_T_K == 10)].N_f1_no_efficiency_factors.iloc[0]) o4 = {} for seasons, tag in ((8, "2010-2017"), (12.5, "2010-2024")): sc = seasons / N_SEASONS_OGLE2 vals = {key: float(np.sum(f4.n18.values * k4[key][:, -1]) * sc) for key in (k10, ("R", 100), k50, k30s)} o4[tag] = vals add_fc(f"OGLE-IV bulge {tag}, I<18", float(f4.n18.sum()), 0.5 * seasons, vals[k10] * eff_mid, vals[k50] * eff_mid, vals[k30s] * eff_mid, "yes", "no", "Mroz et al. 2019 star counts; OGLE-II duty cycle and efficiency factors assumed") v21 = float(np.sum(f4.n21.values * k4[k10][:, -1]) * 8 / N_SEASONS_OGLE2) add_fc("OGLE-IV bulge 2010-2017, I<21 (upper bound)", float(f4.n21.sum()), 4.0, v21, v21 * br50, v21 * br30s, "yes, at low S/N", "no", "all sources; no efficiency applied: most are too faint or blended") add_fc("MACHO bulge 1993-1999 (Drake & Cook 2003)", 3.32e7, 1.71e8 / 3.32e7, float(1e5 * np.sqrt(ratios.min() * ratios.max()) * ob[10] / np.interp(7.0, ob.index.values, ob.values)), np.nan, np.nan, "probably", "partly: dust-like dips; <= 83 unvetted", "their detected-event number for KBS clouds x model rate ratio") add_fc("MACHO LMC 1992-1999", 1.29e7, 9.38e7 / 1.29e7, 9.38e7 * rate("LMC (MW disc + self)"), 9.38e7 * rate("LMC (MW disc + self)") * fr50, 9.38e7 * rate("LMC (MW disc + self)") * fr30s, "probably", "partly (same paper)", "exposure x model rate, no efficiency") nv = 1e8 rv = rate("VVV-WIT-08 position") add_fc("VVV bulge 2010-2019, Ks < 15-16", nv, 5.0, nv * 5.0 * rv, nv * 5.0 * rv * br50, nv * 5.0 * rv * br30s, "no (detection-based)", "no (one found by luck)", "Smith et al. 2021: 1.5e6 sources/deg2 x 200 deg2, a third usable, 10 yr at 50%") add_fc("Gaia DR3 sources with G < 17.65", 2.2e8, 5.5, 2.2e8 * 5.5 * rate("disc star at 3 kpc, l=30 b=2"), 2.2e8 * 5.5 * rate("disc star at 3 kpc, l=30 b=2") * fr50, 2.2e8 * 5.5 * rate("disc star at 3 kpc, l=30 b=2") * fr30s, "no (missing transits only)", "no", "all taken as disc stars at 3 kpc (generous); 66 months of DR4") add_fc("ZTF FGK dwarfs (Tzanidakis et al. 2025)", 6.3e7, 5.0, 6.3e7 * 5 * rate("field star at 1 kpc, l=90 b=10"), 6.3e7 * 5 * rate("field star at 1 kpc, l=90 b=10") * fr50, 6.3e7 * 5 * rate("field star at 1 kpc, l=90 b=10") * fr30s, "no (gap)", "dips only", "stars at 1 kpc, b = 10 deg") add_fc("ASAS-SN 13= 0.3 mag: none total", "stars at 1 kpc, b = 10 deg") add_fc("Kepler prime mission", 2e5, 4.0, 2e5 * 4 * rate("Kepler field star at 1 kpc"), 2e5 * 4 * rate("Kepler field star at 1 kpc") * fr50, 2e5 * 4 * rate("Kepler field star at 1 kpc") * fr30s, "yes", "yes (none)", "stars at 1 kpc") fc = pd.DataFrame(fc) for col in ("N_f1_R10au", "N_f1_allMS_50K", "N_f1_sunlike_30K"): for fv in (1e-3, 1e-4, 1e-5): fc[f"{col.replace('N_f1_', 'N_')}_f{fv:g}"] = fc[col] * fv fc.to_csv(os.path.join(DATA, "forecast_events.csv"), index=False, float_format="%.4g") print(fc[["data_set", "stars", "effective_years", "N_f1_R10au", "N_f1_allMS_50K", "N_f1_sunlike_30K"]].to_string( index=False, float_format=lambda x: f"{x:.3g}")) # ---- static "vanished star" test (two-epoch catalogue comparisons) R172 = float(radius_au(1.0, 30.0)) r100 = occ[(occ.sightline == fld) & (occ.R_au == 100)].iloc[0] vas = {"R_au": R172, "rate_per_star_yr_f1": float(r100.rate_per_star_yr_f1 * R172 / 100), "median_duration_yr": float(r100.dur_d_50 * R172 / 100 / 365.25), "prob_star_vanished_between_two_epochs_60yr_apart_f1": float(r100.rate_per_star_yr_f1 * R172 / 100 * min(60.0, r100.dur_d_50 * R172 / 100 / 365.25)), "note": "per star at 1 kpc, b = 30 deg; solar-luminosity 30 K enclosures only"} # ---- follow-up: thermal flux of the occulter itself at the lens distance fol = {f"L={L:g} Lsun, T={T:g} K, d={d:g} pc": {"F_25um_mJy": bnu_jy(25.5, T, L, d) * 1e3, "F_100um_mJy": bnu_jy(100.0, T, L, d) * 1e3} for L, T, d in ((1.0, 100.0, 4000.0), (1.0, 50.0, 4000.0), (0.01, 100.0, 4000.0), (1.0, 30.0, 4000.0))} # ---- lensing shoulder: magnification of the major image at the moment of disappearance sh = {} for rho in (1.2, 1.5, 2.0, 3.0, 5.0): u = rho - 1 / rho A = (u * u + 2) / (u * np.sqrt(u * u + 4)) sh[f"R/r_E={rho:g}"] = float(0.5 * (A + 1)) hid = pd.read_csv(os.path.join(DATA, "wise_hidden_temperature.csv")) getR = lambda R, col: float(limR[limR.R_au == R][col].iloc[0]) getT = lambda p, T, col: float(limT[(limT.population == p) & (limT.T_K == T)][col].iloc[0]) summary = { "ogle2": {"n_stars_I_lt_18": int(nstar), "n_candidates": summ["n_candidates"], "poisson_ul_0": float(ul0), "T_eff_yr_by_duration": {cond: dict(zip([f"{10 ** x:.0f} d" for x in base[cond]["logd"]], [float(v) for v in base[cond]["all"]])) for cond in ("any", "ref195", "ref185", "ingress")}, "efficiency_factors": {case: {"inclusion_and_centroid_by_bin": INCL[case], "norm": cal[CASES[case]["norm"]], "eps_blend_all_I_lt_18": blend[(case, "baseline")]["all_I_lt_18"], "exposure_table": CASES[case]["teff"], "ramp_4d_applied_to_I_lt_15": CASES[case]["ramp_bright"], "eps_fg": (1 - summ["frac_I18_bluer_than_clump_median_fields_with_V"]) if CASES[case]["fg"] == "blue" else 1.0} for case in CASES}, "inclusion_measured_on_deep_eclipsing_binaries_by_bin": eb_frac, "inclusion_first_version_conservative_by_bin": first_version, "eclipsing_binary_test": {k: ebs[k] for k in ("n_deep_I_lt_18", "frac_deep_I_lt_18_in_catalogue", "frac_deep_I_lt_18_in_catalogue_and_pass_C2")}, "overall_efficiency_R10": {case: float(exp[case][k10] / ev[(ev.kind == "radius") & (ev.R_au_or_T_K == 10)].N_f1_no_efficiency_factors.iloc[0]) for case in CASES}, "systematic_note": "Limits should be read to one significant figure. Beyond the conservative-optimistic bracket, " "which is dominated by the modelled blending correction, an overall systematic of about a factor " "of 2 is attached: catalogue inclusion is measured on eclipsing binaries and by injection into " "catalogue light curves, not by injection into images."}, "f95_vs_R": {f"R={R} au": {c: getR(R, f"f95_{c}") for c in list(CASES) + ["conservative_ingress_only"]} for R in R_GRID}, "N_f1_vs_R": {f"R={R} au": {c: getR(R, f"N_f1_{c}") for c in CASES} for R in R_GRID}, "f95_vs_T": {p: {f"T={T} K": {c: getT(p, T, f"f95_{c}") for c in list(CASES) + ["conservative_ingress_only"]} for T in T_GRID} for p in POPS}, "drake_cook_2003_conditional": dc7, "ogle4_forecast_f95_conservative_efficiency": {tag: {f"{k[0]} {k[1]}": float(ul0 / (v * eff_cons)) for k, v in vals.items()} for tag, vals in o4.items()}, "ogle4_expected_events_f1_no_efficiency": {tag: {f"{k[0]} {k[1]}": v for k, v in vals.items()} for tag, vals in o4.items()}, "static_vanished_star_test": vas, "follow_up_thermal_flux": fol, "lensing_shoulder_magnification_at_disappearance": sh, "comparison": { "gravity_oort_limit_fraction_of_local_stellar_mass": 0.42, "T_hidden_from_AllWISE_K": {f"L={r.L_Lsun:g} d={r.d_pc:g} pc": (float(r.T_hidden_K) if np.isfinite(r.T_hidden_K) else None) for r in hid.itertuples()}, }, } with open(os.path.join(DATA, "limits_summary.json"), "w") as f: json.dump(summary, f, indent=1) print(json.dumps({k: summary[k] for k in ("ogle2", "drake_cook_2003_conditional", "ogle4_forecast_f95_conservative_efficiency", "static_vanished_star_test", "follow_up_thermal_flux", "lensing_shoulder_magnification_at_disappearance")}, indent=1)) # ---- figure 4 ---------------------------------------------------------------------------------- from matplotlib.ticker import FixedLocator, FixedFormatter fig, ax = plt.subplots(1, 3, figsize=(19, 5.8)) a = ax[0] a.fill_between(limR.R_au, limR.f95_optimistic, limR.f95_conservative, color="#d62728", alpha=0.35, label="OGLE-II, this work: 0 events\n(optimistic to conservative efficiency)") a.loglog(limR.R_au, limR.f95_conservative, color="#d62728", lw=1.8) a.loglog(limR.R_au, limR.f95_conservative_ingress_only, color="#d62728", lw=1.0, ls=":", label="conservative, counting only events whose\ningress falls inside the survey") d2 = limR.dropna(subset=["N_f1_low"]) a.fill_between(d2.R_au, d2.f95_if_all_83_are_candidates_strong, d2.f95_if_all_83_are_candidates_weak, color="#1f77b4", alpha=0.3, label="MACHO (Drake & Cook 2003), conditional:\nall 83 unvetted leftovers counted") for tag, ls in (("2010-2017", "--"), ("2010-2024", ":")): v = o4[tag] a.loglog([10, 100], [ul0 / (v[("R", 10)] * eff_cons), ul0 / (v[("R", 100)] * eff_cons)], ls, color="k", marker="s", ms=4, label=f"reach of the same search in OGLE-IV {tag}") a.axhline(0.42, color="#2ca02c", lw=1.2) a.text(900, 0.5, "gravity (local mass budget): f < 0.42", color="#2ca02c", fontsize=8, ha="right") a.axvspan(0.8, 2.0, color="0.85", zorder=0) a.text(0.86, 1.3e-6, "R < Einstein\nradius: lensing,\nno occultation", fontsize=7.5, color="0.3") a.set_xlim(0.8, 1000) a.set_ylim(1e-6, 2) a.set_xlabel("enclosure radius R [au]") a.set_ylabel("95% upper limit on enclosed fraction of stars, f") a.set_title("(a) Limit on f for opaque spheres of fixed radius", fontsize=10) a.legend(fontsize=7.5, loc="lower left", bbox_to_anchor=(0.2, 0.0)) a.grid(alpha=0.3, which="both") top = a.twiny() top.set_xscale("log") top.set_xlim(a.get_xlim()) tt = [300, 200, 100, 50, 30, 20, 15] top.xaxis.set_major_locator(FixedLocator([(393.6 / t) ** 2 for t in tt])) top.xaxis.set_major_formatter(FixedFormatter([str(t) for t in tt])) top.xaxis.set_minor_locator(FixedLocator([])) top.set_xlabel("blackbody temperature of a sphere of this radius for L = 1 L$_\\odot$ [K]", fontsize=8) a = ax[1] for p, c in zip(POPS, ("#d62728", "#9467bd")): s = limT[limT.population == p] a.fill_between(s.T_K, s.f95_optimistic, s.f95_conservative, color=c, alpha=0.3) a.loglog(s.T_K, s.f95_conservative, color=c, lw=1.8, label=("every main-sequence star equally likely" if p.startswith("all") else "Sun-like stars only (f = fraction of them)")) a.axvspan(100, 700, color="#ff7f0e", alpha=0.12, zorder=0) a.text(104, 1.3e-5, "mid-IR excess searches\n(Suazo et al. 2022, 2024:\n100-700 K, d < 300 pc)", fontsize=7.5, color="#b35900", va="bottom") for L, lab in ((1.0, "1 L$_\\odot$"), (0.01, "0.01 L$_\\odot$")): th = float(hid[(hid.L_Lsun == L) & (hid.d_pc == 100)].T_hidden_K.iloc[0]) a.axvline(th, color="0.4", ls=":", lw=1) a.text(th * 1.03, 1.5e-5, f"colder: hidden from AllWISE at 100 pc ({lab})", fontsize=6.8, color="0.3", rotation=90, va="bottom") a.axhline(0.42, color="#2ca02c", lw=1.2) a.text(21, 0.5, "gravity: f < 0.42", color="#2ca02c", fontsize=8) a.set_xlim(20, 320) a.set_ylim(1e-5, 2) a.xaxis.set_major_locator(FixedLocator([20, 30, 50, 75, 100, 150, 200, 300])) a.xaxis.set_major_formatter(FixedFormatter(["20", "30", "50", "75", "100", "150", "200", "300"])) a.xaxis.set_minor_locator(FixedLocator([])) a.set_xlabel("enclosure temperature T [K]") a.set_ylabel("95% upper limit on f") a.set_title("(b) OGLE-II limit when every enclosure has R = sqrt(L) (394 K / T)$^2$ au", fontsize=10) a.legend(fontsize=7.5, loc="upper left", bbox_to_anchor=(0.0, 0.88)) a.grid(alpha=0.3, which="both") a = ax[2] f2 = fc.copy() y = np.arange(len(f2))[::-1] col = ["#2ca02c" if r.startswith("yes (this") else ("#1f77b4" if q.startswith("yes") or q.startswith("prob") else "0.6") for r, q in zip(f2.searched_for_total_occultations, f2.pipeline_records_vanished_star)] a.barh(y, f2.N_f1_R10au * 1e-3, color=col, height=0.36, align="edge") a.barh(y - 0.36, f2.N_f1_allMS_50K.fillna(0) * 1e-3, color=col, height=0.36, align="edge", alpha=0.45) a.set_xscale("log") a.set_yticks(y) a.set_yticklabels([d.replace(" (", "\n(") for d in f2.data_set], fontsize=7.2) a.axvline(1.0, color="k", lw=0.8) a.axvline(ul0, color="k", lw=0.8, ls="--") a.set_xlim(1e-5, 1e3) a.set_xlabel("expected total occultations if f = 10$^{-3}$\n(solid: R = 10 au; pale: all stars at 50 K)", fontsize=9) a.set_title("(c) Where the events would be\ngreen: searched here; blue: pipeline records a vanished star; grey: it does not", fontsize=9) a.grid(alpha=0.3, which="both", axis="x") fig.tight_layout() fig.savefig(os.path.join(FIGS, "fig4_limits.png"), dpi=130) plt.close(fig) if __name__ == "__main__": main()