aartfaac_RT_calib_v1.tex~

% \documentclass[referee]{aa}
\documentclass{aa}
\usepackage[varg]{txfonts}
\usepackage{color}
\usepackage{graphicx}
%% \usepackage{subfig} %% For generating a grid of figures.
\usepackage{caption}
\usepackage{subcaption}
%% \usepackage{float}
%% \usepackage{epstopdf}
%% \newsavebox\mybox   %% For generating figure captions which wrap at figure edges.
%% \newlength\myboxlen
%% \newcommand{\figcap}[2]
%% {%
%%   \sbox\mybox{#1}
%%   \settowidth{\myboxlen}{\usebox{\mybox}}
%%   \centering
%%   \usebox\mybox
%%   \hskip \textwidth
%%   \parbox{\myboxlen}{#2}
%% }

\usepackage{bm}
\bibpunct{(}{)}{;}{a}{}{,} % to follow the A&A style

%% From Wijnholds, for better figure placement
\renewcommand\floatpagefraction{1.0}
\renewcommand\topfraction{1.0}
\renewcommand\bottomfraction{1.0}
\renewcommand\textfraction{0.0}

\begin{document}

\title{Real-time calibration of the AARTFAAC array for transient detection}

\author{P. Prasad    \inst{1,2} 
 \and   S. J. Wijnholds  \inst{3}  
 \and   F. Huizinga    \inst{4}  
 \and    R.A.M.J Wijers    \inst{5}}  
 \institute 
 {Anton    Pannekoek    Institute,  University of  Amsterdam, 
  Postbus  94249, 1090 GE  Amsterdam, The Netherlands\email{prasad@astron.nl} 
\and
  ASTRON, The  Netherlands Foundation for  Radio Astronomy, 
  PO Box 2,  7990 AA Dwingeloo, The  Netherlands\email{prasad@astron.nl}
\and
  ASTRON, The  Netherlands Foundation for  Radio Astronomy, 
  PO Box 2,  7990 AA Dwingeloo, The  Netherlands\email{wijnholds@astron.nl}
\and
  Anton    Pannekoek    Institute,  University of  Amsterdam, 
  Postbus  94249, 1090 GE  Amsterdam, The Netherlands \email{F.Huizinga@uva.nl}
\and
  Anton    Pannekoek    Institute,  University of  Amsterdam, 
  Postbus  94249, 1090 GE  Amsterdam, The Netherlands \email{r.a.m.j.wijers@uva.nl}
}

%\author{P. Prasad    \inst{1,2} \email{prasad@astron.nl}
% \and   S.J. Wijnholds  \inst{2}  
% \and   F. Huizinga    \inst{1}  
% \and    R.A.M.J Wijers    \inst{1}}  
% \institute 
% {Anton    Pannekoek    Institute,  University of  Amsterdam, 
%  Postbus  94249, 1090 GE  Amsterdam, The Netherlands \label{1}
%\and
%  ASTRON, The  Netherlands Foundation for  Radio Astronomy, 
%  PO Box 2,  7990 AA Dwingeloo, The  Netherlands \label{2}
%}

\date{Received <date> / Accepted <date>}

\abstract{ The  search for transient phenomena  at low radio  frequencies is now
  coming of age with the development  of  radio sky monitors with a large field of view, which are
  made  feasible by new  developments in  calibration algorithms  and computing.
  However, accurate calibration of such arrays is challenging, especially within
  the latency requirements of near real-time transient monitors, and is the main
  cause of  limiting their sensitivities.   This paper describes a  strategy for
  real-time, wide-field direction  dependent calibration of the Amsterdam-ASTRON
  Radio Transients  Facility and  Analysis Center (AARTFAAC)  array, which  is a
  sensitive, continuously  available all-sky monitor based on  the Low Frequency
  Array (LOFAR). The monitor operates in a zenith pointing, snapshot imaging
  mode for  image plane detection  of bright radio  transients.  We show  that a
  tracking calibration approach with solution propagation satisfies our latency,
  computing,   and  calibration  accuracy   constraints.   We   characterize  the
  instrument and  verify the calibration  strategy under a variety  of observing
  conditions.  This brings out several phenomena which can bias the calibration.
  The real-time  nature of  the application further  imposes strict  latency and
  computational  constraints.  We  find that  although ionosphere  induced phase
  errors  present a  major  impediment  to accurate  calibration,  these can  be
  corrected  in the  direction of  the  brightest few  sources to  significantly
  improve  image  quality.  Our  real-time  calibration pipeline  implementation
  processes a  single spectral channel of  a snapshot observation  in $\sim$ 0.2
  seconds  on  test hardware, which is well  within  its  latency budget.   Autonomously
  calibrating and imaging one second  snapshots, our approach leads to a typical
  image noise  of $\sim$  10 Jy for  a $\sim$  90 kHz channel,  reaching dynamic
  ranges of $\sim$ 2000:1.  We  also show that difference imaging allows thermal
  noise  limited transient  detection,  despite the  instrument being  confusion
  noise limited.}

\keywords{Instrumentation:Interferometers - telescopes - atmospheric effects - methods: observational}

\maketitle

\section{\label{sec:Introduction}Introduction}
The detection  and characterization of low  frequency radio transients  is a Key
Science  Project \citep{fender2006lofar}  for  the Low  Frequency Array  (LOFAR)
\textbf{\citep{vanhaarlem2013lofar}}  radio  telescope.   In this  context,  the
Amsterdam-ASTRON    Radio    Transients    Facility    and    Analysis    Center
(AARTFAAC) \footnote{see  www.aartfaac.org }  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  ($\sim$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   \ref{tab:AARTFAAC_specs}  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 $\gamma$-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 \citep{bower2007submillijansky,  bower2011search}, follow-up studies of
discovered  transients \citep{chandra2012radio}, and  a few  dedicated observing
programs  \citep{bannister2012limits,  katz2003survey}.   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) \citep{lonsdale2009murchison}  and  the Long  Wavelength
Array (LWA) \citep{ellingson2009long}.

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  \citep{lorimer2007bright,  thornton2013population}  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
\mbox{\citep{wijers2014aartfaac}}.
\begin{table*}[tbhp]
\caption{Specifications of contemporary radio transient detection machines.}
\label{tab:AARTFAAC_specs}
\center{%
\begin{tabular}{c c c c c}
\hline \hline
Parameter                   & AARTFAAC LBA            & MWA        & LWA1 & Comment\\ 
\hline 
Array elements              & 288 inverted V antennas & 128 tiles  & 256 wiregrid bowtie & Dual polarized elements\\ 
Frequency range (MHz)       & 30-80                   & 80-300     & 10-88 & \\ 
Field of view (sr)          & $\pi$                   & 0.06$\pi$  & $\pi$ & FWHM of beam\\ 
%% Effective area ($m^{2})$    & 10                      & 21         & 2.5   & per element\\ 
Total Effective area ($m^{2})$    & 2617\tablefootmark{a} & 2752\tablefootmark{b}         & 1393\tablefootmark{c}   & \\ 
%% $T_{sys}$($\nu^{-2.55}$K)              & $\sim2000$              & $\sim350$  & $\sim1500$ & \\
$T_{sys}$($\nu^{-2.55}$K)              & ~3600              & ~350  & ~2100 & \\
Angular resolution (arcmin) & ~60              & ~3  & ~120 & \\ 
Spectral resolution (kHz)   & 16                      & 40         & 75       & \\ 
Processed Bandwidth (MHz)   & 13                      & 30.72      & 0.075    & LWA: Narrow band Imaging \\
Temporal resolution (sec)   & 1                       & 8          & 5        & \\ 
%% Transient detection FoM $(m^2sr K^{-1})$ & 2.3        & 1.3        & 0.830    & Per unit BW and unit time \\ \hline 
Transient detection FoM\tablefootmark{d} $(m^2sr K^{-1})$ & 2.28        & 1.48        & 2.08    & Per unit BW and unit time \\ \hline 
\end{tabular}}

\tablefoot{The  parameters are  estimated  at the  highest  sensitivity of  each
  instrument, viz. 60 MHz  for AARTFAAC, 150 MHz for the MWA,  and 74 MHz for the
  LWA, which  observes continuously only in  its 'Narrow band'  mode. The zenith
  pointing $A_{eff}$ is quoted. The $T_{sys}$ ( assumed to be sky noise dominated)
  is  also evaluated  at the  same frequencies from \citep{bregman2000concept}.
  \tablefoottext{a}{From    \citep{wijnholds2011situ}.}    
  \tablefoottext{b}{From    \citep{ellingsonLWA1}.} 
  \tablefoottext{c}{From \citep{tingay2013murchison}.} 
  \tablefoottext{d}{See Sec. \ref{sec:The-AARTFAAC-All} for FoM definition}.}
\end{table*}


\textbf{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^{12}$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  fact 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 \citet{law2012all}.   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 SNR 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 all-sky monitors  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  \citep  {nemiroff2003tile}.  Imaging  all-sky
monitors are hence matched for this application.


\textbf{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   Section
\ref{sec:The-AARTFAAC-All}. We lay out  the problem of wide-field calibration in
the  context of transient  detection in  Section \ref{sec:Array-calibration-for}
while describing  the approaches taken by  similar instruments.  We  then give a
detailed   description    of   our   approach   to    calibration   in   Section
\ref{sec:An-Optimal,-tracking},   especially  in   the  context   of  autonomous
operation.   Section  \ref{sec:Performance-of-tracking}  presents  the  observed
performance  of the calibration  approach on  commissioning data,  while Section
\ref{sec:Computational-performance} elaborates  on the computational performance
of our strategy.  Section \ref{sec:Challenges-to-tracking} elaborates on some of
the challenges faced  by the AARTFAAC ASM under  real observing conditions after
which we present  our conclusions. \textbf {The results  presented in this paper
  have been  obtained on test observations  carried out by  using existing LOFAR
  system          infrastructure,           while          the          AARTFAAC
  Uniboard    \footnote{http://www.radionet-eu.org/uniboard}   based   piggyback
  recording system was being built.}


\section{\label{sec:The-AARTFAAC-All}The AARTFAAC All Sky Monitor System
overview} 

\begin{figure*}[htbp]
\centering
\includegraphics[width=1\textwidth]{Figs/AFAAC_blkdia_latest.png}
\caption{AARTFAAC  calibration  and imaging  pipeline  within  the full  all-sky
  monitor flow and leading onto the transient detection pipeline (TraP).}
\label{fig:The-AARTFAAC-calibration}
\end{figure*}


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 $\sim$$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
$\sim$$100$  MHz sampled band  into $\sim$$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   \footnote{http://www.radionet-eu.org/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  PSF.   The estimation  of  the  resulting $\sim$$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 $\sim$$66$
subbands ($\sim$$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 4x4
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 $\sim$$20$ sq. deg,
and the  arbitrary pointing of the  analog beam of  this system, its use  is not
discussed further in this paper.

Table   \ref{tab:AARTFAAC_specs}  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
\ref{fig:The-AARTFAAC-calibration}  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.


\subsection{\label{sub:AARTFAAC-for-Transient}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
$\sim$0.1 to  1 second 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
\citet{bhat2004multifrequency}. This implies a DM upper limit of $\sim$100 at 60
MHz when  lines of sight  are within the  Galaxy and when scatter  broadening is
restricted to 1 second.

The temporal pulse  broadening in microseconds over a  bandwidth $\Delta\nu$ MHz
at $\nu$ GHz for a dispersion of DM units is given by Eq. \ref{eq:dm_broad}.
\begin{equation}
\Delta t_{DM}=8.3\,DM\,\Delta\nu_{MHz}\,\nu_{GHz}^{-3}\, \mu sec\label{eq:dm_broad}
\end{equation}
Imaging at the  highest available spectral resolution of  16 kHz and restricting
dispersion broadening to within 1 second  at the same time allows detection of a
source with  a DM of up to  $\sim$1600.  This allows probing  along any galactic
line of sight \citep{cordes2002ne2001}, although  at a sensitivity a factor that
is $\sim$28 poorer than that with the full 13 MHz band.

%% Assuming all subbands  are chosen to be contiguous  around the peak instrumental
%% sensitivity while the  dispersion broadening is limited to 1  second, the 13 MHz
%% continuum sensitivity of  the AARTFAAC would limit detection  of transients with
%% DMs of only a few.  However, restricting dispersion smearing to \mbox{<1} second
%% under  the upper limit  of DM  100 imposed  by scattering  would allow  a higher
%% spectral integration, before the sum  of the two broadening processes exceed the
%% imaging   cadence.   

\textbf{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  SNR for any pulses  with DMs greater
  than  a few.  This  is because  dispersion broadening  would spread  the pulse
  beyond the 1 second 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 second. 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
\citet{hassall2013detecting},  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  \citep{cordes2004dynamic} from  a
telescope  is  $A\Omega\left(\frac{T}{\Delta  t}\right)$  and can  be  used  for
comparison between instruments.  Here, A  is the total collecting area; $\Omega$
is the solid angle coverage; T is the total duration of observation; and $\Delta
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\Omega/T_{sys}$.  Among  currently  planned ASMs,  the  AARTFAAC's FoM  is
comparable  to   the  LWA,   which  has  similar   field  of  view   (see  Table
\ref{tab:AARTFAAC_specs}). 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  Section
\ref{sub:AARTFAAC-for-Transient}.

\citet{backer1999pers} has  shown that the  signal-to-noise ratio achieved  in a
continuous  wide-field  search scales  with  the  array  filling factor  $f$  as
$\sqrt{f}$. 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 $\sim$$58$  MHz, the PSF has
high roll-off,  low side  lobes $\sim$$15$  dB below peak  with a  resolution of
$\sim$$0.8$ sq.  deg.   \textbf{Assuming the confusion limit to  be reached with
  the  presence  of  one  source  per  ten  PSFs  and  the  source  counts  from
  \citet{bregman2012system},   results  in  a   classical  confusion   noise  of
  $\sim$$10$ Jy.}   The large number of  antennas also make the  AARTFAAC a very
sensitive instrument  with a  combined $A_{eff}/T_{sys}$ of  $\sim$0.7 $m^{2}/K$
\citep{wijnholds2011situ} and  an instantaneous thermal noise of  $\sim$8 Jy for
the standard imaging mode of 1sec/16  kHz.  This implies that the sensitivity of
typical snapshot images are limited by confusion noise and not thermal noise.

The  array is  2-D to  high accuracy  ($\sim$1.4 cm  rms in  the  z-axis); thus,
imaging  the  zenith  region  can  be  done  via  a  simple  \mbox{2-D}  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
the  spatial coherence sampling  function that  is confined  to a  plane passing
through  the origin,  keeping the  2-D relationship  between the  visibility and
imaging   domains  \citep{cornwell1992radio}.    Thus,  wide-field   imaging  is
simplified  as  no 3-D  transforms  or  their  approximations like  W-projection
\textbf{\citep{cornwell2008noncoplanar}} 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 all-sky monitor.

\section{\label{sec:Array-calibration-for}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
\citep{taylor1999synthesis}.   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
light curves  of existing sources, or  by detecting new sources  referenced to a
database   of   lightcurves.     The   Transient   Detection   pipeline   (TraP)
\citep{swinbank2014trap}  component  of  the  AARTFAAC ASM  carries  out  source
extraction  from individual  time slices via  thresholding and  fitting \mbox{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.

\begin{figure*}[tbh]
\centering
\includegraphics[width=0.7\textwidth]{Figs/Fig2_blkdiag.png}
\caption{Block diagram depicting the snapshot
iterative calibration control flow in the AARTFAAC optimal calibration
scheme. }
\label{fig:Block-diagram-depicting}
\end{figure*}

\subsection{\label{sub:All-Sky-cal-iono}Ionospheric effects}

The  antenna complex  gain has  contributions primarily  from the  gains  of the
electronics, receiver characteristics, and telescope geometry. These usually vary
slowly ($\sim$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
$\left(\sim10x\right)$  in large scale  Total Electron  Content (TEC)  with slow
temporal variations,  as well as  fluctuations on relatively small  temporal and
spatial scales  ($\sim0.1\%$ variations over a  few kilometers and  over tens of
secs).   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
\citep{wijnholds2010calibration}.

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.

\subsection{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.  \textbf{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
\citep{noordam2004peel,vdTol2007selfcallofar}, in which self-calibration toward
sources within  the field of  view is carried  out in decreasing order  of their
brightnes,  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  \citep{boonstra2003gain}.   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  use
mostly Peeling  based approaches. 

The  MWA   \citep{lonsdale2009murchison},  located  in   Western  Australia,  has
developed an  algorithm for carrying out real-time,  full polarization direction-
dependent  calibration of snapshots  of $\sim$8 seconds  via a  Peeling mechanism
\citep{mitchell2008real}.   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  \citep{ellingsonLWA1}  is  co-located  with  the  VLA  in  New  Mexico
(USA). Along with forming beams in real-time (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  \citep{cottona2004beyond}  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  minutes, 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  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. % , as can be seen in Figure \ref{fig:Light-curves-of}.

\section{\label{sec:An-Optimal,-tracking}An Optimal tracking calibration
scheme for transient detection}

Here,  we present our  approach for  latency bound,  precise calibration  of the
AARTFAAC  with  a  low  computational  footprint.  We  use  a  model-sky  based,
multi-source   self   calibration   approach   \citep{wijnholds2009multisource}.
Calibration of an  array of identical, wide-field antennas  entails deriving the
direction  dependent response  of every  antenna, each  of which  is  excited by
multiple,  simultaneously   present  calibration  sources.    The  system  noise
contribution   is  also   estimated,  and   subtracted  while   calibrating  the
visibilities.\\    \textbf{Algorithmic    Summary:}    As   seen    in    Figure
\ref{fig:Block-diagram-depicting},  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 direction independent effects.
\begin {enumerate}
\item {The  input to the algorithm  is a running timeseries  of measured spatial
  correlations  between  all  antennas  (visibilities organized  into  an  Array
  Covariance  Matrix  (ACM)) for  a  given time  and  frequency  slice, and  the
  corresponding instantaneous sky-model.}
\item {The Minor cycle within each Major cycle iteratively estimates the complex
  gains        of       all        antennas       using        the       StEFCal
  \citep{salvini2014stef,salvini2014stefursi}   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.}
\item {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.}
\item {The time  axis is used to feed forward parameter  estimates as an initial
  guess for the next timeslice, leading to faster convergence.}
\item {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.}
\end{enumerate}

 In the  following,  we first  describe  the calibration  of an  ACM
corresponding to  a single  timeslice, 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, e.g.,  via CLEAN
components. This is beneficial for a streaming calibration approach.

Note that 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
\citep{zatman1998narrow}.  This is appropriate as  the maximum time of flight on
the longest baseline of 300 m  (\textasciitilde{}1 $\mu s$) is much less than
the  inverse of  the typical  calibration channel  bandwidth of  $16$  kHz.



\subsection{Instantaneous calibration}

Array  calibration  can  be  formulated  as a  non-linear  parameter  estimation
problem, which can be  solved using statistically efficient estimators generated
via a  Non-linear Least Squares approach.  Our approach to  calibrating a single
incoming  timeslice closely  follows  that of  \citet{wijnholds2009multisource},
which establishes  asymptotically efficient estimators for  the model parameters
via a  Maximum-Likelihood (ML)  formulation. Given a  set of $p$  antennas whose
locations are known  precisely and $q$ calibration sources  with nominally known
positions and fluxes within the field  of view, the data model developed for the
instantaneous ACM $\mathbf{R}$ in their scheme is given by:
\begin{equation}
\mathbf{R=GAG_{0}\Sigma_{s}G_{0}^{H}A^{H}G^{H}+\Sigma_{n}}\label{eq:datamodel}
\end{equation}
where    \textbf{$\mathbf{G=}diag(\left[g_{1},\ldots,g_{p}\right])$    }is   the
diagonal matrix of DIE gain of each antenna, $\mathbf{A}$ is the matrix of array
steering vectors,  generated using  the known positions  of the  array antennas,
$\mathbf{G_{0}=}diag(\left[g_{01},\ldots,g_{0q}\right])$ is  the diagonal matrix
of   DDE   gain    in   the   direction   of   each    of   the   $q$   sources,
$\mathbf{\Sigma_{s}}=diag\left(\left[\sigma_{s1}^{2},\ldots,\sigma_{sq}^{2}\right]\right)$
is the known flux of the  sources in the sky model, and $\mathbf{\Sigma_{n}}$ is
the diagonal noise covariance matrix assuming an uncorrelated, but non-identical
system noise contribution from each antenna.

The model parameters to solve for are of the form:
\begin{equation}
\bm{\theta}=[\gamma_{1},\ldots\gamma_{p},\phi_{2}...,\phi_{p,}\sigma_{s1}^{2},...\sigma_{sq}^{2},\mathbf{l_{1}}^{T},...\mathbf{l_{q}}^{T},\sigma_{n1},...,\sigma_{np}],\label{eq:estparam}
\end{equation}
 where $\mathbf{\gamma_{i}}$  is the direction independent gain  of each antenna
 and $\phi_{i}$  is the  associated phase. $\sigma_{si}$  and $\mathbf{l_i}^{T}$
 are  the   estimated  fluxes  and   positions  of  the  model   sources,  while
 $\sigma_{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:
\begin{equation}
\bm{\hat{\theta}}=argmin\left(ln|\bm{R(\theta)}|+tr(\mathbf{R}^{-1}(\bm{\theta})\mathbf{\widehat{R}})\right)\label{eq:normeq}
\end{equation}
 where  $\bm{R(\theta)}$  is  the  model  covariance  matrix  as  a
 function of  the parameters $\bm{\theta}$,  and $\mathbf{\widehat{R}}$ 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 \citep{ottersten1998covariance}.  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).


\subsubsection{\label{sub:Model-generation}Model generation}
\begin{table*}[tbh]
\caption{Details of model sources (the A-team) used for All-sky self-calibration}
\label{tab:Details-of-model}
\center{%
\begin{tabular}{c c c c c}
\hline \hline
Src & $S_{60}(Jy)$ & 74 MHz Size  & \% Flux contribution  & Comment \\ 
\hline 
3C461 (CasA) & 19605.31 & 7' & $\sim$50 & Spectra from (1), Decay from (2)\\
3C405 (CygA) & 19201.65 & 4.5' & $\sim$42 & Spectra from (1)\\
3C144 (TauA) & 2417.31 & 8' & $\sim$11 & Spectra from (1)\\ 
3C274 (VirA) & 1656.46 & 14' & $\sim$8& Spectra from (2)\\
\hline 
\end{tabular}}
\tablebib{
 (1)\citet{baars1977absolute}; (2)\citet{helmboldt2008radio}
}
\tablefoot{  $S_{60}$   has  been  generated  using  the   spectral  indices  of
  \citet{helmboldt2008radio}, while  a secular decrease  in the flux of  CasA of
  0.8\% has been assumed. The  fractional flux contribution of each model source
  to  the  flux  in the  observed  snapshot  visibilities  assumes the  flux  of
  unmodeled sources  adds incoherently due  to a random sky  distribution, along
  with a source density of $\sim66 sr^{-1}$ at 178 MHz.}
\end{table*}

Forming the  normal equations of  Eq. \ref{eq:normeq} for  minimization requires
modeling    the    observed   visibilities    using    the    data   model    in
Eq. \ref{eq:datamodel} 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 \textbf{Northern} sky, viz.,
CasA, CygA, VirA, TauA and the  Sun, whose nominal positions are determined from
the revised third Cambridge (revised)  catalog (3CR).  The model parameters used
are summarized in Table  \ref{tab:Details-of-model}.  The efficacy of this model
for  parameter  estimation  is  elaborated  upon  using  real  data  in  Section
\ref{sec:calib-accuracy}.  Of  these model  sources, CasA 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
$\left(\mathbf{\Sigma_{s}}\right)$ 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  as well as 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  \ref{sub:Direction-dependent-gain} 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 \citep{ottersten1998covariance}.  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 such  effects, and  is
discussed in more detail in Section \ref{sub:System-noise-estimation}.
\begin{comment}
 Values  of  A-team   and  galactic  plane  contribution  is   based  on  actual
 contribution to visibilities in night time data (for CasA/CygA and Galaxy), and
 day  time data  for  VirA,TauA and  the  Sun.  The  total  power of  calibrated
 visibilities  =  sum(abs(acccalgal(:)).^2),   where  acccalgal  =  (gainsol'  *
 gainsol) .* acc; The Galactic  contribution = sum (abs(sigman(:)).^2), but also
 includes the autocorrs,  which is a sum  of Tsys and Tsky, and  which cannot be
 separated without external  measurements of Tsys. This explains  the high value
 of  galactic contribution, against  an expected  40% contribution.   The RAteam
 contribution  is  estimated  by   generating  model  visibilities  from  fluxes
 estimated by LSimaging,  ie, the apparent flux, and  so is elevation dependent.
 NOTE that  we dont  have observations in  which all  Ateam sources are  at high
 elevations,  in which  case we  multiply their  visibility contribution  by the
 approximate beam attenuation.   SNR was removed from the table,  as it does not
 provide very much extra information  as compared to the total flux contribution
 from a source.
\end{comment}

\subsubsection{\label{sub:Direction-independent-gain}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
\citep{salvini2014stef,salvini2014stefursi} algorithm for gain estimation.  This
algorithm  has  a  smaller  memory   footprint,  and  its  computing  scales  as
$\mathcal{O}(N^{2})$  for  N  parameters,  as against  the  $\mathcal{O}(N^{3})$
scaling of traditional solvers.  This significantly reduces the computing during
AARTFAAC calibration due  to the large number of  parameters ($\sim$$590$) being
estimated.  The StEFCal implementation led to a measured increase in calibration
throughput  of   a  factor  $\sim$30  \citep{salvini2014stef},   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.   These,  although  computationally
costly, are  essential to other solvers like  the Levenberg-Marquardt minimizer.
StEFCal algorithm has been shown to be unbiased with rapid convergence.

\textbf{Gain  Constraints:} Since  we simultaneously  solve for  apparent source
fluxes as  well as 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  instead use a  constraint on the  average gain amplitude  over all
antennas,  setting it to  unity.  This  reflects the  behaviour 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  $\chi^{2}$ error  surface towards  an identifiable
global minimum, whose convergence is identified  by a slow enough rate of change
of the estimated parameters.

\subsubsection{\label{sub:Direction-dependent-gain}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,  parameterized  in the  apparent  fluxes  and  positions of  the  model
sources. \textbf{The cycle begins with a one-time estimation of the model source
  positions via the WSF algorithm.}


\textbf{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)
\citep{viberg1991detection}.   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  $\lambda/D$ (D  being the
array size)  possible for sources with  a high SNR. 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.

%%   An  example  of  estimation  with  the  WSF  on  test
%% observations,   and   its   spatial    resolution   is   presented   in   Figure
%% \ref{fig:Estimated-flux-of}.



\begin{table*}[tbh]
\caption{Flux ratios of bright sources within snapshot from Obs. 3.}
\label{tab:Flux-ratios-of-1}
\center{%
\begin{tabular}{c c c c c c c c}
\hline \hline
Source & $S_{60}$ & $\sigma_{S_{60}}$ & Calibrated & $\sigma$  & Ratio & Comment\\
 & (Jy) &(Jy)  & flux (Jy) & (Jy) &  & \\ %\tabularnewline
\hline 
3C10    & 367.9 & 18.81 & 397.02 & 67 & 1.08 & VLSS $N_g$=4, $\alpha$=-0.7\\ 
3C452   & 152.2 & 10.99 & 184.52 & 28 & 1.21 & VLSS $N_g$=2, $\alpha$=-0.7\\
3C390.3 &  99   &  6.00 & 121.89 & 35 & 1.23 & VLSS $N_g$=3, $\alpha$=-0.7\\
3C380   & 156.2 &  2.99 & 156.2  & 38 & 1    & \citep{scaife2012broad}\\
3C295   & 134.5 &  8.31 & 95.23  & 37 & 0.71 & \citep{scaife2012broad}\\
%% 3C405 & 19202 &  & 26927 & 78 & 1.402 & Kassim et. al.\tabularnewline
%% \hline 
%% 3C461 & 16145 &  & 31758 & 89 & 1.967 & Kassim et. al.\tabularnewline
\hline 
\end{tabular}}
\tablefoot{The $S_{60}$  flux has been  derived from the best  available modeled
  spectra, or -0.7 otherwise. $N_g$ is  the number of Gaussian components in the
  source.  The  calibrated flux  ratio is derived  after imaging  the calibrated
  visibilities  and generating  light-curves of  peak flux  from sources  in the
  snapshots.}
\end{table*}


\textbf{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.

\textbf{Model    source   flux    estimation:}   As    explained    in   Section
\ref{sub:Model-generation},    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
\citep{wijnholds2008fundamental}  in an approach  called Least  Squares Imaging.
Here, the  array configuration is  assumed 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.

%% Figure  \ref{fig:Estimated-flux-of} shows  the extracted  flux of  the brightest
%% visible model  sources, and  clearly shows  the effect of  the element  beam and
%% ionospheric scintillation.

  The model  source flux estimator uses a  sparse weighting matrix
due to the modeling of mutual  coupling via a non-diagonal noise matrix (Section
\ref{sub:System-noise-estimation}), which prevents a bias in estimation.

\subsubsection{\label{sub:System-noise-estimation}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 \citep{agrawal1972mutual}.   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
$<-17dB$ from simulations \citep{wijnholds2011situ}.

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 poorer  sensitivity of the  PSF formed by  including these
baselines as a function of pointing  direction and away from the zenith compared
to  the   sensitivity  expected  after  pattern  multiplication,   as  found  by
\citet{ellingson2011sensitivity}.   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 few tens
of wavelengths separation  between antennas forming a baseline,  and so 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 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   \citet{wijnholds2010self}   by  an   additive
non-diagonal  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 non-diagonal 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, as well
as  the receiver  temperatures of  each antenna  and achieve  a PSF  with higher
sensitivity,    although    estimates    of    the   flux    can    be    biased
(\ref{sub:Direction-dependent-gain}).   Hence,  all  our estimators  ignore  the
short baselines (defined  as $<10\lambda$) to avoid the  interaction between the
excess flux at the shortest baselines and the calibration parameters.

%% In practice, it has been seen that the estimated values of the non-diagonal
%% covariance matrix are a substantial fraction of the total received
%% flux on the selected short baselines. Thus, the calibrated visibilities
%% have an inherent taper or 'hole' in the visibility plane which can
%% lead to ringing in the image domain, and needs to be tapered. 

%% The oversampled UV plane of the AARTFAAC allows implementing a spatial
%% filter with high roll-off and low sidelobes (put numbers?). In practice,
%% a filter with a $10-20\lambda$ cutoff has been found to be effective.
%% The filter taper is implemented in the form of a surface of revolution
%% of a gaussian to reduce ringing in the image domain.

\subsubsection{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 $\sim$$192$ kHz,
which  are  further filtered  to  give  the  correlator spectral  resolution  of
$\sim$$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  and typically
over a MHz. However, the LBA bandpass can have variations of up to 80\% from the
peak  resonance at  $\sim$$58$ MHz  \citep{vanhaarlem2013lofar}.  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  multi-frequency  synthesis mode  of
imaging.


\subsubsection{Flux density calibration}

The low resolution of the AARTFAAC allows  the use of a number of bright sources
for calibrating  visibility amplitudes.  Recently,  \citet{scaife2012broad} 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 \ref{tab:Flux-ratios-of-1} 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. Light curves  are then generated for the  various observations. A mean
integrated flux in  pixel counts is estimated by first  segmenting the data into
100 second 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.

%% while the $S_{60}$  for CasA and CygA has  been taken using the VLA  74 MHz flux
%% from \citet{kassim200774},  with a spectral  index (and the  frequency dependent
%% secular  flux decrease)  published by  Baars et.  al.  \citep{baars1977absolute}
%% spectral index  fits. The deviations from  published fluxes are  acceptable to a
%% transient detection system as long as systematics are absent. This is seen to be
%% true due to the limited range of fit errors on the calibrated sources.

\subsubsection{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 SNR 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, and so their Signal to Noise Ratio 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  $\sim48$  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 second.  Thus, in  our standard autonomous mode
of  operation, we  choose a  resolution  of 16kHz  and 1  seconds 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.


\subsection{Tracking calibration}

\begin{figure}[tbh]
\resizebox{\hsize}{!}{\includegraphics[width=1\columnwidth]{Figs/tracking_cal_schematic_new.png}}

\caption{\label{fig:trackcalSchematic}Schematic view of tracking calibration
for bounded latency precise calibration.}
\end{figure}

An  Alternating   Least  Squares  approach   is  already  known  to   have  high
computational  efficiency   over  closed  form  estimators,  as   well  as  fast
convergence  \citep{boonstra2003gain}.   Thus   this  approach  is  particularly
suitable for  our case due to  our latency and compute  budgets.  However, 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 timeslice.   Thus,  considering  the  known stability  of  the
instrument   (see  Sec.    \ref{sub:Stability}),  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 timeslice  to just one. This  is supplemented by  propagating solutions from
previous timeslices 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,  we choose to  use the  time axis due  to the possibility  of the
AARTFAAC subbands being chosen to sparsely sample the full bandwidth in order to
optimize  detection of  steep spectrum  transients.  Note  that computationally,
StEFcal  iterations are  relatively inexpensive  as compared  to e.g.,  WSF (see
Table  \ref{tab:Performance-results}).  This  implies that  their number  can be
increased  in order  to handle  data with  higher 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.

The feedback of solutions as initial estimates makes our approach susceptible to
errors in 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  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 Figure
\ref{fig:trackcalSchematic}.

The  StEFCal  algorithm already  utilizes  the  gain  estimate from  a  previous
iteration. Although  it is seen to converge  quite rapidly even for  a very poor
estimate, initialization with a gain close to the truth leads to even more rapid
convergence.   This can  be seen  in  Figure \ref{fig:Convergence-behaviour-of},
which  compares the  number of  iterations and  the residuals  for  an identical
timeslice, with no  initial gain estimate, and one  determined from the previous
timeslice.

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:
\begin{enumerate}
\item Carry  out a  calibration to convergence  over the  first time slice  in a
  chosen temporal  window at the specified spectral  resolution. This determines
  the baseline of parameter estimates.
\item Generate a mean over the short term history of solutions.
\item For every time and spectral bin within a window:
  \begin{enumerate}
    \item Precalibrate incoming ACM with the mean solution.
    \item Estimate positions of model sources via DOA estimation techniques.
    \item Carry out a single iteration of all minimization procedures within the
      algorithm, using  the previously generated baseline  parameters as initial
      estimates.
    \item  Compare  updated  parameter   estimates  for  temporal  and  spectral
      smoothness.  If a deviant solution  is found, reject; else update solution
      history.
  \end{enumerate}
\end{enumerate}

A  computational profile  of  the  various algorithmic  components  of both  the
convergent  and   the  tracking  calibration   schemes  is  depicted   in  Table
\ref{tab:Performance-results}.  The calibration solution memory pool is utilized
by  the  goodness  of  solution  block  to  decide  on  the  propagation  of  an
instantaneous solution.

\section{\label{sec:Performance-of-tracking}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.

\begin{table}[h]
\caption{Details of commissioning observations carried out with the AARTFAAC}
\label{tab:Details-of-commissioning}
\center{%
\begin{tabular}{ c c c c c c}
\hline  \hline
Obs. & Epoch & $\Delta t$ & $\Delta\nu$& Nchan & $\nu_{obs}$\\
     & UTC   & sec  & kHz & & MHz\\
\hline 
1.&2011-Sep-21,1200-1500 & 5 &3&64{*}5&60 \\
2.&2012-Jul-12,2222-2258 & 1 &3&64{*}5&60 \\
3.&2012-Jul-12,0140-0200 & 1 &3&64{*}5&60 \\
4.&2012-Jul-12,0340-0430 & 1 &3&64{*}5&60 \\
\hline 
\end{tabular}}
\end{table}

\textbf{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 sub-banded
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
\ref{tab:Details-of-commissioning}. 




\subsection{\label{sec:calib-accuracy}Calibration accuracy on a single timeslice}

\begin{figure}[tbh]
\resizebox{\hsize}{!}{\includegraphics[width=1\columnwidth]{Figs/major_minor_cycle_conv.eps}}

\caption{\label{fig:Convergence-behaviour-of}Typical   major  and   minor  cycle
  convergence rates on  a single timeslice.}
\end{figure}

In this section, we elaborate on the convergence behaviour of the calibration on
a typical data  segment, while commenting on the aspects  of weighting for model
source flux estimation  and gain constraints.  A timeslice 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 \citep{offringa2013lofar}, and does
not affect our conclusions.  Figure \ref{fig:Convergence-behaviour-of} 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 timeslice.  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.



\begin{figure*}[tbh]
\subfloat
%% [Model error as a function  of point source model complexity (abscissa)
%%   for a  typical dataset.  The model complexity  increases with the  addition of
%%   bright   sources  in   decreasing  order   of  expected   flux  contribution.]
  {\includegraphics[width=1\columnwidth]{Figs/model_efficacy.eps}}
%%  \hspace{10pt}
\subfloat
%% [Phase from the ratio of calibrated and model visibilities of a single 
%%   timeslice, plotted as a function of increasing calibrated visibility amplitude.
%%   An increase in phase error on higher amplitude visibilities can be seen to the
%%   right, but lies well within 1 radian. This is expected to be 0 for an ideal 
%%   model. The phase error is due to flux from unmodeled sources.]
  {\includegraphics[width=1\columnwidth]{Figs/selfcal_ph_behaviour_new.eps}}

\caption{\label{fig:The-model-amplitude}The  sky model  efficacy evaluated  on a
  single  timeslice.(Left) Effect  of increasing  sky model  complexity. (Right)
  Phase behaviour of residual visibilities after model subtraction.}
\end{figure*}

\textbf{Model efficacy:  } We utilize a  simple point source sky  model for self
calibration. The efficacy of this model can be evaluated via the modelling error
after calibration.   The parameters of the  brightest sources making  up the sky
model are given in Table.  \ref{tab:Details-of-model}, 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  Figure
\ref{fig:The-model-amplitude}  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 \ref{tab:Details-of-model}.  The
right panel of Figure \ref{fig:The-model-amplitude} 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, i.e., 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 SNR,  thus having a
higher calibration error.  Visibilities with high SNR 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.
%%  to  points towards The low phase residue with
%% 90\% of the visibilities having <0.5  rad phase error implies that the vectorial
%% contribution to a visibility from the residual sources of similar flux range add
%% up  incoherently,  as  they should  in  the  absence  of a  dominant  systematic
%% contribution from a bright source. 

%%%%%%% Subfigure arrangement: Doesn't work
%% \begin{figure}[tbh]
%%   \centering
%%         \begin{subfigure}[a]{0.5\textwidth}
%%           \resizebox{\hsize}{!}{\includegraphics[width=1\textwidth]{Figs/model_efficacy.eps}}
%%                 \caption{A gull}
%%                 \label{fig:gull}
%%         \end{subfigure} %%
%%         \begin{subfigure}[b]{0.5\textwidth}
%%           \resizebox{\hsize}{!}{\includegraphics[width=1\textwidth]{Figs/selfcal_phase_behaviour.png}}
%%                 \caption{B gull}
%%         \end{subfigure}
%% \end{figure}
\begin{figure}[tbh]
\resizebox{\hsize}{!}{\includegraphics[width=1\columnwidth]{Figs/SB000_ch30-35_5sec_3hr_9_convcalsol_bin_iter_majmin.eps}}

\caption{\label{fig:Major-and-minor}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}
\end{figure}


Figure   \ref{fig:Major-and-minor}   shows   the   longer  term   behaviour   of
instantaneously calibrating  every timeslice  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.

\begin{figure*}[t]
\includegraphics[width=1\textwidth]{Figs/uncal_cal_deconv_night.eps}
\caption{AARTFAAC zenith projection  all sky images at ~60  MHz, with 90 kHz
  and  1   second  integration.   (Left)  uncalibrated  image,   (Middle)  after
  calibration, and (Right) after subtraction of the brightest sources.}
\label{fig:uncal_cal_deconv_night}
\end{figure*}

Figure  \ref{fig:uncal_cal_deconv_night}   shows  an  all  sky   image  with  an
integration of  1 second  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  CasA, CygA and the  Galactic center at
the South west corner. The middle panel shows 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.

 
\subsection{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.
\begin{figure}[tbh]

\resizebox{\hsize}{!}{\includegraphics[width=1\columnwidth]{Figs/cmpcalsol.eps}}
\caption{\label{fig:The-error-timeseries}The error timeseries of tracking
calibration solutions as compared to those obtained after iterating
to convergence, for Obs. 3.}
\end{figure}

Figure  \ref{fig:The-error-timeseries}  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, probably due
to  the recombining 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   Section
\ref{sub:iono-effect-on-calib}.   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.
\subsection{\label{sub:Stability}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  non-convergence 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 behaviour of each parameter.
To  this  end,  Figure  \ref{fig:gain-Temporal-stability}  shows  the  estimated
complex gain  timeseries from a  few antennas over \textasciitilde{}3  hrs.  The
variance on the estimated phase of the instrument is \textasciitilde{}5 degrees,
while  the amplitude  has variations  of 2\%  on average.   In  conjunction with
Figure   \ref{fig:Estimated-flux-of},   Figure  \ref{fig:gain-Temporal-stability}
demonstrates  the resilience  of  the algorithm  to observational  perturbances,
while  showing  the  segregation  of  sources of  error  into  instrumental  and
observational.

\begin{figure*}[tbh]
\subfloat[Obs. 3]
{\includegraphics[width=1\columnwidth]{Figs/SB002_LBA_OUTER_SPREAD_1ch_8_convcalsol_bin_gains.eps}}
\subfloat[Obs. 4]
{\includegraphics[width=1\columnwidth]{Figs/combined_SB002_2_convcalsol_bin_gains.eps}}
\caption{\label{fig:gain-Temporal-stability}Temporal  stability  of the  complex
  gain of  a randomly chosen  set of 6  antennas from the 288  AARTFAAC antennas
  over  a  period  of \textasciitilde3  hrs.  The  outliers  are caused  due  to
  instances of significant RFI.}
\end{figure*}

\subsubsection{Light curve 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
Figure  \ref{fig:Light-curves-of}.  One  can  see that  the  dominant source  of
variation  in  the extracted  flux  is due  to  systematic  variations on  short
timescales, most probably caused by  the ionosphere.  This limits the estimation
accuracy on the extracted light-curve.

\begin{figure}[tbh]
\resizebox{\hsize}{!}{\includegraphics[width=1\columnwidth]{Figs/SB002_LBA_OUTER_SPREAD_1ch_6_convcal_7_el_fftimg_bin_save_mat_tseries.eps}}

\caption{\label{fig:Light-curves-of}Light curves of a few bright
field sources from images of Obs. 3 of the calibrated visibilities.}
\end{figure}
The position offset  of the extracted sources as  compared to cataloged position
shows  no  significant  bias,  leading  to  the conclusion  that  there  are  no
significant systematic calibration errors.


\subsection{\label{sub:iono-effect-on-calib}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 introducing
spatially  variant refraction and  propagation delays,  which can  require rapid
calibration cycles  to over sample  ionospheric phase coherency  timescales (see
Sec. \ref{sub:All-Sky-cal-iono}).

Recent results \citep{intema2009ionospheric} have indicated towards 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 \mbox{AARTFAAC}.  These are
expected to  be larger than the  typical scale size  of ionospheric fluctuations
when  projected  onto  ionospheric  heights,  resulting  in  \emph{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
\citep{vdTol2007selfcallofar}.

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,     as    well     as     position    wander.      Figure
\ref{fig:Estimated-flux-of} shows the effect  of the ionosphere on the estimated
flux of  the model  sources, as  well as the  position wander  of CasA  over the
course of a night.

\begin{figure*}[tbh]
\subfloat[Obs.2]
{\includegraphics[width=\columnwidth]{Figs/SB000_LBA_OUTER_2222-2258_1ch_1_convcalsol_bin_sigmas.eps}}
\subfloat[Obs. 3]
{\includegraphics[width=\columnwidth]{Figs/SB002_LBA_OUTER_SPREAD_1ch_8_convcalsol_bin_sigmas.eps}}\par
\subfloat[Obs. 4]
{\includegraphics[width=\columnwidth]{Figs/combined_SB002_2_convcalsol_bin_sigmas.eps}}
\subfloat[Obs. 3: CasA WSF estimated position variation from catalog ]
{\includegraphics[width=\columnwidth]{Figs/CasAsrcposoff.eps}}

\caption{\label{fig:Estimated-flux-of}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 CasA
  and CygA  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
  \textasciitilde{}1000 sec, as  compared to its catalog position.  Shown at two
  frequencies separated by \textasciitilde{}30MHz.}
\end{figure*}


The variation  of the position of  model sources from catalog  locations under a
typical  ionosphere is  demonstrated  by Figure  \ref{fig:Estimated-flux-of}(d).
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  behaviour 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 Figure \ref{fig:Estimated-flux-of}.

Such a direction dependent correction  is estimated and carried out only towards
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 light  curves  of field
sources 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 light  curves
extracted from calibrated images.


\subsection{\label{sub:Enhancing-the-transient}Enhancing the transient detection
sensitivity via difference imaging} 

\begin{figure}[h]
\resizebox{\hsize}{!}{\includegraphics[width=1\columnwidth]{Figs/diffimg.eps}}
\caption{\label{fig:Reduction-in-image}Reduction in image noise as a function of
  integration  time over  raw  and differenced  consecutive  image pairs.   Each
  member of the pair being differenced is generated by averaging raw images over
  increasing  cadences.}
\end{figure}

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
Figure  \ref{fig:Reduction-in-image}, which  shows  the decrease  in noise  with
integration  time for  a snapshot  image as  well as  a difference  image,  as a
function of integration time.  The difference image I$^N_{diff}$ with N units of
time integration is generated as:
\begin{equation}
I^N_{diff} = \sum\limits_{i=0}^{N-1} I_{i+N} - \sum\limits_{i=0}^{N-1} I_i \label{eq:diff_img}
\end{equation}
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
$\sqrt{N}$ enhancement on the integration of  N timeslices, and we see the noise
in a  difference image closely  following this curve  out to about  ~30 seconds,
beyond which residuals start getting  affected by earth rotation.  In principle,
short term image differences  (\textasciitilde{}secs) are enough to reduce noise
amplitude by 20-30\%.

\section{\label{sec:Computational-performance}Pipeline architecture and computational performance}
\begin{figure}[tbh]
\resizebox{\hsize}{!}{\includegraphics[width=1\columnwidth]{Figs/pipeline-design}}
\caption{\label{fig:pipeline}  High  level  schematic  of  the  calibration  and
  imaging pipeline design depicting the two stages of parallelism on the subband
  and channel axis.}
\end{figure}

This  section  describes the  implementation  of the  all  sky  monitor and  its
performance. The  upstream correlator implementation generates a  time series 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
timeslice,   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
Figure \ref{fig:pipeline}.   The second stage splits each  subband into channels
which  are  individually  flagged   and  calibrated,  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.

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   \citep{salvini2011pelican}   framework  for   data
distribution,  in   which  a  server  process  provides   high  throughput  data
distribution to  a collection  of worker threads  on a first-come,  first served
basis. This scheduling maintains 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 Figure.   \ref{fig:pipeline}) to implement  a server for distribution  of the
low bandwidth calibration solutions  to the corresponding channels. We implement
the  individual  estimators using  the  Eigen3  \citep{eigenweb} 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.


\subsection{Performance results}


The system processes a channel per CPU core i.e. per thread. We hence define the
performance   as   the   processing   time   of   a   single   channel.    Table
\ref{tab:Performance-results} 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, imaging) of a
subband with $N$  channels using $N$ cpus  to be $250$ ms ($\sigma  = 15$ms) for
nighttime data and  $1700$ ms ($\sigma = 100$ms) for daytime  data on a Intel
Core i7-2600 CPU at 3.40GHz, which is similar to our target ASM hardware.
\begin{table*}[tbh]
\caption{Functional profile of real-time calibration and imaging code}
\label{tab:Performance-results}
\center{%
\begin{tabular}{c c c c c}
\hline \hline
Function & Time measured & Time spent & Time measured & Time spent\tabularnewline
\hline 
Initialization & Night & 5\%  & Day & 3\%\tabularnewline
% \hline 
Model source pos est. & Night & 83\%  & Day & 92\%\tabularnewline
% \hline 
Major cycle & Night & 7\%  & Day & 4\%\tabularnewline
% \hline 
Minor cycle & Night & 5\%  & Day & 1\%\tabularnewline
\hline 
\end{tabular}}

\tablefoot{  Overview of  the  time  spent in  calibration  function blocks  for
  various  datasets  (daytime  and  nighttime).  We measured  the  overall  time
  (flagging, calibrating,  imaging) of an ACM  to be $250$ ms  ($\sigma = 15$ms)
  for nighttime  data and $1700$  ms ($\sigma =  100$ms) for daytime  data using
  current generation  hardware, and calibrating to convergence.  The day latency
  is seen to be  much worse due to the presence of a  very active Sun leading to
  significant modeling errors.}
\end{table*}

\section{\label{sec:Challenges-to-tracking}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 because of  a larger error  on the
estimation  parameters on  data  with poor  SNR  due to  the  limited number  of
iterations.   These effects  tend to  lower the  sensitivity of  each timeslice.
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.

\subsection{Radio Frequency Interference (RFI)}

RFI mimics  a coherent source across the  array, and is usually  very bright and
sporadic,  both in 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 \mbox{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 such  timeslices  which are
significantly different from  the maintained time history, and  the timeslice is
effectively rejected by the filters on the generated solutions.

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  \mbox{AARTFAAC}.  The  fast  recalibration  allows  us to  distinguish
individual  RFI corrupted  timeslices and  discard them  from  further analysis,
while weak  RFI which does  not affect calibration  will have to be  filtered by
higher layers like the TraP.

\subsection{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  \mbox{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
\mbox{AARTFAAC},  thus invalidating  the point  source assumption  of  the Solar
model. A  multi-component point  source Solar model  which varies with  time and
frequency would then be required for adequate modeling the solar contribution to
observed visibilities.  Such  a model can be extracted from  a snapshot via high
resolution  beamforming  towards the  Sun  and  the  calibration model  updated.
However,  our  latency  bounds  preclude  such an  operation  currently.   Thus,
timeslices 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
\mbox{AARTFAAC} images after tracking calibration is $\sim2200:1$, and so images
in the presence of the quiet Sun can be considered for transient detection.


\subsection{Creeping bias in tracking calibration}

Since both  the positions  of the  model sources as  well as  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  timeslice  as  an  initial  estimate,  which  can  be  biased  due  to
non-convergence  of  tracking calibration  on  data  with  poor SNR.   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  Figure \ref{fig:The-error-timeseries},  being subsequently  brought  back on
track via a convergent calibration.

%% \begin{figure*}[tbh]
%% \center{\includegraphics[width=0.7\paperwidth,height=7cm]{Figs/cmpcalsol}}\caption{\label{fig:The-error-timeseries-1}The error timeseries of tracking
%% calibration solutions as compared to those obtained after iterating
%% to convergence. }
%% \end{figure*}


%% \section{\label{sec:Discussion}Discussion}

%% The tracking  multisource calibration has  been found to  function adequately
%% within the  compute and latency  constraints of then \mbox{AARTFAAC}  ASM.  However,
%% several  other factors  can  influence its  efficiency.  These are  discussed
%% below.
%% \subsection{\label{sub:Confusion,-peeling-and}Confusion, peeling and the depth
%% of deconvolution }

%% The  current   algorithm  determines  model  parameters   for  the  brightest
%% unresolved sources  in the  sky to carry  out their  effective deconvolution.
%% These sources  account for  almost TODO fraction  of the received  power, and
%% hence their  deconvolution increases the  image dynamic range to  about TODO.
%% However,  the  effective  noise  in  the generated  instantaneous  images  is
%% $\sim20$ Jy,  accounting for the  contribution by the ionosphere.   Thus, the
%% limitation  on  achieved dynamic  range  is  constrained  by the  ionospheric
%% errors,   rather  than   by  sidelobe   confusion  of   lower   level  bright
%% sources. Further, does  the estimation of the model  source fluxes get biased
%% due to confusion  from background sources in the relatively  broad PSF of the
%% \mbox{AARTFAAC}?  What about  the  simultaneous  peeling we  carry  out, instead  of
%% peeling  individual sources  in decreasing  order  of flux?\textcolor{red}{{}
%% <TODO> }


%% \subsection{Effect of the presence of bright transients on model based calibration}

%% The model based Selfcal.  can be argued to have a high  dependency on the sky
%% model in order to estimate the calibration parameters. A bright enough source
%% not present in the skymodel could  potentially lead to convergence to a local
%% minimum, resulting in degraded sensitivity.  However, the relative brightness
%% of such a source needs to  be substantial to affect the calibration. We carry
%% out an  empirical test for  the flux limits  which affect calibration  to the
%% level  of masking  a bright  transient source.  However, the  brightest model
%% sources dominate the visibility amplitudes  by far, thus usually only RFI can
%% cause a model error. \textcolor{red}{<TODO?>}

\section{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 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 to be  immune to the
rapid  flux and  position  variations of  model  sources in  the  presence of  a
diffractive ionosphere.

Among  contemporary low  frequency  arrays, the  \mbox{AARTFAAC}  array is  well
suited  to  perform   as  a  real-time  monitor  for   bright  and  short  radio
transients.The \mbox{AARTFAAC} tracking calibration is optimized for the precise
2-D  nature  of  the ASM,  its  snapshot  imaging  mode  of operation,  and  its
sensitivity  and  field  of  view.  Other  established  calibration  algorithms,
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
\mbox{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  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 $\sim$ 0.25 seconds on a single time and frequency timeslice 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 \mbox{AARTFAAC} all-sky monitor.

\begin {acknowledgements}
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  \mbox{AARTFAAC}.   We  thank  \textbf{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.
\end{acknowledgements}
\bibliographystyle{aa}
\bibliography{ref}

\end{document}