A&A 674, A119 (2023)

https://doi.org/10.1051/0004-6361/202245124 Astronomy
&
© The Authors 2023
Astrophysics

LOFAR Deep Fields: Probing faint Galactic polarised

emission in ELAIS-N1⋆
Iva Šnidarić1 , Vibor Jelić1 , Maaijke Mevius2 , Michiel Brentjens2 , Ana Erceg1 , Timothy W. Shimwell2,3 ,
Sara Piras4 , Cathy Horellou4 , Jose Sabater5 , Philip N. Best5 , Andrea Bracco6 , Lana Ceraj1 , Marijke Haverkorn7 ,
Shane P. O’Sullivan8 , Luka Turić1 , and Valentina Vacca9

1
Rud̄er Bošković Institute, Bijenička cesta 54, 10000 Zagreb, Croatia
e-mail: isnidar@irb.hr; vibor@irb.hr
2
ASTRON, Netherlands Institute for Radio Astronomy, Oude Hoogeveensedijk 4, 7991 PD Dwingeloo, The Netherlands
3
Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands
4
Department of Space, Earth and Environment, Chalmers University of Technology, Onsala Space Observatory, 43992 Onsala,
Sweden
5
Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK
6
Laboratoire de Physique de l’École Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université de Paris,
75005 Paris, France
7
Department of Astrophysics/IMAPP, Radboud University, PO Box 9010, 6500 GL Nijmegen, The Netherlands
8
School of Physical Sciences and Centre for Astrophysics & Relativity, Dublin City University, Glasnevin D09 W6Y4, Ireland
9
INAF-Osservatorio Astronomico di Cagliari, Via della Scienza 5, 09047 Selargius (CA), Italy

Received 3 October 2022 / Accepted 6 March 2023

ABSTRACT

We present the first deep polarimetric study of Galactic synchrotron emission at low radio frequencies. Our study is based on
21 observations of the European Large Area Infrared Space Observatory Survey-North 1 (ELAIS-N1) field using the Low-Frequency
Array (LOFAR) at frequencies from 114.9 to 177.4 MHz. These data are a part of the LOFAR Two-metre Sky Survey Deep Fields Data
Release 1. We used very low-resolution (4.3′ ) Stokes QU data cubes of this release. We applied rotation measure (RM) synthesis to
decompose the distribution of polarised structures in Faraday depth, and cross-correlation RM synthesis to align different observations
in Faraday depth. We stacked images of about 150 h of the ELAIS-N1 observations to produce the deepest Faraday cube at low radio
frequencies to date, tailored to studies of Galactic synchrotron emission and the intervening magneto-ionic interstellar medium. This
Faraday cube covers ∼36 deg2 of the sky and has a noise of 27 µJy PSF−1 RMSF−1 in polarised √ intensity. This is an improvement in
noise by a factor of approximately the square root of the number of stacked data cubes (∼ 20), as expected, compared to the one in
a single data cube based on five-to-eight-hour observations. We detect a faint component of diffuse polarised emission in the stacked
cube, which was not detected previously. Additionally, we verify the reliability of the ionospheric Faraday rotation corrections esti-
mated from the satellite-based total electron content measurements to be of ∼0.05 rad m−2 . We also demonstrate that diffuse polarised
emission itself can be used to account for the relative ionospheric Faraday rotation corrections with respect to a reference observation.
Key words. ISM: general – ISM: structure – ISM: magnetic fields – radio continuum: ISM – techniques: interferometric –
techniques: polarimetric

1. Introduction The deep fields are selected in regions covered by a wealth

of multi-wavelength data and the first data release includes
The LOw-Frequency ARray (LOFAR, van Haarlem et al. 2013)
the Boötes, Lockman Hole, and European Large Area Infrared
is a radio interferometer utilising a new-generation phased-array
Space Observatory Survey-North 1 (ELAIS-N1) fields. These
design to explore the low-frequency radio sky (10–240 MHz) in
data sets make it possible to probe a new, fainter population
the Northern Hemisphere. This is carried out through a number
of radio sources dominated by star-forming galaxies and radio-
of key science projects, including a series of ongoing LOFAR
quiet active galactic nuclei (Smolčić et al. 2017; Novak et al.
surveys. The LOFAR Two-meter Sky Survey (LoTSS, Shimwell
2018; Kondapally et al. 2021). In addition, they enable stud-
et al. 2017, 2019, 2022) and LOFAR Low Band Antenna Sky Sur-
ies of the under-explored polarised radio source population at
vey (LoLSS, de Gasperin et al. 2021) are the wide-area surveys
sub-milijansky flux densities at low frequencies.
at 120–168 MHz and 42–66 MHz, respectively. These surveys
Herrera Ruiz et al. (2021) did an initial analysis of polarised
are complemented by a few deeper fields, known as the LoTSS-
radio sources in the LoTSS Deep Fields using six 8-h observa-
Deep Fields (Tasse et al. 2021; Sabater et al. 2021) and the
tions of the ELAIS-N1 field. They used low-resolution (20′′ )
LoLSS-Deep Fields (de Gasperin et al. 2021).
polarimetric images (Sabater et al. 2021) and successfully
⋆
The Faraday cube is only available at the CDS via anonymous demonstrated the feasibility of stacking LOFAR data. While they
ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https:// detected three polarised sources in a single observation, this
cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/674/A119 number increased by more than a factor of three for the stacked

A119, page 1 of 15
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0),
which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
 A&A 674, A119 (2023)

data; they reported detection of seven additional sources, yield- rotation corrections is needed, as well as, an assessment of the
ing a surface density of polarised sources of one per 1.6 deg2 . method for LOFAR High Band Antenna (HBA) observations and
This work is being extended through stacking of 19 eight-hour direction-dependent effects.
LoTSS Deep Field observations of ELAIS-N1 re-imaged at The six ELAIS-N1 observations analysed by Herrera Ruiz
higher angular resolution (6′′ ) to further decrease the detection et al. (2021) were corrected for the ionospheric Faraday rotation
threshold and increase the number of detected polarised sources by the satellite-based TEC measurements (Sabater et al. 2021).
and probe the polarised source counts in the sub-mJy regime To check how well they were corrected relatively to each other,
(Piras et al., in prep.). Herrera Ruiz et al. (2021) compared the observed Faraday depth
When stacking, polarisation data need to be first properly of the bright reference source in each observation and found a
calibrated and corrected for the Faraday rotation in the Earth’s relative difference varying from −0.12 to +0.05 rad m−2 . Then
ionosphere (Murray & Hargreaves 1954; Hatanaka 1956). they calculated the difference in the observed polarisation angle,
Ionospheric Faraday rotation is a time- and direction-dependent corrected each observation accordingly, and stacked the data.
propagation effect proportional to the integral along the line of A complementary method to check for a relative alignment
sight (LOS) of the product of the total electron content (TEC) between the observations concerning the Faraday rotation in the
of plasma in the ionosphere and a projection of the geomag- ionosphere is based on using the polarised diffuse Galactic syn-
netic field, Bgeo , to the LOS towards the observed field of view chrotron emission (Lenc et al. 2016; Brentjens 2018). This type
(FoV). It is characterised by the ionospheric rotation measure of emission is ubiquitous at low radio frequencies (e.g. Erceg
(RMion ), which, in the thin-shell model, can be approximated as et al. 2022, and references therein) and allows analysis over
(e.g. Sotomayor-Beltran et al. 2013) a larger portion of the FoV compared to using a single refer-
ence polarised source. Ionospheric Faraday rotation corrections
RMion TECLOS Bgeo,LOS obtained in such a way should improve the accuracy of cor-
= 0.26 , (1) rections and allow the analysis of differential variations across
[rad m ]
−2 [TECU] [G]
the field.
where TECLOS is the total electron content, measured in TEC In this work, we used the polarised diffuse synchrotron emis-
units (1 TECU = 1016 electrons m−2 ), at the ionospheric pierc- sion to study the ionospheric Faraday rotation corrections in 21
ing point of the LOS. A typical RMion is 0.5–2 rad m−2 LOFAR observations of the ELAIS-N1 field. We also stacked
(Sotomayor-Beltran et al. 2013; Jelić et al. 2014, 2015) at moder- very low-resolution images (4.3′ ) to study the faint compo-
ate geographical latitudes during nighttime. Daytime values are nent of the diffuse polarised emission in the ELAIS-N1 field,
higher due to solar irradiation and an increase in TEC. The TEC whose bright component was observed in the commissioning
decreases after the sunset due to recombination of plasma in the phase of the LOFAR (Jelić et al. 2014). The paper is organ-
ionosphere. ised as follows. LOFAR observations and related data products
Given that ionospheric Faraday rotation changes the polar- are described in Sect. 2. Section 3 presents the analysis of the
isation angle θ of the observed emission on timescales smaller ionospheric Faraday rotation corrections. Section 4 describes
than the total integration time of observation, the observed methodology for stacking the very low-resolution images. The
polarised emission may be incoherently added during the syn- final stacked Faraday cube is presented and analysed in Sect. 5.
thesis, resulting in partial, or, in exceptional cases, full depo- The newly detected faint polarised emission is discussed in
larisation. Ionospheric depolarisation effects are mostly relevant Sect. 6. Summary and conclusions are presented in Sect. 7.
at lower radio frequencies, as Faraday rotation is inversely
proportional to a square of the frequency (∆θ ∼ RMion ν−2 ).
At 150 MHz, a change in the ionospheric Faraday rotation of 2. Data and processing
∼0.8 rad m−2 results in a 180◦ rotation of the polarisation vector
and therefore full depolarisation. In this section, we describe the LoTSS-Deep Fields observations
The LOFAR observations are usually corrected for the and the derived data products used in this paper. We also give an
ionospheric Faraday rotation in a direction-independent manner overview of the rotation measure (RM) synthesis technique and
by combining global geomagnetic field models with Global its parameters used to create Faraday cubes.
Navigation Satellite System (GNSS) observations of the iono-
spheric TEC (Sotomayor-Beltran et al. 2013; Mevius 2018). This 2.1. LoTSS-Deep Fields observations and very low-resolution
was first tested on the LOFAR commissioning observations images
of the ELAIS-N1 field (Jelić et al. 2014), and since then, it is
widely used in polarisation studies with LOFAR (e.g. Jelić et al. The ELAIS-N1 data analysed in this paper are part of the LoTSS-
2015; Van Eck et al. 2017; Turić et al. 2021; Erceg et al. 2022). Deep Fields Data Release 1 (Sabater et al. 2021). We used 21 out
Depending on the source of the TEC data, the estimated uncer- of 27 observations, which were of good quality (10 observations
tainty in the calculated ionospheric Faraday rotation is within a from Cycle 2 and 11 observations from Cycle 4, IDs 009–018,
factor of a few of 0.1 rad m−2 at time intervals of 15 min to 2 h. 020–024, 026–028, 030–032 in table 1 in Sabater et al. 2021).
Recently, de Gasperin et al. (2018) showed that LOFAR Low The data were taken with the LOFAR HBA from May 2014 to
Band Antenna (LBA) station-based gain phase can be decom- August 2015 (under project codes LC2_024 and LC4_008), cov-
posed into a few systematic effects related to clock delays and ering the frequency range from 114.9 to 177.4 MHz dived into
ionospheric effects and used directly to obtain independent mea- 320 frequency sub-bands. The observing time of each observa-
surements of the absolute TEC. The LOFAR measured TEC tion was between 5 and 8 h, taken during night-time and symmet-
values are within 10% of the satellite-based measurements and ric around transit. The array was used in the HBA DUAL INNER
have two orders of magnitude better time resolution. This has configuration (van Haarlem et al. 2013). The HBA antennas of
enabled a new, efficient, unified calibration strategy for LOFAR each core station are clustered in two groups of 24 tiles of 16
LBA (de Gasperin et al. 2019). However, further detailed anal- dual-polarised antennas. Each cluster of 24 tiles was then treated
ysis of systematic uncertainties related to ionospheric Faraday as an independent HBA core station. The remote stations have
A119, page 2 of 15
 Šnidarić., I., et al.: A&A proofs, manuscript no. aa45124-22

