# Extend one monochromatic PBH limit to a lognormal mass distribution

This folder runs on its own with Python 3 and the standard library. It was
written and checked with Python 3.9.6 on 2026-11-12.

## What it computes

A published primordial black hole (PBH) limit usually assumes every black hole
has the same mass `M`. It gives `f_max(M)`, the largest fraction of the dark
matter they may make up. A formation model gives a spread of masses instead:
`psi(M)`, the dark-matter fraction per natural-log mass interval.

For a probe whose signal is a sum of independent contributions from each black
hole, Carr et al. 2017 (arXiv:1705.05567v3, Eq. 12) show that the spread is
allowed when

    R = Int dln M  psi(M) / f_max(M)  <=  1.

`extend_constraint.py` evaluates `R` for a lognormal `psi` against one table
and reports the largest total fraction `f = 1 / R(f = 1)` that the table allows.

## Inputs

| File | Content | SHA-256 |
| --- | --- | --- |
| `input/SUB.dat` | Subaru HSC microlensing of M31 (Smyth et al. 2020, arXiv:1910.01285) | `251f0ed1457a8846b55ccae5c79ce93c2b63e6456aa574a668f075e146344860` |
| `input/LIGO.dat` | LIGO-Virgo-KAGRA O3 merger limit (Andrés-Carcasona et al. 2024, arXiv:2405.05732), included only to show the refusal | `7ef737e1139eeefb34ba4b9500257b786948619d1b5ad47a75772f004a401f38` |
| `input/LICENSE-PBHconstraints-MIT.txt` | MIT license of the digitized tables | `79e29afbfe33a2fd4a42c3363b69d411ee58f5ddc6d520d40edbee2c09e2239a` |

The tables are unchanged copies from github.com/vianvask/PBHconstraints at
commit `ef2fb03aa0f6c3fab58e9b23ee64331e88c75a86` (MIT, Copyright (c) 2026
Ville Vaskonen). Each row is `log10(M / M_sun)  log10 f_max`. The table-to-paper
mapping comes from the repository's notebook; the papers themselves were read
only as arXiv abstracts for this tutorial.

## Units

- Mass: solar masses, stored as `log10 M`.
- `sigma`: lognormal width in natural-log mass. One `sigma` is `sigma / ln 10`
  decades; `sigma = 2` spans about 0.87 dex.
- `psi`: dimensionless per unit `ln M`, normalized so `Int psi dln M = f`.

## Run

    python3 test_tutorial.py
    python3 extend_constraint.py --cases --output rerun/results.json
    python3 extend_constraint.py --table input/SUB.dat --central -8 --width 0.5

`expected/results.json` holds the five tutorial cases; a rerun should match it
to about 1e-12 in relative terms on any IEEE double platform.

## Expected results

| `log10 M_c` | `sigma` | max `f` from SUB | monochromatic `f_max(M_c)` |
| --- | --- | --- | --- |
| -8 | 0.005 | 0.0094879 | 0.0094880 |
| -8 | 0.5 | 0.0096371 | 0.0094880 |
| -8 | 1.0 | 0.0098154 | 0.0094880 |
| -8 | 2.0 | 0.010949 | 0.0094880 |
| -6 | 0.5 | 0.15997 | 0.16890 |

The narrow case reproduces the table value to 2.2e-6. A wider spread is weaker
at `M_c = 1e-8` and stronger at `M_c = 1e-6`: whether spreading the mass helps
or hurts depends on how `1 / f_max` curves around `M_c`, so neither direction
can be assumed.

## Numerical error

Trapezoid rule on a uniform `log10 M` grid, 400 points per `sigma` over
`+-9 sigma`. Halving or doubling the grid changes `R` by less than 7e-8 for
`sigma = 0.5` and `2.0`. The truncation at `9 sigma` drops a Gaussian tail of
about 2e-19. These errors are far below the uncertainty of the digitized
curves, which this tutorial does not propagate.

## Refusals

`expected/failures.json` lists four commands that must fail:

| Input | Exit code | Reason |
| --- | --- | --- |
| `LIGO.dat` | 3 | A merger rate depends on pairs of black holes, so it is not linear in `psi` and Eq. 12 does not apply. |
| `--width 0` or `nan` | 2 | A lognormal needs a positive finite width. |
| `--central -3` | 5 | Less than half of the distribution lies inside the table; SUB has no data there and cannot bound it. |

A table whose masses do not strictly increase also fails, with exit code 4.

## Assumptions that change the answer

1. The probe is linear in `psi`. Clustering, mergers and accretion that grows
   the black holes break this (Carr et al. 2026, arXiv:2601.06024v2, Sec. VI).
2. `log10 f_max` is linear between tabulated points.
3. Masses outside the table are dropped. The digitized curves stop where
   `f_max` reaches 1, so this can only lower `R`; the bound is weaker, never
   stronger, than the full curve would give.
4. One table at a time. Combining several probes in quadrature (Carr et al.
   2017, Eq. 13) needs extra assumptions and is not done here. Taking the
   lower envelope of all curves and comparing `psi` with it is wrong:
   `f_max` limits a total, while `psi` is a density per log mass.