LOFAR 150-MHz observations of SS 433 and W 50

J.W. Broderick, R.P. Fender, J.C. A. Miller-Jones, S.A. Trushkin, A.J.

Stewart, G.E. Anderson, T.D. Staley, K.M. Blundell, M. Pietka, S. Markoff, et al.

► To cite this version:

J.W. Broderick, R.P. Fender, J.C. A. Miller-Jones, S.A. Trushkin, A.J. Stewart, et al. LOFAR 150-MHz observations of SS 433 and W 50. Monthly Notices of the Royal Astronomical Society, 2018, 475 (4), pp.5360-5377. 10.1093/mnras/sty081. hal-01730032

HAL Id: hal-01730032

https://hal.science/hal-01730032

Submitted on 13 Mar 2018

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

L'archive ouverte pluridisciplinaire **HAL**, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
<figure><img src="image_1.png" /></figure>

# LOFAR 150-MHz observations of SS 433 and W 50

J. W. Broderick,¹,²,³* R. P. Fender,¹,² J. C. A. Miller-Jones,⁴ S. A. Trushkin,⁵,⁶

A. J. Stewart,¹,² G. E. Anderson,¹,⁴ T. D. Staley,¹,² K. M. Blundell,¹ M. Pietka,¹

S. Markoff,⁷ A. Rowlinson,³,⁷ J. D. Swinbank,⁸ A. J. van der Horst,⁹ M. E. Bell,¹⁰,¹¹

R. P. Breton,¹² D. Carbone,⁷ S. Corbel,¹³,¹⁴ J. Eislöffel,¹⁵ H. Falcke,¹⁶,³

J.-M. Grießmeier,¹⁷,¹⁴ J. W. T. Hessels,³,⁷ V. I. Kondratiev,³,¹⁸ C. J. Law,¹⁹

G. J. Molenaar,⁷,²⁰ M. Serylak,²¹,²² B. W. Stappers,¹² J. van Leeuwen,³,⁷

R. A. M. J. Wijers,⁷ R. Wijnands,⁷ M. W. Wise³,⁷ and P. Zarka²³,¹⁴

*Affiliations are listed at the end of the paper*

Accepted 2017 December 21. Received 2017 December 21; in original form 2017 May 25

## ABSTRACT

We present Low-Frequency Array (LOFAR) high-band data over the frequency range 115–189 MHz for the X-ray binary SS 433, obtained in an observing campaign from 2013 February to 2014 May. Our results include a deep, wide-field map, allowing a detailed view of the surrounding supernova remnant W 50 at low radio frequencies, as well as a light curve for SS 433 determined from shorter monitoring runs. The complex morphology of W 50 is in excellent agreement with previously published higher frequency maps; we find additional evidence for a spectral turnover in the eastern wing, potentially due to foreground free-free absorption. Furthermore, SS 433 is tentatively variable at 150 MHz, with both a debiased modulation index of 11 per cent and a $\chi^2$ probability of a flat light curve of $8.2 \times 10^{-3}$. By comparing the LOFAR flux densities with contemporaneous observations carried out at 4800 MHz with the RATAN-600 telescope, we suggest that an observed ~0.5–1 Jy rise in the 150-MHz flux density may correspond to sustained flaring activity over a period of approximately 6 months at 4800 MHz. However, the increase is too large to be explained with a standard synchrotron bubble model. We also detect a wealth of structure along the nearby Galactic plane, including the most complete detection to date of the radio shell of the candidate supernova remnant G 38.7–1.4. This further demonstrates the potential of supernova remnant studies with the current generation of low-frequency radio telescopes.

**Key words:** stars: individual: SS 433 – ISM: individual objects: W50 – ISM: jets and outflows – ISM: supernova remnants – radio continuum: ISM – radio continuum: stars.

## 1 INTRODUCTION

The W 50 nebula (Westerhout 1958) is a large, non-thermal, Galactic radio source, at the centre of which is the X-ray binary system SS 433 (Stephenson & Sanduleak 1977), a known source of relativistic jets (Abell & Margon 1979; Fabian & Rees 1979; Liebert et al. 1979; Margon et al. 1979a,b; Milgrom 1979). These jets, with mean velocity 0.26c, precess about an axis with a period of 162.4 d, where the half-opening angle of the precession cone is 21° and the precession axis is inclined at 78° to the line of sight (Eikenberry et al. 2001; also see Margon 1984). This creates a corkscrew-shaped

trace on the sky as the jet components move ballistically outwards (Hjellming & Johnston 1981a,b).

At radio frequencies, W 50 takes the form of a circular shell of diameter 58 arcmin with lateral extensions (the ‘wings’ or ‘ears’) of 47 and 25 arcmin to the east and west, respectively (e.g. 327.5- and 1465-MHz maps in Dubner et al. 1998, henceforth D98). Originally classified as a supernova remnant (SNR) by Holden & Caswell (1969), it was suggested with the discovery of the jets in SS 433 that, given the remarkable alignment between the jet axis and the major axis of W 50, they are responsible for the elongated (~2° × 1°) shape of the nebula (Begelman et al. 1980). Hydrodynamical simulations by Goodall, Alouani-Bibi & Blundell (2011a) suggested that W 50 is about 20 000 yr old, and that several distinct episodes of jet activity, with different precession characteristics

<footnote>* E-mail: broderick@astron.nl</footnote>
to the current jets, can explain the morphology of W 50; this is
consistent with the kinematical investigation by Goodall, Blun-
dell & Bell Burnell (2011b). The asymmetry in the sizes of the
ears is usually ascribed to the known density gradient associated
with the Galactic plane; SS 433 is located at Galactic coordinates
(l, b) = (39:7, -2:2), and the major axis of W 50 is oriented at 20°
from the normal to the plane (e.g. Downes, Salter & Pauls 1981a;
Goodall et al. 2011a).

In observations at frequencies ≳1 GHz, flaring activity from
SS 433 has been found to occur on time-scales ≲100–200 d (e.g.
Seaquist et al. 1982; Fiedler et al. 1987; Vermeulen et al. 1993;
Trushkin, Bursov & Nizhelskij 2003; Pal et al. 2006; Trushkin,
Nizhelskij & Tsybulev 2012, 2014; Trushkin et al. 2016). During
flaring, flux densities of ~0.5–4 Jy are measured, on average a
factor of a few above the quiescent baseline value. Flaring has also
been detected at 408 and 843 MHz (e.g. Bonsignori-Facondi et
al. 1986; Vermeulen et al. 1993); at 408 MHz, the source is typically
about 1 Jy brighter than during quiescence (3.4 Jy versus 2.3 Jy;
Bonsignori-Facondi et al. 1986).

At low radio frequencies (<300 MHz), however, SS 433, and indeed W 50 as well, have been relatively poorly characterized. A brief summary of previous observations is as follows. First, using 80- and 160-MHz data from the Culgoora Circular Array (CCA), as well as near-simultaneous higher frequency measurements, Seaquist et al. (1980) suggested that the quiescent radio spectrum of SS 433 turns over at approximately 300 MHz. From two epochs, spread 3 weeks apart, the average 160-MHz flux density was 1.9 Jy, and the average 5σ upper limit at 80 MHz was 1.6 Jy. Secondly, Pandey et al. (2007) monitored SS 433 at 235 MHz with the Giant Metrewave Radio Telescope (GMRT) over the period 2002 July to 2005 January (observational cadence range 2–270 d). The average 235-MHz flux density was 2.71 Jy and the modulation index 21 per cent. Using simultaneous 610-MHz GMRT observations, the two-point spectral index<fnref n="1" /> $\alpha_{235}^{610}$ varied between −1.05 and −0.69, where the uncertainty per individual measurement is much less than this range. This is possibly due to optical depth variations at different stages of flaring (discussion in Pandey et al. 2007). Lastly, Miller-Jones et al. (2007) detected both SS 433 and W 50 with the Very Large Array (VLA) at 74 MHz; the root-mean-square (rms) noise level in their map was 192 mJy beam⁻<fnref n="1" />, with an angular resolution of about 100 arcsec. The flux density of SS 433 was 1.9 ± 0.2 Jy (Miller-Jones, private communication). Using these data, as well as previously published observations, including 83- and 102-MHz measurements from the Pushchino Radio Astronomy Observatory (Kovalenko, Pynzar' & Udal'Tsov 1994), Miller-Jones et al. (2007) demonstrated that the integrated radio spectrum of W 50 can be modelled with a single power law between 74 and 4750 MHz, where $\alpha = −0.51 ± 0.02$ (also see D98 and Gao et al. 2011).

With the advent of a new generation of novel, low-frequency ra-
dio telescopes, including the Low-Frequency Array (LOFAR; van
Haarlem et al. 2013), the Murchison Widefield Array (MWA; Tin-
gay et al. 2013) and the Long Wavelength Array (LWA; Ellingson
et al. 2013), there is now an excellent opportunity to embark on a
much more thorough exploration of the low-frequency synchrotron
radio emission from jet-producing sources such as SS 433. Topics
of particular interest include constraining models for intrinsic vari-
ability over wide frequency ranges (including frequency-dependent
time lags), investigating absorption processes that result in a spec-

tral turnover, as well as the presence of a low-frequency cut-off in
the synchrotron spectrum. In this paper, we present the results of a
LOFAR high-band observing campaign of SS433 and W50, carried
out over the frequency range 115–189 MHz during the period 2013
February to 2014 May. In Section 2, we describe the observations,
and how the data were calibrated and imaged. A deep continuum
map is presented in Section 3, while in Section 4 we combine
our data with previous observations at lower and higher frequen-
cies to investigate the radio spectra of the components of W50.
Our low-frequency SS 433 monitoring observations are analysed in
Section 5, including a comparison with contemporaneous data at
4800 MHz from the RATAN-600 telescope. Then, in Section 6, we
discuss the detection of a nearby candidate SNR, G 38.7–1.4, in
our deep image. Finally, we conclude in Section 7.

## 2 LOFAR OBSERVATIONS AND DATA REDUCTION

A detailed description of LOFAR can be found in van Haarlem
et al. (2013). In this section, we provide a summary of our obser-
vations, carried out with the high-band antennas (HBA; full operat-
ing frequency range 110–240 MHz). An observing log is given in
Table 1.

2.1 Deep continuum observation and monitoring runs

Our first HBA observation of SS 433 and W 50 was a deep 4-h run carried out with the 'HBA Dual' mode, centred on transit (elevation 42°). The frequency range was 115–163 MHz, covered by a total of 244 sub-bands, each with a bandwidth of 195.3 kHz. We used the Dutch stations only: 24 stations in the core, each with two closely spaced sub-stations, and 13 remote stations. The projected baselines range from about 29 m to 79 km. The primary beam full width at half-maximum (FWHM) at 150 MHz is 3:80 for the core sub-stations and 2:85 for the remote stations.

The observation was carried out in blocks of 15 min: 2 min on the primary flux density calibrator, 3C 380, followed by 11 min on the target. Although LOFAR is a software telescope, rapid electronic beam switching was not possible when our study took place: the remaining 2 × 1 min in the block was needed to switch between the calibrator and target (and vice versa). By the end of the run, 16 snapshots of SS 433 and W 50 had been obtained, corresponding to a total integration time of 176 min.

The primary goal of the shorter runs was to monitor the flux
density of SS 433. The setup was very similar, except that the
frequency coverage was extended to 115–189 MHz, spread over
380 sub-bands. In general, 2 × 11- or 2 × 13-min snapshots,
near transit, were obtained per run, along with 2 × 2-min scans
of 3C 380. However, the data from 2013 November 16 were ob-
tained with a slightly different setup: one 19-min snapshot, approx-
imately centred on transit, preceded by a 1.5-min observation of
3C 380.

## 2.2 Pre-processing, calibration, and imaging

For all observations, data were recorded with a time-step of 2 s
and 64 channels per sub-band. Pre-processing was carried out us-
ing standard methods (e.g. van Haarlem et al. 2013). First, radio-
frequency interference was removed using AOFLAGGER (Offringa
et al. 2010; Offringa, van de Gronde & Roerdink 2012), and the
two edge channels at both ends of each sub-band were completely
flagged, reducing the bandwidth per sub-band to an effective value

