Dirty beams and nonunique images

Sum the measured waves back toward each pixel

Let $N$ be the number of distinct complex samples, $y_k$ the calibrated visibility and $(l,m)$ a trial sky position. Multiply each measured signal by the opposite phase and average its real part. If the trial position matches a bright point, its contributions line up. Elsewhere they partly cancel. The result is the dirty image $D$ .

Set every visibility to one and do the same sum. That gives the dirty beam, the response to a unit point at the origin. Substituting the forward model into the dirty-image sum shows that every true point is replaced by a shifted copy of this beam. Gaps between sampled frequencies create sidelobes, including negative ones.

D(l,m)=\frac1N\sum_k\operatorname{Re}\!\left[y_k e^{2\pi i(u_kl+v_km)}\right],\qquad B(l,m)=\frac1N\sum_k\cos[2\pi(u_kl+v_km)]

Uniform weighting and the factor one over N give B at the origin equal to one. Explicitly adding the conjugate samples gives the same real result after accounting for their doubled weight.

D(l,m)=\sum_p x_p B(l-l_p,m-m_p)

This finite-source sum explains the name dirty image. Each source pixel contributes its flux times a shifted dirty beam.

Reading: NRAO Essential Radio Astronomy, Chapter 3, aperture synthesis · EHT Collaboration, First M87 Results IV, arXiv:1906.11241v1

A beam value and a pixel value mean different things

A central $0.60\,\mathrm{Jy}$ point has a dirty-image peak of $0.60\,\mathrm{Jy}$ per unit-peak beam. Summing all dirty-image pixels generally does not recover $0.60\,\mathrm{Jy}$ . The sidelobes and chosen image window affect that sum. A reconstructed pixel, in contrast, stores its own integrated flux in jansky, and those pixel fluxes are constrained to sum to the measured total.

The default even-sized grid has its origin between four pixels. Its sampled beam maximum can therefore be slightly below one although the continuous beam has $B(0,0)=1$ . Every displayed map uses its own peak for color stretch and prints that peak. Read each map’s printed peak to translate its colors into the stated physical units.

Build two genuinely different skies with identical signals

Suppose every baseline has $v=0$ . Move a point north without changing its east coordinate or flux. Every exponent in the forward sum stays the same. Two positive images at different north positions then produce exactly the same sampled signals, even with perfect noiseless measurements. Longer east–west baselines do not resolve this particular ambiguity.

For the default $24\times24$ image there are $576$ real unknown pixel fluxes. Fifty-four complex measurements provide at most $108$ real linear equations. An independently known total flux adds one. Positivity can narrow the possibilities, but sparse sampling has left a large space of directions in which the data change little or not at all.

Extended emission is particularly easy to miss if no short baselines sample it. A broad smooth component can contain substantial total flux yet vary too slowly to contribute much at long-baseline frequencies. Supplying a total flux helps; the reconstruction still has to distribute it among the uncertain low-frequency structures.

Reading: NRAO Essential Radio Astronomy, Chapter 3, aperture synthesis · EHT Collaboration, First M87 Results IV, arXiv:1906.11241v1

State what each reconstruction prefers

Write the discrete Fourier operator as $A$ and the pixel flux vector as $x$ . We minimize the mean squared complex signal residual plus a regularizer $R$ , subject to nonnegative pixels and the separately measured total flux. The regularizer is an explicit preference used to choose an image when measurements alone leave choices.

The smooth option penalizes horizontal and vertical neighbor differences. Each neighbor pair is counted once. The compact option adds a radius-weighted flux penalty and a small squared-flux penalty. Here $\rho_p$ is radius in units of the full field width, measured from the center. The numerical coefficients below are exactly those used by the browser, with fluxes expressed in jansky.

J(x)=\frac1{2N}\sum_k|(Ax)_k-y_k|^2+R(x),\qquad x_p\geq0,\quad\sum_px_p=F_{\mathrm{total}}

The data term has units of jansky squared. The saved synthetic signals contain no added random noise; the reference uncertainty is five millijansky per complex sample and is not used as an extra weighting factor in this objective.

R_{\mathrm{smooth}}=\frac{0.05}{2}\sum_{\langle p,q\rangle}(x_p-x_q)^2

Neighbor differences make rapid pixel-to-pixel changes expensive. A broad structure can be preferred even if finer structure also predicts the observations.

R_{\mathrm{compact}}=(0.004\,\mathrm{Jy})\sum_p\rho_p^2x_p+\frac{0.03}{2}\sum_px_p^2

The first term favors emission near the field center. The second discourages putting all flux in very few pixels. Together these differ from the neighbor-smoothness preference.

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

Derive a step that preserves positivity and total flux

For a real pixel perturbation, differentiating the squared complex residual gives $\operatorname{Re}[A^*(Ax-y)]/N$ . Add the regularizer gradient, take a trial downhill step and project it onto the positive fixed-flux simplex. The projection subtracts one common threshold, clips negative entries to zero and chooses the threshold so that the remaining flux sums to $F_{\mathrm{total}}$ .

Why one threshold? Minimize the squared distance between a trial vector and the allowed image. Differentiating with a multiplier for total flux subtracts the same multiplier from every free pixel. Pixels pushed below zero sit on the nonnegativity boundary. Sorting the trial values identifies which pixels remain free and therefore determines the threshold.

The implementation takes ninety accelerated projected-gradient steps. Backtracking checks that a quadratic local model bounds the trial objective; momentum restarts when the full objective increases. This is a finite numerical calculation, so compare the objective history and predicted signal residuals as well as the image. Increasing iteration count and changing the prior answer different questions.

x^{\mathrm{next}}=\Pi\!\left(x-\eta\nabla J(x)\right),\qquad [\Pi(z)]_p=\max(0,z_p-\tau),\quad \sum_p\max(0,z_p-\tau)=F_{\mathrm{total}}

The projection enforces both constraints together. Clipping negatives and then renormalizing the positives is generally a different operation.

Read the saved default before changing anything

The saved default has fifty-four samples, total flux $0.600\,\mathrm{Jy}$ , a three-hour track and the two-component finite-pixel source. After ninety steps the smooth reconstruction has complex-signal RMS residual $7.041\,\mathrm{mJy}$ ; the compact preference gives $5.934\,\mathrm{mJy}$ . The RMS is the square root of the mean squared complex modulus, not a per-pixel image error.

These nonzero residuals come from the finite optimization and the stated priors acting on noiseless data. They are not estimated astrophysical noise. Both inversions read the saved calibrated signals, grid and independent total flux. The source reference is passed only to its separately labeled display. Exporting the observation preserves those inputs so the calculation can be rerun after the source generator has been removed.

\mathrm{RMS}_{V}=\sqrt{\frac1N\sum_k|(Ax)_k-y_k|^2}

A lower signal residual alone does not prove that every image feature is correct. Compare alternative priors and test how a known synthetic source is recovered under the same sampling.

Image grid and field	$24\times24,\quad120\,\mu\mathrm{as}$
Central Greenwich hour angle	$H=105^{\circ}$
Total flux	$F=0.600\,\mathrm{Jy}$
Recorded complex samples	54
First directed closure phase	$0.287631^{circ}$
Smooth-prior signal RMS	$7.0408,mathrm{mJy}$
Compact-prior signal RMS	$5.9340,mathrm{mJy}$