On the variance of radio interferometric calibration solutions: Quality-based Weighting Schemes

This paper investigates the possibility of improving radio interferometric images using an algorithm inspired by an optical method known as"lucky imaging", which would give more weight to the best-calibrated visibilities used to make a given image. A fundamental relationship between the statistics of interferometric calibration solution residuals and those of the image-plane pixels is derived in this paper. This relationship allows us to understand and describe the statistical properties of the residual image. In this framework, the noise-map can be described as the Fourier transform of the covariance between residual visibilities in a new differential Fourier plane. Image-plane artefacts can be seen as one realisation of the pixel covariance distribution, which can be estimated from the antenna gain statistics. Based on this relationship, we propose a means of improving images made with calibrated visibilities using weighting schemes. This improvement would occur after calibration, but before imaging - it is thus ideally used between major iterations of self-calibration loops. Applying the weighting scheme to simulated data improves the noise level in the final image at negligible computational cost.


Introduction
Interferometers sample Fourier modes of the sky brightness distribution corrupted by instrumental and atmospheric effects rather than measuring the sky brightness directly. This introduces two problems for astronomers to invert: calibration and imaging. Both of these problems are ill-conditioned.
The problem of imaging consists of correcting for the incomplete uv-coverage of any given interferometer by deconvolving the instrument's Point-Spread Function (PSF) from images. Its poor conditioning comes from our limited a priori knowledge of the sky brightness distribution, combined with large gaps in our uv-coverage, which prevents us from placing strong constraints on image deconvolution. It can be better-conditioned in different ways, including through the use of weighting schemes (see Briggs 1995;Yatawatta 2014, and references therein) to improve image fidelity at the start of deconvolution. When inverting the imaging problem, we often assume that the sky is stable within the domain (i.e. is constant in time and frequency). There are exceptions, such as wide-band deconvolution algorithms (e.g. Rau & Cornwell 2011) that explicitly take into account the sky's frequency-dependence, but still assume that the sky brightness distribution does not vary with time.
The problem of calibration is what concerns us in this paper. It consists of estimating and correcting for instrumental errors (which includes effects such as antenna pointing errors, but also the phase-delays caused by ionospheric activity, troposphere, etc). Calibration consists of solving for gain estimates, where a gain models the relationship between the electromagnetic field of an astrophysical source and the voltage that an antenna measures for this source. Because measurements are noisy, calibration often involves some fine-tuning of solution intervals, to ensure that the solutions are well-constrained while the solution intervals stay as small as signal-to-noise allows. The calibration inverse problem involves three competing statistical effects: thermal noise in the measurements, true gain variability, and sky model incompleteness. If gain solutions are sought individually for each measurement, then calibration estimates will be dominated by thermal noise, and will not adequately describe the actual gains. Similarly, if a single gain estimate is fitted to too many measurements, the intrinsic gain variability will be "av-eraged out"; for example, a choice of time and frequency interval that is too large will cause the solver to estimate a constant gain while the underlying function varies quickly, thereby missing much of the gain structure. This will introduce error which will be correlated in time and frequency. This occurs, for example, when solving for ionospheric phase delays: in the most extreme case, where the solution interval is significantly larger than the scale of ionospheric fluctuations, its varying phase can average out to zero over the interval in time and frequency. Finally, if the model being fitted is incomplete, unmodeled physical flux will likely be absorbed unpredictably into both the gain solutions and the residual visibilities: this absorption of physical flux into gain solutions is known as source suppression (see Grobler et al. 2014;Kazemi & Yatawatta 2013, and references therein).
In practice, it is reasonable to assume that gain variation is generally slower than some given scale: we can then reduce the noise of our gain estimates by finding a single gain solution for a small number of measurements, assuming that the underlying gain variation is very small and stable over short intervals. This is generally a valid hypothesis, but the specific value for the variation scale can be contentious. Indeed, while the noise level can often be treated as constant throughout an observation, the gain variability itself is generally not constant: there will be time periods where the gains will tend to remain constant for longer, and others where variability will be very quick. This means that, for any choice of calibration interval, some gain estimates will be better than others, and almost all could have been improved (at a cost to others) by a different choice of time (and frequency) intervals.
Since we have measurements which are bettercalibrated than others (in that better estimates for their gains were obtained through chance alone), we could, in principle, take inspiration from "lucky imaging" (an optical-domain method for making good images: for more details, see Fried 1978, and references therein) to weigh our visibilities according to their calibration quality. Those weights would in effect be an improvement of currently existing methods such as clipping noisy residual visibilities: in the extreme case where all visibilities are equallywell calibrated except a few which are extremely noisy, it should be equivalent to clipping. Otherwise, the weights should show at least a slight improvement over clipping.
The key finding of the present paper is a fundamental relationship between the covariance matrices of residual visibilities and the map of the covariance in the imageplane: the "Cov-Cov relationship" between visibility covariance to image-plane covariance. We show that the pixel statistics in the image-plane are determined by a "noise-PSF", convolved with each source in the sky (modeled or not). This noise-PSF is the product of the Fourier transforms of the gain covariance matrix with each cell mapped not from uv space to lm coordinates but rather between their respective covariance spaces -from a new differential Fourier plane (henceforth "(δuδv)"-plane) to the image-plane covariance space δlδm. This image-plane covariance space describes the variance in each pixel and the covariance between pixels 1 . It describes the expected calibration artefacts and thermal noise around each source, does not vary as a function of direction, and is convolved with each source in the field to yield the final error map. Because all unwanted (in our case, unphysical) signal can be thought of as noise, we will refer to the pixel variance map as the "noise-map".
The notion of a (δuδv)-plane arises organically from the framework of radio interferometry: we are associating a correlation between visibilities to coordinates in covariance space, just as we associate the visibilities themselves to the uv-domain. The (δuδv)-plane is the natural domain of these correlations. As previously stated, even if all sources in the field are perfectly known and modeled, a poor choice of calibration interval can introduce correlated noise in the residuals, which would then introduce larger variance near sources in the noise-map. Conversely, if calibration is perfect, the noise-map should be completely flat (i.e. same variance for all pixels), as there would be no noise-correlation between pixels.
The main result of this paper consists of describing a new adaptive, quality-based weighting scheme based on this insight. Using the Cov-Cov relationship, we can create a new weighting scheme by estimating the residual visibility covariance matrix in a given observation. By weighting visibilities so as to change their covariance matrix, one can change the shape of the noise-PSF and thus improve the final image: this manifests as either decreased noise or decreased calibration artefacts. Note that this weighting is applied after calibration, but before image deconvolution: applying it will therefore not only improve the residual noise in the image, and thus the sensitivity achievable with a given pipeline, but will also improve deconvolution by minimising calibration artefacts in the field: it should thus effectively remove spurious, unphysical emission from final data products. Estimating the covariance matrix is the main difficulty of our framework: we do not know the underlying covariance matrix, and the conditioning of our estimation thereof is limited by the number of measurements within each solution interval. As such, we have no guarantee that our estimate of the corrected visibility covariance matrix is accurate. This problem can be alleviated, for example by estimating the covariance matrix for the antenna gains themselves, and use it to build the visibility covariance matrix: this effectively improves conditioning (cf Sec. 4).
This paper is split into four main sections. In Section 2, we derive the Cov-Cov relationship. With its newfound insights, we propose quality-based weighting schemes with which to improve radio interferometric images in Section 3. We follow in Section 4 by showing how to estimate, from real data, the covariance matrix from which the quality- Baseline selection matrix, which picks out 1 visibility out of the full set. Size npix × n b C b npix × npix convolution kernel that defines the PSF. F bb ′ Convolution matrix mapping one δuδv to δlδm. The set of all F bb ′ determines the noise-PSF. Size npix × npix. based weights are derived. Our approach seems to give good results. Finally, we close the paper on a discussion of the applicability and limitations of the quality-based weighting scheme.