<footnote n="1">¹ In this paper, we use the convention $S_v \propto v^\alpha$, where $S_v$ is the flux density at frequency $v$, and $\alpha$ is the spectral index.</footnote>
<table><caption>**Table 1.** LOFAR high-band observations of SS 433 and W 50. MJD is the Modified Julian Date at the halfway point of each observation. The prefix 'M' in column 1 refers to the monitoring observations. See Section 5 for a description of how the SS 433 flux densities were obtained; uncertainties are ±1σ. The on-source integration time is also given for each run, as well as the LOFAR observation IDs.</caption><thead><tr><th>Run</th><th>Date</th><th>MJD</th><th>S<sub>SS433</sub><br>(Jy)</th><th>t<sub>int</sub><br>(min)</th><th>No. stations<br>(core, remote, int.)</th><th>IDs<br>(range)</th></tr></thead><tbody><tr><td>Deep</td><td>2013 February 18</td><td>56341.38</td><td>2.08 ± 0.21</td><td>176</td><td>24, 13, 0</td><td>L95041–L95072</td></tr><tr><td>M1</td><td>2013 March 20</td><td>56371.28</td><td>2.09 ± 0.21</td><td>22</td><td>24, 13, 0</td><td>L106054–L106057</td></tr><tr><td>M2</td><td>2013 April 25</td><td>56407.22</td><td>1.65 ± 0.17</td><td>22</td><td>16, 13, 8<sup>a</sup></td><td>L128512–L128515</td></tr><tr><td>M3</td><td>2013 May 23</td><td>56435.12</td><td>1.69 ± 0.17</td><td>22</td><td>22, 13, 0</td><td>L136070–L136073</td></tr><tr><td>M4</td><td>2013 July 25/26</td><td>56499.00</td><td>1.49 ± 0.15</td><td>22</td><td>19, 13, 0</td><td>L166060–L166063</td></tr><tr><td>M5</td><td>2013 August 22</td><td>56526.87</td><td>1.60 ± 0.16</td><td>22</td><td>24, 14, 0</td><td>L169277–L169280</td></tr><tr><td>M6</td><td>2013 September 20</td><td>56555.86</td><td>1.74 ± 0.18</td><td>22</td><td>24, 14, 0</td><td>L175340–L175343</td></tr><tr><td>M7</td><td>2013 September 28</td><td>56563.83</td><td>1.64 ± 0.17</td><td>22</td><td>24, 14, 0</td><td>L178910–L178913</td></tr><tr><td>M8</td><td>2013 November 16</td><td>56612.63</td><td>2.14 ± 0.21</td><td>19<sup>b</sup></td><td>24, 14, 0</td><td>L188650–L188651</td></tr><tr><td>M9</td><td>2013 December 11</td><td>56637.55</td><td>2.13 ± 0.21</td><td>26</td><td>24, 14, 0</td><td>L193989–L193992</td></tr><tr><td>M10</td><td>2014 February 11</td><td>56699.39</td><td>2.62 ± 0.26</td><td>26</td><td>24, 11, 0</td><td>L204224–L204227</td></tr><tr><td>M11</td><td>2014 March 12</td><td>56728.30</td><td>2.06 ± 0.21</td><td>26</td><td>24, 14, 0</td><td>L207137–L207140</td></tr><tr><td>M12</td><td>2014 April 24</td><td>56771.25</td><td>1.89 ± 0.19</td><td>26</td><td>24, 14, 0</td><td>L215295–L215298</td></tr><tr><td>M13</td><td>2014 May 10</td><td>56787.10</td><td>1.89 ± 0.19</td><td>26</td><td>24, 14, 0</td><td>L228277–L228280</td></tr></tbody></table>

Notes. Additional monitoring observations on 2013 September 16 and 2014 January 8 were unusable due to various observational and data quality issues.

<sup>a</sup>The observation on 2013 April 25 included eight international stations in place of an equivalent number of core stations, potentially allowing very-high-resolution images to be generated, but we only consider the Dutch stations in this paper.

<sup>b</sup>This snapshot was obtained as part of a second deep observation, but of lower data quality. The 'HBA Dual Inner' configuration was used, in which the primary beam FWHM of the remote stations is matched to the core sub-stations by using the inner 24 tiles of the remote stations only.

of 183.1 kHz. Secondly, because the very bright 'A-team' source Cygnus A is only 17.8 from 3C 380, the 'demixing' algorithm (van der Tol, Jeffs & van der Veen 2007) was used to subtract the response of the former from the visibilities of the latter. Lastly, for practical reasons concerning data volume, and the computing time required for calibration and imaging, the data were averaged in time and frequency: to 10 s per time-step and either four frequency channels per sub-band (deep run) or one frequency channel per sub-band (monitoring runs). For the baseline ranges considered in this paper, the effects of bandwidth and time-average smearing (e.g. Bridle & Schwab 1999) in the averaged data are well within the calibration uncertainty (Section 2.3), even at the edges of the field.

Calibration and imaging were also carried out using standard practices. Each 3C 380 sub-band was calibrated using a model of the source defined in Scaife & Heald (2012); the absolute flux density scale is that of Roger, Bridle & Costain (1973). The gain amplitudes were transferred to the target sub-bands in the same observing block. We then performed phase-only calibration on the target field using data from the global sky model developed by Scheers (2011), within a radius of 10° from the position of SS 433. The basis for our model of the field was the 74-MHz VLA Low-Frequency Sky Survey (VLSS(r); Cohen et al. 2007; Lane et al. 2014), with spectral index information being obtained by cross-correlating the relevant VLSS catalogue entries with the 1.4-GHz NRAO VLA Sky Survey (NVSS; Condon et al. 1998). To increase the signal-to-noise and hence the accuracy of the phase-only calibration step, solutions were computed for at most 20 sub-bands at a time (Δν ≈ 4 MHz). Furthermore, we restricted the projected baseline range over which the phase solutions were determined to 250–5000 λ (i.e. 500 m to 10 km at 150 MHz), so as to avoid the shortest and longest baselines where our sky model will generally be an inadequate representation of the true sky brightness distribution.

Prior to imaging each run, the sub-bands were combined into six, approximately evenly split bands, and the data were visually inspected; any further errant data points were manually flagged.

A primary-beam-corrected image was then generated for each band with AWIMAGER (Tasse et al. 2013); we used a robust weighting parameter (Briggs 1995) of 0. At the time of analysis, the AWIMAGER did not fully account for variations of flux density with frequency in the deconvolution process (e.g. Sault & Wieringa 1994); by making multiple images across the full bandwidth, we kept the fractional bandwidth of the data being imaged to <10 per cent in all cases, which helps to minimize the possible reduction in image quality due to effects associated with this issue.

For the deep run, a maximum projected baseline of 4000 λ was used when imaging; for the monitoring runs, the baseline range was 100–3000 λ. In the latter case, the lower cut-off reduced the contribution from the diffuse emission from W 50, while both higher cut-offs gave the most reliable images for the given uv coverage and relatively simple calibration procedure. A multiscale CLEAN option had not yet been implemented in the AWIMAGER; it will be possible in the future to more accurately image the significant amount of extended emission in this field.

To make a final map for each run, the six separate images were first convolved to a common angular resolution using the task CONVOL in MIRIAD (Sault, Teuben & Wright 1995). We then stacked and averaged these resulting maps together with the MIRIAD tasks IMCOMB and AVMATHS. Equal weights were used; the effective frequency is about 140 MHz for the deep run and 150 MHz for the monitoring runs.

Lastly, after a first pass at imaging, we included an additional round of phase-only self-calibration with a model generated from the LOFAR image. The source finder PYBDSF (Mohan & Rafferty 2015) was used to generate the model; it was restricted to sources sufficiently compact such that they could be well described with either a single Gaussian or a moderate number (≤6) of Gaussians. While we subsequently used these self-calibration solutions (determined over the same projected baseline range as before), we found that they did not make a significant difference to the image properties.
## 2.3 In-band flux density accuracy

At the time of analysis, an issue existed with the standard LOFAR HBA beam model such that the in-band spectral index could be artificially steep. The magnitude of the effect is dependent on the elevation of the target in a given observation, and the distance between the primary calibrator and the target. Unfortunately, this currently precludes an investigation of how the radio properties of SS 433 and W50 potentially vary within the LOFAR high band. However, to ensure that the average HBA flux density scale is as accurate as possible, approximate corrections can be made based on a preliminary analysis of full simulations of the electromagnetic properties of the dipoles in the array, including the effects of mutual coupling and grating lobe attenuation (Heald, private communication).

For the deep observation, such corrections were applied to each of the six images before they were convolved to a common resolution and averaged together; we divided by factors ranging from 1.6–0.8 between 115–163 MHz. To first order, this process removes an artificially steep in-band spectral index of approximately –2.5 in this case. Unfortunately, the simulations are currently incomplete above 163 MHz, meaning that we are not able to fully assess the effect on the monitoring runs. However, based on an extrapolation of the available data between 115 and 163 MHz, we estimate that the average correction across a bandwidth of 115–189 MHz should, to first order, fall within the general calibration uncertainty: for all runs, based on the flux stability of sources in the target field, along with 3C 380 itself, we estimate that the internal flux density calibration uncertainty is ~10 per cent. Hence, we have made no flux density scale adjustments to the images from the monitoring runs.

## 2.4 Absolute flux density calibration

We also cross-checked the band-averaged LOFAR flux densities of the SNRs 3C 396 and 3C 397 (also known as G 039.2–00.3 and G 041.1–00.3, respectively; e.g. Green 2014), as well as 4C +06.66, with flux densities measured at 160 MHz with the CCA (Slee 1995). Using the information presented in Slee (1995) and Scaife & Heald (2012), both our LOFAR data and the data from Slee (1995) are tied to absolute flux density scales that agree to within ~3 per cent. We measured the LOFAR flux densities from our deep image, which has the more complete uv and baseline coverage, and, additionally, we interpolated the 160-MHz CCA flux densities to 140 MHz using the 80-MHz CCA flux densities and two-point spectral indices $\alpha_{80}^{160}$ that were also reported in Slee (1995).

We find that our LOFAR flux densities for 3C 396 and 4C +06.66 agree with the interpolated 140-MHz flux densities from Slee (1995) to within approximately ±2 per cent, giving us confidence in the absolute flux density calibration. However, the LOFAR measurement for 3C 397 is about 25 per cent lower. One possibility to explain at least some of this discrepancy is that the Slee (1995) measurement at 160 MHz has recorded more diffuse emission from the source, and/or there is significant confusion in this case from the Galactic plane, owing to the lower angular resolution, which is roughly twice as coarse as in our study.

To investigate further whether this finding for 3C 397 is a potential outlier, we conducted a similar comparison for 3C 398 (G 043.3–00.2), which we also detect, albeit well into a sidelobe of the primary beam; this source is 4:1 to the north of SS 433. Despite the fact that a flux density measured so far from the phase centre may inherently be unreliable, interestingly the agreement with the adjusted measurement from Slee (1995) is well within the ±2 per cent margin described above for 3C 396 and 4C +06.66.

Given the above results for 3C 396 and 4C +06.66 (and 3C 398), we conclude that our data have been correctly tied to the Roger et al. (1973) absolute flux density scale, and that the main calibration uncertainty is an internal one, as outlined above in Section 2.3. Moreover, we will show throughout this paper that there is a good agreement between our LOFAR flux densities and previously published flux density measurements. Of our key conclusions (Section 7), further evidence for spectral curvature in the eastern wing of W50 (Section 4), as well as the corresponding fit parameters, may be affected if the calibration uncertainty is larger than we have calculated.

Generally speaking, it is important to note that absolute flux density calibration is well known to be challenging at low frequencies with aperture arrays such as LOFAR, particularly near and along the Galactic plane. Due to the excellent sensitivity to diffuse extended emission, future planned releases of Galactic plane catalogues from the LOFAR Multifrequency Snapshot Sky Survey (MSSS; Heald et al. 2015), and the GaLactic and Extragalactic All-sky Murchison Widefield Array survey (GLEAM; Hurley-Walker et al. 2017), will be highly valuable comparison data sets.

## 3 A LOFAR MAP OF SS 433 AND W50

Our deep continuum map is presented in Fig. 1, with a zoomed-in view of SS 433 and W50 shown in Fig. 2. This is the most detailed image of W50 made thus far at frequencies below 300 MHz. The angular resolution is 78 arcsec × 55 arcsec (beam position angle BPA = 20°; throughout this paper we use the convention that the BPA is measured north through east), comparable to that of the 327.5- and 1465-MHz VLA images from D98; in particular, the resolution of the latter map from D98 is 56 arcsec × 54 arcsec (BPA −61°:4). Indeed, at a similar angular resolution, there is a clear consistency between the LOFAR and D98 maps: we also observe the three main components of the nebula (i.e. the central shell, as well as the eastern and western wings) and the complex structure contained within. A similar correspondence can be seen when comparing the LOFAR data with the recent study by Farnes et al. (2017), who conducted a full-Stokes broad-band analysis of SS 433 and W50, centred at about 2.3 GHz, with the Australia Telescope Compact Array (ATCA); the angular resolution of their map is 2.5 arcmin × 1.9 arcmin (BPA 75°:9). Within the wide field of view, we have also detected a wealth of structure near and along the Galactic plane; we return to this topic in Section 6.

In the vicinity of W50, the rms noise level is ~4–6 mJy beam⁻¹. Source and sidelobe confusion so near the Galactic plane, together with the presence of undersampled diffuse emission, are responsible for the noise level being much higher than the theoretical value of ~0.08 mJy beam⁻¹. Moreover, there are clear variations off-source (henceforth referred to as the ‘background’), including ‘negative bowls’, due to the undersampling of extended emission.

### 3.1 Flux densities

A Gaussian fit (using `IMFIT` in `MIRAD`) to SS 433 yields a peak flux density of 2.01 ± 0.20 Jy beam⁻¹ and an integrated flux density of 2.25 ± 0.23 Jy. Note that the 10 per cent calibration uncertainty is the dominant error term. There is very marginal evidence of source extension (at about the 1.6σ level along both the major and minor axes of the fitted Gaussian), but the deconvolved position angle, approximately 30°, is not consistent with the jet axis and precession cone. Intriguingly, however, this position angle is consistent with the general direction of the low-surface-brightness ‘ruff’ emission
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Deep LOFAR 140-MHz high-band continuum map, with dimensions 6.0° × 7.9°, centred on SS 433 and W 50. The angular resolution is 78 arcsec × 55 arcsec (BPA 20°), and the pixels are 10 arcsec × 10 arcsec in size. We have imaged beyond the half-power points of the core and remote station primary beams, which are 4.1 and 3.0 FWHM, respectively, for an average frequency of 140 MHz. The colour bar range does not include all pixel values (maximum 10.4 Jy beam⁻¹); the contrast level has been chosen to show the significant amount of extended emission in this field, and many sources, including SS 433, are saturated using this scheme. SS 433 and W 50 are labelled, as well as several other sources that are discussed in the text.</figcaption></figure>

observed at GHz frequencies (e.g. Blundell et al. 2001), albeit on much smaller angular scales (tens of mas) than in our map.

Assuming that W 50 has a central shell with radius 29 arcmin (D98), we calculated the integrated flux densities of the various components of the nebula by summing the pixels in the respective regions (indicated in Fig. 2) and making corrections for the non-zero background levels. Each background level correction value was determined by measuring the mean flux density in a series of boxes around the component and then taking the average of this set of means; moreover, the standard deviation of the set was used as the uncertainty for the correction. The mean offsets that were

