© ESO 2020

1. Introduction

The large sizes (up to few megaparsecs) of extended extragalactic radio sources allow us to use them in tracing the history of galactic nuclear activity over hundreds of millions of years. Since their discovery, it has been revealed that most powerful radio jets have highly linear morphologies (e.g., Miley 1980). In classical models of radio jets, the orientation is associated with the spin axis of a supermassive black hole (SMBH) in the nucleus of the host galaxy. The alignment of kpc and Mpc-scale radio emission with pc-scale jets (e.g., Fomalont & Miley 1975) has demonstrated that the collimated jets hold a “memory” of their directions for more than 10⁸ years. Our understanding of the accretion processes by which the SMBHs are “fed” or the mechanisms that determine the orientation of their spin axes is still incomplete.

An intriguing question concerns whether there could be some connection between the orientations of the SMBH spin axes and properties of the cosmic filaments in which the radio sources and their host galaxies are found. The possibility of such a connection has been suggested in recent evidence for non-uniformity in radio-source position angles over large regions of the sky found by Taylor & Jagannathan (2016) and Contigiani et al. (2017).

If the radio sources are indeed aligned with respect to the large-scale structure in which they are found, a possible cause could be attributed to angular momentum transfer during the early stages of galaxy formation. The tidal torques imparted on the collapsing halos is found to influence the spin and shape of galaxies in N-body simulations (e.g., White 1984; Codis et al. 2012, 2018; Laigle et al. 2015; Kraljic et al. 2020). However, the angular momentum vector of the active galactic nucleus (AGN) and the host galaxy are found to be misaligned and generally uncorrelated (Hopkins et al. 2012), indicating that this explanation is incorrect or incomplete.

Furthermore, there is substantial evidence that AGNs are associated with mergers, based on both observations and simulations (e.g., Chiaberge et al. 2015; Croton et al. 2006). If these mergers occur preferentially along the filaments of the large-scale structure, these could orient the central SMBHs in a particular way, resulting in a preferential alignment of the extended radio sources. Hence, if the alignment of radio sources on large scales is confirmed, this would have significant implications for models of the formation of galaxies and active galactic nuclei.

Additional evidence that there may be a connection between the orientation of the spin axes of SMBHs that power active galactic nuclei and the cosmic filaments in which they lie comes from observations of large-scale statistical alignments in the optical polarization position angles of quasars (e.g., Hutsemekers 1998; Hutsemékers & Lamy 2001; Jain et al. 2004). Evidence has also been found for the polarization angle of quasars to be either parallel or perpendicular to the large-scale structures they inhabit (e.g., Hutsemékers et al. 2014; Pelgrims & Hutsemékers 2016).

A more extensive investigation of the large-scale distribution of radio source orientations is warranted. Surveys with the Low-Frequency Array (LOFAR) High Band Array (HBA; van Haarlem et al. 2013) are especially suited for carrying out such studies, because they (i) are conducted at sufficiently low frequencies to detect steep-spectrum extended synchrotron radio structures, (ii) have sufficient angular resolution, with a ∼6″ half-power beam width (HPBW), to resolve 50 (100) kpc-sized radio sources out to redshift ∼1 (> 6), and (iii) have the sensitivity and dynamic range needed to detect and measure orientations for an unprecedented number of sources.

Here, we describe such an investigation using position angles of radio sources from the LOFAR Two-metre Sky Survey Data Release I (LoTSS-DR1; Shimwell et al. 2019). We first describe the data in Sect. 2. The criteria we used to select sources with well-defined position angles from the 318 520 radio sources from the survey are discussed in Sect. 3. The statistical methods we used to explore non-uniformity in the source alignments are explained in Sect. 4. Our results are given in Sect. 5, where we report evidence for non-uniformity in the source alignments. Finally, in Sects. 6 and 7, we discuss the robustness and implications of the results.

Throughout this paper, we adopt the Planck 15 cosmology (Planck Collaboration XIII 2016). This cosmology is defined by the following relevant parameters: H₀ = 67.8 kms⁻¹Mpc⁻¹, Ω_m = 0.308, Ω_Λ = 0.692.

2. The data

Our sample is taken from the LoTSS, a sensitive low-frequency (120–168 MHz) survey that will ultimately cover the entire northern sky. The first data release comprises 2% of the whole survey (424 square degrees) in the HETDEX Spring Field region (right ascension 10h45m to 15h30m and declination 45° to 57°; Shimwell et al. 2019). It contains more than 300 000 radio sources that have a signal to noise ratio (S/N) of > 5. The images have a HPBW resolution of ∼6″, a median sensitivity of 71 μJy beam⁻¹, and a positional accuracy better than ∼0.2″.

The data used were taken from the “value-added” radio + optical catalog of Williams et al. (2019) of 318,520 LoTSS sources, which includes, where possible, identifications and redshifts of the optical counterparts. The optical identifications were made using either a likelihood ratio method or by human visual classification through the LOFAR Galaxy Zoo¹. Spectroscopic redshifts in the added-value catalog were taken, where available, from the Sloan Digital Sky Survey Data Release 14 (Abolfathi et al. 2018). Otherwise photometric redshifts were estimated using a hybrid methodology based on traditional template fitting and machine learning (see Duncan et al. 2019).

3. Source selection

For the alignment uniformity analysis, our goal was to select double-lobed radio sources with clearly defined position angles from the LoTSS value-added catalog. To identify such sources, we used the following method:

First, we filter the catalog to contain only high S/N extended sources. We define sources as extended if they have a major axis that is larger than five times the restoring beam size. The adopted selection criteria are:

$$\begin{array}{c}\hfill S_{\mathrm{peak}}/N>10\mathrm{and}a>{30}^{″},\end{array}$$

