© ESO, 2013

1. Introduction

A major epoch in the history of the Universe yet to be understood in detail is its dark ages and the epoch of reionization (EoR). Observational evidence for this era can be gathered with high probability by studying the fluctuations of the redshifted neutral hydrogen at redshifts corresponding to 6 < z < 12. Therefore, there are numerous experiments becoming operational and already collecting data, especially in the frequency range from 115 MHz to 240 MHz to reach this goal.

At the forefront of such experiments is the LOw Frequency ARray (LOFAR; van Haarlem et al., in prep.). Similar EoR experiments using other radio telescopes are already underway. For instance, Paciga et al. (2011) provide a new lower bound for the statistical detection threshold of HI fluctuations using the Giant Metrewave Radio Telescope (GMRT). While Murchison Widefield Array (MWA) and Precision Array to Probe the EoR (PAPER) are not in full hardware deployment yet, there are still relevant results being produced. In Ord et al. (2010) and Williams et al. (2012), initial widefield images of the southern sky using 32 MWA Tiles are presented. In Jacobs et al. (2011), full sky images and source catalogs using PAPER are presented.

In preparing for the LOFAR EoR project we have conducted several pilot experiments with the Low Frequency Frontends on the WSRT in a relevant frequency range: 138–157 MHz. The results of these observations, and a discussion of their limitations, have been described by Bernardi et al. (2009, 2010). For LOFAR in its commissioning phase we have adopted a multi-faceted observing strategy, building on the experience gained from the WSRT data. The rationale behind this is described in more detail in de Bruyn et al. (in prep.). A brief summary follows:

Using LOFAR in its commissioning phase we have observed and processed three very diverse windows. One window contains a very bright compact source, 3C196, which allows exquisite absolute calibration as well as a study of the systematics at very high spectral and image dynamic range. At a declination of only 48 degrees it will allow a study of elevation dependent effects. The many bright compact field sources around 3C196 also allow a study of ionospheric non-isoplanaticity. The results of these observations, with emphasis on all those topics, will be described in a separate paper Labropoulos et al. (in prep.). The second window was chosen to ascertain possible damaging effects of faint signals due to instrumental leakage of bright polarized Galactic foreground signals. These results, on the Elais N1 window, will be described by Jelic et al. (in prep.) in the second paper in this series.

The third window, without bright sources and with relatively faint diffuse linearly polarized emission from our Galaxy, is accessible throughout the year. This window is centered on the north celestial pole (NCP). It has been used to experiment with various calibration approaches as well as conduct a thorough analysis of the noise levels attainable with the current, still incomplete, LOFAR array. These results are the subject of the first paper in the series.

The three windows described above have thus far been observed using a single digital beam with about 48 MHz bandwidth. Between these three windows we expect to address most of the issues that will affect much longer observations with LOFAR, which should go one order of magnitude deeper in noise level. The analysis of the results obtained in these three windows, as they pertain to EoR signal levels, will be discussed in more detail in subsequent publications.

In this paper, we present results of LOFAR observations pointed at the NCP in the frequency range 115 MHz to 163 MHz. The NCP was previously observed using the Westerbork Synthesis Radio Telescope (WSRT) in a similar frequency range, albeit with limited integration time and resolution. As reported by Bernardi et al. (2010), the WSRT observations are mainly limited by broadband (and low level) radio frequency interference (RFI), and classical confusion (due to having limited <3 km longest baseline) that prevents reaching the theoretical noise level.

LOFAR provides significant challenges as well as advantages over conventional low frequency radio telescopes in terms of calibration. Unlike the WSRT, antennas on the ground are much less susceptible to broadband RFI. On the other hand, calibration of a LOFAR observation is challenging due to many reasons including spatially and temporally varying beam shapes with wide fields of view as well as mild to severe ionospheric distortions. Therefore, it is paramount that we test and demonstrate the feasibility of LOFAR for EoR observations before starting the long, dedicated observing campaign.

The results reported in this paper are based on integrations of the NCP which consisted of 3 nights with 6 h each nighttime observing. We provide details of the calibration and imaging that lead us to almost reach the expected theoretical noise level (within a factor of 1.4). We also provide details of current limitations and what we expect with the current, still incomplete LOFAR. Based on this result, we see no show stoppers for the launch of the dedicated LOFAR EoR observations on this window which will last several hundred hours.

This paper is organized as follows: in Sect. 2, we give an extensive overview of the observational setup. Section 3 describes the data processing pipeline. In Sect. 4, we present initial results with deep very widefield images, new sources, and the noise behaviour. We give an analytical explanation for the noise behavior in Sect. 5 that considers the excess noise due to sources far away from the phase center. Finally, we draw our conclusions in Sect. 6.

Notation (mostly in Sect. 5): we use bold lowercase letters for vectors and bold uppercase letters for matrices. The matrix Frobenius norm is given by ∥ . ∥ . The matrix transpose, Hermitian transpose and pseudoinverse are given by (.)^T, (.)^H and (.)^†, respectively. The trace of a matrix is given by trace(.). The real part of a complex number is denoted by Re(.). The statistical expectation operator is denoted by E {.}.

2. Observational setup

In this section we provide details of the LOFAR stations used in the NCP observations. We also provide the motivation behind observing the NCP.

2.1. LOFAR stations

We give a brief overview of LOFAR hardware and a complete overview can be found in van Haarlem et al. (in prep.) and de Bruyn et al. (in prep.). Each LOFAR HBA station consists of multiple elements (dipoles) with dual, linear polarized receivers. For a core station (CS), there are 384 dipoles, arranged in 24 tiles that have dipoles on a 4 × 4 grid. For a remote station (RS), there are 768 dipoles, arranged in 48 tiles. The signals of each dipole in a tile are coherently added (or beamformed) to form a narrow field of view (FOV) along a given direction in the sky. The effective beam shape is the product of the array (beamformer) gain with the dipole beam shape. It should be noted that the dipole beam shape is strongly polarized, and along some directions in the sky the polarization could be as much as 20% of the total intensity.

The nominal LOFAR FOV at around 150 MHz at the NCP (from null to null) is about 11 degrees in diameter for a CS and about 8 degrees in diameter for a RS. Therefore, the effective FOV is about 10 degrees in diameter. There is also a complicated low level sidelobe pattern surrounding the FOV, and the sidelobes change with time and frequency, as the beamformer weights change. In order to minimize the cumulative effect of the sidelobes, each LOFAR station is given a different rotation in its dipole layout (but keeping the dipoles parallel). Due to this reason, each LOFAR station will have a unique beam shape, that varies in time, frequency as well as according to the direction being pointed at. For a widefield image, this naturally leads to beam variations that depend on time, frequency as well as the direction in the sky.

2.2. The NCP and its surroundings

It is of significant importance that the NCP FOV lies on a relatively cold (i.e., having low sky temperature) spot in the Galactic halo in order to minimize the effects of Galactic foregrounds. In Fig. 1, we show the NCP window (or FOV), which is overlaid on a full sky image observed at 50 MHz. The full sky image was made using an array of 16 LOFAR lowband dipoles, with the longest baseline of 450 m.

The Galactic plane lies along an arc joining Cassiopeia A (CasA) and Cygnus A (CygA) in Fig. 1 but it is resolved. CasA is 32 degrees away from the pole while CygA is about 40 degrees away. The closest 3C source to the pole is 3C61.1 which is 4 degrees away, with a total flux about 35 Jy (peak about 7 Jy) at 150 MHz. It is fully resolved and an image made from a previous LOFAR observation is shown in Fig. 2. There are several other bright 3C sources in the vicinity including 3C 390.3 (11 degrees away).

Due to the station beams, many of the strong sources outside the main beam are suppressed significantly and the brightest (apparent) source in the FOV is NVSS J011732+892848, about 30′ away from the pole. This source is unresolved in our observations and has a peak flux of 5.4 Jy at 352 MHz (Rengelink et al. 1997) and a peak flux of 5.3 Jy at 74 MHz (Cohen et al. 2007) as reported in the WENSS and VLSS surveys, respectively. Therefore, we assume this source to have a flat spectrum within the observing band as in Fig. 3, which makes it a suitable candidate for absolute flux calibration and noise estimation. Because it is very close to the pole, the nominal station beam gain is equal to the beam gain along the direction of the NCP (which is unity) and constant in time. However, the element (dipole) gain varies and this is taken into account during calibration.

