"""Handle 3: microlensing. (i) Geometry: a lens surrounded by an opaque sphere of radius R_s. With rho = R_s / r_E: rho < 1 : the major image (outside the Einstein ring) is never covered; the minor image is covered when u > 1/rho - rho. For rho < 0.618 the event is normal inside u < 1. rho > 1 : both images are covered for u < u_c = rho - 1/rho -> the source disappears. There is no microlensing event; there is a total occultation with weak lensing "shoulders" (magnification of the major image at u = u_c). (ii) Optical depth and event rate (OGLE-IV, MOA-II): what a dark lens population would do. (iii) Lenses with measured masses and high-resolution imaging: luminous versus dark. (iv) Forecast of the occultation rate by cold enclosures, normalised in two independent ways (the measured microlensing event rate, and the Galaxia-based rates of Smith et al. 2021). The event table is compiled by hand from the abstracts of the cited papers (archived by 01_download.py in RAW/arxiv/). Revision after the referee report: * the lens-light fraction is INDICATIVE, not a limit on stars in general: the events are a hand-compiled, non-representative set of published follow-ups (mostly planet hosts). The conservative counts are now the lead numbers, and the tension with the ~20% remnant fraction expected among bulge lenses (Gould 2000) is computed. * the optical-depth "ceiling" of the first pass is withdrawn: it was the MOA-II/OGLE-IV ratio times 1.1, not a measurement of room above a model. The ratio is kept only as a remark. * the occultation forecast is always quoted as the range between its two normalisations and before any detection efficiency; its survey inputs are marked as assumed or sourced. The rate calculation for dark occulters at stellar density is due to Smith et al. (2021). """ import json import os import sys import numpy as np import pandas as pd from scipy.stats import beta, binom 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, SIGMA_SB, L_SUN, AU, KAPPA_MAS T1AU = (L_SUN / (4 * np.pi * SIGMA_SB * AU ** 2)) ** 0.25 # 393.6 K # event, M1, M2, D_L [kpc], status, arXiv, basis # status: L = lens light detected and consistent with a main-sequence lens of the stated mass # D = main-sequence lens excluded by a non-detection (dark lens) # A = ambiguous (no lens light attributed, mass not independently fixed; a faint # M dwarf is allowed) -> counted as dark only in the conservative variant # N = listed for completeness, not counted (no high-resolution imaging, or selected # for a long time-scale, i.e. biased towards massive dark lenses) # 'luminous' entries whose identification is weaker than a resolved lens (referee): unresolved # excess flux, a 'candidate' host, or a disputed identification WEAK_LUMINOUS = ["KMT-2018-BLG-0029", "OGLE-2019-BLG-0960", "MOA-2011-BLG-262", "OGLE-2016-BLG-1195"] EVENTS = [ ("OGLE-2005-BLG-169", 0.69, np.nan, 4.1, "L", "1507.08661; 1507.08914", "HST + Keck, lens resolved"), ("OGLE-2005-BLG-071", 0.426, np.nan, np.nan, "L", "1909.04740", "Keck"), ("OGLE-2006-BLG-109", 0.51, np.nan, np.nan, "L", "0911.2706", "mass from parallax + finite source; Keck flux agrees"), ("MOA-2007-BLG-192", 0.28, np.nan, 2.16, "L", "2403.12118", "Keck + HST"), ("MOA-2007-BLG-400", 0.69, np.nan, 6.89, "L", "2009.02329", "Keck"), ("MOA-2008-BLG-379", 0.434, np.nan, 3.44, "L", "2311.00627", "Keck + HST"), ("MOA-2009-BLG-319", 0.514, np.nan, 7.0, "L", "2009.08461", "Keck"), ("OGLE-2011-BLG-0950", 1.12, 0.47, 6.70, "L", "2206.03502", "Keck"), ("MOA-2011-BLG-293", 0.86, np.nan, 7.72, "L", "1310.3706", "Keck"), ("MOA-2011-BLG-262", 0.19, np.nan, 7.49, "L", "2410.09147", "Keck (candidate lens host, K = 22.3)"), ("OGLE-2012-BLG-0026", 1.06, np.nan, 4.0, "L", "1601.01846", "Keck"), ("OGLE-2012-BLG-0563", 0.801, np.nan, 5.46, "L", "2412.03651", "HST + Keck"), ("OGLE-2012-BLG-0950", 0.58, np.nan, 2.19, "L", "1809.02654", "HST + Keck"), ("MOA-2013-BLG-220", 0.88, np.nan, 6.72, "L", "1909.04444", "Keck"), ("OGLE-2013-BLG-0132", 0.495, np.nan, 3.48, "L", "2401.17549", "Keck"), ("OGLE-2014-BLG-0124", 0.89, np.nan, 3.6, "L", "1709.00806", "Spitzer parallax + Keck"), ("OGLE-2014-BLG-0676", 0.60, np.nan, 1.88, "L", "2607.18408", "Keck"), ("OGLE-2014-BLG-1050", 1.05, 0.38, 3.43, "L", "2112.08029", "Spitzer parallax; Magellan AO flux agrees"), ("OGLE-2016-BLG-1195", 0.62, np.nan, 7.45, "L", "2302.01168; 2303.08876", "Keck; identification disputed by Gould 2023"), ("OGLE-2017-BLG-1434", 0.234, np.nan, 0.86, "L", "2107.06510", "Keck"), ("Kojima-1 (TCP J0507+2447)", 0.581, np.nan, 0.505, "L", "1909.11802; 1912.00038", "spectroscopy + AO; VLTI theta_E"), ("MOA-2010-BLG-328", 0.35, np.nan, 3.5, "L", "2504.06347", "Keck + HST: 0.2 Msun at 2-3 kpc or 0.5 Msun at 4-5 kpc"), ("KMT-2018-BLG-0029", 0.70, np.nan, np.nan, "L", "2508.18343", "Keck excess consistent with a main-sequence host"), ("OGLE-2019-BLG-0960", 0.40, np.nan, np.nan, "L", "2508.18343", "Keck excess consistent with a main-sequence host"), ("OGLE-2007-BLG-349", 0.41, 0.30, np.nan, "L", "1609.06720", "HST; too faint for one 0.7 Msun star, fits two M dwarfs"), ("OGLE-2006-BLG-284", 0.43, 0.12, np.nan, "L", "2610.07133", "Keck + HST, lens resolved"), ("MOA-2010-BLG-477", 0.53, np.nan, 2.0, "D", "2110.07934", "Keck non-detection; white dwarf"), ("KMT-2020-BLG-0414", np.nan, np.nan, np.nan, "D", "2409.02157", "Keck non-detection; white dwarf"), ("MOA-2008-BLG-310", 0.21, np.nan, np.nan, "A", "1703.06947", "excess flux is not the lens; lens I > 22.44; mass is a Bayesian estimate"), ("MOA-2016-BLG-227", np.nan, np.nan, np.nan, "A", "1704.01724", "excess flux probably not mainly the lens"), ("MACHO-97-BLG-28", 0.24, 0.07, 7.0, "A", "1905.08540", "no significant excess; faint M dwarf allowed"), ("OGLE-2015-BLG-1285", 1.35, np.nan, np.nan, "A", "1508.06636", "M1 > 1.35 Msun at 80%; main sequence excluded by blue-light limit (no AO)"), ("Gaia16aye", 0.57, 0.36, 0.78, "N", "1901.07281", "two main-sequence stars; no high-resolution imaging quoted"), ("Gaia19bld", 1.13, np.nan, 5.52, "N", "2112.01613", "likely main sequence, 'yet to be verified'"), ("OGLE-2014-BLG-0221", np.nan, np.nan, np.nan, "N", "2310.13066", "late-type star or remnant; no AO yet"), ("OGLE-2011-BLG-0462", 7.1, np.nan, 1.58, "N", "2201.13296; 2308.03302", "black hole (long time-scale selection)"), ] def cp_upper(k, n, cl=0.95): """Clopper-Pearson one-sided upper limit on a binomial fraction.""" return 1.0 if k >= n else float(beta.ppf(cl, k + 1, n - k)) def r_einstein_au(M, DL, DS=8.0): return np.sqrt(KAPPA_MAS * M * (1.0 / DL - 1.0 / DS)) * DL def r_enclosure_au(L, T): return np.sqrt(L) * (T1AU / T) ** 2 def a_total(u): return (u * u + 2) / (u * np.sqrt(u * u + 4)) def lightcurve_opaque(u, rho): """Flux (source = 1) of a point source lensed by a point mass inside an opaque disc of radius rho (Einstein units): images at theta+- = (sqrt(u^2+4) +- u)/2 are removed if < rho.""" u = np.asarray(u, float) A = a_total(u) ap, am = 0.5 * (A + 1), 0.5 * (A - 1) tp = 0.5 * (np.sqrt(u * u + 4) + u) tm = 0.5 * (np.sqrt(u * u + 4) - u) return np.where(tp > rho, ap, 0.0) + np.where(tm > rho, am, 0.0) def main(): ev = pd.DataFrame(EVENTS, columns=["event", "M1_Msun", "M2_Msun", "DL_kpc", "status", "arxiv", "basis"]) ev.to_csv(os.path.join(DATA, "microlens_lens_light.csv"), index=False) nL, nD, nA = int((ev.status == "L").sum()), int((ev.status == "D").sum()), int((ev.status == "A").sum()) lum = ev[ev.status == "L"] nW = int(ev.event.isin(WEAK_LUMINOUS).sum()) assert nW == len(WEAK_LUMINOUS) p_rem = 0.20 # remnant fraction of bulge lenses estimated by Gould (2000) lens = { "status": "indicative only: hand-compiled, non-representative sample of published follow-ups (mostly planet hosts)", "n_luminous": nL, "n_dark": nD, "n_ambiguous": nA, "n_not_counted": int((ev.status == "N").sum()), "n_luminous_weak_identification": nW, "dark_fraction_observed": nD / (nL + nD), # nominal: ambiguous events dropped "f95_dark_nominal": cp_upper(nD, nL + nD), # conservative (lead number): ambiguous events and the one disputed identification counted as dark "n_dark_conservative": nD + nA + 1, "n_total_conservative": nL + nD + nA, "f95_dark_conservative": cp_upper(nD + nA + 1, nL + nD + nA), # most conservative: every weak 'luminous' identification also counted as dark "n_dark_most_conservative": nD + nA + nW, "n_total_most_conservative": nL + nD + nA, "f95_dark_most_conservative": cp_upper(nD + nA + nW, nL + nD + nA), "luminous_mass_range": [float(lum.M1_Msun.min()), float(lum.M1_Msun.max())], "luminous_DL_range_kpc": [float(lum.DL_kpc.min()), float(lum.DL_kpc.max())], "expected_remnant_fraction_of_bulge_events_Gould2000": p_rem, # probability of seeing this few dark lenses if 20% of lenses were remnants: the sample is # probably biased towards luminous lenses (follow-up is published when lens light is found) "P_le_observed_dark_if_20pct_remnants_nominal": float(binom.cdf(nD, nL + nD, p_rem)), "P_le_observed_dark_if_20pct_remnants_conservative": float(binom.cdf(nD + nA + 1, nL + nD + nA, p_rem)), } # ---- optical depth and event rate ------------------------------------------------------------ tau_ogle, tau_moa = 1.36e-6, 1.75e-6 # at b = -3 deg (Mroz+2019; Nunota+2024) gam_ogle, gam_moa = 13.4e-6, 16.08e-6 # per star per year at b = -3 deg sys_src = 0.10 # OGLE-IV systematic from source counts # REMARK ONLY (no limit is derived). A hidden lens population (fraction f of stellar mass, # R_s < r_E) would raise tau by 1/(1-f) over a model normalised to star counts. No published # data/model ratio with an uncertainty exists: Mroz et al. (2019) say their maps "agree # remarkably well" with the Besancon-based simulator without quoting a ratio; Nunota et al. # (2024) say MOA-II is consistent with OGLE-IV yet "inconsistent with a prediction based on # Galactic models"; and bulge models are partly normalised dynamically, so hidden stellar mass # would already be partly inside them. The survey-to-survey ratio below shows only that the # measurements themselves differ by ~30% at b = -3 deg. tau = {"tau_OGLE_IV_b-3": tau_ogle, "tau_MOA_II_b-3": tau_moa, "rate_OGLE_IV": gam_ogle, "rate_MOA_II": gam_moa, "ratio_tau_MOA_over_OGLE": tau_moa / tau_ogle, "ratio_rate_MOA_over_OGLE": gam_moa / gam_ogle, "OGLE_IV_source_count_systematic": sys_src, "status": "remark only; NOT a limit (the first-pass 'ceiling' of 0.3 is withdrawn)"} # ---- occultation forecast -------------------------------------------------------------------- # Survey inputs. N_SRC and T_YR: Mroz et al. (2019), Sec. 2: 121 fields, "over 160 deg2 ... # over 400 million sources in OGLE databases", 2010-2017 (upper bounds: database objects, and # not every field has 8 seasons). GAMMA is ASSUMED: the measured rate is 13.4e-6 per star per # year at b = -3 deg and falls as exp(-0.49 (|b| - 3)); 1e-5 is an order-of-magnitude mean. N_SRC, T_YR, GAMMA = 4.0e8, 8.0, 1.0e-5 n_true = GAMMA * N_SRC * T_YR V_REL = 150.0 / 4.74 # assumed lens-source relative velocity at the lens, au/yr # Smith+2021: 0.80 events deg^-2 yr^-1 (duration > 10 d, depth > 0.5) per 1e6 sources deg^-2 for # 0.5 au occulters with the space density of stars; rates scale linearly with occulter radius SMITH_PER_SOURCE_YR = 0.80e-6 VVV_SRC_YR = 200.0 * 1.5e6 * 10.0 * 0.5 # deg2 x sources/deg2 x yr x duty cycle rows = [] for (M, L, tag) in ((1.0, 1.0, "G dwarf"), (0.6, 0.08, "K/M dwarf"), (0.3, 0.011, "M dwarf")): for T in (300.0, 150.0, 100.0, 50.0, 30.0, 20.0): rE = r_einstein_au(M, 6.0) Rs = r_enclosure_au(L, T) rho = Rs / rE uc = rho - 1 / rho if rho > 1 else 0.0 dur = 0.5 * np.pi * uc * rE / V_REL * 365.25 if uc > 0 else 0.0 if rho < 0.618: regime = "normal event; lens dark" elif rho < 1: regime = "event with minor image clipped" else: regime = "total occultation, no event" rows.append({"lens": tag, "M_Msun": M, "L_Lsun": L, "T_K": T, "R_s_au": Rs, "r_E_au_DL6kpc": rE, "rho": rho, "regime": regime, "u_c": uc, "mean_duration_d": dur, "major_image_flux_at_edge": float(0.5 * (a_total(uc) + 1)) if uc > 0 else np.nan, "N_occ_OGLEIV_per_unit_f (rate norm.)": uc * n_true, "N_occ_OGLEIV_per_unit_f (Smith norm.)": SMITH_PER_SOURCE_YR * (Rs / 0.5) * N_SRC * T_YR if uc > 0 else 0.0, "N_occ_VVV_per_unit_f (Smith norm.)": SMITH_PER_SOURCE_YR * (Rs / 0.5) * VVV_SRC_YR if uc > 0 else 0.0}) occ = pd.DataFrame(rows) occ.to_csv(os.path.join(DATA, "microlens_occultation_forecast.csv"), index=False, float_format="%.4g") ref = occ[(occ.lens == "M dwarf") & (occ.T_K == 50.0)].iloc[0] forecast = { "N_lensing_events_true_OGLEIV_assumed": n_true, "reference_case": "0.3 Msun, 0.011 Lsun, 50 K, D_L = 6 kpc", "reference_R_s_au": float(ref.R_s_au), "reference_rho": float(ref.rho), "reference_duration_d": float(ref.mean_duration_d), "reference_N_per_unit_f_rate_norm": float(ref["N_occ_OGLEIV_per_unit_f (rate norm.)"]), "reference_N_per_unit_f_smith_norm": float(ref["N_occ_OGLEIV_per_unit_f (Smith norm.)"]), # fraction that a search with 10% efficiency and no detections would exclude at 95% "f95_if_zero_found_eff10pct_rate_norm": 3.0 / (0.1 * float(ref["N_occ_OGLEIV_per_unit_f (rate norm.)"])), "f95_if_zero_found_eff10pct_smith_norm": 3.0 / (0.1 * float(ref["N_occ_OGLEIV_per_unit_f (Smith norm.)"])), "smith2021_expected_VVV_events_0.5au_at_stellar_density": 1.2e-2 * 200 * 10 * 0.5 * 1.5, "N_per_unit_f_range_before_efficiency": sorted([float(ref["N_occ_OGLEIV_per_unit_f (rate norm.)"]), float(ref["N_occ_OGLEIV_per_unit_f (Smith norm.)"])]), "inputs": {"N_sources": N_SRC, "years": T_YR, "mean_rate_per_star_per_yr_ASSUMED": GAMMA, "v_rel_km_s_ASSUMED": 150.0, "efficiency_ASSUMED": 0.1, "N_sources_and_years_source": "Mroz et al. 2019 Sec. 2 (upper bounds)"}, "status": ("forecast, not a limit: no efficiency-calibrated search for chance total occultations in bulge " "survey data has been published (Smith et al. 2021 computed the rate for dark occulters at " "stellar density and noted that no concerted search exists)"), } summary = {"lens_light": lens, "optical_depth_remark": tau, "occultation_forecast": forecast} with open(os.path.join(DATA, "microlensing_summary.json"), "w") as f: json.dump(summary, f, indent=1) print("[microlensing]", json.dumps(summary, indent=1)) pd.set_option("display.width", 250) print(occ.to_string(index=False, float_format=lambda x: f"{x:.3g}")) # ---- figure --------------------------------------------------------------------------------- fig, ax = plt.subplots(1, 2, figsize=(12, 4.5)) t = np.linspace(-4.5, 4.5, 1801) u0 = 0.3 u = np.sqrt(u0 ** 2 + t ** 2) for rho, c in ((0.0, "k"), (0.8, "#1f77b4"), (1.5, "#ff7f0e"), (3.0, "#d62728")): lab = "bare lens (ordinary event)" if rho == 0 else f"opaque enclosure, R$_s$ = {rho} r$_E$" ax[0].plot(t, lightcurve_opaque(u, rho), color=c, lw=1.4, label=lab) ax[0].set_xlabel("time from closest approach [Einstein-radius crossing times]") ax[0].set_ylabel("flux of the background star (baseline = 1)") ax[0].set_ylim(-0.1, 3.8) ax[0].legend(fontsize=8, loc="upper right") ax[0].set_title("A lens inside an opaque sphere (impact parameter 0.3 r$_E$)", fontsize=10) Tg = np.logspace(np.log10(15), np.log10(600), 200) for (M, L, tag, c) in ((1.0, 1.0, "1 M$_\\odot$, 1 L$_\\odot$", "#d62728"), (0.6, 0.08, "0.6 M$_\\odot$, 0.08 L$_\\odot$", "#ff7f0e"), (0.3, 0.011, "0.3 M$_\\odot$, 0.011 L$_\\odot$", "#1f77b4")): rho = r_enclosure_au(L, Tg) / r_einstein_au(M, 6.0) uc = np.where(rho > 1, rho - 1 / rho, np.nan) ax[1].loglog(Tg, uc * n_true, color=c, label=tag) ax[1].axhline(n_true, color="k", ls=":", lw=0.8) ax[1].text(16, n_true * 1.15, "microlensing events occurring in the same data (assumed)", fontsize=7.5) ax[1].axvspan(15, 36, color="0.85", zorder=0) ax[1].text(16.5, 2.2e3, "below the AllWISE 22 um\nflux limit beyond 100 pc\n(1 L$_\\odot$)", fontsize=7.5) ax[1].set_xlabel("enclosure temperature [K]") for (M, L, c) in ((1.0, 1.0, "#d62728"), (0.6, 0.08, "#ff7f0e"), (0.3, 0.011, "#1f77b4")): rs = r_enclosure_au(L, Tg) rho = rs / r_einstein_au(M, 6.0) ax[1].loglog(Tg, np.where(rho > 1, SMITH_PER_SOURCE_YR * (rs / 0.5) * N_SRC * T_YR, np.nan), color=c, ls="--", lw=1) ax[1].plot([], [], "k-", label="normalised to the microlensing rate") ax[1].plot([], [], "k--", lw=1, label="normalised to Smith et al. (2021)") ax[1].set_ylabel("total occultations in OGLE-IV per unit enclosed fraction\n(before detection efficiency)") ax[1].set_ylim(1e3, 1e7) ax[1].set_xlim(15, 600) ax[1].set_xticks([20, 30, 50, 100, 200, 500]) ax[1].set_xticks([], minor=True) ax[1].set_xticklabels(["20", "30", "50", "100", "200", "500"]) ax[1].legend(fontsize=8, title="enclosed star (lens at 6 kpc)", title_fontsize=8) ax[1].set_title("Forecast (not a limit): cold enclosures occult instead of lensing", fontsize=10) ax[1].grid(alpha=0.3, which="both") fig.tight_layout() fig.savefig(os.path.join(FIGS, "fig4_microlensing_occultation.png"), dpi=150) plt.close(fig) if __name__ == "__main__": main()