Building the Noise Map
In this section, we derive our first fundamental result: the Cov-Cov relationship, Eq. 25, which describes how the statistics of residual visibilities (and thus the antenna calibration solutions, henceforth "gains") relate to the statistics of the image plane, i.e. of images made using the associated visibilities. The dimensions of the matrices (denoted by boldface capital letters) and vectors (denoted by boldface lowercase letters) used in this paper are given in Table 1, along with the scalar numbers used to denote specific dimensions. All other variables are scalars.

The Cov-Cov Relationship in the δuδv plane
Let us begin by defining visibility gains. Using the Radio Interferometry Measurement Equation formalism for a sky consisting of a single point source (Hamaker et al. (1996), Smirnov (2011), and companion papers), we can write the following relation between the sky and the signal as measured by a single baseline at time t and frequency ν: All the quantities above are 2 × 2 matrices. Eq. 1, implies a linear relationship between the coherency matrix B d ν and the visibilities recorded by a given baseline (V tν pq ), with the addition of a thermal noise matrix N , which is also of shape 2 × 2 and contains different complex-valued realisations of the noise in each cell. Since electric fields are additive, the sky coherency matrix can be described as the sum of the contributions from individual sources in directions d in the sky. We also assume that the sky does not vary over time, i.e. that B d ν is not a function of time. The Jones matrices (J d ...,tν ) contain the antenna gain information in matrix form, while K d ...,tν is the Fourier kernel. Let us limit ourselves to the scalar case, which corresponds to assuming that emission is unpolarised. We assume that B d = s d I, where s is the flux of our single point source and I is the 2 × 2 identity matrix. We also assume that J p,tν = g p,tν I, where g tν p is the complex-valued gains of antenna p at time t and frequency ν. This means that we assume that the gains are direction-independent, and so J d ...,tν becomes J ...,tν . Similarly, K d p,tν = k d p,tν I, the Fourier kernel in the direction of the source, d. N has 1 realisation of in each cell, where: where σ is the variance of the thermal noise. Let us denote each (t, ν) pair by τ , and ignore the sky's frequencydependence. The following scalar formulation is then equivalent to Eq. 1: Calibration is the process of finding an accurate estimate of g τ p for all antennas p, at all times t and frequencies ν. Since we are in a direction-independent regime, the quality of our calibration then determines the statistical properties of the residual visibilities (and the image-plane equivalent, the residual image). The residual visibilities associated with calibration solutions are defined as our measured visibilities minus the gain-corrupted model visibilities.ĝ τ p then denotes our calibration estimate for g τ p .
We now begin to limit the generality of our framework by assuming that all sufficiently bright sources have been modeled and subtracted: unmodeled flux is then negligible. We can then write the residual visibilities as: The flux values in the image-plane pixels 2 are the Fourier transform of the visibility values mapped onto each pixel. This can be written as follows: where lm are the directional cosine positions of a given pixel, and k lm pq,τ = k d p,τ k d q,τ the Fourier coefficient mapping a point in Fourier space to a point on the image-plane. ω pq,τ is the weight associated to a given visibility. Let us now write this using a matrix formalism. The contribution of a single visibility b = (pq, τ ) to the image-plane residuals can be written as: where is a vector of the of Eq. 2 andỹ is a vector of size n pix , with and w b is the scalar weight associated with each visibility. By default, w b = 1 n b : all visibilities then have the same weight, andỹ then becomes the average of allỹ b .γ is a vector of allγ b , and thus of size n b . F is the Fourier kernel, of size n pix × n pix . S b is a matrix of size n pix × n b : its purpose is to encode the uv-coverage. Each S b contains only a single non-zero cell, different for different S b . The height (number of rows) of S b is determined the size of the uv-grid, and its length (number of columns) by the number of visibilities.
The order of operations is thus: each residual visibility (κ bγ + n) is assigned some weight w b and its uvcoordinates are set by S b . The inverse Fourier transform (F H ) is then applied to this grid, and so we recover its image-plane fringe. By averaging over all fringes, we recover the dirty image.
The residual image will thus depend on three quantities: the residual gains, the flux in the image, and the weighting scheme. Let us consider the relationship between the statistics of residual visibilities and the variance at a given point in the corresponding residual image.

