"""Part 1b: occultation optical depth, event rate and duration distribution per monitored star, for bulge, disc-plane, field-star and LMC sightlines, with validation of the Galactic model against (i) the OGLE-IV microlensing optical depth and event rate field by field (Mroz et al. 2019), (ii) the chance-occultation rate of Smith et al. (2021), and (iii) the cold-cloud event numbers of Kerins, Binney & Silk (2002) as used by Drake & Cook (2003). Outputs: data/model_vs_ogle4_microlensing.csv, data/occultation_rates.csv, data/occultation_rates_population.csv, data/rate_sanity_checks.json, figures/fig2_rates.png """ import json import os import numpy as np import pandas as pd import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from astropy.coordinates import SkyCoord import astropy.units as u from common import DATA, FIGS, RAW import galmodel as g import popmodel as pm R_GRID = [0.3, 1, 2, 3, 5, 10, 20, 50, 100, 200, 500, 1000] T_GRID = [20, 30, 50, 75, 100, 150, 200, 300] N_POP, N_MASS = 60000, 24 def wq(t, w, p): o = np.argsort(t) c = np.cumsum(w[o]) return float(np.interp(p * c[-1], c, t[o])) if c[-1] > 0 else np.nan def validate_microlensing(): """Model versus OGLE-IV field-by-field optical depth and event rate.""" f = pd.read_csv(os.path.join(RAW, "mroz2019_fields.csv")).rename(columns={"tE-mean": "tE_mean"}) f = f[(f.Nev >= 20) & (f.tau > 0)].copy() rows = [] for r in f.itertuples(): sl = g.Sightline(r.GLON, r.GLAT, n=40000, seed=int(r.recno)) rows.append({"field": r.Field, "l": r.GLON, "b": r.GLAT, "Nev": r.Nev, "tau_obs_1e-6": r.tau, "e_tau": r.e_tau, "tau_model_1e-6": sl.ml_tau() * 1e6, "gamma_obs_1e-6": r.Gam / 10.0, "e_gamma": r.e_Gam / 10.0, "gamma_model_1e-6": sl.ml_rate() * 1e6, "tE_obs_d": r.tE_mean, "tE_model_d": sl.ml_mean_te(), "column_ms_stars_pc2": sl.column()}) v = pd.DataFrame(rows) v["ratio_tau"] = v["tau_obs_1e-6"] / v["tau_model_1e-6"] v["ratio_gamma"] = v["gamma_obs_1e-6"] / v["gamma_model_1e-6"] v.to_csv(os.path.join(DATA, "model_vs_ogle4_microlensing.csv"), index=False, float_format="%.4g") return v def calibration(v): """Global observed/model microlensing event-rate ratio (event-weighted).""" w = v.Nev.values return {"ratio_gamma_weighted": float(np.sum(w * v.ratio_gamma) / w.sum()), "ratio_gamma_median": float(v.ratio_gamma.median()), "ratio_gamma_16_84": [float(v.ratio_gamma.quantile(0.16)), float(v.ratio_gamma.quantile(0.84))], "ratio_tau_weighted": float(np.sum(w * v.ratio_tau) / w.sum()), "ratio_tau_median": float(v.ratio_tau.median()), "n_fields": int(len(v))} def sightlines(): c = SkyCoord(270.1975825599437 * u.deg, -30.602547973764167 * u.deg).galactic # VVV-WIT-08 (SIMBAD) return { "bulge l=1 b=-3": dict(l=1.0, b=-3.0), "bulge l=0 b=-1.5": dict(l=0.0, b=-1.5), "bulge l=0 b=-6": dict(l=0.0, b=-6.0), "VVV-WIT-08 position": dict(l=float(c.l.deg), b=float(c.b.deg)), "disc plane l=320 b=-0.5": dict(l=320.0, b=-0.5), "disc star at 3 kpc, l=30 b=2": dict(l=30.0, b=2.0, d_src=3.0), "field star at 1 kpc, l=90 b=10": dict(l=90.0, b=10.0, d_src=1.0), "field star at 1 kpc, l=90 b=30": dict(l=90.0, b=30.0, d_src=1.0), "field star at 0.3 kpc, l=90 b=30": dict(l=90.0, b=30.0, d_src=0.3), "Kepler field star at 1 kpc": dict(l=76.3, b=13.5, d_src=1.0), } def main(): v = validate_microlensing() cal = calibration(v) print("[rates] model vs OGLE-IV microlensing:", json.dumps(cal)) CAL = cal["ratio_gamma_weighted"] rows, prow, dur_store = [], [], {} for name, kw in sightlines().items(): sl = g.Sightline(n=200000, seed=11, **kw) is_bulge = name.startswith("bulge") or name.startswith("VVV") scale = CAL if is_bulge else 1.0 # empirical normalisation only where microlensing calibrates it for R in R_GRID: t, w = sl.durations(R, 0.3, True) rows.append({"sightline": name, "R_au": R, "column_ms_stars_pc2": sl.column() * scale, "rate_per_star_yr_f1": w.sum() * scale, "rate_per_star_yr_f1_no_lensing": sl.occ_rate(R, lensing=False) * scale, "rate_per_star_yr_f1_lens_1Msun": sl.occ_rate(R, 1.0) * scale, "tau_f1": sl.occ_tau(R, 0.3) * scale, "dur_d_10": wq(t, w, 0.1), "dur_d_50": wq(t, w, 0.5), "dur_d_90": wq(t, w, 0.9), "frac_dur_gt_10d": float(w[t > 10].sum() / max(w.sum(), 1e-300)), "mean_vperp_kms": float(np.average(sl.vperp, weights=sl.w_col))}) if name in ("bulge l=1 b=-3", "field star at 1 kpc, l=90 b=30") and R in (3, 10, 50, 200): dur_store[(name, R)] = (t, w) if name in ("bulge l=1 b=-3", "field star at 1 kpc, l=90 b=30", "disc plane l=320 b=-0.5"): slp = g.Sightline(n=N_POP, seed=12, **kw) # smaller sample: the mass loop multiplies the cost for T in T_GRID: for tag, lo, hi in (("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)): rate, t, w = g.population_rate(slp, T, n_m=N_MASS, lo=lo, hi=hi) prow.append({"sightline": name, "population": tag, "T_K": T, "rate_per_star_yr_f1": rate * scale, "dur_d_10": wq(t, w, 0.1), "dur_d_50": wq(t, w, 0.5), "dur_d_90": wq(t, w, 0.9)}) # LMC for R in R_GRID: mw, se, sl = g.lmc_rate(R) vse = np.sqrt(np.pi) * g.LMC["sigma_v"] rows.append({"sightline": "LMC (MW disc + self)", "R_au": R, "column_ms_stars_pc2": sl.column() + g.LMC["sigma_n_self"], "rate_per_star_yr_f1": mw + se, "rate_per_star_yr_f1_no_lensing": np.nan, "rate_per_star_yr_f1_lens_1Msun": np.nan, "tau_f1": np.nan, "dur_d_10": np.nan, "dur_d_50": 0.5 * np.pi * R / (vse * 0.21095) * 365.25 * 0.8, "dur_d_90": np.nan, "frac_dur_gt_10d": np.nan, "mean_vperp_kms": vse}) tab = pd.DataFrame(rows) tab.to_csv(os.path.join(DATA, "occultation_rates.csv"), index=False, float_format="%.4g") pop = pd.DataFrame(prow) pop.to_csv(os.path.join(DATA, "occultation_rates_population.csv"), index=False, float_format="%.4g") pd.set_option("display.width", 250) print(tab[tab.R_au.isin([1, 3, 10, 100])][["sightline", "R_au", "column_ms_stars_pc2", "rate_per_star_yr_f1", "tau_f1", "dur_d_10", "dur_d_50", "dur_d_90", "mean_vperp_kms"]] .to_string(index=False, float_format=lambda x: f"{x:.3g}")) print(pop[pop.sightline == "bulge l=1 b=-3"].to_string(index=False, float_format=lambda x: f"{x:.3g}")) # ---- sanity checks --------------------------------------------------------------------------- c = sightlines()["VVV-WIT-08 position"] sl = g.Sightline(n=300000, seed=21, **c) t, w = sl.durations(0.5, 0.3, lensing=False) # Smith et al. assume massless 0.5 au occulters smith = {"position_l_b": [c["l"], c["b"]], "rate_all_per_source_yr": float(w.sum()), "rate_gt10d_per_source_yr": float(w[t > 10].sum()), "rate_gt100d_per_source_yr": float(w[t > 100].sum()), "smith2021_gt10d_depth_gt_0.5": 0.80e-6, "smith2021_gt100d_depth_gt_0.9": 1.2e-8, "note": "This model: point sources, model source distances, no magnitude cut; Smith et al. use Galaxia " "with Ks<17 sources and finite source radii."} # KBS disc clouds toward MACHO-like bulge sightlines kbs = {} for (l, b) in ((1.0, -2.5), (2.0, -3.0), (3.0, -4.0), (4.0, -6.0)): k = g.Sightline(l, b, n=200000, seed=31, lens=g.kbs_disc_lens(), src_pop="bar") s = g.Sightline(l, b, n=200000, seed=31) rk = k.occ_rate(g.KBS_R_AU, lensing=False) tk, wk = k.durations(g.KBS_R_AU, lensing=False) kbs[f"l={l:g} b={b:g}"] = {"kbs_rate_per_star_yr": rk, "kbs_raw_events_in_1.71e8_star_yr": rk * 1.71e8, "kbs_median_duration_d": wq(tk, wk, 0.5), "implied_DC03_efficiency_if_1e5_detected": 1e5 / (rk * 1.71e8), "enclosure_rate_R7au_f1_per_star_yr": s.occ_rate(7.0, 0.3) * CAL, "ratio_enclosure_R7_to_kbs": s.occ_rate(7.0, 0.3) * CAL / rk} checks = {"microlensing_calibration": cal, "smith2021": smith, "kbs_drake_cook": kbs, "model_at_l1_b-3": {"tau_ml": g.Sightline(1, -3, n=200000, seed=2).ml_tau(), "gamma_ml_per_yr": g.Sightline(1, -3, n=200000, seed=2).ml_rate(), "tau_obs_fit": 1.36e-6, "gamma_obs_fit": 13.4e-6}} with open(os.path.join(DATA, "rate_sanity_checks.json"), "w") as f: json.dump(checks, f, indent=1) print(json.dumps(checks, indent=1)) # ---- figure 2 -------------------------------------------------------------------------------- fig, ax = plt.subplots(1, 3, figsize=(16.5, 4.8)) ax[0].errorbar(v["gamma_model_1e-6"], v["gamma_obs_1e-6"], yerr=v["e_gamma"], fmt="o", ms=3, lw=0.6, color="#1f77b4") lim = [0.3, 40] ax[0].loglog(lim, lim, "k--", lw=0.8) ax[0].loglog(lim, [CAL * x for x in lim], "r-", lw=0.8, label=f"observed = {CAL:.2f} x model") ax[0].set_xlim(lim) ax[0].set_ylim(lim) ax[0].set_xlabel("model microlensing event rate [10$^{-6}$ per star per yr]") ax[0].set_ylabel("OGLE-IV measured rate (Mroz et al. 2019)") ax[0].set_title(f"Validation: {len(v)} OGLE-IV fields", fontsize=10) ax[0].legend(fontsize=8) ax[0].grid(alpha=0.3, which="both") show = ["bulge l=1 b=-3", "bulge l=0 b=-6", "disc plane l=320 b=-0.5", "disc star at 3 kpc, l=30 b=2", "field star at 1 kpc, l=90 b=30", "LMC (MW disc + self)"] cols = ["#d62728", "#ff7f0e", "#8c564b", "#2ca02c", "#1f77b4", "#9467bd"] for name, c in zip(show, cols): s = tab[tab.sightline == name] ax[1].loglog(s.R_au, s.rate_per_star_yr_f1, "-o", ms=3, color=c, label=name) s = tab[tab.sightline == "bulge l=1 b=-3"] ax[1].loglog(s.R_au, s.rate_per_star_yr_f1_no_lensing, ":", color="#d62728", lw=1) ax[1].axhline(13.4e-6, color="0.4", ls="--", lw=0.8) ax[1].text(0.32, 16e-6, "measured microlensing rate, bulge b = -3", fontsize=7.5, color="0.3") ax[1].set_xlabel("enclosure radius R [au]") ax[1].set_ylabel("total occultations per monitored star per yr, if f = 1") ax[1].set_title("Occultation rate (0.3 M$_\\odot$ lens; dotted: light bending ignored)", fontsize=10) ax[1].legend(fontsize=7.5) ax[1].grid(alpha=0.3, which="both") bins = np.logspace(0, 5, 51) for (name, R), ls in zip([("bulge l=1 b=-3", 3), ("bulge l=1 b=-3", 10), ("bulge l=1 b=-3", 50), ("bulge l=1 b=-3", 200), ("field star at 1 kpc, l=90 b=30", 10)], ["-", "-", "-", "-", "--"]): t, w = dur_store[(name, R)] h, _ = np.histogram(t, bins, weights=w) lab = f"R = {R} au, " + ("bulge" if name.startswith("bulge") else "field star at 1 kpc") ax[2].semilogx(np.sqrt(bins[1:] * bins[:-1]), h / h.sum(), ls, label=lab) for x, lab in ((365.25, "1 yr"), (3652.5, "10 yr")): ax[2].axvline(x, color="0.6", lw=0.7) ax[2].text(x * 1.05, 0.002, lab, fontsize=7.5, color="0.4") ax[2].set_xlabel("occultation duration [days]") ax[2].set_ylabel("fraction of events per logarithmic bin") ax[2].set_title("Duration distributions", fontsize=10) ax[2].legend(fontsize=7.5) ax[2].grid(alpha=0.3, which="both") fig.tight_layout() fig.savefig(os.path.join(FIGS, "fig2_rates.png"), dpi=150) plt.close(fig) if __name__ == "__main__": main()