"""Handle 2: the local mass budget (the "Oort limit"). A fully enclosed star still gravitates. Vertical stellar kinematics measure the total mass of the local disc; star counts and gas surveys measure the visible baryons. The difference bounds any optically invisible population with a stellar-like scale height. All numbers are published values, read at source (the arXiv identifier is given for each; the abstracts are archived by 01_download.py, and the table values of McKee et al. 2015 were read from the full text). Nothing here is re-measured from Gaia: this is a re-interpretation of the published budgets. The baseline rests ONLY on the 1991-2014 dynamical determinations as re-analysed by McKee et al. (2015); no Gaia-era dynamical determination enters any limit. Two Gaia-era numbers are listed for context (Buch et al. 2019; the de Salas & Widmark 2021 review). Revision (referee point on the Bovy & Rix row): the disc component found by the dynamical authors is now compared with the baryons exactly as McKee et al. do in their Section 7.2, i.e. with the ionised gas (scale height 1590 pc, counted as "halo" by the dynamics) removed: Sigma_b,disc = 44.4 (whole disc) or 42.1 (within 1.1 kpc). A determination quoted for the whole disc (Bovy & Rix 2013: McKee Table 4 footnote d; Bienayme et al. 2014; Holmberg & Flynn 2004) is compared with 44.4, one quoted within 1.1 kpc (Kuijken & Gilmore; Zhang et al.) with 42.1. This reproduces McKee's published values 0 +- 5.3 (Bienayme), 13.1 +- 6.9 (Zhang) and the published average of four, 6.4. The first pass had subtracted 43.8 (all baryons within 1.1 kpc) from the whole-disc Bovy & Rix value. The referee's alternative (subtract all 47.1 baryons) is kept as a variant row. Estimators (Msun/pc^2 or Msun/pc^3): S1 Sigma_enc = Sigma_disc,dyn(<1.1 kpc) - Sigma_baryon,visible(<1.1 kpc) (halo subtracted by the dynamical authors); baseline S2 rho_enc = rho_0,dyn - rho_baryon,visible - rho_DM,halo (mid-plane) S3 Sigma_enc = Sigma_total(<1.1 kpc) - Sigma_baryon,visible (no dark halo at all near the Sun: the extreme upper bound) For each, the 95% upper limit is quoted (i) as the one-sided Gaussian limit and (ii) for a flat prior restricted to Sigma_enc >= 0, and is converted into a fraction of all stars, f = Sigma_enc / (Sigma_enc + Sigma_visible stars), by Monte Carlo over the uncertainties of both terms. """ import json import os import sys import numpy as np import pandas as pd sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) from common import DATA RNG = np.random.default_rng(20261012) NMC = 400000 # ---- compiled literature values ------------------------------------------------------------- # quantity, value, 1-sigma, unit, source, arXiv id, note COMPILATION = [ # baryons (McKee, Parravano & Hollenbach 2015, ApJ 814, 13) ("Sigma visible stars (MS + giants)", 27.0, 2.7, "Msun/pc2", "McKee+2015 eq. 36", "1509.05334", "25.4 within 1.1 kpc (Table 1)"), ("Sigma white dwarfs", 4.9, 0.6, "Msun/pc2", "McKee+2015 abstract", "1509.05334", "0.8 quoted in their Sec. 3"), ("Sigma brown dwarfs", 1.2, np.nan, "Msun/pc2", "McKee+2015 abstract", "1509.05334", "approximate"), ("Sigma stars + remnants", 33.4, 3.0, "Msun/pc2", "McKee+2015 Table 1", "1509.05334", "31.2 within 1.1 kpc"), ("Sigma gas", 13.7, 1.6, "Msun/pc2", "McKee+2015 abstract", "1509.05334", "12.6 within 1.1 kpc"), ("Sigma baryons", 47.1, 3.4, "Msun/pc2", "McKee+2015 abstract", "1509.05334", "43.8 +- 3.4 within 1.1 kpc; 43.1 +- 3.3 within 1.0 kpc"), ("rho visible stars (mid-plane)", 0.036, np.nan, "Msun/pc3", "McKee+2015 Table 1", "1509.05334", "no uncertainty printed for this entry"), ("rho stars + remnants (mid-plane)", 0.043, 0.004, "Msun/pc3", "McKee+2015 eq. 38", "1509.05334", "printed as +-0.04; the text says 'we have estimated the error to be 10%', i.e. +-0.004"), ("Sigma baryons in the disc, ionised gas removed", 44.4, 3.4, "Msun/pc2", "McKee+2015 Sec. 7.2", "1509.05334", "42.1 within 1.1 kpc; uncertainty taken from Sigma baryons"), ("rho baryons (mid-plane)", 0.084, 0.012, "Msun/pc3", "McKee+2015 eq. 46", "1509.05334", ""), ("rho dark-matter halo", 0.013, 0.003, "Msun/pc3", "McKee+2015 abstract", "1509.05334", ""), ("rho main-sequence stars (mid-plane)", 0.040, 0.002, "Msun/pc3", "Bovy 2017", "1704.05063", "Gaia DR1 star counts"), # dynamical: disc ("baryon-like") component after the authors' own halo subtraction (McKee+2015 Table 4) ("Sigma_disc dyn (<1.1 kpc) Kuijken & Gilmore 1991", 48.0, 8.0, "Msun/pc2", "McKee+2015 Table 4", "1509.05334", ""), ("Sigma_disc dyn Bovy & Rix 2013", 51.0, 4.0, "Msun/pc2", "McKee+2015 Table 4; Bovy & Rix 2013", "1309.0809", "38 +- 4 stars + remnants in the original"), ("Sigma_disc dyn (<1.0 kpc) Zhang+2013", 55.0, 5.0, "Msun/pc2", "McKee+2015 Table 4", "1509.05334", ""), ("Sigma_disc dyn (<1.0 kpc) Bienayme+2014", 43.0, 4.0, "Msun/pc2", "McKee+2015 Table 4", "1509.05334", "44.4 +- 4.1 total disc"), ("Sigma_disc dyn Holmberg & Flynn 2004", 56.0, 6.0, "Msun/pc2", "Holmberg & Flynn 2004 abstract", "astro-ph/0405155", "53 visible in their own budget"), # dynamical: total column ("Sigma_total (<1.1 kpc) Kuijken & Gilmore 1991", 73.5, 6.0, "Msun/pc2", "McKee+2015 Table 4", "1509.05334", ""), ("Sigma_total (<1.1 kpc) Bovy & Rix 2013", 68.0, 4.0, "Msun/pc2", "Bovy & Rix 2013 abstract", "1309.0809", ""), ("Sigma_total (<1.0 kpc) Zhang+2013", 69.0, 6.0, "Msun/pc2", "McKee+2015 Table 4", "1509.05334", ""), ("Sigma_total (<1.0 kpc) Bienayme+2014", 70.5, 1.0, "Msun/pc2", "McKee+2015 Table 4", "1509.05334", ""), ("Sigma_total (<1.1 kpc) Holmberg & Flynn 2004", 74.0, 6.0, "Msun/pc2", "Holmberg & Flynn 2004 abstract", "astro-ph/0405155", ""), # dynamical: mid-plane density (the Oort limit proper) ("rho_0 dyn Holmberg & Flynn 2000", 0.102, 0.010, "Msun/pc3", "Holmberg & Flynn 2000 abstract", "astro-ph/9812404", "0.095 visible in their own budget"), ("rho_0 dyn Bienayme+2014", 0.091, 0.006, "Msun/pc3", "McKee+2015 Table 5", "1509.05334", ""), ("rho_0 dyn Garbari+2012", 0.120, 0.018, "Msun/pc3", "McKee+2015 Table 5", "1509.05334", "+0.016/-0.019; strongly coupled to the assumed baryons"), # published dark-disc results quoted for context ("Sigma disc dark matter (average of four)", 6.4, 6.8, "Msun/pc2", "McKee+2015 Sec. 7.2", "1509.05334", "Bienayme 0 +- 5.3; Zhang 13.1 +- 6.9"), ("Sigma thin dark disc, h ~ 10 pc (95% excl.)", 6.0, np.nan, "Msun/pc2", "Schutz+2018 abstract", "1711.03103", "context only: thin disc"), ("Sigma thin dark disc, h < 100 pc (95% excl.)", 12.0, np.nan, "Msun/pc2", "Buch+2019 abstract", "1808.05603", "5-12 depending on scale height; context only"), ("Sigma thin dark disc, non-equilibrium analysis", 14.0, np.nan, "Msun/pc2", "Kramer & Randall 2016 abstract", "1604.01407", "context only"), ("Sigma surplus near the plane (|z| < 60 pc)", 7.0, 2.0, "Msun/pc2", "Widmark 2019 abstract", "1811.07911", "5-9; Gaia DR2; context only"), ("rho_DM local, Gaia DR2 A stars, |z| < 200 pc", 0.016, 0.010, "Msun/pc3", "Buch+2019 abstract", "1808.05603", "Gaia-era; context only; F stars give a significantly higher value"), ("rho_DM, range of local analyses", 0.0135, 0.0025, "Msun/pc3", "de Salas & Widmark 2021 abstract", "2012.11477", "0.011-0.016; mid-point and half-range"), ("rho_DM, range of global (rotation-curve) analyses", 0.0105, 0.0025, "Msun/pc3", "de Salas & Widmark 2021 abstract", "2012.11477", "0.008-0.013; mid-point and half-range"), ] # visible baryons within the height at which each dynamical value is quoted SIG_B = {"1.1": (43.8, 3.4), "1.0": (43.1, 3.3)} # all baryons (used for S3, total column) SIG_B_DISC = {"whole": (44.4, 3.4), "1.1": (42.1, 3.4)} # ionised gas removed (McKee Sec. 7.2), for S1 # visible stars within 1.1 kpc: 25.4 (McKee Table 1); the source gives 27.0 +- 2.7 for the whole # column (eq. 36) and no separate error for 25.4, so the same 10 per cent is ASSUMED here SIG_VIS_11 = (25.4, 2.54) # visible stars, mid-plane: 0.036 (McKee Table 1, no error printed); 10 per cent ASSUMED, as McKee # et al. estimate for the stars + remnants density RHO_VIS = (0.036, 0.004) def ul_gauss(mu, sig, cl=0.95): from scipy.stats import norm return mu + norm.ppf(cl) * sig def ul_truncated(mu, sig, cl=0.95): """Upper limit for a Gaussian likelihood with a flat prior on x >= 0.""" from scipy.stats import norm lo = norm.cdf(0, mu, sig) return norm.ppf(lo + cl * (1 - lo), mu, sig) def frac_limit(mu_enc, sig_enc, vis, cl=0.95): """95% limit on f = x/(x + visible) with x >= 0 (rejection sampling of the truncated Gaussian).""" x = RNG.normal(mu_enc, sig_enc, NMC) v = RNG.normal(vis[0], vis[1], NMC) ok = (x >= 0) & (v > 0) return float(np.quantile(x[ok] / (x[ok] + v[ok]), cl)) def main(): comp = pd.DataFrame(COMPILATION, columns=["quantity", "value", "sigma", "unit", "source", "arxiv", "note"]) comp.to_csv(os.path.join(DATA, "oort_compilation.csv"), index=False) rows = [] def add(label, kind, mu, sig, vis, note=""): rows.append({"estimator": label, "kind": kind, "excess": mu, "sigma": sig, "UL95_gaussian": ul_gauss(mu, sig), "UL95_positive_prior": ul_truncated(mu, sig), "visible_stars": vis[0], "f95_enclosed_fraction_of_stars": frac_limit(mu, sig, vis), "note": note}) # S1: disc component minus disc baryons (ionised gas removed), determination by determination, # following McKee et al. (2015) Sec. 7.2. (name, disc value, error, which baryon figure) s1 = [("Kuijken & Gilmore 1991", 48.0, 8.0, "1.1"), ("Bovy & Rix 2013", 51.0, 4.0, "whole"), ("Zhang+2013", 55.2, 6.0, "1.1"), ("Bienayme+2014", 44.4, 4.1, "whole"), ("Holmberg & Flynn 2004", 56.0, 6.0, "whole")] ex = [] for name, v, s, h in s1: b = SIG_B_DISC[h] mu, sig = v - b[0], float(np.hypot(s, b[1])) ex.append(mu) add(f"S1 disc - baryons: {name}", "surface density", mu, sig, SIG_VIS_11, f"dynamical {v} +- {s} minus {b[0]} +- {b[1]}") ex = np.array(ex) mean4 = float(ex[:4].mean()) assert abs(mean4 - 6.4) < 0.06, mean4 # reproduces the published average of four # combination: mean of the determinations; error = the larger of the scatter and the typical # individual error (the determinations share tracers and the baryon model, so no 1/sqrt(N)) sig_comb = max(float(np.std(ex, ddof=1)), float(np.median([r["sigma"] for r in rows]))) add("S1 disc - baryons: mean of five (this work)", "surface density", float(ex.mean()), sig_comb, SIG_VIS_11, "error = max(scatter, median individual error)") add("S1 as published: McKee+2015 average of four", "surface density", 6.4, 6.8, SIG_VIS_11, "BASELINE") # variant asked for by the referee: Bovy & Rix whole-disc value minus ALL baryons (47.1 +- 3.4) add("S1 variant: Bovy & Rix 2013 minus all baryons (47.1)", "surface density, variant", 51.0 - 47.1, float(np.hypot(4.0, 3.4)), SIG_VIS_11, "referee's subtraction; includes ionised gas the dynamics assigns to the halo") # S2: mid-plane densities for name, v, s in (("Holmberg & Flynn 2000", 0.102, 0.010), ("Bienayme+2014", 0.091, 0.006), ("Garbari+2012", 0.120, 0.018)): mu = v - 0.084 - 0.013 sig = float(np.sqrt(s ** 2 + 0.012 ** 2 + 0.003 ** 2)) add(f"S2 rho_0 - baryons - halo: {name}", "mid-plane density", mu, sig, RHO_VIS) # Gaia-era context (not used for any limit) add("S2 Gaia-era context: Buch+2019 rho_DM(A stars) - halo", "mid-plane density, context", 0.016 - 0.013, float(np.sqrt(0.010 ** 2 + 0.003 ** 2)), RHO_VIS, "Gaia DR2; depends on their baryon model; context only") add("S4 Gaia-era context: local - global dark-matter density", "mid-plane density, context", 0.0135 - 0.0105, float(np.hypot(0.0025, 0.0025)), RHO_VIS, "de Salas & Widmark 2021 ranges; NOT a limit: local analyses fit the disc mass, which would absorb enclosed stars") # S3: no local dark halo at all for name, v, s, h in (("Bovy & Rix 2013", 68.0, 4.0, "1.1"), ("Bienayme+2014", 70.5, 1.0, "1.0"), ("Kuijken & Gilmore 1991", 73.5, 6.0, "1.1")): b = SIG_B[h] add(f"S3 total - baryons (no halo): {name}", "surface density, extreme", v - b[0], float(np.hypot(s, b[1])), SIG_VIS_11, "assumes zero dark-matter halo within 1.1 kpc") res = pd.DataFrame(rows) res.to_csv(os.path.join(DATA, "oort_limits.csv"), index=False, float_format="%.4g") base = res[res.note == "BASELINE"].iloc[0] s1_rows = res[res.estimator.str.startswith("S1 disc") | res.estimator.str.startswith("S1 as published")] summary = { "baseline": {k: (float(v) if isinstance(v, (int, float, np.floating)) else v) for k, v in base.items()}, "S1_f95_range": [float(s1_rows.f95_enclosed_fraction_of_stars.min()), float(s1_rows.f95_enclosed_fraction_of_stars.max())], "S2_f95_range": [float(res[res.estimator.str.startswith("S2 rho_0")].f95_enclosed_fraction_of_stars.min()), float(res[res.estimator.str.startswith("S2 rho_0")].f95_enclosed_fraction_of_stars.max())], "S1_UL95_range_Msun_pc2": [float(s1_rows.UL95_positive_prior.min()), float(s1_rows.UL95_positive_prior.max())], "S3_UL95_range_Msun_pc2": [float(res[res.estimator.str.startswith("S3")].UL95_positive_prior.min()), float(res[res.estimator.str.startswith("S3")].UL95_positive_prior.max())], "BovyRix_McKee_procedure": {k: float(v) for k, v in res[res.estimator == "S1 disc - baryons: Bovy & Rix 2013"].iloc[0][ ["excess", "sigma", "UL95_positive_prior", "f95_enclosed_fraction_of_stars"]].items()}, "BovyRix_referee_variant": {k: float(v) for k, v in res[res.estimator.str.startswith("S1 variant")].iloc[0][ ["excess", "sigma", "UL95_positive_prior", "f95_enclosed_fraction_of_stars"]].items()}, "gaia_era_context": {r_["estimator"]: {"excess": float(r_["excess"]), "sigma": float(r_["sigma"]), "f95": float(r_["f95_enclosed_fraction_of_stars"])} for _, r_ in res[res.kind.str.contains("context")].iterrows()}, "mean_of_four_reproduced": mean4, "gaia_era_determinations_in_baseline": 0, "S3_f95_range": [float(res[res.estimator.str.startswith("S3")].f95_enclosed_fraction_of_stars.min()), float(res[res.estimator.str.startswith("S3")].f95_enclosed_fraction_of_stars.max())], # number of enclosed solar-mass-equivalents per visible star mass "baseline_enclosed_to_visible_mass_ratio_95": float(base.UL95_positive_prior / SIG_VIS_11[0]), } with open(os.path.join(DATA, "oort_summary.json"), "w") as f: json.dump(summary, f, indent=1) pd.set_option("display.width", 250) print(res.drop(columns=["note"]).to_string(index=False, float_format=lambda x: f"{x:.4g}")) print("[oort]", json.dumps(summary, indent=1)) if __name__ == "__main__": main()