"""Added in the revision: what the AllWISE catalogue actually contains at the position in colour space where a complete, thermally radiating enclosure would sit. The first pass compared the gravitational limits with the infrared limits of Suazo et al. (2022). The referee pointed out that those limits concern PARTIAL spheres around stars that have a Gaia counterpart, and say nothing about a complete enclosure, which has no optical source at all. A complete blackbody enclosure of 30-300 K is dark at 3.4 and 4.6 um and bright at 22 um. This script counts every AllWISE source of that kind in 144 non-overlapping 5-degree cones at |b| > 31 deg and converts the count into a CEILING on the fraction of stars that can be such enclosures: f(T) <= N95[ N_obs(T) ] / ( C * N_stars(T) ) N_obs(T) : AllWISE sources with W4 < 7.0 mag (13.3 mJy) whose W1, W2 and W3 measurements or upper limits are compatible, within +-1 mag, with a blackbody of temperature T. EVERY such source is counted as a possible enclosure, whatever it really is (most are expected to be dusty galaxies, the "W1W2-dropouts" of Eisenhardt et al. 2012, and catalogue artefacts), so this is an upper limit and not a search. N_stars(T) : number of main-sequence stars in the same cones whose enclosure would be brighter than W4 = 7.0, from the local luminosity function measured in the Gaia DR3 parent sample of this study (d < 100 pc) and the Juric et al. (2008) thin + thick disc. C : completeness, ASSUMED 0.8 (about 0.85 for chance blends with an unrelated 3.4-um source inside the beam, about 0.95 for catalogue completeness 0.7 mag above the 5-sigma limit). Fields contaminated by a bright moving object or scattered light (count of 22-um-bright, 3.4-um-dark sources above three times the median field) are dropped, with their area. Two less conservative variants are reported but not used for the headline ceiling: (i) sources within 2 arcmin of one another are counted once (about two thirds of the counted sources are fragments of extended objects, e.g. the planetary nebula NGC 246; independent enclosures cannot cluster like that at these surface densities); (ii) only isolated, point-like sources (ext flag 0, no other counted source within 2 arcmin). """ import json import os import sys import numpy as np import pandas as pd from scipy import stats sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) from common import RAW, DATA from mslib import MainSequence import wiselib W4_LIM = 7.0 # mag; 13.3 mJy, 2.7 times the nominal AllWISE 5-sigma limit TOL = 1.0 # mag tolerance on the blackbody colours C_COMPLETE = 0.8 # assumed completeness (see docstring) T_GRID = [30, 35, 40, 45, 50, 60, 75, 100, 125, 150, 200, 250, 300] R0 = 8000.0 # pc DISC = {"Juric+2008 thin+thick": dict(h1=300.0, L1=2600.0, h2=900.0, L2=3600.0, ft=0.12), "thin disc only, h = 250 pc": dict(h1=250.0, L1=2600.0, h2=900.0, L2=3600.0, ft=0.0)} CLASSES = {"all 0.2-1.5 Msun": (0.2, 1.5), "Sun-like 0.8-1.2 Msun": (0.8, 1.2), "M dwarfs 0.2-0.6 Msun": (0.2, 0.6)} def poisson_ul(n, cl=0.95): return stats.chi2.ppf(cl, 2 * (n + 1)) / 2.0 def ms_offset(ms, bprp, MG): order = np.argsort(ms._bprp) return MG - np.interp(bprp, ms._bprp[order], ms._c_MG[order]) def local_luminosity_function(ms): """Number density [pc^-3] of main-sequence stars in cells of absolute magnitude, from the Gaia DR3 parent sample within 100 pc (complete to G = 17.5, i.e. M_G = 12.5 at 100 pc).""" p = pd.read_csv(os.path.join(RAW, "parent_cells_300pc.csv")) p = p[p.dbin < 4].copy() # d < 100 pc p["MG"] = (p.mgbin + 0.5) * 0.25 p["bprp"] = (p.cbin + 0.5) * 0.1 p["dMS"] = ms_offset(ms, p.bprp.values, p.MG.values) p["M"] = ms.mass_of_MG(p.MG.values) p = p[(p.dMS > -1.3) & (p.dMS < 0.7) & (p.M >= 0.2) & (p.M < 1.5)] g = p.groupby("mgbin", as_index=False).agg(n=("n", "sum"), MG=("MG", "first"), M=("M", "first")) g["n0_pc3"] = g.n / (4.0 / 3.0 * np.pi * 100.0 ** 3) g["L"] = 10 ** ms.logL_of_mass(g.M.values) return g def density_ratio(r, l_deg, b_deg, m): """Stellar density relative to the solar neighbourhood at heliocentric distance r [pc].""" l, b = np.radians(l_deg), np.radians(b_deg) x = r * np.cos(b) R = np.sqrt(R0 ** 2 + x ** 2 - 2 * R0 * x * np.cos(l)) z = np.abs(r * np.sin(b)) thin = (1 - m["ft"]) * np.exp(-(R - R0) / m["L1"] - z / m["h1"]) thick = m["ft"] * np.exp(-(R - R0) / m["L2"] - z / m["h2"]) return thin + thick def effective_volume(fields, m): """V(d) = sum over fields of Omega * Int_0^d nu(r)/nu_0 r^2 dr, tabulated on a distance grid.""" r = np.linspace(0.0, 12000.0, 6001) V = np.zeros_like(r) for _, f in fields.iterrows(): omega = 2 * np.pi * (1 - np.cos(np.radians(f.radius_deg))) integrand = density_ratio(r, f.l, f.b, m) * r ** 2 V += omega * np.concatenate([[0.0], np.cumsum(0.5 * (integrand[1:] + integrand[:-1]) * np.diff(r))]) return r, V def compatible(d, T): """Is each source compatible with a blackbody of temperature T (see docstring)?""" ok = np.ones(len(d), bool) for k, band in ((1, "W1"), (2, "W2"), (3, "W3")): pred = d.W4mag.values + wiselib.colour(band, "W4", T) mag = d[f"W{k}mag"].values snr = d[f"snr{k}"].values det = np.isfinite(mag) & (snr >= 3) ok_det = np.abs(mag - pred) <= TOL ok_lim = ~np.isfinite(mag) | (pred >= mag - TOL) # prediction not brighter than the limit ok &= np.where(det, ok_det, ok_lim) return ok def n_groups(d, link_arcmin=2.0): """Number of friends-of-friends groups (linking length 2 arcmin) and number of isolated members.""" if len(d) == 0: return 0, np.zeros(0, bool) from astropy.coordinates import SkyCoord, search_around_sky import astropy.units as u c = SkyCoord(d.RAJ2000.values * u.deg, d.DEJ2000.values * u.deg) i, j, _, _ = search_around_sky(c, c, link_arcmin * u.arcmin) parent = np.arange(len(d)) def find(a): while parent[a] != a: parent[a] = parent[parent[a]] a = parent[a] return a for a, b in zip(i, j): if a < b: ra, rb = find(a), find(b) if ra != rb: parent[rb] = ra roots = np.array([find(a) for a in range(len(d))]) size = np.bincount(roots, minlength=len(d))[roots] return int(len(np.unique(roots))), size == 1 def main(): ms = MainSequence() src = pd.read_csv(os.path.join(RAW, "w4dark_sources.csv")) fields = pd.read_csv(os.path.join(RAW, "w4dark_fields.csv")) n_rows = len(src) # ---- field hygiene ----------------------------------------------------------------------------- per = src.groupby("field").size().reindex(fields.field, fill_value=0) fields["n_w4lt7p7_w1dark"] = per.values med = float(np.median(per.values)) fields["contaminated"] = fields.n_w4lt7p7_w1dark > 3 * med clean = fields[~fields.contaminated] omega_field = 2 * np.pi * (1 - np.cos(np.radians(fields.radius_deg.iloc[0]))) * (180 / np.pi) ** 2 area = len(clean) * omega_field s = src[src.field.isin(clean.field) & (src.W4mag < W4_LIM)].copy() s["w1dark"] = ~(np.isfinite(s.W1mag) & (s.snr1 >= 3)) s["w2dark"] = ~(np.isfinite(s.W2mag) & (s.snr2 >= 3)) s["w3dark"] = ~(np.isfinite(s.W3mag) & (s.snr3 >= 3)) s["artefact_flag_w4"] = s.ccf.astype(str).str[3].isin(list("DPHO")) # is the W4 < 7.0 sample safely above the catalogue's SNR >= 5 threshold? edge = src[src.field.isin(clean.field) & (src.W4mag > W4_LIM - 0.15) & (src.W4mag < W4_LIM)] snr_edge = [float(x) for x in np.percentile(edge.snr4, [1, 5, 50])] if len(edge) else [np.nan] * 3 # cross-check against the published "W1W2-dropout" surface density: Eisenhardt et al. (2012) # select W1 > 17.4 and (W4 < 7.7 and W2 - W4 > 8.2, or W3 < 10.6 and W2 - W3 > 5.3) and quote about # 1000 objects over the sky outside the Galactic plane and bulge (about 0.03 per deg2) AFTER # their artefact vetting. Only the 22-um branch can be reproduced from this download. a_ = src[src.field.isin(clean.field)] w1faint = ~np.isfinite(a_.W1mag) | (a_.W1mag > 17.4) drop = w1faint & (a_.W4mag < 7.7) & (~np.isfinite(a_.W2mag) | (a_.W2mag - a_.W4mag > 8.2)) n_drop = int(drop.sum()) n_drop_clean = int((drop & ~a_.ccf.astype(str).str[3].isin(list("DPHOdpho")) & (a_.ex == 0)).sum()) # ---- stars --------------------------------------------------------------------------------------- lf = local_luminosity_function(ms) lf.to_csv(os.path.join(DATA, "w4dark_local_luminosity_function.csv"), index=False, float_format="%.5g") flim = wiselib.F0_JY["W4"] * 10 ** (-0.4 * W4_LIM) vols = {name: effective_volume(clean, m) for name, m in DISC.items()} rows = [] for T in T_GRID: comp = compatible(s, T) n_obs = int(comp.sum()) ng, iso = n_groups(s[comp]) n_iso_pt = int((iso & (s[comp].ex.values == 0)).sum()) row = {"T_K": T, "N_obs_compatible": n_obs, "N95": poisson_ul(n_obs), "N_obs_without_artefact_flag": int((comp & ~s.artefact_flag_w4).sum()), "N_groups_2arcmin": ng, "N_isolated_pointlike": n_iso_pt, "W3-W4_blackbody": wiselib.colour("W3", "W4", T), "dmax_pc_1Lsun": wiselib.dmax_pc("W4", flim, T, 1.0), "dmax_pc_0.01Lsun": wiselib.dmax_pc("W4", flim, T, 0.01)} dmax = np.array([wiselib.dmax_pc("W4", flim, T, L) for L in lf.L.values]) for dname, (r, V) in vols.items(): veff = np.interp(dmax, r, V) for cname, (mlo, mhi) in CLASSES.items(): inc = (lf.M.values >= mlo) & (lf.M.values < mhi) nstar = float((lf.n0_pc3.values * veff)[inc].sum()) tag = f"{cname} | {dname}" row[f"N_stars | {tag}"] = nstar row[f"f95 | {tag}"] = poisson_ul(n_obs) / (C_COMPLETE * nstar) if nstar > 0 else np.nan if dname.startswith("Juric"): row[f"f95 groups counted once | {cname}"] = poisson_ul(ng) / (C_COMPLETE * nstar) if nstar > 0 else np.nan row[f"f95 isolated point-like only | {cname}"] = poisson_ul(n_iso_pt) / (C_COMPLETE * nstar) if nstar > 0 else np.nan rows.append(row) tab = pd.DataFrame(rows) tab.to_csv(os.path.join(DATA, "w4dark_ceiling_by_T.csv"), index=False, float_format="%.4g") fields.to_csv(os.path.join(DATA, "w4dark_fields.csv"), index=False, float_format="%.6g") keep = ["field", "AllWISE", "RAJ2000", "DEJ2000", "W1mag", "snr1", "W2mag", "snr2", "W3mag", "snr3", "W4mag", "snr4", "ccf", "ex", "qph", "covW4"] cold = s[compatible(s, 100) | compatible(s, 50)] cold[keep].to_csv(os.path.join(DATA, "w4dark_sources_compatible_50K_or_100K.csv"), index=False) base = "Juric+2008 thin+thick" alt = "thin disc only, h = 250 pc" def pick(T, cname, dname=base): return float(tab.loc[tab.T_K == T, f"f95 | {cname} | {dname}"].iloc[0]) summary = { "selection": {"W4_lim_mag": W4_LIM, "W4_lim_mJy": flim * 1e3, "colour_tolerance_mag": TOL, "completeness_assumed": C_COMPLETE}, "n_fields": int(len(fields)), "n_fields_contaminated": int(fields.contaminated.sum()), "contaminated_fields_l_b": [[float(a), float(b)] for a, b in fields[fields.contaminated][["l", "b"]].values], "area_deg2": float(area), "sky_fraction": float(area / 41252.96), "n_rows_downloaded": int(n_rows), "median_rows_per_field": med, "n_W4lt_lim_clean_fields": int(len(s)), "n_W4lt_lim_w1dark_w2dark": int((s.w1dark & s.w2dark).sum()), "n_W4lt_lim_w1dark_w2dark_w3dark": int((s.w1dark & s.w2dark & s.w3dark).sum()), "surface_density_w1w2dark_per_deg2": float((s.w1dark & s.w2dark).sum() / area), "snr4_percentiles_1_5_50_within_0.15mag_of_limit": snr_edge, "w1w2_dropout_like_22um_branch": {"n": n_drop, "per_deg2": n_drop / area, "n_no_W4_artefact_flag_and_pointlike": n_drop_clean, "per_deg2_no_flag_pointlike": n_drop_clean / area, "published_all_branches_after_vetting_per_deg2": 1000.0 / 32000.0}, "local_density_pc3": {c: float(lf.n0_pc3[(lf.M >= a) & (lf.M < b)].sum()) for c, (a, b) in CLASSES.items()}, "ceiling_f95": {c: {str(T): pick(T, c) for T in T_GRID} for c in CLASSES}, "ceiling_f95_thin_disc_only": {c: {str(T): pick(T, c, alt) for T in T_GRID} for c in CLASSES}, "N_obs_by_T": {str(int(r_.T_K)): int(r_.N_obs_compatible) for _, r_ in tab.iterrows()}, "N_groups_by_T": {str(int(r_.T_K)): int(r_.N_groups_2arcmin) for _, r_ in tab.iterrows()}, "N_isolated_pointlike_by_T": {str(int(r_.T_K)): int(r_.N_isolated_pointlike) for _, r_ in tab.iterrows()}, "ceiling_f95_groups_counted_once": {c: {str(int(r_.T_K)): float(r_[f"f95 groups counted once | {c}"]) for _, r_ in tab.iterrows()} for c in CLASSES}, "ceiling_f95_isolated_pointlike_only": {c: {str(int(r_.T_K)): float(r_[f"f95 isolated point-like only | {c}"]) for _, r_ in tab.iterrows()} for c in CLASSES}, "N_stars_sunlike_by_T": {str(int(r_.T_K)): float(r_[f"N_stars | Sun-like 0.8-1.2 Msun | {base}"]) for _, r_ in tab.iterrows()}, "N_stars_mdwarf_by_T": {str(int(r_.T_K)): float(r_[f"N_stars | M dwarfs 0.2-0.6 Msun | {base}"]) for _, r_ in tab.iterrows()}, "status": ("ceiling from source counts, not a search: every colour-compatible AllWISE source is counted as an " "enclosure; assumes an isotropic single-temperature blackbody, the stated disc model and completeness"), } with open(os.path.join(DATA, "w4dark_summary.json"), "w") as f: json.dump(summary, f, indent=1) pd.set_option("display.width", 250) cols = ["T_K", "N_obs_compatible", "N_groups_2arcmin", "N_isolated_pointlike", "W3-W4_blackbody", "dmax_pc_1Lsun", "dmax_pc_0.01Lsun"] cols += [c for c in tab.columns if c.endswith(base)] print(tab[cols].to_string(index=False, float_format=lambda x: f"{x:.3g}")) print("[w4dark]", json.dumps({k: v for k, v in summary.items() if k not in ("N_stars_sunlike_by_T", "N_stars_mdwarf_by_T")}, indent=1)) if __name__ == "__main__": main()