subtracted (per pixel) range from −15.5 to 6 mJy beam⁻¹. Twisted-plane fits were also investigated, but the results were deemed to be less reliable. Henceforth in this section, as well as Sections 3.2 and 3.3, all quoted integrated flux densities and surface-brightness values are background corrected, unless stated otherwise.

The integrated flux densities for the eastern wing, central shell and western wing are 48 ± 13, 121 ± 25, and 38 ± 8 Jy, respectively. The flux density errors were calculated by combining the 10 per cent calibration uncertainty and the uncertainty in the background level in quadrature; there is also a much smaller statistical uncertainty from the noise level. The 2.25-Jy contribution from SS 433 was
<figure><img src="image_3.png" /></figure>

**Figure 2.** The map from Fig. 1, zoomed in on the region containing SS 433 and W 50. We define the central shell as in D98: a circle, radius 29 arcmin, centred on SS 433. Note, however, that the true centre of W 50 has been determined to be approximately 5 arcmin to the east–south–east of SS 433, along the jet axis (Watson et al. 1983; Lockman, Blundell & Goss 2007; Goodall et al. 2011a). The eastern and western wings are labelled; the dotted lines approximately enclose the regions from which we measured the flux densities (discussion in Section 3.1). We also indicate the positions of a selection of filaments, the spur (or ‘chimney’), two artefacts, and the H II region S 74 (discussion in Sections 3.2 and 3.3). The white cross refers to the 327.5-MHz position of a variable point source (Section 5.4). As in Fig. 1, the contrast level is chosen to show the significant amount of extended emission from W 50.

subtracted from the central shell, but we do not make any corrections for possibly unassociated compact sources in the three regions; their effects are essentially negligible in comparison to the other flux density uncertainties.

Summing the integrated flux densities for the three components stated in the previous paragraph, and taking into account the relative contributions of systematic and statistical uncertainties when propagating the errors, our total integrated flux density for W 50 at 140 MHz is $207 \pm 33$ Jy. This value is consistent with the previously calculated integrated radio spectrum (see discussion in Section 1), giving us further confidence in our calibration and background subtraction procedures. However, depending on the relative contributions of the negative bowls and diffuse Galactic emission in Figs 1 and 2, one possibility is that the mean background corrections should be more negative in value (i.e. the flux densities of the three components of W 50 should be higher), and that the respective uncertainties are overestimated. As will become apparent in Section 4.1, this would move the flux densities of the eastern wing and central shell upwards towards the extrapolated single power-law radio spectra determined by D98 at higher frequencies. In the former case, this would reduce the evidence for spectral curvature at low frequencies. However, it seems very unlikely that an incorrect background correction alone could explain the deviation from a single power-law radio spectrum for the eastern wing, as there is insufficient evidence in Figs 1 and 2 for the significant negative bowl (mean $\approx -46$ mJy beam$^{-1}$) that would be needed in the vicinity of

this component. Furthermore, if the background correction for the western wing were to be a larger negative value, then the corrected 140-MHz flux density would quickly become overestimated.

Various second-order effects may result from the use of the Clark (1980) CLEAN algorithm in the AWIMAGER on a target such as W 50, including residual deconvolution sidelobes and artefacts, as well as possible CLEAN bias (e.g. Condon et al. 1998). We plan to re-image our data with multiscale CLEAN and maximum entropy algorithms once they are implemented in the AWIMAGER; an alternative would be to use sparse image reconstruction methods (e.g. Dabbech et al. 2015; Garsden et al. 2015; Girard et al. 2015; Pratley et al. 2018). Other techniques, such as the one outlined in Braun & Walterbos (1985), could also be considered as part of future work.

## 3.2 Low-frequency morphology of W 50

The arcs and filaments observed at higher frequencies are also prominent in our map. These include the thin filament at a position angle of 80° (labelled in Fig. 2), discussed in D98 as being positionally associated with the precession cone of SS 433, with a 140-MHz surface brightness ranging from about 50 to 90 mJy beam$^{-1}$. In addition, we see the set of north–south filaments which are suggestive of large-scale traces of the precessing jets; the brighter of these are also labelled in Fig. 2, with hints of similar fainter features in between. The peak flux densities of the eastern and western wings are approximately 150 and 170 mJy beam$^{-1}$, respectively; they occur
<figure><img src="image_4.png" /><figcaption>Figure 3. The spur, or ‘chimney’, at 140 MHz. To outline its shape, we have overlaid a radio contour at the 30 mJy beam⁻¹ level (not corrected for the non-zero background).</figcaption></figure>

in the brightest north-south filament and near the edge-brightened western boundary, respectively.

We detect the ‘chimney’ (Fig. 3), previously observed at 1465 MHz by D98 and also seen by Farnes et al. (2017). This is a spur-like feature to the north of the bright filament in the eastern wing, near the position RA 19h 16m 18s, Dec. +04° 59′ 30″ (J2000). D98 suggested that it may be analogous in morphology to a spur, approximately an order of magnitude smaller, at the northern edge of the Crab nebula (e.g. Fesen & Gull 1986; Rudie, Fesen & Yamada 2008; Black & Fesen 2015). The approximate dimensions at 140 MHz are very similar to those at 1465 MHz: ~6 arcmin in length and ~7 arcmin in width at the base. This corresponds to a projected linear size of about 10 pc × 11 pc.<fnref n="2" /> The average surface brightness at 140 MHz is ~15–20 mJy beam⁻¹. Using the 1465-MHz image from D98, convolved to the resolution of our data (see Section 4.2 for further details, where we construct a two-point spectral index α₁₄₆₅ map of W 50), we estimate a two-point spectral index α₁₄₀ ~ −0.5, with an error of about 30–40 per cent due to the faintness of the emission and the uncertainties in estimating the background level at both frequencies.

Farnes et al. (2017) identified a potential second chimney, 20 arcmin to the west of the chimney described above. However, it is not visible in Fig. 3. Our non-detection would be consistent with the suggestion of Farnes et al. (2017) that this feature is likely to be a low-level artefact from the deconvolution process.

Related to the previous paragraph, there are, however, some generally low-level artefacts in our LOFAR map, likely due to residual calibration and/or deconvolution errors. For example, 3.7 arcmin to the north-north-west of SS 433, at a position angle of −20°, there is a spurious feature consisting of two components, with a peak flux density of approximately 80 mJy beam⁻¹. A similar structure can be seen below the central shell of W 50 near the compact radio source NVSS J191252+043104.

### 3.3 The H II region S 74

S 74 (Crampton, Georgelin & Georgelin 1978), a foreground H II region (distance ~2–3 kpc; e.g. Forbes 1989; Brand & Blitz 1993; Lockman et al. 2007; Alves et al. 2010) to the north-west of W 50, is also detected at 140 MHz. The peak flux density is ~50 ± 20 mJy beam⁻¹, with the quoted 1σ error dominated by the uncertainty in the background level. In the convolved 1465-MHz image from D98, as mentioned above in Section 3.2, we estimate a peak flux density of 75 ± 10 mJy beam⁻¹; thus, a representative two-point spectral index α₁₄₀³⁵ at the position of peak intensity is 0.2 ± 0.2. Although H II regions can become optically thick at low frequencies (e.g. Kassim et al. 1989; Copetti & Schmidt 1991; Subrahmanyan & Goss 1995; Omar, Chengalur & Anish Roshi 2002), our estimated spectral index for S 74 is consistent with the canonical spectral index of ~0.1 in the optically thin regime: the difference is only 1.5σ. The radio spectrum may in fact be less inverted than our estimate depending on the relative contribution to the brightness temperature at 140 MHz from the Galactic emission behind S 74 (e.g. formulae in Nord et al. 2006; also see sky maps in e.g. Haslam et al. 1982; Guzmán et al. 2011; Remazeilles et al. 2015; Zheng et al. 2017).

In a 30.9-MHz Galactic plane survey, conducted by Kassim (1988) with the Clark Lake Teepee-Tee Telescope, S 74 appears as a negative absorption hole against the Galactic background (we refer the reader to Nord et al. 2006 and Su et al. 2017a,b for other similar cases). There is evidence for such behaviour at 74 MHz as well (Miller-Jones, private communication). Although further measurements at low frequencies are needed to establish exactly where the spectral turnover occurs, we can make some predictions based on an estimated turnover frequency. First, the linear emission measure EM can be calculated approximately using the expression

$$ \tau_v \approx 3.28 \times 10^{-7} \left( \frac{T_e}{10^4 \text{ K}} \right)^{-1.35} \left( \frac{\nu}{\text{GHz}} \right)^{-2.1} \left( \frac{\text{EM}}{\text{pc cm}^{-6}} \right), \quad (1) $$

where τ_v is the optical depth at frequency ν, and T_e is the electron temperature of the ionized gas in the H II region (e.g. Mezger & Henderson 1967). Taking a turnover frequency of ~100 MHz, where the optical depth is unity, and assuming T_e ~ 6000–7000 K for S 74 (e.g. see derived relationships between electron temperature and Galactocentric distance for H II regions in Shaver et al. 1983; Paladini, Davies & De Zotti 2004; Alves et al. 2012), then EM ~ 1.2–1.5 × 10⁴ pc cm⁻⁶.

Secondly, due to the diffuse Galactic emission in the vicinity of the source, it is somewhat difficult to assess the full extent of S 74 in our LOFAR map, but a rough estimate is ~7 arcmin × 13 arcmin (~5.1 pc × 9.5 pc at a distance of 2.5 kpc; also see e.g. Altenhoff et al. 1970; Geldzahler, Pauls & Salter 1980; Kuchar & Clark 1997 for other size estimates). Then, using the simple approximation that EM ≈ n_e<fnref n="2" /> L, where n_e is the electron number density in S 74, and L is the path length through the source along the line of sight, n_e ~ 40–45 cm⁻³ assuming that the path length is roughly the average of the projected dimensions on the sky (i.e. L ~ 7.3 pc). Referring to the sub-classes in e.g. Kurtz (2005), our very crude analysis suggests that S 74 is most consistent with ‘classical’ H II regions.

## 4 RADIO SPECTRA

As mentioned in Section 1, Miller-Jones et al. (2007) found that the integrated radio spectrum of W 50 is well described by a single power law between 74 and 4750 MHz (α = −0.51 ± 0.02). Using their data at 74 MHz, as well as the previously published data from

<footnote n="2">² In this paper, we assume that the distance to SS 433 and W 50 is 5.5 kpc (Blundell & Bowler 2004; Lockman et al. 2007). At this distance, an angular size of 1 arcmin corresponds to a projected linear size of 1.6 pc.</footnote>
<table><caption>**Table 2.** Flux densities and two-point spectral indices for the different components of W 50. The 74-MHz flux densities (Miller-Jones, private communication) are on the same flux density scale as the LOFAR data (see Section 4.1); this scale is consistent with the one used in D98, to an accuracy within the typical uncertainties quoted in the table. The 327.5-MHz flux densities were inferred using the 1465-MHz flux densities given in D98, as well as the $\alpha_{327.5}^{1465}$ values, which were determined by D98 using a linear regression method. Flux density uncertainties at 327.5 and 1465 MHz have been estimated based on the uncertainties reported by D98 for the entire W 50 nebula. At 327.5 MHz, our total inferred flux density (177 ± 21 Jy) is consistent with the value reported for the entire nebula in D98 (160 ± 20 Jy). Apart from $\alpha_{327.5}^{1465}$, two-point spectral indices and the associated uncertainties were deduced from the integrated flux densities.</caption><thead><tr><th>Component</th><th>S<sub>74</sub><br>(Jy)</th><th>S<sub>140</sub><br>(Jy)</th><th>S<sub>327.5</sub><br>(Jy)</th><th>S<sub>1465</sub><br>(Jy)</th><th>α<sub>74</sub><sup>140</sup></th><th>α<sub>140</sub><sup>327.5</sup></th><th>α<sub>327.5</sub><sup>1465</sup></th><th>α<sub>140</sub><sup>1465</sup></th></tr></thead><tbody><tr><td>Eastern wing</td><td>61 ± 5</td><td>48 ± 13</td><td>51 ± 6</td><td>15.0 ± 1.5</td><td>-0.38 ± 0.44</td><td>0.07 ± 0.35</td><td>-0.82 ± 0.05</td><td>-0.50 ± 0.12</td></tr><tr><td>Central shell</td><td>245 ± 16</td><td>121 ± 25</td><td>99 ± 12</td><td>43 ± 4</td><td>-1.11 ± 0.34</td><td>-0.24 ± 0.28</td><td>-0.56 ± 0.05</td><td>-0.44 ± 0.10</td></tr><tr><td>Western wing</td><td>54 ± 4</td><td>38 ± 8</td><td>27 ± 3</td><td>13.0 ± 1.3</td><td>-0.55 ± 0.35</td><td>-0.40 ± 0.28</td><td>-0.49 ± 0.09</td><td>-0.46 ± 0.10</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Broadband radio spectra for the various components of W 50, using the flux densities from Table 2. The solid lines are either a fit based on a free-free absorption model (left-hand panel; see Section 4.1 for further details) or linear fits in log (S<sub>v</sub>)–log (ν) space (middle and right-hand panels). Inverse-variance weighting was used. Dashed lines represent radio spectra based on the two-point spectral indices $\alpha_{327.5}^{1465}$ from D98. For the eastern wing, we report the approximate turnover frequency, where the spectral index is zero, while we give the fitted spectral indices for the central shell and western wing. The full set of parameters can be found in Table 3. Error bars are $\pm 1\sigma$.</figcaption></figure>

D98 at 327.5 and 1465 MHz, they also demonstrated that the only component of W 50 that shows evidence for significant spectral turnover at low frequencies is the eastern wing.

