© ESO, 2014

1. Introduction

The detection and characterization of low frequency radio transients is a key science project (Fender et al. 2006) for the LOw Frequency ARray (LOFAR; van Haarlem et al. 2013) radio telescope. In this context, the Amsterdam-ASTRON Radio Transients Facility and Analysis Center (AARTFAAC)¹ has initiated building an all-sky monitor (ASM) to survey most of the locally visible sky for transients and variable sources. The ASM will be an image plane transient detector at low radio frequencies with images generated via aperture synthesis. It will carry out continuous, autonomous and near real-time sky monitoring for detecting the brightest (~tens of Jy) transients occurring at a variety of cadences (seconds to minutes). Its main aim is the rapid detection of low frequency (LF) transients and the production of their rough position estimates for immediate follow up at high sensitivity and resolution with the full LOFAR telescope. With its very wide field of view, low guaranteed latency of calibration and imaging, it will enable true, real-time all-sky monitoring of the low-frequency radio sky. The AARTFAAC would be the most sensitive of the coming breed of LF Radio Sky Monitors. Table 1 compares various wide-field instruments with transient searches as an explicit goal.

The discoveries of transients are currently biased toward objects emitting at higher energies in the X-ray or γ-ray regime, primarily due to the success of various satellite-based ASMs operating at those energy levels. At radio frequencies, searches for transients have usually been limited to mining archives (Bower et al. 2007; Bower & Saul 2011), follow-up studies of discovered transients (Chandra & Frail 2012), and a few dedicated observing programs (Bannister et al. 2012; Katz et al. 2003). However, the dynamic nature of the sky at low radio frequencies has not been studied in detail yet, mostly due to the restricted fields of view, and the low availability of existing instruments for carrying out dedicated searches for transients. This is set to change with the advent of a number of sensitive, wide field of view instruments at the lowest frequencies – most notably the LOFAR, the Murchison Widefield Array (MWA; Lonsdale et al. 2009) and the Long Wavelength Array (LWA; Ellingson et al. 2009).

These instruments should be sensitive to a variety of phenomena, including Pulsar giant pulses, bright radio flares from brown dwarfs and active stars, emission from OH masers, Jovian bursts, and, most importantly, an unknown population of bright bursters. Recent serendipitous discoveries of short bursts of bright radio emission with unknown origin, and at relatively low radio frequencies (Lorimer et al. 2007; Thornton et al. 2013) have shown the potential variety of sources that might fill the unnaturally empty discovery phase space of LF radio transients. Thus, the development of wide field monitors, like AARTFAAC, is timely and crucial for detecting and characterizing such sources. The science case for the AARTFAAC is detailed in Wijers (in prep.).

Transient detection approaches: a variety of approaches exist for the detection of transient phenomena, depending on factors like the luminosity distribution, spatial distribution, and location of sources (primarily affecting dispersion smearing and scatter broadening), their timescales of variation, emission mechanism, and the spatial and temporal behavior of background noise. These include beam-formed observations, deep imaging (staring), rapid shallow imaging (tiling) etc.

Traditionally, image domain transient detection has been considered suitable only for incoherent sources of emission, whose timescales of transience are slow and whose brightness temperatures are expected to be <10¹² K. Coherent sources of emission have usually been detected using timeseries analysis of beamformed data, although imaging can provide a higher spatial resolution and better discrimination due to the coherent collection of power into a single pixel (compared to incoherent beamforming) and a wider field of view (compared to coherent beamforming). This dichotomy is mostly due to the timescales of coherent emission, which are expected to be short from light travel time arguments and the observation that beamforming can provide high temporal resolutions and sensitivities, although over a reduced field of view.

Another promising alternative for detecting short term transients is the bispectrum method of Law & Bower (2012). This uses interferometric closure quantities from triplets of time differenced visibilities, which makes it independent of instrumental phase errors. However, for an array with a large number of elements and low signal-to-noise ratio (S/N) per baseline, the methods’ absolute sensitivity is lower than that of an imaging array. It is also sensitive to residual flux after the subtraction of constant emission from visibilities. This can arise due to changes in the visibility fringe over the subtraction period and imposes an upper limit on the duration of candidate transients.

The usually sparse UV coverage in snapshot observations from imaging arrays result in low instantaneous sensitivities and image dynamic ranges, mandating long integration. This makes them unsuitable for short-duration transient detection. Another contributor to the lack of imaging detectors at short timescales has been the time-consuming step of calibration and imaging. Recent arrays are being constructed using a large number of small, low sensitivity elements. These provide a perfect base for ASMs due to their very good instantaneous UV coverage, high collective sensitivity and wide fields of view. The challenge of real-time calibration and imaging from these instruments is becoming feasible due to developments in high-speed calibration algorithms and advances in computing resources. Thus, image plane detection of short duration bright transients is now feasible. Further, for transients with an isotropic spatial distribution, unknown duration, and luminosity distributions, a tiling strategy of shallow but rapid observations would result in a larger number of discovered transients (Nemiroff 2003). Imaging ASMs are hence matched for this application.

Fig. 1 AARTFAAC calibration and imaging pipeline within the full ASM flow and leading onto the transient detection pipeline (TraP).

Brief overview of main results: a multisource, model sky-based self-calibration scheme has been found adequate for calibration of the full field of view of the AARTFAAC array. Due to the received source power being dominated by a few, very bright sources and the coarse array resolution, a low-complexity point source sky model is adequate for calibration. Observational constraints, especially due to the ionosphere, require the calibration to be carried out in close to real-time for achieving the necessary high imaging dynamic range. We have developed a calibration scheme capable of handling rapid, direction dependent variations of system parameters, which results in achieving close to thermal-noise limited imaging. This is demonstrated using test data taken under a variety of observational conditions. The scheme has been incorporated into a calibration and imaging pipeline, and an optimized implementation of the algorithm allows the instrument to carry out near real-time calibration and imaging with high dynamic ranges of ~2000:1. To verify that the proposed calibration scheme allows thermal noise limited transient detection, we have analyzed the difference between consecutive snapshot images for which the systematic confusion noise is expected to cancel out. This analysis shows the expected reduction in noise being consistent with a thermal noise limited observation. The achieved noise limit of difference images is found to be ~2–3 Jy over tens of seconds.

This paper is organized as follows. We begin by describing the AARTFAAC array and its suitability as a transient search machine in Sect. 2. We lay out the problem of wide-field calibration in the context of transient detection in Sect. 3 while describing the approaches taken by similar instruments. We then give a detailed description of our approach to calibration in Sect. 4, especially in the context of autonomous operation. Section 5 presents the observed performance of the calibration approach on commissioning data, while Sect. 6 elaborates on the computational performance of our strategy. Section 7 elaborates on some of the challenges faced by the AARTFAAC ASM under real observing conditions after which we present our conclusions. The results presented in this paper have been obtained on test observations carried out by using existing LOFAR system infrastructure, while the AARTFAAC Uniboard² based piggyback recording system was being built.

2. The AARTFAAC all-sky monitor system overview

The AARTFAAC all-sky monitor is an aperture array of 288 sky-noise limited dual polarized antennas (the low band antenna or LBA), which are shared with the LOFAR telescope. These operate between 30 and 80 MHz. They are organized as six stations spread over ~300 m, with each being a random array of 48 antennas, which can be spatially organized in one of a limited number of configurations. The array shares its analog elements with the LOFAR telescope, which thus controls the station subarray configuration. The received analog signal after basic analog conditioning is baseband sampled at 200 MHz with a 12-bit quantization. All dipoles are sampled with a single clock, eliminating differential delay errors between antennas due to such errors as clock misalignment or drift. A hardware digital polyphase filter bank splits the ~100 MHz sampled band into ~192 kHz subbands. An 8-bit complex representation of the subband timeseries has been found to be adequate to maintain linearity in the presence of local RFI. A coupled data path routes a configurable subset of subbands to the Uniboard based hardware FX correlator. This makes the AARTFAAC a continuously available instrument, independent of ongoing LOFAR observations. No delay compensation (apart from the fixed cable delays) is carried out between antennas, resulting in a zenith pointing point spread function (PSF). The estimation of the resulting ~1.6e5 instantaneous visibilities makes this correlator the largest in the world in terms of processed input streams. The I/O restrictions put an input limit of ~66 subbands (~13 MHz) to the correlator, which produces visibilities at 16 kHz and one second resolution. Note that the LOFAR also has a high band antenna (HBA) component operating between 110 and 240 MHz with each element being a 4 × 4 tile equipped with an analog beamformer. This implies that data is available to the AARTFAAC during HBA observations, which are mutually exclusive to LBA observations. However, due to the limited field of view of ~20 sq. deg, and the arbitrary pointing of the analog beam of this system, its use is not discussed further in this paper.

