On the reliability of polarization estimation using Rotation Measure Synthesis
Abstract
We benchmark the reliability of the Rotation Measure (RM) synthesis algorithm using the 1005 Centaurus A field sources of Feain et al. (2009). The RM synthesis solutions are compared with estimates of the polarization parameters using traditional methods. This analysis provides verification of the reliability of RM synthesis estimates. We show that estimates of the polarization parameters can be made at lower S/N if the range of RMs is bounded, but reliable estimates of individual sources with unusual RMs require unconstrainted solutions and higher S/N.
We derive from first principles the statistical properties of the polarization amplitude associated with RM synthesis in the presence of noise. The amplitude distribution depends explicitly on the amplitude of the underlying (intrinsic) polarization signal. Hence it is necessary to model the underlying polarization signal distribution in order to estimate the reliability and errors in polarization parameter estimates. We introduce a Bayesian method to derive the distribution of intrinsic amplitudes based on the distribution of measured amplitudes.
The theoreticallyderived distribution is compared with the empirical data to provide quantitative estimates of the probability that an RM synthesis solution is correct as a function of S/N. We provide quantitative estimates of the probability that any given RM synthesis solution is correct as a function of measured polarized amplitude and the intrinsic polarization amplitude compared to the noise.
1 Introduction
The birefringence of the magnetized plasma that pervades interstellar and intergalactic space causes the observed linear polarization properties of most astrophysical radio sources to be strongly frequency dependent. In cold plasmas the effect of Faraday rotation causes the position angle of the linear polarization vector to increase by an amount . The magnitude of the rotation measure (RM) ranges between rad m for typical sources observed through the interstellar medium of our Galaxy (JohnstonHollitt, Hollitt & Ekers 2004; Schnitzeler 2010) to rad m for the cores of some radio sources to as high as rad m for the compact radio source at the centre of the Milky Way, Sgr A (Macquart et al. 2006, Maronne et al. 2007).
It is often desirable to derive the rotation measure along a particular line of sight and recover the linear polarization vector intrinsic to the radio source at the point of emission. This requires observations over a sufficiently large frequency range that the effect of Faraday rotation is measurable but with fine enough spectral resolution that Faraday rotation does not cause appreciable rotation of the polarization vector across individual spectral channels (bandwidth smearing). However, in many situations of astrophysical interest this tradeoff results in a situation where the signaltonoise ratio of the polarization measurement in each spectral channel is too low to compute either the rotation measure or intrinsic polarization vector without employing a method that simultaneously utilizes the measurements across the entire spectral band. This is particularly pertinent to observations made with nextgeneration backends, such as that of the EVLA (Perley et al. 2011) and the Compact Array Broadband Backend (Wilson et al. 2011).
The technique of rotation measure synthesis was developed by (Killeen et al. 1999) to search for polarization from Sgr A in situations where the large expected RM required such high frequency resolution that the signaltonoise ratio (S/N) per channel was too low. This technique, which was already implicit in the treatment of Burn (1966), was independently discovered by Brentjens & de Bruyn (2005) to handle the large Galactic rotations seen at low radio frequency. Over the past few years RM synthesis has emerged as the method of choice for the derivation of Faraday rotation measures from radio polarimetric data. RM synthesis utilizes both the amplitude and position angle of the polarization vector at all observed frequencies, and is capable of disentangling emission at multiple rotation measures (or Faraday depths) combined within the telescope beam. Because all the observed frequency dependent information is used without bandwidth smearing and without the nonlinear operations needed to calculate the variation in polarization position angle with frequency, it should also be the optimum method to extract other polarization information from noisy data. This technique is fundamentally different from previous methods which derived RMs by fitting position angle estimates at a usually small number of often widely spaced frequencies.
Despite its extensive current use, a number of questions remain concerning the practical application of RM synthesis to real datasets. This includes debate on the confidence of the results of the technique as a function of S/N (e.g. Law et al. 2011), and on the reliability of the technique when there are two emission regions with closelyspaced RM and polarization vectors within the same telescope beam (Farnsworth, Rudnick & Brown 2011). There is also confusion regarding the appropriate means of treating the effects of polarization bias, associated with the fact that the mean magnitude of the polarization vector is nonzero even if no polarized signal is present (Simmons & Stewart 1985).
The objective of this paper is to resolve these issues. The technique and nomenclature of RM synthesis is summarized in §2. In in §3 we proceed to examine several specific facets of the performance of rotation measure synthesis in detail. We utilize the 1005 source dataset of a field surrounding Centaurus A (Feain et al. 2009) to examine its performance in estimating polarization parameters.
In §4 we derive the distribution of polarization amplitudes expected from RM synthesis data as a function of the amplitude of the input polarization signal. We use this to demonstrate the derivation of the intrinsic distribution of polarization amplitudes within a sample. This is applied to the Cen A dataset to examine the polarization flux density distribution down to the submJy level; this presents a means of measuring the polarization distribution to well below the detection threshold of any given source, in a manner somwhat akin to that of “stacking” used in other surveys to achieve statistical detections of populations at detection thresholds below that of an individual source.
In §5 we examine the Cen A dataset for sources with unusual polarization behaviour, including sources with extreme or multiple RMÕs. Analysis of the changes in the goodness of the fit (as quantified by the reduced statistic) is used to identify the small number of sources with multiple RMs. The conclusions are presented in §6.
2 RM synthesis: terminology and summary
For the purposes of the following sections it is useful to briefly recapitulate how RM synthesis determines the rotation measure and polarization of a source subject to Faraday rotation. The input data for RM synthesis is a stream of and Stokes parameters, and where there are spectral channels corresponding to wavelengths . The polarization data can be equivalently represented as a series of complex numbers where . (Throughout the text we use to denote the complex polarization and to denote its amplitude.) RM synthesis works by winding up the measurements of and for a range of different trial rotation measures, . If the weights corresponding to the spectral channels are denoted , then the approximation to the Faraday dispersion function is
(1)  
where  (2) 
where and it is convenient (but not necessary) to choose as the weighted average of the channel values of over the observing band, i.e. . The approximation to the Faraday dispersion function is the convolution of the actual Faraday dispersion function, , with a function, , which describes the sensitivity of an observation to emission at each Faraday depth, :
(3) 
The function is known as the rotation measure transfer function (RMTF) and is computed as follows:
(4) 
In the simple case in which all channel weights are equal, the above relations simplify to
(5)  
(6) 
Now suppose a source contains polarization emission at one or more Faraday depths, , and associated with each depth is a polarization vector . It is clear from eq. (5) that for a dataset of polarization measurements for a source with Faraday depths , the sum over channels results in a coherent summing of the polarization vectors, and the amplitude of the function will be large. The tendency for this sum to add coherently for values of that match (i.e. correctly derotate) the polarization data means that peaks at the locations of the rotation measure(s) in the source, with the quantity giving the amplitude of the polarization vector that is rotating at each Faraday depth, . Thus the problem of measuring the Faraday depths and rotation polarization vectors that are significant in a given dataset is one of locating those peaks in the function that exceed some given threshold.
2.1 RM synthesis deconvolution implementation (CLEAN)
It is often necessary to recover the Faraday dispersion function, , by deconvolving from the RMTF, particularly if the source contains emission at multiple Faraday depths and if the S/N of the observations is high
The RM synthesis code was supplemented by several additional features. Because the functions and were computed on a grid of Faraday depths much finer than the width of , it is sometimes desirable to group together CLEAN components that are very tightly clustered together and identify them as members of a single component. By grouping components in this way we could obtain the estimate for each independent component. We introduced a Faraday depth bunching tolerance, , such that any two components with depths and were merged if . A tolerance rad m was found to be suitable for the dataset examined in the present paper.
The code output a list of clean components ranked in order of amplitude, together with their corresponding Faraday depth and polarization vector evaluated at the mean square wavelength . The code also computed the reduced associated with the polarization solution as a function of the number of clean components used. To be specific, for measurements with standard deviations and in and respectively, and clean components with Faraday depths and associated polarization vectors , the reduced associated with a fit to the data by including the first clean components is,
(7) 
where
(8) 
This statistic was found to be an excellent estimator of the number of significant clean components required to model the polarization behaviour of each source. It was also useful to identify cases in which RM synthesis solution differed significantly from the polarization behaviour of the source; this highlighted cases in which the polarization is not noiselike but where its behaviour is chiefly governed by a mechanism other than Faraday rotation (§5 discusses specific examples of this).
3 RM synthesis of Centaurus A data
In this section we compare the performance of RM synthesis against other techniques of determining the polarization and rotation measures of radio sources. As an initial step, we investigate the accuracy of the rotation measures reported by RM synthesis as a function of signaltonoise ratio. We then compare the rotation measures derived by RM synthesis against a technique that mimics the technique used to derive the rotation measures extracted from the NRAO VLA Sky Survey (Condon et al. 1998; Taylor, Stil & Sunstrum 2009) based on measurements of the polarization in two closely separated bands at 1.4 GHz. We then compare the performance of RM synthesis against a more sophisticated technique that performs a leastsquares fit to the polarization data.
This analysis is based on the ATCA observations of 1005 background sources contained in a deg field centred on the radio galaxy Centaurus A (Feain et al. 2009). The data have a bandwidth of 192 MHz divided into 24 8 MHz channels between 1296 and 1480 MHz. Because of the effect of the reweight option in Miriad a 5% correlation was introduced between channels resulting in 22.75 independent channels. The total bandwidth and spectral resolution of these data mean that the FWHM in Faraday rotation for a single polarization component is 280 rad m, and that it is sensitive to RMs with magnitudes less than 3500 rad m. The relatively low RM resolution in these data avoids the potential complexity in the interpretation of any sources with closelyspaced multiple RMs, as discussed by Farnsworth, Rudnick & Brown (2011). We allowed our RM synthesis code to search to RMs of rad m. The average rootmeansquare error in U and Q is 0.089 mJy and the corresponding average sigma in (which is the measure of noise used throughout this paper) is 0.12 mJy.
3.1 Confidence versus S/N based on RM synthesisderived parameters
We performed RM synthesis on the 1005 sources in the Cen A field to investigate the reliability of the derived Faraday depths as a function of S/N ratio. Several recent papers adopt a S/N of as the threshold below which the results of RM synthesis are regarded as unreliable and are not reported (Feain et al. 2009; Law et al. 2011). It is therefore propitious to investigate whether this threshold represents a sharp barrier to the believability of RM synthesis results, or whether results below this level only decline gradually in accuracy. The large size of our sample enables us to empirically derive the level of confidence associated with RM synthesis solutions over a large range of S/N ratios.
Figures 1 & 2 show the distribution of derived RMs for varying S/N thresholds. For the purposes of this analysis we have used only the single most significant Faraday depth derived by RM synthesis; the small number of sources in which there are multiple significant RMs has a negligible influence on the statistics (see §5). The signal was derived from the value of at the most significant Faraday depth. Some RFI had to be flagged in the time based data resulting a small variations in S/N between channels. This also varied from source to source. The noise, , is the taken to be the average rms noise per channel for each source divided by the squareroot of the weighted number of channels used. To simplify the theoretical analysis in the following sections we have assumed equal noise per channel, but it was then essential to correct for the effect of this nonuniform noise in the data by including the nonuniform weighting in the S/N estimates and in the estimates of the number of independent channels. We estimate this average rms noise per channel from the observations using the fact that the rms errors add in quadrature. Feain et al. (2009) provides specifics on the estimation of errors for the 8MHz and channels which were interpolated from the correlator output. The effect of RFI flagging which resulted in nonuniform noise varying from source to source had to be included in the analysis. Even though these effects are small and subtle it is essential to account for them when analysing the polarization statistics.
It is evident that all the observed RMs for high S/N detections are in the range to rad m, and we assume that the lower S/N detections are drawn from the same distribution. Polarization detections with RMs outside this range are most likely to be due to noise. At the lower S/N regime of interest here, the expected distribution is broadened by noise so we increase the false RM range by and use the range to rad m which is appropriate down to signal strengths of . If Faraday depths associated with false detections are evenly distributed over the depth search range – as appears to be the case based on Figures 1 & 2 – there is only a 4% chance that any given Faraday depth found by RM synthesis will fall in the range of expected real detections by accident. Table 1 gives the distribution of “true” and “false” polarization detections as a function of S/N.
Some RMs outside this range could be real (see §4.5) and while this does not significantly affect the statistical analysis, individual sources with high intrinsic RM are of interest (see §5). This population is rare in the flux density range of our sample since only one of the 309 sources detected above possesses an intrinsic RM magnitude larger than 200 rad m.
S/N  all  
Approx P(mJy)  
Total number  239  288  352  488  775  1003  1005 
Number with good RM  239  287  346  423  499  527  528 
Number with false RM  0  1  6  65  276  476  477 
Differential number  49  64  136  287  228  
Differential number with  48  59  77  76  28  
good RM  
Differential number with  0  1  5  59  211  200  
false RM 
The one high RM source in the 6S/N bin, 132446460218, merits comment. This source has a clear double peak in the RM synthesis spectrum. The two peaks have nearly equal polarization: one is at rad m but the second at rad m does fall within the expected RM range. The S/N difference in amplitude is %. As discussed at the end of §5, the statistical treatment of a second RM component in the presence of a significant polarized signal is well beyond the analysis in this paper and this may be a correct but rare detection.
3.2 Comparison with traditional RM methods
The traditional method to determine the RM has been to determine the position angle of polarization in each frequency channel and then to make a scalar least squares fit to the position angles. Figure 3 (top) shows the result of this analysis when applied to sources in the Cen A polarization catalogue. When compared to the RM synthesis solution for the same sources (bottom panel in Figure 3) it can be seen that a small number of spurious high RM values are obtained. Only those sources with integrated polarization signals S/N7, as identified by RM synthesis, were used in the analysis. It is evident that there is a tail in the distribution of the Faraday depths derived by fitting polarization position angles. Of the 251 sources with S/N7 used in the analysis, 16 of these fall outside the range rad m. This can be explained by noting that with the lower S/N in each channel the least squares fit for weak sources can find solutions which are randomly distributed over the RM space included in the least square solution. This high RM tail is reminiscent of the tail of high RMs sometimes seen in existing RM catalogues e.g. see the comparison of the NVSS catalogs and the Kronberg & NewtonMcGee compilation in Pshirkov et al. (2011). The small number of genuine high RMs seen in our dataset (%) indicates that this class of radio source is rare in surveys at 1.4 GHz, presumably because sources with high intrinsic RM are more likely to be strongly depolarized.
Taylor et al. (2009) report the rotation measures of 37,543 sources from the NVSS based on VLA snapshot images in two 42 MHz wide bands centred on 1364.9 MHz and 1435.1 MHz. It is opportune to investigate the accuracy of this catalog given that it currently represents one of the largest resources of rotation measures, and is rapidly being adopted as a de facto standard in analyses of the Galactic and extragalactic magnetic fields (HarveySmith, Madsen & Gaensler 2011; Law et al. 2011; Opperman et al. 2011; Govoni et al. 2010; McClureGriffiths et al. 2010; Schnitzeler 2010; Stasyszyn et al. 2010).
We used the Cen A polarization catalog to investigate the accuracy of this technique used to derive NVSS rotation measures. The data were obtained at a comparable frequency, and we have collapsed them into two spectral channels with similar centre frequencies of 1364 and 1436 MHz. As for the NVSS sample, the range of Faraday depths in our sample is sufficiently small that phase wrap ambiguities will not occur
3.3 Comparison of RM synthesis and least squares fits
Although RM synthesis outperforms traditional frequencychannelbased positionangle fitting algorithms such as the ones discussed above, it remains to be seen whether it makes the best use of the polarization information available. We address this by considering the output of a brute force minimum least squares search for detected polarization at any RM on all the sources, subject to the assumption that the source contains only a single polarized component that is assumed to have a flat spectrum. To be specific, we implemented a code that seeks the parameters , and (polarization amplitude, position angle and Faraday depth) that minimize the penalty function,
(9) 
where the model for the polarization of the th channel is,
(10) 
Figure 4 compares the solutions from the RM synthesis with the brute force least squares fit approach. The two methods agree very well in both RM and polarization amplitude. At high S/N the solutions are identical. The increased scatter at lower S/N is not surprising as the two methods find different solutions, but the gap (along the Faraday depth axis) surrounding the identical solution indicates that both methods converge to the identical solution if the solutions are close. There is no systematic change in the ratio of signal amplitude between the two methods with S/N, and the average ratio for the two methods is 1.0.
The behaviour of the solutions for the sources with a “false” detection based on the RM has a number of undesirable properties. The brute force solution has a much narrower distribution of RMs at low S/N ( 500 rad m). The amplitudes of the bad solutions are lower and the reduced worse indicating that it does not find the global minimum for noisedominated signals.
In §4.2 it is demonstrated that the noise does fall as for RM synthesis. The naïve expectation would be that it falls similarly for traditional methods, since the estimated standard deviation in the slope of a regression line scales as . However, in practice this dependence is more complicated for traditional methods in a way that is hard to analyse; phase ambiguities in the position angle introduce additional degrees of freedom into the fitting problem, and the fit error depends in detail on how the channels are spaced.
4 The statistics of polarized sources
Determination of the distribution of polarization amplitudes intrinsic to a set of sources involves disentangling the contribution of noise from the polarization amplitudes deduced by RM synthesis. This is most important for sources whose intrinsic polarization is near the noise level.
We are motivated to examine the noise statistics in detail with a view to accounting for its effects in the distribution of polarization amplitudes measured in a set of sources. Specifically, the objective is to evaluate the quantity , namely the probability of obtaining an observed polarization amplitude conditioned on the amplitude of the true intrinsic polarization, . This distribution function is the basis for deriving the probability distribution of intrinsic (i.e. “noisecorrected”) polarization amplitudes based on a set of measured polarization amplitudes determined by RM synthesis. We conclude this section by applying the formalism to investigate the distribution of polarization amplitudes in the Cen A dataset at low polarization amplitude levels.
4.1 The polarization distribution of noise
Here we derive the distribution of polarization amplitudes recovered by RM synthesis from pure noise. The input measurements and are therefore assumed to possess the properties of thermal noise, and are independent and normally distributed with zero mean and noise per spectral channel .
Recall that RM synthesis derotates the channelized measurements of and for a range of different trial rotation measures, . The resultant summed (complex) number is
(11) 
The polarization vector derived by RM synthesis is the value of for which the choice of gives a maximum in the value of . The task here is to determine the distribution of the amplitude of . If is given it is straightforward to recover the polarization distribution. However, an additional complication arises because the RM synthesis algorithm is allowed the extra degree of freedom to chose any value of which maximises . The freedom to rotate the measured polarization vectors over a range of trial values of ensures that the polarization vector amplitude returned by RM synthesis is at least as large as that found in the case. This is because the algorithm may derotate the measured polarization vectors with any value of which improves the coherence of the sum in eq. (11).
There is a direct analogy between the noise characteristics of RM synthesis and those encountered in interferometric radio imaging, which involves deriving the distribution of noise amplitudes from (complex) visibilities (see Section 9.3 of Thomson, Moran & Svenson 2001). In both cases one is attempting to “windup” complex valued quantities whose real and imaginary components are both normally distributed. In the case of RM synthesis, they are wound up for various values of Faraday depth. In the interferometric imaging case the complexvalued interferometric visibilites are searched across a range of fringe rotation values, corresponding to the possible position of the source. Thus the two treatments are mathematically identical. We elaborate on this connection in the analysis below.
Constant
First consider the simple, wellknown case in which we explicitly hold the rotation measure at a fixed value. Since the input data is purely noiselike the actual value of is unimportant: any noiselike polarization vector which is transformed by derotation to any nonzero value of still possesses the same statistical noiselike properties. Here, for the purpose of simplicity, we set . The polarization found by RM synthesis is,
(12) 
The statistics of the sum of independent, normally distributed random variates is also normally distributed, with standard deviation . Thus also follows a normal distribution with zero mean and standard deviation . An identical relation holds for . It is then straightforward to show that the probability distribution for the polarization amplitude when there is no signal present and is just a Rayleigh distribution,
(13) 
where we define . This result is well known in the analogous interferometric imaging case pertaining to noise in VLBI observations, where it corresponds to the statistics of the fringe amplitude for a specific fringe rotation when no signal is present (Thomson, Moran & Swenson 2001, p318).
Unconstrained
Now consider the general case in which RM synthesis performs a fit to the polarization data that is unconstrained in Faraday depth, . This consists of two separate steps. The first is to find the specific value of which maximizes the value of for a given dataset. The second step is to find the distribution of corresponding to these maxima.
The square of the polarization amplitude is
(14)  
The second line follows because the sine function has odd symmetry, so the imaginary part of every pair in the sum is cancelled by every pair .
The rotation measure should obey the following condition at the maximal value of ,
(15) 
Since the imaginary terms are antisymmetric with respect to exchange of the indices and , we see that imaginary part of this sum is zero for any value of . So only the real part of the sum yields a nontrivial constraint on :
(16) 
The foregoing relations show that the distribution of exhibits a detailed dependence on the spectral coverage (i.e. the values of ). This indicates that the exact location and amplitude of the noise peak depends on the , and that the noise distribution for any given spectral setup is best found empirically.
It is possible to construct a simple approximate analytic formulation independent of the detailed spectral coverage. This result has already been derived in the analogous case of the noise properties in VLBI images, where one is concerned with the statistics of the fringe amplitude across (Thomson, Moran & Swenson 2001, pp319211). We briefly summarize the argument here and place it in the context of RM synthesis.
The statistics of the polarization amplitude at any one value of obey a Rayleigh distribution. At any given value of the probability that we detect a polarization amplitude lower than some prescribed threshold amplitude, say , is the cumulative probability distribution corresponding to eq. (13), namely, . If the RM synthesis search is conducted over a large range of possible Faraday depths, , and the width of the amplitude of the RMTF is approximately , there are independent Faraday depth “pixels”, or trials, over which the RM search is being performed. The probability that we do not detect a polarization amplitude above some amplitude after these trials is . The probability of finding is just the culmulative probability distribution function of the amplitude of the polarization vector,
(17) 
where the subscript emphasises that the distribution is associated with a purely noiselike signal. The corresponding probability density function is,
(18) 
This expression describes the expected distribution of the polarization amplitude when the input data are purely noiselike and in which the rotation measure is searched over approximately independent Faraday depth pixels. It is identical the characteristics of the noise in a search for fringes from independent samples (Thomson, Moran & Swenson 2001, eq 9.57). The exact noise distribution pertaining to the Cen A dataset is discussed at the conclusion of this section. We preempt this by noting that we find empirically that the noise distribution for the Cen A observations (shown in Fig. 6) differs slightly from the analytic formula derived here.
The behaviour of the noise distribution has a number of important implications for the statistics of polarization amplitudes determined by RM synthesis. Equation (18) shows that the distribution of polarization amplitudes derived from noise narrows as increases (for constant ), and the peak of the distribution shifts to increasing values of . The location of the peak is determined by the solution to the transcendental equation,
(19) 
This has solutions and respectively for and .
Thus the threshold signaltonoise, , for the believability of a polarization signal extracted using RM synthesis exhibits an important, albeit weak, dependence on the number of independent trials, . The distribution peaks sharply on the low side of the distribution, implying that one expects very few polarization amplitudes with signaltonoise values less than once exceeds . For the analysis in this paper .
4.2 The polarization distribution of a nonzero polarization signal containing noise
There is an important distinction between the polarization distribution obtained from a purely noiselike signal from RM synthesis and the polarization distribution derived when an intrinsic (nonzero) polarization signal is present. The distinction arises because, in the absence of an intrinsic polarization signal, the RM synthesis algorithm picks out the polarization vector associated with the Faraday depth that maximises the length of the “Faradayderotated” noise vector. However, this freedom to optimise the length of the polarization noise vector is removed provided the intrinsic signal is sufficiently large that the RM synthesis code finds the highest amplitude at the correct Faraday depth of the intrinsic signal (i.e. there is no other Faraday depth at which a purely noiselike signal yields a higher polarization amplitude). In this instance the RM synthesis code is constrained to derotate the signal at the Faraday depth of the intrinsic polarization signal.
We compute the probability distribution of the polarization signal subject to the assumption that the intrinsic polarization is nonzero and that the RM synthesis code detrotates the signal at the Faraday depth associated with the intrinsic signal
(20) 
where and (as defined previously). The polarization vector, , therefore obeys the statistics of a random walk with a constant offset . Thus the polarization amplitude in the presence of an intrinsic signal follows the Rician distribution,
(21) 
where is the modified Bessel function of the first kind of order zero and the subscript denotes that the distribution is valid only when a nonzero polarization signal is present so that the signal is derotated at the Faraday depth corresponding to the intrinsic signal, . The associated cumulative probability distribution function (i.e. the probability of ) is
(22) 
where is Marcum’s Q function. This result is analogous to the probability distribution of the interferometric fringe amplitude in which one measures a nonzero signal in the presence of noise (Thomson, Moran & Swenson 2001, eq. 9.37).
Some remarks about the treatment of polarization noise bias are in order. The polarization amplitude distribution in eq. (21) demonstrates that the measured amplitude is, on average, biased by noise. In particular one has,
(23) 
Thus the appropriate noise bias to subtract from a single measurement of derived on the basis of RM synthesis is twice the channelintegrated noise variance, . While correct on average, bias subtraction is clearly inexact for any given source; in as much as the noise can assume a range of values, so is there a range of values of intrinsic that are consistent with a given measured . Our analysis shows that it is entirely inappropriate to simply remove a noise bias from the data and treat the resulting polarization amplitudes as if they were distributed with a Gaussian error distribution.
4.3 Likelihood of detection
In the two forgegoing subsections we have derived the probability distributions of the polarization amplitude both when a signal is absent, and when one is present and its amplitude is sufficiently large that it dominates the noise (i.e. provided exceeds some threshold).
Here we deduce the probability distribution of the polarization amplitude that is valid for all . The key to this derivation is that there is a nonzero probability that the value of drawn from the signal distribution is lower than the value that would be drawn from the noise distribution . Under this circumstance, the value of picked by the RM synthesis algorithm is the greater of the two values, namely that drawn from the noise distribution.
The probability of obtaining a polarization amplitude of is the sum of (i) the probability of drawing from the noise distribution and requiring that its value exceeds the amplitude drawn from the signal distribution, and (ii) the probability that was instead drawn from the signal distribution and requiring its value exceeds the amplitude drawn from the noise distribution. Thus, one obtains the combined probability distribution,
(24) 
This distribution is valid for any value of . It obeys the expected property that at low and the distribution is dominated by the noise distribution, , because is almost zero in this region, while it is dominated by at high because is almost zero in this domain.
By way of illustration, Fig. 5 compares the combined probability distribution for some specific choices of against and .
One can extend the foregoing logic to determine the likelihood that a signal found by RM synthesis at a given amplitude, , is a genuine signal. The probability that the signal has been derotated to the correct value of is,
(25) 
It is apparent that, even for small values of and , there is a nonzero probability that the RM synthesis code deduces the correct Faraday depth of the signal.
This is borne out by the Cen A data: close inspection of the lower panel in Fig. 1 indicates that even at S/N the RM synthesis algorithm preferentially finds signals within the range of physically acceptable values of .
4.4 Deduction of unbiased polarization amplitude source counts
The foregoing formalism provides an obvious prescription to determine the distribution of intrinsic polarization amplitudes, based on the measured distribution of polarization amplitudes as deduced by RM synthesis, namely :
(26) 
The solution of eq. (26) to find is formally an illposed problem, but is it of a form that is often encountered in numerical optimization theory.
It is instructive to demonstrate how the solution of eq.(26) is readily ammenable to numerical solution. The distribution of intrinsic polarization amplitudes, must be solved by sampling onto a grid with resolution, say, . We represent this distribution by the vector . The corresponding distribution of measured polarization amplitudes is discretized onto a grid in with identical resolution, and this is represented by the vector . If one similarly discretizes into a matrix whose th element corresponds to the parameters and , with resolution in both dimensions, the integral in eq.(26) can be cast in the form,
(27) 
This class of problems is ammenable to solution by iterative deconvolution and expectationmaximization algorithms.
A simpler but more robust approach is to forwardmodel the expected distribution to determine whether it is consistent with the observed amplitude distribution.
4.5 Application to Centaurus A data
As an illustration of the foregoing formalism, we apply it to determine the distribution of intrinsic polarization amplitudes in the Cen A dataset.
In order to model the effect of noise in this dataset precisely, we determined the noise distribution empirically based on the specific spectral coverage of the ATCA observations. The same RM synthesis code used to process the Cen A data was employed. We synthesised purely noiselike sources where, for each source, synthetic and data for each spectral channel were drawn from a normal distribution of zero mean and unit standard deviation, and assigned to spectral channels with frequencies identical to those used in the Cen A observations.
For each source the “signal” was extracted at the value of at which the amplitude in the Faraday dispersion spectrum was a maximum. The statistics of the polarization amplitude at a fixed Faraday depth, , were also recorded. The resulting probability distributions are shown in Fig. 6. There is excellent agreement between the distribution for the case and the Rayleigh distribution discussed in §4.1.1. The lower panel of Fig. 6 shows the amplitude distribution determined from the synthetic noise dataset when the RM synthesis algorithm is allowed to find a source anywhere over the Faraday depth range, rad m. One expects naïvely, on the basis of the width of the RMTF and the range of the search in Faraday depth, that there are independent Faraday depth trials over the search space. There is good agreement in the shape of the theoretically and empiricallyderived distributions for this value of . However, it is evident that the theoretical distribution peaks at a lower S/N than the empiricallyderived distribution. This cannot be accounted in terms of a change in since, although this would shift the peak of the distribution to higher S/N, it would also narrow the distribution to a value lower than is observed. We note that, from a purely empirical viewpoint, the two distributions agree closely if the S/N axis is scaled such that the recovered S/N values from the simulated sources are systematically higher than those expected in the theoretical distribution by 8.5%. The distribution corresponding to this is also shown in Fig. 6. However, the reason for such a S/N discrepancy is not obvious; the lower right panel of Fig. 4 shows that there is no systematic bias in the S/N of sources recovered by the RM synthesis technique and a bruteforce fit to the polarization data.
One can directly compare the polarization amplitude distribution for the sources in the Cen A field with “false” RMs (see §3.1) against the distribution, as determined by the synthetic purenoise sources. This comparison is shown in Fig. 7. In this plot we see excellent agreement between the polarization amplitude distribution of sources with false RMs and the expected noise distribution. However, there is an excess of sources in the region 4S/N5 which deserves further investigation. A likely explanation for this excess is that a number of these sources are mistakenly identified as false detections. It is possible that several contain real signal, but their S/N is sufficiently low that their Faraday depth is poorly localized by RM synthesis (i.e. they are only localized to about the rad m width of the RM transfer function), thus placing the detection RM outside the nominal range rad m considered to constitute a “good” detection. The excess of sources can be seen in Fig. 2 as the concentration of sources in the bin with RMs in the range rad m. The results of §4.2 (and Fig. 5) indicate that when some small signal is present the polarization amplitude distribution shifts to higher S/N. This suggests that some sources with supposedly “false” RMs do indeed contain intrinsic signal.
A related statistic is , the expected fraction of sources whose Faraday depths are identified correctly as a function of for a given . We compute this quantity using the numericallyderived distribution, . Figure 8 displays this likelihood. This gives a direct estimate that the probability is correct at any given RM. At S/N the probability of retrieving the correct Faraday depth is %, and it approaches 100% for S/N values between 5 and 6, depending on the value of the intrinsic polarization value, ,
It is interesting to compare the differential source counts from the Cen A data in Table 1 against this curve, despite uncertainty in the intrinsic polarization amplitudes that contribute to the Cen A data. At S/N values of 3, 4 and 5, the measured fractions of good RMs are 26%, 57% and 92% respectively, while the predicted fractions at these S/N values are in the ranges 2557%, 6193% and 94%100% respectively. In estimating the range of fractions, we have assumed that the S/N of the intrinsic polarization, , is between 2.5 and 5.
Forward Modelling of Polarization Source Counts
It is beyond the scope of this paper to implement a complete iterative deconvolution scheme to recover the intrinsic polarization source counts, , using eq.(26) in conjunction with .
We have instead developed a simple iterative forwardmodelling scheme in which the intrinsic polarized source count distribution is taken to be of the form,
(28) 
which is a broken power law distribution with index that turns over at a S/N of and scales with an index below this S/N break, and where is an overall normalization
(29) 
where we take the errors to be (and set when ). The quantity was searched for a minimum in terms of the parameters , , and . The minimum occurs at for , and the bestfitting solution is plotted over the measured polarization distribution in Figure 9. However, the minimum formally supports a range of indices in the range , and we note that the values of and and are all highly correlated. To illustrate this, Table 2 below lists the derived fit parameters for a range in the break S/N, .
0.0  1.12  0.96  –  119 
1.0  0.93  1.53  4  1.56 
2.0  0.88  1.90  1.01  8.09 
3.0  0.86  2.04  0.05  1.48 
3.72  0.85  2.17  0.30  2.60 
4.0  0.85  2.22  0.37  3.27 
5.0  0.87  2.40  0.54  7.39 
6.0  0.91  2.54  0.66  1.45 
7.0  0.95  2.64  0.74  2.50 
The bestfitting index of the slope of the differential counts in at high S/N, , is very close to the value of obtained by fitting just the stronger welldetected sources in the sample with in the range 0.6 to 6 mJy. This indicates that the source counts in continue with the same slope down to at 0.35 mJy, and possibly lower. Note that this forward modelling method does not require any assumptions about polarization amplitude bias at low S/N.
Following the technique used by Tucci et al. (2004) we can now estimate the change in fractional polarization as a function of flux density by comparing the source counts in and at the same source density. For both FIRST (White et al. 1997) and these Centaurus A field sources the differential counts in have a slope of in this flux range. The fractional polarization is 3.0% for the strongest sources (mJy) and this increases to 3.5% for sources below 20 mJy. This can be compared with the mean fractional polarization of ATLSB sources ( 0.4 mJy) of 4.3%. Taylor et al. (2007) make a detailed analysis of polarization in the Elais N1 field. They have 786 sources with 83 polarization detections. We have 1005 sources and at least 346 polarization detections in the same flux density range. We agree well in the source density and polarized fraction at the higher flux density but our fractional polarization for sources below 20 mJy (3.7%) is less than the increase to 4.8% seen by Taylor et al. (2007).
5 Sources with complex polarization
To evaluate the reality of sources which have multiple RMs we evaluated the reduced after RM synthesis identified each new CLEAN component

