The LOFAR Two Meter Sky Survey: Deep Fields, I -- Direction-dependent calibration and imaging

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 Bootes and LH fields at 150 MHz, totaling $\sim80$ and $\sim100$ hours of integration respectively and reaching unprecedented noise levels at these low frequencies of $\lesssim30$ and $\lesssim23$ $\mu$Jy/beam in the inner $\sim3$ deg$^2$. This approach is also being used to reduce the LoTSS-wide data for the second data release.


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 building tiered extragalactic surveys with LOFAR, of different depth and areas, and at frequen-cies ranging from ∼ 30 to ∼ 200 MHz.Specifically the LO-FAR LBA Sky Survey (LoLSS) aims at surveying the northern hemisphere using the LOFAR LBA antennas while the LOFAR Two Meter Sky Survey (LoTSS) uses the High Band Antennas (HBA).Its widest component (LoTSS-wide) has been described by Shimwell et al. (2017aShimwell et al. ( , 2019)), and aims at reaching noise levels of 100 µJy.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 target noise levels of ultimately 10 µJy.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.Fig. 1 is inspired by that of Smolvcić et al. (2017a) and shows a sensitivity and surveyed area comparison between various existing and future surveys.These include TGSS (Intema et al. 2017), FIRST (Becker et al. 1995), NVSS (Condon et al. 1998), VLA-COSMOS (Schinnerer et al. 2004;Smolvcić et al. 2017b), VLASS (Lacy et al. 2020), EMU (Norris 2010), 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), LOTSS-DR1 (Shimwell et al. 2017b), HDF-N (Richards 2000), WENSS (Rengelink et al. 1997), GLEAM (Wayth et al. 2015), and SKA (Prandoni & Seymour 2015a).
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 Square Kilometer Array (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 2015b).In order to construct the LoTSS-Deep Fields, we have selected the Boötes, Lockman Hole, and ELAIS-N1 fields, together with the North Ecliptic Pole (NEP).Each of them is covered by a wealth of multiwavelength data, necessary to derive photometric redshifts and low frequency radio luminosities, thereby providing an efficient way to estimate Star Formation Rate (SFR hereafter) in galaxies for example.Together, these four multiwavelength fields allow us to probe a total sky area of 30 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 (dde), such as the ionospheric distortions, and 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 were able to estimate and compensate for the dde, and thus use LOFAR to produce thermal noise limited maps from 8-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 ∼ 30 µJy.beam −1 noise levels, while being more robust against absorbing faint unmodeled extended emission and dynamic range issues.In Sec. 2 we introduce the 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 Sec. 3 (for completion in appendix A we describe ddf-pipeline-v1 that was presented in detail in Shimwell et al. 2019).In Sec. 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. 2020, in prep.), the multiwavelength cross matching (Kondapally et al. 2020, in prep.) and the photometric redshifts and galaxy characterisation (Duncan et al. 2020, in prep.).