Table 1 lists the current specifications of the instrument and compares it with other low frequency wide field instruments with transient detection as an explicit goal. Figure 1 shows the functional components of the AARTFAAC monitor with details of the near real-time transient detection pipeline. The calibration component of the pipeline is the subject of this paper.

2.1. AARTFAAC for transient searches

Temporally pulsed emission from a transient source interacts with an ionized intervening medium and is distorted due to dispersion and scattering. These effects can be parameterized by the medium’s dispersion measure, the integrated electron density along the line of sight (DM, pc/cm^-3). The DM is also used as a proxy for distance. Dispersion causes an arrival time delay of the pulse at different frequencies and leads to pulse broadening if the signal is not dedispersed before spectral integration. However, this broadening can be completely corrected if the DM of the source is known. Scattering is caused by electron density fluctuations along the line of sight and leads to pulse broadening due to arrival of delayed pulses via multiple paths. This effect, however, cannot be corrected. Both effects are enhanced at large DMs and low frequencies and can limit the detection of pulses, instead of the limitation coming from the inherent luminosity and distance of the sources.

Since pulse attenuation is proportional to broadening, transient detectors routinely increase the pulse detection sensitivity by carrying out trial DM searches. However, the AARTFAAC does not carry out de-dispersion before imaging. Thus, given a sensitivity threshold and a temporal resolution of ~0.1 to 1 s of the AARTFAAC, the transient sensitivity achieved via dispersed spectral collapse limits the DM range of a detectable source, while scatter broadening limits the DM range irrespective of spectral collapse.

An empirical relationship between the scatter broadening time due to the interstellar medium and the DM has been derived by Bhat et al. (2004). This implies a DM upper limit of ~100 at 60 MHz when lines of sight are within the Galaxy and when scatter broadening is restricted to 1 s.

The temporal pulse broadening in microseconds over a bandwidth Δν MHz at ν GHz for a dispersion of DM units is given by Eq. (1). (1)Imaging at the highest available spectral resolution of 16 kHz and restricting dispersion broadening to within 1 s at the same time allows detection of a source with a DM of up to ~1600. This allows probing along any galactic line of sight (Cordes & Lazio 2002), although at a sensitivity that is a factor ~28 times poorer than that with the full 13 MHz band.

The AARTFAAC imager does not carry out DM measure searches via de-dispersion on trial DMs. Thus, attempts to increase the transient detection sensitivity of AARTFAAC by spectral integration over its full 13 MHz band would lead to a lowered detected S/N for any pulses with DMs greater than a few. This is because dispersion broadening would spread the pulse beyond the 1 s imaging cadence, restricting the transient search spatial radius. However, as stated above, scattering along Galactic lines of sight broadens pulses with DM ~100 to beyond 1 s. Thus, some level of spectral integration can be afforded by this limit before scattering, and dispersion broadens the pulse beyond the AARTFAAC cadence. Finally, as shown by Hassall et al. (2013), de-dispersion is not beneficial for highly scattered sources. The scatter broadening of the pulse increases more rapidly with DM than the dispersive delay. At high DMs, this makes the inherent pulse width (due to scatter broadening) larger than the dispersion broadening. This would allow the AARTFAAC to probe the high DM domain with full spectral sensitivity.

The ability to spread subbands over the full analog bandwidth allows us to search for ultra steep spectrum sources and can aid in distinguishing between terrestrial and celestial transient sources. The large spectral baseline also help in carrying out image-based dedispersion in the future for enhancing the transient detection sensitivity of AARTFAAC.

A figure of merit for transient detection (Cordes et al. 2004) from a telescope is and can be used for comparison between instruments. Here, A is the total collecting area; Ω is the solid angle coverage; T is the total duration of observation; and Δt is the time resolution. Due to the limited duration of fast transients, it is relevant to add an instantaneous sensitivity parameter to the figure of merit, as represented by the system temperature T_sys. Large bandwidths increase transient sensitivity only if de-dispersion with the correct DM value is carried out. Since this requires a computationally expensive DM space search, we consider a narrow band figure of merit. Assuming continuous availability and the simultaneous imaging of the entire field of view, a per second figure of merit is AΩ /T_sys. Among currently planned ASMs, the AARTFAAC’s FoM is comparable to the LWA, which has similar field of view (see Table 1). However, the AARTFAAC’s higher spectral and temporal resolution, as well as its larger processed bandwidth would give it a higher sensitivity toward scatter broadened transients, as elaborated in Sect. 2.1.

Backer (1999) has shown that the S/N achieved in a continuous wide-field search scales with the array filling factor f as . Thus, the highly sampled aperture of the AARTFAAC dense array would be an advantage for our application. The resulting zenith pointing PSF is stable and repeatable in time. At its peak sensitivity at ~58 MHz, the PSF has high roll-off, low side lobes ~15 dB below peak with a resolution of ~0.8 sq. deg. Assuming the confusion limit to be reached with the presence of one source per ten PSFs and the source counts from Bregman (2012), results in a classical confusion noise of ~10 Jy. The large number of antennas also make the AARTFAAC a very sensitive instrument with a combined A_eff/T_sys of ~0.7 m²/K (Wijnholds & van Cappellen 2011) and an instantaneous thermal noise of ~8 Jy for the standard imaging mode of 1 s/16 kHz. This implies that the sensitivity of typical snapshot images are limited by confusion noise and not thermal noise.

The array is 2D to high accuracy (~1.4 cm rms in the z-axis); thus, imaging the zenith region can be done via a simple 2D Fourier transformation of the calibrated visibilities. Lines of sight away from the phase center (the zenith) have a neglected phase term due to the w-component of the observed visibilities. However, the co-planarity of the array results in a spatial coherence sampling function that is confined to a plane passing through the origin, keeping the 2D relationship between the visibility and imaging domains (Cornwell & Perley 1992). Thus, wide-field imaging is simplified as no 3D transforms or their approximations like W-projection (Cornwell et al. 2008) are required. Simple faceting recovers off-axis lines of sight with high fidelity, making latency bound imaging feasible. The requirement of snapshot imaging implies that the sky observed by every baseline remains constant within a visibility set, thus making the calibration parameters amenable to estimation by self calibration. Since all dipoles are physically aligned parallel to each other, the instrumental polarization can be determined via a bi-scalar polarization calibration, and corrections can be applied in the image domain.

Dedicated hardware is used to achieve low latencies of transient detection. These include a hardware correlator based on the Uniboards for estimating the visibilities in real-time, and a software solution running on general purpose hardware for the calibration, imaging, and transient detection pipelines. Thus, the AARTFAAC array is well suited for its role as an ASM.

3. All-sky calibration

Calibrating a radio array telescope entails estimating the instrumental and (possibly environmental) contribution to observed visibilities. These can be parameterized partially by per antenna complex gains and noise. Observations of sources with accurate models allow the estimation of these parameters and can enhance the sensitivity of the instrument substantially (Taylor et al. 1999). An advantage of all-sky instruments is the availability of enough bright sources within the beam for in-beam calibration, making dedicated observations of calibrator sources unnecessary.

Calibration of the full field of view (i.e. beyond the half power beam width) is usually necessary for accurate estimation of the flux contributed by bright sources to any pointing of the synthesized beam via its far sidelobes, aiding in their accurate deconvolution. A wide-field image plane transient detector typically searches image timeseries for transients by extracting and analyzing lightcurves of existing sources, or by detecting new sources referenced to a database of lightcurves. The Transient Detection pipeline (TraP; Swinbank, in prep.) component of the AARTFAAC ASM carries out source extraction from individual time slices via thresholding and fitting 2D Gaussian models to regions of emission. Hence, high point-source sensitivity and accurate flux recovery across the full field of view is a fundamental requirement of calibration for transient searches, while a slight error in source positions can be tolerated by the pipeline via spatial source association. The flux estimation of any source via imaging directly depends on the flux calibration of the instantaneous visibilities, which in turn depends on a proper division of the uncalibrated flux on a visibility between element gain, system noise, and source flux. Due to the dominance of emission from our own Galaxy on system temperatures, the system noise can vary by almost 40%, depending on the pointing of the array. Thus, system gain and temperature need to be continuously estimated via calibration.

An important requirement for real-time transient detection is a hard deadline on calibration time and a limited computational budget in view of the large number of parameters to be estimated. These dictate quick and timebound convergence of the calibration parameter estimators. Further, the streaming nature of the application and the large number of visibilities makes their buffering or multiple passes of processing unfeasible.

Finally, the calibration should be robust against sources of terrestrial transients due to, for example, radio frequency interference (RFI), which can cause the algorithm to converge to a sub-optimal solution. However, the most dominant of these effects is due to the ionosphere, which is discussed next.

Fig. 2 Block diagram depicting the snapshot iterative calibration control flow in the AARTFAAC optimal calibration scheme.

3.1. Ionospheric effects