2.3. Motivation for observing the NCP

The NCP is one of several observational windows for LOFAR EoR observations (de Bruyn et al., in prep.). The reasons behind choosing the NCP are numerous although this does not imply it to be the optimal choice. We list some of the positive and some of the potentially negative aspects of this choice:

• +

  The geographical location of LOFAR makes the NCP windowobservable at night time, throughout the year, at high elevation(53 degrees).

• +

  Due to minimum projection effects of the uv tracks, we get almost circular uv coverage and therefore, we get an almost circular point spread function (PSF).

• ±

  The strongest source in the FOV, 3C61.1 is attenuated by about 70% and the strongest (apparent) source is about 5 Jy in peak flux. Not having a strong source in the FOV has both advantages and disadvantages. First, not having a strong source means a not so high signal-to-noise ratio in calibration. However, there are less artefacts resulting from deconvolution residuals of strong sources.

• ±

  Geostationary RFI that manages to escape flagging routines, which work on high-noise samples, may end up near the NCP. Nonetheless, this provides a sensitive diagnostic of the presence of any faint, stationary, undetected RFI.

• −

  Finally, the NCP is located at a Galactic latitude of only 38 degrees. The overall system noise is therefore higher than the coldest regions near the north Galactic pole.

2.4. Observational parameters

We used 40 core stations and 7 remote stations in our observations. The shortest baseline is 60 m and the longest baseline is about 30 km. The observing frequency range is from 115 MHz to 163 MHz. There are about 240 subbands in each observation, within this observing frequency range. Each subband has 256 channels, covering a bandwidth of 195 kHz. The monochromatic uv coverage for a typical 6 h NCP observation is shown in Fig. 4.

With this uv coverage, we get a resolution of about 12′′ at 150 MHz. For comparison, in previous NCP observations using the WSRT (Bernardi et al. 2010), the longest baseline used was only 2.7 km yielding an angular resolution of only 120′′. The correlator integration time is set at 2 s. We used data taken on 3 different nights for the results presented in this paper. In Table 1, we summarize the observational parameters for these 3 different nights.

3. Data reduction

In this section, we describe the major steps taken to calibrate the NCP observations. Apart from the initial processing, the data was completely processed in a CPU/GPU¹ cluster dedicated for LOFAR EoR computing needs. The processing of these observations also enabled us to fine tune the software used in various processing steps.

3.1. Initial processing

The LOFAR correlator (Romein et al. 2010) outputs data at a very fine resolution (2 s and 0.78 kHz), mainly to facilitate RFI mitigation. However, the data volume makes it cumbersome for further processing. Therefore, the data of each subband (having 256 channels at 2 s integration) is flagged using the aoflagger (Offringa et al. 2010, 2012) and averaged to 15 channels (after removing the 8 channels at both subband edges, mainly to remove edge effects from the polyphase filter). This significantly reduces the size of the data that has to be processed in the following stages. However, with improvements in software, we intend to process data at a finer frequency resolution, once regular EoR observing has begun.

3.2. Sky model

Because there is no single bright source in the FOV, the sky model used for initial calibration of the NCP data contains about three hundred discrete sources, spread across the FOV of about 10 × 10 square degrees. The most complex source in this region is 3C61.1 as shown in Fig. 2. In order to efficiently model this source, we use a model including shapelets (Yatawatta 2011) and point sources. The rest of the sky model is modeled as a set of discrete sources, having multiple point source components. The brightest source (NVSS J011732+892848) is modeled as a single point source with a flat spectrum and a peak flux of 5.3 Jy. We like to emphasize here that we diverge from the traditional “clean component” based sky model construction in order to minimize the number of components used without the loss of accuracy (Yatawatta 2010). In order to automate this process, we have developed custom software (buildsky) that creates a sky model with the minimum number of source components required. The principle behind buildsky is to select the simplest model for a given source by choosing the correct number of degrees of freedom (Yatawatta 2011). While a point source has only one degree of freedom (for its shape), a double source has two and so on. There are additional degrees of freedom due to its position and flux. We use information theoretic criteria as given in Yatawatta (2011) to select the optimum number of degrees of freedom for any given source. All the sources in the sky model are unpolarized. The sky model was updated using two calibration and imaging cycles.

3.3. Calibration

The aim of calibration of LOFAR EoR observations is twofold: (i) correction for instrumental and ionospheric errors in the data; (ii) removal of strong foreground sources from the data such that specialized foreground removal algorithms (e.g. Harker et al. 2009) can be applied. The basic description of the LOFAR EoR data model used in calibration is given by Labropoulos (2010). We use an enhanced version of the LOFAR calibration pipeline (Pandey et al. 2009) for the EoR data calibration.

3.3.1. Data correction

Major steps in our calibration pipeline are as follows:

• 1.

  We first calibrate for clock errors as well as small time scaleionospheric errors along the center of the FOV. This is the so calleduv plane or direction independent calibration (Labropoulos 2010) and is performed using the Black Board Selfcal (BBS) package (Pandey et al. 2009). At this stage, each subband has 15 channels at 2 s integration time. We determine the calibration solutions for every 10 s, and one solution per subband. Since we do not have a dominant source at the center of the FOV, the solutions thus obtained correspond to small time scale ionospheric phase fluctuations common to the full FOV plus the clock errors.

• 2.

  Once uv plane calibration is done, the data is corrected for these errors. The data is also corrected for the element beam gain along the center of the FOV. As discussed previously, the dipole beam of LOFAR is strongly polarized and we use an element beam model based on numerical simulations. For an area about 10 degrees in diameter in the sky, the variation of the dipole beam shape is assumed to be small and correction for the center of the FOV is considered accurate enough for the full FOV.

• 3.

  The corrected data is averaged to 183 kHz (one channel per subband) and 10 s integration time. The data is also flagged by clipping any spikes present in the data after correction.

3.3.2. Source subtraction

LOFAR has a very wide FOV and along each direction, the errors present in the data are different due to varying beam shape and ionospheric effects (Koopmans 2010). Therefore, source subtraction is not a simple deconvolution problem for LOFAR observations. Even for a simple deconvolution, it is better to subtract the sources directly from the visibilities (Yatawatta 2010). In the case of LOFAR, this subtraction has to be done with the appropriate gain corrections along each direction.

In order to efficiently and accurately solve the multi-source calibration problem, we have developed algorithms and software based on Expectation Maximization (Yatawatta et al. 2009, 2012; Kazemi et al. 2011b). We have implemented these algorithms (SAGECal) with accelerated processing using graphics processing units (GPUs). In the NCP window, there are about 500 bright discrete sources (note that we subtract more sources than what we use for the uv plane calibration, for which we only use about 300) that are subtracted from calibrated visibilities with the correct directional gains. We have “clustered” (Kazemi et al. 2011a) these sources into about 150 different directions. Thus we effectively determine the errors along 150 directions during the source subtraction. An example of clustering is shown in Fig. 5. In the left panel of Fig. 5, we show two sources that are about 5′ apart and apparently having identical error patterns. Therefore, instead of calibrating along each source individually, we can cluster them into one complex source and determine the common errors. The right panel in Fig. 5 shows the result after subtracting the cluster and restoring the model.

There are still some errors remaining in the right panel of Fig. 5 mainly because of errors in the sky model and due to the effect of surrounding sources (that were not included in the sky model) and also due to short time scale ionospheric errors. It should also be mentioned that while most of the sources subtracted lie within the FOV, we have also subtracted strong sources far away from the NCP as shown in Fig. 1, for example CasA.

3.4. Imaging

We make images at different stages during calibration. All images are made using CASA². In order to update the sky model, we make images of the calibrated and source subtracted data. We keep the highest available resolution in order to create accurate source models. For the results presented in this paper, we have a resolution of about 12′′ and we choose a pixel size of 4′′ with uniform weighting. Even though the nominal FOV is about 10 degrees, we make images that have an FOV of about 13 degrees, to detect sources at the edge of the beam. We restore the subtracted sources onto these images, after convolving with the nominal (Gaussian) PSF. Afterwards, we use Duchamp (Whiting 2012) to select areas with positive flux and update the sky model using buildsky.

More relevant for EoR signal detection are images made at low resolution using the short baselines of LOFAR. Therefore, we also make images using baselines less than 1200 wavelengths at 35′′ pixel size. The image size is chosen to be about 65 degrees so that we can see any contributions from the Galactic plane, which is about 30 degrees away from the NCP.