where S_peak is the peak flux density of the LOFAR source at 144 MHz and a is the size of the major axis of the source. The major and minor axes of some sources are not directly provided for sources that have been processed by the LOFAR Galaxy Zoo (LGZ). Instead, an equivalent “LGZ_Size” and “LGZ_Width” parameter is provided. The construction of the source dimensions from the LGZ data is described in Williams et al. (2019). Throughout this paper, we set the major and minor axes of the sources processed by LGZ as the “LGZ_Size” and “LGZ_Width”, respectively. Additionally, uncertainties for the LGZ shape parameters (source size, width, and position angle) are not provided by the value-added catalog. We discuss any effects due to uncertainties in the position angles in Sect. 6.

Next, we keep only the sources with a double lobed structure. We enforce this criterion by imposing the condition that sources must be fitted by multiple Gaussian components by the initial source finder PyBDSF (Mohan & Rafferty 2015). This is indicated by the “S_Code” of the source in the catalog. It is also possible that the source is a bright resolved nearby galaxy and these are identified with the “ID_flag” code where the first digit is 2. We remove these sources as well, using:

$$\begin{array}{c}\hfill \text{``}\mathrm{S}\_\mathrm{Code}=M"\mathrm{and}\text{``}\mathrm{ID}\_\mathrm{flag}"\ne 2.\end{array}$$

Imposing these criteria results in a reduction of the sample of 318 520 sources to a sample of 7688 bright extended linear sources. We check the catalog for sources that might have been identified multiple times by examining the distance from every source to its nearest neighbor. We investigate all sources that have a nearest neighbor within ten synthesized beams (1 arcmin). If the source has a different optical identification from its nearest neighbor, we can be reasonably sure that it is not a duplicate entry. When the source has an optical identification while the nearest neighbor does not, or both the source and the nearest neighbor lack an optical identification, we cannot be certain that these entries are not duplicates. To err on the side of caution, we remove all sources from our sample that have a nearest neighbor within ten synthesized beams, unless they have a different optical identification from their nearest neighbor. We find that 165 sources have a nearest neighbor within ten synthesized beams, and 32 of these have a different optical identification from their nearest neighbor. We expect that removing the other 133 entries would not impact the strength of a possible alignment effect since radio source alignments have been claimed on scales of at least a degree (Taylor & Jagannathan 2016; Contigiani et al. 2017) and these source separations are on a smaller angular scale than this. Thus, the final sample contains 7555 selected sources.

4. Statistical methods

To determine the departure from uniformity of the alignment of radio sources on the sky, an appropriate statistical method must be used that accounts for effects due to the geometry of the celestial sphere. We shall do this by introducing the concepts of “parallel transport” and “dispersion measure”.

4.1. Parallel transport

The position angle in the LoTSS catalog is defined as the angle of the major axis of a source measured east of the local m (north) direction. To have a consistent definition of the position angle across all pointings, we translated the position angles to be measured east of the direction of the north celestial pole.

Because the position angle is defined with respect to the local meridian, the vectors corresponding to the position angles on different points of the celestial sphere cannot be compared directly. These vectors must be transported along the great circle joining these points. Following Jain et al. (2004) and Contigiani et al. (2017), we use the parallel transport method, by which the radio source “vectors” can be transported to a different position on the celestial sphere. This method is described below for completeness.

We parametrize the celestial sphere with local unit vectors (u_r, u_δ, u_α) which point, respectively, to the center of the sphere, north along the local meridian and eastwards on the sphere. We wish to compare the position angles, θ_p1 and θ_p2, of sources 1 and 2 with positions, P₁ and P₂, on the celestial sphere (Fig. 1). The vector resulting from the position angle, θ_p1, of source 1 at location, P₁, is given, in terms of the local basis, by

$$\begin{array}{c}\hfill {\mathbf{v}}_1=cos{\theta }_{p_1}{\mathbf{u}}_{{\delta }_1}+sin{\theta }_{p_1}{\mathbf{u}}_{{\alpha }_1}.\end{array}$$(1)

Fig. 1. Illustration of parallel transport. Vector v1 corresponding to a position angle θp1 and vector v2 corresponding to position angle θp2 are shown. In order to compare v1 to v2, v1 must be parallel-transported along the great circle indicated by the curve from location P1 to location P2. The transported vector is indicated by $v_1^{\prime }$ and the local basis vectors are denoted by (uδ, uα). In parallel transport, the angle α between the vector tangent to the sphere ut and the vector v remains fixed. Figure adapted from Jain et al. (2004).

To define a coordinate-invariant inner product we parallel transport the vector, v₁, to the position, P₂, to obtain vector, ${\mathbf{v}}_{\mathbf{1}}{}^{\prime }$. Vector ${\mathbf{v}}_{\mathbf{1}}{}^{\prime }$ then makes an angle, ${\theta }_p^{\prime }$, with respect to the local north-pointing vector, u_δ2. To find the transported angle, ${\theta }_p^{\prime }$, let u_s be the unit vector perpendicular to the plane containing the two radial vectors u_r1 and u_r2. Thus, u_s is found by

$$\begin{array}{c}\hfill {\mathbf{u}}_s=\frac{{\mathbf{u}}_{r_1}\times {\mathbf{u}}_{r_2}}{|{\mathbf{u}}_{r_1}\times {\mathbf{u}}_{r_2}|}.\end{array}$$(2)

Consider now the unit vectors, u_t1 and u_t2, at points, P₁ and P₂, and tangent to the great circle connecting P₁ and P₂. These vectors are given by:

$$\begin{array}{c}\hfill {\mathbf{u}}_{t_{1,2}}={\mathbf{u}}_s\times {\mathbf{u}}_{r_{1,2}}.\end{array}$$(3)

In terms of the local basis, these vectors can be written as:

$$\begin{array}{c}\hfill {\mathbf{u}}_{t_1}={\mathbf{u}}_{{\delta }_1}·{\mathbf{u}}_{t_1}{\mathbf{u}}_{{\delta }_1}+{\mathbf{u}}_{{\alpha }_1}·{\mathbf{u}}_{t_1}{\mathbf{u}}_{{\alpha }_1},\end{array}$$(4)

where

$$\begin{array}{c}\hfill {\mathbf{u}}_{{\delta }_1}·{\mathbf{u}}_{t_1}=\frac{-sin{\delta }_{1}cos{\delta }_2+cos{\delta }_{1}sin{\delta }_{2}cos({\theta }_{p_1}-{\theta }_{p_2})}{\sqrt[]{1-{({\mathbf{u}}_{r_1}·{\mathbf{u}}_{r_2})}^2}},\end{array}$$(5)

$$\begin{array}{c}\hfill {\mathbf{u}}_{{\alpha }_1}·{\mathbf{u}}_{t_1}=\frac{sin{\delta }_{2}sin({\alpha }_2-{\alpha }_1)}{\sqrt[]{1-{({\mathbf{u}}_{r_1}·{\mathbf{u}}_{r_2})}^2}},\end{array}$$(6)

$$\begin{array}{c}\hfill {\mathbf{u}}_{{\delta }_2}·{\mathbf{u}}_{t_2}=\frac{-sin{\delta }_{2}cos{\delta }_1+cos{\delta }_{2}sin{\delta }_{1}cos({\theta }_{p_1}-{\theta }_{p_2})}{\sqrt[]{1-{({\mathbf{u}}_{r_1}·{\mathbf{u}}_{r_2})}^2}},\end{array}$$(7)

$$\begin{array}{c}\hfill {\mathbf{u}}_{\alpha ,2}·{\mathbf{u}}_{t_2}=\frac{-sin{\delta }_{1}sin({\alpha }_1-{\alpha }_2)}{\sqrt[]{1-{({\mathbf{u}}_{r_1}·{\mathbf{u}}_{r_2})}^2}}.\end{array}$$(8)

As v₁ is parallel-transported along the great circle to position, P₂, with its angle with respect to the tangent of the great circle remaining fixed. Thus, to determine the angle by which the vector has turned due to this transport, we consider the orientation of u_t1 and u_t2 with respect to the local basis at the two points where the sources lie. We call ξ₁ the angle between u_t1 and u_α1, and ξ₂ the angle between u_t2 and u_α2. These angles are given, per definition of the inner product, by

$$\begin{array}{c}\hfill {\xi }_{1,2}=arccos({\mathbf{u}}_{{\alpha }_{1,2}}·{\mathbf{u}}_{t_{1,2}}).\end{array}$$(9)

Thus, the transported ${\mathbf{v}}_{\mathbf{1}}{}^{\prime }$ makes an angle ${\theta }_{p,1}^{\prime }={\theta }_{p,1}+({\xi }_2-{\xi }_1)$ defined with respect to the local coordinates in P₂. Hence we can now define the generalized dot product between v₁ and v₂ as the dot product between the transported vector ${\mathbf{v}}_{\mathbf{1}}{}^{\prime }$ and v₂:

$$\begin{array}{c}\hfill {\mathbf{v}}_1\odot {\mathbf{v}}_2={\mathbf{v}}_{\mathbf{1}}^{\prime }·{\mathbf{v}}_2=cos({\theta }_{p_1}-{\theta }_{p_2}+{\xi }_2-{\xi }_1).\end{array}$$(10)

Equation (10) can generally be used in any problem that considers angles on a sphere. In particular, when comparing the difference between position angles, it makes sense to redefine the generalized inner product between two position angles as

$$\begin{array}{c}\hfill ({\theta }_{p_1},{\theta }_{p_2})=cos[2({\theta }_{p_1}-{\theta }_{p_2}+{\xi }_2-{\xi }_1)],\end{array}$$(11)

where, since the position angles range from 0 to π, it assumes values of ∈ ( − 1, 1) and where +1 expresses the perfect alignment between θ_p1 and θ_p2 and −1 indicates perpendicular orientations.

4.2. Statistical test

To test the significance of a possible alignment in source position angles, we use the dispersion measure (Jain et al. 2004; Contigiani et al. 2017). We briefly repeat the definition of the dispersion measure here for completeness.

The dispersion depends only on the differences between neighboring position angles and it is, therefore, a suitable choice when testing for alignment on different scales. The dispersion measure of source, i, as a function of a position angle, θ, is defined as

$$\begin{array}{c}\hfill d_{i,n}(\theta )=\frac{1}{n}{\displaystyle \sum _{k=1}^n(\theta ,{\theta }_k),}\end{array}$$(12)

where n is the number of nearest neighbors that are considered around source i, including the source itself, and θ_k is the position angle of the respective neighbors. The generalized inner product (θ, θ_k) is defined by Eq. (11).

The position angle θ that maximizes the dispersion around source i is analogous to the definition of the mean position angle of source i and its n nearest neighbor s. The magnitude of d_i, n|_max is, then, a measure of the dispersion around this mean. The dispersion can take a maximum value of 1, which corresponds to perfect alignment of all n nearest neighbors. To find the value of θ that maximizes the dispersion, we take the derivative of Eq. (12) with respect to θ and, after some intermediate steps, we arrive at the following expression for d_i, n|_max:

$$\begin{array}{c}\hfill d_{i,n}{|}_{max}=\frac{1}{n}{{\displaystyle [{(\sum _{k=1}^{n}cos{\theta }_k)}^2+{(\sum _{k=1}^{n}sin{\theta }_k)}^2]}}^{1/2}.\end{array}$$(13)

The statistic, so that we may test for the non-uniformity of alignment in a sample of N sources, is then defined as:

$$\begin{array}{c}\hfill S_n=\frac{1}{N}{\displaystyle \sum _{i=1}^N{d}_{i,n}{|}_{max},}\end{array}$$(14)

which is simply the average of the maximum dispersion for a number of nearest neighbors, n, calculated over all N sources in the sample. This statistic thus measures the strength of a local alignment signal in the full sample of N sources while considering the n nearest neighbors of every source.

The significance level for rejecting the null hypothesis that a sample of sources is randomly oriented is then given by comparing the statistic of the dataset, S_n, to the distribution of the statistic for simulated samples that are randomly oriented. It is found through a one-tailed significance test, expressed as:

$$\begin{array}{c}\hfill \mathrm{SL}=1-\mathrm{\Phi }\left(\frac{S_n-\langle S_n{|}_{\mathrm{MC}}\rangle }{{\sigma }_n}\right),\end{array}$$(15)

where Φ is the cumulative standard normal distribution function. Here, < S_n|_MC> and σ_n are, respectively, the expectation value and standard deviation of S_n in the absence of alignment. These values can be found through Monte Carlo simulations of randomly oriented sources.

Jain et al. (2004) verified that for randomly oriented samples of sources, S_n is normally distributed if N ≫ n ≫ 1 is satisfied. With the dispersion measure and the resulting statistic, the significance level at which the hypothesis of uniformity in the position angles should be rejected can be calculated on a local scale by probing different numbers of nearest neighbors. Since the number of nearest neighbors can be translated to fixing apertures with angular radii extending to the nth nearest neighbor around all sources, S_n can be used to probe the significance of alignment on different angular scales. We note that different S_n are not independent since the dispersion is an average of n neighbors. This statistic thus probes alignment up to scales corresponding to n and once a signal is detected for some n, a preferentially positive signal is expected for larger n.

If the redshifts of the sources are known, this method can be extended to probing nearest neighbors in 3D space. In this way, the dependence of a possible alignment effect and S_n as a function of physical scale can be probed.

5. Results

We first tested the uniformity of the LoTSS radio source position angles over the complete 424 square degrees of the available survey to give an indication of possible systematic effects. The distribution of position angles is given in Fig. 2. We expect the position angles to be uniformly distributed over this relatively large patch of the sky, if no systematic effects are present. From Fig. 2, we can see that no major systematic effects are present, although the distribution is not quite uniform. To check if the distribution is consistent with a uniform distribution of sources, we applied the Kolmogorov-Smirnov (K–S) test (e.g., Frommert et al. 2012). The K–S test resulted in a p-value of 0.030 per cent. This is strong evidence for rejecting the null hypothesis that the distribution of position angles over the complete sample is uniform, which indicates some systematic (survey-wide) bias in our sample. Still, the local alignment signal might be stronger or weaker depending on the nature of the effect that is causing the alignment.

Fig. 2. Position angle distribution of the complete sample of selected sources.

5.1. Two-dimensional analysis

To determine whether the hypothesis of uniformity in position angles on different angular scales should be rejected and if so, at what significance level, we compared results for the observed LoTSS sample with those for 1000 simulated randomly distributed position angle samples. These samples were generated by randomly shuffling the position angles among the sources to maintain the same global position angle distribution and source positions.

The sample was checked for local alignment by probing the statistic, S_n, for different numbers of nearest neighbors. To express the statistic in terms of angular scale, a circular aperture with a radius extending to the nth neighbor of every source is drawn. We translated the number of nearest neighbors to an approximate corresponding angular scale by taking the median angular radius of all these apertures. This dependency is shown in Fig. 3.

Fig. 3. Median angular radius corresponding to drawing a circular aperture around every source with an angular radius bound by the nth neighbor around that source.

The significance level at which the null hypothesis should be rejected (of the position angles being uniformly distributed) is given as a function of the number of nearest neighbors (or corresponding angular scale) in Fig. 4. There is strong evidence that the hypothesis of uniformity in radio source position angles should be rejected on angular scales of about four degrees, with a significance level of < 10⁻⁵.

Fig. 4. Logarithm of the significance level for which uniformity in position angles as a function of the number of nearest neighbors n should be rejected for the sample of 7555 selected sources. The conversion to angular scale is shown in Fig. 3.

To investigate the effect further, we split our sample into four equal frequency flux density bins to have the maximum number of sources in every bin, as given in Table 1. For each bin, this table includes the median flux density, the median redshift, the median source angular size, and the maximum significance level at which the null-hypothesis of position angle uniformity should be rejected, taken from Fig. 5.

Fig. 5. Logarithm of the significance level for the sample of sources split into four equal frequency total flux density bins, as defined in Table 1.

Figure 5 shows the significance level of position angle alignment for the four flux density bins as a function of angular distance. Interestingly, the highest flux density bin shows very strong evidence for alignment, up to scales of roughly ten degrees, but most significantly around four degrees, while all other bins are consistent with uniform distributions. This shows that the effect seen in the total sample is caused by the highest flux density sources only.

5.2. Three-dimensional analysis