The antenna complex gain has contributions primarily from the gains of the electronics, receiver characteristics, and telescope geometry. These usually vary slowly (~hours), are direction independent effects (DIEs), and are estimated via observation of calibration sources or using self-calibration. At low frequencies, the disturbed ionosphere can act as a phase screen with temporally and spatially varying refractive and diffractive effects. Ionization of the ionosphere during the day by the solar wind and photoionization is balanced by recombination during the night. These can cause bulk variations in large scale Total Electron Content (TEC) with slow temporal variations, as well as fluctuations on relatively small temporal and spatial scales (~0.1% variations over a few kilometers and over tens of s). If the rapid TEC fluctuation sets up an instantaneous spatial phase gradient across the array in the direction of a source, an apparent position shift of the source is seen. However, if the phase structure function varies from a gradient, a deformation of the source structure may result due to defocusing, resulting in a reduction of the recovered source peak flux.

Since the arrival time delay per antenna depends on the TEC along the line of sight to the source, these are termed direction dependent effects (DDEs). A phase correction varying over the field of view then needs to be applied, and a single phase correction per antenna (as obtained by traditional self calibration) is insufficient. Uncalibrated phase errors cause scattering of source power into sidelobes, affecting the reliability of the extracted lightcurves. Further, a change in the apparent source shape leads to an increase in residual sidelobes after deconvolution, as the mean source model deviates from the apparent instantaneous sky, and the source subtraction is incomplete. These residual sidelobes increase the background noise level and can introduce spurious structure into the image.

During DDE parameter estimation, we make the simplification that the effects are identical for all antennas. This is a reasonable assumption for a compact array, whose size projected to ionospheric heights are expected to be smaller than the typical isoplanatic patch. Thus, identical lines of sight from individual antennas are expected to traverse an ionospheric patch with the same refractive properties. The estimation of a complex gain in the direction of each source is required, but this can be assumed to be identical for all antennas. Thus, a multisource, model-based calibration approach is appropriate for the calibration of wide-field, compact arrays (Wijnholds et al. 2010).

Although instrumental parameters are expected to vary on longer timescales (hours to days for LOFAR), their precise calibration also mandates temporal oversampling. If uncorrected, these effects can raise the image noise floor and contribute to variations in the recovered flux of sources, leading to false positives. Thus, real-time calibration of both DIEs and DDEs is necessary for maximizing transient detection sensitivity of autonomously operating, near real-time instruments.

3.2. Current approaches to all-sky calibration

Typical calibration schemes, such as those implemented in the CASA package, address a different problem than what is applicable to the AARTFAAC due to the restricted field of view, much longer baselines, and an order of magnitude fewer stations in typical arrays. Their algorithms are further complicated by the need to compensate for higher order effects like the temporal variation of the primary beam during long synthesis due to changing array geometry, or the beam rotation of altaz mounts, etc., which are moot in our application. Snapshot, zenith pointing imaging from a co-planar array simplifies the calibration and imaging from the ASM to a certain extent.

A successful approach to wide-field calibration is the Peeling algorithm (Noordam 2004; van der Tol et al. 2007), in which self-calibration toward sources within the field of view is carried out in decreasing order of their brightness with removal of the brightest remaining source’s calibrated contribution to all visibilities. This source is selected via rephasing the array to its direction, and the contribution of other sources is attenuated by averaging over their unphased visibility contribution (fringe washing). Thus, a single source calibration approach is adequate, and a least squares estimation of the antenna gains in the direction of the source can be carried out using one of several methods (Boonstra & van der Veen 2003). The solution yields a set of time-variable, antenna-based phase corrections per source, and a source model. However, the algorithm requires the storage of the estimated visibilities with multiple calibration and imaging passes through the data before the iterative peeling can be concluded and calibrated images produced. This, along with the increased computational load, makes it unsuitable for real-time operation. Other contemporary instruments with similar specifications and goals like ours mostly use Peeling based approaches.

The MWA (Lonsdale et al. 2009), located in Western Australia, has developed an algorithm for carrying out real-time, full polarization direction- dependent calibration of snapshots of ~8 s via a Peeling mechanism (Mitchell et al. 2008). They use the sequential deconvolution of several bright sources in decreasing order of apparent flux to iteratively fit apparent angular offsets induced by ionospheric phase shifts and antenna gains toward these sources. This is suitable for the much reduced field of view of the MWA as compared to the AARTFAAC, allowing them to fit a “rubber sheet” model of the ionosphere over their field of view. The main departure of our approach to theirs is in the estimation of image wander and the deconvolution of sources affected by the ionosphere. The sequential nature of peeling ensures that the sidelobe contamination of the brighter source does not affect the flux measurement of the weaker source and hence, its deconvolution. In our case, the apparent fluxes of the brightest sources are simultaneously estimated, and their effects subtracted together. This is efficient computationally and is aided by the filled nature of the AARTFAAC UV plane, which results in low sidelobe contamination. Further, the ionospheric effect is estimated by fitting a phase ramp across the available bandwidth in the MWA, which requires synchronization between the spectrally separated calibration threads. Our estimation uses a direction of arrival (DoA) algorithm independently for each spectral channel and does not require large bandwidths for fitting phase ramps. This allows a loosely coupled parallel architecture for our implementation. The AARTFAAC calibration has a stricter latency bound, and its lowered sensitivity and better PSF reduce the benefits of an iterative deconvolution.

The LWA (Ellingson et al. 2013) is colocated with the VLA in New Mexico (USA). Along with forming beams in realtime (its primary mode of operation), it can also operate in an “all-sky” mode for a single tuning of 70 kHz. The calibration of antenna phases utilizes the method of fringe-rate filtering using narrow band data to localize one of multiple point source calibrators within the field of view. Here, the time-varying phases of visibilities due to contribution by discrete sources is converted into a fringe rate spectrum by taking a Fourier transform of the correlation, where the sources are apparent as localized components. The fringe of the brightest source is filtered, and antenna delays extracted by fitting a model of the presumed delay response to the measured phase response. Such an approach provides a phase calibration solution in the unique direction of the brightest source. However, fitting a fringe requires an extended period of time to generate the fit baseline, which is not suitable for our streaming, real-time application. Further, given the calibration approach, the presence of a disturbed ionosphere would add a random noise on the parameter estimates, while any bias due to ionospheric systematics would probably be averaged out over the time period over which the fringe is fitted. This could lead to a bias in the snapshot calibrated visibilities.

Field-based calibration (Cotton et al. 2004) is another technique for direction-dependent ionospheric phase calibration and was developed for the VLA 74 MHz observations. Over a time interval of 1–2 min, it measures and converts the apparent position shifts of 5–10 detectable bright sources within the field of view into ionospheric phase gradients over the array. Subsequently, an independent phase screen (based on a 5-term basis of Zernike polynomials) is fitted onto the observations and is used to predict phase gradients in arbitrary viewing directions to image the full field of view. Such a technique is applicable to a limited field of view of a high-sensitivity system, which implies the presence of bright sources that are close enough to interpolate the phases attributed to the ionosphere. The AARTFAAC wide field of view and lower sensitivity preclude the interpolation between the widely spaced phase calibrator sources. This implies that sources in locations intermediate to model sources may be affected by ionospheric amplitude and position variations, which result in the presence of structured variations in a lightcurve for such sources.

4. An optimal tracking calibration scheme for transient detection

Here, we present our approach for latency-bound and precise calibration of the AARTFAAC with a low computational footprint. We use a model-sky based, multisource self calibration approach (Wijnholds & van der Veen 2009a). Calibration of an array of identical, wide-field antennas entails deriving the direction dependent response of every antenna, of which each is excited by multiple, simultaneously present calibration sources. The system noise contribution is also estimated and subtracted while calibrating the visibilities.

Algorithmic Summary: as seen in Fig. 2, the calibration of the system for each incoming time and frequency slice is carried out iteratively in two parts: a major cycle, which carries out calibration of direction dependent effects, and a minor cycle, which carries out calibration of DIEs.

• 1.

  The input to the algorithm is a running timeseries of measuredspatial correlations between all antennas (visibilities organizedinto an array covariance matrix (ACM)) for a given time,frequency slice, and the corresponding instantaneoussky-model.

• 2.

  The minor cycle within each major cycle iteratively estimates the complex gains of all antennas using the StEFCal (Salvini & Wijnholds 2014) algorithm. This is done by iteratively fitting a parameterized model to the observed visibilities with the model itself being generated by the encompassing major cycle.

• 3.

  The major cycle then estimates the direction dependent effects (modeled via the apparent flux and positions of model sources) and system noise using the latest antenna gains and updates the model with these higher accuracy parameters.

• 4.

  The time axis is used to feed forward parameter estimates as an initial guess for the next time slice, leading to faster convergence.

• 5.

  The frequency axis is compensated for the known spectral response of the dipoles during bandpass calibration, allowing spectral averaging of the calibrated visibilities. Finally, an appropriate flux scale is provided by scaling the recovered flux of a few model sources to their actual, known fluxes.

