# The LOFAR Two-meter Sky Survey: Deep Fields Data Release 1

C. Tasse, T. Shimwell, M. J. Hardcastle, S. P. O'sullivan, R. van Weeren, P.

N. Best, L. Bester, B. Hugo, O. Smirnov, J. Sabater, et al.

► To cite this version:

C. Tasse, T. Shimwell, M. J. Hardcastle, S. P. O'sullivan, R. van Weeren, et al. The LOFAR Two-meter Sky Survey: Deep Fields Data Release 1: I. Direction-dependent calibration and imaging. Astronomy and Astrophysics - A&A, 2021, 648, pp.A1. 10.1051/0004-6361/202038804. hal-03192293

HAL Id: hal-03192293

https://hal.science/hal-03192293

Submitted on 7 Apr 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# The LOFAR Two-meter Sky Survey: Deep Fields Data Release 1

## I. Direction-dependent calibration and imaging

C. Tasse<sup>1,2,3</sup>, T. Shimwell<sup>4,5</sup>, M. J. Hardcastle<sup>6</sup>, S. P. O'Sullivan<sup>7</sup>, R. van Weeren<sup>5</sup>, P. N. Best<sup>8</sup>, L. Bester<sup>9,2</sup>, B. Hugo<sup>9,2</sup>, O. Smirnov<sup>2,9</sup>, J. Sabater<sup>8</sup>, G. Calistro-Rivera<sup>10</sup>, F. de Gasperin<sup>11</sup>, L. K. Morabito<sup>12</sup>, H. Röttgering<sup>5</sup>, W. L. Williams<sup>5</sup>, M. Bonato<sup>13,14,15</sup>, M. Bondi<sup>13</sup>, A. Botteon<sup>5,13</sup>, M. Brüggen<sup>11</sup>, G. Brunetti<sup>13</sup>, K. T. Chyży<sup>16</sup>, M. A. Garrett<sup>17,5</sup>, G. Gürkan<sup>18</sup>, M. J. Jarvis<sup>19,20</sup>, R. Kondapally<sup>8</sup>, S. Mandal<sup>5</sup>, I. Prandoni<sup>13</sup>, A. Repetti<sup>21</sup>, E. Retana-Montenegro<sup>22</sup>, D. J. Schwarz<sup>23</sup>, A. Shulevski<sup>24</sup>, and Y. Wiaux<sup>21</sup>

Received 1 July 2020 / Accepted 10 November 2020

## ABSTRACT

The Low Frequency Array (LOFAR) is an ideal instrument to conduct deep extragalactic surveys. It has a large field of view and is sensitive to large-scale and compact emission. It is, however, very challenging to synthesize thermal noise limited maps at full resolution, mainly because of the complexity of the low-frequency sky and the direction dependent effects (phased array beams and ionosphere). In this first paper of a series, we present a new calibration and imaging pipeline that aims at producing high fidelity, high dynamic range images with LOFAR High Band Antenna data, while being computationally efficient and robust against the absorption of unmodeled radio emission. We apply this calibration and imaging strategy to synthesize deep images of the Boötes and Lockman Hole fields at ~150 MHz, totaling ~80 and ~100 h of integration, respectively, and reaching unprecedented noise levels at these low frequencies of ~30 and ~23 µJy beam⁻<fnref n="1" /> in the inner ~3 deg<fnref n="2" />. This approach is also being used to reduce the LOTSS-wide data for the second data release.

**Key words.** techniques: interferometric – techniques: image processing – surveys – galaxies: active – galaxies: starburst – radio continuum: galaxies

## 1. Introduction

With its low observing frequency, wide fields of view, high sensitivity, large fractional bandwidth, and high spatial resolution, the Low Frequency Array (LOFAR, see van Haarlem et al. 2013)

is well suited to conduct deep extragalactic surveys. The LOFAR Surveys Key Science Project is using LOFAR to build tiered extragalactic surveys of different depths and areas and at frequencies ranging from ~30 to ~200 MHz. Specifically, the LOFAR LBA Sky Survey (LOLSS) aims at surveying the

A1, page 1 of 18

<footnote n="1">¹ GEPI, Observatoire de Paris, CNRS, Université Paris Diderot, 5 place Jules Janssen, 92190 Meudon, France
e-mail: cyril.tasse@obspm.fr</footnote>

<footnote n="2">² Centre for Radio Astronomy Techniques and Technologies, Department of Physics and Electronics, Rhodes University, Grahamstown 6140, South Africa</footnote>

<footnote n="3">³ USN, Observatoire de Paris, CNRS, PSL, UO, Nançay, France</footnote>

<footnote n="4">⁴ ASTRON, Netherlands Institute for Radio Astronomy, Oude Hoogeveensedijk 4, 7991 PD, Dwingeloo, The Netherlands</footnote>

<footnote n="5">⁵ Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands</footnote>

<footnote n="6">⁶ Centre for Astrophysics Research, University of Hertfordshire, College Lane, Hatfield AL10 9AB, UK</footnote>

<footnote n="7">⁷ School of Physical Sciences and Centre for Astrophysics & Relativity, Dublin City University, Glasnevin, D09 W6Y4, Ireland</footnote>

<footnote n="8">⁸ SUPA, Institute for Astronomy, Royal Observatory, Blackford Hill, Edinburgh, EH9 3HJ, UK</footnote>

<footnote n="9">⁹ South African Radio Astronomy Observatory, Observatory 7925, Cape Town, South Africa</footnote>

<footnote n="10">¹⁰ European Southern Observatory, Karl-Schwarzchild-Strasse 2, 85748 Garching bei München, Germany</footnote>

<footnote n="11">¹¹ Hamburger Sternwarte, University of Hamburg, Gojenbergsweg 112, 21029 Hamburg, Germany</footnote>

<footnote n="12">¹² Centre for Extragalactic Astronomy, Department of Physics, Durham University, Durham, DH1 3LE, UK</footnote>

<footnote n="13">¹³ INAF – Istituto di Radioastronomia, Via P. Gobetti 101, 40129 Bologna, Italy</footnote>

<footnote n="14">¹⁴ INAF – Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, 35122, Padova, Italy</footnote>

<footnote n="15">¹⁵ Italian ALMA Regional Centre, Via Gobetti 101, 40129 Bologna, Italy</footnote>

<footnote n="16">¹⁶ Astronomical Observatory, Jagiellonian University, ul. Orla 171, 30-244 Kraków, Poland</footnote>

<footnote n="17">¹⁷ Jodrell Bank Centre for Astrophysics, University of Manchester, Alan Turing Building, Oxford Road, M13 9PL, UK</footnote>

<footnote n="18">¹⁸ CSIRO Astronomy and Space Science, PO Box 1130, Bentley WA 6102, Australia</footnote>

<footnote n="19">¹⁹ Astrophysics, Department of Physics, Keble Road, Oxford, OX1 3RH, UK</footnote>

<footnote n="20">²⁰ Department of Physics & Astronomy, University of the Western Cape, Private Bag X17, Bellville, Cape Town 7535, South Africa</footnote>

<footnote n="21">²¹ Institute of Sensors, Signals and Systems, Heriot-Watt University, Edinburgh EH14 4AS, UK</footnote>

<footnote n="22">²² Astrophysics & Cosmology Research Unit, School of Mathematics, Statistics & Computer Science, University of KwaZulu-Natal, Durban 3690, South Africa</footnote>

<footnote n="23">²³ Fakultät für Physik, Universität Bielefeld, Postfach 100131, 33501 Bielefeld, Germany</footnote>

<footnote n="24">²⁴ Anton Pannekoek Institute for Astronomy, University of Amsterdam, Postbus 94249, 1090 GE Amsterdam, The Netherlands</footnote>
<figure><img src="image_1.png" /><figcaption>**Fig. 1.** Effective noise in the LOTSS-Deep Fields continuum maps as compared to other existing and future surveys. A spectral index of $\alpha = -0.7$ has been used to convert flux densities to the 1.4 GHz reference frequency.</figcaption></figure>

northern hemisphere using the LOFAR Low Band Antenna (LBA) antennas, while the LOFAR Two-metre Sky Survey (LOTSS) uses the High Band Antennas (HBAs). Its widest component (LOTSS-wide) has been described by Shimwell et al. (2017, 2019) and aims at reaching noise levels of $\le 100 \ \mu\text{Jy}\ \text{beam}^{-1}$ over the whole northern hemisphere. While the bright sources identified in LOTSS-wide are largely radio loud active galactic nuclei (AGN), the population of faint sources consists of star forming galaxies and radio quiet AGN (see Padovani 2016, and references therein). The LOTSS-Deep Fields targets noise levels of ultimately $\le 10 \ \mu\text{Jy}\ \text{beam}^{-1}$, thereby entering a new fainter, higher redshift regime where star forming galaxies and radio quiet AGN will outnumber the population of radio loud AGN, and thereby probing the evolution of those populations over cosmic time. Figure 1 is inspired by that of Smolčić et al. (2017a) and shows a sensitivity and surveyed area comparison between various existing and future surveys. These include the TIFR Giant Metrewave Radio Telescope Sky Survey (TGSS, Intema et al. 2017), the Faint Images of the Radio Sky at Twenty Centimeters (FIRST, Becker et al. 1995), the NRAO Very Large Array Sky Survey (NVSS, Condon et al. 1998), Very Large Array-COSMOS (VLA-COSMOS, Schinnerer et al. 2004; Smolčić et al. 2017b), the Very Large Array Sky Survey (VLASS, Lacy et al. 2020), Evolutionary Map of the Universe (EMU, Norris 2010), the Very Large Array-Spitzer Wide-Area Infrared Extragalactic Survey (VLA-SWIRE, Owen & Morrison 2008), SSA13 (Fomalont et al. 2006), Stripe 82 (Heywood et al. 2016), XXL (Butler et al. 2018, and references therein), DEEP2 (Mauch et al. 2020), the LOTSS first data release (DR1) (LOTSS-DR1, Shimwell et al. 2017), Hubble Deep Field - North (HDF-N, Richards 2000), the Westerbork Northern Sky Survey (WENSS, Rengelink et al. 1997), the GaLactic and Extra-Galactic All-Sky MWA (GLEAM, Wayth et al. 2015), and the Square Kilometer Array (SKA, Prandoni & Seymour 2015).

The depth of the LOTSS-Deep Fields is unlikely to be routinely surpassed at these low frequencies even into the era of the first phase of the SKA (Dewdney et al. 2009) because, although the SKA will have the raw sensitivity to easily reach such depths, the confusion noise of the SKA-low will likely increase the image rms to values exceeding the target depth of the LOTSS-deep images (see e.g., Zwart et al. 2015; Prandoni & Seymour 2015).

In order to construct the LOTSS-Deep Fields, we have selected the Boötes, Lockman Hole, and European Large Area Infrared Space Observatory Survey-North 1 (ELAIS-N1) fields, together with the North Ecliptic Pole (NEP). Each of them is covered by a wealth of multiwavelength data, which are necessary to derive photometric redshifts and low frequency radio luminosities, thereby providing an efficient way to estimate star formation rates (SFRs) in galaxies, for example. Together, these four multi-wavelength fields allow us to probe a total sky area of $\ge 30 \ \text{deg}^2$ in order to probe all galaxy environments at $z > 1$.

It is, however, quite challenging to make thermal noise limited images at low frequencies because of the presence of direction-dependent effects (DDEs), such as the ionospheric distortions, and because of the complex primary beam shapes of phased arrays. We have shown (Shimwell et al. 2019) that, using a novel set of calibration and imaging algorithms developed by Tasse (2014a), Smirnov & Tasse (2015), and Tasse et al. (2018), we are able to estimate and compensate for the DDEs and thus use LOFAR to produce thermal noise limited maps from eight-hour LOFAR observations in a systematic and automated way, while keeping the computational efficiency high enough to be able to deal with the high LOFAR data rates.

