2.1. Prior work

There are multiple ways to estimate Faraday complexity, including detecting non-linearity in $\chi(\lambda^2)$ (Goldstein & Reed Reference Goldstein and Reed1984), change in fractional polarisation as a function of frequency (Farnes, Gaensler, & Carretti Reference Farnes, Gaensler and Carretti2014), non-sinusoidal variation in fractional polarisation in Stokes Q and U (O’Sullivan et al. Reference O’Sullivan2012), counting components in the FDF (Law et al. Reference Law2011), minimising the Bayesian information criterion (BIC) over a range of simple and complex models (called ‘QU-fitting’; O’Sullivan et al. Reference O’Sullivan, Purcell, Anderson, Farnes, Sun and Gaensler2017), the method of Faraday moments (Anderson et al. Reference Anderson, Gaensler, Feain and Franzen2015; Brown Reference Brown2011), and deep convolutional neural network classifiers (CNNs; Brown et al. Reference Brown2018). See Sun et al. (Reference Sun2015) for a comparison of these methods.

The most common approaches to estimating complexity are QU-fitting (e.g. O’Sullivan et al. Reference O’Sullivan, Purcell, Anderson, Farnes, Sun and Gaensler2017) and Faraday moments (e.g. Anderson et al. Reference Anderson, Gaensler, Feain and Franzen2015). To our knowledge, there is currently no literature examining the accuracy of QU-fitting when applied to complexity classification specifically, though Miyashita et al. (Reference Miyashita, Ideguchi, Nakagawa, Akahori and Takahashi2019) analyse its effectiveness on identifying the structure of two-component sources. Brown (Reference Brown2011) suggested Faraday moments as a method to identify complexity, a method later used by Farnes et al. (Reference Farnes, Gaensler and Carretti2014) and Anderson et al. (Reference Anderson, Gaensler, Feain and Franzen2015), but again no literature examines the accuracy. CNNs are the current state-of-the-art with an accuracy of 94.9% (Brown et al. Reference Brown2018) on simulated ASKAP Band 1 and 3 data, and we will compare our results to this method.

2.2. Assumptions on Faraday dispersion functions

Before we can classify FDFs as Faraday complex or Faraday simple, we need to define FDFs and any assumptions we make about them. An FDF is a function that maps Faraday depth $\phi$ to complex polarisation. It is the distribution of Faraday depths in an observed polarisation spectrum. For a given observation, we assume there is a true, noise-free FDF F composed of at most two Faraday screens. This accounts for most actual sources (Anderson et al. Reference Anderson, Gaensler, Feain and Franzen2015) and extension to three screens would cover most of the remainder—O’Sullivan et al. (Reference O’Sullivan, Purcell, Anderson, Farnes, Sun and Gaensler2017) found that 89% of their sources were best explained by two or less screens, while the remainder were best explained by three screens. We model the screens by Dirac delta distributions:

(1) \begin{equation} F(\phi) = A_0 \delta(\phi - \phi_0) + A_1 \delta(\phi - \phi_1). \end{equation}

$A_0$ and $A_1$ are the polarised flux of each Faraday screen, and $\phi_0$ and $\phi_1$ are the Faraday depths of the respective screens. With this model, a Faraday simple source is one which has $A_0 = 0$ , $A_1 = 0$ , or $\phi_0 = \phi_1$ . By using delta distributions to model each screen, we are assuming that there is no internal Faraday dispersion (which is typically associated with diffuse emission rather than the mostly compact sources we expect to find in wide-area polarised surveys). F generates a polarised spectrum of the form shown in Equation (2):

(2) \begin{equation} P(\lambda^2) = A_0 e^{2i\phi_0\lambda^2} + A_1 e^{2i\phi_1\lambda^2}. \end{equation}

Such a spectrum would be observed as noisy samples from a number of squared wavelengths $\lambda^2_j, j \in [1, \dots, D]$ . We model this noise as a complex Gaussian with standard deviation $\sigma$ and call the noisy observed spectrum $\hat P$ :

(3) \begin{equation} \hat P(\lambda_j^2) \sim \mathcal N(P(\lambda^2_j), \sigma^2). \end{equation}

The constant variance of the noise is a simplifying assumption which may not hold for real data, and exploring this is a topic for future work. By performing RM synthesis (Brentjens & de Bruyn Reference Brentjens and de Bruyn2005) on $\hat P$ with uniform weighting we arrive at an observed FDF:

(4) \begin{equation} \hat F(\phi) = \frac{1}{D} \sum_{j = 1}^D \hat P(\lambda^2_j) e^{-2i\phi\lambda^2_j}. \end{equation}

Examples of F, $\hat F$ , P, and $\hat P$ for simple and complex observations are shown in Figures 1 and 2, respectively. Note that there are two reasons that the observed FDF $\hat F$ does not match the groundtruth FDF F. The first is the noise in $\hat P$ . The second arises from the incomplete sampling of $\hat P$ .

We do not consider external or internal Faraday dispersion in this work. External Faraday dispersion would broaden the delta functions of Equation (1) into peaks, and internal Faraday dispersion would broaden them into top-hat functions. All sources have at least a small amount of dispersion as the Faraday depth is a bulk property of the intervening medium and is subject to noise, but the assumption we make is that this dispersion is sufficiently small that the groundtruth FDFs are well-modelled with delta functions. Faraday thick sources would also invalidate our assumptions, and we assume that there are none in our data as Faraday thickness can be consistent with a two-component model depending on the wavelength sampling (e.g. Ma et al. Reference Ma, Mao, Stil, Basu, West, Heiles, Hill and Betti2019; Brentjens & de Bruyn Reference Brentjens and de Bruyn2005). Nevertheless some external Faraday dispersion would be covered by our model, as depending on observing parameters Faraday thick sources may appear as two screens (Van Eck et al. Reference Van Eck2017).

To simulate observed FDFs we follow the method of Brown et al. (Reference Brown2018), which we describe in Appendix E.

3.1. Features

Our features are based on a simple idea: all simple FDFs look essentially the same, up to scaling and translation, while complex FDFs may deviate. A noise-free peak-normalised simple FDF $\hat F_{\mathrm{simple}}$ has the form

(5) \begin{align} \hat F_{\mathrm{simple}}(\phi;\ \phi_s) &= R(\phi - \phi_s), \end{align}

where R is the rotation measure spread function (RMSF), the Fourier transform of the wavelength sampling function which is 1 at all observed wavelengths and 0 otherwise. $\phi_s$ traces out a curve in the space of all possible FDFs. In other words, $\hat F_{\mathrm{simple}}$ is a manifold parametrised by $\phi_s$ . Our features are derived from relating an observed FDF to the manifold of simple FDFs (the ‘simple manifold’). We measure the distance of an observed FDF to the simple manifold using distance measure $D_f$ , that take all values of the FDF into account:

(6) $${\varsigma _f}(\hat F) = {\min _{{\phi _s} \in }}{D_f}(\hat F(\phi )||{\hat F_{{\rm{simple}}}}(\phi ;\;{\phi _s})).$$

We propose two distances that have nice properties:

• invariant over changes in complex phase,

• translationally invariant in Faraday depth,

• zero for Faraday simple sources (i.e. when $A_0 = 0$ , $A_1 = 0$ , or $\phi_0 = \phi_1$ ) when there is no noise,

• symmetric in components (i.e. swapping $A_0 \leftrightarrow A_1$ and $\phi_0 \leftrightarrow \phi_1$ should not change the distance),

• increasing as $A_0$ and $A_1$ become closer to each other, and

• increasing as screen separation $|\phi_0 - \phi_1|$ increases over a large range.

Our features are constructed from this distance and its minimiser. In other words we look for the simple FDF $\hat{F}_{\mathrm{simple}}$ that is ‘closest’ to the observed FDF $\hat{F}$ . The minimiser $\phi_s$ is the Faraday depth of the simple FDF.

While we could choose any distance that operates on functions, we used the 2-Wasserstein ( $W_2$ ) distance and the Euclidean distance. The $W_2$ distance operates on probability distributions and can be thought of as the minimum cost to ‘move’ one probability distribution to the other, where the cost of moving one unit of probability mass is the squared distance it is moved. Under $W_2$ distance, the minimiser $\phi_s$ in Equation (6) can be interpreted as the Faraday depth that the FDF $\hat F$ would be observed to have if its complexity was unresolved (i.e. the weighted mean of its components). The Euclidean distance is the square root of the least squares loss which is often used for fitting $\hat{F}_{\mathrm{simple}}$ to the FDF $\hat F$ . Under Euclidean distance, the minimiser $\phi_s$ is equivalent to the depth of the best-fitting single component under assumption of Gaussian noise in $\hat F$ . We calculated the $W_2$ distance using Python Optimal Transport (Flamary & Courty Reference Flamary and Courty2017), and we calculated the Euclidean distance using scipy.spatial.distance.euclidean (Virtanen et al. Reference Virtanen2020). Further intuition about the two distances is provided in section 3.2.