Our LOFAR 140-MHz map not only bridges the frequency gap in the Miller-Jones et al. (2007) analysis, but has sufficient dynamic range and angular resolution to permit us to study changes in the spectral index across the W 50 nebula, either using the integrated flux densities for the specific components, or pixel by pixel. We now discuss our findings using a combination of both approaches.

## 4.1 Integrated flux densities

Table 2 lists the integrated flux densities for the eastern wing, central shell, and western wing at 74, 140, 327.5, and 1465 MHz. While our LOFAR data are on the Roger et al. (1973) flux density scale, the 74-MHz observations used in the Miller-Jones et al. (2007) study were calibrated with the Baars et al. (1977) scale. Hence, we took the 74-MHz flux densities for the three components (Miller-Jones, private communication), scaling each measurement up by 10 per cent and propagating the corresponding uncertainty term of 6 per cent with the previous uncertainty; these values were found to be appropriate in an analogous rescaling of the 74-MHz VLSSr (Lane et al. 2014).

In Table 2, we also give two-point spectral indices calculated from the integrated flux densities, although we note that the $\alpha_{327.5}^{1465}$ values are from D98 and were calculated with a linear regression method using the pixel values in specified regions of the source (further details in Section 4.2). In addition, broad-band radio spectra are presented in Fig. 4.

As can be seen in Fig. 4, our LOFAR flux densities are consistent with the previously published measurements, in particular strengthening the evidence for a spectral turnover in the eastern wing. We have fitted each spectrum with either a free-free absorption model (eastern wing), or a linear fit in log (S<sub>v</sub>)–log (ν) space (central shell and western wing). The best-fitting parameters are reported in Table 3, and the fits are overlaid in the panels in Fig. 4. The reduced $\chi^2$ value for the western wing (0.035) suggests that our uncertainties are possibly overestimated in this case.

In the case of the eastern wing, our model is of the form

$$S_v = S_{140,u} \left( \frac{\nu}{140 \text{ MHz}} \right)^\alpha \exp \left[ -(\tau_{140}) \left( \frac{\nu}{140 \text{ MHz}} \right)^{-2.1} \right], \quad (2)$$

where $\tau_{140}$ is as before, $S_{140,u}$ is the unabsorbed flux density at 140 MHz for optically thin emission, and $\alpha$ is the optically thin spectral index (e.g. Weiler et al. 1986). We assume that $\alpha = -0.82$,
<table><caption>**Table 3.** Values of the best-fitting parameters to the radio spectra in Fig. 4. The units of *S* and *ν* are Jy and MHz, respectively.</caption><tbody><tr><td colspan="4">Free-free absorption fit:</td></tr><tr><td>S<sub>ν</sub> = S<sub>140,u</sub> (<sup>ν</sup>/<sub>140 MHz</sub>)<sup>-0.82</sup> exp [-(τ<sub>140</sub>) (<sup>ν</sup>/<sub>140 MHz</sub>)<sup>-2.1</sup>]</td><td>S<sub>140,u</sub> (Jy)</td><td>τ<sub>140</sub></td><td>Reduced χ<sup>2</sup></td></tr><tr><td>Eastern wing</td><td>99 ± 8</td><td>0.27 ± 0.03</td><td>2.57</td></tr><tr><td colspan="4">Linear fit:</td></tr><tr><td>log (S<sub>ν</sub>) = a + b log (ν)</td><td>a</td><td>b (≡ α)</td><td>Reduced χ<sup>2</sup></td></tr><tr><td>Central shell</td><td>3.46 ± 0.09</td><td>-0.58 ± 0.04</td><td>1.23</td></tr><tr><td>Western wing</td><td>2.62 ± 0.10</td><td>-0.48 ± 0.04</td><td>0.035</td></tr></tbody></table>

using the $α_{327.5}^{1465}$ value from D98. The turnover frequency (at about 120 MHz, where the fitted integrated flux density is approximately 80 Jy) is near the bottom of the LOFAR high band, although more low-frequency measurements are needed to confirm this. Further discussion can be found in Section 4.4.

## 4.2 W 50 spectral index map

We generated a two-point spectral index map for W 50 using the 140-MHz image shown in Figs 1 and 2, as well as the 1465-MHz image from D98; see Fig. 5. First, to account for the fact that the angular resolutions are slightly different, we convolved the 1465-MHz image (resolution 56 arcsec × 54 arcsec, BPA −61.4°; figs 1a and 1b in D98), to the resolution of our LOFAR map (78 arcsec × 55 arcsec, BPA 20°). Secondly, for consistency with the 1465-MHz map, we regridded the LOFAR image from 10 to 15 arcsec pixel⁻¹. Lastly, we calculated the spectral index values: log (S<sub>140</sub>/S<sub>1465</sub>)/log (140/1465) per pixel, with average background corrections of −4.5 and 3 mJy beam⁻¹ at 140 and 1465 MHz, respectively, applied first.

Generally, there is a very good agreement between the astrometry of the two maps, although there are some examples of non-systematic offsets of order 1–2 pixels for moderately bright, compact sources in the field of view. Investigating this further, it appears as though our astrometry is generally more consistent with, for example, the 1.4-GHz NVSS. Given the limitations in our analysis, which we describe below, we have not attempted to make differential alignment corrections to account for these relatively small offsets.

One difficulty with the method described above for calculating the spectral index map is that varying, non-zero background levels will cause biases; in particular, this is an issue for our LOFAR image. Hence, the results are best interpreted in a qualitative sense. As in the $α_{327.5}^{1465}$ spectral index map (fig. 4 in D98), we too observe regions of flatter spectral indices towards the south-west of the source, either side of the steeper, bright filament. However, in general, the spectral index gradients are not as pronounced in our map; moreover, while the apparent steepening from west to east across the source is consistent with the D98 spectral index map, this is, at least in part, likely to result from biases related to the background levels (see below). The brightest filamentary structure in the eastern wing may slightly flatten immediately to the east and west of the central region, although higher resolution maps are needed to confirm this. The very flat spectral indices along some of the boundary regions are artificial: the negative residual backgrounds at 140 MHz result in underestimated 140-MHz flux densities.

There is a hole-like feature approximately centred at RA 19h 10m 58s.1, Dec. +04° 39′ 57.8 (J2000), in the south-western region of the central shell below the bright filament, that has an inverted spectrum (maximum value ~0.2). In the LOFAR band, it is ~7–8 times fainter than the diffuse emission surrounding it; the feature is labelled in Fig. 2. In D98, Figs 1b and 2 suggest possibly similar behaviour at 327.5 and 1465 MHz, but not as pronounced in the latter case (also see fig. 1 in Farnes et al. 2017). This may be a significantly absorbed region, although a corresponding foreground H II region at or near these coordinates could not be found in the SIMBAD (Wenger et al. 2000) and VizieR (Ochsenbein, Bauer & Marcout 2000) online data bases.<fnref n="3" /> Moreover, Farnes et al. (2017) presented a continuum-corrected H α mosaic of the W 50 region, using data from the Isaac Newton Telescope Photometric H α Survey of the Northern Galactic Plane (IPHAS; Drew et al. 2005); however, an excess of emission at the position of the feature is not observed (fig. 6 in Farnes et al. 2017).

As mentioned above, generating a spectral index map using the standard formula, pixel by pixel, can potentially give highly misleading quantitative results. Therefore, following Anderson & Rudnick (1993), as well as D98 (also see e.g. Leahy & Roger 1991 and Leahy et al. 1998), we investigated a linear regression method. Various box sizes (5–35 pixels, i.e. 75 arcsec to 8.75 arcmin, per side) were slid over both maps in increments of one pixel in both directions. We then used the average slope from the regressions of S<sub>1465</sub> versus S<sub>140</sub> and S<sub>140</sub> versus S<sub>1465</sub> (taking the inverse of the calculated slope in the latter case) to determine the spectral index for each area, which was assigned to the central pixel. Regions off-source were masked out. Rather than considering the ratio of the flux densities at 140 and 1465 MHz, this approach, also often referred to as the T–T plot method (where T is the brightness temperature), removes any bias which is constant over the box size, but averages over the spectra of small-scale variations within the box (e.g. structure within filaments). Larger box sizes will have a range of faint and bright emission to provide more accurate constraints when determining each slope.

After using the technique described in the previous paragraph, the apparent $α_{140}^{1465}$ spectral index gradient visible in Fig. 5 is significantly reduced. Each of the three components of W 50 has about the same mean spectral index, ~−0.4, a value consistent with the two-point spectral indices $α_{140}^{1465}$ derived from the integrated flux densities in Table 2. This mean is not affected significantly by the choice of the dimensions of the box, apart from the typical 1σ error, which decreases from ~0.3 to ~0.1 with increasing box size. The effects of correlated noise, which can, for example, be reduced by only considering pixels at regularly spaced intervals (e.g. Green 1990), are second order here. Unlike between 327.5 and 1465 MHz, an approximately constant average spectral index $α_{140}^{1465}$ across W 50 simply follows from the fact that there is spectral turnover in the eastern wing at low frequencies (Table 2 and Fig. 4), and so $α_{140}^{1465}$ is significantly flatter than $α_{327.5}^{1465}$ for the eastern wing.

## 4.3 Spectral index asymmetry between the wings

D98 suggested that the density gradient due to the nearby Galactic plane can explain the difference in the radio spectra of the two wings between 327.5 and 1465 MHz: the western wing has a flatter spectrum because the jet is expanding into a denser region, leading

<footnote n="3">³ http://simbad.u-strasbg.fr/simbad/ and http://vizier.u-strasbg.fr/viz-bin/VizieR</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Two-point spectral index map between 140 and 1465 MHz. Pixels with flux densities $S_{1465} < 5$ mJy beam$^{-1}$ have been masked out. LOFAR contours have also been overlaid to outline the shape of the nebula; the lowest contour is at 30 mJy beam$^{-1}$ and each successive contour is a factor of two higher in surface brightness (up to a maximum of 1.92 Jy beam$^{-1}$). In addition, as in Fig. 2, we show the central shell. The pixel-by-pixel spectral indices across W50 range from about -0.9 to 0.2, but, to ensure a suitable contrast range in the plot, we use a subset of this range (indicated in the grey-scale bar). Average background levels of -4.5 and 3 mJy beam$^{-1}$ at 140 and 1465 MHz, respectively (see Section 3.1 and D98) were subtracted before the spectral index map was generated. However, note that the contours plotted are the raw values.</figcaption></figure>

to stronger shocks and higher levels of compression. However, our results in Section 4.2 are less consistent with this interpretation, albeit with the possible complication of absorption effects at low frequency. On the other hand, in Fig. 5, the steepening of the spectral index from west to east across W50 is a relatively smooth. Some of this is highly likely to be due to the corresponding gradient in the background level at 140 MHz, but we cannot determine with certainty whether there is a residual gradient after accounting for the varying background levels: the uncertainties in the spectral index maps generated with the linear regression method are too large.

We now explore an alternative explanation for any spectral index asymmetry between the wings. Expanding into the higher density environment closer to the Galactic plane, the western wing has been less affected by adiabatic expansion, and so its magnetic field will consequently be stronger than that in the more expanded eastern wing of the nebula. Electrons emitting synchrotron radiation at a specific frequency $\nu$ in an ambient magnetic field of strength $B$ have a Lorentz factor $\gamma$ given by

$$ \gamma = \left( \frac{2\pi m_e \nu}{eB} \right)^{1/2}, \qquad (3) $$

where $e$ is the elementary charge and $m_e$ is the electron rest mass. Therefore, in the western wing, where the magnetic field is stronger, we would probe lower Lorentz factors at a fixed radio frequency. Moreover, for a curved distribution of the underlying electron Lorentz factors, we should also measure a flatter spectral index (e.g. Scheuer & Williams 1968).

If this hypothesis is correct, then we may expect the equipartition magnetic field $B_{eq}$ to be stronger in the western wing. However, X-ray observations of W50 suggest similar ranges of $B_{eq}$ in the eastern and western wings: Safi-Harb & Ögelman (1997) calculated values of 7–60 µG (eastern wing) and 7–49 µG (western wing; also see Moldowan et al. 2005). The equipartition magnetic fields can also be estimated from radio data, although the calculation can involve significant uncertainties; revisions to the classical approach (e.g. Longair 1994) have been suggested (e.g. Beck & Krause 2005; Arbutina et al. 2012, 2013; Duffy & Blundell 2012). In particular,

Arbutina et al. (2012, 2013) developed modified equipartition calculations for the particular case of SNRs, although it is possible that further modifications would be needed for W50, given both its distinct morphology and the relatively flat spectral index of the western wing (their model is valid for $-1 \le \alpha \le -0.5$).

While it is beyond the scope of this paper to carry out a detailed comparison between the different methods, we now briefly consider the formalism developed by Duffy & Blundell (2012). Despite the fact that the focus of their study is the lobes of radio galaxies and quasars, the framework described is based upon fundamental synchrotron physics. The key quantities in equation 26 of Duffy & Blundell (2012) are as follows: (i) a spectral curvature term (that is intrinsic curvature and not from absorption effects); (ii) the frequency of the peak of the radio spectrum; and (iii) the emissivity at this frequency. All terms vary slowly, however, due to power-law dependencies of 1/7 or 2/7. We assume that the majority of the radio emission in the eastern wing can be modelled as a cylinder, height 47 arcmin and diameter 22 arcmin; similarly, we model the western wing as a cone with height 25 arcmin and base diameter 27 arcmin. Some straightforward calculations can be made to show that, very tentatively, and in a qualitative sense, $B_{eq}$ may be at least slightly stronger in the western wing if the associated spectral curvature is milder than for the eastern wing. However, the comparison will be affected by spectral curvature that is caused by foreground free-free absorption (Section 4.4).