In this first paper of a series, we present an improved strategy to reach thermal noise limited images after hundreds of hours of integration on the Boötes and Lockman Hole extragalactic fields, reaching $\sim 30 \ \mu\text{Jy}\ \text{beam}^{-1}$ noise levels, while being more robust against absorbing faint unmodeled extended emission and dynamic range issues. In Sect. 2, we introduce the direction-dependent (DD) calibration and imaging problem together with the existing software that we use to tackle it. We describe our DDE calibration and imaging strategy (DDF-PIPELINE-V2) in Sect. 3 (to expand on this, in Appendix A we describe DDF-PIPELINE-V1, which was presented in detail in Shimwell et al. 2019). In Sect. 4, we use DDF-PIPELINE-V2 to synthesize deep images over the Boötes and Lockman Hole extragalactic fields, and present these deep low frequency images. The subsequent papers in this series will present the deeper ELAIS-N1 data products (Sabater et al. 2021), the multiwavelength cross matching (Kondapally et al. 2021), and the photometric redshifts and galaxy characterization (Duncan et al. 2021).

## 2. LOTSS and the third generation calibration and imaging problem

Calibration and imaging techniques have greatly evolved since the first radio interferometers became operational. First generation calibration is commonly referred to as direction-independent (DI) calibration, where calibration solutions are transferred to the target from an amplitude and/or phase calibrator field. Second generation calibration is the innovation, beginning in the mid-1970s, of using the target field to calibrate itself (self-calibration; Pearson & Readhead 1984). Third generation calibration and imaging consists in estimating and compensating for DDEs.

As mentioned above, it is challenging to synthesize high resolution thermal noise limited images with LOFAR (van Haarlem et al. 2013). Specifically, LOFAR (i) operates at very low frequencies ($\nu < 250 \ \text{MHz}$), (ii) has very large fields of view (full width at half maximum of 2–10 degrees), and (iii) combines short ($\sim 100 \ \text{m}$) and long ($\sim 2000 \ \text{km}$) baselines to provide sensitivity to both the compact and extended emission. Because of the presence of the ionosphere and the usage of phased array beams, the combination of (i) and (ii) makes the calibration problem direction-dependent by nature. In Sect. 2.1, we introduce
the mathematical formalism used throughout this paper, while in Sects. 2.2 and 2.3 we describe the algorithms and software used for DI and DD calibration and imaging.

## 2.1. Measurement equation formalism

The radio interferometry measurement equation (RIME, see Hamaker et al. 1996) describes how the underlying electric field coherence (the sky model), as well as the various direction-independent and direction-dependent Jones matrices (denoted $G$ and $J$ respectively), map to the measured visibilities. In the following, we consider the electric field in linear notation (along the X and Y axes) at frequency $v$ in direction $s = [l, m, n = \sqrt{1 - l^2 - m^2}]^T$ (where $T$ is the matrix transpose) and write the $4 \times 1$ sky coherency matrix as $x_{sv} = [XX, XY, YX, YY]^T_{sv}$. If $G_b = G_{qtv}^* \otimes G_{ptv}$ and $J_b^s = J_{qtv}^* \otimes J_{ptv}^s$ are the direction-independent and direction-dependent $4 \times 4$ Mueller matrices<fnref n="1" /> on a baseline $b \leftrightarrow \{pq,t\} \rightarrow [u, v, w]^T$ between antenna $p$ and $q$ at time $t$, then the four-visibility $v_b$ is given by

$$v_b = G_b \int_s \mathcal{J}_b^s \mathcal{B}_b^s x_{sv} k_b^s ds + n_b, \quad (1)$$

$$\text{with } k_b^s = \exp\left(-2i\pi\frac{v}{c}(\mathbf{b}_{pq,t}^T(s - s_0))\right) \quad (2)$$

$$\text{and } b_{pq,t} = \begin{bmatrix} u_{pq,t} \\ v_{pq,t} \\ w_{pq,t} \end{bmatrix} = \begin{bmatrix} u_{p,t} \\ v_{p,t} \\ w_{p,t} \end{bmatrix} - \begin{bmatrix} u_{q,t} \\ v_{q,t} \\ w_{q,t} \end{bmatrix} \quad (3)$$

$$\text{and } s = \begin{bmatrix} l \\ m \\ n \end{bmatrix} \text{ and } s_0 = \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix}, \quad (4)$$

where $c$ is the speed of light in the vacuum and $n_b$ is a $4 \times 1$ random matrix following a normal distribution $N(0, \sigma_b)$. Depending on the context, in the rest of this paper we will either make use of the antenna-based Jones matrices or the baseline-based Mueller matrices.

The elements of $G_b$ describe the direction-independent effects, such as the individual station electronics (the bandpass) or their clock drifts and offsets. The $J_b^s$ models the DDES, including the ionospheric distortion (phase shift, Faraday rotation, and scintillative decoherence) and phased array station beam that depend on time, frequency, and antenna. Importantly, in order to estimate the intrinsic flux densities, we used a description of the LOFAR station primary beam that is built from semi-analytic models<fnref n="2" /> and write it as $B_b^s$ in Eq. (1).

Solving for the third generation calibration and imaging problem consists of estimating the terms on the right-hand side of Eq. (1), namely the Mueller matrices $G_b$ and $J_b^s$ and the sky model $x_v$ from the set of visibilities $v_b$. Due to the bilinear structure of the RIME, instead of estimating all these parameters at once, the inversion of Eq. (1) is split into two steps. In the first step, the sky term $x_v$ is assumed to be constant and the Jones matrices are estimated. This step is referred to as “calibration” or as the DD-C-RIME system later in this text (or simply C-RIME, depending on the context). In the second step, the Jones matrices are assumed to be constant and the sky term $x_v$ is estimated. This step is commonly called “imaging” and is referred to as

solving the DD-I-RIME system later in the text. The C-RIME and I-RIME problems constitute two sub-steps in inverting the RIME system. We will later describe the idea of alternating between DD-C-RIME and DD-I-RIME as DD self-calibration.

While the vast majority of modern developments in the field of algorithmic research for radio interferometry aim at addressing either the DD-C-RIME (DD calibration, see Yatawatta et al. 2008; Kazemi et al. 2011; Tasse 2014a; Smirnov & Tasse 2015) or DD-I-RIME (direction-dependent imaging, see Bhatnagar et al. 2008; Tasse et al. 2010, 2018), in this article we aim at developing a robust approach using existing DD-C-RIME and DD-I-RIME algorithms to tackle the complete RIME inversion problem.

## 2.2. Direction-independent calibration

The standard LOTSS HBA observations consist of a ten-minute scan on a bright calibrator source (in general 3C 196 or 3C 295) before observing the target field for eight hours. On both calibrator and target fields, the visibilities of the 240 subbands (SBs) are regularly distributed across the 120–168 MHz bandpass, with 64 channels per 195.3 kHz SB and 1 sec integration time. The data are first flagged using AOFLAGGER<fnref n="3" /> (Offringa et al. 2012) and are averaged to 16ch/sb and 1s.

The interferometric data taken on the calibrator field are used to estimate the direction-independent Jones matrices $G$ that are, to the first order, the same in the target and calibrator fields. These include (i) the individual LOFAR station electronics and (ii) the clock offsets and drifts. This first pass of calibration is conducted using the PREFACTOR software package<fnref n="4" /> (de Gasperin et al. 2019). Specifically, as the calibrator field essentially consists of a single bright source, the measurement equation is direction-independent and the visibilities are modeled as

$$\mathbf{v}_b^{\text{cal}} = \mathcal{G}_b^0 \mathbf{v}_b^{\text{model}}, \quad (5)$$

where $\mathbf{v}_b^{\text{model}} = \int_s x_{sv} k_b^s ds$ is the sky model of the calibrator. We cannot just use $\mathcal{G}_b^0$ to calibrate the target field since the ionosphere is different for the calibrator and target fields. Instead, we want to extract (i) the bandpass and (ii) the clock offsets from the calibrator field, these being valid for the target field. The effective Mueller matrix of a given baseline $\mathcal{G}_b^0$ can be decomposed as the product of a direction-independent and direction-dependent term

$$\mathcal{G}_b^0 = \mathcal{G}_b^0 \mathcal{J}_b^0, \quad (6)$$

$$\text{with } \mathcal{G}_{pv} = \mathcal{A}_{pv} \exp(i\nu\Delta_p^t) \quad (7)$$

$$\text{and } \mathbf{J}_{pv}^0 = \exp(iK\nu^{-1}\Delta_p^T), \quad (8)$$

where $K = 8.44 \times 10^9 \text{ m}^3 \text{ s}^{-2}$ and $A_{pv}$, $\Delta_p^t$, and $\Delta_p^T$ are real-valued and represent, respectively, the bandpass, the clock, and the ionospheric total electron content (TEC) offset with respect to a reference antenna. The terms $\Delta_p^t$ and $\Delta_p^T$ can be disentangled from the frequency-dependent phases because the former is linear with $v$ while the latter is nonlinear.

Assuming the clocks and bandpass are the same for the calibrator and for the target field, the corrected visibilities $\mathbf{v}_b^c$ for the target field can be built from the raw data $v_b$ as $\mathbf{v}_b^c = \mathcal{G}_b^{-1} \mathbf{v}_b$. In order to calibrate for the remaining phase errors, the target field

<footnote n="1">¹ As described by Hamaker et al. (1996), the Mueller matrices and the Jones matrices can be related to each other using the Vec operator. In the context of the RIME, if $J_p$ and $J_q$ are $2 \times 2$ Jones matrices of antenna $p$ and $q$, and $X$ is the $2 \times 2$ source’s coherency matrix, then we have $\text{Vec}(J_p X J_q) = (J_q \otimes J_p) \text{Vec}(X)$, where $\otimes$ is the Kronecker product.</footnote>

<footnote n="2">² https://github.com/lofar-astron/LOFARBeam</footnote>

<footnote n="3">³ https://sourceforge.net/p/aoflagger/wiki/Home</footnote>

<footnote n="4">⁴ https://github.com/lofar-astron/prefactor</footnote>
is DI calibrated against the TGSS catalogs (Intema et al. 2017) and visibilities are averaged to 2 ch/sb and 8s.

## 2.3. Direction-dependent calibration and imaging

As discussed by Tasse (2014b), there are two families of calibration algorithms. “Physics-based” solvers aim at estimating the underlying Jones matrices, whose product gives the effective $G_{ptv}^s$ and $J_{ptv}^s$. Depending on the observing frequency and instrumental setup, these can be the product of the ionospheric Faraday rotation matrix, the scalar phase shift, or the individual station primary beams. This approach has the great advantage of constraining the free parameters to a low number, but it requires one to analytically model the physics of the various effects over the {$sptv$} space to be able to disentangle them. The second family of algorithms directly estimates the effective Jones matrices over piecewise constant domains in {$sptv$} space. These “Jones-based” solvers have the advantage of not requiring any physical modeling of the DDEs. However, this makes the number of degrees of freedom increase dramatically, typically by a few orders of magnitude. These additional degrees of freedom can often make the inverse problem ill-posed<fnref n="5" />. This means, in practice, that the DD solvers can overfit the data, leading to the unmodeled sky flux being absorbed by the calibration solutions. This happens differently at different scales and has a greater effect on the extended emission, which is measured only by the less numerous shorter baselines. Similar to linear problems, the situation depends on the sizes of the parameter space, as well as on the shape of the neighboring domains in the {$sptv$} spaces. Furthermore, as explained by Shimwell et al. (2019), experience shows that we need to split the sky model into a few tens of directions (“facets”) to be able to properly describe the spatial variation of the DD Jones matrices. This effect is amplified by the difficulty of properly modeling the extended emission itself. Indeed, even in the absence of calibration errors, the deconvolution problem consisting of inverting Eq. (1) by estimating $x_{sv}$ for given $G$, $J$, and $v$ is ill-posed. Furthermore, the situation is more severe when the point spread function (PSF) is less point-like (i.e., when the uv plane is not well sampled). While the true measured visibilities are described by Eq. (1), the (“model”) visibilities $\tilde{v}_b$ that are estimated<fnref n="6" /> over the piecewise constant domains $p, \Omega_\varphi, \Delta t$, and $\Delta v$ can be written as

$$ \tilde{v}_b = \tilde{G}_b \tilde{v}_b^\Sigma, \quad (9) $$

$$ \text{with } \tilde{v}_b^\Sigma = \sum_\varphi \tilde{v}_b^\varphi \quad (10) $$

$$ \text{and } \tilde{v}_b^\varphi = \tilde{J}_b^\varphi \mathcal{B}_b^\varphi \int_{s \in \Omega_\varphi} \tilde{x}_{sv} k_b^s \, ds, \quad (11) $$