In the following, we first describe the calibration of an ACM corresponding to a single time slice, termed “Instantaneous calibration”. The latency and computational constraints for real-time calibration suggest utilizing the available time axis for rapid convergence. This is implemented in the tracking component of our algorithm, which implements a temporal feed-forward of filtered solutions as initial estimates. It may be noted that our calibration algorithm is limited to operating in the visibility domain and does not use the image domain to determine source structure, such as via CLEAN components. This is beneficial for a streaming calibration approach.

Since the telescope gains are frequency dependent, the received band is divided into multiple channels, and the parameter estimation is carried out separately per channel. A fundamental assumption is that the bandwidth over which the ACM is formed is narrow enough for phase rotations to compensate for temporal delays between the reception of the signal at different antennas (Zatman 1998). This is appropriate as the maximum time of flight on the longest baseline of 300 m (~1 μs) is much less than the inverse of the typical calibration channel bandwidth of 16 kHz.

4.1. Instantaneous calibration

Array calibration can be formulated as a nonlinear parameter estimation problem, which can be solved using statistically efficient estimators generated via a nonlinear least squares approach. Our approach to calibrating a single incoming timeslice closely follows that of Wijnholds & van der Veen (2009a), which asymptotically establishes efficient estimators for the model parameters via a maximum-likelihood (ML) formulation. Given a set of p antennas whose locations are known precisely, q calibration sources with nominally known positions and fluxes within the field of view, the data model developed for the instantaneous ACM R in their scheme is given by (2)where is the diagonal matrix of DIE gain of each antenna, A is the matrix of array steering vectors which are generated using the known positions of the array antennas, is the diagonal matrix of DDE gain in the direction of each of the q sources, is the known flux of the sources in the sky model, and Σ_n is the diagonal noise covariance matrix assuming an uncorrelated but nonidentical system noise contribution from each antenna.

The model parameters to solve for are of the form (3)where γ_i is the direction independent gain of each antenna and φ_i is the associated phase. The variables σ_si and are the estimated fluxes and positions of the model sources, while σ_ni is the system noise, assumed to be independent for each antenna.

Since all signals are assumed to be independent, identically distributed Gaussian random variables, the normal equations for this estimation problem are generated by minimizing the negative log-likelihood function: (4)where R(θ) is the model covariance matrix as a function of the parameters θ, and is the sample covariance matrix.

Due to the difficulty in solving this minimization problem in closed form, a numerical optimization based on a weighted least-squares, covariance-matching estimation technique (COMET) is utilized (Ottersten et al. 1998). This is known to lead to estimates that are equivalent to ML estimates for a large number of samples, thus being asymptotically efficient and reaching the theoretical limit of the error on the estimator, the Cramer-Rao lower bound (CRLB).

The estimation of the instantaneous calibration parameters is carried out by partitioning the parameters into subsets, which are then estimated using a least squares approach. This approach of alternating between parameter subsets using the best current estimate of a subset in an iteration to estimate the other subsets of parameters is called weighted alternating least squares (WALS).

4.1.1. Model generation

Forming the normal equations of Eq. (4) for minimization requires modeling the observed visibilities using the data model in Eq. (2) and the current estimates of the parameters. Due to the limited resolution of the array, most of the brightest objects (apart from the Galactic plane) are unresolved by the AARTFAAC. Thus, the sky can be modeled by a collection of points at the nominal locations of the bright sources. We use a sky model composed of the brightest sources in the Northern sky, viz. Cas A, Cyg A, Vir A, Tau A, and the Sun, whose nominal positions are determined from the revised third Cambridge (revised) catalog (3CR). The model parameters used are summarized in Table 2. The efficacy of this model for parameter estimation is elaborated upon using real data in Sect. 5.1. Of these model sources, Cas A is visible throughout the year from the location of the AARTFAAC.

Galactic emission, having a steep spectrum, is extremely bright at LOFAR frequencies and also difficult to model due to the detailed structure resolvable at the AARTFAAC’s resolution and sensitivity. Hence, we filter out the low spatial frequencies to suppress this emission during calibration. The model can then account for a large fraction of the received flux on the filtered visibilities, and is constrained by the fluxes and positions of the model sources, as well as the noise model.

The model building requires an estimate of the actual instantaneous flux of the model sources, along with the current estimate of the direction-independent gain solutions and the noise. Typically, an accurate sky model provides the model fluxes as well. However, significant flux variation of the model sources can be introduced by the ionosphere over the AARTFAAC array and by the direction dependent gain. Hence, estimating the model source fluxes via an efficient estimator results in the generation of a consistent sky model. Section 4.1.3 summarizes this approach. Note that a relative flux stability is more relevant for a transient detection instrument and can be achieved to a higher accuracy than an absolute flux extraction.

A COMET based estimator for the noise covariance has been shown to be well suited if a linearly parameterized model of the noise covariance matrix is used for this problem (Ottersten et al. 1998). However, antennas in a crowded array like the AARTFAAC can mutually couple with each other, resulting in coloring of the noise. Our noise model also includes these effects and is discussed in more detail in Sect. 4.1.4.

4.1.2. Direction independent gain estimation (minor cycle)

Once an appropriate model for the observed visibilities has been created, direction-independent gains can be estimated using an appropriate estimator. The AARTFAAC calibration is latency bound; hence, we implement the StEFCal (Salvini & Wijnholds 2014) algorithm for gain estimation. This algorithm has a smaller memory footprint, and its computing scales as for N parameters, compared to the scaling of traditional solvers. This significantly reduces the computing during AARTFAAC calibration due to the large number of parameters (~590) being estimated. The StEFCal implementation led to a measured increase in calibration throughput of a factor ~30 (Salvini & Wijnholds 2014) compared to a traditional estimator and makes real-time calibration of every time and frequency slice feasible. This speed gain is primarily due to linearizing the cost function, allowing the solution to be computed analytically and hence avoiding the need for matrix inversions. Although computationally costly, these are essential to other solvers like the Levenberg-Marquardt minimizer. The StEFCal algorithm has been shown to be unbiased with rapid convergence.

Gain constraints: since we simultaneously solve for apparent source fluxes and antenna gain amplitudes, a constraint on either the source fluxes or the gain amplitudes are needed. However, the variability of model source fluxes due to ionospheric scintillation prevents their use to constrain the self-calibration. Doing so couples the estimated gain solutions to variability in the model fluxes, thus transferring the scintillation to calibrated visibilities and hence, over the entire field of view. To break this coupling, we use instead a constraint on the average gain amplitude over all antennas, setting it to unity. This reflects the behavior of antenna gain amplitudes, which were observed to be stable to a few percent across all antennas and over long times during commissioning. The phases are referenced to the first antenna in the array as an additional constraint. These constraints guide the traversal of the χ² error surface toward an identifiable global minimum, whose convergence is identified by a slow enough rate of change of the estimated parameters.

4.1.3. Direction dependent gain estimation (major cycle)

The major cycle creates an updated model sky using the higher accuracy parameter estimates from the minor cycle. It uses this for estimating direction dependent effects, which are parameterized in the apparent fluxes and positions of the model sources. The cycle begins with a one-time estimation of the model source positions via the WSF algorithm.

Pre-calibration model source position estimation: the tracking calibration approach also incorporates a direction of arrival (DoA) algorithm using a vector subspace technique called weighted subspace fitting (WSF; Viberg et al. 1991). This allows us to estimate the instantaneous position of a model source, thus compensating for the image wander of the sources due to ionospheric refraction. The WSF is an algorithm for determining the directions of multiple narrow-band, far-field emitters, whose uncorrelated signals are received by an arbitrary (but known) array geometry. It is a high-resolution algorithm with a resolution better than λ/D (D being the array size), which is possible for sources with a high S/N. It operates by partitioning the received ACM into a signal and a noise subspace based on its eigen-structure. It then confines the signal subspace corresponding to the true signal parameters to the subspace spanned by the array manifold vectors, which represent the functional dependence of the array response vector to a source over the region of interest in parameter space. Then, the signal parameters are estimated by finding the best least-squares fit of the two subspaces. A crucial assumption is that the array manifold is known, implying that the array is calibrated. To satisfy this requirement, our approach generates initial calibration solutions based on the catalog positions of the model sources, calibrates the array, and then carries out a WSF estimation.

Bias of the WSF estimator: the requirement for the majority of flux in the ACM to be contributed by the model sources implies that the algorithm can produce biased solutions in the presence of unmodeled flux. Tests of WSF show that spatial filtering of diffuse emission is adequate to prevent this biasing.

Thus, WSF along with LS imaging allows us to accommodate the effect of the ionosphere into the model visibilities, leading to lower modeling errors, and hence, better estimates of the calibration parameters. Calibrating the data then requires applying the DIE solutions to the incoming data, while the DDEs are corrected in the image plane via a beam model, which incorporates the effects.