one group of 48 HBA tiles of 16 dual-polarised antennas. They 5

009 016 024
were reduced to inner 24 tiles, to have the same shape and num- 010 017 026
ber of tiles as the dual core HBA stations. This provided a uni- 4 011 018 027
form general shape of the primary beam over the entirety of the 012 020 028

RMion [rad m 2]
LOFAR stations in the Netherlands. The phase centre of the main 013 021 030
3 014 022 031
target field was at RA 16h 11m 00s and Dec +55◦ 00′ 00′′ (J2000). 015 023 032
Cycle 2 data were taken and pre-processed jointly with the
LOFAR Epoch of Reionisation Key Science Project team in a 2
slightly different way than Cycle 4 data. This created a difference
in frequency configurations of the final data products of the two 1
cycles. Here we give a brief overview of the main processing
steps and relevant differences for each cycle, while details are
provided in Sabater et al. (2021). 0
0 100 200 300 400 500
The pre-processing of the data included averaging in time Time [min]
and frequency. Before averaging, the Cycle 2 data were auto-
matically flagged for radio-frequency interference (RFI) using Fig. 1. Calculated RMion corrections given at 30-min intervals for dif-
AOFlagger (Offringa et al. 2012). The first two and the last two ferent observations using the satellite-based TEC measurements and the
frequency channels were then removed from each 64-channel global geomagnetic field model. The observed decrease of RMion dur-
sub-band to minimise the band-pass effects. The remaining ing each nighttime observation is due to recombination of plasma in
60 channels were averaged to 15 channels per sub-band. The the ionosphere, which happens after the sunset and decreases the TEC
Cycle 4 data were originally averaged by the observatory to throughout night.
16 channels per sub-band, without discarding the channels at
Cycle 2
edges of each sub-band. After that, they were flagged for the 20.0
RFI. The data from both cycles were averaged in time to 2 s.
The direction-independent calibration was done using the 17.5
PREFACTOR pipeline (van Weeren et al. 2016; de Gasperin 15.0
]
1

et al. 2019), which corrects for the polarisation phase offset

NoiseQU mJy PSF

introduced by the station calibration table, the instrumental time 12.5

delay associated with clocks in the remote stations, the ampli- 10.0
tude band-pass, and ionospheric direction-independent delays
[

and Faraday rotation. 7.5

The ionospheric Faraday rotation corrections were done by 5.0
RMextract (Mevius 2018), which combines the satellite-based
TEC measurements and the global geomagnetic field models 2.5
to predict the corrections. The ionosphere was modelled as a 0.0
single-phase screen above the array, taking into account its spa- 120 130 140 150 160 170 180
tial structure. The RMion corrections were calculated for each Frequency (MHz)
LOFAR station separately; however, the model does not allow Cycle 4
for direction-dependent corrections within the FoV. Moreover, 20.0
the model uses a thin-screen approximation. The contribution
from the plasmasphere to the total integrated electron content
17.5
along the LOS can be significant (up to 40% at moderate geo- 15.0
]
1

graphical latitudes, Yizengaw et al. 2008). Since the GNSS data

NoiseQU mJy PSF

include the full integrated electron density, including the plasma-