where $\Omega_\varphi$ is the set of directions $s$ for a facet $\varphi$, $\tilde{x}_{sv}$ is the estimated underlying sky, and $\tilde{G}_b$ and $\tilde{J}_b^\varphi$ are the DI and DD Mueller matrices for baseline $b$, built from the corresponding estimated Jones matrices in $p, \Delta t$, and $\Delta v$.

Specifically, in order to solve for the DDEs, the size and shape of the domains are critical. Intuitively, if the domains are too small, not enough data points are used and the solutions are subject to ill-conditioning. On the other hand, if the domains are too large, the true Jones matrices can vary within the domain and the

piecewise constant approximation cannot account for the physics that underlies the measurement. Due to (i) the nonlinear nature of Eq. (1) and (ii) the complexity of the background radio sky, optimizing over the shape of these piecewise constant domains is a difficult problem (and is non-differentiable to some extent).

The faceted Jones-based approach is to find sky $\tilde{x}_{sv}$ as well as the DI $\tilde{G}_{ptv}$ and the DD piecewise constant $\tilde{J}_{ptv}^\varphi$ for all {$s\varphi ptv$} such that $\tilde{v}_b \sim v_b$ (we remain intentionally vague here since the cost function that is minimized depends on the specific DD algorithm). In practice, inverting Eq. (1) (estimating the Jones matrices and sky terms) is done by (i) solving for the Jones matrices assuming the sky is known (the calibration step) and (ii) solving for the sky-term assuming the Jones matrices are given (imag-ing step). Using this sky model and repeating steps (i) and (ii) is known as self-calibration, but in a third-generation approach we must explicitly model the direction-dependent aspects.

Since the computing time evolves as $\sim n_d^3$, where $n_d$ is the number of directions in the DD solvers, the problem of DD calibration has in general been tackled using direction alternating peeling-like techniques. Major breakthroughs have been made in the field of DD-C-RIME solvers over the past decade by Yatawatta et al. (2008); Kazemi et al. (2011), making this DD calibration computationally affordable. In addition, Tasse (2014a) and Smirnov & Tasse (2015) have described an alternative way to write the Jacobian of the cost function by using the Wirtinger differentiation method. The Jacobian and Hessian harbor different structures, and shortcuts can be taken to invert the calibration problem. The net gain over the classical method is not trivial: It can be as high as $n_d^2$ (Smirnov & Tasse 2015), where $n_d$ is the number of antennas in the interferometer. This Jones-based approach is therefore fast but is still subject to the same flaws as any other Jones-based solver.

Only a very few CI-RIME algorithms using a full DD self-calibration loop have been described and implemented. They include pointing self-calibration (Bhatnagar & Cornwell 2017) and peeling-based techniques such as MFIMAGE (implemented in the OBIT<fnref n="7" /> package) and FACTOR (van Weeren et al. 2016, see also Sect. 4.3). Similar to peeling, and developed for reducing LOFAR data, FACTOR is sequential along the direction axis. Looping over the different facets, it consists of (i) subtracting all sources, with the exception of calibration sources, in that one facet and (ii) DI self-calibrating in that direction. In addition to the ill-conditioning issues discussed above regarding DD-C-RIME and DD-I-RIME solvers, an expensive computational problem arises when estimating the $\tilde{J}_{ptv}^\varphi$.

The approach presented by Shimwell et al. (2019; also described in detail in Appendix A and referred to as DDF-PIPELINE-V1 in the following) is based on the kMS DD-C-RIME solver (Tasse 2014a; Smirnov & Tasse 2015) and DDFacet DD-I-RIME imager (Tasse et al. 2018) and is algebraically simultaneous in directions. The direction-dependent pipeline DDF-PIPELINE<fnref n="8" /> is a high level wrapper that mainly calls DDFacet<fnref n="9" /> and kMS<fnref n="10" /> for DD self-calibration. This type of algorithm has a number of advantages. Specifically, the interaction terms between the different directions are properly taken into account within the DD-C-RIME solver, that is to say, the DD affected sidelobes leaking from any facet to any other facet are accounted for within the algebraic operations of the algorithm. Another advantage over the FACTOR approach is that the data only need to be read, as

<footnote n="5">⁵ Linear and nonlinear problems can be ill-posed, meaning that nonunique solutions can be found.</footnote>

<footnote n="6">⁶ Throughout this paper, the notation $\tilde{x}$ should be read as “the estimate of $x$.”</footnote>

<footnote n="7">⁷ https://www.cv.nrao.edu/~bcotton/Obit.html</footnote>

<footnote n="8">⁸ https://github.com/mhardcastle/ddf-pipeline</footnote>

<footnote n="9">⁹ https://github.com/saopicc/DDFacet</footnote>

<footnote n="10">¹⁰ https://github.com/saopicc/killMS</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 2.** Some of the images produced at different steps in the DD self-calibration loop implemented as Algorithm 1. The maps correspond (from left to right, top to bottom) to Steps 1.1, 1.8, 1.12, and 1.18, respectively. The white lines show the facet locations. The color scale is the same on all panels and is displayed using an inverse hyperbolic sine function to render both the low level artifacts and some bright sources.</figcaption></figure>

opposed to modified, making the DDF-PIPELINE more efficient in terms of memory and disk access.

As explained in Shimwell et al. (2019), the DDF-PIPELINE-V1 is, however, affected by a number of issues. These include (i) artifacts and dynamic range limits around the brightest radio sources, (ii) artificial and diffuse haloes around moderately bright radio sources, and (iii) unmodeled flux absorption.

### 3. Calibration and imaging robustness

In this section, we describe in detail a DD calibration and imaging algorithm that aims to make the overall RIME imaging and calibration solver more robust against artifacts around the brightest sources (Sect. 3.1) and unmodeled flux absorption (Sects. 3.2 and 3.3). An overview of this approach is shown in Algorithm 1

(the implementation of which is referred to as DDF-PIPELINE-V2), and the corresponding DD self-calibration loop is displayed in Fig. 2. The DDF-PIPELINE-V2 products are described in Sect. 3.4. In Sect. 3.5, we show that DDF-PIPELINE-V2 produces improved images as compared to those previously made with DDF-PIPELINE-V1 (see Appendix A and Shimwell et al. 2019). In Sect. 3.6, we discuss the DDF-PIPELINE-V2 computing efficiency.

#### 3.1. *Dynamic range issue*

With LOFAR's very large field of view, it is common to observe bright sources within the station's primary beam. Tests we conducted on fields containing bright sources (such as 3C 295, which has a flux density of ~100 Jy at 150 MHz) have shown
**Algorithm 1:** Overview of the algorithm implemented in DDF-PIPELINE-V2. The function $\tilde{I}$ represents the imaging step and takes as input the visibility vector **v** together with the beam model $\mathbf{B}_{\Omega_n}$ and kMS-estimated Jones matrices $\mathbf{J}_{\Omega_n}$ at locations $\Omega_n$. The function $\mathcal{K}$ abstracts the DD calibration step, and takes as arguments the visibilities **v**, the skymodel $\tilde{\mathbf{x}}_v$, a solver mode (estimating for either scalar or full Jones matrices), a time-frequency solution interval (in min and MHz), and a set of directions $\Omega_n$ in which to solve for. The extra functions $\mathcal{C}$, $\mathcal{B}$, and $\mathcal{F}$ represent the clustering, bootstrapping and smoothing steps respectively.

**Data:** Visibilities **v** calibrated from DI effects using PREFACTOR.

/* On 60 LOFAR HBA subbands */

/* DI initial deconv and clustering */

1.1 $\tilde{\mathbf{x}}_v \leftarrow \mathcal{I}(\mathbf{v}_6, \mathbf{J}_{\Omega_r} = \mathbf{1}, \mathbf{B}_{\Omega_r});$

1.2 $\Omega_n \leftarrow C(\tilde{\mathbf{x}}_v);$ /* DI calibration and imaging */

1.3 $\tilde{\mathbf{x}}_v \leftarrow \mathcal{I}(\mathbf{v}_6, \mathbf{B}_{\Omega_n});$

1.4 $\mathbf{v}_6^c \leftarrow \mathcal{K}(\mathbf{v}_6, \tilde{\mathbf{x}}_v, \mathbf{B}_{\Omega_n}|\text{full}, \delta t_0, \delta v_0, \Omega_0);$

1.5 $\tilde{\mathbf{x}}_v \leftarrow \mathcal{I}(\mathbf{v}_6^c, \mathbf{B}_{\Omega_n});$ /* Bootstrapping the flux density scale */

1.6 $\mathbf{v}^c \leftarrow \mathcal{B}(\mathbf{v}_6^c);$ /* Phase only DD calibration and imaging */

1.7 $\tilde{\mathbf{J}} \leftarrow \varphi \circ \mathcal{F} \circ \mathcal{K}(\mathbf{v}_6^c, \tilde{\mathbf{x}}_v, \mathbf{B}_{\Omega_n}|\text{scalar}, 1\text{min}, 2\text{MHz}, \Omega_n);$

1.8 $\tilde{\mathbf{x}}_v \leftarrow \mathcal{I}(\mathbf{v}_6^c, \tilde{\mathbf{J}}\mathbf{B}_{\Omega_n});$ /* DD calibration and imaging */

1.9 $\tilde{\mathbf{J}} \leftarrow \mathcal{F} \circ \mathcal{K}(\mathbf{v}_6^c, \mathbf{B}_{\Omega_n}, \tilde{\mathbf{x}}_v|\text{scalar}, 1\text{min}, 2\text{MHz}, \Omega_n);$

1.10 $\tilde{\mathbf{x}}_v \leftarrow \mathcal{I}(\mathbf{v}_6, \tilde{\mathbf{J}}\mathbf{B}_{\Omega_n});$ /* DI calibration and imaging */

1.11 $\mathbf{v}_6^c \leftarrow \mathcal{K}(\mathbf{v}_6, \tilde{\mathbf{J}}\mathbf{B}_{\Omega_n}, \tilde{\mathbf{x}}_v|\text{full}, \delta t_0, \delta v_0, \Omega_0);$

1.12 $\tilde{\mathbf{x}}_v \leftarrow \mathcal{I}(\mathbf{v}_6^c, \tilde{\mathbf{J}}\mathbf{B}_{\Omega_n});$

/* On 240 LOFAR HBA subbands */

/* DD calibration */

1.13 $\tilde{\mathbf{J}} \leftarrow \mathcal{F} \circ \mathcal{K}(\mathbf{v}_{24}, \mathbf{B}_{\Omega_n}, \tilde{\mathbf{x}}_v|\text{scalar}, 1\text{min}, 2\text{MHz}, \Omega_n);$

/* DI calibration */

1.14 $\mathbf{v}_{24}^c \leftarrow \mathcal{K}(\mathbf{v}_{24}, \tilde{\mathbf{J}}\mathbf{B}_{\Omega_n}, \tilde{\mathbf{x}}_v|\text{full}, \delta t_0, \delta v_0, \Omega_0);$

/* DD imaging */

1.15 $\tilde{\mathbf{x}}_v \leftarrow \mathcal{I}(\mathbf{v}_{24}^c, \tilde{\mathbf{J}}\mathbf{B}_{\Omega_n});$

/* DD calibration */

1.16 $\tilde{\mathbf{J}} \leftarrow \mathcal{F} \circ \mathcal{K}(\mathbf{v}_{24}^c, \mathbf{B}_{\Omega_n}, \tilde{\mathbf{x}}_v|\text{scalar}, 1\text{min}, 2\text{MHz}, \Omega_n);$

/* Slow DD calibration */

1.17 $\tilde{\mathbf{J}}_s \leftarrow \mathcal{K}(\mathbf{v}_{24}^c, \tilde{\mathbf{J}}\mathbf{B}_{\Omega_n}, \tilde{\mathbf{x}}_v|\text{scalar}, 43\text{min}, 2\text{MHz}, \Omega_n);$

/* Final imaging steps */

1.18 $\tilde{\mathbf{x}}_v \leftarrow \mathcal{I}(\mathbf{v}_{24}^c, \tilde{\mathbf{J}}_s\tilde{\mathbf{J}}\mathbf{B}_{\Omega_n});$

1.19 Facet-based astrometric correction (see Shimwell et al. 2019, for details);

that the related residual errors create powerful artifacts that largely dominate the thermal noise, thereby driving a dynamic range limit. As explained in Sect. 2.2, the initial phase calibration is done against TGSS at 150 MHz (Intema et al. 2017).

