Recompute point-lens amplification and test a lens search

Read the exact amplification factor and the matched-filtering chapter first. Nothing here uses detector strain. All signals and noise are synthetic.

Get the files

Put these eight files in one empty folder:

• README.md, run instructions and failure cases
• point_lens_reference.py, the closed form and the two-image limit
• independent_fresnel_integral.py, the second route to $F$
• signal_injection.py, noise, chirp, injections and search
• reproduce.py, runs everything and compares
• injection-inputs.json, every injection setting and seed
• expected-outputs.json, the numbers quoted on this page
• provenance.json, sources, constants and SHA-256 checksums

Then run:

python3 -m pip install numpy scipy mpmath
python3 reproduce.py

The script prints the largest change from the saved outputs and exits with status 0 when every amplification value agrees to $2\times10^{-6}$ and every injection median to $10^{-6}$. We ran it with Python 3.9.6, numpy 2.0.2, scipy 1.13.1 and mpmath 1.3.0. Other versions may change the last digits of the random noise and fail the strict comparison; if so, python3 reproduce.py --write saves your own outputs, and the medians should stay close to ours.

Units and conventions

The solar mass parameter is the IAU 2015 nominal value $1.3271244\times10^{20}\,\mathrm{m^3\,s^{-2}}$ divided by $c^3$.

Step 1: the closed form

point_lens_reference.py evaluates Eq. (17) of Takahashi and Nakamura (2003) at 60 significant digits with mpmath, and the two-image formula, their Eq. (18). We use Eq. (17) exactly as printed. Takahashi and Nakamura 2003, Sec. 2

Step 2: a second route that never uses the closed form

independent_fresnel_integral.py starts again from the lens-plane integral. Group the points of the lens plane by their arrival time. If $A(\tau)$ is the area where $T<\tau$, then

The second form uses $\int_0^\infty e^{iw\tau}\mathrm d\tau=i/w$. $I(\tau)$ is the response the lens would give to a sharp pulse, normalized to 1 without a lens. It is the time-domain picture in Eq. (16) of the 2003 paper.

The script computes $I(\tau)$ in four steps:

1. For each direction from the lens, find where $T=\tau$ by bisection in $\ln r$, on both sides of the minimum along that ray.
2. Integrate over direction with 240-point Gauss–Legendre quadrature, after a change of variable that removes the square-root singularity where the contour closes.
3. Above the saddle time, cluster the nodes near the far side of the lens, where the contour pinches.
4. Fourier-transform the tabulated $I(\tau)$ exactly for a piecewise-linear function, then add a tail beyond $\tau=400$ fitted as $c_1/\tau+c_2/\tau^2$.

Because this route never touches the gamma or hypergeometric function, agreement in both real and imaginary parts tests the phase convention of Eq. (17) as well as its magnitude.

Result. On $y\in\{0.1,0.3,1,3\}$ and $w\in\{0.1,1,3,10,30\}$, the worst relative difference is $4.2\times10^{-4}$, at $y=0.3$, $w=30$. The fitted tail coefficient is $c_1=0.499$ for every $y$.

These four rows are a sample. All 20 are in expected-outputs.json.

The arrival-time constant matters here. If you write $T$ with $+\phi_m(y)$, as Eq. (12) of the 2003 paper is printed, the Fresnel route gives $F$ times an extra phase $e^{2iw\phi_m}$ and no longer matches Eq. (17). The paper's text says the minimum arrival time is set to zero, which needs $-\phi_m$. We use that and make no claim about the paper beyond this.

Step 3: the two limits

These tests use the closed form, now trusted from step 2.