Statistical Analysis
In the following analysis, we treat our gain solutions and thermal noise as random variables in order to compute the covariance matrix of our residual image, Cov{ỹ}. The diagonal of this matrix gives the variance for each pixel, while the wings give the covariance between pixels. Using the property that Cov{Ax} = ACov{x}A H , we can apply the Cov{} operator to Eq. 9 to write: So far, we have only applied definitions. The net effect of is to encode where a given baseline samples the uv-plane, and map one cell at matrix coordinates (b, b ′ ) from the correlation matrix Cov{γ} onto the visibility grid. S b is not the gridding kernel, but rather the sampling matrix, which determines where we have measurements and where we do not. We can thus write that bb ′ is the value from the appropriate cell and 1 is the vector-of-ones of appropriate length (here, n b ). This allows us to write: Here, C b is a Toeplitz matrix, i.e. a convolution matrix, associated with baseline b. The set of all C b defines the convolution kernel which characterises the Point-Spread Function (henceforth PSF) associated with a given uvcoverage, of size n pix × n pix . F bb ′ , meanwhile, is not generally Toeplitz. Its cells can be written as: Let us investigate how the sky brightness distribution (i.e. d-dependence) affects the noise-map. We can write the sum over bb ′ as two sums: one over b, b ′ = b and one over b, b ′ ≠ b. Thus: Note that the only direction-dependent terms in the above are s d and κ d b , which are both inside φ bb ′ (for both b = b ′ and b ≠ b ′ ). By making the approximation that the Fourier kernels of different sources are orthogonal (i.e. that κ d Note that F bb = C b , since those are the coordinates along the diagonal: for these values, the matrix-of-ones at the centre of F bb ′ becomes the identity matrix. Note also that We can then write Eq. 20 as: where we have now limited our formalism to the case where the sky is dominated by distant point-like sources. This is our fundamental result: assuming unpolarised emission coming from distant point sources and normallydistributed thermal noise, it gives a direct relationship between the covariance of the residual visibilities and the covariance of the residual image-pixel values. We thus call it the Cov-Cov relationship. It describes the statistical properties of the image-plane as the result of a convolution process changing an average noise level at different points in the image-plane, allowing us to describe the behaviour of variance and covariance in the image. Let us focus on the first. By applying the Diag {} operator (which returns the diagonal of an input matrix as a vector) to both sides of Eq.25, we can find an expression for the variance map in the image-plane: In Eq 27, we have: For b ≠ b ′ , the diagonals of F bb ′ are the Fourier kernels mapping δuδv space to δlδm. F bb ′ can then be thought of as a Fourier transform. It is not a diagonal matrix. It behaves as a covariance fringe, allowing us to extend standard interferometric ideas to covariance space: each fringe can be thought of as a single "spatial filter" applied to the pixel covariance matrix. Just as a given baseline has coordinates in uv-space, a given correlation between baseline residual errors has coordinates in uv correlation space, which we will henceforth refer to as δuδv-space. This δuδv space warrants further discussion: Fig. 1 shows, for a given uv-track (Fig. 1b), both the corresponding δuδv domain (Fig. 1c) and point-spread function (Fig.  1a). The symmetric, negative uv-track is treated as a separate track, and thus ignored. This means that we do not fully constrain the noise-PSF (since the covariance matrix of the symmetric track is simply the Hermitian of the first), but we do not seek to constrain it in this section, but rather to show that our results hold.. We can see that the δuδv-tracks are symmetrical about the origin. The δuδv space corresponding to a given uv-track can thus be most concisely described as a "filled uv-track", with its boundaries defined by the ends of the uv-track. The set of F bb ′ , each of which maps one value of the covariance matrix to a fringe in the image-plane, would then characterise a PSF equivalent for the noise distribution, which we refer to as the noise-PSF. In our formalism, the only source of statistical effects in the field are calibration errors and thermal noise. The average variance in all pixels will be given by the diagonal of the covariance matrix and the thermal noise, provided that it truly follows a normal distribution. The only effects which will cause the variance in the imageplane to vary from one pixel to the next are those mapped onto the covariance fringes, i.e. such position-dependent variance fluctuations will be caused by correlated gain errors, which are spurious signal introduced by erroneous gain estimates. Assuming all sources in the field are pointlike and distant, then these variance fluctuations will follow a specific distribution, convolved to every source in the field. Since the variance fluctuations act as tracers for calibration artefacts, artefacts in the image can be understood as one realisation of the variance map, which is characterised by an average level determined by the variance in the gains and thermal noise, and a noise-PSF convolved with the sky brightness distribution. The actual artefacts in the image will still be noisy, as a realisation of the true variance map. For the same reason, in the absence of correlated gain errors, [Cov{γ}] b≠b ′ are all zero and Cov{γ} is a diagonal matrix. We then recover a "flat" noise-map: the variance will be the same in all pixels, as the noise-PSF is absent. In the ideal case, were we to recover the true value of the gains for all times and frequencies, this becomes pure thermal noise.

Noise Map Simulations
We have shown in Eq. 25 that there exists an analytical relationship between residual visibility statistics and image-plane residual statistics. This section gives details of simulations we have performed to support our claims on this "Cov-Cov" relationship. Specifically, we simulate residual visibilities for a single baseline by generating a set of correlated random numbers with zero mean and a distribution following a specified covariance matrix C. It contains a periodic function of period T along the diagonal, which is then convolved with a Gaussian of width σ τ corresponding to the characteristic scale of correlation. The values of these parameters are chosen arbitrarily. A small constant term is added on the diagonal, the net value of which is strictly positive. This simulates a low thermal Fig. 2: Example of a non-stationary covariance matrix, which can be used to simulate Cov{γ}.The colourbar units are Jy 2 . The correlation scale σ τ is 40 cells, and the variability period is 500 cells. The matrix is made positive semi-definite (and therefore a covariance matrix) through SVD decomposition. The maximum size of the "bubbles" is determined by σ τ .
noise. Finally, singular value decomposition is used to ensure that this matrix is Hermitian positive semi-definite. The net effect is a non-stationary correlation: some residuals are correlated with their neighbours, and uncorrelated with others. An example of this covariance matrix for arbitrary parameter values is shown in Fig. 2. We see that, for any given point, correlation is stronger with some neighbours than others (as determined by σ τ ). Samples are drawn as follows: C is built following user specifications as described above. We then find its matrix square root C 0 so as to apply it to random numbers generated from a normal distribution. We generate 2000 realisations i of our random variables r i : and x ←N (0, 1) whereỹ is a matrix of dimensions n realis × n b . Since x follows a normal distribution, Cov{x} = I and the covariance matrix of eachỹ i is, by construction, C = C 0 C H 0 . The covariance matrix ofỹ is thus also C.
As for the uv-track, our simulations read a single one from a specified dataset. In this case, we read an 8-hour uv-track for a baseline between two arbitrary LOFAR stations (specifically, CS001HBA0 and RS310HBA) in an observation of the Bootes extragalactic field. The effective baseline length varies between 37.9km and 51.8km. The dataset included 20 channels, each with a spectral width of 97.7 kHz; the central observing frequency is 139 MHz. The temporal resolution is 1 measurement per second.
We compare the measured variance map V m y , built by measuring the variance across realisations at each pixel in the image-plane, with the predicted variance map Var{ỹ}, built using the Cov-Cov relationship (Eq. 25). Since we are only interested in the variance map, rather than the covariance between pixels, we compute only the diagonal terms.
where the thermal noise is already incorporated into the diagonal of C and Eq. 36 is merely the diagonal operator applied to the Cov-Cov relationship.

Simulation with a single point source
We model our sky as containing a single 1 Jy point source at phase centre: we thus have φ bb = w 2 b . The source as seen through the set of uv-tracks used in our simulation, along with their corresponding (δu, δv) space, are shown in Fig. 1. The simulated noise-map is calculated by drawing a large sample (n realis = 2000) of random numbers from the correlated distribution, thereby creating 2000 sets of residual visibilities. By Fourier-transforming the visibilities to the image-plane and taking the variance of the values for each image pixel (i.e. each l, m pair) as per Eq. 35, we can estimate Var{ỹ}. The predicted noise-map, meanwhile, was found by assigning each cell of C to the appropriate point in the (δu, δv) plane and Fourier transforming from this plane into the image-plane, as per Eq. 36. We compare the outcome of simulating a large number n realis of realisations and taking the variance across these realisations for each pixel with mapping the covariance matrix onto the δuδv-plane and using the Cov-Cov relationship.The results of our simulations are shown sideby-side in Fig. 4: the predicted and simulated noise-PSFs match. The peak-normalised predicted noise-PSF is less noisy, as shown in Fig. 3 for different correlation scales. This is expected, since it is calculated directly from the underlying distribution, rather than an estimate thereof. As n realis → ∞, we expect the two methods to fully converge. As the maximum characteristic correlation length σ τ increases, the variance becomes ever more sharply peaked.
Since our simulated sky consists of a 1Jy source at phase centre, there is only one noise-PSF to modulate the average noise-map, and it lies at phase centre. Let us test our formalism further by considering a model with multiple point sources.

Simulation with 3 point sources
We wish to test our prediction that the noise-map can be described as a convolutional process modulating an average noise level. We thus perform another simulation, this time with three 1Jy point sources. The associated dirty image is shown in Fig. 5d. This dirty image simply consists of performing a direct Fourier transform (i.e. without using a Fast Fourier Transform algorithm) on simulated visibilities corresponding to these three point sources. We now perform a similar test as above on this field. Firstly, we "apply" gain errors to these visibilities by multiplying our model with our residual gain errors. This allows us to find 2000 realisations of residual visibilities, and find the variance for each pixel across these realisations. This gives us the simulated noise-map, shown in Fig. 5a. Secondly, we perform a DFT from the differential Fourier plane to the (l, m) plane as before, assigning one cell of Cov{γ} to each point of the differential Fourier plane. This time, however, φ d bb ′ is not simply unity for all points in the differential Fourier plane. Instead, it is calculated for the three-point-source model, and applied for each point. This gives us the predicted noise-map in Fig. 5b. Finally, Fig. 5c shows the absolute value of the difference between the two images. We see that there is some structure present in these residuals: this is expected, as the PSF of the sources in the dirty imags clearly overlap. We are thus not quite in the regime where emission is fully spatially incoherent. Nevertheless, our predictions hold to better than 5% accuracy.
It bears repeating that, for correlated noise, this map can be understood as a distribution map for calibration artefacts: the amount of spurious correlated emission seen by each baseline will determine the noise-map, and the true image-plane artefacts will then be one set of realisations of this underlying distribution.

Adaptive Quality-based Weighting Schemes
As discussed in Section 1, some intervals of an observation will have lower gain variability. These will show up in the gain covariance matrix as intervals with lower variance. Similarly, those with larger intrinsic gain variability will have greater error in their gain estimate. By giving greater weights to the former, and lower weights to the latter, we expect to be able to improve image reconstruction. We thus talk of adaptive quality-based weighting, as the weights will adapt based on the calibration quality.
The pixel variance is determined by the visibility covariance matrix, as shown in Eq. 27. The diagonal of the visibility covariance matrix will add a flat noise to all pixels, while its wings will determine the calibration artefact distribution, which will be convolved to the sky brightness distribution. We thus have two sources of variance in the image-plane. Minimising the far-field noise (i.e. the variance far from sources) in an image would involve downweighting noisier calibration intervals while up-weighting the more quiescent ones, without taking noise-correlation between visibilities into account. This is because the farfield noise will be dominated by the diagonal component of the covariance matrix (cf. Eq. 27). By the same token, minimising calibration artefacts would involve downweighting measurements with strongly-correlated noise, and up-weighting the less-correlated. This would not, however, minimise the diagonal component: in fact, it will  likely exaggerate its up-weighting and down-weighting. As such, it will increase the constant level of the noisemap, but flatten the noise-PSF's contribution. There are thus two competing types of noise that we seek to minimise: uncorrelated noise, which corresponds to δuδv = 0 (i.e. the diagonal components of the gain covariance matrix), and correlated noise, which corresponds to δuδv ≠ 0 (i.e. its wings). Minimising the first will minimise far-field noise without optimally reducing artefacts, while minimising the last will minimise noise near sources at a cost to far-field noise. In the following sub-sections, we will discuss weighting schemes used to accomplish this.

Optimising sensitivity
The Cov-Cov relationship (Eq. 25) tells us that, far from any sources, the variance map (Eq. 27) is dominated by a constant term: the contribution from thermal noise and the diagonal of the residual visibility covariance matrix. Maximising sensitivity far from sources therefore implies minimising Diag {Cov{γ}}. This is equivalent to assigning visibilities weights inversely proportional to their variance: For each baseline, those times with larger variance in the residuals will be down-weighted, and those with smaller variance will be up-weighted; this scheme does not require information about the underlying gains, only the error on our solutions. Since we are treating σ 2 n as a constant for all antennas and all times, those times where our gains estimate is closer to the true gains will be upweighted, and those moments where they are farther from the actual gains will be down-weighted: hence the term "adaptive quality-based weighting". Note that the diagonal of the weighted residuals' covariance matrix should therefore become constant: this weighting scheme explicitly brings the residuals closer to what is expected in the case of perfect calibration, assuming uncorrelated noise. For the remainder of this paper, we will refer to these weights as sensitivity-optimal weighting.

Minimising Calibration Artefacts
Minimising calibration artefacts -i.e. optimising the sensitivity near bright sources -means flattening the noisemap. Since the noise-map can be understood as a noise-PSF convolved with all the modeled sources in the sky modulating the background variance level, it will be flattest when its peak is minimised. From the Cov-Cov relationship (Eq. 25), we can see that, at the peak of the noise-PSF (which would be the variance at the pixel where l = m = 0), the Fourier kernel is unity: the variance for that pixel is thus the sum of all the cells in the covariance matrix. By accounting for normalisation, we can write the variance at the centre of the noise-PSF as: Our optimality condition is then, after some algebra: (40) We find that one w which satisfies the above is: where 1 is a vector of ones. These weights depend only on calibration quality: badly-calibrated cells will include spurious time-correlated signal introduced by trying to fit the noise n on visibilities. Down-weighting these cells helps suppress artefacts in the field, at the cost of far-field sensitivity. This weighting scheme is thus only a function of the relative quality of calibration at different times, boosting better-calibrated visibilities and suppressing poorlycalibrated visibilities. For the remainder of this paper, we will refer to these weights as artefact-optimal weighting.

Estimating the Covariance Matrix
In our simulations, we have worked from a known covariance matrix and shown that our predictions for the residual image's behaviour hold. With real data, however, we do not have access to this underlying covariance matrix. Since our weights are extracted from said matrix, estimating it as accurately as possible remains a challenge: this is in turn limited by the number of samples which can be used for each cell. Each cell in the covariance matrix is built by averaging a number of measurements, or samples. The more samples are available, the better our estimate becomes: once we have more samples than degrees of freedom, we say that our estimation is well-conditioned. Otherwise, it is poorlyconditioned. In this section, we will discuss ways in which we can improve the conditioning of the covariance matrix estimation.

Baseline-based Estimation
One way to improve the conditioning of our covariance matrix estimation is to make the same hypothesis as the calibration algorithm: we can treat the underlying gains as constant within each calibration interval. Provided this interval is known, this allows us to find a single estimate for each interval block of the covariance matrix, turning a n b ×n b matrix into a smaller n intervals ×n intervals equivalent, where n intervals is the number of solution intervals used for to find the gain solutions. We then improve our conditioning by a factor of n int , which is the number of samples in a calibration interval. The estimate Cov{γ}of the covariance matrix Cov{γ} is built by applying the covariance operator: where the ⟨⋯⟩ operator denotes taking the average over the full vector. If the calibration solver's gain estimates are unbiased (i.e. E{ĝ} = g) and the model of the sky is sufficiently complete, this quantity should be zero. Having createdĈγ, which will be of size n b × n b , its cells can now be averaged over blocks of n int × n int . This allows us to estimate the weights for each baseline and each time.
Mathematically, we retrace the steps of Section 2. In the absence of direction-dependent effects, we define the residual visibilities as before, and use them to define the normalised residual visibilities ρ b : We then organise the residuals in cells: R corresponds to a matrix containing all the residual visibilities within one calibration cell C, i.e. for b ∈ C where g C = const. It is therefore of size n intervals ×n C , where n intervals is the number of calibration intervals in the observation. Normalising the residual visibilities by k b allows us to recover the underlying covariance matrix by multiplying the residual visibility matrix R with its Hermitian conjugate:Ĉγ Note that we have divided the noise term by the flux model S b , which can be very small in some cells. As such, care must be taken not to cause the relative thermal noise contribution to explode: those cells where this would occur are dominated by thermal noise, and information on the covariance matrix cannot be recovered from them. In this framework, we simply treat the index C as containing all the times and frequencies, for individual baselines, corresponding to a single calibration interval.Ĉγ is then an estimate of the residual visibility covariance matrix.

Antenna-based Estimation
In the subsection above, we assumed that finding one solution per interval will give us strong enough constraints to make the problem of estimating the covariance matrix well-conditioned: this may not be true in all cases. Conditioning may then need to be improved further: in this subsection, we show one way in which this can be done. There are others, e.g. using the rank of the matrix itself to find better-conditioned estimates of the covariance matrix at a lower resolution (i.e. a single estimate for a greater number of cells). They will not be presented in this paper, but are a possible avenue future work on this topic.
In estimating the covariance matrix for each baseline and each calibration cell, we are severely limited by the small number of samples in each cell. One way to overcome this problem is to find estimates for the variance of antenna gains, and use these to return to the baseline variances. In this formalism, we extend C to include all visibilities pointing at a single antenna at a given time. Let us begin by writing an expression for the gain vector, which contains the gains for all antennas and all calibration cells:ĝ and the variance on each antenna in each calibration cell is then: As we can see, Eq. 50 is simply a vector form of Eq. 12. The residual gains of Eq. 12 can now be understood as random samples of the covariance between the gains for antennas p and q at a given time, assuming complete skymodel subtraction. We can thus define the variance sample matrix as an estimate of the variance matrix : We define the residual matrix as: where we explicitly place ourselves in the limits of our formalism, i.e. that we do not have direction-dependent gains. We now see that at the core of Eq. 53 liesVĝ τ , where ∑ τ ∈CVĝ τ = V C . The K-matrix is defined as follows: Since the residual matrix depends on the gains, we define the residual visibility vectors as: r C corresponds to a matrix containing all the residual visibilities within one calibration cell C, i.e. for τ ∈ C where g C = const. Let us define n C as the number of elements in each calibration cell. The residual variance sample matrix can now be built by multiplying the residual visibility matrix with its Hermitian conjugate: Note that we do this because it allows us to turn a single noise realisation into a statistical quantity σ. We can relate ⩔ to the variance of individual antenna gains: To reach this point, in Eq. 23, we made the hypothesis that the sky brightness distribution is dominated by spatially incoherent emission. Applying this hypothesis again here, we can make the approximation that the cross-terms in the sum over d, d ′ average to zero: ∑ d,d ′ ≠d ≈ 0. We then have: where ○ denotes the Hadamard or entrywise product and k = Diag {K}. Thus, ⩔ allows us to estimate the variance of each antenna and for each calibration cell by using all the visibilities pointing to that antenna within that calibration cell. With this information, we can then rebuild the baseline-dependent matrix, having improved our sampling by a factor of n ant .

Applying the Correction to Simulated Data
In this section, we show the impact of our weighting schemes on a noise-map made from arbitrarily stronglycorrelated residuals. Here, we assume that our sky contains only a single point source at phase centre: there is thus only a single instance of the noise-PSF, placed at phase centre, to modulate the average variance level. We sample this instance by taking a cross-section from (l, m) = 0 to a large l, keeping m constant. The only difference between these cross-sections is the weighting scheme applied: unit weights for all visibilities ("Uncorrected", blue), sensitivity-optimal weights (green), and artefactoptimal weights (red). We plot both the result predicted by the Cov-Cov relationship (solid line) and the variance estimated across 2000 realisations (dashed line): the result is shown in Fig. 7. The two remain in such agreement throughout the cross-section as to be nearly indistinguishable.
There are a few significant points to note on this figure. Firstly, most of the power in the matrix lies along the diagonal: both weighting schemes thus give good improvements in variance across the map. The artefact-optimal weights, while decreasing the peak further, as expected, also increases the noise far from sources: this is due to the fact that the artefact-optimal weights are in a sense more "selective" than the sensitivity-optimal weights: they upweigh and down-weigh more severely, and will only result in a constant covariance matrix if this matrix is zero everywhere outside of the diagonal. In effect, the noise-map becomes flatter, but much broader.

Applying the Correction to Real Data
In this section, we show the effect of adaptive qualitybased weighting on real data. The dataset used in this section is a single sub-band from an 8-hour LOFAR observation centred on the Extended Groth Strip (α=14:19:17.84,δ=52:49:26.49). The observation was performed on August 28th, 2014. The subband includes 8 channels of width 24.414 kHz each, for a total bandwidth ranging from 150.2 to 150.5 MHz. The data have been averaged in time to 1 data point per second. The data was calibrated using Wirtinger calibration (see Tasse 2014; Smirnov & Tasse 2015, and references therein) and a sky model consisting only of a nearby calibrator source, 3C295. A reference image (a cutout of which is shown in Fig. 6a) was made by calibrating the data according to best practice for LOFAR survey data: 1 calibration solution per 8 seconds and per 4 channels. The residual data was then corrected by the gain solutions and imaged using Briggs weighting (robust=0), pixel size of 1.5 ′′ , and deconvolved using the default devoncolution algorithm in DDFacet (Tasse et al., in press).
The data was then time-averaged to create a new, 2.4 GB dataset with 1 data point per 8 seconds. Deliberately poor calibration was then performed on this dataset, solving for 1 calibration solution per 2 minutes (caeteris paribus). The resulting corrected residual data was imaged using the same imaging parameters as the reference image, and a cutout of the result is shown in Fig. 6b. As expected, the very long calibration intervals introduce calibration artefacts in the image. The brightest sources are still visible, but much of the fainter emission is buried under these artefacts. We are then in a case where our residual visibilities are dominated by calibration error rather than sky model incompleteness.
Weights were then calculated based on the badlycalibrated residual visibilities. Fig. 6c was made using the same visibilities as Fig. 6b and applying baselinebased, sensitivity-optimal weight. Similarly, Fig. 6d used the poorly-calibrated residual visibilities with the application of baseline-based, artefact-optimal weighting. These weights are likely to be poorly-conditioned. In both cases, all other parameters were conserved.
Note that applying antenna-based sensitivity-optimal weighting to the badly-calibrated data (not shown here) allows us to recover the reference image with only a very small increase in rms (increased by a factor of 1.14). Further testing on complex field simulations will be required to ascertain the usefulness of artefact-optimal weighting: it is likely that it fails to correct the image fully due to the poor conditioning of the covariance matrix used here.
The pixel histograms show us that the weights do not completely mitigate the poor calibration interval choice, but certainly give a dramatic improvement over the unweighted, poorly-calibrated residuals. This is compatible with our statement that the weights give similar residuals in the image with a dramatic improvement in time at some cost in sensitivity. It is interesting to note that while Fig.  6d looks noisier than Fig 6c, its residual flux histogram is actually closer to that of Fig. 6a.
As for performance, the weights used for Fig. 6d took 8 hours of computing time on a single core 4 , working on a 29 GB dataset, which is not particularly large for LOFAR data. Since the problem is massively parallel, this cost can be alleviated. The main bottleneck is likely due to very poor code optimization. As for the weights used for Fig.  8b, they are computed in 1min6s on the same single core.
(a) Well-calibrated, unweighted restored image of the sky near the centre of the Extended Groth Strip. Used for comparison with the other images. Units of colourbar are Janskys. This image was made with data calibrated following best practice (solution intervals of 8 seconds, half the bandwidth). rms=5.87mJy/beam (b) Poorly-calibrated, unweighted restored image of the sky near the centre of the Extended Groth Strip. Units of colourbar are Janskys. This image was made with the same data as for Fig. 6a, but averaged in time and calibrated using larger gain solution intervals: 2 minutes and half the bandwidth. rms=86.4mJy/beam (c) Image made using the same imaging parameters and corrected visibilities as Fig. 6b, using sensitivityoptimal weighting. Units of colourbar are Janskys. rms=9.69mJy/beam (d) Image made using the same imaging parameters and corrected visibilities as Fig. 6b, using artefactoptimal weighting. Units of colourbar are Janskys. rms=15.8mJy/beam Fig. 6: Restored images of the centre of the Extended Groth Strip, as seen with an 8-hour observation using the full LOFAR array. Fig. 6a shows an image of the field made with good calibration intervals. Fig. 6b shows an image of the field made with poor calibration intervals. Fig. 6c shows image made with the same visibilities and imaging parameters, but with the application of the sensitivity-optimal weighting scheme. Fig. 6d, similarly, differs from Fig  6c only in that artefact-optimal weights, rather than sensitivity-optimal weights, were used. The histograms of pixel values in each image have 1000 flux bins ranging from -0.16 Jy to 0.16 Jy. Their ordinates are in log scale. Pixel size is 1.5" in all images.

Discussion
This paper began by investigating the use of an algorithm inspired by "lucky imaging" to improve images made using radio interferometric data. By investigating the statistics of residual visibilities, we have made the following findings: -A relationship between the statistics of residual visibilities and residual pixel values (the "Cov-Cov relationship"). -A description of the noise-map in the image plane as a constant variance level modulated by a noise-PSF convolved with the sources in the field. This gives the variance in the flux of the image as a function of distance from the sources in the sky for a given calibration. Fig. 7: The sensitivity-optimal (green) and artefactoptimal (red) weights both give improvements over the unweighted noise-map (blue).
-An Adaptive Quality-Based weighting scheme, which reduces the noise in the image (and the presence of calibration artefacts) by minimising either the constant noise term or the noise-PSF.
While our results are not a panacea for poor calibration, they show that we can not only improve images made with well-calibrated data, but also mitigate the worst effects of poorly-calibrated visibilities in otherwise well-calibrated datasets. Provided that the gain variability timescale is long enough at certain points of the observation, we can effectively get images of similar quality using both the "standard" best-practice calibration interval for LOFAR survey data (calibration solution interval of 8 seconds) and a significantly larger solution interval of 2 minutes (frequency interval unchanged). Of course, if no such stable interval exists, there will be no good intervals to upweigh, and we will be left only with equally-poor data chunks. This means that, in the right conditions, net pipeline time can be sped up by a factor of nearly three, at a slight cost in sensitivity. This increase will be greater than what could be achieved with existing comparable methods such as "clipping".
We emphasize that the adaptive quality-based weighting schemes work because the noise-map describes the underlying noise distribution, of which calibration artefacts are one single realisation. To fully characterise the artefacts, the correlation between different pixels (i.e. offdiagonal elements of Cov{ỹ}) must be computed; this has not been done in this paper. Nevertheless, lesser constraints on the spatial distribution of artefacts can be found using only the diagonal elements of Cov{ỹ}. The weighting schemes merely seek to minimise this spatial distribution as much as possible: the end result is fewer artefacts, which can be distributed across a much larger area. This is the source of the dramatic improvement from 6b to Fig. 6c. We have simply down-weighted those visibilities where spurious signal was introduced by the calibration solutions, and up-weighted those visibilities where such signal was lesser. Since this spurious signal is the source of calibration artefacts, downweighting the associated visibilities reduces it dramatically. The poor improvement from Fig. 6b to 6d is likely due to limits in the conditioning of our estimation of the covariance matrix.
The work presented here can be improved upon, notably by working on improving the conditioning of our covariance matrix estimate: for real observations, it is impossible to have more than one realization of each gain value for all antennas. By treating each visibility within a calibration interval as a realization of the true distribution, we can better estimate the covariance matrix per baseline, and thus reach a better estimate of the variance in the image-plane. Of course, in practice, we can never access to the true, underlying time-covariance matrix for each baseline. Significant hurdles remain: -The impact of sky model incompleteness (since calibration requires a sky model) is ignored in this paper; we start by assuming that we have a complete sky model. In practice, of course, acquiring a complete sky model is often a key science goal in and of itself. The impact of this hypothesis therefore ought to be investigated in future work. -The conditioning of our covariance matrix estimation remains a concern. By using an antenna-based approach, we can improve conditioning by a factor of n ant , but this is only one approach among many. Further work is needed to investigate which method, if any, proves optimal.
(a) Well-calibrated, unweighted restored image of the sky near the centre of the Extended Groth Strip. Used for comparison with the other images. Calibration solution intervals used were 8 seconds, half the bandwidth. rms=5.87mJy/beam (b) Image made using the same imaging parameters and corrected visibilities as Fig. 6b, with the application of antenna-based, sensitivity-optimal weighting. Solution interval of 2 minutes, half the bandwidth. rms=6.69mJy/beam Fig. 8: Comparison between the well-calibrated image (i.e. the same image as Fig 6a) and antenna-based sensitivityoptimal weights. Units of both colourbars are Jansky. We see that we recover a very similar image, despite the fact that the data used for the weighted image are averaged by a factor of 8 compared to those used for the unweighted image.