12.5
spheric contribution, and the magnetic field contribution from 10.0
the higher layers is significantly smaller, the derived RM val-
[

ues using a thin screen model are likely an overestimate of the

7.5
ionospheric Faraday rotation. Figure 1 shows calculated values 5.0
for the LOFAR station CS002, as an example. The curves are
given for different nights as a function of the observing time at 2.5
30-min intervals. The absolute RMion values are between 0.5 and 0.0 120 130 140 150 160 170 180
3 rad m−2 , while their relative variations during observations are
on average 0.9 ± 0.3 rad m−2 .
Frequency (MHz)
The final step of processing included the direction-dependent Fig. 2. Noise in Stokes QU data cubes for observations from Cycle 2
calibration done by the DDF pipeline (Tasse et al. 2021) and (top plot) and Cycle 4 (bottom plot) as a function of frequency. Cycle 4
imaging of the data in full Stokes parameters (I, Q, U, and V). data are much more affected by broad RFI than Cycle 2 data due to
In this work, we used very low (vlow) resolution (4.3′ ) Stokes DABs and DVBs.
QU data cubes (Sabater et al. 2021). They are split in 800 or
640 frequencies of 73.24 kHz or 97.66 kHz width in the case of as a standard deviation in the corner of the primary-beam-
Cycle 2 or 4 data, respectively, due to their different frequency uncorrected image (farthest out of the primary beam), where no
configurations. polarised emission is present. The noise in Stokes Q and U is
Figure 2 shows the noise in Stokes QU data cubes as a func- comparable. Cycle 4 data are much more affected by broad RFI
tion of frequency. The noise at each frequency was calculated (Offringa et al. 2013) at frequencies around 134, 151, 167, and

A119, page 3 of 15
 A&A 674, A119 (2023)

174 MHz than Cycle 2 data. These RFIs are due to man-made 1.0 CYCLE 2
wireless applications such as digital audio or video broadcasts CYCLE 4
(DABs or DVBs). Over the observed frequency range, a typical
noise at frequencies not affected by the RFI is ∼3.3 mJy PSF−1 in 0.8
Cycle 2 data and ∼2.7 mJy PSF−1 in Cycle 4 data. A small differ-
ence between the two cycles arises from their different frequency 0.6

RMSF
configurations and hence frequency channel widths.
0.4
2.2. RM synthesis and Faraday depth cubes
We created Faraday data cubes of the ELAIS-N1 deep field 0.2
observations for our analysis. They were produced by apply-
ing the RM synthesis technique (Burn 1966; Brentjens & de 0.0
Bruyn 2005) to Stokes QU frequency data cubes. This technique 15 10 5 0 5 10 15
decomposes the observed polarised emission by the amount
of Faraday rotation of its polarisation angle, θ, experienced at
[radm 2]

wavelength λ: Fig. 3. RMSF for ten observations in Cycle 2 (dashed) and 11 observa-
2 tions in Cycle 4 (solid line) of the ELAIS-N1 field.
∆θ Φ λ
= . (2)
[rad] [rad m−2 ] [m2 ]
The quantity Φ is called Faraday depth, and it is defined as of each observation. The frequency channels with noise
>7.5 mJy PSF−1 were flagged. The resulting Faraday cubes cov-
Φ B∥ dl
Z d
ne ered Faraday depths from −50 to +50 rad m−2 in 0.25 rad m−2
= 0.81 , (3)
[rad m ] 0 [cm ] [µG] [pc] steps, given the expected Faraday depth range of the observed
−2 −3

where ne is the density of thermal electrons and B∥ is the mag- emission in this field (from −10 to +13 rad m−2 , Jelić et al.
netic field component parallel to the LOS. The integral is taken 2014). The resolution in Faraday depth was δΦ = 0.9 rad m−2
over the path length dl from the source (l = 0) to the observer for all observations. The side lobes of the RMSF in Cycle 4
(l = d). If the magnetic field component is pointing towards the data were higher than in Cycle 2 data (see Fig. 3) due to the
observer, the value of the Faraday depth is positive and vice gaps at frequencies contaminated by the broad RFIs (see Fig. 2).
versa. Equation (3) and the sign convention related to the mag- Because the resolution in Faraday depth is comparable to the
netic field component along the line of sight are in agreement maximum detectable Faraday scale (∆Φscale = 1.1 rad m−2 ), we
with the correct sense of Faraday rotation discussed by Ferrière are only sensitive to Faraday-thin structures (λ2 ∆Φscale ≪ 1) or
et al. (2021). the edges of Faraday-thick structures (λ2 ∆Φscale ≫ 1 Brentjens
For a given location in the sky, the RM synthesis gives us the & de Bruyn 2005).
distribution of the observed polarised emission in Faraday depth. The noise in the Faraday cubes for the different observations
This so-called Faraday spectrum is the Fourier transform of the is given in Table 1. The noise was estimated as the standard
complex polarisation of the observed signal, P(λ2 ) = Q(λ2 ) + deviation of an image given in the polarised intensity√at Faraday
iU(λ2 ), from λ2 - to Φ-space (Brentjens & de Bruyn 2005): depth of −50 rad m−2 and multiplied by a factor of 2. At this
Z +∞ Faraday depth, we do not observe any polarised √ emission, and
1 2
the image is dominated by noise. The factor 2 addresses the
F(Φ) = P(λ2 )P∗ (λ2 ) exp−i2Φλ dλ2 , (4)
W(λ2 ) −∞ Rician distribution of the noise in the polarised intensity, which
roughly corresponds to a normally distributed noise in Stokes
where W(λ2 ) is the non-zero-weighting function, usually taken Q and U (e.g. Brentjens & de Bruyn 2005). A mean value of the
to be 1 at λ2 where measurements are taken and 0 elsewhere. noise in Cycle 2 observations is 91 ± 10 µJy PSF−1 RMSF−1 and
If the RM synthesis is applied over a sky area, we can study in Cycle 4 observations is 121 ± 26 µJy PSF−1 RMSF−1 . Higher
the morphology of the observed polarised emission at different noise in Faraday cubes of Cycle 4 data is due to a larger number
Faraday depths, to perform the so-called Faraday tomography. of frequency channels in this cycle affected by RFI (see Fig. 2).
Characteristics of the λ2 distribution constrain scales in Faraday Observation 014 has the lowest noise among both Cycle 2 and
depth that we can probe when performing the RM synthesis. 4 observations, and observation 021 has the lowest noise among
A resolution in Faraday depth is inversely
√ proportional to the Cycle 4 observations. Hence, we choose the 014 observation as
spectral bandwidth (∆λ2 ) as δΦ ≈ 2 3/∆λ2 and corresponds a reference for Cycle 2. For Cycle 4 we take for consistency the
to the width of the rotation measure spread function (RMSF, same reference observation (024) as in Herrera Ruiz et al. (2021),
Brentjens & de Bruyn 2005). The maximum detectable Faraday which is our second-best observation in terms of the noise in this
scale is inversely proportional to the smallest (λ2min ) measured λ2 cycle. We cannot choose the same reference observation for both
as ∆Φscale ≈ π/λ2min . cycles because of their different frequency configurations.
We used the publicly available code rm-synthesis1 and We used publicly available code rmclean3d from
applied it to Stokes Q and U images, which had comparable RM-Tools3 (Purcell et al. 2020) to deconvolve the Faraday
noise levels (<7.5 mJy PSF−1 )2 in the frequency data cube cubes for the side lobes of the RMSF. The code is based on
1
RM-CLEAN algorithm described in Heald et al. (2009). We
https://github.com/brentjens/rm-synthesis used a threshold of five times the noise in the Faraday cube
2 The noise threshold of 7.5 mJy PSF−1 is estimated based on the noise
during the RM-CLEAN process.
characteristics in Stokes QU datacubes of Cycle 2 observations. It cor-
responds to the mean value of it plus six times its variations measured
by the standard deviation at frequencies not affected by the RFI. 3 https://github.com/CIRADA-Tools/RM-Tools

A119, page 4 of 15
 Šnidarić., I., et al.: A&A proofs, manuscript no. aa45124-22

58° 3.5 58° 8

3.0 6
4
56° 2.5 56°

[mJy PSF 1 RMSF 1]

2
2.0

[rad m 2]
(J2000)
(J2000)

0
1.5
54° 54° 2
1.0
4
0.5 6
52° 52°
0.0 8
248° 246° 244° 242° 240° 238° 248° 246° 244° 242° 240° 238°
(J2000) (J2000)
Fig. 4. Image of the highest peak of the Faraday depth spectrum in the polarised intensity (left) and a corresponding image of a Faraday depth
of the highest peak (right) for the observation with the lowest noise (014, the reference observation for Cycle 2). The blue circle in the left image
marks a randomly chosen location for which a Faraday spectrum is presented in Fig. 5.

Table 1. Calculated noise in the Faraday cubes given for different obser- Noise in a Faraday cube of the commissioning observa-
vations and their relative shift in Faraday depth (∆Φshift ) with respect to tion (a single 8h synthesis) was 300 µJy PSF−1 RMSF−1 (Jelić
the reference observation (calculated in Sect. 3). et al. 2014). This is around 3.6 times higher than the noise
in the individual Faraday cubes presented in this work. The
ID Cycle Noise ∆Φshift difference arises from the limited available frequency band-
(µJy PSF−1 RMSF−1 ) (rad m−2 ) width during the commissioning phase of LOFAR and the use
of a simpler calibration strategy that addressed only direction-
009 2 106 −0.159 ± 0.007 independent effects.
010 2 94 −0.065 ± 0.005 The commissioning observation of the ELAIS-N1 field
011 2 88 >1 revealed polarised diffuse emission over a wide range of Fara-
012 2 84 0.053 ± 0.007 day depths ranging from −10 to +13 rad m−2 (Jelić et al. 2014)
013 2 88 0.002 ± 0.006 given a resolution of 1.75 rad m−2 in Faraday depth. The most
014 2 82 ‘reference’ prominent features of that emission are seen in the left image of
015 2 84 0.010 ± 0.005 Fig. 7 in Jelić et al. (2014), showing the highest peak value of
016 2 86 0.003 ± 0.005 the Faraday depth spectrum at each pixel (RA, Dec). The mean
017 2 85 0.010 ± 0.006 surface brightness of that emission is 2.6 mJy PSF−1 RMSF−1 .
018 2 112 0.052 ± 0.005 The same figure also shows the Faraday depth of each peak in an
020 4 107 −0.033 ± 0.005 image presented on the right.
021 4 91 −0.021 ± 0.004 We constructed the same images for the observations anal-
022 4 133 0.020 ± 0.004 ysed in this work. The images for the observation that has the
023 4 111 0.011 ± 0.004 lowest noise level (014) are presented in Fig. 4 as an example.
024 4 102 ‘reference’ Images for all other observations are very similar to these. The
026 4 112 0.033 ± 0.005 observed diffuse emission in the left image of Fig. 4 shows mor-
027 4 112 0.017 ± 0.004 phological similarity with the one detected in the commissioning
028 4 179 0.045 ± 0.006 observation (see left image in Fig. 7 in Jelić et al. 2014). The
030 4 164 0.114 ± 0.005 observed morphological features appear much sharper despite
031 4 103 0.080 ± 0.004 comparable angular resolution in both observations. This is due
032 4 111 0.101 ± 0.004 to almost two times better resolution in Faraday depth than in
the commissioning observation. As a consequence, the observed
Notes. An ID of each observation corresponds to the one given in emission suffers less from depolarisation, as is the case, for
Table 1 in Sabater et al. (2021). example, for a filamentary structure oriented north-south in the
central part of the image. The filament is depolarised in the
commissioning observation, while it is visible in observations
2.3. Comparison with a previous LOFAR commissioning
presented in work. Due to a better signal-to-noise ratio, there
observation
is also more emission visible towards the edges of the image,
where the emission is attenuated by the LOFAR primary beam.
The ELAIS-N1 field was observed previously with LOFAR The mean surface brightness of the observed emission in the
during its commissioning phase (Jelić et al. 2014). That obser- central part of the image is 3.0 mJy PSF−1 RMSF−1 , which
vation was done in a limited frequency range from 138 MHz to is a bit brighter than in the commissioning observation. The
185 MHz. Here we make a comparison between that observa- emission appears in a range of Faraday depths from −16 to
tion and observations used in this work. The comparison is done +14 rad m−2 , starting at slightly smaller and ending at slightly
using Faraday cubes in the polarised intensity. larger Faraday depths than in the commissioning observation.

A119, page 5 of 15
 A&A 674, A119 (2023)

Fig. 5. Example of a Faraday spectrum given in the polarised intensity for the 009 observation (blue line) and the reference (014) observation at
a randomly chosen location marked with a blue circle in Fig. 4 (RA 242◦ 18′ 03.60′′ and Dec 56◦ 08′ 16.80′′ ), given on the left panel. Calculated
modulus of the complex cross-correlation function |ζ| (red line) for the given Faraday spectra, fitted with a Gaussian (blue dashed line) to estimate
the misalignment between the two observations (black vertical line) at this specific location, given on the right panel.

Further discussion on characteristics of the observed emission

are in Sect. 5.2.
| < > | [ Jy PSF 1 RMSF 1]
5

3. Analysis of a relative shift in Faraday depth 4

between different observations
Data analysed in this work were corrected for the ionospheric 3
Faraday rotation by the satellite-based TEC measurements (see
Sect. 2.1). To check how well the data are corrected, we make a 2
relative comparison between each observation and the reference
observation by cross-correlating Faraday cubes. 1
Instead of explicitly cross-correlating Faraday cubes of two
observations (a and b) as in Jelić et al. (2015), it is compu- 0
4 2 0 2 4
tationally more efficient to use the Fourier transform’s cross-
correlation property. We calculated the cross-correlation func- [radm 2]
tion by effectively performing RM-synthesis on Pa (λ2 )P∗b (λ2 ), Fig. 6. Calculated modulus of the averaged complex cross-correlation
as proposed by Brentjens (2018) and implemented in the above- function for observations 009 and 014 (solid magenta line) and vari-
mentioned publicly available code rm-synthesis. Using this ations of the cross-correlation function across the Fov as measured
code, we evaluate the cross-correlation function for Faraday by a standard deviation (dashed cyan line). A misalignment between
depths between −5 and +5 rad m−2 in 0.01 rad m−2 steps. two observations is determined by fitting a Gaussian to the peak (solid
We expect the relative shift between the observations to be yellow line).
≲1 rad m−2 .
We illustrate the applied method in Fig. 5 by giving examples
of Faraday spectra in the polarised intensity for two obser- Figure 6 shows the calculated modulus of the averaged
vations (009 and 014; left image) and the modulus of their complex cross-correlation function for observations 009 and
evaluated complex cross-correlation function ζ (right image). 014 (magenta solid line). The same figure also gives variations
Faraday spectra are taken for a random (RA, Dec) pixel in the of the cross-correlation function across the FoV as measured by
cube (marked with a blue circle in Fig. 4), where the observed a standard deviation (cyan dashed line). We then fit a Gaussian to
emission is relatively bright. To find a Faraday depth of the cross- the peak and find that the 009 observation is shifted by −0.159 ±
correlation function’s peak, we fitted a Gaussian to the peak. 0.007 rad m−2 with respect to the reference observation.
At this specific position in the image, the Faraday spectrum of Calculated relative shifts in Faraday depth for all other obser-
the 009 observation is shifted by −0.047 ± 0.006 rad m−2 with vations are given in Table 1. There is no significant difference in
respect to a reference (014) observation. misalignment between the two cycles. The observations are on
To find a common shift in Faraday depth across the FoV, we average misaligned by ±(0.046 ± 0.042) rad m−2 with respect to
averaged the complex cross-correlation functions for all pixels the reference observation. The only exception is 011 observa-
(RA, Dec), where the observed emission in a reference observa- tion, which shows a misalignment larger than 1 rad m−2 and is
tion has the highest peak value of the Faraday depth spectrum at analysed in detail in Appendix A.
least ten times larger than the noise (≥82 mJy PSF−1 RMSF−1 ). The estimated misalignments are comparable to the one
This will improve the signal-to-noise ratio and therefore the loca- found in the analysis of five LOFAR observations of the 3C 196
tion of the main peak of the cross-correlation function. When field (0.1 ± 0.08 rad m−2 ; Jelić et al. 2015). This verifies the
averaging, we assume that the variation of shifts in Faraday depth reliability of ionospheric Faraday rotation corrections estimated
across the FoV is much smaller than the width of the main peak using the satellite-based TEC measurements. The related uncer-
of the RMSF (≪0.9 rad m−2 ). tainties are mostly connected to daily systematic biases in the
A119, page 6 of 15
 Šnidarić., I., et al.: A&A proofs, manuscript no. aa45124-22

TEC measurements of ∼1 TEC unit, translating to an error in the Cycle2

ionospheric rotation measure of ∼0.1 rad m−2 . The misalignment 10
for observations 020, 027, 028, 030, and 031 are in agreement 8
with the one estimated by Herrera Ruiz et al. (2021) within the 6

Number of observations
errors. Herrera Ruiz et al. (2021) based their analysis using a 4
single Faraday spectrum at the location of the peak pixel of the 2
0
reference polarised source, while we used all pixels that show Cycle4
bright polarised diffuse emission. Therefore, estimated errors are 10
∼5 times smaller in our work than in their work. 8
6
4. Stacking very low-resolution data 4
2
To stack images of different observations together, we first need 0
to ‘de-rotate’ the observed polarisation angle of each observation 120 130 140 150 160 170 180
Frequency (MHz)
by its estimated shift with respect to the reference observation
(∆Φshift , see Table 1). We multiplied the complex polarisation Fig. 7. Number of observations per frequency channel used in the
Pi (λ2 ) = Qi (λ2 ) + iUi (λ2 ) given at each wavelength (frequency) stacked data cube from Cycle 2 (upper plot) and Cycle 4 (lower plot)
2
by exp−i2∆Φshift λ : data. Although the data cover the same frequency range, they have
different frequency configurations and hence a different number of fre-
2 quency channels (see Sect. 2.1).
P̃i (λ2 ) = Q̃i (λ2 ) + iŨi (λ2 ) = Qi (λ2 ) + iUi (λ2 ) exp−i2∆Φshift λ .
 

(5) Sect. 3. The resulting cross-correlation functions were then

averaged, and their common peak was fitted with a Gaussian.
This way, the correction is applied to the whole Faraday spec- The two cubes are aligned in Faraday depth within the error of
trum simultaneously. the fit and can be combined directly to the final Faraday cube.
We then stack all corrected images of each observing We combined the stacked Faraday cubes of two cycles by
cycle by calculating the weighted average at each wavelength calculating the weighted average in Faraday depth:
(frequency):
Φ Φ
i Qi (Φ)wi i U i (Φ)w
P P
Q̃i (λ2 )wiQ̃ (λ2 ) Ũi (λ2 )wŨ 2
i (λ )
Pcombined (Φ) = +i P Φ i , (7)
P P
P Φ
i wi i wi
2 i i
Pcombined (λ ) = +i , (6)
i wi (λ )
P Ũ 2
i wi (λ )
P Q̃ 2
where wΦ i is a Faraday depth independent weight for each
where wiQ̃,Ũ (λ2 )
is a wavelength- (frequency) dependent weight Faraday cube defined as the inverse of the variance of the noise in
for each observation defined as the inverse of the variance of Stokes Q and U Faraday cubes. The noise in the Faraday cube is
the noise in Stokes Q and U images. We recall that the noise in estimated as the standard deviation of an image given in Stokes Q
Stokes Q and U were comparable and were calculated in the cor- and U at −50 rad m−2 . This is the Faraday depth, where we do
ner of each image where the polarised emission was not present. not observe any polarised emission and the image is dominated
We are not able to stack the data from two cycles directly because by noise.
of their different frequency channel widths (see Sect. 2.1). They
are combined at a later stage in Faraday depth. Figure 7 shows a 5. The final stacked Faraday cube
number of images per frequency channel used in the final stacked
The final stacked Faraday cube combines images of 20 ELAIS-
data cube for Cycle 2 and Cycle 4. There are on average nine
N1 LoTSS-Deep Fields observations, ∼150 h of data in total.
images added per frequency channel in Cycle 2 and 11 images in
The cube covers Faraday depths from −50 to +50 rad m−2
Cycle 4.
in 0.25 rad m−2 steps. The resolution in Faraday depth is
Images of observation 011 are not used for the stacked data
0.9 rad m−2 , as defined by the resolution of the stacked
cube of Cycle 2 because of their relatively poor quality compared
Faraday cubes of the two observing cycles.
√ The final image noise
to the images of all other observations. This choice does not have
any significant impact on the final result. Once the data of each is 27 µJy PSF−1 RMSF−1 , which is ∼ 20 smaller than the mean
observing cycle were stacked, we applied the RM synthesis. value of noise in Faraday cubes of every individual observation
The noise in the stacked Faraday cube of Cycle 2 in the two cycles, as expected. In the following two subsections,
data is 32 µJy PSF−1 RMSF−1 and of Cycle 4 data is we first cross-check if we detect the radio sources presented in
40 µJy PSF−1 RMSF−1 . In both cases, this is ∼3 times the catalogue of Herrera Ruiz et al. (2021), and then we analyse
less than the mean value of noise in Faraday cubes and discuss the observed diffuse Galactic polarised emission.
of individual observations (91 ± 10 µJy PSF−1 RMSF−1 and 5.1. Cross-checking the detection of the radio sources
121 ± 26 µJy PSF−1 RMSF−1 , respectively). The noise in the
stacked Faraday cube is reduced by the square root of the number We used the catalogue of the polarised sources provided by
of observations that are stacked, as expected. Herrera Ruiz et al. (2021) to check how many of them we detect
We calculated the cross-correlation between the stacked in our final stacked Faraday cube. The purpose of this compari-
Faraday cubes of two cycles as a function of the displacement in son is only to verify our stacking method on very low-resolution
Faraday depth to check for their alignment. We only considered data. It is not meant to provide an in-depth analysis of the
Faraday spectra that have a peak flux in the polarised intensity polarised sources. This is done in Herrera Ruiz et al. (2021) and
≥82 mJy PSF−1 RMSF−1 , the same limit as the one used in in a follow-up work using the high-resolution data (20′′ and 6′′ ,

A119, page 7 of 15
 A&A 674, A119 (2023)

= 18.43 rad m 2 = 18.43 rad m 2

58° 0.25 58° 0.25

0.20 0.20

]
56° 56°
Declination (J2000)

Declination (J2000)
1

1
mJy PSF 1 RMSF

mJy PSF 1 RMSF

0.15 0.15

54° 0.10 54° 0.10

0.05 0.05

[
52° 52°
0.00 0.00
248° 246° 244° 242° 240° 238° 248° 246° 244° 242° 240° 238°
Right Ascension (J2000) Right Ascension (J2000)
= 4.75 rad m 2 = 4.75 rad m 2
58° 0.8 58° 0.8
0.7 0.7
0.6 0.6
]

]
56° 56°
Declination (J2000)