In addition, despite the limitations mentioned above, given that the radio spectrum of the western wing could have a spectral index slightly steeper than -0.5 within the $1\sigma$ uncertainties (Table 3), then the online calculator<fnref n="4" /> for the Arbutina et al. (2012, 2013) model also possibly suggests a slightly stronger value of $B_{eq}$ for the western wing, if the non-spherical geometries of both wings are compensated for, to first order. Potentially significant differences in the depths and filling factors of the two wings are additional sources of uncertainty, both for this model and the Duffy & Blundell (2012) formalism. More complete radio spectra are also needed.

<footnote n="4">⁴ http://poincare.matf.bg.ac.rs/~arbo/eqp/</footnote>
## 4.4 Low-frequency spectral curvature

Below 327.5 MHz, it is the eastern wing and not the western wing that has a flatter spectrum; this is contrary to our suggestion in Section 4.3. One possibility is that the underlying electron energy spectra in the two wings are different, with perhaps shock acceleration in the denser ambient environment of the western wing accounting for the lack of spectral curvature, at least over the frequency range considered in this paper. Alternatively, the radio spectrum of the eastern wing may correspond to enhanced low-frequency free-free absorption along the line of sight to this component of the nebula. This hypothesis is supported by the fact that Farnes et al. (2017) found a global anticorrelation for W 50 between the linearly polarized flux density at 2.3 GHz and the strength of the diffuse H$\alpha$ emission, suggesting that this is due to a foreground Faraday screen of warm, ionized plasma (also see the polarization measurements in Downes, Salter & Pauls 1981b, and Downes, Pauls & Salter 1986). Inspecting the H$\alpha$ mosaic from this study (their fig. 6), bright emission is coincident with the position of the eastern wing.

In Section 4.1, we fitted the radio spectrum of the eastern wing with a free-free absorption model. Taking the fitted value $\tau_{140} = 0.27 \pm 0.03$ (Table 3), and using equation (1) with the assumption that $T_e \sim 10^4$ K, then the corresponding linear emission measure is EM $\sim 1.3 \times 10^4$ pc cm$^{-6}$. However, as mentioned before, more low-frequency data points are needed (e.g. within the LOFAR high band), particularly in light of the relatively poor reduced $\chi^2$ value of 2.57.

The potential effects of a synchrotron ageing break also need to be considered. Radiative losses from the emitting electrons result in a break in the synchrotron spectrum, and, for continuous particle injection, the spectrum is steepened from $v^\alpha$ to $v^{\alpha-0.5}$. The break frequency $v_B$ is given by

$$v_B = \frac{1610^2 \text{ GHz}}{(B/\mu\text{G})^3(t_s/\text{Myr})^2}, \quad (4)$$

where $t_s$ is the 'synchrotron age' of the distribution, i.e. the time since it was a power law out to infinite frequency (e.g. Carilli et al. 1991). If the age of W 50 is approximately 20 000 yr (Goodall et al. 2011a), then a break frequency of e.g. $\sim 100-300$ MHz would require a magnetic field of strength $\sim 2800-4000$ µG. These values are 1.7-2.8 dex larger than the equipartition magnetic fields that were discussed above in Section 4.3. Therefore, assuming that the lobes do not deviate significantly from equipartition, the effects of synchrotron ageing can be neglected.

In the 30.9-MHz Galactic plane survey by Kassim (1988), two thirds of a sample of 32 SNRs exhibit low-frequency spectral turnovers (Kassim 1989a,b). W 50 was detected, although it is not part of the study presented in Kassim (1989a,b); the eastern half of the source falls outside of the survey area, due to the Galactic latitude cut-off imposed. The brightest part of the western wing and the northern part of the central shell are clearly detected, with corresponding flux densities of 35.2 and 23.7 Jy, respectively (flux density uncertainty approximately 20 per cent; angular resolution 13.0 arcmin $\times 11.1$ arcmin). However, the more southerly regions of the nebula are at best marginally detected. Given that the 30.9-MHz flux densities are on the Baars et al. (1977) scale, we adjusted the 74-MHz flux densities in Table 2 back to this scale (dividing by a factor of 1.1; see Section 4.1), as well as making an additional correction for the 74-MHz flux density of the central shell: about 25 per cent of the flux density originates from the same spatial region as at 30.9 MHz (Miller-Jones, private communication). We therefore obtain approximate spectral indices $\alpha_{30,9}$ for the

western wing and central shell of $\sim 0.4$ and $\sim 1.0$, respectively. These spectral indices, coupled with the distinct change in morphology, provide further evidence for a low-frequency spectral turnover in W 50, possibly due to free-free absorption along the line of sight. LOFAR low-band measurements, which cover the frequency range 30-80 MHz, would allow a far more detailed investigation to be made.

Finally, we briefly consider whether a spectral turnover at low frequencies could also be related to a low-frequency cut-off in the electron energy spectrum. In this case we would expect a $S_v \propto v^{0.3}$ dependence (e.g. Scheuer & Williams 1968). Using equation (3) and assuming the equipartition magnetic field strengths calculated by Safi-Harb & Ögelman (1997), $\gamma \sim 420-1240$ at 30 MHz. It is not entirely clear whether a low-frequency cut-off would therefore be discernible with LOFAR, particularly given the upper end of this range (e.g. discussion in Harris 2005). Further measurements between 30 and 80 MHz would also potentially allow contributions from free-free absorption, a low-frequency cut-off and intrinsic spectral curvature to be disentangled.

# 5 LOW-FREQUENCY VARIABILITY OF SS 433

## 5.1 150-MHz light curve

A light curve of SS 433 from all 14 observations (Table 1) is shown in Fig. 6. For the deep run, we have not used the SS 433 flux density reported in Section 3. Instead, we made another averaged map from the same baseline range as that used when imaging the monitoring runs, and over a similar, more restricted uv coverage. We then divided by a correction factor of 1.2 to shift the effective average frequency from 140 to 150 MHz, bringing the flux density scale into line with the monitoring observations (see Section 2.3 for further details).

For the monitoring runs and the regenerated map from the first run described above, the angular resolution is about 140 arcsec $\times 100$ arcsec (median BPA $\approx 60^\circ$). Each Gaussian fit to the emission from SS 433 was point-like. Note that the new

<figure><img src="image_7.png" /><figcaption>Figure 6. LOFAR 150-MHz light curve for SS 433 over the period 2013 February to 2014 May, as well as the RATAN-600 4800-MHz light curve from 2012 November to 2014 May. We show $\pm 1\sigma$ error bars for each LOFAR data point; these are dominated by the $\pm 10$ per cent calibration uncertainty. The uncertainty for each 4800-MHz data point is $\pm 5$ per cent.</figcaption></figure>
measured flux density from the deep run (as given in Table 1) agrees with the value given in Section 3 within the 1σ uncertainties.

In Fig. 6, SS 433 shows tentative evidence for variability in the LOFAR band. To quantify its magnitude and significance, we used two statistics. First, the modulation index *m* is given by

$$m = \frac{\sigma_S}{\bar{S}}, \qquad (5)$$

where $\sigma_S$ is the standard deviation of the combined set of flux densities from Table 1, and $\bar{S}$ is the mean flux density. We find that *m* is 16 per cent, which is comparable to the modulation index at 235 MHz (21 per cent; Pandey et al. 2007). If we take into account the errors on each data point and debias the modulation index (as, for example, is described in Bell et al. 2014), then *m* is reduced to 11 per cent.

Secondly, we calculated the $\chi^2$ probability that the flux density remained constant over the period covered by our observations (e.g. Gaensler & Hunstead 2000; Bell et al. 2014). The quantity

$$x^2 = \sum_{i=1}^{n} \frac{(S_i - \bar{S})^2}{\sigma_i^2}, \qquad (6)$$

where $S_i$ is the ith flux density measurement, $\bar{S}$ is the weighted mean flux density, and $\sigma_i$ is the uncertainty for the ith flux density measurement, follows a $\chi^2$ distribution with $n-1$ degrees of freedom, assuming that the uncertainties are drawn from a normal distribution. The weighted mean is calculated using the formula

$$\tilde{S} = \sum_{i=1}^{n} \left( \frac{S_i}{\sigma_i^2} \right) / \sum_{i=1}^{n} \left( \frac{1}{\sigma_i^2} \right). \qquad (7)$$

The assumed 10 per cent calibration uncertainty is the dominant error term. Background corrections off-source, including from W50, are second-order effects because of the uv range chosen, and hence we do not consider the small contributions from the associated uncertainties.

Using equations (6) and (7) for the *n* = 14 observations in Table 1, we find that the probability *P* of $x^2$ exceeding our calculated value by chance [i.e. $P(x^2 > 28.3)$ for 13 degrees of freedom] is $8.2 \times 10^{-3}$. This probability value may be too large to be regarded as a significant indicator of variability (e.g. the criteria discussed in Gaensler & Hunstead 2000), but we suggest that it is sufficiently small that it could be interpreted as being at least marginally significant. Furthermore, in Fig. 6, there are hints of a general consistency in the calibration for several of the more closely spaced observations, suggesting that the calculated $\chi^2$ probability does not simply result from a light curve with large scatter, but rather genuine variability. Reducing the calibration uncertainty further would be needed to better quantitatively assess the magnitude of the intrinsic variability of the source.

## 5.2 Comparison with 4800-MHz data

In Fig. 6, we have also plotted 4800-MHz flux densities from contemporaneous monitoring observations (angular resolution 5 arcmin × 50 arcsec) carried out with the RATAN-600 telescope.<fnref n="5" /> The same broad trends are observed at both 150 and 4800 MHz: an initial decrease in flux density to a relatively constant baseline level (approximately 1.6 and 0.44 Jy at 150 and 4800 MHz, respectively), followed by increased activity at later epochs, before a second decline in flux density at the end of the sampling period.

### 5.2.1 Individual flares

At 4800 MHz, there were bright flares on MJD 56442.99 and 56724.22, which peaked at 1.40 and 1.28 Jy, respectively. However, neither of them have obvious counterparts at 150 MHz. There was also activity before and shortly after our LOFAR observing campaign began, with a maximum 4800-MHz flux density of 1.40 Jy on MJD 56264.48.

We now discuss whether we would have expected to see the outbursts in the LOFAR band. First, if the plasmons associated with a flare are well described by, for example, the van der Laan (1966) synchrotron bubble model (also see e.g. Hjellming & Johnston 1988, and Ball & Vlassis 1993), then the rise time of an outburst is longer at lower frequencies (owing to opacity effects), and the peak flux density is often lower (relative to the quiescent baseline level). The ratio of the peak flux densities at 150 and 4800 MHz is given by

$$\frac{S_{\text{max},150}}{S_{\text{max},4800}} = \left( \frac{4800}{150} \right)^{\frac{-(7p+3)}{(4p+6)}}, \qquad (8)$$

where $p = 1 - 2\alpha$ is the negative of the power-law index of the electron energy (*E*) distribution (i.e. $N(E) dE \propto E^{-p} dE$). We use a fiducial value of 1 Jy for the 4800-MHz flares (after subtracting the quiescent component), and, for simplicity, assume that the spectral index is in the range -1 and 0, which is broadly consistent with long-term monitoring data between 960 MHz and 21.7 GHz from RATAN-600 (Trushkin et al. 2003). Then, for the resulting values of *p* (1–3), we would expect the peak flux density at 150 MHz to be about 1–3 per cent of the 4800-MHz value, i.e. approximately 10–30 mJy. This is much less than the typical calibration uncertainty (Table 1). Conversely, a hypothetical flare of, say, 500 mJy at 150 MHz, marginally detectable given our current calibration uncertainty, would be ~tens of Jy at 4800 MHz. Such a bright outburst has not been observed before from SS 433.

An alternative approach is as follows. First, Trushkin et al. (2003) empirically derived an average power-law relationship between the maximum flux density $S_{max}$ of a flare (after the quiescent component has been subtracted) and the observing frequency $\nu$, such that

$$S_{\text{max}} = 1.3\nu^{-0.46}, \qquad (9)$$

where $S_{max}$ is in units of Jy, and $\nu$ is in GHz. Equation (9) is inconsistent with the predictions of the van der Laan (1966) model for the values of *p* that we considered above; this equation only agrees with equation (8) for a significantly inverted spectral index $\alpha \approx 0.83$. However, a negative power-law dependence is, for example, expected in the model developed by Martí, Paredes & Estalella (1992), at frequencies where the outburst is initially optically thin; this model assumes particle injection into twin jets that are adiabatically expanding in the lateral direction, exponentially at first.

From equation (9), a hypothetical 1-Jy flare at 4800 MHz is brighter than what we would expect on average at this frequency (0.63 Jy). None the less, assuming that the average power-law index still holds, a 1-Jy flare at 4800 MHz would be about 2 Jy at 960 MHz, the lowest RATAN-600 frequency used in Trushkin et al. (2003). Moreover, in the absence of a spectral turnover, an extrapolation to lower frequencies would suggest a flux density of around 5 Jy in the LOFAR high band. Despite our limited sampling, there is no evidence of such a large flux density increase at 150 MHz in Fig. 6.

However, if we assume that a flare from SS 433 will initially be optically thick below ~1 GHz, then we can use equations (29) and (30) in Martí et al. (1992) to estimate the flux density ratio of the flare between 150 and 960 MHz. These equations are

$$S_{\text{max}} \propto \nu^{(2.78p+3.46)/9} \qquad (10)$$

<footnote n="5">⁵ http://www.sao.ru/cats/cgi-bin/ss433.cgi</footnote>
in the case of free-free absorption, and

$$S_{\max} \propto \nu^{(13p+2)/(7p+8)} \quad (11)$$