We denote by $\phi_w$ and $\phi_e$ , the Faraday depth of the simple FDF that minimises the respective distances (2-Wasserstein and Euclidean).

\begin{align*} \phi_w &= \underset{\phi_w}{\mathrm{argmin}}\ D_{W_2}{(\hat F(\phi)||}{\hat F_{\mathrm{simple}}(\phi;\ \phi_w))},\\ \phi_e &= \underset{\phi_e}{\mathrm{argmin}}\ D_E{(\hat F(\phi)||}{\hat F_{\mathrm{simple}}(\phi;\ \phi_e))}. \end{align*}

These features are depicted on an example FDF in Figure 3. For simple observed FDFs, the fitted Faraday depths $\phi_w$ and $\phi_e$ both tend to be close to the peak of the observed FDF. However, for complex observed FDFs, $\phi_w$ tends to be at the average depth between the two major peaks of the observed FDF, being closer to the higher peak. For notation convenience, we denote the Faraday depth of the observed FDF that has largest magnitude as $\phi_a$ , i.e.

\begin{equation*} \phi_a = \underset{\phi_a}{\mathrm{argmax}}\ |\hat F(\phi_a)|.\\ \end{equation*}

Note that in practice $\phi_a \approx \phi_e$ . For complex observed FDFs, the values of Faraday depths $\phi_w$ and $\phi_a$ tend to differ (essentially by a proportion of the location of the second screen). The difference between $\phi_w$ and $\phi_a$ therefore provides useful information to identify complex FDFs. When the observed FDF is simple, the 2-Wasserstein fit will overlap significantly, hence the observed magnitudes $\hat F(\phi_w)$ and $\hat F(\phi_a)$ will be similar. However, for complex FDFs $\phi_w$ and $\phi_a$ are at different depths, leading to different values of $\hat F(\phi_w)$ and $\hat F(\phi_a)$ . Therefore, the magnitudes of the observed FDFs at the depths $\phi_w$ and $\phi_a$ indicate how different the observed FDF is from a simple FDF.

Figure 3. An example of how an observed FDF $\hat F$ relates to our features. $\phi_w$ is the $W_2$ -minimising Faraday depth, and $\phi_a$ is the $\hat F$ -maximising Faraday depth (approximately equal to the Euclidean-minimising Faraday depth). The remaining two features are the $W_2$ and Euclidean distances between the depicted FDFs.

In summary, we provide the following features to the classifier:

• $\log |\phi_w - \phi_a|$ ,

• $\log \hat F(\phi_w)$ ,

• $\log \hat F(\phi_a)$ ,

• $\log D_{W_2}{(\hat F(\phi)}||{\hat F_{\mathrm{simple}}(\phi;\ \phi_w))}$ ,

• $\log D_{E}{(\hat F(\phi)}||{\hat F_{\mathrm{simple}}(\phi;\ \phi_e))}$ ,

where $D_E$ is the Euclidean distance, $D_{W_2}$ is the $W_2$ distance, $\phi_a$ is the Faraday depth of the FDF peak, $\phi_w$ is the minimiser for $W_2$ distance, and $\phi_e$ is the minimiser for Euclidean distance.

3.2. Interpreting distances

Interestingly, in the case where there is no RMSF, Equation (6) with $W_2$ distance reduces to the Faraday moment already in common use:

(7) \begin{align} \zeta_{W_2}(F) &= \min_{\phi_w \in \mathbb{R}} D_{W_2}{(F(\phi)}||{F_{\mathrm{simple}}(\phi;\ \phi_w))}\end{align}

(8) \begin{align} &= \left(\frac{A_0A_1}{(A_0 + A_1)^2} (\phi_0 - \phi_1)^2\right)^{1/2}. \end{align}

See Appendix A for the corresponding calculation. In this sense, the $W_2$ distance can be thought of as a generalised Faraday moment, and conversely an interpretation of Faraday moments as a distance from the simple manifold in the case where there is no RMSF. Euclidean distance behaves quite differently in this case, and the resulting distance measure is totally independent of Faraday depth:

(9) \begin{align} \zeta_{E}(F) &= \min_{\phi_e \in \mathbb{R}} D_E{(F(\phi)}||{F_{\mathrm{simple}}(\phi;\ \phi_e))}\end{align}

(10) \begin{align} &= \sqrt{2} \frac{\min(A_0, A_1)}{A_0 + A_1}. \end{align}

See Appendix B for the corresponding calculation.

3.3. Classifiers

We trained two classifiers on simulated observations using these features: logistic regression (LR) and extreme gradient boosted trees (XGB). These classifiers are useful together for understanding Faraday complexity classification. LR is a linear classifier that is readily interpretable by examining the weights it applies to each feature, and is one of the simplest possible classifiers. XGB is a powerful off-the-shelf non-linear ensemble classifier, and is an example of a decision tree ensemble which are widely used in astronomy (e.g. Machado Poletti Valle et al. Reference Machado Poletti Valle, Avestruz, Barnes, Farahi, Lau and Nagai2020; Hložek et al. Reference Hložek2020). We used the scikit-learn implementation of LR and we use the XGBoost library for XGB. We optimised hyperparameters for XGB using a fork of xgboost-tuner Footnote ^a as utilised by Zhu, Ong, & Huttley (Reference Zhu, Ong and Huttley2020). We used $1\,000$ iterations of randomised parameter tuning and the hyperparameters we found are tabulated in Table C.1. We optimised hyperparameters for LR using a fivefold cross-validation grid search implemented in sklearn.model_selection.GridSearchCV. The resulting hyperparameters are tabulated in Table D.1 in Appendix C.

4.1. Data

4.1.1. Simulated training and validation data

Our classifiers were trained and validated on simulated FDFs. We produced two sets of simulated FDFs, one for comparison with the state-of-the-art method in the literature and one for application to our observed FDFs (described in Section 4.1.2). We refer to the former as the ‘ASKAP’ dataset as it uses frequencies from the Australian Square Kilometre Array Pathfinder 12-antenna early science configuration. These frequencies included 900 channels from 700 to $1\,300$ and $1\,500$ to $1\,800$ MHz and were used to generate simulated training and validation data by Brown et al. (Reference Brown2018). We refer to the latter as the ‘ATCA’ dataset as it uses frequencies from the 1 to 3 GHz configuration of the ATCA. These frequencies included 394 channels from 1.29 to 3.02 GHz and match our real data. We simulated Faraday depths from $-50$ to 50 rad m^–2 for the ‘ASKAP’ dataset (matching Brown) and $-500$ to 500 for the ‘ATCA’ dataset.

For each dataset, we simulated $100\,000$ FDFs, approximately half simple and half complex. We randomly allocated half of these FDFs to a training set and reserved the remaining half for validation. Each FDF had complex Gaussian noise added to the corresponding polarisation spectrum. For the ‘ASKAP’ dataset, we sampled the standard deviation of the noise uniformly between 0 and $\sigma_{\max} = 0.333$ , matching the dataset of Brown et al. (Reference Brown2018). For the ‘ATCA’ dataset, we fit a log-normal distribution to the standard deviations of O’Sullivan’s data (O’Sullivan et al. Reference O’Sullivan, Purcell, Anderson, Farnes, Sun and Gaensler2017) from which we sampled our values of $\sigma$ :

(11) \begin{equation} \sigma \sim \frac{1}{0.63 \sqrt{2 \pi} \sigma} \exp \left(-\frac{\log\left(50 \sigma - 0.5\right)^2}{2 \times 0.63^2}\right). \end{equation}

4.2. Results on ‘ASKAP’ dataset

The accuracy of the LR and XGB classifiers on the ‘ASKAP’ testing set was 94.4 and 95.1%, respectively. The rates of true and false identifications are summarised in Table 1. These results are very close to the CNN presented by Brown et al. (Reference Brown2018), with a slightly higher true negative rate and a slightly lower true positive rate (recalling that positive sources are complex, and negative sources are simple). The accuracy of the CNN was 94.9, slightly lower than our XGB classifier and slightly higher than our LR classifier. Both of our classifiers therefore produce similar classification performance to the CNN, with faster training time and easier interpretation.

4.3. Results on ‘ATCA’ dataset