4. Results

In this section, we present results mainly to highlight important stages in the calibration and finally, to present the noise limits that we have reached using LOFAR. The results are based on all three datasets given in Table 1, unless stated otherwise, and continuum images are made using inverse variance weighted averaging of the 240 subband images.

4.1. The performance of SAGECal in directional calibration

The effects due to beam shapes and the ionosphere are major causes of errors in LOFAR images. Therefore, directional calibration is essential. We present a few images to highlight the performance of SAGECal.

In Fig. 6, we present the area around the brightest source NVSS J011732+892848. The image on the left is before multi-directional calibration and source subtraction, with only a deep CLEAN based deconvolution. The image on the right in Fig. 6 is after running SAGECal and after restoring the sky model onto the residual image. It is clear that the errors due to beam variations and ionospheric variations have largely been eliminated in the right panel of Fig. 6. Some errors are still present, mainly caused by inaccurate source models used in multi-directional calibration. We emphasize that longer baselines are needed to construct accurate models for such complex sources.

In Fig. 7, we give another example for the effect of the time interval chosen in SAGECal for multi-directional calibration. The image on the left is without any multi-directional calibration and significant errors due to beam shape and ionosphere are visible. The image in the middle of Fig. 7 shows the image obtained after running SAGECal with directional calibration performed every 20 min. The beam variations, which are slower, are completely eliminated by this procedure. However, the ionospheric variations, that could have time scales much less than 20 min are still present.

The right panel of Fig. 7 shows the result after running SAGECal with a “hybrid” solution scheme. In this case, we solve for bright sources once every 5 min and for fainter sources once every 20 min. Most of the ionospheric errors present in the middle panel have been removed in this figure. There are still errors due to inaccurate source models (the source was assumed to be a perfect point source, but this is not accurate enough) and also due to ionospheric phase variations with a time scale smaller than 5 min.

The time interval of 20 min chosen for obtaining the solutions gives us 120 time samples (each sample is of 10 s duration). For each time sample, we have 990 baselines with 45 stations. Therefore, we have about 8 × 120 × 990 = 1 million real constraints to obtain a solution. The number of real parameters in a solution is 45 × 150 × 8 = 54   000 for 150 directions in the sky. Therefore, the ratio between the number of constraints and the number of parameters is about 18 which is more than sufficient to obtain a reliable solution. This can be further improved by using data points at different frequencies as proposed by Bregman (2012).

4.2. Widefield images

We present widefield images obtained for the full dataset given in Table 1. First, in Fig. 8, we present the image obtained after calibration as described in Sect. 3.3.1, but before running SAGECal. Therefore, no source subtraction is performed and only traditional CLEAN based deconvolution has been applied. The circle indicates an area of diameter 10 degrees. The peak flux of this image is about 5.3 Jy and the noise level is about 400 μJy/PSF. The complex source 3C61.1 is at the bottom left hand corner. The striking features in this image are the artefacts surrounding almost every source. As described previously, there are three major reasons for these artefacts: (i) varying LOFAR beam shapes which are different for each station; (ii) ionospheric phase errors; and (iii) classical deconvolution errors due to having partially resolved sources. For instance for the case of 3C61.1, all three of the aforementioned causes create artefacts, which are clearly visible close to the bottom left hand corner.

The only way to improve the image in Fig. 8 is multi-directional calibration as described in Sect. 3.3.2. We have shown the image obtained after running SAGECal in Fig. 9. The circle indicates an area of diameter 10 degrees. Comparison of Figs. 8 and 9 shows that most significant artefacts in Fig. 8 have been eliminated in Fig. 9. The prominent artefacts that still remain are due to the fact that CS and RS beam shapes have different FOVs and also due to frequency smearing. We now reach a noise level of about 100 μJy/PSF at the outskirts of Fig. 9, while the peak value in the image is about 5.3 Jy. This corresponds to a formal dynamic range of 50 000:1.

In Fig. 10, we give the image made only with the short baselines (<1200 wavelengths) using the same data of Fig. 9. In this image, the circle shows an area with 10 degrees in diameter and the density of the sources close to this circle is clearly less than in other areas of the image. We also see a significant number of sources away from the FOV that are seen through sidelobes of the beam. Most of these sources have not been included in our multi-directional calibration and hence, they have significant artefacts. The noise level in this image is about 300 μJy/PSF with the peak flux of about 5.3 Jy. Note, however, that the noise level is a strong function of the distance from the field centre.

4.3. New sources

Since we reach a noise limit of about 100 μJy/PSF, we detect a large number of sources that have not been detected in previous observations, even at higher frequencies. In Fig. 11, we present a small area (0.6 × 1.0 degrees) of Fig. 9 to compare to an image from WENSS.

We present an area close to the NCP in Fig. 12. The left panel in Fig. 12 shows an image made with all baselines which gives a PSF of 12′′. The right panel in Fig. 12 shows an image made using only the core baselines and has a PSF of about 150′′. An important lesson that can be drawn from Fig. 12 is the total absence of artefacts close to the pole. If there would be any residual geostationary RFI, their effects would accumulate near the pole. However, we see no unexplained artefacts.

4.4. Effects of bright sources at large angular distances

As shown in Fig. 1, there are a few bright sources in the neighborhood of the NCP. We have already mentioned 3C61.1 which is still well inside the FOV. The other source that has a significant effect is CasA, which is about 30 degrees away from the NCP. In Fig. 13, we present the images around CasA, made with only the core station baselines. The images with baselines using core stations only are more affected by CasA than images that include remote stations. There are at least four reasons that contribute to this. First, core stations have wider station beams (compared to a remote stations) and therefore, CasA is less attenuated on core-core baselines. Secondly, the core stations have more short baselines than remote stations, hence see more flux from CasA, which is heavily resolved at baselines longer than 1000λ. Thirdly, time and frequency smearing lead to a significant attenuation of the visibilities of distant sources. Fourthly, ionospheric effects such as non-isoplanaticity rapidly increase with the length of the interferometer baseline.

In Fig. 13 left panel, we show an image where we have run SAGECal while ignoring the effect of CasA. In other words, we did not include CasA in our sky model and neither did we solve along the direction of CasA. CasA is at the bottom of this image and is heavily distorted due to beam and ionospheric errors as well as time and frequency smearing. The ringlike structures centered at CasA and spreading throughout the image is clearly visible. However, there is an area shaped like a cone that is directed towards the NCP where there are no ripples. This is due to the fact that multi-directional calibration that ignores CasA (mainly close to the NCP) has absorbed the effect of CasA directed towards the NCP. Similar effects can be seen in sequential source subtraction schemes such as “peeling”.

In the right panel of Fig. 13, we have included CasA in our multi-directional calibration using SAGECal. Most of the ripples in the left panel are eliminated in the right panel of Fig. 13. Furthermore, the noise level is reduced by a few percent after including CasA in the calibration. There are still some errors close to the location of CasA. This is mainly due to errors in the CasA source model used in the subtraction and we expect to get better results with an updated source model (Yatawatta et al., in prep.).

Because CasA is a strong source, we clearly see its effect as we have just described. Conversely, even faint sources would have a similar effect, albeit at a low magnitude. We perform a statistical analysis of this effect due to faint sources in Sect. 5.

4.5. Noise