<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Amplitude of the diagonal (black) and off-diagonal (gray) terms of the estimated Jones matrices for a remote station using the LOFAR observation synthesized 6" image as the sky model (self-calibration). If the initial DI calibration and correction by PREFACTOR on the lower resolution TGSS sky model were good enough, the calibration solutions found here would be the unity matrix at all times and frequencies. Therefore, in DDF-PIPELINE-V2, we carry out a few full-Jones DI-only self-calibration steps. This approach seems very efficient in increasing the dynamic range around the brightest radio sources.</figcaption></figure>

However, since LOTSS resolution is much higher than that of TGSS (6" × 6" vs. 25" × 25", respectively), small spatial uncertainties on how the individual bright sources are modeled lead to large Jones matrix errors. Specifically, this effect can be severe when the true point sources are erroneously found to be resolved by TGSS as this leads to large calibration errors for the long baselines. In these situations, our experience shows that the initial DI calibration is not good enough to start doing a DD calibration (due to ill-conditioning issues, it has to be done on larger time-frequency solution intervals). In the following, we study the DI calibration solutions and assess whether they need to be recomputed using a high angular resolution sky model. When using the PREFACTOR DI-calibrated LOTSS data and associated imaging products, this amounts to doing a round of DI self-calibration at the beginning of DDF-PIPELINE-V2.

Solution time and frequency variability is, however, hard to interpret. Indeed, because the RIME formalism is subject to unitary ambiguity (see Hamaker 2000, for a detailed discussion), the off-diagonal or absolute phase terms found by a solver are not meaningful. Instead, they are given with respect to a reference antenna. When Jones matrices are scalar, this amounts to zeroing the phases $\varphi_0$ of the reference antenna by subtracting $\varphi_0$ from all phases of all antennas. To do this in the general case of non-diagonal Jones matrices, we used a polar decomposition on the Jones matrix $\mathbf{J}_0$ of the reference antenna such that $\mathbf{J}_0 = \mathbf{U}\mathbf{P}_0$, where $\mathbf{U}$ is a unitary matrix$^{11}$. We then applied $\mathbf{U}$ to all Jones matrices as $\mathbf{J}_p \leftarrow \mathbf{U}^H\mathbf{J}_p$. Intuitively, when the Jones matrices are all scalar, the unitary matrix $\mathbf{U}$ is simply $\exp(i\varphi_0)\mathbf{I}$, and that step makes the phases of all $\mathbf{J}_p$ relative to the reference antenna (and specifically zeros the phases of $\mathbf{J}_0$). In the case of nontrivial $2 \times 2$ Jones matrices, finding and applying $\mathbf{U}$ has the effect of removing a common rotation from all Jones matrices, and it orthogonalizes them.

We apply this in Fig. 3, where we show the typical DI Jones matrices we can estimate at Step 1.3 for a given remote station

<footnote>$^{11}$ The unitary matrix $\mathbf{U}$ is found by doing a singular value decomposition $\mathbf{J}_0 = \mathbf{W}\Sigma\mathbf{V}$ and is then built as $\mathbf{U} = \mathbf{W}\mathbf{V}^H$.</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 4. Simulations conducted to test the robustness of the algorithm described in Sect. 3 and implemented in DDF-PIPELINE-V2. We have simulated an unmodeled extended emission (left-hand panel). The emission is absorbed by the DD-calibration step (central panel), while it can be partially recovered (right-hand panel) by decreasing the effective size of the unknown solutions space (Sects. 3.2 and 3.3).</figcaption></figure>

and frequency, with respect to a reference station in the LOFAR core. They were estimated using kMS and the sky model synthesized by DDFacet from the visibilities corrected by PREFACTOR. Since the polar transform has been applied, the variations of the amplitude of the off-diagonal Jones matrices are genuine. These are interpretable in terms of differential Faraday rotation: The rotation of the electric field polarization changes across the LOFAR array. This demonstrates the need to conduct a full-Jones calibration on the PREFACTOR-calibrated LOTSS data.

Therefore, in Step 1.4 the visibilities are calibrated against the modeled visibilities generated by DDFacet in Step 1.3. As the sky is mostly unpolarized, we assumed $Q = U = V = 0$ Jy in this full-Jones DI calibration step (see Sect. 3.4 for a discussion of polarization related data products). The solution intervals $\delta t_0$ and $\delta v_0$ along time and frequency are determined such that $n_b \propto (T/\langle|\mathbf{x}_v|\rangle)^2 \text{Var}\{\mathbf{n}\}$, where $n_b$ is the number of points in the $\delta t_0 \times \delta v_0$ time-frequency domain, $T$ is the target solution S/N, and $\text{Var}\{\mathbf{n}\}$ is the variance of the visibilities' noise (see Sob et al. 2020 for a justification).

We note that after the initial DD calibration solutions are obtained in Steps 1.7 and 1.9, a more accurate DI calibration can be performed. Specifically, in the DI calibration Steps 1.11 and 1.14, the model visibilities $\tilde{\mathbf{v}}_b$ (Eq. (10)) are predicted on any baseline $b$ based on the previously estimated DD-Jones matrices $\tilde{\mathbf{J}}$ (as is done by Smirnov 2011).

## 3.2. Regularization

The absorption of unmodeled flux by calibration is a well-known issue connected to the calibration of DDES. Intuitively speaking, when real flux is missing from the modeled sky $\tilde{\mathbf{x}}_i$ of $\mathbf{x}$ at step $i$, and since the RIME inversion is often ill-posed, the estimates $\tilde{\mathbf{J}}_i$ of $\mathbf{J}$ can be biased in a systematic way. Experience and simulations show that building a new estimate $\tilde{\mathbf{x}}_{i+1}$ from $\tilde{\mathbf{J}}_i$ can be biased in that the unmodeled emission is not and never will be recovered (Fig. 4). This effect is especially severe when the extended emission is poorly modeled or unmodeled since it is detected only by the shortest baselines. Effectively, during the inversion of the RIME system of equations, the DD self-calibration algorithm falls into the wrong (local) minimum.

One way to address this problem is to reduce the effective number of free parameters used to describe the Jones matrices in

the ${\text{pdtv}}$-space (see for example Tasse 2014a; Yatawatta 2015; van Weeren et al. 2016; Repetti et al. 2017; Birdi et al. 2020). Forcing the estimated Jones matrix shape to look like that of the real underlying ones improves the conditioning of the inverse problem. In Algorithm 1 (implemented in DDF-PIPELINE-V2), we replaced that normalization method with a smoothing of the kMS-estimated Jones matrices. This function $\mathcal{F}$ updates the Jones matrices $\mathbf{J} \leftarrow \mathcal{F}(\mathbf{J})$ in each direction independently by imposing on them a certain behavior in the time-frequency space (see below), effectively reducing the size of the unknown stochastic process. This can be thought of as a regularization. This is done independently on the phases and amplitudes on the scalar Jones matrices generated in Steps 1.9, 1.13, 1.16. The updated Jones matrices take the analytical form

$$ \tilde{\mathbf{J}}_{pd,tv} = \tilde{a}_{pd,t} \mathcal{P}(t, \tilde{\theta}_{pd,v}) \exp(iKv^{-1} \Delta_p^T) \mathbf{I}, \quad (12) $$

where $\tilde{\Delta}_p^T$ is the differential TEC (see also Sect. 2.2 and Eq. (8)), $\mathcal{P}$ is a polynomial parametrized by the coefficients in $\theta_{pd,v}$ (of size 10), and $\tilde{a}_{pd,t}$ is a scalar meant to describe the loss of correlation due to ionospheric scintillation, as seen in the top left-hand panel of Fig. 5. Typically, for the ~8-h integration of LOTSS pointings and solving every 30 s and 2 MHz, this parametrization of the Jones matrices reduces the number of free parameters by a factor of $\ge 20$.

To assess the recovery of unmodeled flux in DDF-PIPELINE-V2, a series of simulations were conducted in which faint simulated sources of various fluxes and extents were injected into real LOFAR data that had been fully processed with the DDF-PIPELINE-V2 strategy. The properties of the injected sources were chosen to be typical for large extragalactic objects such as radio haloes of galaxy clusters. After the injection of the artificial extended sources, Steps 1.16 and 1.17 were repeated using the sky model derived in Step 1.18 prior to the injection of the sources. These simulations will be discussed further by Shimwell et al. (in prep.), but briefly, in each simulation the recovered flux of the completely unmodeled emission exceeded 60%. Examples of the injected and recovered emission are shown in Fig. 4.

As suggested by the results of the simulations, the effect on real data is in general very satisfactory and allows us to recover a good fraction of the unmodeled extended emission even when it is quite faint and extended. This is shown in Fig. 6 for a typical LOTSS observation. Here, the extended emission is about 10'
<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Amplitude and phase (top and bottom, respectively) of a scalar Jones matrix for a given station in a given direction in the example observation. **Left-hand panel:** solution as estimated by the kMS solver. **Right-hand panel:** regularized solution as updated by the *F* function. The amplitude color scale ranges from 0 to 1.5.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 6.** Effects of two calibration steps on imaging artifacts and flux absorption. **Left-hand panel:** shows that a faint and unmodeled extended emission (on the level of ~0.7σ here) can be totally absorbed. While regularizing the DD calibration solutions can help in recovering the unmodeled emission (typically after Step 1.16), it can also produce negative imaging artifacts and “holes” around bright sources (**central panel**). **Right-hand panel:** shows that solving the residuals on longer time intervals (Step 1.17) corrects for this issue.</figcaption></figure>

across, with a mean flux density at the peak of only 70% of the
local standard deviation.

**3.3. Conditioning and solution interval**

The additional issue of arcminute-scale negative haloes appear-
ing around bright compact sources (at a level of ≤1% or the
peak) could, however, be seen in ~10–20% of the LOTSS point-
ings processed with DDF-PIPELINE-V1. As shown in Fig. 6, we
believe this to be connected to the solution regularization itself.
This issue is hard to understand in detail because of the non-
linearity in the C-RIME inversion, but it is likely due to the
pointings that show these issues being more severely affected
by the incompleteness of the sky model. Specifically, by con-
ducting several experiments we were able to observe that the
situation was improved by deconvolving more deeply or by tak-
ing into account sources outside the synthesized image field of
view.

An additional way to improve the conditioning of the prob-
lem is to increase the amount of data used to constrain the
estimated Jones matrices. For the DD calibration steps pre-
sented in Algorithm 1, we used solution intervals of 0.5–1 min.

Following van Weeren et al. (2016), we added an extra calibra-
tion, Step 1.17, where the visibilities are modeled using the latest
available sky model together with the smoothed Jones matrices
estimated in Step 1.16 and which are defined over the finer time
and frequency mesh. Intuitively, since the negative haloes are
produced by some systematic effects, the idea is to calibrate for
a slowly varying differential effect. The time interval is set to
~43 min in DDF-PIPELINE-V2, giving 11 solution intervals in
the eight-hour LOTSS pointings. This interval has to be long
enough to reach a good conditioning for the C-RIME inversion
and short enough to sample the Jones matrices’ remaining phys-
ical variations. As shown in the right-hand panel of Fig. 6, this
method is very efficient in reducing the negative haloes.

3.4. Data products

3.4.1. Unpolarized flux

Once the estimated DD-Jones matrices and sky model $\hat{\mathbf{x}}_v$
have been obtained at the highest available spatial resolution
following the DI and DD self-calibration steps presented in
Algorithm 1, additional data products are formed. Users can
adapt the weighting scheme depending on the scientific exploitation of the interferometric data they want to make use of. This is very much tied to how the calibration and deconvolution algorithms work, and concurrent effects take place along the self calibration loop. Extended emission is hard to properly model since the deconvolution problem is more ill-posed in these cases (more pixels are nonzero). To tackle this issue, the PSF can be modified to make the convolution matrix more diagonal and the deconvolution problem correspondingly better conditioned. This is done at the cost of a lower sensitivity, which in turn can drive systematic errors in the calibration solution estimates because extended emission is poorly modeled on the shorter baselines.