Declination (J2000)
1

1
0.5 0.5
mJy PSF 1 RMSF

mJy PSF 1 RMSF

0.4 0.4
54° 0.3 54° 0.3
0.2 0.2
[

[
0.1 0.1
52° 52°
0.0 0.0
248° 246° 244° 242° 240° 238° 248° 246° 244° 242° 240° 238°
Right Ascension (J2000) Right Ascension (J2000)

Fig. 8. Example of a successful (upper panels) and an unsuccessful detection of a polarised source (lower panels) in the presented Faraday cubes
(sources with ID 10 and 07 in Herrera Ruiz et al. 2021, respectively). Polarised intensity images in the reference (014, left images) and in the final
stacked Faraday cube (middle images) are given at the closest available Faraday depth, such as that of the source. The location of the source in
each image is marked with the red circle, while the corresponding Faraday spectra are given in plots on the right. A reported Faraday depth of the
sources by Herrera Ruiz et al. (2021) are marked with vertical red lines.

respectively), which are better suited for such analysis than the 5.2. Faint diffuse Galactic polarised emission
very low-resolution data (4.3′ ) used in this work.
We extracted the Faraday spectra and inspect the images We detect diffuse polarised emission in the final stacked Faraday
in our final polarised intensity cube at locations of polarised cube over a range of Faraday depths from –16 up to +18 rad m−2
sources provided in the catalogues. We have a clear detection (see Appendix B). Its brightest and prominent morphological
of nine out of ten radio sources from Herrera Ruiz et al. (2021, features were already detected by Jelić et al. (2014), but over a
Table 2, ID 01–06 and 08–10), while one of them (ID 07) is smaller Faraday depth range, ranging from –10 to +13 rad m−2 ,
difficult to identify due to the presence of the diffuse polarised and with a poorer resolution of 1.75 rad m−2 . Here we give a
emission in our Faraday cube. Two examples are given in Fig. 8 description of all morphological features observed in our final
for sources with IDs 10 and 07. In the first example, the source stacked cube.
is not contaminated by diffuse emission. There is a clear sig- From −16 to −4 rad m−2 there is a northwest to the south-
nature of it in the Faraday spectrum of the stacked data. This east gradient of emission. It starts as a small-scale feature in the
source is, however, difficult to detect in the reference observa- northwest part of the image, and then it grows diagonally across
tion due to a poorer signal-to-noise ratio than in the stacked the centre of the image to an extended northeast-southwest struc-
data. In the second example, we don’t find the signature of the ture. Its mean surface brightness is 3.1 µJy PSF−1 RMSF−1 .
source, either in the stacked data or in the reference observation, From −4 to −0.5 rad m−2 there is diffuse emission whose mor-
due to contamination by diffuse polarised emission that dom- phology is more patchy, but it spreads over the full FoV. It has
inates the image and the Faraday spectrum at the location of a mean surface brightness of 3.5 µJy PSF−1 RMSF−1 . A con-
the source. spicuous, stripy morphological pattern of diffuse emission with
The rotation measures of successfully detected sources in north-to-south orientation dominates in the eastern part of the
our final cube are in agreement with the values provided in the image from +0.5 up to +4 rad m−2 . Its mean surface bright-
catalogue, taking into account a resolution in Faraday depth of ness reaches 4.3 µJy PSF−1 RMSF−1 . Towards higher Faraday
0.9 rad m−2 and a difference in angular resolution of the used depths, structures become very patchy, emission gets fainter, and
data. The polarised radio source catalogue is based on high- then it disappears completely at +18 rad m−2 . The mean surface
resolution LoTSS data (20′′ ), while in our work we used very brightness of this faint emission is 0.4 µJy PSF−1 RMSF−1 .
low-resolution LoTSS data (4.3′ ). Therefore, morphologies of We constructed Faraday moments to make a comparison
polarised sources are mostly not resolved in our data. If a source between the observed diffuse emission in the final stacked cube
was unresolved in our data, while in reality it has, for example, and the reference (014) observation. Faraday moments provide
two lobes (see a source with ID 07, Fig. 7 in Herrera Ruiz et al. a statistical description of Faraday tomographic cubes, as intro-
2021) whose RMs do not differ more than a resolution of the duced by Dickey et al. (2019). The zeroth Faraday moment, M0 ,
data in Faraday depth, we observed its rotation measure as an is the polarised intensity PI(Φ) integrated over the full Faraday
averaged value of the two lobes and additionally weighted by depth range, given in units of mJy PSF−1 RMSF−1 rad m−2 . It
their relative brightness. gives the total polarised brightness of the emission in the Faraday

A119, page 8 of 15
 Šnidarić., I., et al.: A&A proofs, manuscript no. aa45124-22

Fig. 9. Moments of reference (014),

√ (upper images) and final stacked Faraday cube (lower images). The left images give M0 , the middle ones show
M1 , and the right images give M2 .

cube. The first Faraday moment, M1 , is the polarised intensity emission, which is the same in both cases. However, mean val-
weighted mean of Faraday spectra in units of rad m−2 . It mea- ues of the M1 are −0.85 rad m−2 and −0.65 rad m−2 for reference
sures a mean Faraday depth at which the brightest emission and stacked cubes, respectively, indicating that there is on aver-
is observed. Finally, the second Faraday moment, M2 , is the age more emission at positive Faraday
√ depths in the stacked cube
intensity-weighted variance of Faraday spectra, whose square than in the reference cube. The M2 shows the most noticeable
root gives the spread of the spectrum in units of rad m−2 . Its differences. The measured spread in Faraday depth is on average
square root measures a range of Faraday depths over which 42% larger in the stacked than the reference cube. This is again
the brightest emission is observed. The Faraday moments are due to faint emission at larger Faraday depths, which does not
defined as contribute to the second moment of the reference observation.
n
Examples of the faint emission, which is only clearly
X detected in the stacked cube, are shown in Fig. 10. The images
M0 = PIi · ∆Φ, (8) are given at Faraday depths of +14.5 (upper images) and
i=1
+16.25 rad m−2 (lower images). The images in the first panels
PIi · Φi are shown for the reference cube, while in the second panels for
Pn
M1 = i=1
Pn , (9) the final stacked cube. The third panels show the corresponding
i=1 PIi
Faraday spectra at a location of the red circle in the images. The
and brightness of the faint emission is comparable to the noise in the
reference cube and therefore is not detected there.
PIi · (Φi − M1 )2
Pn
M2 = i=1
Pn , (10)
i=1 PIi 6. Discussion on the faint polarised emission
where ∆Φ is a step in Faraday depth. The Faraday moments are newly detected
calculated only for emission whose brightness is larger than a The diffuse polarised emission detected in the final
defined threshold to exclude noise-dominated areas in the data. stacked Faraday cube has the mean polarised intensity of
We used a threshold of mP + 5σP , where mP is the polarised 10 mJy PSF−1 RMSF−1 rad m−2 , as measured in the central
intensity bias and σP is noise in the polarised intensity. region of the M0 (see bottom left image in Fig. 9). This translates
Figure 9 shows the calculated Faraday moments, both for to a mean brightness temperature of ∼9.5 K4 .
the reference (upper images) and the final stacked Faraday cube We recalculate the M0 of the final stacked Faraday cube
(lower images). There is around 15% more integrated emission by restricting it to Faraday depths ≥ + 13 rad m−2 to estimate
in the stacked cube than in the reference observation, as mea-
sured by the M0 . This is due to a better signal-to-noise ratio in 4 The intensity of 1 mJy PSF−1 RMSF−1 rad m−2 corresponds to a
the stacked cube than the reference cube and the contribution of brightness temperature of ∼0.95 K at 144 MHz, a frequency that
the detected faint emission to the M0 . The first Faraday moments corresponds to the weighted average of the observed λ2 used in RM
do not differ much, as they are mostly driven by the brightest synthesis.

A119, page 9 of 15
 A&A 674, A119 (2023)

= + 14.5 rad m 2 = + 14.5 rad m 2

58° 0.30 58° 0.30

0.25 0.25

]
56° 56°
Declination (J2000)