We carried out an analysis of alignment uniformity using 3D source positions, after removing all sources from our sample that do not have a spectroscopic or photometric redshift tabulated in the value-added catalog (Williams et al. 2019; Duncan et al. 2019). This reduced the size of our sample to 4212 sources. The number of photometric and spectroscopic redshifts is 2311 and 1901, respectively. We emphasize that the statistical method is exactly the same for this analysis. The only difference between the 2D and 3D analysis is that 3D source positions are now used to find the n nearest neighbors for every source. The distribution of position angles of these 4212 sources is shown in Fig. 6. The K–S test indicates a p-value of 1.0 per cent, indicating that for this sample of 4212 sources, there is weak evidence for rejecting the null hypothesis for the uniformity of position angles.

Fig. 6. Distribution of position angles for the sample of 4212 sources that have a redshift measurement.

In the analysis for local alignment, the samples were compared again with 1000 simulated uniformly distributed position angle samples, generated by randomly shuffling the position angles among the sources. To repeat the analysis in three dimensions, each source was assigned a position in 3D space according to their right ascension, α, declination, δ, and comoving distance, r, as follows:

$$\begin{array}{c}\hfill \begin{array}{c}x=rcos\alpha cos\delta ,\hfill \\ y=rsin\alpha cos\delta ,\hfill \\ z=rsin\delta .\hfill \end{array}\end{array}$$(16)

The nearest neighbors were then computed in 3D space according to these positions to probe for alignment on local scales.

Figure 7 shows the significance level at which the hypothesis of uniformity in position angles can be rejected for the 4212 sources that have a redshift, both in a 3D and a 2D analysis. The 3D analysis does not show strong evidence for an alignment effect.

Fig. 7. Logarithm of the significance level at which position angle uniformity should be rejected, as a function of the number of nearest neighbors n for the 4212 sources in that have redshifts available. The dashed line indicates the results of the 2D analysis and the solid line the results of the 3D analysis.

Since the sources with the largest flux densities are the main contributor to the alignment effect in the 2D analysis, we also calculated the significance for the highest flux density sources in the 3D analysis. We split the 4212 sources into four equal frequency total flux density bins, which defines the highest flux density bin as all sources with a total flux density > 108 mJy. This makes the flux cut for the highest flux density bin slightly higher than the equivalent in the 2D analysis, but we decide to use this flux cut to have a fairer comparison between the different flux density bins within the 3D analysis.

The significance level at which position angle uniformity can be rejected for the highest flux density bin in 3D is shown in Fig. 8. This figure shows, interestingly, that the 2D analysis of these 1051 sources still shows strong evidence for alignment up to scales of four degrees. However, this signal is not present in the 3D analysis. No signal was found in the other flux density bins, either in 2D or in 3D.

Fig. 8. Logarithm of the significance level at which position angle uniformity should be rejected as a function of the number of nearest neighbors n for the 1051 sources with total flux density > 108 mJy and a redshift measurement. The dashed line indicates the results of the 2D analysis and the solid line the results of the 3D analysis.

The difference between the 2D and 3D analysis indicates that the 2D alignment effect is due to some unknown systematic effect, since a physical effect would invariably cause stronger alignment in the 3D analysis than in the 2D analysis. Additionally, we inspected whether the most radio luminous sources are also the most aligned sources, which would be expected from the similar median redshift per flux density bin. However, no alignment signal was found in either the 2D or 3D analysis of the 1,051 highest radio power sources.

Although it reduces the sub-sample sizes even further, we also tested if the results depend on whether the redshifts were photometric or spectroscopic. Figure 9 shows the results for the 523 sources that have a spectroscopic redshift. The figure shows that in both the 2D analysis and 3D analysis of these subsets no significant signal is present. This is not surprising given the small number of sources in the spectroscopic subsample.

Fig. 9. Logarithm of the significance level at which position angle uniformity should be rejected as a function of the number of nearest neighbors n for the 523 sources with total flux density > 108 mJy and a spectroscopic redshift measurement. The dashed line indicates the results of the 2D analysis and the solid line the results of the 3D analysis.

6. Discussion

6.1. Robustness of the results

The robustness of the results depends on the uncertainties of the position angles that were fit to the sources. The position angles of sources in the LOFAR value-added source catalog do not all include uncertainties. The subset of sources that was classified in the LGZ project do not come with position angle uncertainties. To examine the position angle uncertainties, we are thus restricted to using the sources that are classified by the source finder PyBDSF only. To examine the position angle uncertainties, we plot the 1σ uncertainties as given by the catalog for these sources. These are shown in Fig. 10. The figure shows that 81% of the sources have position angle uncertainties smaller than ten degrees. Thus, we can approximate the uncertainty in the final significance level by assuming every source in our sample has a 1σ uncertainty of 10. Since most sources have smaller uncertainties, this assumption is likely to overestimate the uncertainty in the fitted position angles and, thus, in the final significance level.

Fig. 10. 1σ uncertainties on the position angles of sources in our sample that are not classified by the LOFAR galaxy zoo, but by the source finder only.

To approximate the error in the final significance level as a function of a 1σ error of ten degrees in the position angle, we must propagate this error through the statistical analysis of Sect. 4. However, there is no straightforward procedure to define the general error on the extracted significance level as a function of the error on the measured position angles. Simple error propagation can be applied to the calculation of S_n, but it becomes complicated when a one-tailed significance level is extracted. This is due to the dependence of the significance level on the position of S_n in the distribution of S_n|_MC (Eq. (15)). If S_n lies far from the mean of the normal distribution, a given change in S_n will lead to a smaller change in significance level than when S_n lies near the mean of the distribution of S_n|_MC. This is a direct effect of the cumulative normal distribution function being steepest near the mean and flattest near the edges. Moreover, considering that for every sample, S_n|_MC is found by simulating 1000 random datasets by randomly shuffling the position angles of the sources, the distribution of S_n|_MC will be unique for every sample that we have considered. Therefore, we can only approximate the error on concrete results and cannot give a general 1σ confidence level that will apply for a range of samples.