Model source flux estimation: as explained in Sect. 4.1.1, instantaneous estimates of the apparent (deconvolved) model source flux are required for model generation. This can be framed as a parameter estimation problem, for which a closed form estimator based on the COMET approach has been developed in (Wijnholds & van der Veen 2008) in an approach called least squares imaging. Here, the array configuration is assumed to be known to a high accuracy, while the number of sources requiring flux estimation is assumed to be smaller than the number of resolution elements in the field of view. This approach allows the deconvolved flux of a model source to be estimated due to the incorporation of a deconvolution operator as part of the estimator.

The model source flux estimator uses a sparse weighting matrix due to the modeling of mutual coupling via a nondiagonal noise matrix (Sect. 4.1.4), which prevents a bias in estimation.

4.1.4. System noise estimation in the presence of Galactic emission and antenna mutual coupling

Both Galactic emission and coupling between closely spaced antennas can cause correlated noise on the visibility formed from a pair of antennas. Mutual coupling due to the antennas being placed closer than half a wavelength at the lower frequencies can result in fluctuation of the PSF gain and sidelobe levels as a function of pointing direction (Agrawal & Lo 1972). This implies that antenna pattern multiplication (assuming that a single element of an array behaves similarly to it being isolated) may not apply. Mutual coupling is a short baseline effect and can be ignored between the antennas of the different stations making up the AARTFAAC. The impact of mutual coupling to the antenna primary beam in the typical LBA_OUTER mode of observation is expected to be <−17 dB from simulations (Wijnholds & van Cappellen 2011).

The very bright synchrotron background of the Galactic plane at low frequencies also poses problems. Firstly, the total system noise temperature being dominated by Galactic emission also leads to correlated noise between antennas. This, in turn, again leads to a poorer sensitivity of the PSF formed by the inclusion of these baselines as a function of pointing direction and away from the zenith, when compared to the sensitivity expected after pattern multiplication, as found by Ellingson (2011). They also found that the effect of the Galactic correlated noise dominates that due to mutual coupling and is a function of the location of the Galactic plane within the field of view. The Galactic noise correlation is expected to reduce significantly after a separation of a few tens of wavelengths between antennas forming a baseline so it can also be termed a short baseline effect. Secondly, the Galactic plane presents a hard to model source structure due to its extended nature (tens of degrees around the Galactic equator). Even if modeled, calibration would be complicated by a model with a much larger number of parameters than a point source model. Finally, for high sensitivity transient detection, its contribution to the image brightness needs to be subtracted due to the limited dynamic range of the generated images.

The coupling effects and the correlation of Galactic noise have been phenomenologically modeled in Wijnholds & van der Veen (2009b) by an additive nondiagonal noise covariance matrix (defined over a limited set of short baselines) in the measurement equation. A WALS approach is taken to estimate the parameters in the nondiagonal noise covariance matrix. This allows a simple point source model to be used for calibration, leading to a lower bias in calibration parameter estimates. Further, the calibrated visibilities can subtract the contribution due to mutual coupling, correlated sky noise, and the receiver temperatures of each antenna and achieve a PSF with higher sensitivity, although estimates of the flux can be biased (4.1.3). Hence, all our estimators ignore the short baselines (defined as <10λ) to avoid the interaction between the excess flux at the shortest baselines and the calibration parameters.

4.1.5. Bandpass calibration

The ASM operates in a narrow-band mode by default. The incoming raw complex voltages from each antenna are organized as several subbands of ~192 kHz, which are further filtered to give the correlator spectral resolution of ~16 kHz. These subbands can be spread anywhere within the 100 MHz baseband bandwidth via the configurable hardware and, as such, sample different components of the receiver bandpass. The ASM generates channel-collapsed calibrated visibilities for imaging spanning multiple subbands which are typically over a MHz. However, the LBA bandpass can have variations of up to 80% from the peak resonance at ~58 MHz (van Haarlem et al. 2013). Thus, bandpass calibration is necessary, and the calibrated visibilities are scaled by the known bandpass response of the LBA unilaterally. This is adequate to spectrally collapse narrow bandwidths and allow a multifrequency synthesis mode of imaging.

4.1.6. Flux-density calibration

The low resolution of the AARTFAAC allows the use of a number of bright sources for calibrating visibility amplitudes. Recently, Scaife & Heald (2012) have modeled the low frequency spectra of several bright and compact sources to provide a low-frequency, broadband flux scale. Their set of six sources span the full RA range, and a few of these sources are typically visible in every AARTFAAC snapshot due to its very wide field of view. Hence, the AARTFAAC fluxes are tied to this scale. The dipole primary beam response is derived from a full-polarization simulation incorporating all antennas of the AARTFAAC with individual ground planes. The simulation incorporates the effects of mutual coupling between closely spaced antennas. An average response over all antennas is then used as the primary beam response, which is applied as a correction to the full field of view in the image domain.

Table 3 presents the ratio of the integrated fluxes of several bright sources in test observations to the expected value from literature. The integrated fluxes are generated by fitting 2D Gaussians to sources detected in a snapshot image and taking the volume of the fitted Gaussian. Lightcurves are then generated for the various observations. A mean integrated flux in pixel counts is estimated by first segmenting the data into 100 s blocks, finding the median of each block, and then taking the mean of the medians. The calibrated fluxes of several bright sources within the field of view are then compared with those derived from literature, as indicated in the comment column. The dominant source of error to the flux ratios is ionospheric scintillation, which creates systematic variations at timescales of tens of seconds. Another contributor is the error in modeling the primary beam. The fluxes and errors from VLSS have been scaled to 60 MHz assuming a spectral index of −0.7.

4.1.7. Time/spectral bin for calibration

The time and spectral binning of the data has important consequences on the range of science parameter space explored, as well as on the computational load of the calibration pipeline. High spectral resolutions allow better bandpass calibration and RFI discrimination, while higher time resolutions allow probing a larger DM range. Both, however, have a computational cost. Adequate S/N of the model sources on each baseline is required for maintaining reasonable errors on the parameter estimates during calibration. At LOFAR frequencies, model sources sufficiently dominate the visibilities to allow calibration using only a few thousand time samples, so their S/N does not constrain the resolution. Considering the wide field of view, bandwidth smearing effects can affect the flux of sources at the edge of the field. Keeping a limit of 10% flux attenuation at the edge of the field for a fractional bandwidth of 10^-3 puts a limit of ~48 kHz on the channel width. The time resolution is constrained by the feasible latency of the calibration and imaging pipeline and is chosen to be 1 s. Thus, in our standard autonomous mode of operation, we choose a resolution of 16 kHz and 1 s for real-time calibration. The calibrated visibilities or images can be further integrated over a larger time or frequency span in accordance with transient detection requirements.

4.2. Tracking calibration

An alternating least squares approach is already known to have high computational efficiency over closed-form estimators and fast convergence (Boonstra & van der Veen 2003). Thus, this approach is particularly suitable for our case due to our latency and compute budgets. However, the least squares algorithms are quite sensitive to an initial estimate, which can otherwise lead to convergence to local minima. Due to the high temporal sampling of the calibration parameters, a very good initial estimate can be the solution from a previous time slice. Thus, considering the known stability of the instrument (see Sect. 5.3), and the required frequent re-estimation of the calibration parameters, we optimize the computational load of the algorithm by reducing the number of iterations of the iterative solver per time slice to just one. This is supplemented by propagating solutions from previous time slices as initial estimates for the model parameters, which keeps the algorithm on track.

Spectral stability can also be utilized to track solutions in the frequency axis; however, non-contiguous subbands may be chosen during operations to optimize detection of steep spectrum sources by sparsely sampling the full bandwidth. The nature of solutions is then sufficiently different to prevent tracking along the spectral axis. Thus, we choose to use the time axis for solution tracking. Note that computationally, StEFcal iterations are relatively inexpensive as compared to, for example, WSF (see Table 5). This implies that their number can be increased to handle data with larger errors with little effect on latency bounds. Thus, the results from a single iteration show the worst case calibration (and best case computational) performance of our approach.

Fig. 3 Schematic view of tracking calibration for bounded latency precise calibration.

The feedback of solutions as initial estimates makes our approach susceptible to errors in the case of bad data leading to incorrect solutions. We prevent this by maintaining a short time window of solutions against which all new estimates are compared before being used for feedback. In the case of an inconsistent solution, an average gain over the window is used. This approach works well to eliminate biases in estimates due to incorrect tracking. A schematic view of the tracking calibration and its major components as described is shown in Fig. 3.

The StEFCal algorithm already utilizes the gain estimate from a previous iteration. Although it is seen to converge quite rapidly for even a very poor estimate, initialization with a gain that is close to the truth leads to even more rapid convergence. This can be seen in Fig. 4, which compares the number of iterations and the residuals while calibrating a time slice, first with no initial gain estimate (convergent calibration) and then with a gain estimate determined from the previous time slice (tracking calibration).