One important question that needs to be answered is whether we have reached the theoretical noise limits in our images and if not, provide plausible explanations for the difference. The theoretical noise can be calculated as follows (Nijboer et al. 2009): ${\mathrm{noise}}_{\mathrm{subband}}\mathrm{=}\frac{{\mathit{w}}_{\mathrm{imaging}}}{\sqrt{\mathrm{4}{\mathrm{\Delta }}_{\mathit{f}}{\mathrm{\Delta }}_{\mathit{t}}\left(\frac{{\mathit{N}}_{\mathrm{c}}\mathrm{(}{\mathit{N}}_{\mathrm{c}}\mathrm{-}\mathrm{1}\mathrm{)}}{\mathrm{2}{\mathit{S}}_{\mathrm{c}}^{\mathrm{2}}}\mathrm{+}\frac{{\mathit{N}}_{\mathrm{c}}{\mathit{N}}_{\mathrm{r}}}{{\mathit{S}}_{\mathrm{c}}{\mathit{S}}_{\mathrm{r}}}\mathrm{+}\frac{{\mathit{N}}_{\mathrm{r}}\mathrm{(}{\mathit{N}}_{\mathrm{r}}\mathrm{-}\mathrm{1}\mathrm{)}}{\mathrm{2}{\mathit{S}}_{\mathrm{r}}^{\mathrm{2}}}\right)}}\mathrm{·}$(1)In (1), the noise per subband (bandwidth 183 kHz) is given by noise_subband while the system equivalent flux densities for core and remote stations are given by S_c and S_r, respectively. We assume S_c = 3360 Jy and S_r = 1680 Jy, taking the proximity to the Galactic plane into account. The total bandwidth and integration time are considered as Δ_f = 183 kHz, and Δ_t = 6 × 3600 s. The total number of core and remote stations are given by N_c and N_r. The scale factor w_imaging is used to take imaging weights into account and we take this equal to 2 for uniform weighting. For the observation with highest sensitivity from Table 1 (L24560), with N_c = 38, N_r = 7, we get noise_subband = 1.46 mJy. Moreover, for the full observation, with 229 subbands, and taking the variation of S_c and S_r with frequency into account, we get a theoretical noise value of about 93 μJy. As seen in Table 1, we are only a factor of 1.4 from reaching the theoretical noise limit. An important question is how to estimate the noise in the images, in other words, which area of the image to use to calculate the noise. In Fig. 14, we show the variation of the noise across the 13 × 13 degrees image shown in Fig. 9. As expected, the noise close to the NCP is higher and is about 180 μJy due to additional contamination by unsubtracted compact sources as well as by diffuse foregrounds. However, with the subtraction of many more sources as well as foregrounds, we expect to reach the noise that we see at the edge of the image in Fig. 14. Therefore, we study the behavior of the noise seen at the edge of this image as this is what we intend to reach. We emphasize that the achieved continuum noise levels in low resolution images, i.e. made using only the core stations, will be limited by classical confusion noise in the inner parts of the images. Bernardi et al. (2009) and Pizzo et al. (2011) estimated this to be about 0.6 mJy for the WSRT at a frequency of 140 MHz. Because the LOFAR core has a similar extent as the WSRT we expect this to be the asymptotic limit of the continuum noise level in core-only LOFAR observations.

To study the noise behavior in more detail, we provide several figures where we give the noise (or the standard deviation) of a small rectangular area in the images (about 4 degrees away from the NCP) at different frequencies and using different observations listed in Table 1. In Fig. 15, we present noise levels per subband (standard deviation), determined with images made using the three nights with imaging parameters as in Fig. 9. Additionally, we also show the noise estimated at the NCP in Fig. 15. The main conclusion that can be drawn from Fig. 15 is the variability of the noise (or the sensitivity) of LOFAR from night to night. This is due to some stations not working or not working at full sensitivity at different nights. We expect this to be much more stable before the commencement of dedicated EoR observations.

The best night in Fig. 15 is for the observation number L24560 of Table 1 which has a noise level of about 2 mJy per subband at the high frequency end. There is a steep rise in the noise level at frequencies below 130 MHz. This we attribute to the rising contribution of Galactic background noise, increased sidelobe noise from an increasing number of faint background sources (due to both a wider primary beam and steeper source spectra), emission from the Galactic plane and other very bright sources like CygA (see Sect. 5). There are also some spikes and dips in the noise curves where RFI removal has flagged significant amounts (more than 30%) of data. We have discarded such images from further analysis.

In Fig. 16, we give the image difference noise estimates for the best night (L24560) of Table 1. The image difference noise is calculated by (i) subtracting two images adjacent in frequency; (ii) estimating the standard deviation of the subtracted image; (iii) multiplying the standard deviation by $\mathrm{1}\mathit{/}\sqrt{\mathrm{2}}$. We have also plotted the absolute noise level in the same figure. It is clear from Fig. 16 that the absolute and differential noise curves are almost identical. This suggests that our absolute noise estimate does not contain any frequency independent systematic effects, except systematic effects that are uncorrelated between adjacent images. Future processing would take the difference at narrower frequencies to verify this.

We have also plotted in Fig. 16 image difference noise curves where instead of taking the difference of adjacent subbands, we have averaged B images (or subbands) together and calculated the B difference between images. We have done this for B = 2,4,8 in Fig. 16. As B increases, the differential noise should be more affected by systematic effects due to the large bandwidth of averaging. This, in turn should be reflected by the differential noise not decreasing as $\mathrm{1}\mathit{/}\sqrt{\mathit{B}}$. However, we do see a decrease by $\mathrm{1}\mathit{/}\sqrt{\mathit{B}}$ in the plots in Fig. 16, even when B = 8.

We also show the noise comparison between the images made with all baselines and with the images made with core only (<1200 wavelengths) baselines in Fig. 17. The corresponding continuum image is given in Fig. 10 where the imaging parameters are also given. We have estimated the noise of image with core only baselines at a location about 5 degrees away from the center of Fig. 10.

In Fig. 18, we give the image difference noise plots for images made only with core baselines. We have also plotted the absolute noise curve which agrees well with the image difference noise. Similar to Fig. 16, we have also averaged adjacent subbands of groups 2, 4 and 8 and found the differential noise between averaged images.

In Fig. 19, we show the noise plots for all 4 Stokes images (I,Q,U,V) for one night (L24560) of observed data. The noise level of Stokes I is significantly higher which we attribute to unsubtracted sources. We give a detailed analysis of this in Sect. 5.

4.6. Linear polarization

Due to the fact that we correct for the element (dipole) beam polarization along the direction of the NCP during calibration, we do not see substantial instrumental polarization in our Stokes Q, U and V images. Note also that we expect most of the sources to be intrinsically unpolarized in this frequency range. Even though we only correct for the element beam gain along the direction of the NCP, the relative variation within the full FOV of the dipole beam shape is very little (therefore this correction is satisfactory within the full FOV). However, we expect to see an increase in instrumental polarization at the edge of the FOV but the station beam attenuation makes this instrumental polarization less obvious. During the multi-directional calibration phase using SAGECal, we also solve for a full Jones matrix and therefore, what remains of this instrumental polarization is correctly subtracted.

In order to detect any diffuse Galactic foregrounds that might appear in polarized images, we have performed rotation measure (RM) synthesis (Brentjens & de Bruyn 2005) using data from one of the nights listed in Table 1 (L24560). We show the total polarized intensity image for RM = 0 rad m^-2 in Fig. 20. The noise in this image is about 110 μJy/PSF and apart from many weakly instrumentally polarized discrete sources, we detect only very faint diffuse structure, at any value of RM. However, we have seen that subtraction of compact sources using SAGECal would suppress the unmodeled diffuse foregrounds. To alleviate this from happening in future processing, we would ignore the contribution from short baselines while running SAGECal for calibration along multiple directions.

5. Effect of outlier sources in image difference noise

In this section, we present an analysis of the contribution of sources far away from the field center to the (differential) noise of images made at the NCP. In fact, this analysis can be extended to any interferometric observation. We show that our ignorance of these sources indeed act as an additional source of noise.

