"""A standard microlensing-style Galactic model (double-exponential disc + triaxial bar) giving, along any sightline, the number density and kinematics of main-sequence stars, and from them the occultation optical depth, event rate and duration distribution for opaque spheres. The same machinery returns the microlensing optical depth and event rate, which are compared with the OGLE-IV measurements (Mroz et al. 2019) as a validation and used as an empirical normalisation. Geometry: Galactic centre at the origin, Sun at (-R0, 0, 0), y towards l = 90 deg, z towards the north Galactic pole. Units: kpc, km/s, Msun, au. """ import numpy as np from common import (AU_PC, KMS_TO_PC_PER_YR, KMS_TO_AU_PER_YR, RE2_AU2_PER_MSUN_KPC, G_NEWT, M_SUN, C_LIGHT, PC) import popmodel as pm R0 = 8.2 # kpc (GRAVITY) V_CIRC = 233.0 # km/s SUN_PEC = (11.1, 12.24, 7.25) # Schoenrich et al. 2010 (U towards GC, V, W) # ---- disc: Zheng et al. (2001)-type profile used by Han & Gould (2003) -------------------------- DISC_HR = 2.75 # kpc DISC_H1, DISC_H2, DISC_BETA = 0.156, 0.439, 0.381 N0_MS = 0.085 # main-sequence stars per pc^3 at the Sun (10-pc census, Reyle et al. 2021) # ---- bar: Dwek et al. (1995) G2, rescaled to R0 -------------------------------------------------- BAR_X0, BAR_Y0, BAR_Z0 = (v * R0 / 8.5 for v in (1.58, 0.62, 0.43)) BAR_ANGLE = np.radians(25.0) BAR_MASS = 1.6e10 # Msun, stars + remnants + brown dwarfs (Portail et al. 2017; Valenti et al. 2016) BAR_OMEGA = 40.0 # km/s/kpc pattern-like mean streaming (solid body), capped at 100 km/s BAR_SIG = (110.0, 82.5, 66.3) # km/s along bar axes x', y', z' (Han & Gould 1995) DISC_SIG = (35.0, 25.0, 20.0) # km/s (R, phi, z) DISC_LAG = 10.0 # km/s asymmetric drift _POP = pm.ms_stars_per_msun() N_MS_PER_MSUN = _POP["n_ms_per_msun_total"] # 2.54 for Kroupa 0.08-1 Msun def _bar_norm(): # volume integral of exp(-rs^2/2) in units of x0*y0*z0 u = np.linspace(-6, 6, 241) U, V, W = np.meshgrid(u, u, u, indexing="ij") rs = (((U ** 2 + V ** 2) ** 2 + W ** 4)) ** 0.25 return np.exp(-0.5 * rs ** 2).sum() * (u[1] - u[0]) ** 3 _BAR_K = _bar_norm() BAR_RHO0 = BAR_MASS / (_BAR_K * BAR_X0 * BAR_Y0 * BAR_Z0 * 1e9) # Msun / pc^3 def xyz(l, b, D): l, b = np.radians(l), np.radians(b) return -R0 + D * np.cos(b) * np.cos(l), D * np.cos(b) * np.sin(l), D * np.sin(b) def n_disc(X, Y, Z): R = np.hypot(X, Y) prof = (1 - DISC_BETA) / np.cosh(Z / DISC_H1) ** 2 + DISC_BETA * np.exp(-np.abs(Z) / DISC_H2) return N0_MS * np.exp(-(R - R0) / DISC_HR) * prof # stars / pc^3 def _bar_coords(X, Y, Z): xp = -X * np.cos(BAR_ANGLE) + Y * np.sin(BAR_ANGLE) yp = X * np.sin(BAR_ANGLE) + Y * np.cos(BAR_ANGLE) return xp, yp, Z def rho_bar(X, Y, Z): xp, yp, zp = _bar_coords(X, Y, Z) rs = (((xp / BAR_X0) ** 2 + (yp / BAR_Y0) ** 2) ** 2 + (zp / BAR_Z0) ** 4) ** 0.25 return BAR_RHO0 * np.exp(-0.5 * rs ** 2) # Msun / pc^3 def n_bar(X, Y, Z): return N_MS_PER_MSUN * rho_bar(X, Y, Z) def rho_disc_mass(X, Y, Z): return n_disc(X, Y, Z) / N_MS_PER_MSUN # same mass function assumed def sample_velocity(comp, X, Y, Z, rng): """Galactocentric Cartesian velocity (km/s) of stars of component comp (0 disc, 1 bar).""" n = X.size R = np.maximum(np.hypot(X, Y), 1e-3) v = np.zeros((n, 3)) d = comp == 0 if d.any(): vr = rng.normal(0, DISC_SIG[0], d.sum()) vp = V_CIRC - DISC_LAG + rng.normal(0, DISC_SIG[1], d.sum()) vz = rng.normal(0, DISC_SIG[2], d.sum()) er = np.c_[X[d] / R[d], Y[d] / R[d]] ep = np.c_[Y[d] / R[d], -X[d] / R[d]] v[d, 0] = vr * er[:, 0] + vp * ep[:, 0] v[d, 1] = vr * er[:, 1] + vp * ep[:, 1] v[d, 2] = vz bsel = comp == 1 if bsel.any(): nb = bsel.sum() vstream = np.minimum(BAR_OMEGA * R[bsel], 100.0) ep = np.c_[Y[bsel] / R[bsel], -X[bsel] / R[bsel]] s = rng.normal(0, 1, (nb, 3)) * np.array(BAR_SIG) ca, sa = np.cos(BAR_ANGLE), np.sin(BAR_ANGLE) # bar axes: x' = (-cos a, sin a, 0), y' = (sin a, cos a, 0) v[bsel, 0] = vstream * ep[:, 0] + s[:, 0] * (-ca) + s[:, 1] * sa v[bsel, 1] = vstream * ep[:, 1] + s[:, 0] * sa + s[:, 1] * ca v[bsel, 2] = s[:, 2] return v V_SUN = np.array([SUN_PEC[0], V_CIRC + SUN_PEC[1], SUN_PEC[2]]) # Sun at (-R0,0,0): U is +x, rotation +y class Sightline: """Monte Carlo sample of (source, lens) pairs along a sightline. Each sample carries a weight w such that, for an opaque sphere of effective radius R_eff [au], rate per source per year = f * sum( w_rate * R_eff ) (R_eff may depend on the pair) optical depth = f * sum( w_tau * R_eff^2 ) with f the fraction of stars that are enclosed. """ def __init__(self, l, b, n=200000, seed=1, d_src=None, src_pop="model", dmax=14.0, beta=1.0, lens=None, v_src_override=None): rng = np.random.default_rng(seed) self.l, self.b = l, b el = np.array([-np.sin(np.radians(l)), np.cos(np.radians(l)), 0.0]) eb = np.array([-np.sin(np.radians(b)) * np.cos(np.radians(l)), -np.sin(np.radians(b)) * np.sin(np.radians(l)), np.cos(np.radians(b))]) # ---- sources if d_src is None: Dg = np.linspace(0.05, dmax, 2800) Xg, Yg, Zg = xyz(l, b, Dg) wd, wb = n_disc(Xg, Yg, Zg), n_bar(Xg, Yg, Zg) if src_pop == "bar": wd = wd * 0 ws = (wd + wb) * Dg ** (2 - 2 * beta) cdf = np.cumsum(ws) Ds = np.interp(rng.uniform(0, cdf[-1], n), cdf, Dg) else: Ds = np.full(n, float(d_src)) if np.isscalar(d_src) else rng.choice(np.asarray(d_src, float), n) Xs, Ys, Zs = xyz(l, b, Ds) if d_src is None or v_src_override is None: pb = n_bar(Xs, Ys, Zs) / np.maximum(n_bar(Xs, Ys, Zs) + n_disc(Xs, Ys, Zs), 1e-30) if src_pop == "bar": pb = np.ones(n) comp_s = (rng.uniform(size=n) < pb).astype(int) vs = sample_velocity(comp_s, Xs, Ys, Zs, rng) if v_src_override is not None: vs = np.tile(np.asarray(v_src_override, float), (n, 1)) # ---- lenses: uniform in distance along the line of sight, weighted by density Dl = rng.uniform(0, 1, n) * Ds Xl, Yl, Zl = xyz(l, b, Dl) if lens is None: nd, nb = n_disc(Xl, Yl, Zl), n_bar(Xl, Yl, Zl) ntot = nd + nb comp_l = (rng.uniform(size=n) < nb / np.maximum(ntot, 1e-30)).astype(int) vl = sample_velocity(comp_l, Xl, Yl, Zl, rng) else: # custom occulter population: dict(n=fn(X,Y,Z), vel=fn(X,Y,Z,rng)) ntot = lens["n"](Xl, Yl, Zl) nd, nb = ntot, 0 * ntot comp_l = np.zeros(n, int) vl = lens["vel"](Xl, Yl, Zl, rng) x = Dl / Ds vrel = vl - (x[:, None] * vs + (1 - x[:, None]) * V_SUN) vperp = np.hypot(vrel @ el, vrel @ eb) # km/s at the lens plane self.Ds, self.Dl, self.vperp, self.comp_l = Ds, Dl, vperp, comp_l self.n_lens = ntot # stars / pc^3 self.path_pc = Ds * 1e3 # MC volume factor (pc) self.re2_per_msun = RE2_AU2_PER_MSUN_KPC * Dl * (Ds - Dl) / Ds # au^2 per Msun # weights: integral over D_l -> mean over samples times path length self.w_col = ntot * self.path_pc / n # stars / pc^2 per sample self.w_rate = 2 * AU_PC * self.w_col * vperp * KMS_TO_PC_PER_YR # per au of R_eff, per yr self.w_tau = np.pi * AU_PC ** 2 * self.w_col # per au^2 of R_eff^2 self.rho_mass = (nd + nb) / N_MS_PER_MSUN # Msun / pc^3 (all stellar mass) self.n = n self.rng = rng # ---- generic numbers def column(self): return float(self.w_col.sum()) # foreground MS stars per pc^2 per source def r_eff(self, R, M): """Radius of the region in which the source is totally hidden: the major image must fall inside the sphere, theta_+ < R/r_E, i.e. impact parameter < R - r_E^2/R (zero if R < r_E).""" re2 = self.re2_per_msun * M return np.clip(R - re2 / R, 0.0, None) def occ_rate(self, R, M=0.3, lensing=True): reff = self.r_eff(R, M) if lensing else np.full(self.n, float(R)) return float(np.sum(self.w_rate * reff)) # per star per year, for f = 1 def occ_tau(self, R, M=0.3, lensing=True): reff = self.r_eff(R, M) if lensing else np.full(self.n, float(R)) return float(np.sum(self.w_tau * reff ** 2)) def durations(self, R, M=0.3, lensing=True): """Event durations [days] and rate weights (sum of weights = rate for f = 1).""" reff = self.r_eff(R, M) if lensing else np.full(self.n, float(R)) u = self.rng.uniform(0, 1, self.n) chord = 2 * reff * np.sqrt(1 - u ** 2) t = chord / np.maximum(self.vperp * KMS_TO_AU_PER_YR, 1e-9) * 365.25 return t, self.w_rate * reff # ---- microlensing validation (all lenses: BD + MS + remnants, Kroupa mass function) def ml_tau(self): k = 4 * np.pi * G_NEWT * M_SUN / C_LIGHT ** 2 / PC # (pc) per (Msun/pc^3 * pc^2) -> dimensionless x = self.Dl / self.Ds return float(np.sum(k * self.rho_mass * (self.Dl * 1e3) * (1 - x) * self.path_pc / self.n)) def ml_rate(self): mlens, wnum = pm.lens_mass_function() # number per Msun of stellar mass mean_sqrt_m_per_msun = np.sum(wnum * np.sqrt(mlens)) # sum over lenses of sqrt(M) per Msun re1 = np.sqrt(self.re2_per_msun) # au for 1 Msun w = 2 * AU_PC * self.rho_mass * self.path_pc / self.n * self.vperp * KMS_TO_PC_PER_YR return float(np.sum(w * re1) * mean_sqrt_m_per_msun) def ml_mean_te(self): """Rate-weighted mean Einstein-radius crossing time (days) for the validation.""" mlens, wnum = pm.lens_mass_function() a = np.sum(wnum * mlens) / np.sum(wnum * np.sqrt(mlens)) # / w = self.rho_mass * self.path_pc * np.sqrt(self.re2_per_msun) * self.vperp te1 = np.sqrt(self.re2_per_msun) / np.maximum(self.vperp * KMS_TO_AU_PER_YR, 1e-9) * 365.25 return float(np.sum(w * te1) / np.sum(w) * a) def population_rate(sl, T, lensing=True, n_m=60, lo=pm.M_MS_LO, hi=pm.M_MS_HI): """Rate per source per year for f = 1 when every main-sequence star in [lo, hi] is equally likely to be enclosed at temperature T; also returns pooled duration samples.""" m, w, R = pm.population_radii(T, n=n_m, lo=lo, hi=hi) rate, ts, ws = 0.0, [], [] for mi, wi, Ri in zip(m, w, R): t, wt = sl.durations(Ri, mi, lensing) rate += wi * wt.sum() ts.append(t) ws.append(wi * wt) return rate, np.concatenate(ts), np.concatenate(ws) # ---- the cold-cloud disc model of Kerins, Binney & Silk (2002) as used by Drake & Cook (2003) ---- KBS_N0 = 0.03 / 1e-3 # clouds per pc^3 at the Sun: 0.03 Msun/pc^3 in 1e-3 Msun clouds KBS_R_AU = 7.0 def kbs_disc_lens(): def n(X, Y, Z): R = np.hypot(X, Y) return KBS_N0 * np.exp(-(R - R0) / 2.5) / np.cosh(Z / 0.19) ** 2 def vel(X, Y, Z, rng): R = np.maximum(np.hypot(X, Y), 1e-3) v = rng.normal(0, 25.0, (X.size, 3)) v[:, 0] += 220.0 * Y / R v[:, 1] += -220.0 * X / R return v return {"n": n, "vel": vel} # ---- Large Magellanic Cloud sightline: Milky Way disc foreground plus LMC self-occultation -------- LMC = {"l": 280.47, "b": -32.89, "D": 50.0, "sigma_n_self": 190.0, # MS stars / pc^2 in front of a typical source in the central fields "sigma_v": 22.0, # km/s 1-D dispersion of LMC disc stars "dls_kpc": 0.5} # typical lens-source separation for self-events def lmc_rate(R, M=0.3, lensing=True, n=100000, seed=3): """Occultation rate per LMC source per year for f = 1: (MW disc part, LMC self part).""" # LMC systemic motion relative to the Sun is irrelevant for nearby lenses (x -> 0). sl = Sightline(LMC["l"], LMC["b"], n=n, seed=seed, d_src=LMC["D"], v_src_override=(0.0, 0.0, 0.0)) mw = sl.occ_rate(R, M, lensing) re2 = RE2_AU2_PER_MSUN_KPC * M * LMC["dls_kpc"] if lensing else 0.0 reff = max(R - re2 / R, 0.0) vmean = np.sqrt(np.pi) * LMC["sigma_v"] # mean |v1 - v2| in two dimensions self_rate = 2 * reff * AU_PC * LMC["sigma_n_self"] * vmean * KMS_TO_PC_PER_YR return mw, self_rate, sl