for synchrotron self-absorption. For the value of p corresponding to equation (9), i.e. $p = 1.92$, we would therefore expect the ~2-Jy flare at 960 MHz to have a peak flux density of either ~340 or ~200 mJy at 150 MHz. Repeating the entire alternative approach for $p = 1-3$, including modifying the spectral index from equation (9) accordingly, the expected peak flux density at 150 MHz varies between ~160 and 440 mJy for a 1-Jy flare at 4800 MHz. The flare, although now brighter than was predicted using equation (8), would still not be detectable with the current calibration uncertainty.

Also note that as a rough consistency check, if one were to start with the average maximum flux density at 960 MHz from equation (9), i.e. 1.3 Jy, and do the same analysis as described in the previous paragraph, then the expected peak flux density at 408 MHz would be either ~560 or ~440 mJy. On average, flares at 408 MHz have an excess flux density of $1110 \pm 410$ mJy above the quiescent baseline level (Bonsignori-Facondi et al. 1986), which overlaps with our calculated values within at most $1.6\sigma$ (also see Vermeulen et al. 1993). If the initial transition between optically thin and optically thick flares occurs below ~1 GHz, then the flares at 408 MHz would generally be brighter than predicted above using equations (10) and (11), and the chances of a detection at 150 MHz would also increase.

The fact that the 150- and 4800-MHz observations were independent observing programmes also needs to be considered. Trushkin et al. (2003) empirically derived a power-law relationship between the time at which the maximum flux density occurs during a flare, $t_{Smax}$, and the observing frequency, such that

$$t_{Smax} = 5.1\nu^{-0.43}, \quad (12)$$

where $t_{Smax}$ is in units of days, and $\nu$ is in GHz. Equation (12) suggests a delay of about 9 d between 4800 and 150 MHz, assuming that it is still valid at low frequencies. The first LOFAR observation after the 4800-MHz flare on MJD 56442.99 occurred 56 d afterwards, while there was a LOFAR observation 4.1 d after the flare on MJD 56724.22. While in theory we may have been able to detect the flare rise with LOFAR in the latter case, the extended 4800-MHz activity before this particular flare may complicate the analysis at low frequencies (Section 5.2.2). Closely coordinating observations across a wide range of radio frequencies will be the most robust approach for studying future flaring activity from SS 433.

### 5.2.2 Extended flaring activity

The onset of 4800-MHz activity at MJD ~56560 is accompanied by a sustained ~0.5-1 Jy rise in flux in the LOFAR band. There is also a hint of a similar flux density increase when comparing the 150- and 4800-MHz data taken at times up to and including the epoch corresponding to the second LOFAR data point. It is possible that the numerous separate high-frequency bursts have become blended together at 150 MHz, owing to the longer associated rise timescales at low frequencies. This is also very tentatively suggested by the apparent difference between the 150-MHz flux densities at the start and end of the 4800-MHz flaring activity beginning at MJD ~56560: approximately 1.63 and 1.89 Jy, respectively.

The median peak brightness of these flares at 4800 MHz is about 360 mJy, after correcting for the quiescent levels on either side of the period of activity. A few tens of flares, spaced ~5 d apart, would not replicate the observed rise in the LOFAR flux density if the van

der Laan (1966) model is assumed: a plateau can be modelled, but it is over an order of magnitude too faint and decays too rapidly. Although beyond the scope of this paper, if individual flares are instead much brighter at 150 MHz, as, for example, was suggested above in Section 5.2.1 when considering the Martí et al. (1992) model, then a sustained rise of order several hundred mJy or more is a possibility. More complex scenarios previously outlined in the literature (e.g. Vermeulen et al. 1993; Blundell, Schmidtbreick & Trushkin 2011; Jeffrey et al. 2016) may also be relevant.

### 5.2.3 Possible effects of refractive interstellar scintillation

Another potential explanation for the relatively slow variability that we observe at 150 MHz is refractive interstellar scintillation (RISS; e.g. Rickett 1986). Pandey et al. (2007) investigated whether RISS could explain the variability of SS 433 in their GMRT data at 235 and 610 MHz, concluding that the source is most likely affected by this phenomenon. Taking their calculated time-scales for RISS and making a standard assumption that this time-scale is proportional to $\lambda^{2.2}$ (where Kolmogorov turbulence is assumed), then the associated time-scale at 150 MHz is ~2.3 yr. This value is roughly double the length of our monitoring campaign, and so it seems unlikely that we are seeing a corresponding effect in our data. A much longer time baseline at 150 MHz would be beneficial.

## 5.3 Quiescent flux densities

Considering the data points between MJD 56407 and 56560, before the onset of the extended flaring activity, the median 4800-MHz flux density, 0.44 Jy, is close to the long-term average quiescent value of 0.43 Jy (Trushkin et al. 2003; also note that the median flux density at the end of the flaring activity, 0.40 Jy, is very similar). In addition, the average flux density of the five LOFAR observations over this period is $1.63 \pm 0.07$ Jy, which is consistent with both 160 MHz quiescent flux density measurements from Seaquist et al. (1980): $1.6 \pm 0.2$ Jy (1979 July 1) and $2.2 \pm 0.3$ Jy (1979 July 21 and 22). These data were obtained at a similar angular resolution to the LOFAR observations ($2.3$ arcmin $\times$ $1.9$ arcmin; BPA $0^{\circ}$) and also have an absolute flux density scale that is consistent with our data (Roger et al. 1973; Seaquist et al. 1980; Scaife & Heald 2012).

Seaquist et al. (1980) suggested that the quiescent radio spectrum of SS 433 turns over near 300 MHz, although there was a gap in their data between 160 and 408 MHz. In Table 4, we list the average flux densities from this study, along with more recent measurements, including our LOFAR data. We have estimated the quiescent flux densities at 235 and 610 MHz from the monitoring data presented by Pandey et al. (2007); given the possibility of flaring activity and/or the effects of interstellar scintillation during this observing campaign (discussion in Pandey et al. 2007), our estimate in each case is the range between the minimum flux density and the median flux density.

Fig. 7 shows the quiescent radio spectrum at low and mid frequencies. We have also plotted an extrapolated empirical power-law fit that was derived by Trushkin et al. (2003) using higher frequency RATAN-600 data, where the spectral index $\alpha = -0.60 \pm 0.14$. The spectrum still appears to be consistent with a turnover at around 300 MHz, as suggested by Seaquist et al. (1980), and our LOFAR data further suggest that the spectrum has truly turned over at 150-160 MHz. However, there may be discrepancies resulting from the different angular resolutions of the various data sets, which would affect the relative level of contamination from W 50. Again,
<table><caption>**Table 4.** Quiescent flux density estimates for SS 433 at low and mid frequencies. We have calculated average values at 80 and 160 MHz based on the individual measurements presented in Seaquist et al. (1980). In addition, the ranges we have calculated at 235 and 610 MHz are based on the minimum and median flux densities in Pandey et al. (2007). The 843-MHz value was determined from the tabulated flux densities in Vermeulen et al. (1993). Note that the upper limit at 80 MHz is approximately consistent with the 74-MHz flux density from Miller-Jones et al. (2007), although there are not enough higher frequency data available to fully establish the extent of activity from SS 433 when the 74-MHz data were obtained (Miller-Jones, private communication).</caption><thead><tr><th>Frequency (MHz)</th><th>S<sub>SS433</sub> (Jy)</th><th>Reference(s)</th></tr></thead><tbody><tr><td>80</td><td>&lt;1.6<sup>a</sup></td><td>(1)</td></tr><tr><td>150</td><td>1.63 ± 0.07</td><td>(2)</td></tr><tr><td>160</td><td>1.9 ± 0.2</td><td>(1)</td></tr><tr><td>235</td><td>1.73–2.74</td><td>(3)</td></tr><tr><td>408</td><td>1.1–2.6</td><td>(4.1)</td></tr><tr><td>408</td><td>2.30 ± 0.32</td><td>(5)</td></tr><tr><td>610</td><td>0.73–1.23</td><td>(3)</td></tr><tr><td>843</td><td>1.13 ± 0.09</td><td>(6)</td></tr></tbody></table>

References: (1) Seaquist et al. (1980); (2) this paper; (3) Pandey et al. (2007); (4) Spencer (1979); (5) Bonsignori-Facondi et al. (1986); (6) Vermeulen et al. (1993).

<sup>a</sup>5σ upper limit.

<figure><img src="image_8.png" /><figcaption>**Figure 7.** The low- and mid-frequency radio spectrum for SS 433 during quiescence, using the flux densities from Table 4. The dashed line is the RATAN-600 quiescent flux density empirical power-law fit from Trushkin et al. (2003), i.e. $S_v = 1.1v^{-0.6}$ (where $S_v$ and $v$ are in units of Jy and GHz, respectively), extrapolated to lower frequencies. Error bars indicate ±1σ uncertainties, the solid lines between asterisks indicate flux density ranges at 235, 408, and 610 MHz, and the triangle represents the average 5σ upper limit at 80 MHz.</figcaption></figure>

in-band LOFAR flux densities could be used to obtain a more complete picture of the low-frequency properties of the quiescent radio spectrum of SS 433. It may also be particularly useful to obtain new HBA measurements at the upper end of the operating frequency range (210–240 MHz) and GMRT data at 325 MHz.

## 5.4 Variability of nearby sources

We searched for transients, as well as other variable sources in the field, by running all of our LOFAR maps through the TRANSIENTS PIPELINE (TRAP; Swinbank et al. 2015). No transients were found over the course of our monitoring campaign. Moreover, SS 433 has the most significant variability statistics of all the detected sources. Given the discussion in Section 5.1, our data are therefore not sufficiently precise to allow us to draw conclusions concerning the potential low-frequency variability of other sources in the field.