Consider an elementary interferometer. The visibility V(u_p,v_p,w_p) at coordinates u_p,v_p,w_p on the uv plane is $\begin{array}{ccc}\mathit{V}\mathrm{\left(}{\mathit{u}}_{\mathit{p}}\mathit{,}{\mathit{v}}_{\mathit{p}}\mathit{,}{\mathit{w}}_{\mathit{p}}\mathrm{\right)}\mathrm{=}& & \\ & & \mathrm{\int }\mathrm{\int }\mathit{S}\mathrm{\left(}\mathit{l,m}\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}\mathit{j}\mathrm{2}\mathit{\pi }\frac{\mathit{f}}{\mathit{c}}\left({\mathit{u}}_{\mathit{p}}\mathit{l}\mathrm{+}{\mathit{v}}_{\mathit{p}}\mathit{m}\mathrm{+}{\mathit{w}}_{\mathit{p}}\mathrm{(}\sqrt{\mathrm{1}\mathrm{-}{\mathit{l}}^{\mathrm{2}}\mathrm{-}{\mathit{m}}^{\mathrm{2}}}\mathrm{-}\mathrm{1}\mathrm{)}\right)}\frac{\mathrm{d}\mathit{l}\mathrm{d}\mathit{m}}{\sqrt{\mathrm{1}\mathrm{-}{\mathit{l}}^{\mathrm{2}}\mathrm{-}{\mathit{m}}^{\mathrm{2}}}}\end{array}$(2)where S(l,m) is the sky flux density and l,m are the direction cosines. The frequency of the observation is f while the speed of light is c. We assume the sky to consist of a set of discrete sources, and we arrive at $\mathit{V}\mathrm{\left(}{\mathit{u}}_{\mathit{p}}\mathit{,}{\mathit{v}}_{\mathit{p}}\mathit{,}{\mathit{w}}_{\mathit{p}}\mathrm{\right)}\mathrm{=}\sum _{\mathit{q}}\mathit{I}\mathrm{\left(}{\mathit{l}}_{\mathit{q}}\mathit{,}{\mathit{m}}_{\mathit{q}}\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}\mathit{j}\mathrm{2}\mathit{\pi }\frac{\mathit{f}}{\mathit{c}}\left({\mathit{u}}_{\mathit{p}}{\mathit{l}}_{\mathit{q}}\mathrm{+}{\mathit{v}}_{\mathit{p}}{\mathit{m}}_{\mathit{q}}\mathrm{+}{\mathit{w}}_{\mathit{p}}\left(\sqrt{\mathrm{1}\mathrm{-}{\mathit{l}}_{\mathit{q}}^{\mathrm{2}}\mathrm{-}{\mathit{m}}_{\mathit{q}}^{\mathrm{2}}}\mathrm{-}\mathrm{1}\right)\right)}$(3)where I(l,m) is the sky intensity.

Considering a bandwidth of Δ for smearing, we have the smeared value of V(u_p,v_p,w_p) around frequency f₀ as $\mathit{V̅}\mathrm{\left(}{\mathit{u}}_{\mathit{p}}\mathit{,}{\mathit{v}}_{\mathit{p}}\mathit{,}{\mathit{w}}_{\mathit{p}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{\Delta }}{\mathrm{\int }}_{{\mathit{f}}_{\mathrm{0}}\mathrm{-}\mathrm{\Delta }\mathit{/}\mathrm{2}}^{{\mathit{f}}_{\mathrm{0}}\mathrm{+}\mathrm{\Delta }\mathit{/}\mathrm{2}}\mathit{V}\mathrm{\left(}{\mathit{u}}_{\mathit{p}}\mathit{,}{\mathit{v}}_{\mathit{p}}\mathit{,}{\mathit{w}}_{\mathit{p}}\mathrm{\right)}\mathrm{d}\mathit{f}$(4)and assuming I(l_q,m_q) variation is small over this bandwidth, this reduces (3) to $\begin{array}{ccc}\mathit{V̅}\mathrm{\left(}{\mathit{u}}_{\mathit{p}}\mathit{,}{\mathit{v}}_{\mathit{p}}\mathit{,}{\mathit{w}}_{\mathit{p}}\mathrm{\right)}\mathrm{=}& & \\ & & \sum _{\mathit{q}}\mathit{I}\mathrm{\left(}{\mathit{l}}_{\mathit{q}}\mathit{,}{\mathit{m}}_{\mathit{q}}\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}\mathit{j}\mathrm{2}\mathit{\pi }\frac{\mathit{f}}{\mathit{c}}\left({\mathit{u}}_{\mathit{p}}{\mathit{l}}_{\mathit{q}}\mathrm{+}{\mathit{v}}_{\mathit{p}}{\mathit{m}}_{\mathit{q}}\mathrm{+}{\mathit{w}}_{\mathit{p}}\mathrm{(}\sqrt{\mathrm{1}\mathrm{-}{\mathit{l}}_{\mathit{q}}^{\mathrm{2}}\mathrm{-}{\mathit{m}}_{\mathit{q}}^{\mathrm{2}}}\mathrm{-}\mathrm{1}\mathrm{)}\right)}\\ & & \mathrm{\times }\mathrm{sinc}\left(\mathit{\pi }\frac{\mathrm{\Delta }}{\mathit{c}}\left({\mathit{u}}_{\mathit{p}}{\mathit{l}}_{\mathit{q}}\mathrm{+}{\mathit{v}}_{\mathit{p}}{\mathit{m}}_{\mathit{q}}\mathrm{+}{\mathit{w}}_{\mathit{p}}\left(\sqrt{\mathrm{1}\mathrm{-}{\mathit{l}}_{\mathit{q}}^{\mathrm{2}}\mathrm{-}{\mathit{m}}_{\mathit{q}}^{\mathrm{2}}}\mathrm{-}\mathrm{1}\right)\right)\right)\mathit{.}\end{array}$Let M denote the number of samples in the uv plane and N denote the number of sources in the sky. We represent the visibilities at all points on the uv plane in vector form as $\begin{array}{}{\mathrm{􏽥}}^{}\\ {b}\end{array}$ (size M × 1) and the intensities of all discrete sources in the sky as b (size N × 1) where $\begin{array}{}\mathrm{􏽥}\\ {b}\end{array}\mathrm{=}\left[\begin{array}{c}\\ \mathit{V̅}\mathrm{\left(}{\mathit{u}}_{\mathrm{1}}\mathit{,}{\mathit{v}}_{\mathrm{1}}\mathit{,}{\mathit{w}}_{\mathrm{1}}\mathrm{\right)}\\ \mathit{V̅}\mathrm{\left(}{\mathit{u}}_{\mathrm{2}}\mathit{,}{\mathit{v}}_{\mathrm{2}}\mathit{,}{\mathit{w}}_{\mathrm{2}}\mathrm{\right)}\\ \mathit{...}\\ \mathit{V̅}\mathrm{\left(}{\mathit{u}}_{\mathit{M}}\mathit{,}{\mathit{v}}_{\mathit{M}}\mathit{,}{\mathit{w}}_{\mathit{M}}\mathrm{\right)}\end{array}\right]\mathit{,}{b}\mathrm{=}\left[\begin{array}{c}\\ \mathit{I}\mathrm{\left(}{\mathit{l}}_{\mathrm{1}}\mathit{,}{\mathit{m}}_{\mathrm{1}}\mathrm{\right)}\\ \mathit{I}\mathrm{\left(}{\mathit{l}}_{\mathrm{2}}\mathit{,}{\mathit{m}}_{\mathrm{2}}\mathrm{\right)}\\ \mathit{...}\\ \mathit{I}\mathrm{\left(}{\mathit{l}}_{\mathit{N}}\mathit{,}{\mathit{m}}_{\mathit{N}}\mathrm{\right)}\end{array}\right]\mathit{.}$(7)We can relate $\begin{array}{}{\mathrm{􏽥}}^{}\\ {b}\end{array}$ and b as $\begin{array}{}\mathrm{􏽥}\\ {b}\end{array}\mathrm{=}{\mathbf{T}}_{{\mathit{f}}_{\mathrm{0}}}{b}$(8)where the elements in the matrix T_f0 (size M × N) are given as $\begin{array}{ccc}\mathrm{[}{\mathbf{T}}_{{\mathit{f}}_{\mathrm{0}}}{\mathrm{]}}_{\mathit{pq}}& \mathrm{=}& {\mathrm{e}}^{\mathrm{-}\mathit{j}\mathrm{2}\mathit{\pi }\frac{{\mathit{f}}_{\mathrm{0}}}{\mathit{c}}\left({\mathit{u}}_{\mathit{p}}{\mathit{l}}_{\mathit{q}}\mathrm{+}{\mathit{v}}_{\mathit{p}}{\mathit{m}}_{\mathit{q}}\mathrm{+}{\mathit{w}}_{\mathit{p}}\mathrm{(}\sqrt{\mathrm{1}\mathrm{-}{\mathit{l}}_{\mathit{q}}^{\mathrm{2}}\mathrm{-}{\mathit{m}}_{\mathit{q}}^{\mathrm{2}}}\mathrm{-}\mathrm{1}\mathrm{)}\right)}\\ & & \mathrm{\times }\mathrm{sinc}\left(\mathit{\pi }\frac{\mathrm{\Delta }}{\mathit{c}}\left({\mathit{u}}_{\mathit{p}}{\mathit{l}}_{\mathit{q}}\mathrm{+}{\mathit{v}}_{\mathit{p}}{\mathit{m}}_{\mathit{q}}\mathrm{+}{\mathit{w}}_{\mathit{p}}\left(\sqrt{\mathrm{1}\mathrm{-}{\mathit{l}}_{\mathit{q}}^{\mathrm{2}}\mathrm{-}{\mathit{m}}_{\mathit{q}}^{\mathrm{2}}}\mathrm{-}\mathrm{1}\right)\right)\right)\mathit{.}\end{array}$(9)Let D be the number of pixels of the image in which the noise is calculated. Consider the construction of the image pixels given by $\begin{array}{}{\mathrm{􏽢}}^{}\\ {b}\end{array}$ (size D × 1) $\begin{array}{}\mathrm{􏽢}\\ {b}\end{array}\mathrm{=}\left[\begin{array}{c}\\ \mathit{I}\mathrm{\left(}\begin{array}{}\mathrm{˜}\\ {\mathit{l}}_{\mathrm{1}}\end{array}\mathit{,}\begin{array}{}\mathrm{˜}\\ {\mathit{m}}_{\mathrm{1}}\end{array}\mathrm{\right)}\\ \mathit{I}\mathrm{\left(}\begin{array}{}\mathrm{˜}\\ {\mathit{l}}_{\mathrm{2}}\end{array}\mathit{,}\begin{array}{}\mathrm{˜}\\ {\mathit{m}}_{\mathrm{2}}\end{array}\mathrm{\right)}\\ \mathit{...}\\ \mathit{I}\mathrm{\left(}\begin{array}{}\mathrm{˜}\\ {\mathit{l}}_{\mathit{D}}\end{array}\mathit{,}\begin{array}{}\mathrm{˜}\\ {\mathit{m}}_{\mathit{D}}\end{array}\mathrm{\right)}\end{array}\right]$(10)from the observed data $\begin{array}{}{\mathrm{˜}}_{}\\ {b}\end{array}$. We can write this as $\begin{array}{}\mathrm{􏽢}\\ {b}\end{array}\mathrm{=}\begin{array}{}\mathrm{􏽥}\\ {\mathbf{T}}_{{\mathit{f}}_{\mathrm{0}}}^{†}\end{array}\begin{array}{}\mathrm{􏽥}\\ {b}\end{array}$(11)where the elements in the matrix $\begin{array}{}{\mathrm{􏽥}}_{}\\ {\mathbf{T}}_{{\mathit{f}}_{\mathrm{0}}}\end{array}$ (size M × D) are given as $\mathrm{[}\begin{array}{}\mathrm{􏽥}\\ {\mathbf{T}}_{{\mathit{f}}_{\mathrm{0}}}\end{array}{\mathrm{]}}_{\mathit{pq}}\mathrm{=}{\mathrm{e}}^{\mathrm{-}\mathit{j}\mathrm{2}\mathit{\pi }\frac{{\mathit{f}}_{\mathrm{0}}}{\mathit{c}}\left({\mathit{u}}_{\mathit{p}}{\mathit{l}}_{\mathit{q}}\mathrm{+}{\mathit{v}}_{\mathit{p}}{\mathit{m}}_{\mathit{q}}\mathrm{+}{\mathit{w}}_{\mathit{p}}\mathrm{(}\sqrt{\mathrm{1}\mathrm{-}{\mathit{l}}_{\mathit{q}}^{\mathrm{2}}\mathrm{-}{\mathit{m}}_{\mathit{q}}^{\mathrm{2}}}\mathrm{-}\mathrm{1}\mathrm{)}\right)}$(12)and in the reconstruction, no smearing is assumed. Note also that the sets ℒ and $\begin{array}{}{\mathrm{􏽥}}_{}\\ \mathrm{ℒ}\end{array}$ that denote the positions of the outlier sources and the positions of the pixels $\mathrm{ℒ}\mathrm{=}\mathrm{\left\{}\mathrm{\right(}{\mathit{l}}_{\mathrm{1}}\mathit{,}{\mathit{m}}_{\mathrm{1}}\mathrm{\left)}\mathit{,}\mathit{...}\mathit{,}\mathrm{\right(}{\mathit{l}}_{\mathit{N}}\mathit{,}{\mathit{m}}_{\mathit{N}}\mathrm{\left)}\mathrm{\right\}}\mathit{,}\begin{array}{}\mathrm{􏽥}\\ \mathrm{ℒ}\end{array}\mathrm{=}\mathrm{\left\{}\mathrm{\right(}\begin{array}{}\mathrm{˜}\\ {\mathit{l}}_{\mathrm{1}}\end{array}\mathit{,}\begin{array}{}\mathrm{˜}\\ {\mathit{m}}_{\mathrm{1}}\end{array}\mathrm{\left)}\mathit{,}\mathit{...}\mathit{,}\mathrm{\right(}\begin{array}{}\mathrm{˜}\\ {\mathit{l}}_{\mathit{D}}\end{array}\mathit{,}\begin{array}{}\mathrm{˜}\\ {\mathit{m}}_{\mathit{D}}\end{array}\mathrm{\left)}\mathrm{\right\}}$(13)have no relation to each other. In general the pixels coordinates $\begin{array}{}{\mathrm{􏽥}}_{}\\ \mathrm{ℒ}\end{array}$ where we calculate the differential noise, are on a regular grid. Using (8) and (11), we get $\begin{array}{}\mathrm{􏽢}\\ {b}\end{array}\mathrm{=}\begin{array}{}\mathrm{􏽥}\\ {\mathbf{T}}_{{\mathit{f}}_{\mathrm{0}}}^{†}\end{array}{\mathbf{T}}_{{\mathit{f}}_{\mathrm{0}}}{b}$(14)and the difference image at frequencies f₁ and f₀ can be given as ${e}\mathrm{=}\mathrm{(}\begin{array}{}\mathrm{􏽥}\\ {\mathbf{T}}_{{\mathit{f}}_{\mathrm{1}}}^{†}\end{array}{\mathbf{T}}_{{\mathit{f}}_{\mathrm{1}}}\mathrm{-}\begin{array}{}\mathrm{􏽥}\\ {\mathbf{T}}_{{\mathit{f}}_{\mathrm{0}}}^{†}\end{array}{\mathbf{T}}_{{\mathit{f}}_{\mathrm{0}}}\mathrm{)}{b}$(15)assuming the variation of b is negligible within this bandwidth. The noise variance in the difference image is proportional to ∥ e ∥ ² and the average noise variance per pixel is ∥ e ∥ ²/D and the standard deviation is $\sqrt{\mathrm{|}{e}{\mathrm{\parallel }}^{\mathrm{2}}\mathit{/}\mathit{D}}$.