The initial sample of sources rejected uniformity at a significance level of < 10⁻⁵ (Fig. 4). The signal was found with a number of nearest neighbors between 467 and 916, corresponding to an angular scale between 3.1 and 4.6 degrees. Figure 11 shows the distribution of the simulated data and the highly significant value of S_n for these two bounds. We calculate the error on S₄₆₇ and S₉₁₆ and translate these errors to bounds on the significance values.

Fig. 11. Distribution of the 1000 simulated values Sn|MC and the highly significant values Sn for n = 467 and n = 916. Plotted for the initial sample of sources.

Assigning for each position angle in our sample a 1σ error of 10 degrees and applying standard error propagation, we find for the resulting values of S₄₆₇ and S₉₁₆, 0.070 ± 0.0025 and 0.065 ± 0.0026, respectively. Taking the 1σ lower and upper bound of S₄₆₇ and calculating the significance level of these two bounds results in the lower and upper bound logarithmic significance levels of −3.37 and −7.03. For S₉₁₆, the same method leads to lower and upper bound logarithmic significance levels of −3.24 and −7.16. Thus, strong evidence to reject uniformity in this sample at scales between 3.1 and 4.5 degrees is still found after applying possible uncertainties in the position angles. We can conclude that assuming a 1σ error of ten degrees on the position angle of all sources, the effect of an uncertainty in the position angles is quite powerful, but the significance level does remain strong enough to reject uniformity. Thus, the resulting significance level of log SL = −5.0 to reject uniformity found in Fig. 4 for angular scales between 3.1 and 4.6 degrees should be stated with the approximate bound of log SL = −5.0 ± 2.0.

As stated previously in this paper, the difference in significance level for the same variation in S_n is dependent on the position of S_n, thus, it is also dependent on the significance level itself. Therefore, we reiterate that the change of two orders of magnitude in significance, found for the subset considered in this section, should not be applied to different subsets. We can apply the same calculation to the 3D analysis of the initial sample (Fig. 7), where no result was found. We chose to investigate n = 500, which corresponds to a significance level of 10^−1.6. This results in $log\mathrm{SL}=-1.7_{-0.97}^{+0.71}$; still without changing the signal to strong (< 10⁻³) evidence for alignment. Repeating the same calculations for the 2D analysis of the highest flux density sources that have an available redshift (Fig. 8) results in the approximate bounds $log\mathrm{SL}=-5.3_{-1.6}^{+1.4}$ for n = 100.

6.2. Interpretation of the results

Our complete sample of 7555 double sources with a well-defined orientation was found to be inconsistent with a uniform distribution with a K–S test significance of 0.030 percent, which already indicates a global systematic effect in the data. However, the analysis of local alignment depends on the contrast between the statistic, S_n, found for our dataset and the statistic, S_n|_MC, found in absence of alignment. The statistic in absence of alignment was generated by randomly shuffling the position angles amongst the sources to maintain the same geometry and global position angle distribution. The advantage of this method over generating position angles from the uniform distribution 𝒰[0, 180) is that it diminishes the effect of a possible global systematic present in our data sample. This is due to a global systematic then also being included in the distribution of the statistic S_n|_MC. Therefore, as long as n ≪ N, where N is the number of sources in the sample that is examined, the effect of the deviation from uniformity of the whole sample will not have considerably impacted the result of the significance of local alignment.

To identify which particular sources are causing the observed signal, we examine which sources show the strongest alignment effect in 2D space. For this, we use the calculated maximum dispersion measure d_i, n|_max (Eq. (13)), which measures the significance of the alignment of a source, i, and its nearest neighbors, n. We plot the maximum dispersion for every source in the initial sample of 7555 sources as a function of right ascension and declination for n = 700 in Fig. 13. From this figure, it becomes apparent that there is not a single region where the alignment is most pronounced, but rather, that there is an alternation between strongly aligned and less strongly aligned regions. This contradicts the observed effect being attributed to a survey-wide systematic effect, as then all sources would have similar maximum dispersion, regardless of their position. Additionally, the scale of the alternation between aligned and non-aligned regions is larger than the typical separation between LOFAR pointings (2.58 degrees; Shimwell et al. 2019), which makes the origin of the systematic effect even more elusive.

We also found that the alignment signal was most significant for sources with the largest flux densities, as indicated by Figs. 5 and 8. However, an analysis of the sources with the highest radio power did not show an alignment effect, either in 2D or in 3D. Thus, it seems that only apparent source properties, rather than physical source properties, are correlated with the alignment effect, which could point towards an intrinsic effect of the survey, although radio power and source brightness are not strongly correlated for radio sources. Most importantly, the fact that the alignment effect is not present when using the 3D positions of the high flux density sources to find the nearest neighbors but is present when using 2D source positions (Fig. 8) may indicate a systematic error in the survey images or overall catalog, which is most noticeable or perhaps only present for the highest flux density sources. However, interpreting this result is not straightforward due to the relatively large uncertainties in the third (redshift) dimension.

To further examine the impact of redshift uncertainties, we investigated whether 2D source positions are a better indicator of physical proximity than 3D source positions given different uncertainties in the photo-z estimates. This was done for the sample in Fig. 8 with the 1051 highest flux density sources that showed a signal in 2D around n = 101 and no signal in 3D. We assumed, for this simulation, that the “true” source positions are given by the spectroscopic redshifts and best available photo-z estimates (i.e., that the photo-z scatters around the true redshift). The goal is to investigate what fraction of nearest neighbors that are found by using 3D positions agrees with the nearest neighbors found using the “true” source positions.