For all these concurring reasons, the faint and extended flux in the highest resolution maps produced by Algorithm 1 is either poorly modeled or not deconvolved at all. Since the pixel values of extended sources are not interpretable in the residual maps, the flux density of the radio sources cannot be measured if they are not deconvolved. We therefore intentionally degraded the resolution of some of the imaging to allow survey users to choose a resolution based on the broad scientific topic that they need to address. Additionally, we stored the sub-space deconvolution (SSD, see Tasse et al. 2018) masks as residual images so the end user can know if any given source has been deconvolved. With this in mind, the following Stokes *I* products are generated: (i) a high resolution (6", 20.000 × 20.000 pixels) wide-bandwidth Stokes *I* image (Step 1.18), (ii) a low resolution (20") wide-bandwidth Stokes *I* image (Step 2b.1), and (iii) a high resolution (6") Stokes *I* image in three frequency chunks spread over the whole HBA bandwidth (Step 2b.2).

The DI-calibrated visibilities as well as the final sky models and DD calibration solutions are stored. This allows for additional postprocessing to be made, such as better calibration toward a particular point on the sky (van Weeren et al. 2020), as well as for reimaging at different resolutions if required.

### 3.4.2. QUV images

The DDFacet DD-imager only deals with *I*-Stokes deconvolution. As discussed by Tasse et al. (2018), estimating the QUV Stokes parameters is complex in the context of DD imaging due to the leakage terms. Indeed, for the problem to be properly addressed, 16 PSFs have to be computed (as there are 16 terms in the quadratic mean of the Mueller matrices). As most of the sources are unpolarized, the leakage terms are properly taken into account in the DD predict (i.e., the forward mapping from sky and Jones matrix estimates to modeled visibilities). Instead of deconvolving the polarized flux, we grid the IQUV residual data. The polarized flux is directly interpretable when the sources are unresolved. Hence, we also generated the following additional products: (i) low resolution (20") spectral Stokes QU cubes (480 planes - Step 2b.3), (ii) very low resolution Stokes QU cubes (480 planes - Step 2b.4), by cutting the baselines >1.6 km, giving an effective resolution of ~3', and (iii) a low resolution (20") wide-bandwidth Stokes V image (Step 2b.5).

The output QU cubes are processed using Faraday rotation measure (RM) synthesis (Brentjens & de Bruyn 2005) to find polarized sources and their RMs with the sensitivity of the full bandwidth. The wide bandwidth (120 to 168 MHz) combined with the narrow channel width (97.6 kHz) provides a resolution in RM space of ~1.1 rad m⁻² and an ability to measure RMs of up to ~450 rad m⁻² (e.g. O'Sullivan et al. 2020).

The 3' QU cubes are sensitive to the large-scale polarized emission from the Milky Way, while the 20'' QU cubes are excellent for finding compact polarized sources. However, detailed

<figure><img src="image_7.png" /><figcaption>**Fig. 7.** Faraday depth spectrum, or Faraday dispersion function (FDF), for a radio galaxy in both the DR1 and DR2 datasets, showing the improvement in the suppression of the instrumental polarization signal. The dashed blue line shows the FDF from the DR1 data with a strong instrumental polarization feature near Faraday depths of $\phi \sim 0$ rad m⁻², while the solid orange line shows the FDF from the DR2 data in which the instrumental feature is suppressed below the noise level. In both cases, the Faraday depth of the real astrophysical signal is the same.</figcaption></figure>

studies of the polarization and RM structure of resolved extra-galactic sources will require the deconvolution of the Q and U data. The DDF-PIPELINE-V2 output provides significantly better performance when correcting for the effect of the instrumental polarization (Fig. 7), which is typically at the level of 1% or less for bright total intensity sources (O'Sullivan et al., in prep.).

There is no absolute polarization angle calibration for each LoTSS observation, meaning that while the RM values of sources in overlapping fields are consistent, the polarization angles are not. Therefore, to avoid unnecessary depolarization for both mosaicing and the deep fields, the polarization angles between the observations need to be aligned. The simplest way to do this is by choosing a reference angle of a polarized source in a single observation and applying a polarization angle correction to all other observations to align them with this reference angle, as presented in Herrera Ruiz et al. (2020). An alternative approach is to use the diffuse polarized emission that is present in the ~3' QU cubes.

Bright polarized sources are rare in the LoTSS data, with only three sources having a polarized intensity greater than 50 mJy beam⁻¹ in the DR1 HETDEX sky area (Van Eck et al. 2018; O'Sullivan et al. 2018). However, in the fields containing these bright polarized sources, the DDF-PIPELINE-V2 output becomes unreliable for polarized sources. This limitation likely arises from assuming $Q = U = V = 0$ Jy for a field in the DI calibration step. While only a few percent of fields are strongly affected, the exact extent of this issue is being investigated further through simulations, where bright polarized sources are inserted into existing LoTSS uv-datasets. Possible solutions will be tested in future pipeline developments.

### 3.5. Comparison between DDF-PIPELINE-V1 and DDF-PIPELINE-V2

Figure 8 shows the comparison between the final high resolution images produced by DDF-PIPELINE-V1 and DDF-PIPELINE-V2 for an eight-hour integration LoTSS pointing (P26HETDEX03). Many processes are involved in the sky reconstruction from radio interferometric data. Imaging and calibration affect the final synthesized maps and introduce complex and systematic residual
<figure><img src="image_8.png" /><figcaption>**Fig. 8.** Differences between the maps produced by Algorithms 0 and 1 from typical eight-hour scans (here the P26HETDEX03 pointing in the Hobby-Eberly Telescope Dark Energy Experiment, HETDEX, field, see Shimwell et al. 2017). The color scale is the same on all panels and is displayed using an inverse hyperbolic sine function to render both the low level artifacts and some bright sources. (a) The central part of the P26Hetdex03 8 hours LOFAR-HBA scan as imaged by Alg. 0. (b) Region (1) as imaged by Alg. 0. (c) Region (1) as imaged by Algorithm 1. (d) The central part of the P26Hetdex03 eight-hour LOFAR-HBA scan as imaged by Algorithm 1. (e) Region (2) as imaged by Algorithm 0. (f) Region (2) as imaged by Algorithm 1.</figcaption></figure>
<figure><img src="image_9.png" /></figure>

**Fig. 9.** Nature and ordering of the different steps of Algorithm 1 and how the computing time is distributed across them. The light and dark gray areas represent the imaging and calibration steps, respectively. The black area represents the miscellaneous tasks (additional data products, see Sect. 3.4) that are done once the DI and DD self-calibration loops have completed. These graphs were created from a DDF-PIPELINE-V2 run on a node equipped with 192 GBytes RAM and 2 Intel Xeon Gold 6130 CPU at 2.10 GHz, giving 32 physical compute cores. The dashed area is a quadrant representing a day, while the inner pie shows the total contributions of the imaging, calibration, and miscellaneous tasks.

errors. It is therefore difficult to find a good and absolute metric to compare the final imaging products.

As discussed in Sect. 3.1, the quality of the initial DI cal-
ibration proved to be quite crucial for the feasibility of the
following DD calibration and imaging steps. DDF-PIPELINE-V1
was indeed failing at imaging certain fields with very bright
sources, while artifacts were present around most moderately
bright ones thereby driving the dynamic range limit in large
areas. In Figs. 8b and 8c, we show a radio source imaged by
DDF-PIPELINE-V1 and DDF-PIPELINE-V2.

Another important issue with the approach we presented in Shimwell et al. (2019) was the presence of a low spatial frequency pattern corresponding to a positive or negative halo around radio sources. Although the effect is complex to analyze, we concluded from various experiments that these systematics were due to the combination of (i) sky model incompleteness, (ii) a uv-distance cut used during the calibration, and (iii) the N normalization function (see Appendix A for details) that we had introduced for DDF-PIPELINE-V1 to be robust against the absorption of extended emission. As shown in Figs. 8e and 8f, the approach developed in Sects. 3.2 and 3.3 to conserve unmodeled extended emission and implemented in DDF-PIPELINE-V2 does not produce any significant low spatial frequency systematics.

3.6. DDF-PIPELINE-V2 robustness and performance

As explained above, DDF-PIPELINE-V2 is a high level script
that interfaces kMS and DDFacet. Both of the underlying soft-
ware packages are efficiently parallelized using a custom version
of the Python multiprocessing package for process-level
parallelism and using the SharedArray<fnref n="12" /> module. As explained

<figure><img src="image_10.png" /></figure>

**Fig. 10.** Sensitivity of the various deep dedicated surveys covering the Boötes (top) and Lockman Hole (bottom) fields as a function of observing frequency. The resolution of the various surveys corresponds to the radius of the black dot, while the diameter of the corresponding surveyed area is encoded in the size of the gray circle. The LOTSS-Deep Fields pointings are marked with a red cross, and the dashed line corresponds to a source that has a spectral index of -0.7.

by Tasse et al. (2018), this pythonic approach minimizes the
process interconnections for both the kMS and DDFacet software.

This paper considers the application of DDF-PIPELINE-V2 to the LOTSS-Deep Fields. The pipeline is also being used to process data from the wider and shallower LOTSS survey. The LOTSS project is presently observing at a rate of up to 1500 h every six-month cycle, which corresponds to approximately two eight-hour pointings (observed simultaneously) each day. The DDF-PIPELINE-V2 compute time is roughly split equally between calibration and imaging tasks (see Fig. 9). The total run time for an eight-hour pointing is ~5 days (on a node equipped with 192 GBytes RAM and 2 Intel Xeon Gold 6130 CPU at 2.10 GHz, giving 32 physical compute cores) and takes an extra ~30% of computing time for completion as compared to DDF-PIPELINE-V1. Hence ten compute nodes are sufficient to keep up with the observing rate. However, in practice, more compute nodes are used because LOTSS has been observing since 2014, and, as of June 1, 2019, over 1000 pointings existed in the archive. Over ~1000 pointings and ~12 Petabytes of averaged and compressed LOFAR data (~40 Petabytes uncompressed) have now been processed with DDF-PIPELINE-V2.

4. LOTSS-Deep Fields data and processing

4.1. *Observations*

LOTSS-Deep Fields observations are being carried out over
the four northern fields with high-Galactic latitude and the

<footnote n="12">¹² https://pypi.python.org/pypi/SharedArray</footnote>
<table><caption>**Table 1.** Overview of the deep field pointings used to synthesize the images on the Boötes and Lockman Hole extragalactic fields.</caption><thead><tr><th>Project ID</th><th>LOFAR Obs. ID</th><th>Obs. Date</th><th>Start time<br>UTC</th><th>Integration<br>time [h]</th><th>ν<sub>min</sub><br>[MHz]</th><th>ν<sub>max</sub><br>[MHz]</th><th>f<sub>flag</sub></th><th>nMS</th></tr></thead><tbody><tr><td colspan="9" style="text-align:center; font-style:italic;">Boötes</td></tr><tr><td>LC2_038</td><td>L236786</td><td>21/07/2014</td><td>10:44:00</td><td>8.0</td><td>120.0</td><td>168.7</td><td>37.8</td><td>25</td></tr><tr><td>LC2_038</td><td>L243561</td><td>15/09/2014</td><td>10:22:42</td><td>8.0</td><td>120.0</td><td>168.7</td><td>19.2</td><td>25</td></tr><tr><td>LC4_034</td><td>L346004</td><td>11/06/2015</td><td>16:04:35</td><td>8.0</td><td>120.2</td><td>167.0</td><td>10.9</td><td>24</td></tr><tr><td>LC4_034</td><td>L373377</td><td>17/09/2015</td><td>10:21:18</td><td>8.0</td><td>120.2</td><td>168.9</td><td>20.1</td><td>25</td></tr><tr><td>LC4_034</td><td>L374583</td><td>24/09/2015</td><td>10:09:57</td><td>8.0</td><td>120.2</td><td>168.9</td><td>10.6</td><td>25</td></tr><tr><td>LC4_034</td><td>L387597</td><td>29/09/2015</td><td>09:13:00</td><td>8.0</td><td>120.2</td><td>168.9</td><td>27.4</td><td>25</td></tr><tr><td>LC4_034</td><td>L387569</td><td>01/10/2015</td><td>09:00:00</td><td>8.0</td><td>120.2</td><td>168.9</td><td>32.3</td><td>25</td></tr><tr><td>LC4_034</td><td>L400135</td><td>10/10/2015</td><td>08:46:22</td><td>8.0</td><td>120.2</td><td>168.9</td><td>13.6</td><td>25</td></tr><tr><td>LC4_034</td><td>L401825</td><td>21/10/2015</td><td>08:00:30</td><td>8.0</td><td>120.2</td><td>168.9</td><td>9.0</td><td>25</td></tr><tr><td>LC4_034</td><td>L401839</td><td>22/10/2015</td><td>07:55:23</td><td>8.0</td><td>120.2</td><td>168.9</td><td>8.0</td><td>25</td></tr><tr><td colspan="9" style="text-align:center; font-style:italic;">Lockman Hole</td></tr><tr><td>LC3_008</td><td>L274099</td><td>08/03/2015</td><td>20:11:00</td><td>8.7</td><td>120.2</td><td>168.9</td><td>12.4</td><td>25</td></tr><tr><td>LC3_008</td><td>L281008</td><td>14/03/2015</td><td>18:26:39</td><td>8.7</td><td>120.4</td><td>169.1</td><td>8.1</td><td>25</td></tr><tr><td>LC3_008</td><td>L294287</td><td>21/03/2015</td><td>19:11:00</td><td>8.7</td><td>120.2</td><td>168.9</td><td>18.7</td><td>25</td></tr><tr><td>LC3_008</td><td>L299961</td><td>24/03/2015</td><td>17:47:20</td><td>8.7</td><td>120.2</td><td>168.9</td><td>12.1</td><td>25</td></tr><tr><td>LC3_008</td><td>L340794</td><td>25/04/2015</td><td>17:08:00</td><td>8.7</td><td>120.2</td><td>168.9</td><td>14.5</td><td>25</td></tr><tr><td>LC3_008</td><td>L342938</td><td>08/05/2015</td><td>14:50:24</td><td>8.7</td><td>120.2</td><td>168.9</td><td>16.6</td><td>25</td></tr><tr><td>LT10_012</td><td>L659554</td><td>10/07/2018</td><td>11:11:00</td><td>8.0</td><td>120.2</td><td>168.9</td><td>9.9</td><td>25</td></tr><tr><td>LT10_012</td><td>L659948</td><td>12/07/2018</td><td>11:08:10</td><td>8.0</td><td>120.2</td><td>168.9</td><td>11.9</td><td>25</td></tr><tr><td>LT10_012</td><td>L664320</td><td>15/08/2018</td><td>08:49:00</td><td>8.0</td><td>120.2</td><td>168.9</td><td>11.0</td><td>25</td></tr><tr><td>LT10_012</td><td>L664480</td><td>19/08/2018</td><td>08:38:46</td><td>8.0</td><td>120.2</td><td>168.9</td><td>11.3</td><td>25</td></tr><tr><td>LT10_012</td><td>L667204</td><td>12/09/2018</td><td>07:06:09</td><td>8.0</td><td>120.2</td><td>168.9</td><td>11.2</td><td>25</td></tr><tr><td>LT10_012</td><td>L667218</td><td>13/09/2018</td><td>07:05:12</td><td>8.0</td><td>120.2</td><td>168.9</td><td>10.3</td><td>25</td></tr></tbody></table>

**Notes.** Columns *f*<sub>flag</sub> and nMS show the fraction of flagged data and the number of measurement sets present in the archives.

highest-quality multidegree-scale ancillary data across the electromagnetic spectrum: the Boötes, Lockman Hole, ELAIS-N1, and NEP fields. The ultimate aim of the LoTSS-Deep Fields project is to reach noise levels of 10–15 µJy beam⁻¹ in each of these fields (requiring ~500 h of integration). The first LoTSS-Deep Fields data release consists of initial observations in three of these fields: Boötes (~80 h) and Lockman Hole (~112 h), presented in the current paper, and ELAIS-N1 (presented by Sabater et al. 2021, for an integration time of ~170 h in Paper II). This first data release also includes an extensive effort of optical-to-IR cross-matching, which has obtained host galaxy identifications for over 97% of the ~80 000 radio sources detected within the ~25 deg² overlap with the high-quality multiwavelength data (Kondapally et al. 2021, Paper III). This is supplemented by high-quality photometric redshifts and a characterization of host galaxy properties (Duncan et al. 2021, Paper IV) and source classification (Best et al., in prep.). In order to put the LoTSS-Deep Fields observations into a wider context, in this section we briefly describe the multiwavelength data available on the Boötes and Lockman Hole fields, focusing on the radio coverage (for a more detailed description see Kondapally et al. 2021, Paper III).

### 4.1.1. Boötes field

The Boötes field is one of the NOAO Deep Wide Field Survey (NDWFS, Jannuzi & Dey 1999) fields covering ~9.2 deg². It contains multiwavelength data, including infrared (Spitzer space telescope, see Ashby et al. 2009; Jannuzi et al. 2010), X-rays (Chandra space telescope, see Murray et al. 2005; Kenter et al. 2005; Masini et al. 2020), and optical data (Jannuzi & Dey 1999;

Cool 2007; Brown et al. 2007, 2008). At radio frequencies, it has been mapped with the Westerbork Synthesis Radio Telescope (WSRT, see de Vries et al. 2002), the VLA (see Croft et al. 2008; Coppejans et al. 2015), the Giant Metrewave Radio Telescope (GMRT, see Intema et al. 2011; Williams et al. 2013), and LOFAR (van Weeren et al. 2014; Williams et al. 2016; Retana-Montenegro et al. 2018) at various depths, frequencies, resolutions, and covered areas (see Fig. 10 for an overview of the available radio data).

The Boötes pointing data that are presented in this paper are centered on (α, δ) = (14h32m00s, +34°30′00″) and were observed with the LOFAR-HBA in HBA_DUAL_inner mode during Cycles 2 and 4, with a bandwidth of 48 MHz (see Table 1). The total integration time of ~80 h is spread over ten scans of eight hours each.

### 4.1.2. Lockman Hole

The Lockman Hole field is also covered by a large variety of multiwavelength data. Specifically, it has been observed by SWIRE (Lonsdale et al. 2003) over ~11 deg², and over 16 deg² by the Herschel Multi-tiered Extragalactic Survey (Oliver et al. 2012). It has also been observed in UV (Martin & GALEX Team 2005), optical (González-Solares et al. 2011), near IR (UKIRT Infrared Deep Sky Survey - Deep Extragalactic Survey; UKIDSS-DXS, see Lawrence et al. 2007), and with the Submillimeter Common-User Bolometer Array (Coppin et al. 2006; Geach et al. 2017). At higher energy, it has been observed with XMM-Newton (Brunner et al. 2008), and Chandra (Polletta et al. 2006). In the radio domain, the Lockman Hole has been observed over the two deep aforementioned X-ray fields over small sub-deg² areas
<figure><img src="image_11.png" /><figcaption>**Fig. 11.** Restored high resolution image toward the center of the Boötes field for the eight-hour image produced with Algorithm 1 (top panel) and the 80-h image produced with Algorithm 2 (bottom panel). Both images are thermal noise limited, with the same color scale used on both.</figcaption></figure>

(de Ruiter et al. 1997; Ciliegi et al. 2003; Biggs & Ivison 2006; Ibar et al. 2009). Wide surveys of the Lockman Hole have been done with GMRT (Garn et al. 2010), VLA (Owen et al. 2009), WSRT (Guglielmino et al. 2012; Prandoni et al. 2018), and LOFAR at 150 MHz (Mahony et al. 2016). Figure 10 presents an overview of the available radio data on the Lockman Hole.

The Lockman Hole observation that we are presenting in this paper consists of 12 pointings of ~8 h centered on $(\alpha, \delta) = (10h47m00s, +58^\circ05'00'')$ and observed from March 2015 (Cycle 3) to November 2018 (Cycle 4). As for the Boötes field observation, we observed in HBA_DUAL INNER with ~48 MHz bandwidth, while the integration time depended on

**Algorithm 2:** Overview of the algorithm implemented in DDF-PIPELINE-V2. The function *I* represents the imaging step and takes as input the visibility vector **v** together with the beam model **B**<sub>Ωn</sub> and kMS-estimated Jones matrices **J**<sub>Ωn</sub> at locations Ω<sub>n</sub>. The function *K* abstracts the DD calibration step, and takes as arguments the visibilities **v**, the skymodel **x**<sub>v</sub>, a solver mode (estimating for either scalar or full Jones matrices), a time-frequency solution interval (in min and MHz), and a set of directions Ω<sub>n</sub> in which to solve for. The extra functions *C*, *B*, and *F* represent the clustering, bootstrapping and smoothing steps respectively.

**Data:** Visibilities **v** calibrated from DI effects using PREFACTOR of $n_p \times 8$ hours observations (each with 240 LOFAR-HBA subbands), as well as the high resolution skymodel built in step 1.18.

**Result:** Deconvolved image **x**<sub>v</sub>

/* On $n_p \times 240$ LOFAR HBA subbands */

/* DD calibration */

2.1 $\tilde{\mathbf{J}} \leftarrow \mathcal{F} \circ \mathcal{K}(\mathbf{v}_{n_p \times 24}, \mathbf{B}_{\Omega_n}, \tilde{\mathbf{x}}_v|\text{scalar}, 1\text{min}, 2\text{MHz}, \Omega_n)$

/* DI calibration */

2.2 $\mathbf{v}_{n_p \times 24}^c \leftarrow \mathcal{K}(\mathbf{v}_{n_p \times 24}, \tilde{\mathbf{J}}\mathbf{B}_{\Omega_n}, \tilde{\mathbf{x}}_v|\text{full}, \delta t_0, \delta v_0, \Omega_0)$

/* DD imaging */

2.3 $\tilde{\mathbf{x}}_v \leftarrow \mathcal{I}(\mathbf{v}_{n_p \times 24}^c, \tilde{\mathbf{J}}\mathbf{B}_{\Omega_n})$

/* DD calibration */

2.4 $\tilde{\mathbf{J}} \leftarrow \mathcal{F} \circ \mathcal{K}(\mathbf{v}_{n_p \times 24}^c, \mathbf{B}_{\Omega_n}, \tilde{\mathbf{x}}_v|\text{scalar}, 1\text{min}, 2\text{MHz}, \Omega_n)$

/* Slow DD calibration */

2.5 $\tilde{\mathbf{J}}_s \leftarrow \mathcal{K}(\mathbf{v}_{n_p \times 24}^c, \tilde{\mathbf{J}}\mathbf{B}_{\Omega_n}, \tilde{\mathbf{x}}_v|\text{scalar}, 43\text{min}, 2\text{MHz}, \Omega_n)$

/* Final imaging steps */

2.6 $\tilde{\mathbf{x}}_v \leftarrow \mathcal{I}(\mathbf{v}_{n_p \times 24}^c, \tilde{\mathbf{J}}_s, \tilde{\mathbf{J}}\mathbf{B}_{\Omega_n})$;
/* Absolute flux density scale correction (see Sabater et al. 2021, for details) */

2.7 $\tilde{\mathbf{x}}_v \leftarrow f_c \tilde{\mathbf{x}}_v$;

2.8 Facet-based astrometric correction (see Shimwell et al. 2019, for details);

the LOFAR cycle (8.7 h in Cycle 3, 8 h in Cycle 10, see Table 1). The total integration time is ~100 h.

## 4.2. Image synthesis

The Lockman Hole and Boötes field data were reduced using Algorithm 2. In this approach, we first built a wide-band DI+DD self-calibrated sky model **x**<sub>v</sub> from a single wide-band ~8 -h observation using Algorithm 1. This model was then used to DI+DD calibrate all the $n_p$ pointings (with $n_p = 10$ and $n_p = 12$ for the Boötes and Lockman Hole datasets, respectively) following Algorithm 2. This amounts to repeating Steps 1.13 to 1.18 of Algorithm 1 on a larger dataset. A comparison between the images synthesized from the 8- and 80-h datasets is presented in Fig. 11. On a single node equipped with ~500 GB of 2.4 GHz RAM and 2 Intel Xeon CPU E5-2660 v4 at 2.00 GHz with 14 physical cores each, Algorithm 2 took ~21 days to process the 80 h of Boötes data. Figures 12 (further discussed in Sect. 4.3) and 13 show the central parts of the of these deep LOFAR Boötes and Lockman Hole observations.
<figure><img src="image_12.png" /></figure>

**Fig. 12.** Comparison between the LOFAR-HBA maps generated at 150 MHz by Retana-Montenegro et al. (2018) and in the current work. The color scale is the same on all panels and is displayed using an inverse hyperbolic sine function to render both the low level artifacts and some bright sources. (a) The central ≥2 deg² part of the Bootes field as imaged by the direction dependent factor algorithm (Retana-Montenegro et al. 2018). (b) Zoom in on region (1) of the map synthesised by Retana-Montenegro et al. (2018). (c) Zoom-in on region (1) of the map synthesised by kms-ddfacet (this work). (d) The same as in 12a, but imaged with Algorithm 2. (e) Zoom-in on region (2) of the map synthesised by Retana-Montenegro et al. (2018). (f) Zoom-in on region (2) of the map synthesised by kms-ddfacet (this work).
<figure><img src="image_13.png" /><figcaption>**Fig. 13.** Central region of the deep LOFAR-HBA maps of the Lockman Hole field generated at 150 MHz. The color scale is the same on all panels and is displayed using an inverse hyperbolic sine function to render both the low level artifacts and some bright sources. The stripy artifact seen in the zoom-out in *panel b* seems to be produced by the residual deconvolution and calibration errors of a few ≥10 mJy beam⁻¹ bright sources that are a few degrees away from the center of the field. (a) The central ≥2 deg² part of the Lockman Hole field as imaged by Alg. 2. (Sect. 4.2). (b) Zoom in on region (1) of the map shown in Fig. 13a. (c) Zoom-in on region (2) of the map shown in Fig. 13a.</figcaption></figure>

Estimating the noise in radio maps is not straightforward since noise is correlated and non-Gaussian. Furthermore, while the covariance matrix should be entirely described by the PSF, the real covariance matrix is hard to estimate due to the calibration artifacts (see Tasse et al. 2018; Bonnassieux et al. 2018, for a detailed discussion). Here, in order to estimate the local noise, we used the statistics of the min{.} estimator (which returns the minimum value of a given sample). Intuitively, while the I-Stokes image max{.} statistics have contributions from both artifacts and real sources, the min{.} only accounts for the artifacts. A min{.} filter with a given box size is therefore run through a restored image, and, depending on the box size<fnref n="13" />, the effective standard deviation is derived.

Figure 14 shows the cumulative distribution of the local noise in the Lockman Hole and Boötes fields maps, reaching ≤23 and ≤30 µJy beam⁻¹, respectively. Taking into account the number of pointings with their respective amounts of flagged data, we get total integration times of ~65 and ~88 h on the Boötes and Lockman Hole fields, respectively, giving a theoretical thermal noise difference of a factor of ~1.16, which is compatible with the observed value of ~1.3. Other factors to be taken into account to compare noise properties include the bootstrapping errors, the individual fields' average elevations, and the Galactic noise differences.

<figure><img src="image_14.png" /><figcaption>**Fig. 14.** Cumulative distribution of the local noise estimates in the various maps discussed here. As shown, we have imaged a larger fraction of LOFAR’s HBA primary beam than the image presented in Retana-Montenegro et al. (2018).</figcaption></figure>

### 4.3. Comparison with deep FACTOR/ image synthesis

The image of the Boötes field based on 55 h of LOFAR-HBA data and presented in Retana-Montenegro et al. (2018) reaches an unprecedented noise level image of ~55 µJy beam⁻¹ at 150 MHz. To achieve such high sensitivity, Retana-Montenegro et al. (2018) applied third generation calibration and imaging to correct for the DDEs using the FACTOR package (developed by

<footnote n="13">¹³ The cumulative distribution $\mathcal{F}$ of $Y = \min\{X\}$ with $X \sim N\{\mu = 0, \sigma = 1\}$ is $\mathcal{F}\{y\} = 1 - \left[ \frac{1}{2} \left( 1 - \operatorname{erf}\left\{ \frac{y}{\sqrt{2}} \right\} \right) \right]^n$, where $n$ is the number of pixels in a given box. Finding $y_\sigma$ such that $\mathcal{F}\{y_\sigma\} = 1/2$ given the box size gives us a conversion factor from the minimum estimate to the standard deviation.</footnote>
van Weeren et al. 2016, see Sect. 2 for more detail). Because the set of LOFAR datasets used by Retana-Montenegro et al. (2018) is different<fnref n="14" /> the comparison can only be approximate. In Fig. 12 we compare the images produced by Retana-Montenegro et al. (2018) and by Algorithm 2. While the noise difference should be on the order of 20%, as shown in Fig. 14 the measured one is on the level of ~60%. Consistently artifacts around bright sources are also much less severe in the maps generated by Algorithm 2 and implemented in DDF-PIPELINE-V2.

### 4.4. Cataloging

In order to extract astrophysical information, we built a catalog of radio sources from the images produced by Algorithm 2 and the data described in Sect. 4. Even in the apparent flux maps, because of the imperfect calibration and imaging, the LOTSS-Deep Fields images have spatially variable noise; to deal with this issue, we used PyBDSF<fnref n="15" /> (Python Blob Detector and Source Finder, see Mohan & Rafferty 2015) since it measures noise locally rather than globally. The sources were detected with 3 and 5σ for the island and peak detection threshold, respectively. The position-dependent noise was estimated using a sliding box algorithm with a size of 40 × 40 synthesized beams, except around bright sources where the box size was decreased to 15 × 15 beams to more accurately capture the increased noise in these regions. The columns kept in the final catalog are the source position, peak and integrated flux density, source size and orientation, the associated uncertainties, the estimated local rms at the source position, and a code describing the type of structure fit by PyBDSF. As described in Sabater et al. (2021), the peak and integrated flux densities of the final catalogs and images are corrected from overall scaling factors of 0.920 and 0.859 for the Lockman Hole and Boötes fields, respectively. These numbers were estimated from the comparison between the LOTSS-Deep Fields flux densities and a variety of radio data available at other frequencies. The full catalogs cover out to 30% of the power primary beam and contain 36 767 entries over 26.5 square degrees and 50 112 over 25.0 square degrees for Boötes and Lockman Hole, respectively. These raw PyBDSF catalogs are available online at the LOFAR survey web page<fnref n="16" />, and a thorough analysis of the source catalogs will be presented by Mandal et al. (2021).

## 5. Conclusion and future plans

Imaging low frequency LOFAR data at high resolution and over wide fields of view is extremely challenging. This is mainly due to the RIME system being complex in this regime: the background wide-band sky is unknown, as are the time-frequency-antenna DD-Jones matrices. Due to the high number of free parameters in that system, and due to the finite amount of data points in the nonlinear RIME system, the inversion can be subject to ill-conditioning and the DD-C-RIME solver can absorb unmodeled extended flux.

In order to address this robustness issue, we have developed a strategy that aims at conserving the unmodeled emission without affecting the final dynamic range. The method we have developed has similarities with those presented by Yatawatta (2015),

van Weeren et al. (2016), Repetti et al. (2017), Birdi et al. (2020) and relies on reducing the effective size of the unknown stochastic process. We show that this allows us to recover most of the faint unmodeled extended emission.

We have applied this third generation calibration and imaging DD algorithm to both the wide-field imaging of the LOTSS survey and to the synthesis of deep 150 MHz resolution images on the Boötes and Lockman Hole fields. The synthesized images are the deepest ever obtained at these frequencies. Detailed analyses of the LOTSS-Deep Fields catalogs (including the source counts of the Lockman Hole, Boötes, and ELAIS-N1 fields) are presented in Mandal et al. (2021). In the future, we plan to continue increasing the depths of these fields: The data needed to double the integration time on each field, with a further aim to increase this to 500 h in each field, are already in hand (or are scheduled).

*Acknowledgements.* This paper is based (in part) on data obtained with the International LOFAR Telescope (ILT). LOFAR (van Haarlem et al. 2013) is the Low Frequency Array designed and constructed by ASTRON. It has observing, data processing, and data storage facilities in several countries, which are owned by various parties (each with their own funding sources), and which are collectively operated by the ILT foundation under a joint scientific policy. The ILT resources have benefitted from the following recent major funding sources: CNRS-INSU, Observatoire de Paris and Université d'Orléans, France; BMBF, MIWF-NRW, MPG, Germany; Science Foundation Ireland (SFI), Department of Business, Enterprise and Innovation (DBEI), Ireland; NWO, The Netherlands; The Science and Technology Facilities Council, UK; Ministry of Science and Higher Education, Poland. This work makes use of kernt astronomical software package (available at https://kernsuite.info and presented in Molenaar & Smirnov 2018). M.B. acknowledges support from INAF under PRIN SKA/CTA FORECaST. M.B. acknowledges the support from the Ministero degli Affari Esteri della Cooperazione Internazionale - Direzione Generale per la Promozione del Sistema Paese Progetto di Grande Rilevanza ZA18GR02. M.J.J. acknowledges support from the UK Science and Technology Facilities Council [ST/N000919/1] and the Oxford Hintze Centre for Astrophysical Surveys which is funded through generous support from the Hintze Family Charitable Foundation. P.N.B. and J.S. are grateful for support from the UK STFC via grant ST/R000972/1. M.J.H. acknowledges support from STFC via grant ST/R000905/1. W.L.W. acknowledges support from the ERC Advanced Investigator programme NewClusters 321271. W.L.W. also acknowledges support from the CAS-NWO programme for radio astronomy with project number 629.001.024, which is financed by the Netherlands Organisation for Scientific Research (NWO). A.B. acknowledges support from the VIDI research programme with project number 639.042.729, which is financed by the Netherlands Organisation for Scientific Research (NWO). I.P. acknowledges support from INAF under the SKA/CTA PRIN "FORECaST" and the PRIN MAIN STREAM "SAuROS" projects M.B. acknowledges support from INAF under PRIN SKA/CTA FORECaST and from the Ministero degli Affari Esteri della Cooperazione Internazionale - Direzione Generale per la Promozione del Sistema Paese Progetto di Grande Rilevanza ZA18GR02. R.K. acknowledges support from the Science and Technology Facilities Council (STFC) through an STFC studentship.

## References

* Ashby, M. L. N., Stern, D., Brodwin, M., et al. 2009, ApJ, 701, 428

* Becker, R. H., White, R. L., & Helfand, D. J. 1995, ApJ, 450, 559

* Bhatnagar, S., & Cornwell, T. J. 2017, AJ, 154, 197

* Bhatnagar, S., Cornwell, T. J., Golap, K., & Uson, J. M. 2008, A&A, 487, 419

* Biggs, A. D., & Ivison, R. J. 2006, MNRAS, 371, 963

* Birdi, J., Repetti, A., & Wiaux, Y. 2020, MNRAS, 492, 3509

* Bonnassieux, E., Tasse, C., Smirnov, O., & Zarka, P. 2018, A&A, 615, A66

* Brentjens, M. A., & de Bruyn, A. G. 2005, A&A, 441, 1217

* Brown, M. J. I., Dey, A., Jannuzi, B. T., et al. 2007, ApJ, 654, 858

* Brown, M. J. I., Zheng, Z., White, M., et al. 2008, ApJ, 682, 937

* Brunner, H., Cappelluti, N., Hasinger, G., et al. 2008, A&A, 479, 283

* Butler, A., Huynh, M., Delhaize, J., et al. 2018, A&A, 620, A3

* Ciliegi, P., Zamorani, G., Hasinger, G., et al. 2003, A&A, 398, 901

* Condon, J. J., Cotton, W. D., Greisen, E. W., et al. 1998, AJ, 115, 1693

* Cool, R. J. 2007, ApJS, 169, 21

* Coppejans, R., Cseh, D., Williams, W. L., van Velzen, S., & Falcke, H. 2015, MNRAS, 450, 1477

<footnote n="14">¹⁴ Out of the sets of 7 and 10 observations used in Retana-Montenegro et al. (2018) and in this work respectively, 4 are common, namely L243561, L374583, L400135, L401825.</footnote>

<footnote n="15">¹⁵ https://www.astron.nl/citt/pybdsf</footnote>

<footnote n="16">¹⁶ https://www.lofar-surveys.org/</footnote>
Coppin, K., Chapin, E. L., Mortier, A. M. J., et al. 2006, MNRAS, 372, 1621

Croft, S., van Breugel, W., Brown, M. J. I., et al. 2008, AJ, 135, 1793

de Gasperin, F., Dijkema, T. J., Drabent, A., et al. 2019, A&A, 622, A5

de Ruiter, H. R., Zamorani, G., Parma, P., et al. 1997, A&A, 319, 7

de Vries, W. H., Morganti, R., Röttgering, H. J. A., et al. 2002, AJ, 123, 1784

Dewdney, P. E., Hall, P. J., Schilizzi, R. T., & Lazio, T. J. L. W. 2009, IEEE Proc., 97, 1482

Duncan, K. J., Kondapally, R., Brown, M. J. I., et al. 2021, A&A, 648, A4 (LoTSS SI)

Fomalont, E. B., Kellermann, K. I., Cowie, L. L., et al. 2006, ApJS, 167, 103

Garn, T. S., Green, D. A., Riley, J. M., & Alexander, P. 2010, Bull. Astron. Soc. India, 38, 103

Geach, J. E., Dunlop, J. S., Halpern, M., et al. 2017, MNRAS, 465, 1789

González-Solares, E. A., Irwin, M., McMahon, R. G., et al. 2011, MNRAS, 416, 927

Guglielmino, G., Prandoni, I., Morganti, R., & Heald, G. 2012, Resolving The Sky - Radio Interferometry: Past, Present and Future, 22

Hamaker, J. P. 2000, A&AS, 143, 515

Hamaker, J. P., Bregman, J. D., & Sault, R. J. 1996, A&AS, 117, 137

Herrera Ruiz, N., O'Sullivan, S. P., Vacca, V., et al. 2020, MNRAS, 463, 2997

Heywood, I., Jarvis, M. J., Baker, A. J., et al. 2016, MNRAS, 460, 4433

Ibar, E., Ivison, R. J., Biggs, A. D., et al. 2009, MNRAS, 397, 281

Intema, H. T., van Weeren, R. J., Röttgering, H. J. A., & Lal, D. V. 2011, A&A, 535, A38

Intema, H. T., Jagannathan, P., Mooley, K. P., & Frail, D. A. 2017, A&A, 598, A78

Jannuzi, B. T., & Dey, A. 1999, ASP Conf. Ser., 193, 258

Jannuzi, B., Weiner, B., Block, M., et al. 2010, BAAS, 42, 513

Kazemi, S., Yatawatta, S., Zaroubi, S., et al. 2011, MNRAS, 414, 1656

Kenter, A., Murray, S. S., Forman, W. R., et al. 2005, ApJS, 161, 9

Kondapally, R., Best, P. N., Hardcastle, M. J., et al. 2021, A&A, 648, A3 (LoTSS SI)

Lacy, M., Baum, S. A., Chandler, C. J., et al. 2020, PASP, 132, 035001

Lawrence, A., Warren, S. J., Almaini, O., et al. 2007, MNRAS, 379, 1599

Lonsdale, C. J., Smith, H. E., Rowan-Robinson, M., et al. 2003, PASP, 115, 897

Mahony, E. K., Morganti, R., Prandoni, I., et al. 2016, MNRAS, 463, 2997

Mandal, S., Prandoni, I., Hardcastle, M. J., et al. 2021, A&A, 648, A5 (LoTSS SI)

Martin, C., & GALEX Team. 2005, IAU Symp., 216, 221

Masini, A., Hickox, R. C., Carroll, C. M., et al. 2020, ApJS, 251, 1

Mauch, T., Cotton, W. D., Condon, J. J., et al. 2020, ApJ, 888, 61

Mohan, N., & Rafferty, D. 2015, PyBDSF: Python Blob Detection and Source Finder

Molenaar, G., & Smirnov, O. 2018, Astron. Comput., 24, 45

Murray, S. S., Kenter, A., Forman, W. R., et al. 2005, ApJS, 161, 1

Norris, R. 2010, AAS Meeting Abstracts, 215, 604.05

Offringa, A. R., van de Gronde, J. J., & Roerdink, J. B. T. M. 2012, A&A, 539, A95

Oliver, S. J., Bock, J., Altieri, B., et al. 2012, MNRAS, 424, 1614

O'Sullivan, S., Brüggen, M., Van Eck, C., et al. 2018, Galaxies, 6, 126

O'Sullivan, S. P., Brueggen, M., & Vazza, F. 2020, MNRAS, 495, 2607

Owen, F. N., & Morrison, G. E. 2008, AJ, 136, 1889

Owen, F. N., Morrison, G. E., Klimek, M. D., & Greisen, E. W. 2009, AJ, 137, 4846

Padovani, P. 2016, A&ARv, 24, 13

Pearson, T. J., & Readhead, A. C. S. 1984, ARA&A, 22, 97

Polletta, M. d. C., Wilkes, B. J., Siana, B., et al. 2006, ApJ, 642, 673

Prandoni, I. & Seymour, N. 2015, Advancing Astrophysics with the Square Kilometre Array (AASKA14), 67

Prandoni, I., Guglielmino, G., Morganti, R., et al. 2018, MNRAS, 481, 4548

Rengelink, R. B., Tang, Y., de Bruyn, A. G., et al. 1997, A&AS, 124, 259

Repetti, A., Birdi, J., Dabbech, A., & Wiaux, Y. 2017, MNRAS, 470, 3981

Retana-Montenegro, E., Röttgering, H. J. A., Shimwell, T. W., et al. 2018, A&A, 620, A74

Richards, E. A. 2000, ApJ, 533, 611

Sabater, J., Best, P. N., Tasse, C., et al. 2021, A&A, 648, A2 (LoTSS SI)

Schinnerer, E., Carilli, C. L., Scoville, N. Z., et al. 2004, AJ, 128, 1974

Shimwell, T. W., Röttgering, H. J. A., Best, P. N., et al. 2017, A&A, 598, A104

Shimwell, T. W., Tasse, C., Hardcastle, M. J., et al. 2019, A&A, 622, A1

Smirnov, O. M. 2011, A&A, 527, A108

Smirnov, O. M., & Tasse, C. 2015, MNRAS, 449, 2668

Smolčić, V., Novak, M., Bondi, M., et al. 2017a, A&A, 602, A1

Smolčić, V., Novak, M., Delvecchio, I., et al. 2017b, A&A, 602, A6

Sob, U. M., Bester, H. L., Smirnov, O. M., Kenyon, J. S., & Grobler, T. L. 2020, MNRAS, 491, 1026

Tasse, C. 2014a, ArXiv e-prints [arXiv:1410.8706]

Tasse, C. 2014b, A&A, 566, A127

Tasse, C., Röttgering, H., & Best, P. N. 2010, A&A, 525, A127

Tasse, C., Hugo, B., Mirmont, M., et al. 2018, A&A, 611, A87

Van Eck, C. L., Haverkorn, M., Alves, M. I. R., et al. 2018, A&A, 613, A58

van Haarlem, M. P., Wise, M. W., Gunst, A. W., et al. 2013, A&A, 556, A2

van Weeren, R. J., Williams, W. L., Tasse, C., et al. 2014, ApJ, 793, 82

van Weeren, R. J., Williams, W. L., Hardcastle, M. J., et al. 2016, ApJS, 223, 2

van Weeren, R. J., Shimwell, T. W., Botteon, A., et al. 2020, ArXiv e-prints [arXiv:2011.02387]

Wayth, R. B., Lenc, E., Bell, M. E., et al. 2015, PASA, 32, e025

Williams, W. L., Intema, H. T., & Röttgering, H. J. A. 2013, A&A, 549, A55

Williams, W. L., van Weeren, R. J., Röttgering, H. J. A., et al. 2016, MNRAS, 460, 2385

Yatawatta, S. 2015, MNRAS, 449, 4506

Yatawatta, S., Zaroubi, S., de Bruyn, G., Koopmans, L., & Noordam, J. 2008, ArXiv e-prints [arXiv:0810.5751]

Zwart, J., Wall, J., Karim, A., et al. 2015, Advancing Astrophysics with the Square Kilometre Array (AASKA14), 172
# Appendix A: LOTSS first data release: overview of
## DDF-PIPELINE-V1

The data processing strategy of the LOTSS first data release has been extensively described by Shimwell et al. (2019). Since addressing the issues described in Sect. 2 involves making improvements relative to this approach, we give a brief description here of the data reduction strategy in DDF-PIPELINE-V1 (the various steps are outlined in Algorithm 0).

**Algorithm 0:** Overview of the algorithm implemented in DDF-PIPELINE-V1 to produce the LOTSS-DR1 images. The function *I* represents the imaging step and takes as input the visibility vector **v** together with the beam model **B**<sub>Ω<sub>n</sub></sub> and kMS-estimated Jones matrices **J**<sub>Ω<sub>n</sub></sub> at locations Ω<sub>n</sub>. The function *K* abstracts the DD calibration step, and takes as arguments the visibilities **v**, the skymodel **x̃**<sub>**v**</sub>, a solver mode (estimating for either scalar or full Jones matrices), a time-frequency solution interval (in min and MHz), and a set of directions Ω<sub>n</sub> in which to solve for. The extra functions *C*, *N*, and *B* represent the clustering, normalisation (see text), and bootstrapping steps respectively.

**Data:** Visibilities **v** calibrated from DI effects using PREFACTOR.

**Result:** Deconvolved image **x̃**<sub>**v**</sub>

/* On 60 LOFAR HBA subbands */

/* DI initial deconv and clustering */

0.1 **x̃**<sub>**v**</sub> ← *I* (**v**<sub>6</sub>, **B**<sub>Ω<sub>n</sub></sub>);

0.2 Ω<sub>n</sub> ← *C* (**x̃**<sub>**v**</sub>);

/* Phase only DD calibration */

0.3 **J** ← φ ◦ *K* (**v**<sub>6</sub>, **x̃**<sub>**v**</sub>, **B**<sub>Ω<sub>n</sub></sub>|scalar, 1min, 2MHz, Ω<sub>n</sub>);

/* Absolute flux density scale bootstrapping */

0.4 **ṽ** ← *B* (**ṽ**<sub>6</sub>);

0.5 **x̃**<sub>**v**</sub> ← *I* (**v**<sub>6</sub>, **JB**<sub>Ω<sub>n</sub></sub>);

/* DD calibration and imaging */

0.6 **J** ← *N* ◦ *K* (**v**<sub>6</sub>, **x̃**<sub>**v**</sub>, **B**<sub>Ω<sub>n</sub></sub>|scalar, 1min, 2MHz, Ω<sub>n</sub>);

0.7 **x̃**<sub>**v**</sub> ← *I* (**v**<sub>6</sub>, **JB**<sub>Ω<sub>n</sub></sub>);

/* On 240 LOFAR HBA subbands */

/* Deep DD calibration and imaging */

0.8 **J** ← *N* ◦ *K* (**v**<sub>24</sub>, **x̃**<sub>**v**</sub>, **B**<sub>Ω<sub>n</sub></sub>|scalar, 1min, 2MHz, Ω<sub>n</sub>);

0.9 **x̃**<sub>**v**</sub> ← *I* (**v**<sub>24</sub>, **JB**<sub>Ω<sub>n</sub></sub>);

0.10 Facet-based astrometric correction (see Shimwell et al. 2019, for details);

As discussed in Sect. 2, the calibration and imaging problem is non-convex and ill-posed. Beyond the computational issues, the great difficulty of the DDE calibration is sky incompleteness because the DD-C-RIME nonlinear system can be subject to ill-conditioning. This is due to the fact that the extended emission (i) is hard to model in the deconvolution step and (ii) is seen by only the shortest baselines and therefore sky incompleteness biases the Jones matrices in the calibration step. Experience shows that this leads to some of the unmodeled extended emission being absorbed when running a DD deconvolution with DDFacet.

To try to compensate for this effect, in DDF-PIPELINE-V1 (Algorithm 0) we introduced an inner uv-distance cut during calibration, as well as a normalization of the Jones matrix. With this, the DDF-PIPELINE-V1 was able to recover some of the unmodeled extended emission. The underlying idea was to assume that the sky incompleteness was generating some baseline-dependent systematic errors. So for every given direction and solution interval in Shimwell et al. (2019) we were trying to find a gain vector **g** such that $gg^H \sim \langle g_{tv}g_{tv}^H \rangle$ (where $A^H$ is the Hermitian transpose of matrix **A**). This amounts to constraining the baseline-dependent error to be solely antenna-dependent. This normalization (described by the function symbol *N* in Algorithm 0) was able to recover some extended emission otherwise absorbed in the calibration solution. However, as shown in Fig. 8a and explained by Shimwell et al. (2019), it also produced large-scale fake haloes centered on extended sources together with artifacts around bright sources. On fields that have a bright $\ge 1$ Jy source within the primary beam (such as 3C sources), DDF-PIPELINE-V1 was not able to converge.