LoTSS and the third generation calibration and imaging problem
Calibration and imaging techniques have greatly evolved since the first radio interferometers have become operational.First generation calibration is commonly refered as directionindependent (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 direction-dependent effects (dde).
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 frequency (ν < 250 MHz), (ii) has very large fields of view (fwhm of 2 -10 degrees), and (iii) combines short (∼ 100 m) and long (∼ 2000 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) make the calibration problem direction-dependent by nature.In Sec.2.1 we introduce the mathematical formalism used throughout this paper, while in Sec.2.2 and 2.3 we describe the algorithms and software used for di and dd calibration and imaging.

Measurement equation formalism
The Radio Interferometry Measurement Equation (rime, see Hamaker et al. 1996) describes how the underlying electric field coherence (the sky model), and the the various directionindependent 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 T (where T is the matrix transpose) and write the 4 × 1 sky coherency matrix as ptν are the directionindependent and direction-dependent 4 × 4 Mueller matrices 1 on a baseline b ↔ {pqt} → [u, v, w] T between antenna p and q at time t, then the 4-visibility v b is given by where c is the speed of light in the vacuum, and n b is a 4 × 1 random matrix following a normal distribution N (0, σ b ).Depending on the context, in the rest of this paper we will either make use of the antenna-based Jones matrices or the baselinebased 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 s b models the dde including the ionospheric distortion (phase shift, Faraday rotation, scintillative decoherence) and phased array station beam that depend on time, frequency, and antenna.Importantly, in order to estimate the intrisic flux densities we use a description of the LO-FAR station primary beam that is built from semi-analytic models 2 , and write it as B s b 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 s b and the sky model x ν from the set of visibilities v b .Due to the bilinear structure of the rime, instead of estimating all these parameters at once, inverting Eq. 1 is split into two steps.In the first step, the sky term x ν is assumed to be constant, and the Jones matrices are estimated.The step is referred as "calibration" and as the dd-Crime 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 ν is estimated.This step is commonly called "imaging" and is referred 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-selfcalibration.
While the vast majority of modern developments in the field of algorithmic research for radio interferometry aim at addressing either the dd-C-rime (direction dependent 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. 2010Tasse et al. , 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.
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 × 2 Jones matrices of antenna p and q, and X is the 2 × 2 source's coherency matrix then we have Vec J p XJ q = J * q ⊗ J p Vec {X}, where ⊗ is the Kronecker product. 2 https://github.com/lofar-astron/LOFARBeam

Direction-independent calibration
The standard LoTSS HBA observations consist of a 10 minute scan on a bright calibrator source (in general 3C 196 or 3C 295) before observing the target field for 8 hours.On both calibrator and target fields, the visibilities of the 240 subbands (SB) are regularly distributed across the 120-168 MHz bandpass, with 64 channels per 195.3 kHz subband and 1 sec integration time.The data are first flagged using AOFlagger3 (Offringa et al. 2012) and 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 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 package4 (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 where v model b = s x sν k s b ds is the skymodel of the calibrator.We cannot just use G 0 b 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 G 0 b can be decomposed as the product of a direction independent and direction dependent term and where K = 8.44 × 10 9 m 3 s −2 , and A pν , ∆ t p and ∆ T p are realvalued and represent respectively the bandpass, the clock and ionospheric Total Electron Content offset with respect to a reference antenna.The terms ∆ t p and ∆ T p can be disentangled from the frequency dependent phases because the former are linear with ν while the latter are non-linear.
Assuming the clocks and bandpass are the same for the calibrator and for the target field, the corrected visibilities v c b for the target field can be built from the raw data In order to calibrate for the remaining phase errors, the target field is di calibrated against the TIFR GMRT Sky Survey (tgss) catalogs (Intema et al. 2017) and visibilities are averaged to 2 ch/sb and 8s.

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 s ptν and J s ptν .Depending on the observing frequency and instrumental setup, these can be the product of the ionospheric Faraday rotation matrix, the scalar phase shift, and 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 model analytically the physics of the various effects over the {sptν} space to be able to disentangle them.The second family of algorithms directly estimate the effective Jones matrices over piecewise constant domains in {sptν} space.These "Jonesbased" solvers have the advantage of not requiring any physical modeling of the DDE.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 to be ill-posed5 .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, and also on the shape of the neighboring domains in the {sptν} spaces.Also, 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 sν 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 v b that are estimated6 over the piecewise constant domains p, Ω ϕ , ∆t, ∆ν can be written as where Ω ϕ is the set of directions s for a facet ϕ, x sν is the estimated underlying sky, and G b and J ϕ b are the di and dd Mueller matrices for baseline b, built from the corresponding estimated Jones matrices in p, ∆t, ∆ν.
Specifically, in order to solve for the dde, 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 non-linear nature of Eq. 1, and (ii) the complexity of the background radio sky, optimising 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 x sν as well as the di G ptν and the dd piecewise constant J ϕ ptν for all {sϕptν} such that v b ∼ v b (we remain intentionally vague here, since the cost function that is minimised 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 (imaging step).Using this skymodel and repeating steps (i) and (ii) is known as self-calibration, but in a third-generation approach we must explicitly model the DD aspects.
Since the computing time evolves as ∼ n 3 d , where n d is the number of directions in the dd-solvers, the problem of ddcalibration has in general been tackled using direction alternating peeling-like techniques.Major breakthroughs have been made in the field of dd-C-rime solvers in 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 a different structure and shortcuts can be taken to invert the calibration problem.The net gain over the classical method is not trivial, but can be as high as n 2 a (Smirnov & Tasse 2015) where n a 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 Jones-based solvers.
Only a very few CI-rime algorithms using a full dd selfcalibration loop have been described and implemented.They include pointing self-calibration (Bhatnagar & Cornwell 2017), or peeling-based techniques such as MFImage (implemented in the obit7 package) and factor (van Weeren et al. 2016, see also Sec. 4.3).Similar to peeling, and developed for reducing LO-FAR data, factor is sequential along the direction axis.Looping over the different facets it consists of (i) subtracting all sources besides calibration sources in that one facet, and (ii) di-selfcalibrating in that direction.In addition to the ill conditioning issues discussed above on dd-C-rime and dd-I-rime solvers, an expensive computational problem arises when estimating the J ϕ ptν .The approach presented by Shimwell et al. (2019) (also described in detail Sec.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-pipeline8 is a high level wrapper that mainly calls ddfacet9 and kms10 for direction dependent 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, i.e. the dd affected sidelobes leaking from any facet to any other facet are accounted for within the algebraic operations of the algorithm.Another advantage compared to the factor approach is that the data need only to be read rather than modified, making the ddf-pipeline more I/O efficient.
As explained in Shimwell et al. (2019), the ddf-pipeline-v1 was however affected by a number of issues including (i) artifacts and dynamic range limit around the brightest radio sources, (ii) artificial and diffuse haloes around moderately bright radio sources and (iii) unmodeled flux absorption.

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 cal-Fig.2: Some of the images produced at different steps in the dd-self-calibration loop implemented as Alg. 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 facets' locations.The colorscale is the same on all panels, and diplayed using an inverse hyperbolic sine function to render both the low level artifacts and some bright sources.ibration solver more robust against artifacts around the brightest sources (Sec.3.1) and unmodeled flux absorption (Sec.3.2 and Sec.3.3).An overview of this approach is shown in Alg. 1 (the implementation of which is referred as ddf-pipeline-v2), and the corresponding dd self-calibration loop is displayed in Fig. 2. The ddf-pipeline-v2 products are described in Sec.3.4.In Sec.3.5 we show that ddf-pipeline-v2 produces improved images as compared to those previously made with ddf-pipeline-v1 (see Sec.A and Shimwell et al. 2019).In Sec.3.6 we discuss the ddf-pipeline-v2 computing efficiency.

Dynamic range issue
With LOFAR's very large field of view, it is quite common to observe bright sources within the station's primary beam.It turns out from tests we conducted on fields containing bright sources (such as 3C 295 which has a flux density of ∼ 100 Jy at 150 MHz) that the related residual errors create powerful artifacts that largely dominate the thermal noise, thereby driving a dynamic range limit.As explained in Sec.2.2 the initial phase calibration is done against tgss at 150 MHz (Intema et al. 2017).However, since LoTSS resolution is much higher than tgss's (6 ×6 against 25 ×25 respectively), small spatial uncertainties Data: Visibilities v calibrated from di effects using PreFactor.
/* On 60 LOFAR HBA subbands */ /* DI initial deconv and clustering */ /* DI calibration and imaging */ /* DI calibration and imaging */ Shimwell et al. 2019, for details); 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 (that due to ill conditioning issues has to be done on larger timefrequency solutions 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.
Fig. 3: This plot shows the amplitude of the diagonal (black) and off-diagonal (gray) terms of the estimated Jones matrices for a remote station, using the LOFAR's 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 would be good enough, the calibration solutions found here would be the unity matrix at all times and frequency.Therefore, in ddf-pipeline-v2 we carry out a few full-Jones dionly self-calibration steps.This approach seems very efficient in increasing the dynamic range around the brightest radio sources.
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, these are given with respect to a reference antenna.When Jones matrices are scalar, this amounts to zeroing the phases ϕ 0 of the reference antenna, by subtracting ϕ 0 from all phases of all antennas.To do this in the general case of non-diagonal Jones matrices, we use a polar decomposition on the Jones matrix J 0 of the reference antenna such that J 0 = UP 0 where U is a unitary matrix11 .We then apply U to all Jones matrices as J p ← U H J p .Intuitively, when the Jones matrices are all scalar, the unitary matrix U is simply exp (iϕ 0 )I, and that step makes the phases of all J p relative to the reference antenna (and specifically zeros the phases of J 0 ).In the case of non-trivial 2 × 2 Jones matrices, finding and applying U has the effect of removing a common rotation from all Jones matrices, and orthogonalises them.
We apply this in Fig. 3 where we show the typical di Jones matrices we can estimate at stage 1.3 for a given remote station and frequency, with respect to a reference station in the LOFAR core.They are estimated using kms and the skymodel 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 polarisation changes across the LO-FAR 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 modeled visibilities generated by ddfacet in Step 1.3.The sky being mostly unpolarised, in this full-Jones di calibration step, we assume Q = U = V = 0 Jy (see Sec.Note that after the initial dd calibration solutions have been 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, on any baseline b the model visibilities v Σ b (Eq.10) are predicted based on the previously estimated dd-Jones matrices J (as is done by Smirnov 2011).

Regularisation
The absorption of unmodeled flux by calibration is a well known issue connected to the calibration of dde.Intuitively speaking, when real flux is missing from the modeled sky x i of x at step i, and since the rime inversion is often ill-posed, the estimates J i of J can be biased in a systematic way.Experience and simulations show that building a new estimate x i+1 from J i can be biased in that the unmodeled emission is not and will never be recovered (Fig. 5).This effect is especially severe when the extended emission is poorly modeled or unmodeled since this is detected only by the shortest baselines.Effectively, during the inversion of the rime system of equations, the dd-self-calibration algorithm has fallen into the wrong (local) minimum.
In order to address this problem, one idea is to reduce the effective number of free parameters used to describe the Jones matrices in the {pdtν}-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 matrices' shape to look like that of the real underlying ones improves the conditioning of the inverse problem.In Alg. 1 (implemented in ddfpipeline-v2) we have replaced that normalisation method by a smoothing of the kms-estimated Jones matrices.This function F updates the Jones matrices J ← F (J) in each direction independently by imposing on them a certain behavior in the timefrequency space (see below), effectively reducing the size of the unknown stochastic process.This can be thought of as a regularisation.This is done independently on the phases and amplitudes on the scalar Jones matrices generated at Steps 1.9, 1.13, 1.16.The updated Jones matrices take the analytical form where ∆ T pd,t is the differential tec (see also Sec. 2.2 and Eq.8), P is a polynomial parametrised by the coefficients in θ pd,ν (of size 10), and a pd,t is a scalar meant to describe the loss of correlation due to ionospheric scintillation as seen in the left panel of Fig. 4. Typically, for the ∼ 8 hours' integration of LoTSS pointings and solving every 30 sec.and 2 MHz, this parametrisation of the Jones matrices reduces the number of free parameters by a factor 20.
To assess the recovery of unmodelled 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 LO-FAR 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 halos of galaxy clusters.After the injection of the artificial extended sources the steps 1.16 and 1.17 were repeated using the sky model derived at step 1.18 prior to the injection of the sources.These simulations will be discussed further by Shimwell et al. (in preparation)  As suggested by the results of 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 across, with a mean flux density at the peak of only 0.7 of the local standard deviation.

Conditioning and solution interval
The additional issue of arcmin-scale negative haloes appearing around bright compact sources (at a level of 1% or the peak) could be seen however in ∼ 10 − 20% of the LoTSS pointings processed with ddf-pipeline-v1.As shown in Fig. 6, we believe this to be connected to the solution regularisation itself.This issue is hard to understand in detail because of the non-linearity in the C-rime inversion, but is likely to be due to the pointings showing these issues being more severely affected by the incompleteness of the sky model.Specifically, conducting several experiments, we were able to observe that the situation was improved by deconvolving deeper or taking into account sources outside the synthesized image field of view.
An additional way to improve the conditioning of the problem is to increase the amount of data used to contrain the estimated Jones matrices.For the dd calibration steps presented in Alg. 1 we use solution intervals of 0.5 − 1 minute.Following van Weeren et al. ( 2016), we add an extra calibration step 1.17, where the visibilities are modeled using the latest available skymodel together with the smoothed Jones-matrices estimated in Step 1.16 and that 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 minutes in ddf-pipeline-v2 giving 11 solution intervals in the 8 hours 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 physical variations.As shown in the right panel of Fig. 6 this method is very efficient in reducing the negative haloes.

Unpolarised flux
Once the estimated dd-Jones matrices and skymodel x ν have been obtained at the highest available spatial resolution following the di/dd-self-calibration steps presented in Alg. 1, additional data products are formed.
Users can adapt the weighting scheme depending on the scientific exploitation they want to make of the interferometric data.This is very much tied to how the calibration and deconvolution algorithms are working, 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 non-zero).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, that can drive, in return, systematic errors in the calibration solutions estimates, be-cause 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 Alg. 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 degrade 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.Also, we store the sub-space deconvolution masks (ssd hereafter, see Tasse et al. 2018)   The di-calibrated visibilities as well as the final skymodels and dd calibration solutions are stored.This allows for additional postprocessing to be made such as better calibration towards a particular point on the sky (van Weeren et al in prep.), and also reimaging at different resolutions if required.

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 psf have to be computed (as there are 16 terms in the quadartic mean of the Mueller matrices).As most of the sources are unpolarised, the leakage terms are properly taken into account in the dd-predict (i.e. the forward mapping from sky and Jones matrices estimates to modeled visibilities).Instead of deconvolving the polarised flux, we grid the IQUV residual data.The polarised flux is directly interpretable when the sources are unresolved.Hence we also generate the following additional products: The output QU cubes are processed using Faraday rotation measure (RM) synthesis (Brentjens & de Bruyn 2005) to find polarised sources and their RM 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 2 and an ability to measure RMs of up to ∼450 rad/m 2 (e.g.O'Sullivan et al. 2020).
The 3 QU cubes are sensitive to the large-scale polarised emission from the Milky Way, while the 20 QU cubes are excellent for finding compact polarised sources.However, detailed studies of the polarisation and RM structure of resolved extragalactic sources will require deconvolution of the Q and U data.The ddf-pipeline-v2 output provides significantly better performance in correcting for the effect of the instrumental polarisation (Fig. 8), 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 polarisation angle calibration for each LoTSS observation, meaning that while the RM values of sources in overlapping fields are consistent, the polarisation angles are not.Therefore, to avoid unnecessary depolarisation for both mosaicing and the deep fields, the polarisation angles between the observations need to be aligned.The simplest way to do this is by choosing a reference angle of a polarised source in a single observation and applying a polarisation angle correction to all other observations to align with this reference angle, as presented in Herrera Ruiz et al. (2020).An alternative approach is to use the diffuse polarised emission that is present in the ∼ 3 QU cubes.
Bright polarised sources are rare in the LoTSS data, with only three sources having a polarised intensity greater than 50 mJy beam −1 in the DR1 HETDEX sky area (Van Eck et al. 2018;O'Sullivan et al. 2018).However, in the fields containing these bright polarised sources the ddf-pipeline-v2 output becomes unreliable for polarised 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 polarised sources are inserted into existing LoTSS uv-datasets.Possible solutions will be tested in future pipeline developments.

Comparison between ddf-pipeline-v1 and ddf-pipeline-v2
Fig. 7 shows the comparison between the final high resolution images produced by ddf-pipeline-v1 and ddf-pipeline-v2 for an 8 hours integration LoTSS pointing (P26Hetdex03).Many pro-  Shimwell et al. 2017a).The colorscale is the same on all panels, and diplayed using an inverse hyperbolic sine function to render both the low level artifacts and some bright sources.cesses are involved in the sky reconstruction from radio interferometric data.Imaging and calibration affect the final synthesized maps and introduce complex and systematic residual errors.It is therefore difficult to find a good and absolute metric to compare the final imaging products.
As discussed in Sec.3.1, the quality of the initial di calibration 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 Fig. 7b and 7c 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) skymodel incompleteness, (ii) a uv-distance cut used during the calibration and (iii) the N normalisation function (see Sec.A for details), that we had introduced for ddf-pipeline-v1 to be robust against the absorption of extended extended emission.As shown in Fig. 7e and 7f, the approach developed in Sec.3.2 and 3.3 to conserve unmodeled extended emission and implemented in ddf-pipeline-v2 does not produce any significant low spatial frequencies systematics.

ddf-pipeline-v2 robustness and performance
As explained above ddf-pipeline-v2 is a high level script interfacing kms and ddfacet.Both of the underlying software packages are efficiently parallelised using a custom version of the Python multiprocessing package for process-level parallelism, and Fig. 9: This pie graph shows the nature and ordering of the different steps of Alg. 1 and how the computing time is distributed across them.The lighter and darker grey areas represent the imaging and calibration steps respectively.The black area are the miscellaneous tasks (additional data products, see Sec. 3.4) that are done once the di and dd self-calibration loops have completed.It has been created from a ddf-pipeline-v2 run on a node equipped with 192 GBytes RAM and 2 Intel Xeon Gold 6130 CPU@2.10GHz,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.using the SharedArray12 module.As explained 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 1,500 hrs every 6 month cycle which corresponds to approximately two 8 hr 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 8 hour pointing is ∼ 5 days (on a node equipped with 192 GBytes RAM and 2 Intel Xeon Gold 6130 CPU@2.10GHz,giving 32 physical compute cores), and takes an extra ∼ 30% of computing time to completion as compared to ddf-pipeline-v1.Hence 10 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 1st 2019 over 1,000 pointings exist in the archive.Over ∼ 1000 pointings and ∼ 12 PB of averaged and compressed LOFAR data (∼ 40 PB uncompressed) have now been processed with ddf-pipeline-v2.

Observations
LoTSS-Deep Fields observations are being carried out over the four northern fields with high-Galactic latitude and the highestquality multi-degree-scale ancillary data across the electromagnetic spectrum: the Boötes field, the Lockman Hole, ELAIS-N1 and the North Ecliptic Pole fields.The ultimate aim of the LoTSS Deep Fields project is to reach noise levels of 10-15 µJy.beam −1 in each of these fields (requiring ∼ 500 hours of integration).The first LoTSS-Deep Fields data release consists of initial observations in three of these fields: Boötes (∼ 80 hrs) and Lockman Hole (∼ 112 hrs) presented in the current paper, and ELAIS-N1 (presented by Sabater et al. 2020, for an integration time of ∼ 170 hrs in paper 2).This first data release also includes an extensive effort of optical/IR cross-matching, which has obtained host galaxy identifications for over 97% of the ∼80,000 radio sources detected within the ∼ 25 deg 2 overlap with the high-quality multi-wavelength data (Kondapally et al. 2020, Paper 3).This is supplemented by high quality photometric redshifts, and characterisation of host galaxy properties (Duncan et al. 2020, Paper 4), and source classification (e.g.star-forming vs AGN: Best et al. 2020, Paper 5).
In order to put the LoTSS-deep observations in a wider context, in this section we briefly describe the multi-wavelength data available on the Boötes and Lockman Hole fields, focusing on the radio coverage (for a more detailed description see Kondapally et al. 2020, Paper 3).
The Boötes pointings 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 Cycle 2 and Cycle 4, with a bandwidth of 48 MHz (see Tab. 1).The total integration time of ∼ 80 hours is spread over 10 scans of 8 hours.

Lockman hole
The Lockman Hole field is also covered by a large variety of multiwavelength data.Lawrence et al. 2007), and with the Submillimetre 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 2 areas (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).Fig. 10 presents an overview of the available radio data on the Lockman Hole.

Image synthesis
The Lockman Hole and Boötes fields data have been both reduced using Alg. 2. In this approach we first build a wide-band di+dd self-calibrated sky model x ν from a single wide band ∼ 8 hours observation using Alg. 1.This model is 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 Alg. 2. This amounts to repeating Steps 1.13 to 1.18 of Alg. 1 on a larger dataset.A comparison between the images synthetised from 8 and 80 hours 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@2.00GHz with 14 physical cores each, Alg. 2 took ∼ 21 days to process the 80 hours of Boötes data.Estimating the noise in radio maps is not straightforward since noise is correlated and non-Gaussian.Also, 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 use the statistics of the min {.} estimator (that returns the minimum value of a given sample).Intuitively, while the I-Stokes image max {.} statistics has 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 13 , the effective standard deviation is derived.
Fig. 14 shows the cumulative distribution of the local noise in the Lockman Hole and Boötes fields maps, reaching 23 and 30 µJy.beam −1 respectively.Taking into account the number of pointings with their respective amount of flagged data, we get total integration times of ∼ 65 and ∼ 88 hours on the Boötes and Lockman Hole fields respectively, giving a theoretical thermal noise difference of a factor ∼ 1.16 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 elevation, and the Galactic noise differences.2018) is different 14 the comparison can only be approximate.In Fig. 12 we compare the images produced by Retana-Montenegro et al. (2018) and by Alg. 2. While the noise difference should be on 13 The cumulative distribution , where n in the number of pixels in a given box.Finding y σ such that F {y σ } = 1/2 given the box size gives us a conversion factor from the minimum estimate to the standard deviation. 14Out 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. 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 Alg. 2 and implemented in ddf-pipeline-v2.

Cataloguing
In order to extract astrophysical information we build a catalogue of radio sources from the images produced by Alg. 2 and the data described in Sec. 4.Even in the apparent flux maps, because of the imperfect calibration and imaging, the LoTSS-deep     The colorscale is the same on all panels, and diplayed using an inverse hyperbolic sine function to render both the low level artifacts and some bright sources.The stripy artifact seen in the zoomout Fig. 13b seems to be produced by the residual deconvolution and calibration errors of a few 10 mJy.beam −1 bright sources that are a few degrees away from the center of the field.

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-frequencyantenna dd-Jones matrices.Due to the high number of free parameters in that system, and to the finite amount of data points in the non-linear 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 both to 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 analysis of the LoTSS-deep catalogues (including the source counts of the Lockman Hole, Boötes and ELAIS-N1 fields) are presented in Mandal et al. (2020).In the future we plan to continue increasing the depth of these fields: data are already in hand, or scheduled, to double the integration time on each field, with a further aim to increase this to 500 hours in each field.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 Ω n and kms-estimated Jones matrices J Ω n at locations Ω n .The function K abstracts the dd calibration step, and takes as arguments the visibilities v, the skymodel x ν , 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 Ω n in which to solve for.The extra functions C, N, and B represent the clustering, normalisation (see text), and bootstrapping steps respectively.

Acknowledgements
Data  The data processing strategy of the LoTSS first data release (DR1) has been extensively described by Shimwell et al. (2019).Since addressing the issues described in Sec. 2 involves making improvements relative to this approach, we give here a brief description of the data reduction strategy in ddf-pipeline-v1 (the various steps are outlined in Alg.0).As discussed in Sec. 2, the calibration and imaging problem is non-convex and ill-posed.Beyond the computational issues, the great difficulty of the calibration of the dde is sky incompleteness, because the dd-C-rime non-linear 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 (Alg.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 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 ∼ g tν g H tν (where A H is the hermitian transpose of matrix A).This amounts to constraining the baselinedependent error to be solely antenna-dependent.This normalization (described by the function symbol N in Alg.0) was able to recover some extended emission otherwise absorbed in the calibration solution.However, as shown in Fig. 7a 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 having a bright 1 Jy source within the primary beam (such as 3C sources), ddf-pipeline-v1 was not able to converge.

Fig. 1 :
Fig. 1: This figure shows the effective noise in the LoTSS-Deep continuum maps as compared to other existing and future surveys.A spectral index of α = −0.7 has been used to convert flux densities to the 1.4 GHz reference frequency.

Algorithm 1 :
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 Ω n and kmsestimated Jones matrices J Ω n at locations Ω n .The function K abstracts the dd calibration step, and takes as arguments the visibilities v, the skymodel x ν , 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 Ω n in which to solve for.The extra functions C, B, and F represent the clustering, bootstrapping and smoothing steps respectively.
3.4 for a discussion of polarisation related data products).The solution intervals δt 0 and δν 0 along time and frequency are determined such that n b ∝ (T/ |x ν | ) 2 Var{n} where n b is the number of points in the δt 0 × δν 0 time-frequency domain, T is the target solution SNR, and Var{n} is the variance of the visibilities' noise (see Mbou Sob et al. in preparation for a justification).
but in each simulation the recovered flux of the completely unmodelled emission exceeded 60%.Examples of the injected and recovered emission are shown in Fig 5.

Fig. 4 :
Fig. 4: This figure shows the amplitude and phase (top/bottom respectively) of a scalar Jones matrix for a given station in a given direction in the example observation.The left panel shows the solution as estimated by the kms solver.The right panel shows the regularised solution, as updated by the F function.The amplitude color scale ranges from 0 to 1.5.

Fig. 5 :
Fig.5: In order to test the robustness of the algorithm described in Sec. 3 and implemented in ddf-pipeline-v2, we have simulated an unmodeled extended emission (left panel).The emission is absorbed by the dd-calibration step (middle), while it can be partially recovered (right panel) by decreasing the effective size of the unknown solutions space (Sec.3.2 and 3.3).
Fig.6: Conserving the unmodeled extended emission while keeping high dynamic range is extremely challenging in the context of dd calibration and imaging.The left panel shows that a faint and unmodeled extended emission (on the level of ∼ 0.7σ here) can be totally absorbed.While regularising 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 (middle panel).The right panel shows that solving the residuals on longer time intervals (Step 1.17) corrects for this issue.
3. High resolution (6 ) Stokes I image in 3 frequency chunks spread over the whole HBA bandwidth (Step 2b.2)

4.
Low resolution (20 ) spectral Stokes QU cubes (480 planes -Step 2b.3) 5. Very low resolution Stokes QU cubes (480 planes -Step 2b.4), by cutting the baselines > 1.6 km, giving an effective resolution of ∼ 3 6.Low resolution (20 ) wide-bandwidth Stokes V image (Step 2b.5) Fig. 7: This figure shows the differences between the maps produced by Alg.0 and Alg. 1 from a typical 8 hour scans (here the P26Hetdex03 pointing in the HETDEX field, seeShimwell et al. 2017a).The colorscale is the same on all panels, and diplayed using an inverse hyperbolic sine function to render both the low level artifacts and some bright sources.

Fig. 8 :
Fig. 8: A plot of the Faraday depth spectrum, or Faraday dispersion function (FDF), for a radio galaxy in both DR1 and DR2 datasets, showing the improvement in the suppression of the instrumental polarisation signal.The blue dashed line shows the FDF from the DR1 data with a strong instrumental polarisation feature near Faraday depths of φ ∼ 0 rad/m 2 , while the orange solid 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.
Fig. 10: The sensitivity of the various deep dedicated surveys covering the Boötes (top) and Lockman Hole fields (bottom) 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 pointings are marked with a red cross, the dashed line corresponding to a source having a spectral index of −0.7.

Fig. 12 (
further discussed in Sec.4.3) and 13 show the central parts of the of these deep LOFAR Boötes and Lockman Hole observations.

4. 3 .
Comparison with deep factor image synthesis The image of the Boötes field based on 55 hours of LOFAR HBA data and presented Retana-Montenegro et al. (2018) reaches an unprecedented noise level image of ∼ 55 µJy.beam −1 at 150 MHz.To achieve such high sensitivity, Retana-Montenegro et al. (2018) have applied third generation calibration and imaging to correct for the dde using the factor package (developped by van Weeren et al. 2016, see Sec. 2 for more detail).Because the set of LOFAR datasets used by Retana-Montenegro et al. (

Fig. 11 :
Fig. 11: This figure shows the restored high resolution image towards the center of the Boötes field for the 8 hours image produced with Alg. 1 (top panel) and the 80 hours image produced with Alg. 2 (bottom panel).Both images are thermal noise limited, with the same colorscale being used on both.

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 Ω n and kms-estimated Jones matrices J Ω n at locations Ω n .The function K abstracts the dd calibration step, and takes as arguments the visibilities v, the skymodel x ν , 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 Ω n 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 × 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 ν /* On n p ×240 LOFAR HBA subbands */ /* DD calibration */ . The sources were detected with a 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 synthesised 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 catalogue are the source position, peak and integrated flux density, source size and orientation, the associated uncertainties, the estimated local rms at the source position, as well as a code describing the type of structure fitted by PyBDSF.As described in Sabater et al. (2020), 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 for Lockman Hole and Boötes fields respectively.These numbers were estimated from the comparison C. Tasse: The LOFAR Two Meter Sky Survey: Deep Fields (a) The central 2 deg 2 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 Alg. 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).

Fig. 12 :
Fig. 12: Comparison between the LOFAR-HBA maps generated at 150 MHz by Retana-Montenegro et al. (2018) and in the current work.The colorscale is the same on all panels, and diplayed using an inverse hyperbolic sine function to render both the low level artifacts and some bright sources.

Fig. 13 :
Fig. 13: This figure shows the central region of the deep LOFAR-HBA maps of the Lockman Hole field generated at 150 MHz.The colorscale is the same on all panels, and diplayed using an inverse hyperbolic sine function to render both the low level artifacts and some bright sources.The stripy artifact seen in the zoomout Fig.13bseems to be produced by the residual deconvolution and calibration errors of a few 10 mJy.beam −1 bright sources that are a few degrees away from the center of the field.

Fig. 14 :
Fig. 14: The cumulative distribution of the local noise estimates in the various maps discussed here.As shown here, we have imaged a larger fraction of LOFAR's HBA primary beam than the image presented in Retana-Montenegro et al. (2018).
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 kern astronomical software package (available at https://kernsuite.info and presented in Molenaar & Smirnov 2018).MB acknowledges support from INAF under PRIN SKA/CTA FORECaST.MB 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.MJJ 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.PNB and JS are grateful for support from the UK STFC via grant ST/R000972/1.MJH acknowledges support from STFC via grant ST/R000905/1.WLW acknowledges support from the ERC Advanced Investigator programme NewClusters 321271.WLW 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).AB acknowledges support from the VIDI research programme with project number 639.042.729, which is financed by the Netherlands Organisation for Scientific Research (NWO).IP acknowledges support from INAF under the SKA/CTA PRIN "FORECaST" and the PRIN MAIN STREAM "SAuROS" projects MB 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.RK acknowledges support from the Science and Technology Facilities Council (STFC) through an STFC studentship.

Table 1 :
Overview of the deep fields pointings used to synthetise the images on the Boötes and Lockman Hole extragalactic fields.Columns f f lag and nMS stand for the fraction of flagged data and number of measurement sets present in the archives.