The redshift of sources with a photo-z estimate was perturbed by a Gaussian with standard deviations of the usual form σ(1 + z) and spectroscopic redshifts were left intact. The n = 101 nearest neighbors in 3D were then found for every source given the perturbed redshifts and the fraction of “correct” nearest neighbors was calculated. What we mean by “correct” here is that a nearest neighbor that was found is also one of the n = 100 nearest neighbors using the “true” source positions, thus, we do not take the ordering into account (as the statistical method for a single value of n does not do either). Figure 12 shows the result of re-sampling the photometric redshifts 200 times and computing the fraction of “correct” nearest neighbors. The figure shows that using 2D coordinates leads to finding 23% of the true physically close sources, while using 3D coordinates leads to finding more than 35% of the physically close sources, even with standard deviations as as large as 0.1 × (1 + z). Our assumed scatter of 0.1 × (1 + z) represents a conservative upper limit on the expected precision of the LoTSS photometric redshift estimates, with the typical scatter for the radio population found to range from 0.03 × (1 + z) for radio sources dominated by stellar emission and 0.08 to 0.1 × (1 + z) for the more difficult quasar and AGN population (see Duncan et al. 2019).

Fig. 12. Percentage of the n = 101 nearest neighbors of every source in the photo-z perturbed sample of 1051 brightest sources that agrees with the nearest neighbors found in the unperturbed sample. See text for more details.

Fig. 13. Scatter plot of the maximum dispersion measure (Eq. (13)) for every source, which indicates the strength of an alignment signal, of the selected sample of 7555 radio sources plotted for n = 700 as a function of right ascension and declination.

Thus, for this sample of sources and n = 101, it would be likely to find a stronger alignment using 3D coordinates if the alignment effect is correlated with physical source positions. However, we are finding stronger alignment using 2D coordinates, which qualitatively implies that the alignment effect is more correlated with observed 2D source positions than it is with 3D source positions.

6.3. Scale of the alignment

The angular scale of the observed alignment effect is substantially larger than that of the two previous radio structure studies. Taylor & Jagannathan (2016) investigated an area of 1.2 square degrees, and were thus limited to finding alignment within this area. Therefore, the angular scale of one degree found in that study might be underestimated and may still be in agreement with the results of this study. Contigiani et al. (2017), however, did not suffer this limitation, as they studied an area of 7000 square degrees and found an effect up to scales smaller than 2.5 degrees, with the maximum alignment signal at 1.5 degrees, while the distribution of source redshifts is not significantly different from that in this study. While the scale of the maximum effect does not agree with the angular scale of larger than three degrees found in this study, Contigiani et al. (2017) limited their search to angular scales below 2.5 degrees, so the signal may perhaps be present on larger scales in the FIRST survey as well. Further research into radio jet alignment at larger angular scales is thus needed.

Should the effect turn out to be physical, it is useful to compute the approximate physical scale corresponding to the effect that is observed. We computed the physical scale corresponding to the angular scale at which the alignment was found in this study by assuming the median redshift of the sample of sources for which a redshift is available (z = 0.56). Converting the angular scale of four degrees to comoving distances yields a corresponding physical scale of 103 h⁻¹ Mpc. Although, as expected due to the limits in angular scales of the previous studies, the physical scale of 100 h⁻¹ Mpc found in this study does not agree with the physical scale of the two previous studies of radio lobe alignment discussed earlier, it is in agreement with physical scales where other studies have found AGN alignment effects. As stated in Sect. 1, several studies have found that the radio polarization of quasars is preferentially aligned either perpendicular or parallel to the major axis of the surrounding large-scale large quasar groups (LQGs). These effects range from distances of the order of 150 Mpc (Tiwari & Jain 2013) to distances larger than 300 h⁻¹ Mpc (Pelgrims & Hutsemékers 2016). The physical scales found in this study agree with the physical scale for alignment with large-scale structures and coincides with the observed first peak of the Baryon Acoustic Oscillations (BAO; Eisenstein et al. 2005), while still abiding by the upper limits of homogeneity of the Universe, found to be on the order of 260 h⁻¹ Mpc (e.g., Yadav et al. 2010).

7. Conclusion

In this study, we analyze the uniformity in the position angles of extended radio sources with well-defined linear double structures from the initial installment of the LOFAR Two-metre Sky Survey (LoTSS). The combination of low frequencies (with sensitivity to extended structures) and the relatively high angular resolution of LOFAR makes it an excellent survey in the search for systematic alignments in the position angles of radio sources.

We extracted 7555 LoTSS-extended sources with well-defined position angles from the 318 520 sources in the radio/optical value-added catalog of LOFAR sources in the HETDEX Spring Field region. To test for the alignment of position angles in this sample, the spherical nature of position angles and the effect of transporting these angles over the celestial sphere were taken into account using statistical methods originally developed to test for the alignment of polarization vectors. We find evidence for alignment in our initial sample of sources. The null hypothesis that the position angles are distributed uniformly can be rejected with a significance level of < 10⁻⁵ for an angular scale of four degrees, with the most non-uniformity present for radio sources with the largest flux densities.

Approximately half of the sources in our final sample have estimated redshifts available, either photometric or spectroscopic. This allows us to analyze the uniformity of radio source position angles in 3D space, but no strongly significant deviation from uniformity was found. We think it is more likely that the effect is caused by systematic effects, given the fact that the 2D analysis of the same reduced sample of sources still show an effect. However, the results are not straightforward to interpret due to the added uncertainties on the photometric redshifts, leaving no indisputable conclusion.