Considering a random distribution of outlier sources, we can also find the expected value of ∥ e ∥ ² as $\mathit{E}\mathrm{\{}\mathrm{\parallel }{e}{\mathrm{\parallel }}^{\mathrm{2}}\mathrm{\}}\mathrm{=}\mathit{E}\mathrm{\left\{}{{e}}^{\mathit{T}}{e}\mathrm{\right\}}\mathrm{=}\mathit{E}\mathrm{\left\{}{{b}}^{\mathit{T}}{\mathbf{T}}_{{\mathit{f}}_{\mathrm{1}}{\mathit{f}}_{\mathrm{0}}}{b}\mathrm{\right\}}$(16)where T_f1f0 (size N × N) is given by ${\mathbf{T}}_{{\mathit{f}}_{\mathrm{1}}{\mathit{f}}_{\mathrm{0}}}\mathrm{=}\mathrm{Re}\mathrm{\left(}\mathrm{\right(}\begin{array}{}\mathrm{􏽥}\\ {\mathbf{T}}_{{\mathit{f}}_{\mathrm{1}}}^{†}\end{array}{\mathbf{T}}_{{\mathit{f}}_{\mathrm{1}}}\mathrm{-}\begin{array}{}\mathrm{􏽥}\\ {\mathbf{T}}_{{\mathit{f}}_{\mathrm{0}}}^{†}\end{array}{\mathbf{T}}_{{\mathit{f}}_{\mathrm{0}}}{\mathrm{)}}^{\mathit{H}}\mathrm{(}\begin{array}{}\mathrm{􏽥}\\ {\mathbf{T}}_{{\mathit{f}}_{\mathrm{1}}}^{†}\end{array}{\mathbf{T}}_{{\mathit{f}}_{\mathrm{1}}}\mathrm{-}\begin{array}{}\mathrm{􏽥}\\ {\mathbf{T}}_{{\mathit{f}}_{\mathrm{0}}}^{†}\end{array}{\mathbf{T}}_{{\mathit{f}}_{\mathrm{0}}}{\mathrm{\left)}}^{\mathrm{\right)}}$(17)and note that T_f1f0 is symmetric positive semi-definite because of the factorization as above. We also consider T_f1f0 to be real because the sky image is real.

