Fit a real EHT closure triangle

Freeze the observation before fitting it

The file SR1_M87_2017_100_lo_hops_netcal_StokesI.csv contains actual low-band EHT-HOPS SR1 Stokes I observations of M87 on 10 April 2017. The channel frequency in its header is 227.0707 GHz. “Network-calibrated” identifies a released calibration stage; it does not say that every station’s atmospheric phase has become a known sky phase.

Our copy has 194,347 bytes. The script checks its SHA-256 against the pinned source and then checks that the three-row excerpt reproduces source lines 1272, 1274 and 1283 byte for byte, with the two header lines retained. All three rows have UTC 04.90972209 hours. The data repository places these measurements under PDDL 1.0; the tutorial reader and fitter were written independently.

python3 fit-eht-triangle.py --data-dir . --output reproduced-eht.json
# Compare reproduced-eht.json with eht-triangle-results.json
# Publisher CSV SHA-256:
# 1dc51e3e31a840ec34de5b5ff4d8a6af0cbb73cd8c6db85e917366ed3ec7924a
# Three-row excerpt SHA-256:
# 62ffea8dadaae26aa7fb3788c0456c4625d0e13eaad61e57529ba20311d86502

Reading: EHT 2019-D01-01 M87 data release, pinned commit 27ae76a

Read one simultaneous triangle

AA is ALMA, AZ is the Submillimeter Telescope in Arizona, and LM is the Large Millimeter Telescope in Mexico. The rows below list UTC hour, first station, second station, eastward spatial frequency, northward spatial frequency, amplitude, phase and thermal amplitude uncertainty. Spatial frequencies are already in wavelengths. Amplitudes and uncertainties are in jansky; phases are in degrees.

The publisher stores AA to AZ, AA to LM, and AZ to LM. To walk AA to AZ to LM to AA, the AA-to-LM edge must be reversed. Conjugating its complex visibility keeps the amplitude and negates the phase. The exported spatial frequencies close to within a few thousand wavelengths, tiny beside the billion-wavelength baselines; retain that rounding residual instead of editing the published numbers.

04.90972209,AA,AZ,3028965632.0000,-4496652800.0000,0.10473071,179.9880,0.00501226
04.90972209,AA,LM,2296739584.0000,-3319742720.0000,0.07996017,59.4402,0.00309308
04.90972209,AZ,LM,-732228224.0000,1176914944.0000,0.51672293,-62.7827,0.04022401

Reading: EHT 2019-D01-01 M87 data release, pinned commit 27ae76a · EHT Collaboration, First M87 Results IV, arXiv:1906.11241v1

Calculate the phase and its thermal approximation

Add the phase on AA–AZ to the phase on AZ–LM, then subtract the stored AA–LM phase. Wrapping the answer into a single cycle gives $57.7651^{\circ}$ . The script independently forms the three complex numbers and takes the argument of their directed product; the two calculations agree.

For high signal-to-noise ratio, a small transverse complex-noise perturbation changes a visibility phase by approximately noise divided by amplitude. Assuming independent thermal errors, add those phase variances. The minimum amplitude signal-to-noise ratio here is $12.85$ ; the resulting closure uncertainty is $5.68545^{\circ}$ . Amplitude-calibration errors, correlated errors and nonclosing systematics are outside this estimate.

\Phi=\operatorname{wrap}(179.9880^{\circ}-62.7827^{\circ}-59.4402^{\circ})=57.7651^{\circ}

The subtraction comes from the edge direction, not from which phase happens to be negative.

\sigma_\Phi\simeq\left[\sum_{e=1}^3\left(\frac{\sigma_e}{A_e}\right)^2\right]^{1/2}=0.09923\,\mathrm{rad}

Use radians in the uncertainty propagation and convert the final result to degrees. This is the high-SNR thermal approximation for this single triangle.

Reading: EHT Collaboration, First M87 Results IV, arXiv:1906.11241v1

Choose an empirical family that can be asymmetric

Use two circular Gaussian components with equal full width at half maximum $w$ . Put the primary at the origin and a secondary with flux ratio $r$ at displacement $(d_l,d_m)$ . Its separation is $d$ ; the angle $\psi$ is measured from east toward north. Write $s=w/\sqrt{8\ln2}$ for the Gaussian standard deviation in radians.