Some residual bad data points were readily identified in this process and these (6) sources do not appear in the list when the bad data is flagged out.

High S/N (3 sources). The high associated with some high S/N sources could be due to real but weak components with different RM which are only visible at high S/N. It may also result from instrumental dynamic range errors.

Multiple RM (8 sources). These are well fit by multiple components. Note that for sources with multiple RMs we are unable to uniquely separate components within our RM resolution of rad m.

Significant spectral change in across the band (4 sources). This is likely to be caused by regions with closely spaced RM and polarization vectors as discussed by Farnsworth et al. (2011), although it may in principle be caused by depolarization effects.
Component no.  RM (rad m)  reduced after inclusion  (mJy)  S/N 

1  120.4  4.5  1.41  15.0 
2  171.2  2.6  0.94  10.0 
3  323.2  1.6  0.54  5.8 
Component no.  RM (rad m)  reduced after inclusion  (mJy)  S/N 

1  74.8  2.7  13.1  66.8 
2  118.4  1.7  0.56  2.9 
3  686.8  1.1  0.45  2.3 
We illustrate these effects with two sources from this list:
131943445904 has three significant and comparable amplitude RM components which are just separated at our RM resolution. The three components are a good fit and reduce from 4.5 to 1.6. See Figure 10 and Table 3.
131713410934 is a high S/N case with a dominant RM components at rad m. It also has a barely resolved component at rad m which causes the slope in and a much weaker (3%) well resolved component at RM rad m. See Figure 11 and Table 4.
It is beyond the scope of this paper to present a detailed discussion of the sources with complex RM structure but we can make a few useful observations. Sources with RM structure less than our resolution (280 rad m) can still cause strong variations in polarization amplitude across the band and the RM clean process can fit this change in amplitude with components separated by approximately the HPBW. While this is still a valid indication of significant RM structure, the “super resolution” of our CLEAN process (i.e. the ability to resolve RM components within the width of the RMTF) is unlikely to find a unique solution for the actual RM components.
Although we do see clear evidence of multiple RMs in 12 of the 359 sources with S/N in polarization amplitude greater than 5, this is only 3% of our sample and the majority of sources are well fit by a single RM component.
To better quantify this fraction we can ask down to what S/N one can we detect a second RM component by examining to what degree the addition of the second component improves the fit to the spectropolarimetric data. We start with a model consisting of one polarized component and compare it against the measurement vector . Let us assume that the measurements in each spectral channel have the same error, . Then, if the vector contains components, one has
(30) 
Now let us introduce a second polarized component, that improves the fit and produces a new reduced of
(31) 
Subtracting the two foregoing expressions and, if the data are well explained by two polarized components, one has we find,
(32) 
where we use the fact that . Recognising that is the average squared amplitude of the second polarized component, say , we find,
(33) 
For our measurements, one has and we can see that, to improve the fit from to , one reqires a component with a S/N of .
We note that the statistical treatment for the detection of a second RM component in the presence of another strongly polarised signal is very different from the analysis of the detection of a polarised signal in the presence of noise. The perturbation that an additional polarised signal generates will be closer to the normal detection statistics so the weak multiple components we find at the S/N level are already significant.
6 Conclusions
Our major conclusions are as follows:

We have assessed the distribution of false and correct RM detections in the 1005source sample of Feain et al. (2009). This yields a quantitative estimate of the likelihood that a given RM, as found by RM synthesis, is correct, as a function of S/N. If S/N, we find the likelihood of finding a correct RM is about 52%. Restricting the range of possible RM solutions to the range 210 to 110 rad m (the range found for the high S/N sources) further decreases the likelihood of a false RM detection; there is only a 4% chance that any given Faraday depth found by RM synthesis in a search over the range rad m will fall in the range of expected real detections by accident.

There is no systematic difference in the polarization amplitude recovered by RM synthesis compared to a leastsquares fitting approach as a function of S/N, and the average ratio of the amplitudes found by the two methods is .

We have investigated the effect of noise in RM synthesis and examined the probability distribution of the polarization amplitude as determined by this algorithm. From this, we have developed a formalism to recover the distribution of intrinsic polarization amplitudes from a measured distribution which is affected by noise.
This enables us to recover the polarization distribution at S/N levels well below the confidence level for any single detection. We have applied this to examine the polarization counts of the sources in the Feain et al. (2009) Cen A catalog; the best fitting polarization source count distribution follows , which is consistent with found by fitting only to S/N detections, and with no turnover at an amplitude above 0.2 mJy.
A related point is that polarization “bias” should not be subtracted from the derived polarization amplitudes when analysing the distribution of . Subtraction of this bias is often at best misleading, since the distribution of is highly nonGaussian.