We find an upper and a lower bound for E {∥ e ∥ ²}  as $\frac{\mathrm{1}}{{\mathit{N}}^{\mathrm{2}}}\mathit{E}\mathrm{\left\{}\mathrm{\right|}{\mathbf{T}}_{{\mathit{f}}_{\mathrm{1}}{\mathit{f}}_{\mathrm{0}}}\mathrm{\left|}\mathrm{\right\}}\mathit{E}\mathrm{\left\{}\mathrm{\right|}{b}{\mathrm{|}}^{\mathrm{2}}\mathrm{\}}\mathrm{\le }\mathit{E}\mathrm{\{}\mathrm{\parallel }{e}{\mathrm{\parallel }}^{\mathrm{2}}\mathrm{\}}\mathrm{\le }\mathit{E}\mathrm{\{}\max \mathrm{\left(}\mathrm{diag}\mathrm{\right(}{\mathbf{T}}_{{\mathit{f}}_{\mathrm{1}}{\mathit{f}}_{\mathrm{0}}}\mathrm{\left)}\mathrm{\right)}\mathrm{\left\}}\mathit{E}\mathrm{\right\{}\mathrm{\left|}{b}{\mathrm{|}}^{\mathrm{2}}\mathrm{\right\}}$(18)where |T_f1f0| is the sum of all elements in T_f1f0 and |b| is the sum of all elements in b. The diagonal entry with the highest magnitude in T_f1f0 is given by max(diag(T_f1f0)). The proof and the underlying assumptions can be found in Appendix A.

We draw a few conclusions from (18): (i) by properly selecting the uv coverage and the frequencies f₁ and f₀ and also imaging weights (in the derivation natural weights are assumed) and image pixel sizes, we can change T_f1f0; (ii) however, |b|² is entirely determined by the intensities of the outlier sources and the only way to minimize |b|² is by suppressing the beam sidelobe level or by subtracting the outlier sources, as we have done for CasA; (iii) we get the lowest value for the noise when all outlier sources have intensities that are equally distributed while we get the highest noise when there is one very bright source (see Appendix A).

In order to relate (18) to NCP observations, we need to estimate |b|² due to the outlier sources in the images. For the particular case of the NCP, we select the short baselines (<1200 wavelengths) from the uv coverage in Fig. 4. The corresponding image of the FOV is shown in Fig. 10. With the same pixel size of 35′′, we have made an image of 6400 × 6400 pixels. From this image, we have selected all pixels that have a flux >1.2 mJy that are outside the FOV and we show the selected pixels in Fig. 21.

We consider the vector b to be formed by the pixels selected in Fig. 21. After forming this vector, we have calculated |b| for different frequencies (using images made at those frequencies) as shown in Fig. 22. We have also fitted a model for |b| (the red line) which is $\mathrm{\left|}{b}\mathrm{\right|}\mathrm{=}\mathrm{89.32}{\left(\frac{\mathit{f}}{\mathrm{150}\mathrm{\times }{\mathrm{10}}^{\mathrm{6}}}\right)}^{-3.687}\mathrm{+}\mathrm{70.}$(19)Using the model for |b| given by (19), we can approximate E {|b|²} . Next, using this approximation, we simulate (18), for various image differencing frequencies. We simulate a sky consisting of N = 500 outlier sources. Note that some of these sources will have intensities below the receiver noise level at any given frequency. But we make an important distinction that what matters is the cumulative effect of all the sources taken together, and not their individual contributions. The pixel grid for the simulation at the center has dimensions 100 × 100, so D = 10   000. The uv coverage is chosen to be the core (<1200 wavelengths) uv coverage from Fig. 4. For a longest baseline of 1200 wavelengths, the resolution is about 3′ and the pixel size is chosen to be large enough (10′) so that we are not dominated by sampling errors (Yatawatta 2010) and that the noise is not correlated between pixels. The geometry of one realization of the simulation is shown in Fig. 23. Note that the positions of the outlier sources are randomly varied in each realization.

We have done simulations for three different values of smearing bandwidth, Δ = 200 kHz, Δ = 20 kHz, and Δ = 2 kHz. The bandwidth where the difference is taken is also equal to Δ. For each value of Δ, we generated 100 different sky realizations as in Fig. 23.

We have shown the upper and lower bounds for the standard deviation of the noise due to outlier sources per pixel in Fig. 24. We assume the noise contribution due to outlier sources at the highest frequency (160 MHz) is zero. As expected, as the value of Δ decreases, the noise due to far away sources is also reduced. So for low values of Δ, the result will be dominated by the receiver noise. In future processing, we intend to use Δ = 20 kHz and lower.

6. Discussion/conclusions

We have presented the first deep imaging results made with LOFAR that reach a noise level of about 100 μJy/PSF both in total intensity and in polarization. Note that for calculation of the noise in polarization, we have only used one night of data while for calculation of the noise in total intensity, we have used all three nights of data listed in Table 1. Using all three nights of data for polarization analysis requires correction for the effects due to differential Faraday rotation which varies from night to night and will be done in future processing. These results act as a precursor of results with more dedicated EoR observations that are about to begin with LOFAR.

We briefly discuss the differences between previous WSRT observations of the NCP (72 h, 14 MHz) and the LOFAR observations (10 h effective, 48 MHz) presented in this paper. Due to having only a 2.7 km longest baseline, the WSRT observations were limited by classical source confusion in the continuum. Moreover, the elevated location of the WSRT receivers, and the interference from the building, made them much more prone to RFI. The WSRT observations also used only 14 stations while LOFAR had about 40 stations Furthermore, the data processing capabilities are significantly better than what was available to process WSRT observations. Therefore, despite the much longer integration time used in the WSRT observations we have obtained images that are much better in terms of resolution as well as noise. In fact, we have presented the deepest interferometric images ever made at this frequency range, to this date. The observed noise levels – in high resolution images, linear polarization and at the edges of low resolution wide field images, are within a factor of 1.4 from the theoretical noise that can be achieved based on the nominal system equivalent flux density and the effective integration time.

Nonetheless, we can see limitations in the current results that can be attributed to the choices made in our processing:

• We have only a partially correct sky model that contains about 500discrete sources. However, we see much more than this in ourimages but at present we are limited by the resolution that is onlyabout 12′′. We therefore intend to use more data with longer baselines (even international LOFAR stations) to update our sky model.

• We have thus far not seen (and ignored) the effect of (spatial) differential Faraday rotation (DFR) due to the ionosphere in our data. Note that even though we calibrate for full Jones matrices using SAGECal, we do not correct the data using the solutions. However, we shall soon incorporate taking DFR into account in our calibration (Yatawatta 2012). We believe this will be crucial when we get data with significantly longer (80 km) baselines.

• In the processing of our data we were limited by computational resources. We have therefore invested significant effort in optimizing the performance of our calibration pipeline using GPUs to accelerate our processing. This will enable us to employ sophisticated algorithms yielding better results. For instance, with a slightly enhanced version of SAGECal (Yatawatta et al. 2012), we intend to process data at a finer spectral resolution, therefore eliminating the bandwidth smearing visible at the edge of the FOV in our images. This will also decrease the effect of outlier sources, as discussed in Sect. 5.

In summary, we do not see any major obstacles to prevent us from going considerably deeper through much longer multiple-nights integration. In forthcoming publications we will investigate the properties of the sources detected in the images, make a detailed study of effects due to the Galactic foreground and determine limits on signals from the epoch of reionization.

---

Acknowledgments

LOFAR, designed and constructed by ASTRON, has facilities in several countries, that are owned by various parties (each with their own funding sources), and that are collectively operated by the International LOFAR Telescope (ILT) foundation under a joint scientific policy. Chiara Ferrari acknowledges financial support by the Agence Nationale de la Recherche through grant ANR-09-JCJC-0001-01.

References

Appendix A: Proof of Eq. (18)

A.1. Mathematical background