• Low frequency. $\lvert F-1\rvert$ should fall almost in proportion to $w$, a little more slowly because of the $w\ln w$ term. At $y=1$ it is $2.2\times10^{-2}$, $3.3\times10^{-3}$ and $4.4\times10^{-4}$ for $w=10^{-2}$, $10^{-3}$ and $10^{-4}$. At $y=0.1$ and $w=10^{-4}$ it is $5.13\times10^{-4}$, and the first-order formula $\pi w/4+(iw/2)[\ln(w/2)+\gamma-2\phi_m]$ gives the same value.
• High frequency. The difference from the two-image formula, divided by $\lvert F\rvert$, should fall as $w$ grows. At $y=1$ it goes from $2.6\times10^{-3}$ at $w=30$ to $2.6\times10^{-5}$ at $w=1000$. At $y=0.3$ it is $3.1\times10^{-2}$, $6.5\times10^{-3}$, $8.7\times10^{-3}$ and $1.1\times10^{-3}$ at $w=30$, 100, 300 and 1000. The rise between 100 and 300 is real: the leftover difference oscillates with $w$.

Step 4: inject and search

signal_injection.py builds an 8 s stretch at 1024 samples per second:

1. Noise. Gaussian, coloured by Eq. (4.7) of Ajith (2011), with standard deviation $\sqrt{S_n/(4\Delta f)}$ in each real and imaginary part of each frequency bin. Ajith 2011, Sec. IV
2. Signal. A $10+10\,M_\odot$ leading-order chirp from 20 Hz to 219.9 Hz, scaled to an unlensed signal-to-noise ratio of 20, and multiplied by the exact $F$ for a lensed injection. Cutler and Flanagan 1994, Eqs. 41–42
3. Unlensed search. 17 component masses from 9.6 to $10.4\,M_\odot$, each maximized over amplitude, arrival time and phase.
4. Lensed search. 36 lens masses from 10 to $1500\,M_\odot$ times 9 values of $y$ from 0.1 to 3, then 48 finer masses between the neighbours of the best mass for each $y$, then a bounded one-dimensional polish of the best two. Lensed templates keep the true component mass.
5. Statistic. $\Delta\ln L$, the best lensed minus the best unlensed log likelihood ratio. The reporting threshold is 4.

Before any noise is added, the script checks that the search finds the injected lens. It does: lens mass within 0.001 percent and $\Delta\ln L$ within $10^{-8}$ of the exact-template value for both lensed cases.

Result.

The clear lens is unmistakable. The faint lens is a coin toss, and the noise-only tail reaches values a faint lens also produces.

Failure cases kept on purpose

• Coarse grid. The first version searched only the 36-mass grid, about 15 percent apart. At $y=2$ the 90 ms delay shifts the fringe phase at 200 Hz by several radians between neighbouring masses, so even without noise the best grid template reached $\Delta\ln L=0.73$ instead of 6.03. Its faint-lens median, 2.79, was mostly noise. The refinement in step 4 fixes this, and the noise-free check would catch it again.
• Coarse amplification table. The first table of $F$ used 1500 points in $w$. At $y=3$ and $w$ near 41 it sampled each fringe about four times, and interpolation errors reached 3.8 percent. 6000 points bring this to 0.24 percent.
• False lenses. 14 of 120 noise-only trials pass the threshold. The largest, seed 3058 with 7.97, prefers $M_{Lz}\approx1150$ and $y=1.5$.
• What is missing. No spins, unequal masses, merger, ringdown, detector antenna pattern, glitches or galaxy around the lens; $y$ is searched only on its grid. These numbers illustrate the method and are not the sensitivity of any real search.

Checklist before you trust a rerun

• python3 reproduce.py exits with status 0, or your own --write outputs differ from ours only in noise-dependent digits.
• The SHA-256 checksums in provenance.json match your files.
• The Fresnel route and Eq. (17) agree in real and imaginary parts, not just in magnitude.
• $\lvert F-1\rvert$ falls about as fast as $w\ln w$ at low frequency, and the two-image error falls overall at high frequency.
• The noise-free search recovers the injected lens before you read any noisy result.
• The noise-only distribution is reported next to every lensed result, with the same threshold.