Preexisting analyses have used a S/N cutoff of as a limit to the believability of results from RM synthesis. We derive here the likelihood that a given RM detection is correct as a function of S/N and the intrinsic polarization amplitude () of the polarization signal. We also derive the criteria for the likelihood of the detection of a second RM component. For instance, a detection at a S/N is sufficient to believe the existence of a second component if the addition of the second component reduces the value of the fit to the polarization data from to .
Footnotes
 affiliation: ICRAR/Curtin University of Technology, Bentley, WA 6845, Australia; J.Macquart@curtin.edu.au
 affiliation: ICRAR/Curtin University of Technology, Bentley, WA 6845, Australia; J.Macquart@curtin.edu.au
 affiliation: CSIRO Astronomy and Space Science, P.O. Box 76, Epping, NSW 1710, Australia
 affiliation: CSIRO Astronomy and Space Science, P.O. Box 76, Epping, NSW 1710, Australia
 affiliation: School of Chemical & Physical Sciences, Victoria University of Wellington, PO Box 600, Wellington 6140, New Zealand
 When the S/N is sufficient, this enables us to “superresolve” and search for multiple peaks in a scale smaller than the width of .
 In the NVSS, the actual rotation measure of a source is related to the deduced value assuming no phase wraps, , by the relation, , where represents the number of integral phase wraps that occurred between the two observing frequencies. Taylor et al. (2009) assume , given that RMs of the magnitude required to effect a phase wrap are rare when the source is located off the Galactic plane and away from the inner Galaxy.
 The contribution of noise, at the Faraday depth of the intrinsic signal may, in principle, slightly alter the depth at which the peak amplitude occurs, and hence the Faraday depth reported by the RM synthesis algorithm. This may, in turn, influence the statistics of . However, this is a small effect in most practical situtations. If the Faraday depth derived by RM synthesis is sufficiently close to the correct depth it will not alter the statistics of ; this is the case provided that the S/N is sufficient that the error in Faraday depth is a small fraction of the width of the rotation measure transfer function, (i.e. the S/N ). Winding up an intrinsic polarization signal at a depth that is wrong by an amount yields an intrinsic signal of amplitude , so the error in the amplitude of the intrinsic polarization vector, , is quadratic in . This is a small quantity since, for S/N the peak is localised to within the width of , one has , so that . We see in the following subsection that the condition S/N is always satisfied.
 Since the distribution is zero for , the fact that diverges for small when does not pose a problem for our modelling. We are effectively insensitive to the behaviour of the distribution at .
 Recall that CLEAN recognised all components that were within 5 rad m of each other as being part of the same component.