Understanding the systematic effect or physical effect that causes the observed alignment in different radio surveys is beyond the scope of this study, but should be investigated further. In particular, these subtle effects will be important for cosmological analyses with radio data, such as weak lensing studies with the Square Kilometer Array (e.g., Harrison et al. 2016; Bonaldi et al. 2016).

The number of sources considered here comprises less than 2% of the complete LoTSS survey. Hence, future studies by LOFAR should result in information about radio source alignments caused by substantially more subtle effects than we are presently able to determine. Additionally, the WEAVE-LOFAR project (Smith et al. 2016) will obtain over a million spectra of radio sources in LoTSS, which will allow for a much more detailed study of alignment in 3D space. This will provide the statistics needed to prove or disprove whether the alignment effect observed in this study is physical.

Acknowledgments

We kindly thank the anonymous referee for the comments and instructive insights. E. O. and R. J. vW. acknowledge support from the VIDI research programme with project number 639.042.729, which is financed by the Netherlands Organisation for Scientific Research (NWO). M. J. H. acknowledges support from STFC [ST/R000905/1]. W. L. W. acknowledges support from the ERC Advanced Investigator programme NewClusters 321271. W. L. W. also acknowledges support from the CAS-NWO programme for radio astronomy with project number 629.001.024, which is financed by the Netherlands Organisation for Scientific Research (NWO). KJD acknowledges support from the ERC Advanced Investigator programme NewClusters 321271. A. P. M. would like to acknowledge the support from the NWO/DOME/IBM programme “Big Bang Big Data: Innovating ICT as a Driver For Astronomy”, project #628.002.001. H. R. acknowledges support from the ERC Advanced Investigator programme NewClusters 321271. LOFAR is the Low Frequency Array designed and constructed by ASTRON. It has observing, data processing, and data storage facilities in several countries, which are owned by various parties (each with their own funding sources), and which are collectively operated by the ILT foundation under a joint scientific policy. The ILT resources have benefitted from the following recent major funding sources: CNRS-INSU, Observatoire de Paris and Université d’Orléans, France; BMBF, MIWF-NRW, MPG, Germany; Science Foundation Ireland (SFI), Department of Business, Enterprise and Innovation (DBEI), Ireland; NWO, The Netherlands; The Science and Technology Facilities Council, UK; Ministry of Science and Higher Education, Poland. This research made use of the Dutch national e-infrastructure with support of the SURF Cooperative (e-infra 180169) and the LOFAR e-infra group. The Jülich LOFAR Long Term Archive and the German LOFAR network are both coordinated and operated by the Jülich Supercomputing Centre (JSC), and computing resources on the Supercomputer JUWELS at JSC were provided by the Gauss Centre for Supercomputing e.V. (Grant CHTB00) through the John von Neumann Institute for Computing (NIC). This research made use of the University of Hertfordshire high-performance computing facility and the LOFAR-UK computing facility located at the University of Hertfordshire and supported by STFC [ST/P000096/1].

References

All Tables

All Figures

Fig. 1. Illustration of parallel transport. Vector v1 corresponding to a position angle θp1 and vector v2 corresponding to position angle θp2 are shown. In order to compare v1 to v2, v1 must be parallel-transported along the great circle indicated by the curve from location P1 to location P2. The transported vector is indicated by $v_1^{\prime }$ and the local basis vectors are denoted by (uδ, uα). In parallel transport, the angle α between the vector tangent to the sphere ut and the vector v remains fixed. Figure adapted from Jain et al. (2004).

In the text

	Fig. 2. Position angle distribution of the complete sample of selected sources.
In the text

	Fig. 3. Median angular radius corresponding to drawing a circular aperture around every source with an angular radius bound by the nth neighbor around that source.
In the text

	Fig. 4. Logarithm of the significance level for which uniformity in position angles as a function of the number of nearest neighbors n should be rejected for the sample of 7555 selected sources. The conversion to angular scale is shown in Fig. 3.
In the text

	Fig. 5. Logarithm of the significance level for the sample of sources split into four equal frequency total flux density bins, as defined in Table 1.
In the text

	Fig. 6. Distribution of position angles for the sample of 4212 sources that have a redshift measurement.
In the text

	Fig. 7. Logarithm of the significance level at which position angle uniformity should be rejected, as a function of the number of nearest neighbors n for the 4212 sources in that have redshifts available. The dashed line indicates the results of the 2D analysis and the solid line the results of the 3D analysis.
In the text

	Fig. 8. Logarithm of the significance level at which position angle uniformity should be rejected as a function of the number of nearest neighbors n for the 1051 sources with total flux density > 108 mJy and a redshift measurement. The dashed line indicates the results of the 2D analysis and the solid line the results of the 3D analysis.
In the text

	Fig. 9. Logarithm of the significance level at which position angle uniformity should be rejected as a function of the number of nearest neighbors n for the 523 sources with total flux density > 108 mJy and a spectroscopic redshift measurement. The dashed line indicates the results of the 2D analysis and the solid line the results of the 3D analysis.
In the text

	Fig. 10. 1σ uncertainties on the position angles of sources in our sample that are not classified by the LOFAR galaxy zoo, but by the source finder only.
In the text

	Fig. 11. Distribution of the 1000 simulated values Sn|MC and the highly significant values Sn for n = 467 and n = 916. Plotted for the initial sample of sources.
In the text

	Fig. 12. Percentage of the n = 101 nearest neighbors of every source in the photo-z perturbed sample of 1051 brightest sources that agrees with the nearest neighbors found in the unperturbed sample. See text for more details.
In the text

	Fig. 13. Scatter plot of the maximum dispersion measure (Eq. (13)) for every source, which indicates the strength of an alignment signal, of the selected sample of 7555 radio sources plotted for n = 700 as a function of right ascension and declination.
In the text