Consider a vector of real numbers x of size N × 1 with x_i denoting the ith element, ${x}\mathrm{=}\mathrm{[}{\mathit{x}}_{\mathrm{1}}\mathit{,}{\mathit{x}}_{\mathrm{2}}\mathit{,}\mathit{...}\mathit{,}{\mathit{x}}_{\mathit{N}}{\mathrm{]}}^{\mathit{T}}\mathit{.}$(A.1)We can order the elements in x in descending order and construct a new vector $\mathrm{\left[}{x}\mathrm{\right]}\mathrm{=}\mathrm{[}{\mathit{x}}_{\mathrm{\left[}\mathrm{1}\mathrm{\right]}}\mathit{,}{\mathit{x}}_{\mathrm{\left[}\mathrm{2}\mathrm{\right]}}\mathit{,}\mathit{...}\mathit{,}{\mathit{x}}_{\mathrm{\left[}\mathit{N}\mathrm{\right]}}{\mathrm{]}}^{\mathit{T}}$(A.2)where x _[1]  ≥ x _[2]  ≥ ... ≥ x _[N] . Given two vectors of real numbers, x and y of size N × 1, we say x is majorized by y, or x ≺ y when $\begin{array}{ccc}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{k}}{\mathit{x}}_{\mathrm{\left[}\mathit{i}\mathrm{\right]}}\mathrm{\le }\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{k}}{\mathit{y}}_{\mathrm{\left[}\mathit{i}\mathrm{\right]}}\mathit{,}\mathrm{1}\mathrm{\le }\mathit{k}\mathrm{\le }\mathit{N,}\mathrm{and}& & \end{array}$(A.3)$\begin{array}{ccc}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{x}}_{\mathrm{\left[}\mathit{i}\mathrm{\right]}}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{y}}_{\mathrm{\left[}\mathit{i}\mathrm{\right]}}& & \end{array}$are satisfied (Marshall et al. 2011).

Consider the quadratic form $\mathit{\phi }\mathrm{\left(}{x}\mathrm{\right)}\mathrm{=}{{x}}^{\mathit{T}}\mathbf{A}{x}$(A.4)where A is an N × N real symmetric matrix. If the elements in A satisfy $\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathit{i}}\mathrm{\left(}\mathrm{\right[}\mathbf{A}{\mathrm{]}}_{\mathit{k,j}}\mathrm{-}\mathrm{\left[}\mathbf{A}{\mathrm{]}}_{\mathit{k}\mathrm{+}\mathrm{1}\mathit{,j}}\mathrm{\right)}\mathrm{\ge }\mathrm{0}\mathit{,}\mathit{i}\mathrm{=}\mathrm{1}\mathit{,}\mathit{...}\mathit{,N,}\mathit{k}\mathrm{=}\mathrm{1}\mathit{,}\mathit{...}\mathit{,N}\mathrm{-}\mathrm{1}$(A.5)then we say φ(x) is Schur-convex (Marshall et al. 2011).

If φ(x) is Schur-convex and with two vectors of positive real numbers, x and y, with x ≺ y, then $\mathit{\phi }\mathrm{\left(}{x}\mathrm{\right)}\mathrm{\le }\mathit{\phi }\mathrm{\left(}{y}\mathrm{\right)}\mathit{.}$(A.6)

A.2. Assumptions

Using conditional mean, we rewrite (16) as $\mathit{E}\mathrm{\left\{}{{b}}^{\mathit{T}}{\mathbf{T}}_{{\mathit{f}}_{\mathrm{1}}{\mathit{f}}_{\mathrm{0}}}{b}\mathrm{\right\}}\mathrm{=}{\mathit{E}}_{\mathrm{\left|}{b}\mathrm{\right|}}\left\{\mathit{E}\mathrm{\{}\frac{{{b}}^{\mathit{T}}}{\mathrm{\left|}{b}\mathrm{\right|}}{\mathbf{T}}_{{\mathit{f}}_{\mathrm{1}}{\mathit{f}}_{\mathrm{0}}}\frac{{b}}{\mathrm{\left|}{b}\mathrm{\right|}}\mathrm{|}{b}{\mathrm{|}}^{\mathrm{2}}\mathrm{|}\mathrm{|}{b}\mathrm{|}\mathrm{\}}\right\}\mathit{,}$(A.7)where E_|b| {.}  is expectation with respect to |b|, i.e. the sum of all elements in b.

We make to following assumptions in order to proceed:

• The source fluxes (elements in b) and their positions are both random variables. We assume that the probability distribution for the fluxes is statistically independent of the probability distribution for the source positions. In other words, the statistical properties of b in (A.7) are independent of the statistical properties of T_f1f0.

• We assume T_f1f0 to satisfy the properties given in (A.5). More precisely, we assume that there exists a permutation of the rows and the columns of T_f1f0 to satisfy this condition. Since we did not impose any ordering on the pixel positions in (13) this permutation is feasible. We do not need to find this permutation since the end result does not depend on it. But we assume that T_f1f0 is permuted such that the largest diagonal entry is at the 1-st row and the 1-st column.

A.3. Derivation of bounds

Let ${z}\mathrm{=}\frac{{b}}{\mathrm{\left|}{b}\mathrm{\right|}}$ and we see that z is a vector of N positive elements that add up to 1, i.e. |z| = 1. Consider the following vectors $\begin{array}{ccc}\underline{{z}}& \mathrm{=}& \mathrm{[}\mathrm{1}\mathit{/}\mathit{N,}\mathrm{1}\mathit{/}\mathit{N,}\mathit{...}\mathit{,}\mathrm{1}\mathit{/}\mathit{N}{\mathrm{]}}^{\mathit{T}}\mathit{,}\mathrm{and}\\ \overline{{z}}& \mathrm{=}& \mathrm{[}\mathrm{1}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}\mathit{...}\mathit{,}\mathrm{0}{\mathrm{]}}^{\mathit{T}}\mathit{.}\end{array}$(A.8)We see that from the properties of (A.3), $\underline{{z}}\mathrm{\prec }{z}\mathrm{\prec }\overline{{z}}\mathit{.}$(A.9)Together with (A.6) and (A.9), and the assumptions on T_f1f0, we get ${\underline{{z}}}^{\mathit{T}}{\mathbf{T}}_{{\mathit{f}}_{\mathrm{1}}{\mathit{f}}_{\mathrm{0}}}\underline{{z}}\mathrm{\le }{{z}}^{\mathit{T}}{\mathbf{T}}_{{\mathit{f}}_{\mathrm{1}}{\mathit{f}}_{\mathrm{0}}}{z}\mathrm{\le }{\overline{{z}}}^{\mathit{T}}{\mathbf{T}}_{{\mathit{f}}_{\mathrm{1}}{\mathit{f}}_{\mathrm{0}}}\overline{{z}}\mathit{.}$(A.10) We can simplify (A.10) as $\frac{\mathrm{1}}{{\mathit{N}}^{\mathrm{2}}}\mathrm{\left|}{\mathbf{T}}_{{\mathit{f}}_{\mathrm{1}}{\mathit{f}}_{\mathrm{0}}}\mathrm{\right|}\mathrm{\le }{{z}}^{\mathit{T}}{\mathbf{T}}_{{\mathit{f}}_{\mathrm{1}}{\mathit{f}}_{\mathrm{0}}}{z}\mathrm{\le }\max \mathrm{\left(}\mathrm{diag}\mathrm{\right(}{\mathbf{T}}_{{\mathit{f}}_{\mathrm{1}}{\mathit{f}}_{\mathrm{0}}}\mathrm{\left)}\mathrm{\right)}$(A.11) where |T_f1f0| is the sum of all elements in T_f1f0 and the diagonal entry with the highest magnitude of T_f1f0 (at row 1 and Col. 1) is given by max(diag(T_f1f0)).

Substituting (A.11) back to (A.7) and considering the independence of b and T_f1f0, we get (18).

One additional point to note is that the lower bound is obtained when ${z}\mathrm{=}\underline{{z}}$, i.e. when the elements in b are uniformly distributed with equal values. Moreover, the upper bound is obtained when ${z}\mathrm{=}\overline{{z}}$, i.e. the total value of the elements in b is concentrated at one element. Reverting back to the fluxes of the sources, the lower bound is obtained when all sources have equal fluxes while the upper bound is obtained when only one source has the total flux (or when there is one very bright source).

All Tables

All Figures