References
 Brentjens, M.A. & de Bruyn, A.G. 2005, A&A, 441, 1217
 Burn, B. J. 1966, MNRAS, 133, 67
 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
 Farnsworth, D., Rudnick, L. & Brown, S. 2011, AJ, 41, 19
 Feain, I.J., Ekers, R.D., Murphy, T., Gaensler, B.M., Macquart, J.P., Norris, R.P., Cornwell, T.J., JohnstonHollitt, M., Ott, J. & Middelberg, E. 2009, ApJ, 707, 114
 Govoni F., Dolag, K., Murgia, M., Feretti, L., Schnidler, S., Giovannini, G., Boschin, W., Vacca, V. & Bonafede, A. 2010, A&A, 522, A105
 Heald, G., Braun, R. & Edmonds, R. 2009, A&A, 503, 409
 HarveySmith, L., Madsen, G.J. & Gaensler, B.M. 2011, ApJ, 736, 83
 Högbom, J.A. 1974, Astronomy & Astrophysics Supp, 15, 417
 JohnstonHollitt, M., Hollitt, C.P. & Ekers, R.D. 2004, in The Magnetized Interstellar Medium, eds. B. Uyaniker, W. Reich & W. Wielebinski, Copernicus GmbH, 13
 Killeen, N.E.B., Fluke, C.J., Zhao, JumHui, Ekers, R.D. 1999, MNRAS, submitted (MZ381)
 Law, C.J., Gaensler, B.M., Bower, G.C., Backer, D.C., Bauermeister, A., Croft, S., Forster, R., GutierrezKraybill, C., HarveySmith, L., Heiles, C., Hull, C., Keating, G., MacMahon, D., Whysong, D. Williams, P.K.G. & Wright, M. 2011, ApJ, 728, 57
 Macquart, J.P., Bower, G.C., Wright, M.C.H., Backer, D.C. & Falcke, H. 2006, ApJ, 646, L111
 Marrone D.P., Moran, J.M., Zhao, J.H. & Rao, R. 2007, ApJ, 654, L57
 McClureGriffiths N. M., Madsen, G. J., Gaensler, B. M., McConnell, D. & Schnitzeler, D. H. F. M. 2010, ApJ, 725, 275
 Oppermann, N., Junklewitz, H., Robbers, G. & Enßlin, T.A., A&A, 530, A98
 Perley, R.A., Chandler, C.J., Butler, B.J. & Wrobel, J.M. 2011, ApJ, 739, L1
 Schnitzeler, D.H.F.M. 2010, MNRAS, 409, L99
 Simmons, J.F.L. & Stewart, B.G. 1985, A&A, 142, 100
 Stasyszyn F., Nuza, S. E., Dolag, K., Beck, R. & Donnert, J. 2010, MNRAS, 408, 684
 Taylor, A.R., Stil, J.M., Grant, J.K., Landecker, T.L., Kothes, R., Reid, R.I., Gray, A.D., Scott, D., Martin, P.G., Boothroyd, A.I., Joncas, G., Lockman, F.J., English, J., Sajina, A. & Bond, J.R. 2007, ApJ, 666, 201
 Taylor, A.R., Stil, J.M. & Sunstrum, C. 2009, ApJ, 702, 1230
 Thompson, A. R., Moran, J. M., & Swenson, G. W., Jr. 2001, “Interferometry and synthesis in radio astronomy” by A. Richard Thompson, James M. Moran, and George W. Swenson, Jr. 2nd ed. New York : Wiley
 Tucci, M., MartnezGonzález, E., Toffolatti, L., GonzálezNuevo, J. & De Zotti, G. 2004, MNRAS, 349, 1267
 White, R.L., Becker, R.H., Helfand, D.J. & Gregg, M.D. 1997, ApJ, 475, 479
 Wilson, W.E. et al. 2011, MNRAS, 416, 832