D98 detected a ~340-mJy point source at 327.5 MHz, close to the bright filaments in the eastern wing of W 50 (coordinates RA 19h 15m 49.8°, Dec. +04° 51' 45"; J2000). We have labelled the position of this source with a white cross in Fig. 2. Given several non-detections or marginal detections at higher frequencies, D98 suggested that it may be variable. We do not detect this source in any of our LOFAR data sets (e.g. 3σ upper limit ~25 mJy beam⁻¹ in the deep observation, after taking into account the background level from W 50). Furthermore, the source appears not to be present in the Farnes et al. (2017) ATCA map. There are also non-detections in four 327.5-MHz VLA maps over the period 2001 January–October, as well as the 1996 August 19 run used in the D98 study; the 3σ upper limit from the map generated using all five observations is ~17 mJy beam⁻¹ (Miller-Jones, private communication). This result suggests that the source may have only been present in the second 327.5-MHz data set that D98 obtained on 1997 August 18. Another possibility may be that the source is an artefact generated during the imaging process described in D98, which involved source subtraction and re-addition steps; see section 2.2 in D98.

## 6 CANDIDATE SNR G 38.7–1.4

In Fig. 1, west-south-west of SS 433 and W 50, there is a diffuse radio source that is coincident with the candidate SNR G 38.7–1.4 (also known as SNR G 38.7–1.3). This candidate was first suggested by Schaudel et al. (2002), who searched the ROSAT all-sky survey (Voges et al. 1999) for extended X-ray sources that could be unidentified SNRs. They found a 12 arcmin × 8 arcmin source that is coincident with an incomplete radio shell observed in the Effelsberg Galactic Plane Survey at 11 cm (Reich et al. 1984, 1990), favouring an SNR interpretation. An independent identification was also made at optical wavelengths in the IPHAS Hα survey by Sabin et al. (2013), who further supported their claims through comparisons with the same radio shell, as seen in other archival radio data sets.

Additional evidence for an SNR identification was found by Huang et al. (2014). Using XMM-Newton and Chandra X-ray Observatory observations, these authors suggested that the X-ray emission is from shock-heated plasma. In addition, their analysis of archival VLA radio observations at 1.4 GHz revealed a coincident broken radio shell, consistent with that reported by Schaudel et al. (2002). Huang et al. (2014) estimated that G 38.7–1.4 is at a distance of ~4 kpc, with an age of ~14 000 yr, and also suggested that the source belongs to the SNR mixed-morphology category (Rho & Petre 1998).

Although previously detected as a broken radio shell, our deep observation of SS 433 and W 50 has sufficient surface-brightness sensitivity and angular resolution to reveal a more complete radio shell morphology (left-hand panel of Fig. 8). In particular, the candidate SNR is significantly larger than had been indicated by the X-ray observations, with a diameter of approximately
<figure><img src="image_9.png" /></figure>

**Figure 8.** *Left:* LOFAR detection of the radio shell of G 38.7–1.4. The north-western (NW) and south-eastern (SE) components of the shell are labelled to assist the discussion in the text, and the knot from which a spectral index was estimated is circled. The positions of PSR J1904+0451 and PSR J1906+0414 are indicated by crosses. An average background level correction of –40 mJy beam–1 has been applied to the map. *Right:* MSX 8-µm image of the G 38.7–1.4 region, overlaid with LOFAR contours indicating surface-brightness levels of 30, 40, 60, 80, 100, 200, and 400 mJy beam–1 (also background corrected). As in the left-hand panel, the knot from which a spectral index was estimated is circled. There is no infrared counterpart to the nebula. Note that compared to the left-hand panel, the right-hand panel shows a slightly more zoomed-in view on G 38.7–1.4.

40 arcmin × 30 arcmin, centred at Galactic coordinates (l, b) = (38.7, −1.1). However, accurate flux density measurements are hampered by the fact that the source sits in a significant negative bowl. After making a first-order background correction and assuming that the compact sources coincident with the shell are unrelated, the average 140-MHz surface brightness of the diffuse emission is ~30 mJy beam−1.

Traditionally, the confirmation of an SNR is achieved by demonstrating that the radio emission is non-thermal and/or linearly polarized. A polarimetric analysis of our LOFAR data set was not possible at the time of writing. Also, while the south-eastern part of the shell (labelled ‘SE’ in the left-hand panel of Fig. 8) falls within the ATCA 2.3-GHz mosaic from the Farnes et al. (2017) study, no associated linearly polarized emission is detected, which is suggested to be due to foreground depolarization. Therefore, the potential non-thermal nature of G 38.7–1.4 can instead be tested by using the selection criteria outlined in Brogan et al. (2006): other than having a shell-like or partial-shell-like morphology, the candidate must have a negative spectral index (where $S_ν ∝ ν^{α}$, as before) and must show a lack of associated mid-infrared 8-µm emission.

The VLA Galactic Plane Survey (VGPS; Stil et al. 2006), conducted at 1420 MHz with 60-arcsec angular resolution, covers the region containing G 38.7–1.4. While only a hint of the south-eastern part of the shell is detected in the VGPS map, the north-western part of the shell (labelled ‘NW’ in the left-hand panel of Fig. 8) is slightly brighter. Therefore, in order to calculate a rough spectral index for G 38.7–1.4, we first measured the flux density from the brightest knot of diffuse radio emission that we assume to be associated with the north-western part of the shell, indicated in the left-hand panel of Fig. 8. Our estimate is $\alpha_{140}^{1420} \sim -0.8 \pm 0.1$, where the $1\sigma$ uncertainty is dominated by the uncertainty in the background subtraction in the LOFAR map. The fainter diffuse emission in the vicinity of this feature has a spectral index ranging from about –1.0 to –0.5, although the $1\sigma$ spectral index uncertainty can be as large as about 0.3. Moreover, using a $3\sigma$ upper limit at 1420 MHz, an approximate spectral index constraint for the south-eastern part of the shell is $\alpha_{140}^{1420} \lesssim -0.6$. Therefore, given these various values,

there is no significant evidence to suggest a discrepancy in relation to the general SNR spectral index distribution (e.g. $\bar{\alpha} \approx -0.48$ and $\sigma_\alpha \approx 0.15$ in the Green 2014 SNR catalogue), and, in addition, our measurements are indicative of non-thermal radiation.

Referring to the third selection criterion from Brogan et al. (2006), at the position of G 38.7–1.4, there is no evidence for associated mid-infrared emission at 8 µm, using Galactic plane survey data collected by the Midcourse Space Experiment (MSX; Price et al. 2001). The MSX map is shown in the right-hand panel of Fig. 8. This result further suggests a non-thermal nature for the shell, strengthening the potential SNR identification.

We also investigated the possibility of nearby pulsar associations by searching the Australia Telescope National Facility Pulsar Catalogue (Manchester et al. 2005).<fnref n="6" /> The two closest known pulsars are PSR J1904+0451 (angular distance ≈21 arcmin from the centre of the nebula), with a 1400-MHz flux density of 0.117 mJy, and PSR J1906+0414 (angular distance ≈26 arcmin), with a 1400-MHz flux density of 0.23 mJy. Their positions are indicated in the left-hand panel of Fig. 8. At 140 MHz, the $3\sigma$ flux density upper limits are ~15–20 mJy beam−1 (background corrected). PSR J1904+0451 is an isolated millisecond pulsar with a spin period of 6.09 ms; the dispersion measure distance and characteristic age are 4.7 kpc (using the Cordes & Lazio 2002 model) and $1.7 \times 10^{10}$ yr, respectively (Scholz et al. 2015a,b). PSR J1906+0414 is also not part of a binary system and has a spin period of 1.04 s; furthermore, the dispersion measure distance and characteristic age are 7.5 kpc and $1.4 \times 10^6$ yr, respectively (Lorimer et al. 2006). While there is certainly a very large discrepancy between the characteristic ages of the two pulsars and the estimated age of the SNR candidate, the potentially significant uncertainties associated with the various age and distance estimates do not allow any strong conclusions concerning a pulsar-SNR association to be made (e.g. following the guidelines in Kaspi 1996). Moreover, assuming an age of roughly 14 000 yr, an association in

<footnote n="6">⁶ http://www.atnf.csiro.au/people/pulsar/psrcat/</footnote>
either case would imply that the pulsar has a projected velocity in the plane of the sky of up to several thousand km s⁻¹, well into the tail of the velocity distribution for isolated pulsars (e.g. Arzoumanian, Chernoff & Cordes 2002).

Brogan et al. (2006) demonstrated through 330-MHz VLA observations of the Galactic plane that deep, mid-frequency surveys can significantly increase the number of known SNRs, helping to solve the long-standing problem of ‘missing’ SNRs in the Galaxy (e.g. Helfand et al. 1989; Brogan et al. 2004). Our detection of G 38.7−1.4 with LOFAR is an encouraging result in this sense: although this source was previously known, we have shown that low-frequency observations with competitive low-surface-brightness sensitivity can also aid in the detection of diffuse emission. This is particularly so given that the surface brightness values will generally increase with decreasing frequency for the SNR spectral index distribution stated above. Galactic plane surveys conducted with the current generation of low-frequency radio telescopes therefore have the potential to reveal a hitherto missing population.

## 7 CONCLUSIONS

In this paper, we have presented the results from a LOFAR high-band observing campaign of the X-ray binary SS 433 and its associated SNR W 50. We have drawn the following conclusions from our study:

(i) Our 140-MHz wide-field image of SS 433 and W 50 is the deepest low-frequency map of this system to date. The morphology of W 50, as well as the flux densities of the various components, are consistent with previous investigations. Foreground free-free absorption may explain the curvature of the radio spectrum of the eastern wing, although more low-frequency data points are needed to confirm this result.

(ii) We have tentatively detected variability from SS 433, with both a debiased modulation index of 11 per cent and a $\chi^2$ probability of a flat light curve of $8.2 \times 10^{-3}$. The ~0.5–1 Jy rise in the 150-MHz flux density may correspond to extended flaring activity over a period of about 6 months, observed at 4800 MHz in contemporaneous RATAN-600 observations. The flux density increase is too large to explain with a standard synchrotron bubble model, although a model along the lines of the one presented in Martí et al. (1992), involving laterally expanding twin jets, may potentially offer an alternative explanation. While LOFAR could be a key trigger for other facilities if there is a prolonged period of activity, low-frequency detections of individual flares will require more accurate calibration, as well as higher cadence sampling.

(iii) Within the large field of view, a significant amount of extended structure has been detected near and along the Galactic plane. We have described the most complete detection thus far of the radio shell of SNR candidate G 38.7−1.4.

Our findings demonstrate the potential of LOFAR, as well as other low-frequency Square Kilometre Array precursors and pathfinders, for studying not only X-ray binaries, but SNRs as well.

## ACKNOWLEDGEMENTS

We thank the anonymous referee for their insightful, detailed comments, which were valuable in refining our analysis and presentation. We also acknowledge the considerable efforts of the ASTRON Radio Observatory, particularly Carmen Toribio, in setting up the observations and pre-processing the data. In addition, we warmly thank George Heald and Christian Kaiser for useful discussions, as

well as Gloria Dubner for kindly providing us with a 1465-MHz map of SS 433 and W 50.

JWB, RPF, AJS, GEA, TDS, and MP acknowledge support from European Research Council Advanced Grant 267697 ‘4 Pi Sky: Extreme Astrophysics with Revolutionary Radio Telescopes’. AR, JDS, AJvdH, DC, GJM, and RAMJW acknowledge support from European Research Council Advanced Grant 247295 ‘AART-FAAC’. JCAM-J is the recipient of an Australian Research Council Future Fellowship (FT140101082). SAT is grateful to the Russian Base Research Foundation for financial support (grant number 12-02-0812a) and acknowledges support through the Russian Government Program of Competitive Growth of Kazan Federal University. RPB acknowledges support from the European Research Council under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 715051; ‘Spiders’). SC acknowledges financial support from the UnivEarthS Labex program of Sorbonne Paris Cité (ANR-10-LABX-0023 and ANR-11-IDEX-0005-02). JWTH acknowledges funding from an NWO Vidi fellowship and European Research Council Starting Grant ‘DRAG-NET’ (337062). JvL acknowledges funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement n. 617199. The financial assistance of the South African SKA Project (SKA SA) towards this research is hereby acknowledged. Opinions expressed and conclusions arrived at are those of the authors and are not necessarily to be attributed to the SKA SA.

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 facilities in several countries, that are owned by various parties (each with their own funding sources), and that are collectively operated by the ILT foundation under a joint scientific policy.

This research has made use of the SIMBAD data base and the VizieR catalogue access tool, operated at CDS, Strasbourg, France. This research has also made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

This research made use of data products from the *MSX*. Processing of the data was funded by the Ballistic Missile Defense Organization with additional support from NASA Office of Space Science. This research has also made use of the NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

## REFERENCES

Abell G. O., Margon B., 1979, Nature, 279, 701

Altenhoff W. J., Downes D., Goad L., Maxwell A., Rinehart R., 1970, A&AS, 1, 319

Alves M. I. R., Davies R. D., Dickinson C., Davis R. J., Auld R. R., Calabretta M., Staveley-Smith L., 2010, MNRAS, 405, 1654

Alves M. I. R., Davies R. D., Dickinson C., Calabretta M., Davis R., Staveley-Smith L., 2012, MNRAS, 422, 2429

Anderson M. C., Rudnick L., 1993, ApJ, 408, 514

Arbutina B., Urošević D., Andjelić M. M., Pavlović M. Z., Vukotić B., 2012, ApJ, 746, 79

Arbutina B., Urošević D., Vučetić M. M., Pavlović M. Z., Vukotić B., 2013, ApJ, 777, 31
Arzoumanian Z., Chernoff D. F., Cordes J. M., 2002, ApJ, 568, 289

Baars J. W. M., Genzel R., Pauliny-Toth I. I. K., Witzel A., 1977, A&A, 61, 99

Ball L., Vlassis M., 1993, PASA, 10, 342

Beck R., Krause M., 2005, Astron. Nachr., 326, 414

Begelman M. C., Hatchett S. P., McKee C. F., Sarazin C. L., Arons J., 1980, ApJ, 238, 722

Bell M. E. et al., 2014, MNRAS, 438, 352

Black C. S., Fesen R. A., 2015, MNRAS, 447, 2540

Blundell K. M., Bowler M. G., 2004, ApJ, 616, L159

Blundell K. M., Mioduszewski A. J., Muxlow T. W. B., Podsiadlowski P., Rupen M. P., 2001, ApJ, 562, L79

Blundell K. M., Schmidtobreick L., Trushkin S., 2011, MNRAS, 417, 2401

Bonsignori-Facondi S. R., Padrielli L., Montebugnoli S., Barbieri R., 1986, A&A, 166, 157

Brand J., Blitz L., 1993, A&A, 275, 67

Braun R., Walterbos R. A. M., 1985, A&A, 143, 307

Bridle A. H., Schwab F. R., 1999, in Taylor G. B., Carilli C. L., Perley R. A., eds, ASP Conf. Ser. Vol. 180, Synthesis Imaging in Radio Astronomy II. Astron. Soc. Pac., San Francisco, p. 371

Briggs D. S., 1995, PhD thesis, New Mexico Institute of Mining and Technology

Brogan C. L., Devine K. E., Lazio T. J., Kassim N. E., Tam C. R., Brisken W. F., Dyer K. K., Roberts M. S. E., 2004, AJ, 127, 355

Brogan C. L., Gelfand J. D., Gaensler B. M., Kassim N. E., Lazio T. J. W., 2006, ApJ, 639, L25

Carilli C. L., Perley R. A., Dreher J. W., Leahy J. P., 1991, ApJ, 383, 554

Clark B. G., 1980, A&A, 89, 377

Cohen A. S., Lane W. M., Cotton W. D., Kassim N. E., Lazio T. J. W., Perley R. A., Condon J. J., Erickson W. C., 2007, AJ, 134, 1245

Condon J. J., Cotton W. D., Greisen E. W., Yin Q. F., Perley R. A., Taylor G. B., Broderick J. J., 1998, AJ, 115, 1693

Copetti M. V. F., Schmidt A. A., 1991, MNRAS, 250, 127

Cordes J. M., Lazio T. J. W., 2002, preprint (arXiv:astro-ph/0207156)

Crampton D., Georgelin Y. M., Georgelin Y. P., 1978, A&A, 66, 1

Dabbech A., Ferrari C., Mary D., Slezak E., Smirnov O., Kenyon J., 2015, A&A, 576, A7

Downes A. J. B., Salter C. J., Pauls T., 1981a, A&A, 97, 296

Downes A. J. B., Salter C. J., Pauls T., 1981b, A&A, 103, 277

Downes A. J. B., Pauls T., Salter C. J., 1986, MNRAS, 218, 393

Drew J. E. et al., 2005, MNRAS, 362, 753

Dubner G. M., Holdaway M., Goss W. M., Mirabel I. F., 1998, AJ, 116, 1842 (D98)

Duffy P., Blundell K. M., 2012, MNRAS, 421, 108

Eikenberry S. S., Cameron P. B., Fierce B. W., Kull D. M., Dror D. H., Houck J. R., Margon B., 2001, ApJ, 561, 1027

Ellingson S. W. et al., 2013, IEEE Trans. Ant. Prop., 61, 2540

Fabian A. C., Rees M. J., 1979, MNRAS, 187, 13P

Farnes J. S., Gaensler B. M., Purcell C., Sun X. H., Haverkorn M., Lenc E., O'Sullivan S. P., Akahori T., 2017, MNRAS, 467, 4777

Fesen R. A., Gull T. R., 1986, ApJ, 306, 259

Fiedler R. L. et al., 1987, AJ, 94, 1244

Forbes D., 1989, A&AS, 77, 439

Gaensler B. M., Hunstead R. W., 2000, PASA, 17, 72

Gao X. Y., Han J. L., Reich W., Reich P., Sun X. H., Xiao L., 2011, A&A, 529, A159

Garsden H. et al., 2015, A&A, 575, A90