Declination (J2000)
0.20 0.20

1
mJy PSF 1 RMSF

mJy PSF 1 RMSF

0.15 0.15
54° 0.10 54° 0.10

[
0.05 0.05
52° 52°
0.00 0.00
248° 246° 244° 242° 240° 238° 248° 246° 244° 242° 240° 238°
Right Ascension (J2000) Right Ascension (J2000)
= + 16.25 rad m 2 = + 16.25 rad m 2
58° 0.30 58° 0.30

0.25 0.25
]

]
56° 56°
Declination (J2000)

Declination (J2000)

0.20 0.20
1

1
mJy PSF 1 RMSF

mJy PSF 1 RMSF

0.15 0.15
54° 0.10 54° 0.10
[

[
0.05 0.05
52° 52°
0.00 0.00
248° 246° 244° 242° 240° 238° 248° 246° 244° 242° 240° 238°
Right Ascension (J2000) Right Ascension (J2000)

Fig. 10. Examples of faint Galactic polarised emission, which is only clearly detected in the final stacked Faraday cube. The images are given in
the polarised intensity at Faraday depths of +14.5 (upper images) and +16.25 rad m−2 (lower images) for the reference (left images) and the final
stacked cube (middle images). The corresponding Faraday spectra at a location of the red circle in the images are given in plots on the right.

the mean brightness temperature of the newly detected faint high-frequency polarisation data (DRAO GMIMS, Dickey et al.
emission at higher Faraday depths. We get its mean polarised 2019) and to the Galactic Faraday Sky map (Hutschenreuter
intensity to be of 0.5 mJy PSF−1 RMSF−1 rad m−2 , which is et al. 2022). The latter compliments the low- and high-frequency
∼0.475 K. Although this faint emission is not contributing more observations, as it represents the total RM yielded from
than ∼5% to the total observed polarised emission, its rele- the Galaxy. It is constructed using the observed RM of a
vance comes from the fact that it is present at Faraday depths large sample of extragalactic polarised sources, including
at which the emission was not observed before. It increases
√ the the one in the LoTSS polarised source catalogue (O’Sullivan
range of Faraday depths, usually characterised by M2 , over et al. 2023).
which the emission is detected in this field with LOFAR. This is Erceg et al. (2022) found a correlation between the Galactic
especially important for the interpretation of the LOFAR obser- Faraday Sky map and the LOFAR first Faraday moment image.
vations regarding an extent of the probed volume along the However, the ratio of the two cannot be explained by a simple
LOS and underlying distribution of synchrotron-emitting and model of a Burn slab (Burn 1966), which seems to be appli-
Faraday-rotating regions. cable to the high-frequency data (Ordog et al. 2019). A Burn
Depolarisation effects associated with Faraday rotation are slab assumes a mixture of uniform synchrotron-emitting and
significant at low radio frequencies. Only a few percent of the Faraday-rotating regions along the LOS and predicts a ratio of
intrinsically polarised synchrotron emission is observed with two between the modelled total Galactic RM and the observed
the LOFAR (Jelić et al. 2014, 2015; Van Eck et al. 2017, 2019; polarised emission. The observed LOFAR Faraday spectra are
Turić et al. 2021). The questions that arise relate to where more complex to understand, highlighting the high level of com-
along the LOS does depolarisation happen and from where plexity of the LOS distribution of synchrotron emission and
does the observed emission originate. The idea is that we Faraday rotation.
observe mostly close-by emission, while far-away emission is We compare the M1 of our stacked Faraday cube with
depolarised in the magneto-ionic medium on the way to us. the Galactic Faraday Sky map (Hutschenreuter et al. 2022).
However, determining this from the LOFAR observations alone Figure 11 shows a cut-out of this map in the area of ELAIS-
is very difficult. A Faraday depth is not necessarily a good N1 field. A visual comparison with the bottom middle panel
proxy for the distance. We need to take into consideration the of Fig. 9 shows that the northwest-to-southeast gradient in the
full-complexity of the magnetic fields, its possible reversals, first moment is also present in the Galactic Faraday sky map.
and the multi-phase nature of the interstellar medium. This is The values are more negative around the northwest corner of
challenging, but it has been attempted recently in a number the image; then, diagonally towards the centre of the image they
of the multi-tracer and -frequency studies of the LOFAR verge towards zero, and then they increase to more positive val-
observations (Zaroubi et al. 2015; Van Eck et al. 2017; Jelić ues towards the southeast corner of the image. This gradient
et al. 2018; Bracco et al. 2020; Turić et al. 2021; Erceg et al. implies a bending magnetic field in a southeast to a northwest
2022) and by using the magneto-hydrodynamical simulations direction. The magnetic field mostly points towards us in the
(Bracco et al. 2022). For example, Erceg et al. (2022) compared southeast corner of the image, in the central part of the image
the Faraday moments of the LOFAR observations of around it is mostly in the plane of the sky, and then in the northwest
3100 square degrees in the high-latitude outer Galaxy to the corner it points away from us.

A119, page 10 of 15
 Šnidarić., I., et al.: A&A proofs, manuscript no. aa45124-22

of 27 µJy PSF−1 RMSF √ in polarised intensity, which is

−1

an improvement of ∼ 20 in comparison with noise in an

individual five-to-eight-hour observation;
3. The rotation measures of successfully detected polarised
sources in our final Faraday cube are in agreement with
the values provided in the catalogue by Herrera Ruiz et al.
(2021), which is based on the higher angular resolution
LoTSS Deep Field images;
4. We detected a faint component of diffuse polarised emission
in the stacked cube at high Faraday depths, ranging from
+13 to +17 rad m−2 , which was not detected in the com-
missioning observation by Jelić et al. (2014). The brightness
temperature of this emission is ∼475 mK, which is almost
an order of magnitude fainter than the brightest emission
observed in this field;
5. The observed northwest to southeast gradient of emission
in Faraday depth we associate with a bending magnetic field
across the FoV. It is probably connected to the radio Loop III,
Fig. 11. Total Galactic RM in the area of the ELAIS-N1 field, as the ELAIS-N1 field is just at the edge of it.
extracted from the publicly available the Galactic Faraday Sky map by The presented stacking technique provides a valuable tool and
Hutschenreuter et al. (2022). gives perspective for the future deep polarimetric studies of
Galactic synchrotron emission at low radio frequencies with
LOFAR, the Square Kilometre Array, and its other precursors.
Comparable negative values towards the northwest corner of
For example, it can be applied to other LoTSS Deep Fields
the image between the two maps means that we are probing the
(Lockman Hole and Boötes), as well as to the Great Observa-
same volume of magneto-ionic medium. On the contrary, larger
tories Origins Deep Survey-North (GOODS-N) field. It is also
positive values of Galactic Faraday Sky towards the southeast
of importance for cosmological studies, where the polarised
corner of the image than in the first moment means that we are
emission of the Milky Way is the foreground contaminant. For
not probing the same volume. Emission that comes from fur-
instance, if one wants to measure magnetic field properties in
ther away is either depolarised in LOFAR observations or it is
the cosmic web (e.g. Carretti et al. 2022), it is crucial to have
Faraday thick, as discussed by Erceg et al. (2022) in a broader
an independent measurement of the Galactic contribution to
discussion of Faraday moments in the high-latitude outer Galaxy.
the total Faraday rotation observed toward extragalactic sources.
Moreover, from the same work (see Figs. 5, 8, 9, and 10 in Erceg
Moreover, successful extraction of the cosmological signal from
et al. 2022), it is clear that the ELAIS-N1 field is just at the edge
the cosmic dawn and epoch of reionisation also relies on good
of a region associated with the polarised emission from the radio
knowledge of the foreground emission, including the Galac-
Loop III. The observed gradient is perpendicular to the shape of
tic polarised emission (e.g. Jelić et al. 2010; Asad et al. 2015;
the loop, and it is probably associated with it. The same is true for
Spinelli et al. 2019).
the emission observed in the surrounding area of the ELAIS-N1
field. Acknowledgements. We thank the anonymous referee for their constructive com-
Faraday depths at which we observe the diffuse emission in ments. V.J., A.E., L.C. and L.T. acknowledge support by the Croatian Science
the ELAIS-N1 field are comparable or smaller than the total Foundation for a project IP-2018-01-2889 (LowFreqCRO) and additionally L.T.
Galactic RM in the same area. Therefore, the observed diffuse and V.J. for the project DOK-2018-09-9169. P.N.B. is grateful for support from
the UK STFC via grant ST/V000594/1. M.H. acknowledges funding from the
emission is probably Galactic, including its newly detected faint European Research Council (ERC) under the European Union’s Horizon 2020
component. research and innovation program (grant agreement No 772663). VV acknowl-
edges support from Istituto Nazionale di Astrofisica (INAF) mainstream project
“Galaxy Clusters Science with LOFAR”, 1.05.01.86.05. This paper is based on
7. Summary and conclusions data obtained with the International LOFAR Telescope (ILT) under project codes
LC0_019, LC2_024 and LC4_008. LOFAR (van Haarlem et al. 2013) is the Low
We used 21 LOFAR HBA observations of the ELAIS-N1 field Frequency Array designed and constructed by ASTRON. It has observing, data
(about 150 h of data) to conduct the deepest polarimetric study of processing, and data storage facilities in several countries, that are owned by
Galactic synchrotron emission at low radio frequencies to date. various parties (each with their own funding sources), and that are collectively
operated by the ILT foundation under a joint scientific policy. The ILT resources
The analysis was performed on very low-resolution (4.3′ ) Stokes have benefited from the following recent major funding sources: CNRS-INSU,
QU data cubes produced as part of the LoTSS-Deep Fields Data Observatoire de Paris and Université d’Orléans, France; BMBF, MIWF-NRW,
Release 1 (Sabater et al. 2021). A stacking technique was devel- MPG, Germany; Science Foundation Ireland (SFI), Department of Business,
oped to improve sensitivity of the data based on diffuse polarised Enterprise and Innovation (DBEI), Ireland; NWO, The Netherlands; The Sci-
ence and Technology Facilities Council, UK; Ministry of Science and Higher
emission. The outcomes of this analysis are the following: Education, Poland.
1. We verified the reliability of the absolute ionospheric
Faraday rotation corrections estimated using the satellite-
based TEC measurements to be of ∼0.05 rad m−2 . We also References
demonstrated that diffuse polarised emission itself can be Asad, K. M. B., Koopmans, L. V. E., Jelić, V., et al. 2015, MNRAS, 451, 3709
used to account for the relative ionospheric corrections with Bracco, A., Jelić, V., Marchal, A., et al. 2020, A&A, 644, L3
respect to some reference observation; Bracco, A., Ntormousi, E., Jelić, V., et al. 2022, A&A, 663, A37
Brentjens, M. A. 2018, in Astrophysics and Space Science Library, 426, 159
2. We showed the feasibility of the developed stacking tech- Brentjens, M. A., & de Bruyn, A. G. 2005, A&A, 441, 1217
nique by combining 20 single-night observations into one Burn, B. J. 1966, MNRAS, 133, 67
150 h data set. The resulting Faraday cube has a noise Carretti, E., Vacca, V., O’Sullivan, S. P., et al. 2022, MNRAS, 512, 945

A119, page 11 of 15
 A&A 674, A119 (2023)

de Gasperin, F., Mevius, M., Rafferty, D. A., Intema, H. T., & Fallows, R. A. Offringa, A. R., de Bruyn, A. G., Zaroubi, S., et al. 2013, A&A, 549, A11
2018, A&A, 615, A179 Ordog, A., Booth, R., Van Eck, C., Brown, J.-A., & Landecker, T. 2019, Galaxies,
de Gasperin, F., Dijkema, T. J., Drabent, A., et al. 2019, A&A, 622, A5 7, 43
de Gasperin, F., Williams, W. L., Best, P., et al. 2021, A&A, 648, A104 O’Sullivan, S. P., Shimwell, T. W., Hardcastle, M. J., et al. 2023, MNRAS, 519,
Dickey, J. M., Landecker, T. L., Thomson, A. J. M., et al. 2019, ApJ, 871, 5723
106 Purcell, C. R., Van Eck, C. L., West, J., Sun, X. H., & Gaensler, B. M. 2020, RM-
Erceg, A., Jelić, V., Haverkorn, M., et al. 2022, A&A, 663, A7 Tools: Rotation measure (RM) synthesis and Stokes QU-fitting, Astrophysics
Ferrière, K., West, J. L., & Jaffe, T. R. 2021, MNRAS, 507, 4968 Source Code Library [record ascl:2005.003]
Hatanaka, T. 1956, PASJ, 8, 73 Sabater, J., Best, P. N., Tasse, C., et al. 2021, A&A, 648, A2
Heald, G., Braun, R., & Edmonds, R. 2009, A&A, 503, 409 Shimwell, T. W., Röttgering, H. J. A., Best, P. N., et al. 2017, A&A, 598, A104
Herrera Ruiz, N., O’Sullivan, S. P., Vacca, V., et al. 2021, A&A, 648, A12 Shimwell, T. W., Tasse, C., Hardcastle, M. J., et al. 2019, A&A, 622, A1
Hutschenreuter, S., Anderson, C. S., Betti, S., et al. 2022, A&A, 657, A43 Shimwell, T. W., Hardcastle, M. J., Tasse, C., et al. 2022, A&A, 659, A1
Jelić, V., Zaroubi, S., Labropoulos, P., et al. 2010, MNRAS, 409, 1647 Smolčić, V., Delvecchio, I., Zamorani, G., et al. 2017, A&A, 602, A2
Jelić, V., de Bruyn, A. G., Mevius, M., et al. 2014, A&A, 568, A101 Sotomayor-Beltran, C., Sobey, C., Hessels, J. W. T., et al. 2013, A&A, 552,
Jelić, V., de Bruyn, A. G., Pandey, V. N., et al. 2015, A&A, 583, A137 A58
Jelić, V., Prelogović, D., Haverkorn, M., Remeijn, J., & Klindžić, D. 2018, A&A, Spinelli, M., Bernardi, G., & Santos, M. G. 2019, MNRAS, 489, 4007
615, L3 Tasse, C., Shimwell, T., Hardcastle, M. J., et al. 2021, A&A, 648, A1
Kondapally, R., Best, P. N., Hardcastle, M. J., et al. 2021, A&A, 648, A3 Turić, L., Jelić, V., Jaspers, R., et al. 2021, A&A, 654, A5
Lenc, E., Gaensler, B. M., Sun, X. H., et al. 2016, ApJ, 830, 38 Van Eck, C. L., Haverkorn, M., Alves, M. I. R., et al. 2017, A&A, 597, A98
Mevius, M. 2018, RMextract: Ionospheric Faraday Rotation calculator, Astro- Van Eck, C. L., Haverkorn, M., Alves, M. I. R., et al. 2019, A&A, 623, A71
physics Source Code Library [record ascl:1806.024] van Haarlem, M. P., Wise, M. W., Gunst, A. W., et al. 2013, A&A, 556, A2
Murray, W. A. S., & Hargreaves, J. K. 1954, Nature, 173, 944 van Weeren, R. J., Williams, W. L., Hardcastle, M. J., et al. 2016, ApJS, 223, 2
Novak, M., Smolčić, V., Schinnerer, E., et al. 2018, A&A, 614, A47 Yizengaw, E., Moldwin, M., Galvan, D., et al. 2008, J. Atmos. Solar-Terrestrial
Offringa, A. R., van de Gronde, J. J., & Roerdink, J. B. T. M. 2012, A&A, 539, Phys., 70, 1541
A95 Zaroubi, S., Jelic, V., de Bruyn, A. G., et al. 2015, MNRAS, 454, L46

A119, page 12 of 15
 Šnidarić., I., et al.: A&A proofs, manuscript no. aa45124-22

Appendix A: Restoring the 011 observation using The polarised emission is now visible in the restored Faraday
Galactic polarised emission cubes. An example of polarised intensity is shown in Fig. A.1
(right image) at Faraday depth of −3.25 rad m−2 . Observed mor-
The 011 observation shows the largest shift with respect to phology of polarised emission in the restored 011 observation
the reference observation (≳ 1 rad m−2 ) among all observations and the reference observation are visually very similar. To
analysed in this work. We inspected the Faraday cube for this quantify this similarity, we calculate the Pearson correlation
observation in the polarised intensity. There is almost no emis- coefficient between the images of the highest peak value of the
sion visible in the Faraday cube in comparison to the reference Faraday depth spectrum in the polarised intensity of the two
observation. An example is given in Fig. A.1 (left image) at a observations and a ratio of their peak intensity distributions.
Faraday depth of −2.25 rad m−2 . The same figure (middle image) We are using in the calculation only the inner 3◦ × 3◦ of the
shows an image of the reference observation, but at a Faraday images. We get a correlation coefficient of 0.95, and find that
depth of −3.25 rad m−2 to account for a relative misalignment polarised emission in the restored cube is (73 ± 14)% of that in
between the two observations in Faraday depth. The lack of the the reference observation. The majority of the emission and its
observed polarised emission shows that RMion corrections were morphology is restored.
not applied properly to the data due to some unfortunate pro- To increase the percentage of the recovered brightness of
cessing error. We confirm this by inspection of the processing the observed emission even further, we would need to address
log files. the depolarisation that happens within ten-minute intervals. To
The calculated RMion correction for this observation is achieve that, we could re-image the eight-hour observation to
2.4 rad m−2 at the beginning of the observation, and then it even smaller time intervals, for example of one minute. How-
decreases to 1.5 rad m−2 within the first 430 minutes. This rel- ever, the signal-to-noise ratio in that case would be a limiting
ative change of 0.9 rad m−2 is enough to fully depolarize the factor for our methodology, making it out of the scope of current
signal at 150 MHz (see Sect. 1), if we do not correct the data for work.
it. In the remaining 50 minutes of the observation, it increases Furthermore, we assess if the estimated shifts can be applied
again to 1.6 rad m−2 . to the high-resolution images (6′′ ) of the same observation. We
To restore the polarised signal in 011 observation, we test if re-imaged a small part of the high-resolution data centred at a
the observed polarised emission itself can be used to account polarised source at RA 16h 24m 32s and Dec +56◦ 52′ 28′′ (Her-
for the ionospheric Faraday rotation correction that should be rera Ruiz et al. 2021, source with ID 01) in ten-minute intervals.
applied to the data. We first re-imaged the eight-hour 011 obser- Then, we ‘de-rotated’ the observed polarisation angle of each
vation by creating 48 Stokes QU images of ten-minute intervals ten-minute interval by the shift estimated using the very-low res-
of the observation. Then, for each ten-minute interval, we found olution data, combined high-resolution ten-minute intervals with
its relative shift in Faraday depth with respect to the full eight- the full eight-hour synthesis frequency cube, and applied the RM
hour reference observation following the methodology described synthesis to it. The resulting primary-beam-uncorrected Faraday
in Sect. 3. spectrum at a location of the polarised source is presented in
Figure A.2 shows the results (thick solid black line), which Fig. A.3 (black line). The source appears at a Faraday depth of
are compared with the corrections calculated using the satellite- ∼ 9.5 rad m−2 , as expected. Its recovered peak polarised flux is
based TEC measurements (black dashed line). The two curves 43% of the value reported by Herrera Ruiz et al. (2021), once
are showing the same trend. A systematic shift of ∼0.3 rad m−2 we take into account the primary beam correction at the location
between the two curves is due to the different nature of these two of the source (a factor of 2.3×). The same figure also shows the
methods. The satellite-based corrections give absolute RMion , Faraday spectrum before the ionospheric corrections are applied
while the corrections based on the observed polarised emis- (cyan line), where the source is fully depolarised. A success-
sion give relative values with respect to the used reference ful detection of the source demonstrates a potential of using the
observation (∆Φshift ). very low-resolution data to correct the high-resolution data. This
We also tested for any angular variations of ∆Φshift across method is computationally more efficient than the one that uses
the FoV. We did this by splitting the frequency cube spatially the high-resolution data only.
into quadrants and then repeating the procedure to find a rela-
tive time varying shift for each quadrant separately. The results
are over-plotted in Fig. A.2 with thin solid coloured lines. The
northwest (NW) and the southeast (SE) quadrants show relative
shifts which are consistent with the result of the full cube. Larger
values are found in the northeast (NE) quadrant; on average, the
shifts are larger by 0.096 rad m−2 than the one from the full cube.
In the southwest (SW) quadrant, they are smaller by 0.1 rad m−2 .
This points to a relative spatial gradient of ∼0.2 rad m−2 in the
northeast-southwest direction across the FoV.
We ‘de-rotated’ the observed polarisation angle of each ten-
minute interval by its estimated shift across the full image with
respect to the reference observation (∆Φshift ). This is done by
multiplying the complex polarisation given at each wavelength
2
(frequency) by exp−i2∆Φshift λ . Then, we combined all corrected
ten-minute intervals to obtain a Stokes QU cube over the full
eight-hour synthesis and used the RM synthesis to obtain the
final restored Faraday cubes of the 011 observation.

A119, page 13 of 15
 A&A 674, A119 (2023)

Observation 011, = 2.25 rad m 2 Reference observation, = 3.25 rad m 2 Corrected 011 observation, = 3.25 rad m 2
58° 1.0 58° 1.0 58° 1.0

0.8 0.8 0.8

]
56° 56° 56°
Declination (J2000)

Declination (J2000)

Declination (J2000)
1

1
mJy PSF 1 RMSF

mJy PSF 1 RMSF

mJy PSF 1 RMSF

0.6 0.6 0.6

54° 0.4 54° 0.4 54° 0.4

0.2 0.2 0.2

[
52° 52° 52°
0.0 0.0 0.0
248° 246° 244° 242° 240° 238° 248° 246° 244° 242° 240° 238° 248° 246° 244° 242° 240° 238°
Right Ascension (J2000) Right Ascension (J2000) Right Ascension (J2000)

Fig. A.1. Example of image in Faraday cube given in the polarised intensity at −2.25 rad m−2 for the 011 observation (left image), which is not
properly corrected for the RMion . There is almost no emission visible in comparison to the reference (014) observation (middle image), whose image
is given at −3.25 rad m−2 to account for a relative misalignment of +1.0 rad m−2 between the two observations. The polarised emission is visible in
the restored Faraday cube of the 011 observation (right image), which is corrected using the estimated ∆Φshift given in Fig. A.2.

2.4 RMion
full cube
or RMion [rad m 2]

2.2 NE quadrant
NW quadrant
2.0 SE quadrant
1.8 SW quadrant

1.6
1.4
shift

1.2
1.0
0 100 200 300 400 500
Time [min]
Fig. A.2. Estimated relative shifts in Faraday depth (∆Φshift ) of each
ten-minute interval of the 011 observation with respect to the full eight-
hour reference (014) observation (thick solid black line). The calculated
RMion corrections based on the satellite TEC measurements are plotted
with a thick dashed black line. The thin solid coloured lines give the
∆Φshift in the field-of-view quadrants.

Fig. A.3. Primary-beam-uncorrected Faraday spectrum at a location of

a polarised source (ID 01 in Herrera Ruiz et al. 2021) in the high-
resolution restored (solid black line) and unrestored (solid cyan line)
Faraday cube of the 011 observation. The source should be present at a
Faraday depth of +(9.44 ± 0.03) rad m−2 (red vertical line), as reported
by Herrera Ruiz et al. (2021).

A119, page 14 of 15
 Šnidarić., I., et al.: A&A proofs, manuscript no. aa45124-22

Appendix B: Supplementary figures

This appendix presents images in the polarised intensity of the final stacked Faraday cube described in Sect. 5.2.

Fig. B.1. Images of the ELAIS-N1 field in the polarised intensity given at Faraday depths of -15.0, -10.0, -5.0, -3.0, -1.0, +1.0, +3.0, +5.0,
+7.0, +10.0, +12.0 and +16.0 rad m−2 of the final stacked Faraday cube. The cube is based on ∼150 hours of the LOFAR observations in the
frequency range from 114.9 to 177.4 MHz. Angular resolution of the images is 4.3′ . These are primary-beam-uncorrected images with the noise of
27 µJy PSF−1 RMSF−1 .

A119, page 15 of 15