Finally, several heuristics have been put in place for handling exceptional data cases which can derail the tracking calibration due to the constant feedback of incongruous calibration solutions to future initial estimates. The tracking calibration algorithm can be summarized as follows:

• 1.

  Carry out a calibration to convergence over the first time slice ina chosen temporal window at the specified spectral resolution.This determines the baseline of parameter estimates.

• 2.

  Generate a mean over the short term history of solutions.

• 3.

  For every time and spectral bin within a window,

  • (a)

    Precalibrate incoming ACM with the mean solution.

  • (b)

    Estimate positions of model sources via DOA estimation techniques.

  • (c)

    Carry out a single iteration of all minimization procedures within the algorithm, using the previously generated baseline parameters as initial estimates.

  • (d)

    Compare updated parameter estimates for temporal and spectral smoothness. If a deviant solution is found; reject or else update solution history.

A computational profile of the various algorithmic components of both the convergent and the tracking calibration schemes is depicted in Table 5. The calibration solution memory pool is utilized by the goodness of solution block to decide on the propagation of an instantaneous solution.

5. Observational performance of tracking calibration

In this section, we present the results of applying tracking calibration on data taken under different observing conditions that an autonomous calibration algorithm might encounter, and evaluate the performance of the various estimators in those conditions. Our metric of performance is two-tiered: We first establish the calibration efficiency of our approach on a single time and frequency slice. During this test, we set high convergence thresholds in order to extract the best performance from the calibration parameter estimators. We refer to this estimator as the “convergent calibration” estimator.

We next compare the calibration efficiency of the tracking approach against the previously established best performance. Here, we carry out a single iteration of both the Major and Minor cycles, consequently providing the worst calibration performance, but with minimum latency and compute load.

Fig. 4 Typical major and minor cycle convergence rates on a single time slice.

Data description: due to the hardware per-dipole correlator still being under development, these data were obtained via a unique mode of operation of the LOFAR. Here, the beamformer weights of all antennas in a station were set to 0, except for one antenna, whose weight was set to 1. Passing the subbanded datastream through the beamformer then resulted in a single dipoles’s data appearing as a beamformed output. These were subsequently recorded to disk. The LOFAR hardware allows the continuous recording of 5 subbands from all 288 dual-polarized AARTFAAC antennas in this mode. These were subsequently correlated offline before further processing and analysis.

Data corresponding to the subbanded timeseries from each antenna of the AARTFAAC array were recorded between 1200–1500 UTC, 21 Sep. 2011, and again throughout the night of 12 July 2012. The stations were in the LBA_OUTER array configuration during all observations to facilitate comparisons. The salient details of the observations are tabulated in Table 4.

5.1. Calibration accuracy on a single time slice

In this section, we elaborate on the convergence behavior of the calibration on a typical data segment, while commenting on the aspects of weighting for model source flux estimation and gain constraints. A time slice from Obs. 3 is used for this analysis, and could be considered typical in terms of data quality and observing conditions. The fraction of data affected by RFI in our observations is minimal, consistent with the findings of (Offringa et al. 2013), and does not affect our conclusions. Figure 4 shows the typical number of major and minor cycle iterations taken for convergence via their convergence curves. Here, the stars indicate calibration with no prior knowledge of gain solutions, as used in convergent calibration, while the diamonds show the convergence rate after pre-initializing the solver with solutions from a previous time slice. Both the major and minor cycles are seen to benefit from the initialization. Their steep convergence rate is a motivation to apply tracking calibration for the AARTFAAC.

Fig. 5 Sky model efficacy evaluated on a single time slice. Left: effect of increasing sky model complexity. Right: phase behavior of residual visibilities after model subtraction.

Model efficacy: we utilize a simple point source sky model for self calibration. The efficacy of this model can be evaluated via the modeling error after calibration. The parameters of the brightest sources making up the sky model are given in Table 2, where the fractional flux contribution of the model sources to the total received flux per snapshot has been estimated using the 3CR catalog. In the left panel of Fig. 5, we plot the modeling error as a function of skymodel complexity, which is increased by the addition of a point source into the sky model in descending order of brightness. One can see that the addition of a larger number of weaker sources to the skymodel does not significantly decrease model error, indicating that a low complexity model is adequate. This is due to the brightness of the model sources, which contribute a significant fraction of the received flux, as seen in Table 2. The right panel of Fig. 5 shows the phase of the ratio between calibrated and model visibilities as a function of calibrated visibility power. The ratios should ideally have the response of a point source, that is, 0 phase across all visibilities. This is seen in the figure, along with an error contribution on all baselines due to unmodeled sources, which constitute the real information from the calibrated image. The large phase scatter at the low residual power end is contributed by visibilities with low S/N, thus having a higher calibration error. Visibilities with high S/N show <0.15 rad phase error, implying that the simple point source sky model containing a large fraction of received flux is adequate for snapshot calibration. This is typical for the AARTFAAC, unless RFI or other unmodeled sources have a large contribution.

Fig. 6 Number of iterations of major and minor during calibration of Obs. 1. The first of the dataset is affected by the flaring Sun, requiring larger number of major and minor iterations for convergence.

Figure 6 shows the longer term behavior of instantaneously calibrating every time slice of the timeseries from Obs. 1 to convergence. An average of close to 40 iterations are taken in the minor cycle (spread over the 3–4 major cycle iterations) in the initial part of the dataset, while this reduces to about 30 iterations in the latter half of the observation with about 3 major cycles. This small reduction in the number of minor cycles indicates the steep convergence of StEFCal, which results in a smaller number of minor iterations for the later major cycles due to a better initial estimate. The convergence criterion is based on the rate of change of estimated parameters reaching a (low) threshold rather than a limit on the fitting error. The slope obtained for the major cycle while calibrating a single time and frequency bin and that for the minor cycles during the final major cycle are also shown in the figure.

Fig. 7 AARTFAAC zenith projection of all sky images at ~60 MHz with 90 kHz and 1 s integration. Left: uncalibrated image (middle) after calibration and (right) after subtraction of the brightest sources.

Figure 7 shows an all-sky image with an integration of 1 s, as generated by the AARTFAAC calibration and imaging pipeline, along with the effect of calibrating the data and subtracting the brightest sources. The axes are in direction cosine units of a tangential projection of the celestial hemisphere on the zenith plane. The boundary of the circular image represents the local horizon. The image on the left shows the uncalibrated field of view dominated by Cas A, Cyg A, and the Galactic center at the southwest corner. The middle panel shows the same data after calibration. Note the coherency of the flux of the dominant point sources and the absence of the bright Galactic center due to spatial filtering. The right panel shows the effect of subtracting the brightest sources (the A-team) using their estimated model fluxes and positions. Note the appearance of secondary bright sources, corresponding to several known 3C sources.

5.2. Accuracy of tracking calibration

We present the performance of our algorithm on Obs. 3. The performance is tuned for maximum computational efficiency by restricting the major and minor cycles to just one. The tracking solutions are compared to those obtained by calibration to convergence. The comparison is carried out on estimated phase parameters, which affects image fidelity the most.

Figure 8 shows the mean absolute deviation of tracking calibration phase solutions from those obtained via convergent calibration, showing the worst case error for a single major and minor cycle iteration. The dataset is characterized by a disturbed ionosphere, which is probably due to the recombination of the disturbed plasma. This results in short-timescale amplitude and position jitter of the model sources. These are effectively tracked by the short cadence calibration, as is shown in Sect. 5.4. The tracking calibration solution error per parameter rises during certain times due to slight variations in model source position estimates between the tracking and convergent calibration, which are affected by the convergence error of tracking calibration.

Fig. 8 Error timeseries of tracking calibration solutions as compared to those obtained after iterating to convergence for Obs. 3.

5.3. Stability of calibration solutions

The stability of the tracking calibration can be gauged from the errors on the estimated parameters. The aggregate phase of a receiver path is more susceptible to disruption than the amplitude, primarily due to the rapid phase distortions introduced by the ionosphere during a calibration cycle. The errors on the estimated parameters are most likely caused by nonconvergence of the model fitting, or biases due to model fitting errors. The former can be due to a complete model mismatch because of the flaring Sun or intense RFI and are usually easily identified by exploiting temporal stability constraints. The latter can arise during observations with a significantly disturbed ionosphere, or due to low level RFI. These are more difficult to isolate, although they can be identified by referring to the typical temporal behavior of each parameter. To this end, Fig. 9 shows the estimated complex gain timeseries from a few antennas over ~3 h. The variance on the estimated phase of the instrument is ~5 degrees, while the amplitude has variations of 2% on average. In conjunction with Fig. 11, Fig. 9 demonstrates the resilience of the algorithm to observational perturbances while showing the segregation of sources of error into instrumental and observational.