The accuracy of the LR and XGB classifiers on the ‘ATCA’ dataset was 89.2 and 90.5%, respectively. The major differences between the ‘ATCA’ and the ‘ASKAP’ experiments are the range of the simulated Faraday depths and the distribution of noise levels. The ‘ASKAP’ dataset, to match past CNN work, only included depths from $-50$ to 50 rad m^–2, while the ‘ATCA’ dataset includes depths from $-500$ to 500 rad m^–2. The rates of true and false identifications are again shown in Table 1.

Table 1. Confusion matrix entries for LR and XGB on ‘ASKAP’ and ‘ATCA’ simulated datasets, and the CNN confusion matrix entries are adapted from Brown et al. (Reference Brown2018)

As we know the true Faraday depths of the components in our simulation, we can investigate the behaviour of these classifiers as a function of physical properties. Figure 4 shows the mean classifier prediction as a function of component depth separation and minimum component amplitude. This is tightly related to the mean accuracy, as the entire plot domain contains complex spectra besides the left and bottom edge: by thresholding the classifier prediction to a certain value, the accuracy will be 100% on the non-edge for all sources with higher prediction values.

Figure 4. Mean prediction as a function of component depth separation and minimum component amplitude for (a) XGB and (b) LR.

Figure 5. Principal component analysis for simulated data (coloured dots) with observations overlaid (black-edged circles). Observations are coloured by their XGB- or LR-estimated probability of being complex, with blue indicating ‘most simple’ and pink indicating ‘most complex’.

4.4. Results on observed FDFs

We used the LR and XGB classifiers, which were trained on the ‘ATCA’ dataset to estimate the probability that our 142 observed FDFs (Section 4.1.2) were Faraday complex. As these classifiers were trained on simulated data, they face the issue of the ‘domain gap’: the distribution of samples from a simulation differs from the distribution of real sources, and this affects performance on real data. Solving this issue is called ‘domain adaptation’ and how to do this is an open research question in machine learning (Zhang Reference Zhang2019; Pan & Yang Reference Pan and Yang2010). Nevertheless, the features of our observations mostly fall in the same region of feature space as the simulations (Figure 5) and so we expect reasonably good domain transfer.

Two apparently complex sources in the Livingston sample are classified as simple with high probability by XGB. These outliers are on the very edge of the training sample (Figure 5) and the underdensity of training data here is likely the cause of this issue. LR does not suffer the same issue, producing plausible predictions for the entire dataset, and these sources are instead classified as complex with high probability.

With a threshold of 0.5, LR predicted that 96 and 83% of the Livingston and O’Sullivan sources were complex, respectively. This is in line with expectations that the Livingston data should have more Faraday complex sources than the O’Sullivan data due to their location near the Galactic Centre. XGB predicted that 93 and 100% of the Livingston and O’Sullivan sources were complex, respectively. Livingston et al. (Reference Livingston, McClure-Griffiths, Gaensler, Seta and Alger2021) found that 90% of their sources were complex, and O’Sullivan et al. (Reference O’Sullivan, Purcell, Anderson, Farnes, Sun and Gaensler2017) found that 64% of their sources were complex. This suggests that our classifiers are overestimating complexity, though it could also be the case that the methods used by Livingston and O’Sullivan underestimate complexity. Modifying the prediction threshold from 0.5 changes the estimated rate of Faraday complexity, and we show the estimated rates against threshold for both classifiers in Figure 6. We suggest that this result is indicative of our probabilities being uncalibrated, and a higher threshold should be chosen in practice. We chose to keep the threshold at 0.5 as this had the highest accuracy on the simulated validation data. The very high complexity rates of XGB and two outlying classifications indicate that the XGB classifier may be overfitting to the simulation and that it is unable to generalise across the domain gap.

Figure 6. Estimated rates of Faraday complexity for the Livingston and O’Sullivan datasets as functions of threshold. The horizontal lines indicate the rates of Faraday complexity estimated by Livingston and O’Sullivan respectively.

Figures D.1 and D.2 show every observed FDF ordered by estimated Faraday complexity, alongside the models predicted by Livingston and O’Sullivan et al. (Reference O’Sullivan, Purcell, Anderson, Farnes, Sun and Gaensler2017), for LR and XGB, respectively. There is a clear visual trend of increasingly complex sources with increasing predicted probability of being complex.