A normalized circular Gaussian transforms to a Gaussian in spatial frequency. Multiplying that envelope by the two-component phasor gives the model below. The total flux $F$ multiplies every amplitude and cancels from closure phase. A common translation also cancels from closure, so fixing the primary at the origin removes an unmeasurable position here.

V(u,v)=F e^{-2\pi^2s^2(u^2+v^2)}\frac{1+r e^{-2\pi i(ud_l+vd_m)}}{1+r},\qquad(d_l,d_m)=d(\cos\psi,\sin\psi)

All angular quantities inside the exponent are in radians. The family has five fitted quantities: separation, angle, flux ratio, common width and total flux. It is an empirical example, not an interpretation that M87 contains two such emitters.

Reading: NRAO Essential Radio Astronomy, Chapter 3, aperture synthesis

Solve the flux analytically, then search the shapes

For a candidate shape let $q_e$ be its unit-total-flux visibility amplitude on edge $e$ . The amplitude residual is $A_e-Fq_e$ . Divide by $\sigma_e$ , square, and add. Differentiating this sum with respect to $F$ gives a single linear equation and the weighted solution below. Thus the numerical search needs only the four shape parameters.

The script searches separations from ten to one hundred microarcseconds, five-degree angles, six flux ratios and five starting widths. It refines eight separated starting candidates by successively smaller coordinate steps. The output lists the search choices so it can be repeated exactly. This finite search demonstrates multiple good fits; it does not guarantee a global optimum.

\widehat F=\frac{\sum_e A_eq_e/\sigma_e^2}{\sum_e q_e^2/\sigma_e^2}

The denominator is positive for a nonzero candidate. This flux is fitted within that source family; it is not an independently measured zero-baseline flux.

Q=\sum_e\frac{(A_e-\widehat Fq_e)^2}{\sigma_e^2}+\frac{\operatorname{wrap}(\Phi_{\mathrm{model}}-\Phi_{\mathrm{obs}})^2}{\sigma_\Phi^2}

The wrapped phase difference avoids penalizing equivalent angles on opposite sides of the cycle boundary. The expression uses the stated approximate thermal error model.

Compare the actual saved fits

One refined fit has separation $60.16\,\mu\mathrm{as}$ , angle $240.00^{\circ}$ , flux ratio $0.596$ , width $23.56\,\mu\mathrm{as}$ and total flux $0.67307\,\mathrm{Jy}$ . It predicts amplitudes $(0.104679,0.080008,0.516147)\,\mathrm{Jy}$ in the stored row order, and closure $57.7359^{\circ}$ . Its approximate objective is $0.000575$ .

A distinctly different fit has separation $97.44\,\mu\mathrm{as}$ , angle $230.00^{\circ}$ , ratio $0.592$ and width $24.24\,\mu\mathrm{as}$ . Its objective is $0.000681$ . Both fit much more closely than the stated thermal errors require. We supplied only three amplitudes and one closure angle for five free quantities, so a tiny objective is expected to leave considerable freedom.

For comparison, imposing zero closure gives a phase-only contribution of about $103.2$ under this thermal approximation. The measured angle carries strong asymmetric-shape information, even though it leaves many asymmetric shapes available. Report the objective together with the parameter count: five free quantities for four data constraints.

Keep the real triangle and synthetic imaging separate

The real-data calculation never inverse-transforms the released raw station phases into an asserted sky image. It fits quantities that survive unknown station phases: amplitudes and a directed closure. A realistic image requires many triangles, more calibration information and a treatment of errors shared between measurements.

The interactive array instead uses a known synthetic finite-pixel sky and saves an independent ideal calibrator and total-power value. That controlled experiment lets us compare dirty images with two regularized inversions and the separate known source. The real triangle and synthetic image answer complementary questions: what the published measurements say, and how sparse sampling affects a reconstruction when the truth is known.

Reading: EHT 2019-D01-01 M87 data release, pinned commit 27ae76a · EHT Collaboration, First M87 Results IV, arXiv:1906.11241v1

What to preserve in the next triangle

Keep the release version, original rows, station ordering, timestamp definition, frequency and uncertainty columns. Derive the directed product before adding phases. Check that every edge belongs to the same integration, propagate the stated uncertainty model and compare more than one candidate shape. Preserve the predictions as well as the fitted parameters: another baseline can then test where the candidates disagree.