Fig. 9 Temporal stability of the complex gain of a randomly chosen set of 6 antennas from the 288 AARTFAAC antennas over a period of ~3 h. The outliers are caused due to instances of significant RFI.

5.3.1. Lightcurve stability

An important aspect of calibration is the flux recovery of sources with known fluxes. This is now demonstrated for our calibration scheme. We carry out DFT imaging of data calibrated with the tracking calibration algorithm and extract out fluxes of some field sources in the observation. These are presented in Fig. 10. One can see that the dominant source of variation in the extracted flux is due to systematic variations on short timescales, which is most probably caused by the ionosphere. This limits the estimation accuracy on the extracted lightcurve.

The position offset of the extracted sources as compared to a cataloged position shows no significant bias, leading to the conclusion that there are no significant systematic calibration errors.

5.4. Effect of the ionosphere on AARTFAAC calibration

The expected strong phase coherence of a plane wave across the AARTFAAC array due to its very short baselines can be mitigated by the ionosphere that introduces spatially variant refraction and propagation delays, which can require rapid calibration cycles to over sample ionospheric phase coherency timescales (see Sect. 3.1).

Recent results (Intema et al. 2009) have indicated the presence of a turbulent layer below the peak TEC, which has more power in the smaller scale fluctuations than in the case of pure Kolmogorov turbulence. This is exacerbated by the very large fields of view of the AARTFAAC. These are expected to be larger than the typical scale size of ionospheric fluctuations when projected onto ionospheric heights, resulting in anisoplanatic conditions where the ionospheric phase error varies over the field of view of each antenna.

Note that short baselines between closely spaced stations were expected to be immune to ionospheric effects and are expected to play a significant role in constraining the calibration of the full LOFAR instrument (van der Tol et al. 2007).

During our test observations, the most common effects seen seem to be due to small scales “bubbles” of turbulent plasma, which affect identical sightlines from different antennas. This results in high-frequency scintillation of the brightest sources, and position wander. Figure 11 shows the effect of the ionosphere on the estimated flux of the model sources and the position wander of Cas A over the course of a night.

The variation of the position of model sources from catalog locations under a typical ionosphere is demonstrated by Fig. 11d. Here, the actual position estimates have been obtained using the WSF algorithm. Systematic variations of up to 20% of a PSF are present and are seen to be a strong function of frequency, as expected due to the ionosphere.

Tracking such behavior is the primary reason for carrying out near real-time calibration, more so because the streaming nature of the application precludes multiple calibration passes through the recorded visibilities. Tracking calibration is able to recover such short term effects quite successfully, as can be seen from Fig. 11.

Fig. 10 Lightcurves of a few bright field sources from images of Obs. 3 of the calibrated visibilities.

Fig. 11 Short term ionospheric effects in the form of amplitude scintillation and position wander being successfully recovered during calibration. Panels a–c show different levels of scintillation on Cas A and Cyg A during trial observations. Note the rising Sun in panel c) Panel d shows the deviation of the WSF estimated position of 3C461 during obs. 2 over ~1000 s, as compared to its catalog position. Shown at two frequencies separated by ~30 MHz.

This direction dependent correction is estimated and carried out only toward the directions of a few model sources. Thus, the instantaneous gains and model source fluxes and positions are estimated such that the flux and positions of model sources are consistent. This does not apply to the other sources within the field of view. Hence, we have observed that the lightcurves of field sources which are only a few degrees away from the model source, can have high frequency flux fluctuations contributed due to the ionosphere. From the test observations, this places another constraint on the achievable accuracy of lightcurves extracted from calibrated images.

5.5. Enhancing the transient detection sensitivity via difference imaging

The snapshot dynamic range of our calibrated images can be lower than expected due to several reasons. Observational conditions like a flaring Sun or significant ionospheric effects can only be addressed in a limited manner due to latency and computing constraints. The effect of spatial filtering of diffuse emission via visibility tapers for ease of modeling can give rise to beam sidelobes ripples from unmodeled flux. Carrying out only a few iterations while calibrating can lead to partial convergence of a snapshot calibration. This, in turn, can also lead to higher sidelobe confusion noise. Further, the ASM is expected to hit its classical confusion limit in a few tens of seconds of integration due to its poor resolution. At short timescales, all these effects are systematic due to the negligible change in the sky brightness distribution and instrumental parameters and hence, do not average down on integration.

Thus, the ASM is expected to be a confusion-noise limited instrument. Here, we examine the temporal difference image domain as a viable alternative to snapshot images for transient searches and quantify the improvement in sensitivity of such images.

A difference image is expected to have a higher sensitivity due to the cancellation of the above mentioned systematic effects with the residual noise expected to better follow Gaussian statistics. Further, steady sources are also cancelled in the difference domain, making it an excellent domain for searches for short-term transients. The sensitivity of the difference image depends on the extent to which systematics cancel, which in turn can be influenced by observing conditions. The effectiveness of this domain is demonstrated in Fig. 12, which shows the decrease in noise with integration time for a snapshot and as a difference image, as a function of integration time. The difference image with N units of time integration is generated as: (5)Here, all operations are carried out per pixel, and the noise is measured in a source free region of the image. The theoretical curve shows the expected enhancement on the integration of N timeslices, and we see the noise in a difference image that closely following this curve out to about ~30 s, beyond which is when residuals is when start getting affected by earth rotation. In principle, short-term image differences (~s) are enough to reduce noise amplitude by 20–30%.

6. Pipeline architecture and computational performance

Fig. 12 Reduction in image noise as a function of integration time over raw and subtracted consecutive image pairs. Each member of the pair being subtracted is generated by averaging raw images over increasing cadences.

This section describes the implementation of the ASM and its performance. The upstream correlator implementation generates a timeseries of ACMs per frequency channel. These are continuous within a subband of 192 kHz, while the subbands may be spread over a large frequency span. For a single time slice, the low coupling between calibration of channels suggests parallelizing along the frequency axis with a subsequent collation for creating a frequency collapsed image. Tracking calibration requires a per-channel, conditioned calibration solution for initialization, which requires maintaining this state information per channel.

We thus utilize a two level parallelization scheme. The first stage splits the visibilities on each subband to be processed by a single pipeline, as shown in Fig. 13. The second stage splits each subband into channels, which are individually flagged and calibrated, as depicted by the stacked rectangles within each pipeline. Finally, all channels within a pipeline are imaged and combined into an image cube for further analysis.

Fig. 13 High level schematic of the calibration and imaging pipeline design depicting the two stages of parallelism on the subband and channel axis.

The iterative and sequential nature of the estimators operating on an individual time/frequency bin is better suited for a CPU, rather than a single algorithm multiple data architecture like a GPU. We thus find that an implementation on commodity multicore CPUs is matched for our application.

We utilize the Pelican (Salvini et al. 2011) framework for data distribution, in where a server process provides high throughput data distribution to a collection of worker threads on a first-come, first served basis. This scheduling maintains a minimal state between threads; however, it does not guarantee delivery of data from a given spectral channel to the same thread. We thus utilize the Pelican service data channel (indicated by the dashed lines in Fig. 13) to implement a server for distribution of the low bandwidth calibration solutions to the corresponding channels. We implement the individual estimators using the Eigen3 (Guennebaud et al. 2010) linear algebra library’s numerical solvers and matrix operator implementations. Errors and statistics from each pipeline are also reported using the Pelican service data channel and can be tapped by network based utilities for monitoring and logging. To minimize data copies between per CPU core caches, we tie a single core to an individual pipeline implementation.

6.1. Performance results

The system processes a channel per CPU core, or per thread. We hence define the performance as the processing time of a single channel. Table 5 gives an overview of the time spent on the various calibration functions. Note that these timings were obtained for data that have a calibration solution, and were obtained for the worst case computational case of calibrating to convergence. When no solution is found, the timings depend on the maximum number of iterations allowed in both the minor and major cycle and the WSF algorithm.

The computational cost is dominated by the model source position estimation, while an individual iteration of DIE calibration is seen to be quite low. We measured the total time (flagging, calibration, and imaging) of a subband with N channels using N cpus to be 250 ms (σ = 15 ms) for nighttime data and 1700 ms (σ = 100 ms) for daytime data on a Intel Core i7-2600 CPU at 3.40 GHz, which is similar to our target ASM hardware.

7. Observational challenges to tracking calibration

This section discusses some observational aspects seen in test data, which can adversely affect the tracking calibration scheme, and our approach to handling them. There are two categories of challenges faced by our approach to real-time calibration. Firstly the instantaneous calibration can become biased due to the presence of unmodeled signal components, or due to a larger error on the estimation parameters on data with poor S/N because of limited number of iterations. These effects tend to lower the sensitivity of each time slice. However, the TraP has several checks to avoid being biased by these inconsistencies. Secondly due to the feedback in tracking calibration introduced via the initial estimation, biases in the solutions can be propagated in time. This can lead to unintended variations in the extracted flux due to under/over estimation of the gain amplitudes and can also lead to a higher deconvolution noise due to position errors on the brightest sources. This effect is more challenging to counter due to the minuteness of the bias but can be countered effectively by introducing periodic convergent calibration cycles in the real-time stream. We first illustrate the former set of challenges.