Geldzahler B. J., Pauls T., Salter C. J., 1980, A&A, 84, 237

Girard J. N., Garsden H., Starck J. L., Corbel S., Woiselle A., Tasse C., McKean J. P., Bobin J., 2015, JINST, 10, C08013

Goodall P. T., Alouani-Bibi F., Blundell K. M., 2011a, MNRAS, 414, 2838

Goodall P. T., Blundell K. M., Bell Burnell S. J., 2011b, MNRAS, 414, 2828

Green D. A., 1990, AJ, 100, 1927

Green D. A., 2014, Bull. Astron. Soc. India, 42, 47

Guzmán A. E., May J., Alvarez H., Maeda K., 2011, A&A, 525, A138

Harris D. E., 2005, in Kassim N., Perez M., Junor M., Henning P., eds, ASP Conf. Ser. Vol. 345, From Clark Lake to the Long Wavelength Array: Bill Erickson's Radio Science. Astron. Soc. Pac., San Francisco, p. 254

Haslam C. G. T., Salter C. J., Stoffel H., Wilson W. E., 1982, A&AS, 47, 1

Heald G. H. et al., 2015, A&A, 582, A123

Helfand D. J., Velusamy T., Becker R. H., Lockman F. J., 1989, ApJ, 341, 151

Hjellming R. M., Johnston K. J., 1981a, Nature, 290, 100

Hjellming R. M., Johnston K. J., 1981b, ApJ, 246, L141

Hjellming R. M., Johnston K. J., 1988, ApJ, 328, 600

Holden D. J., Caswell J. L., 1969, MNRAS, 143, 407

Huang R. H. H., Wu J. H. K., Hui C. Y., Seo K. A., Trepl L., Kong A. K. H., 2014, ApJ, 785, 118

Hurley-Walker N. et al., 2017, MNRAS, 464, 1146

Jeffrey R. M., Blundell K. M., Trushkin S. A., Mioduszewski A. J., 2016, MNRAS, 461, 312

Kaspi V. M., 1996, in Johnston S., Walker M. A., Bailes M., eds, ASP Conf. Ser. Vol. 105, Pulsars: Problems and Progress. Astron. Soc. Pac., San Francisco, p. 375

Kassim N. E., 1988, ApJS, 68, 715

Kassim N. E., 1989a, ApJ, 347, 915

Kassim N. E., 1989b, ApJS, 71, 799

Kassim N. E., Weiler K. W., Erickson W. C., Wilson T. L., 1989, ApJ, 338, 152

Kovalenko A. V., Pynzar' A. V., Udal'Tsov V. A., 1994, AZh, 71, 110

Kuchar T. A., Clark F. O., 1997, ApJ, 488, 224

Kurtz S., 2005, in Cesaroni R., Felli M., Churchwell E., Walmsley M., eds, Proc. IAU Symp. 227, Massive Star Birth: A Crossroads of Astrophysics. Cambridge University Press, Cambridge, p. 111

Lane W. M., Cotton W. D., van Velzen S., Clarke T. E., Kassim N. E., Helmboldt J. F., Lazio T. J. W., Cohen A. S., 2014, MNRAS, 440, 327

Leahy D. A., Roger R. S., 1991, AJ, 101, 1033

Leahy D. A., Xizhen Z., Xinji W., Jiale L., 1998, A&A, 339, 601

Liebert J., Angel J. R. P., Hege E. K., Martin P. G., Blair W. P., 1979, Nature, 279, 384

Lockman F. J., Blundell K. M., Goss W. M., 2007, MNRAS, 381, 881

Longair M. S., 1994, High Energy Astrophysics. Vol. 2: Stars, The Galaxy and The Interstellar Medium, 2nd edn. Cambridge University Press, Cambridge

Lorimer D. R. et al., 2006, MNRAS, 372, 777

Manchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005, AJ, 129, 1993

Margon B., 1984, ARA&A, 22, 507

Margon B., Ford H. C., Katz J. I., Kwitter K. B., Ulrich R. K., Stone R. P. S., Klemola A., 1979a, ApJ, 230, L41

Margon B., Grandi S. A., Stone R. P. S., Ford H. C., 1979b, ApJ, 233, L63

Martí J., Paredes J. M., Estalella R., 1992, A&A, 258, 309

Mezger P. G., Henderson A. P., 1967, ApJ, 147, 471

Milgrom M., 1979, A&A, 76, L3

Miller-Jones J. C. A., Kapinska A. D., Blundell K. M., Stappers B., Braun R., 2007, in Tzioumis A. K., Lazio T. J., Fender R. P., eds, Proc. Sci., Bursts, Pulses and Flickering: wide-field monitoring of the dynamic radio sky. SISSA, Trieste, PoS(Dynamic2007)011

Mohan N., Rafferty D., 2015, Astrophysics Source Code Library, record ascl:1502.007

Moldowan A., Safi-Harb S., Fuchs Y., Dubner G., 2005, Adv. Space Res., 35, 1062

Morrison P., Roberts D., 1985, Nature, 313, 661

Nord M. E., Henning P. A., Rand R. J., Lazio T. J. W., Kassim N. E., 2006, AJ, 132, 242

Offringa A. R., de Bruyn A. G., Biehl M., Zaroubi S., Bernardi G., Pandey V. N., 2010, MNRAS, 405, 155

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

Omar A., Chengalur J. N., Anish Roshi D., 2002, A&A, 395, 227

Ochsenbein F., Bauer P., Marcout J., 2000, A&AS, 143, 23

Pal S., Chakrabarti S. K., Kraus A., Mandal S., 2006, Bull. Astron. Soc. India, 34, 1

Paladini R., Davies R. D., De Zotti G., 2004, MNRAS, 347, 237

Pandey M., Rao A. P., Ishwara-Chandra C. H., Durouchoux P., Manchanda R. K., 2007, A&A, 463, 567
Pratley L., McEwen J. D., d'Avezac M., Carrillo R. E., Onose A., Wiaux Y., 2018, MNRAS, 473, 1038

Price S. D., Egan M. P., Carey S. J., Mizuno D. R., Kuchar T. A., 2001, AJ, 121, 2819

Reich W., Fuerst E., Haslam C. G. T., Steffen P., Reif K., 1984, A&AS, 58, 197

Reich W., Fuerst E., Reich P., Reif K., 1990, A&AS, 85, 633

Remazeilles M., Dickinson C., Banday A. J., Bigot-Sazy M.-A., Ghosh T., 2015, MNRAS, 451, 4311

Rho J., Petre R., 1998, ApJ, 503, L167

Rickett B. J., 1986, ApJ, 307, 564

Roger R. S., Bridle A. H., Costain C. H., 1973, AJ, 78, 1030

Rudie G. C., Fesen R. A., Yamada T., 2008, MNRAS, 384, 1200

Sabin L. et al., 2013, MNRAS, 431, 279

Safi-Harb S., Ögelman H., 1997, ApJ, 483, 868

Sault R. J., Wieringa M. H., 1994, A&AS, 108, 585

Sault R. J., Teuben P. J., Wright M. C. H., 1995, in Shaw R. A., Payne H. E., Hayes J. J. E., eds, ASP Conf. Ser. Vol. 77, Astronomical Data Analysis Software and Systems IV. Astron. Soc. Pac., San Francisco, p. 433

Scaife A. M. M., Heald G. H., 2012, MNRAS, 423, L30

Schaudel D., Becker W., Voges W., Aschenbach B., Reich W., Weisskopf M., 2002, in Slane P. O., Gaensler B. M., eds, ASP Conf. Ser. Vol. 271, Neutron Stars in Supernova Remnants. Astron. Soc. Pac., San Francisco, p. 391

Scheers L. H. A., 2011, PhD thesis, Astronomical Institute Anton Pannekoeck, Univ. Amsterdam

Scheuer P. A. G., Williams P. J. S., 1968, ARA&A, 6, 321

Scholz P. et al., 2015a, ApJ, 800, 123

Scholz P. et al., 2015b, ApJ, 805, 85

Seaquist E. R., Gilmore W., Nelson G. J., Payten W. J., Slee O. B., 1980, ApJ, 241, L77

Seaquist E. R., Gilmore W. S., Johnston K. J., Grindlay J. E., 1982, ApJ, 260, 220

Shaver P. A., McGee R. X., Newton L. M., Danks A. C., Pottasch S. R., 1983, MNRAS, 204, 53

Slee O. B., 1995, Aust. J. Phys., 48, 143

Spencer R. E., 1979, Nature, 282, 483

Stephenson C. B., Sanduleak N., 1977, ApJS, 33, 459

Stil J. M. et al., 2006, AJ, 132, 1158

Su H. et al., 2017a, MNRAS, 465, 3163

Su H. et al., 2017b, MNRAS, 472, 828

Subrahmanyan R., Goss W. M., 1995, MNRAS, 275, 755

Swinbank J. D. et al., 2015, Astron. Comput., 11, 25

Tasse C., van der Tol S., van Zwieten J., van Diepen G., Bhatnagar S., 2013, A&A, 553, A105

Tingay S. J. et al., 2013, PASA, 30, 7

Trushkin S. A., Bursov N. N., Nizhelskij N. A., 2003, Bull. Spec. Astrophys. Obs., 56, 57

Trushkin S. A., Nizhelskij N. A., Tsybulev P. G., 2012, Astron. Telegram, 4484

Trushkin S. A., Nizhelskij N. A., Tsybulev P. G., 2014, Astron. Telegram, 6492

Trushkin S. A., Nizhelskij N. A., Tsybulev P. G., Zhekanis G. V., 2016, Astron. Telegram, 9481

van der Laan H., 1966, Nature, 211, 1131

van der Tol S., Jeffs B. D., van der Veen A.-J., 2007, IEEE Trans. Sig. Process., 55, 4497

van Gorkom J. H., Goss W. M., Seaquist E. R., Gilmore W. S., 1982, MNRAS, 198, 757

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

Vermeulen R. C., McAdam W. B., Trushkin S. A., Facondi S. R., Fiedler R. L., Hjellming R. M., Johnston K. J., Corbin J., 1993, A&A, 270, 189

Voges W. et al., 1999, A&A, 349, 389

Wenger M. et al., 2000, A&AS, 143, 9

Watson M. G., Willingale R., Grindlay J. E., Seward F. D., 1983, ApJ, 273, 688

Weiler K. W., Sramek R. A., Panagia N., van der Hulst J. M., Salvati M., 1986, ApJ, 301, 790

Westerhout G., 1958, Bull. Astron. Inst. Neth., 14, 215

Zheng H. et al., 2017, MNRAS, 465, 2901

This paper has been typeset from a TEX/LATEX file prepared by the author.

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

<footnote n="2">²Physics and Astronomy, University of Southampton, Highfield, Southampton SO17 1BJ, UK</footnote>

<footnote n="3">³ASTRON, the Netherlands Institute for Radio Astronomy, Postbus 2, NL-7990 AA Dwingeloo, the Netherlands</footnote>

<footnote n="4">⁴International Centre for Radio Astronomy Research, Curtin University, GPO Box U1987, Perth, WA 6845, Australia</footnote>

<footnote n="5">⁵Special Astrophysical Observatory RAS, Karachaevo-Cherkassian Republic, Nizhniy Arkhyz 369167, Russian Federation</footnote>

<footnote n="6">⁶Kazan Federal University, Kazan 420008, Russian Federation</footnote>

<footnote n="7">⁷Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, NL-1098 XH Amsterdam, the Netherlands</footnote>

<footnote n="8">⁸Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA</footnote>

<footnote n="9">⁹Department of Physics, The George Washington University, 725 21st Street NW, Washington, DC 20052, USA</footnote>

<footnote n="10">¹⁰University of Technology Sydney, 15 Broadway, Ultimo, NSW 2007, Australia</footnote>

<footnote n="11">¹¹ARC Centre of Excellence for All-sky Astrophysics (CAASTRO), The University of Sydney, NSW 2006, Australia</footnote>

<footnote n="12">¹²Jodrell Bank Centre for Astrophysics, School of Physics and Astronomy, The University of Manchester, Manchester M13 9PL, UK</footnote>

<footnote n="13">¹³Laboratoire AIM (CEA/IRFU – CNRS/INSU – Université Paris Diderot), CEA DSM/IRFU/SAp, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="14">¹⁴Station de Radioastronomie de Nançay, Observatoire de Paris, PSL Research University, CNRS, Univ. Orléans, F-18330 Nançay, France</footnote>

<footnote n="15">¹⁵Thüringer Landessternwarte, Sternwarte 5, D-07778 Tautenburg, Germany</footnote>

<footnote n="16">¹⁶Department of Astrophysics/IMAPP, Radboud University Nijmegen, PO Box 9010, NL-6500 GL Nijmegen, the Netherlands</footnote>

<footnote n="17">¹⁷LPC2E – Université d'Orléans / CNRS, F-45071 Orléans cedex 2, France</footnote>

<footnote n="18">¹⁸Astro Space Centre, Lebedev Physical Institute, Russian Academy of Sciences, Prosfoyuznaya Str. 84/32, Moscow 117997, Russian Federation</footnote>

<footnote n="19">¹⁹Department of Astronomy and Radio Astronomy Lab, University of California, Berkeley, CA 94720, USA</footnote>

<footnote n="20">²⁰Department of Physics and Electronics, Rhodes University, PO Box 94, Grahamstown 6140, South Africa</footnote>

<footnote n="21">²¹Department of Physics and Astronomy, University of the Western Cape, Private Bag X17, Bellville 7535, South Africa</footnote>

<footnote n="22">²²SKA South Africa, 3rd Floor, The Park, Park Road, Pinelands 7405, South Africa</footnote>

<footnote n="23">²³LESIA, Observatoire de Paris, CNRS, UPMC, Université Paris-Diderot, 5 place Jules Janssen, F-92195 Meudon, France</footnote>