7.1. Radio frequency interference (RFI)

RFI mimics a coherent source across the array and is usually very bright and sporadic in both time and frequency. It can thus be very detrimental to a transient machine due to the difficulty of its separation from genuine celestial transients with a corresponding increased false detection rate. Strong terrestrial RFI near the AARTFAAC site has been found to affect only a small fraction of data. The high time and spectral resolution of calibration allows RFI rejection to be carried out within our pipelines via simple sliding window filters operating on the total received power per ACM. However, low power RFI usually cannot be detected by these, resulting in a model versus data inconsistency. During times of reasonably heavy RFI, the estimation does not converge and so is easily distinguished. For tracking calibration, the steep convergence curve usually generates solutions for these time slices, which are significantly different from the maintained time history, and the time slice is effectively rejected by the filters on the generated solutions.

The RFI from a single source usually results in the estimation placing the brightest source in the model sky at the location of the RFI source. This results in significantly different phase solutions, which can also be trapped. This also includes effects due to lightning. An insidious source of RFI are reflections of terrestrial RFI (which is usually at low elevations, and hence trappable) from meteor trails or passing aircraft. These can be dim enough to not bias the skymodel significantly, which makes them difficult to filter out. They can, however, be trapped by higher level filters due to their rapid movement across the field of view.

In conclusion, bright RFI should not be a significant source of false positives for the AARTFAAC. The fast recalibration allows us to distinguish individual RFI corrupted time slices and discard them from further analysis, while weak RFI, which does not affect calibration, has to be filtered by higher layers like the TraP.

7.2. The active Sun

The quiet Sun is a well-known broadband radio source at frequencies of tens of MHz due to black-body radiation from the Corona being in the Rayleigh regime. In spite of the radio extent of the quiet Sun being much larger than the optical disk, simple point or Gaussian models are adequate for its representation with the AARTFAAC PSF. However, during periods of sunspot activity, small regions of emission can extend over several solar radii. These typically have an inverted spectrum, making them the brightest sources in the low frequency sky. During intense activity, these components can be resolved by the AARTFAAC, thus invalidating the point source assumption of the solar model. A multicomponent point source solar model, which varies with time and frequency, would then be required for adequate modeling the solar contribution to observed visibilities. This model can be extracted from a snapshot via high resolution beamforming toward the Sun and the calibration model updated. However, our latency bounds preclude such an operation currently. Thus, time slices with a flaring Sun are either discarded using filters based on the ACM total power or lead to images with significantly higher noise, which are appropriately handled by higher level processors. The dynamic range of AARTFAAC images after tracking calibration is ~2200:1, so images in the presence of the quiet Sun can be considered for transient detection.

7.3. Creeping bias in tracking calibration

Since both the positions of the model sources and the phases of the antenna gains are estimated during the calibration, a biased estimate of one set of parameters can be compensated by a corresponding change in the other during the cost function minimization while estimating calibration parameters. This is especially seen in tracking calibration due to the use of solutions from a previous time slice as an initial estimate, which can be biased due to nonconvergence of tracking calibration on data with poor S/N. A simple solution for this has been found to be the periodic recalibration to full convergence of the incoming data every few minutes. This limits the bias to a small fraction of the estimates. The effect of this can be seen on the phases in Fig. 8, which is subsequently brought back on track via a convergent calibration.

8. Conclusions

The transient detection figure of merit of a fast-transient monitor is a function of its field of view, response time, and sensitivity. At low frequencies, the latter crucially depends on adequate calibration of instrumental and various observing parameters over the short detection timescales.

In this paper, we have described a real-time, wide-field calibration strategy optimized for the detection of such transient signals at low frequencies. We utilize a model-based, direction-dependent multisource self-calibration algorithm with high calibration efficiency and demonstrate that parameter temporal tracking can be utilized to achieve bounded latency. We describe the implementation of the calibration algorithm and demonstrate this on test data taken under typical observing conditions that latency and compute constraints for real-time operation are met in the computationally worst case scenario of calibration to convergence. We utilize the expected slow variations of the instrumental parameters to constrain the self calibration, so it is immune to the rapid flux and position variations of model sources in the presence of a diffractive ionosphere.

Among contemporary low frequency arrays, the AARTFAAC array is well suited to perform as a real-time monitor for bright and short radio transients.The AARTFAAC tracking calibration is optimized for the precise 2D nature of the ASM, its snapshot imaging mode of operation, and its sensitivity and field of view. Other established calibration algorithms, which is typically optimized for earth rotation synthesis, apply poorly to the simplified constraints of the ASM and do not satisfy the latency and compute requirements of the real-time monitor. We have shown under a variety of conditions that both RFI and ionospheric visibility distortions do not significantly affect AARTFAAC’s sensitivity due to its rapid snapshot observing mode. Instead, its sensitivity is limited by classical confusion noise caused by its relatively coarse resolution. We have further shown that the confusion noise limit can be effectively breached via difference imaging over short timescales (tens of seconds) and that difference images approaching the thermal limit can be generated. This makes the instrument more sensitive to short term transients than quiescent sources.

To operate autonomously, effective filtering of inconsistent calibration solutions from the tracking loop is necessary under challenging observing conditions due to RFI, an active Sun, and ionospheric effects. This is done by careful selection of thresholds on the relative error on the timeseries of calibration solutions while feeding the tracking loop. Coupled with a “coasting” strategy, this has been found to be effective for real-time, autonomous operation. With a current C++ implementation, the tracking calibration operates within ~0.25 s on a single time and frequency time slice of a typical observation, which is well within the compute budget.

Thus, tracking calibration is an appropriate strategy for an autonomous, real-time instrument like the AARTFAAC all-sky monitor.

Acknowledgments

This work is funded by the ERC grant C.2320.0056 awarded to Prof. Ralph Wijers, University of Amsterdam. We thank the LOFAR team for building an extremely configurable instrument, allowing the collection of dipole level raw data for commissioning observations for the AARTFAAC. We thank The Netherlands Institute for Radio Astronomy (ASTRON) for support provided in carrying out the commissioning observations. We thank John Swinbank, Alexander van der Horst and Antonia Rowlinson for helpful comments. Finally, we thank the anonymous referee for their comments, which allowed us to clarify a few aspects and generate a more readable paper.

References

All Tables

All Figures

	Fig. 1 AARTFAAC calibration and imaging pipeline within the full ASM flow and leading onto the transient detection pipeline (TraP).
In the text

	Fig. 2 Block diagram depicting the snapshot iterative calibration control flow in the AARTFAAC optimal calibration scheme.
In the text

	Fig. 3 Schematic view of tracking calibration for bounded latency precise calibration.
In the text

	Fig. 4 Typical major and minor cycle convergence rates on a single time slice.
In the text

	Fig. 5 Sky model efficacy evaluated on a single time slice. Left: effect of increasing sky model complexity. Right: phase behavior of residual visibilities after model subtraction.
In the text

	Fig. 6 Number of iterations of major and minor during calibration of Obs. 1. The first of the dataset is affected by the flaring Sun, requiring larger number of major and minor iterations for convergence.
In the text

	Fig. 7 AARTFAAC zenith projection of all sky images at ~60 MHz with 90 kHz and 1 s integration. Left: uncalibrated image (middle) after calibration and (right) after subtraction of the brightest sources.
In the text

	Fig. 8 Error timeseries of tracking calibration solutions as compared to those obtained after iterating to convergence for Obs. 3.
In the text

	Fig. 9 Temporal stability of the complex gain of a randomly chosen set of 6 antennas from the 288 AARTFAAC antennas over a period of ~3 h. The outliers are caused due to instances of significant RFI.
In the text

	Fig. 10 Lightcurves of a few bright field sources from images of Obs. 3 of the calibrated visibilities.
In the text

Fig. 11 Short term ionospheric effects in the form of amplitude scintillation and position wander being successfully recovered during calibration. Panels a–c show different levels of scintillation on Cas A and Cyg A during trial observations. Note the rising Sun in panel c) Panel d shows the deviation of the WSF estimated position of 3C461 during obs. 2 over ~1000 s, as compared to its catalog position. Shown at two frequencies separated by ~30 MHz.

In the text

	Fig. 12 Reduction in image noise as a function of integration time over raw and subtracted consecutive image pairs. Each member of the pair being subtracted is generated by averaging raw images over increasing cadences.
In the text

	Fig. 13 High level schematic of the calibration and imaging pipeline design depicting the two stages of parallelism on the subband and channel axis.
In the text