28 -- 40 GHz variability and polarimetry of bright compact sources in the QUIJOTE cosmological fields
Yvette C. Perrott, Marcos López-Caniego, Ricardo T. Génova-Santos, Jose Alberto Rubiño-Martín, Mark Ashdown, Diego Herranz, Anne Lähteenmäki, Anthony N. Lasenby, Carlos H. López-Caraballo, Frédérick Poidevin, Merja Tornikoski
MMNRAS , 1–17 (2021) Preprint 10th February 2021 Compiled using MNRAS L A TEX style file v3.0
28 – 40 GHz variability and polarimetry of bright compact sources in theQUIJOTE cosmological fields
Yvette C. Perrott , ★ , Marcos López-Caniego ,Ricardo T. Génova-Santos , , Jose Alberto Rubiño-Martín , , Mark Ashdown , ,Diego Herranz , Anne Lähteenmäki , , Anthony N. Lasenby , ,Carlos H. López-Caraballo , , , Frédérick Poidevin , , Merja Tornikoski Astrophysics Group, Cavendish Laboratory, JJ Thomson Avenue, Cambridge CB3 0HE, UK School of Chemical and Physical Sciences, Victoria University of Wellington, PO Box 600, Wellington 6140, New Zealand European Space Agency, ESAC, Camino bajo del Castillo, s/n, Urbanización Villafranca del Castillo, Villanueva de la Cañada, E-28692 Madrid, Spain Instituto de Astrofísica de Canarias (IAC), E-38205 La Laguna, Tenerife, Spain Universidad de La Laguna, Dpto. Astrofísica, E-38206 La Laguna, Tenerife, Spain Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK Instituto de Física de Cantabria (CSIC-Universidad de Cantabria), Avda. de los Castros s/n, E-39005 Santander, Spain Aalto University Metsähovi Radio Observatory, Metsähovintie 114, 02540 Kylmälä, Finland Aalto University Department of Electronics and Nanoengineering, P.O. BOX 15500, FI-00076 AALTO, Finland INFN Sezione di Milano, Via Celoria 16, Milano, Italy
Accepted 2021 February 02. Received 2021 February 01; in original form 2020 February 03
ABSTRACT
We observed 51 sources in the Q-U-I JOint TEnerife (QUIJOTE) cosmological fields which were brighter than 1 Jy at 30 GHzin the
Planck
Point Source Catalogue (version 1), with the Very Large Array at 28 – 40 GHz, in order to characterise theirhigh-radio-frequency variability and polarization properties. We find a roughly log-normal distribution of polarization fractionswith a median of 2%, in agreement with previous studies, and a median rotation measure (RM) of ≈ − with oneoutlier up to ≈ − which is among the highest RMs measured in quasar cores. We find hints of a correlation betweenthe total intensity flux density and median polarization fraction. We find 59% of sources are variable in total intensity, and 100%in polarization at 3 𝜎 level, with no apparent correlation between total intensity variability and polarization variability. Thisindicates that it will be difficult to model these sources without simultaneous polarimetric monitoring observations and they willneed to be masked for cosmological analysis. Key words: (cosmology:) cosmic background radiation – cosmology: observations – radio continuum: general – (galaxies:) quasars: general
The Q-U-I Joint TEnerife (QUIJOTE) experiment (Rubiño-Martínet al. 2010) aims to detect inflationary B-modes using CMB ob-servations at 31 and 42 GHz with supporting observations between10 – 20 GHz to characterise foregrounds. The low angular resolu-tion ( ≈ ◦ ) of the instrument means that polarized extragalactic radiosources, a major contaminant at small angular scales (e.g. Puglisi etal. 2018), cannot be accurately characterised by the experiment andancillary data must be used to either subtract or mask them. Ideallythese sources would be monitored at higher resolution simultaneouslywith the cosmological observations as for the Very Small Array ex-periment (Grainge et al. 2003), however since this was not possiblewe investigated the long-term properties of the bright sources in the ★ E-mail: [email protected] field to assess the level of contamination an inaccurate subtractionwould introduce. The polarization and variability properties of ra-dio sources at these frequencies are not well-studied; some of therelevant existing studies are listed in Table 1. Polarimetric very-long-baseline interferometric monitoring studies exist at 43 GHz, e.g. Parket al. (2018), Algaba, Gabuzda, & Smith (2011), Algaba, Gabuzda,& Smith (2012), Jorstad et al. (2007), Lister (2001); however thesefocus on spatially resolving properties of the sources rather thaninvestigating their integrated properties, which are of relevance toexperiments such as QUIJOTE. At the other end of the resolutionscale, Bonavera et al. (2017) stack
Planck data to statistically detectpolarization from known source positions and constrain average po-larization properties, but are too limited by sensitivity to investigatevariability. Galluzzi & Massardi (2016), Galluzzi et al. (2017) andGalluzzi et al. (2018) detect a large sample of sources between 2 –38 GHz, however they only have repeated observations for a small © a r X i v : . [ a s t r o - ph . C O ] F e b Y. Perrott et al. subsample of sources and focus on broadband polarimetry rather thanvariability. It is clear that there are very few studies addressing theintegrated variability of bright, polarized sources at 30 – 40 GHz; wetherefore observed the brightest sources in the QUIJOTE fields usingthe Very Large Array (VLA) to determine at what level in intensitya source could be safely subtracted in polarization.
MNRAS , 1–17 (2021) GH z v a r i ab ilit y andpo l a r i m e t r y o f b r i gh t c o m pa c t s ou r ce s i n t h e QU IJ O T E c o s m o l og i c a l fi e l d s Table 1.
Characteristics of polarimetric studies of radio sources at frequencies overlapping the 30 – 40 GHz band of this study.Reference Instrument Frequencies Epochs Number ofsources CommentRudnick et al.(1985), Jones etal. (1985) 1.4 – 90 GHz (8 bands) Kitt Peak 10.6 m(for 31 GHz) 6 total, 2 at31 GHz ≈
20 Strong, flat-spectrum sourcesLópez-Caniego etal. (2009) WMAP 23, 33, 41, 61 GHz Single 5-year aver-age 41 (30) at 99%significance at 33(41) GHz Survey maps correlated with total in-tensity source catalogueJackson et al.(2010), Battye etal. (2011) VLA 8, 22, 43 GHz 1 203 Brighter than 1 Jy in the WMAP22 GHz catalogueSajina et al. (2011) VLA 4, 8, 22, 43 GHz 1 at 4, 8, 43 GHz;2 at 22 GHz 159 Radio galaxies with Australia Tele-scope 20 GHz (AT20G) survey fluxdensities >
40 mJy in Atacama Cos-mology Telescope cosmological fieldKravchenko, Cot-ton, & Kovalev(2015) VLA 1.4 – 43 GHz (8bands) 3 7 Bright calibration sourcesHuffenberger et al.(2015) Q/U Imaging Experiment 43, 95 GHz Single 25-monthaverage 13 at S/N > . >
200 mJy at 20 GHz,Southern Ecliptic Pole regionBonavera et al.(2017), PlanckCollaborationXXVI (2016)
Planck
30 – 353 GHz 4 years, averaged 1560 (stacked) PCCS2 30 GHz sources >
427 mJy M N R A S , ( ) Y. Perrott et al.
The paper is organised as follows. In Section 2 we introduce oursample; in Section 3 we describe our observations, data reductionand analysis methods. In Section 4 we summarize our results andinvestigate correlations between the measured quantities. In Section 5we investigate in more detail some interesting or anomalous cases.In Section 6 we discuss implications for QUIJOTE and in Section 7we conclude.Throughout this paper we will use the following conventions.When fitting a power-law as a function of frequency we will usethe convention 𝑆 ∝ 𝜈 𝛼 (i.e. 𝛼 > I , Q , U to denote fluxdensity in total intensity and Stokes Q and U respectively. We willuse 𝑞 = Q/I and 𝑢 = U/I to indicate the fractional quantities andsimilarly P = √︁ Q + U and 𝑝 = P/I to denote the total and frac-tional linear polarization. Polarization angle will be indicated as 𝛷 ,i.e. tan ( 𝛷 ) = 𝑢 / 𝑞 . Circular polarization will be assumed to be zerothroughout. We selected the 54 sources with flux densities greater than 1 Jy at30 GHz in the first
Planck
Catalogue of Compact Sources (PCCS1,Planck Collaboration XXVIII 2014) lying within 30 ◦ of the centresof the three QUIJOTE cosmological fields at RA, 𝛿 = (00h40m,+25 ◦ ), (09h40m, +45 ◦ ) and (16h20m, +50 ◦ ) respectively. We notethat the second, more-complete version of the catalogue (PCCS2,Planck Collaboration XXVI 2016) was not available at the timeof sample selection however completeness is not an issue for thishigh-signal-to-noise-ratio (SNR) sample since at 30 GHz PCCS1 isestimated to be 90% complete at 575 mJy. Had we selected fromPCCS2, ≈
10% of the sample would have been different due to thesource flux densities varying over/under the 1 Jy threshold.Three of the sources in the sample were excluded due to largeangular size, leaving a sample of 51 which were observed with theVLA at frequencies between 28 and 40 GHz, with 41 sources havingtwo epochs of observation. The sources, as expected given the high-frequency selection, are mostly blazars and have identifications inthe 5th edition of the Roma-BZCAT Multifrequency Catalogue ofBlazars (Massaro et al. 2015); these are listed in Table 2 alongwith coordinates and other identifications taken from the SIMBADAstronomical Database . The two sources not in the BZCAT are aLINER-type AGN and a Seyfert-1 galaxy. The sources were observed on the dates listed in Table 3. Observa-tions were carried out with a custom correlator configuration whichallowed simultaneous observing at two frequency bands, 28.5 – 32.4and 35.5 – 39.4 GHz. Each band was then divided into 32 spectralwindows, each with 64 channels of 2 MHz width.We used either 3C 286 or 3C 48 as the flux density calib-rator, both in total intensity and polarization, depending on whichwas closer in the sky to the sources being observed. In thecase of 3C 286, we used J1407+2827 as the polarization leak-age calibrator; in the case of 3C 48 we used 3C 84. These areboth recommended as primary low polarization leakage calibratorsfor the VLA ( https://science.nrao.edu/facilities/vla/ http://simbad.u-strasbg.fr/simbad/ docs/manuals/obsguide/modes/pol ; see Table 7.2.3). Sincethese sources are blazars, their emission at these frequencies is dom-inated by the compact core and they can be treated to a good approx-imation as point sources. We therefore used the targets themselvesas phase calibrators, and performed pointing calibration scans at X-band throughout the observations on targets less than 20 ◦ from thesubsequent scans.The data reduction was performed in casa 5.3.0 and was based onthe standard VLA data reduction pipeline adapted for our observa-tions. As in the pipeline, we use the standard flux calibration sources3C 48 or 3C 286 to calculate delay and bandpass calibrations. De-parting from the pipeline, we then assume a point source model foreach of the target sources and self-calibrate them in amplitude andphase. We use the ‘fluxscale’ task to translate the amplitude solutionsinto a spectral fit for each target source; see Fig. 3 for an example fit).To check that our point source model accurately describes thesources, we subtract the model from the calibrated data and inspectthe residual maps and visibilities. In most cases the residuals are <
1% of the core flux density. In the few cases that the residuals arelarger, we iteratively improve the model by adding in the non-coreflux observed on the map, then solving again for the core flux densityin the updated core + non-core model. Maps and 𝑢𝑣 -plane residualsare shown for all sources in Appendix A (available as supplementarymaterial).We then follow EVLA Memo 201 (Hales 2017a) and Hales(2017b) for the polarization calibration strategy, in particular follow-ing strategy ‘C1’ with one unpolarized calibrator and one calibratorwith known polarization: • Set the full polarization model for the flux density calibrator us-ing fits to polarization fractions and angle as a function of frequencyfrom Perley & Butler (2013). If using 3C 48, which is resolved, flagbaselines where the model I flux is less than 90% of the total model I flux and treat as an unresolved source • Find cross-hand delay (‘kcross’) solutions using calibrator withknown polarization • Solve for polarization leakage per-channel (‘Df’) using the un-polarized calibrator if bright enough; otherwise (observations usingJ1407+2827 only) solve using a combination of the science targets(see below for more details) • Solve for the absolute polarization angle (‘Xf’) using calibratorwith known polarization.As for the total intensity measurements, we solve for the spectralenergy distributions of the sources in polarization in the uv -plane;this will be described in Section 3.2. We also check the image and 𝑢𝑣 -plane for significant residuals from the point-source model andfind none with the exception of the source PCCS1 G182.17+34.17,which displays significant polarization residuals at the positions ofthe intensity peaks seen in Fig. A42. For consistency with the othersources we quote the point-source fit for this source, however clearlywhen observed by a lower-resolution instrument such as QUIJOTEthis source would appear to have a different polarization.In the 2018 – 2019 observations, the polarization leakage cal-ibrator J1407+2827 had become too faint to successfully solve forthe per-channel leakage solutions. We therefore adapted our calibra-tion procedure as follows. We made a first-pass leakage calibrationwith J1407+2827 using per-spectral-window (‘D’) solutions. We thenused all science targets with polarization fractions <
5% (from thefirst-pass analysis), plus 3C286, to make a second-pass solution, asfollows. We used the first-pass solution to insert a model for thepolarized flux for all sources. We then repeated the polarization leak-age calibration using ‘Df+X’ mode to solve simultaneously for the
MNRAS , 1–17 (2021) Table 2.
Sources in the sample, positions and associations taken from SIMBAD, classifications taken from BZCAT (Massaro et al. 2015). Redshifts are fromBZCAT and references therein unless otherwise specified. The three resolved sources are excluded from the VLA sample.
PCCS1030 𝑧 RA Dec Name Class CommentG107.00-50.62 0.089 00:10:31.01 +10:58:29.5 5BZQ J0010+1058 FSRQ Mrk 1501 (Sy1)G115.13-64.85 0.21968 𝑎 𝑏 𝑐 𝑑 𝑒 𝑓 𝑎 Schmidt (1965), 𝑏 Trager et al. (2000), 𝑐 Popescu et al. (1996), 𝑑 Meisner & Romani (2010), 𝑒 Huchra, Vogeley & Geller (1999), 𝑓 de Vaucouleurs et al. (1991) leakage and polarization angle terms. We expect any errors in thefirst-pass solution to leave only a small polarized residual which willtend to cancel out between the different sources since the polariza-tion angles are uncorrelated. We then solve for the polarization angleusing 3C286 as usual.To demonstrate the validity of this method, in Fig. 1 we show acomparison between some sources calibrated using this method andwith the ‘D’ solutions using J1407+2827; it is clear that the resultsare the same but the calibration is more accurate. We also attempteda more standard calibration by merely substituting J1407+2827 withthe source with the lowest polarization fraction, as determined usingthe J1407+2827 ‘D’ solutions. We show an example calibration ofthis type in Fig. 1 as well; although the new leakage calibration source has a very low polarization fraction of 𝑝 = .
4% it isstill detected and using it as the leakage calibrator gives results thatare biased with respect to the J1407+2827 and combined solutions.Our ‘combined’ method is therefore the best way to calibrate thepolarization leakage given the limitations of the data. We apply thismethod also to the 2016 observation using J1407+2827 as the leakagecalibrator to improve the SNR of the leakage solutions.
MNRAS000
MNRAS000 , 1–17 (2021)
Y. Perrott et al.
Table 3.
Dates on which the sources were observed. ‘Num.’ refers to thenumber of target sources observed in each observation. ‘Beam size’ is theaverage synthesized beam major and minor full width at half maximum, inarcsec.Project Date Num. Flux VLA Beamcode calibrator config. size15A-083 10/03/2015 19 3C 48 B 0 . × . . × . . × . . × . . × . . × . 𝑎 . × . . × . 𝑏 . × . . × . 𝑎 Data affected by instrumental issues and not used 𝑏 Repeat of 27/11/2018 observation
Figure 1. 𝑞 and 𝑢 values recovered in each spectral window for two sources inthe 14/12/2018 observation, after full calibration using (i) ‘D’ solutions usingJ1407+2827 (blue dots); (ii) ‘Df+X’ solutions using all science targets with 𝑝 <
5% (black crosses); (iii) ‘Df’ solutions using PCCS1 G174.43+69.81with a low 𝑝 = .
4% (magenta triangles). In the 3C286 plot the black lineis the model. The ‘Df+X’ solutions clearly lead to a more accurate calibrationwhich is unbiased compared to the model and the (less accurate) J1407+2827calibration; the PCCS1 G174.43+69.81 calibration is biased even though thepolarization fraction is very small.MNRAS , 1–17 (2021) The expected calibration accuracy is 5% given that we are using the3-bit samplers and therefore cannot apply the switched power calib-ration . Given the variable nature of the sources, assessment of thecalibration accuracy requires simultaneous monitoring. We use twodatasets for this purpose; the Metsähovi 37 GHz blazar monitoringdata (Terasranta et al. 1992) and the Owens Valley Radio Obser-vatory (OVRO) blazar monitoring data at 15 GHz (Richards et al.2011); since many of our sources are frequently monitored by bothprogrammes we can interpolate the monitoring flux densities to theappropriate time and compare to our VLA measurements.In the case of Metsähovi, we find 48 semi-simultaneous (within 10days) measurements of 30 sources, including the calibration source3C 84. We correct our VLA measurements to 37 GHz using the fit-ted in-band spectral index and compare the measured flux dens-ities; this is shown in Fig. 2. In general we see good agreement.There is a slight overabundance of sources with significantly lowVLA flux density compared to Metsähovi which probably indicatesa resolution effect. The Metsähovi beam is 2.4 arcmin at 37 GHz.The VLA primary beam across our band has a FWHM of at most1.6 arcmin, and the shortest VLA baselines are ≈
10 k 𝜆 measuringstructure on ∼
20 arcsec scales, so both nearby sources outside of theVLA primary beam and flux resolved on scales between 20 arcsecand 2.4 arcmin would be included in the Metsähovi measurementsbut not detectable on our VLA maps. In addition, the Metsähovierrorbars are relatively large (3% systematic error for the brightestsources, and errors up to 20% for the faintest in the overlappingsample), and the lightcurves are often sparsely sampled (around 20day average cadence). The combination of resolution, large errorbarsand sparse sampling makes testing for the expected ∼
5% VLA calib-ration uncertainty difficult using the Metsähovi data. Fitting a straightline to 𝑆 VLA vs 𝑆 Metshovi using orthogonal least-squares regressiongives a slope and offset consistent with 1 and 0 within 2 𝜎 , where 𝜎 are the errorbars of the fit. We also note that the measurements ofthe brightest source, 3C 84, are consistent within 2 𝜎 (without addingadditional calibration errors to the VLA data) and conclude that theVLA and Metsähovi data are consistent within the limitations of thedata.In the case of OVRO, we find 83 semi-simultaneous measurementsof 43 sources. To compare the VLA and OVRO flux measurements,we calculate spectral indices between 15 and 34 GHz. Assuming thatthere is no spectral curvature between 15 and 40 GHz, this shouldbe equal to the spectral index over the VLA bands. Fig. 2 shows thatthere is indeed a very good correspondence between these quantities.We assume offsets from this relationship will be dominated by theVLA systematic calibration uncertainties in 𝛼 due to the muchshorter frequency lever arm and therefore use this comparison to testfor the value of this uncertainty. We do this by fitting a straight lineto the data (excluding the significant outlier PCCS1 G107.00-50.62)using orthogonal least-squares regression. We use the OVRO errorsas given in the database and use error propagation to calculate theerror in the interpolated quantity. We add an extra VLA calibrationerror in 𝛼 to the statistical error calculated from the fit over theVLA band, and increase this error until reduced 𝜒 = See https://science.nrao.edu/facilities/vla/docs/manuals/oss/performance/vla-samplers for an explanation of theexpected calibration accuracy. least-squares regression, with the nominal 5% flux density calibrationuncertainty and spectral index calibration uncertainty of 0.1 addedto the statistical errors in quadrature, gives a best-fit slope consistentwith one (1 . ± .
06) and some indication of an offset ( − . ± . 𝑝 𝜖 as the spurious on-axis fractional polarization producedby a small amount of polarization actually being present:2 ( 𝜎 𝑑 ) ≈ 𝑝 + 𝑁 𝑎 𝐴 𝑝 𝜖 ≈ 𝜎 𝑑 √︂ 𝜋 𝑁 𝑎 (1)(using equations 30 and 28 from Hales 2017b), where 𝜎 𝑑 is theerror in the leakage ‘d-term’ calibration accuracy, 𝑁 𝑎 is the num-ber of antennas, and 𝐴 is the full-array dual-polarization total in-tensity signal to noise of the calibrator within the single spectralchannel of interest. 3C 84 is known to approach 1% true polar-ization at 43 GHz (e.g. ); J1407+2827 has been shown tohave polarization of (cid:47) .
3% at Q-band (Liu et al. 2018) and < . 𝑝 true = 1% for 3C 84and 0.3% for J1407+2827; in the case of the 2018 – 2019 solutionswe conservatively adopt 𝑝 true = 1% as the maximum residual withrespect to the first-pass fit. This gives per-channel 𝑝 𝜖 values < . 𝛷 𝜖 using their publicly available softwarepolcalsims (Hales 2017c) which uses Monte Carlo simulations topredict position angle error as a function of signal-to-noise of thepolarized calibrator. The per-channel estimates range from 0.15 ◦ to1.8 ◦ in the case of 3C 48 and 0.05 ◦ to 0.8 ◦ in the case of 3C 286, againwith the higher errors corresponding to higher frequencies where thesignal-to-noise of the sources is lower. https://github.com/chrishales/polcalsims MNRAS000
3% at Q-band (Liu et al. 2018) and < . 𝑝 true = 1% for 3C 84and 0.3% for J1407+2827; in the case of the 2018 – 2019 solutionswe conservatively adopt 𝑝 true = 1% as the maximum residual withrespect to the first-pass fit. This gives per-channel 𝑝 𝜖 values < . 𝛷 𝜖 using their publicly available softwarepolcalsims (Hales 2017c) which uses Monte Carlo simulations topredict position angle error as a function of signal-to-noise of thepolarized calibrator. The per-channel estimates range from 0.15 ◦ to1.8 ◦ in the case of 3C 48 and 0.05 ◦ to 0.8 ◦ in the case of 3C 286, againwith the higher errors corresponding to higher frequencies where thesignal-to-noise of the sources is lower. https://github.com/chrishales/polcalsims MNRAS000 , 1–17 (2021)
Y. Perrott et al.
Figure 2.
Simultaneous 37 GHz total intensity flux densities measured byMetsähovi and the VLA (top), and simultaneous 15 – 34 GHz spectral indices, 𝛼 calculated using the OVRO blazar monitoring data, compared to spectralindices across the VLA bands, 𝛼 (bottom). The points are coloured byVLA observing epoch. In both cases, the solid black line shows a one-to-onecorrespondence while the black dashed line shows the best straight-line fitusing orthogonal least-squares regression. In the Metsähovi plot, no extrasystematic errors have been added to the VLA datapoints. In the OVROplot, conservative VLA calibration uncertainties of 5% in the reference fluxdensity and 0.1 in the spectral index are included in the errorbars plottedand used for the least-squares fit. The points marked with squares in bothplots are PCCS1 G107.00-50.62 which may have a turnover in its spectrumconsidering the 37 GHz flux agrees well but the spectral index is discrepantwith OVRO. These points are excluded from the fit to the OVRO data. We assume both of these systematic errors are completely correl-ated between channels in a spectral window and so use the averageper-channel estimate in each spectral window as the overall system-atic error for the spectral window estimate. We use error propagationto translate these systematic errors in 𝑝 and 𝛷 in errors in 𝑞 and 𝑢 : Δ 𝑞 = 𝑝 𝜖 cos ( 𝛷 ) + 𝑢 𝛷 𝜖 Δ 𝑢 = 𝑝 𝜖 sin ( 𝛷 ) + 𝑞 𝛷 𝜖 . (2) 𝑞𝑢 -fitting To fit the polarization properties, we take a vector average over timeand baseline of the calibrated RL and LR cross-hand visibility datain each spectral window and solve for Q and U using 𝑉 𝑅𝐿 = Q + 𝑖 U 𝑉 𝐿𝑅 = Q − 𝑖 U . (3)We estimate the uncertainty in 𝑞 and 𝑢 using the standard deviationin the data and adding in quadrature the systematic errors estimatedas in Section 3.1 (the error in I is negligible compared to the errorin Q and U ).We follow a similar 𝑞𝑢 -fitting procedure to O’Sullivan et al. (2012)in that we fit jointly to 𝑞 and 𝑢 , minimizing the total 𝜒 (rather thanfitting to 𝑝 and 𝛷 separately). This guarantees consistency betweenthe derived 𝑝 and 𝛷 properties. Our frequency band does not containenough information to usefully constrain physical parameters, sounlike O’Sullivan et al. (2012), rather than fitting a physical modelwe simply fit a power-law to the polarization fraction as a functionof frequency, and a rotation measure (RM) law, 𝛷 = 𝛷 + RM 𝜆 lawto the polarization angle as a function of wavelength. This gives: 𝑞 = 𝑝 cos ( 𝛷 ) = 𝑝 (cid:18) 𝜈𝜈 (cid:19) 𝛼 𝑃 cos (cid:104) (cid:16) 𝛷 + RM 𝜆 (cid:17)(cid:105) 𝑢 = 𝑝 sin ( 𝛷 ) = 𝑝 (cid:18) 𝜈𝜈 (cid:19) 𝛼 𝑃 sin (cid:104) (cid:16) 𝛷 + RM 𝜆 (cid:17)(cid:105) (4)where our reference frequency 𝜈 is again 34 GHz. We use thescipy.optimize.leastsq least-squares fitting algorithm to find thebest fitting 𝑝 , 𝛼 𝑃 , 𝛷 and RM and estimate their associated errors( Δ 𝑝 , Δ 𝛼 𝑃 , Δ 𝛷 and Δ RM).We have added the systematic errors to the per-spectral window 𝑞 and 𝑢 error estimates before performing the fits to account forsmall differences in systematic error levels across the frequency band.However, in general we expect the errors to be correlated across thefrequency band and therefore the parameter error estimates (producedassuming uncorrelated errors) will be underestimated. To account forthis, we rescale the 𝑞 and 𝑢 errors so that reduced 𝜒 = 𝑝 𝜖 (to Δ 𝑝 ) and 𝛷 𝜖 (to Δ 𝛷 ) to produce the final error.The 𝑞𝑢 -fitting procedure using our simple model provides reason-able fits in all cases and we do not see systematic trends indicating amore sophisticated model would be warranted; some sample fits areshown in Fig. 3. We define as undetected in polarization those sourceswith fitted 𝑝 / Δ 𝑝 <
3. We note that there are no 𝑛𝜋 ambiguities inthe fitted RMs, since our frequency channels are closely spaced. Wenote also that this method circumvents the issue of bias in polarizationfraction (i.e. 𝑝 meas = √︃ ( 𝑞 + 𝜎 𝑞 ) + ( 𝑢 + 𝜎 𝑢 ) > 𝑝 true = √︁ 𝑞 + 𝑢 )since we are fitting a spectrum in 𝑞 and 𝑢 across a large number offrequency measurements, each of which can be displaced in eitherthe positive or negative direction. This is demonstrated in Fig. 4where we show the fit to J1407+2807, the low-SNR source withnegligible polarization. The 𝑞𝑢 -fitting method clearly produces anunbiased result whereas fitting directly to 𝑝 = √︁ 𝑞 + 𝑢 produces abiased result. In Table 4 we show a sample section of our results table, containingflux density and spectral index estimates in total intensity and polar-ization fraction, as well as polarization angles extrapolated to 𝜆 = MNRAS , 1–17 (2021) Figure 3.
Sample fits to the 𝑢𝑣 -plane flux densities: left shows total intensity, centre polarization fraction, and right polarization angle as a function of frequency.The source is PCCS1 G147.84-44.04, observed on 10/03/2015, which shows the most extreme rotation measure in the sample. In all cases, each point representsa spectral window. In the total intensity plot, the errorbars are derived from the scatter within the frequency channels in the spectral window. In the polarizationplots, the errorbars represent the scatter within the frequency channels added in quadrature with the systematic errors estimated as described in Section 3.1. Theerrorbars do not take into account correlations between frequency channels and/or baselines and are not independent. Figure 4.
Sample fits to the polarization spectrum of J1407+2827, a low-SNR source with negligible polarization. From left to right we show 𝑞 , 𝑢 and 𝑝 . Asin Fig. 3, the errorbars are the total error including statistical and systematic contributions. The solid black lines show fits to the datapoints using the 𝑞𝑢 -fittingmethod, while the dashed black line in the 𝑝 plot shows the bias that would result from fitting directly to 𝑝 = √︁ 𝑞 + 𝑢 . Table 4.
A sample section of the results table. Each source has a separate entry for each epoch showing flux density at the central frequency of 34 GHz andspectral index in total intensity; polarization fraction (not percent) at 34 GHz and spectral index, extrapolated polarization position angle and RM, along withassociated errors. The errors in the table include the systematic errors estimated as in Section 3.1. The full table is available as an online-only supplement andat CDS.PCCS1 030 Epoch 𝐼 Δ 𝐼 𝛼 Δ 𝛼 𝑝 Δ 𝑝 𝛼 𝑃 Δ 𝛼 𝑃 𝛷 Δ 𝛷 RM Δ RMMJD Jy Jy deg deg rad m − rad m − G178.26+33.40 57091.06 0.6300 0.0315 -0.360 0.101 0.02972 0.00186 0.080 0.073 108.32 1.18 391 223G182.17+34.17 57091.06 0.7121 0.0356 -0.426 0.101 0.06189 0.00191 0.323 0.079 44.50 0.68 -789 93G200.04+31.88 57091.06 0.7326 0.0366 0.058 0.101 0.01624 0.00195 -0.358 0.360 31.72 1.70 3003 330G143.53+34.42 57091.06 2.8432 0.1422 0.171 0.101 0.01885 0.00133 -1.187 0.139 92.04 0.86 961 162G206.82+35.81 57091.06 5.5528 0.2777 -0.156 0.101 0.02064 0.00186 0.551 0.101 144.15 0.93 160 169G175.72+44.81 57091.06 2.3493 0.1175 0.022 0.100 0.01957 0.00186 -0.793 0.105 76.35 1.53 -1007 301G152.23+41.00 57091.06 0.8449 0.0423 -0.245 0.101 0.04056 0.00192 0.226 0.122 137.42 0.87 -1529 150G198.82+44.43 57091.06 0.9126 0.0456 0.238 0.100 <0.00558 0.00186 - - - - - -G183.71+46.17 57091.06 7.3728 0.3687 -0.301 0.101 0.01772 0.00204 0.096 0.460 137.41 1.35 -195 240G181.02+50.29 57091.06 0.5709 0.0286 -0.455 0.100 0.01440 0.00185 -0.899 0.114 176.76 1.53 4664 300
In Fig. 5 we show the distributions of central flux densities (at34 GHz) and spectral indices, including each epoch of measurementfor each source separately. These show that most (87%) of the sourcesare flat-spectrum with fitted spectral indices 𝛼 > − .
5. This is as ex-pected given the high-frequency selection criterion and in agreementwith the BZCAT classifications. We also see that although the sourceswere selected to be brighter than 1 Jy at 30 GHz, in fact only slightlymore than half (61%) have 34 GHz flux densities greater than 1 Jy; this shows the variability of these sources and highlights the difficultyin selecting a complete sample.
Since we only have two epochs of observation for each source, wefirst use the Metsähovi lightcurves to check whether our observationscan fairly test the variability of our sample. For each of the 31 (outof 54) sources with Metsähovi measurements spanning ≈ MNRAS000
Since we only have two epochs of observation for each source, wefirst use the Metsähovi lightcurves to check whether our observationscan fairly test the variability of our sample. For each of the 31 (outof 54) sources with Metsähovi measurements spanning ≈ MNRAS000 , 1–17 (2021) Y. Perrott et al.
Figure 5.
Distributions of total intensity flux densities at 34 GHz (left) andspectral indices (right) across the VLA observation band. The solid line inthe flux density histogram indicates the selection cut at 1 Jy (at 30 GHz, butwe do not correct for the small frequency shift) and the dashed line in thespectral index histogram indicates the traditional division into steep- andflat-spectrum sources.
Figure 6.
Stacked distributions of fractional variations in I with respect tothe mean for sources with Metsähovi lightcurves. The blue, filled histogramshows the total lightcurve while the red, unfilled histogram shows the dis-tribution when sampling each lightcurve only at two points, mimicking theVLA sampling. urements of (I− ¯ I)/ ¯ I where I are the individual flux measurementsand ¯ I is the mean flux density for each source over the whole light-curve. Then, we replicate our VLA sampling by interpolating thelight curves at points corresponding to the VLA measurements foreach source, which produces 55 individual measurements. We cal-culate the fractional variation for these points as above, with respectto the mean over the whole lightcurve. Fig. 6 shows the resulting twohistograms. They are of course not identical but show very similardistributions, and the standard deviation of each set of measurementsis nearly identical, 𝜎 = .
30 and 0 .
29. Restricting the sample tosources for which we have two VLA measurements (24 sources)does not change the distribution or the standard deviation. We there-fore can be confident that, although we only have two epochs ofmeasurement, our sample of sources is large enough to ensure thatwe are probing the variability statistics of the overall sample fairly.We note that there is a small asymmetry to larger positive values ofthe fractional variation apparent in both histograms, which can beexplained by the fact that most sources spend a large fraction of theirtime near a ‘baseline’ flux density value, with occasional excursionsto very high, flaring states. This results in the mean over a long periodof time being closer to the baseline, and the flares appearing as large,positive fractional deviations.Now considering the VLA measurements, in Fig. 7 we show his- tograms of the variability in total intensity flux density and spectralindex for the 41 sources for which we have two measurements. Sincewe only have two measurements of each source, we define frac-tional variation in flux density as (I , − I , )/I where I isthe weighted mean of the two flux densities, and spectral index vari-ation as 𝛼 − 𝛼 . We attempt to quantify the variability as follows.We have no reason to expect the flux and spectral index to changemore in one direction than the other, given that the sources wereobserved at different epochs and all will have different variabilitytimescales; we do not expect to see the bias toward positive frac-tional deviations in I found in the Metsähovi measurements abovesince our fractional deviations are with respect to the weighted meanof two measurements only, rather than the whole lightcurve. There-fore we symmetrize the distributions by adding the reflection ofeach point to the dataset, i.e. our symmetrized datasets consist of (cid:104) (I , − I , )/I , (I , − I , )/I (cid:105) and [ 𝛼 − 𝛼 , 𝛼 − 𝛼 ] .To fit a model probability distribution to each of the variabilitydistributions, we use a Kolmogorov–Smirnov (KS) test. For a givenmodel for the distribution (either Cauchy or Gaussian) and appropri-ate set of parameters for the model, we calculate the KS ‘D’-statisticbetween the symmetrized data and the model which gives a measureof the maximum difference between the cumulative probability dis-tributions of the two. We then vary the parameters in the model untilthe smallest ‘D’-statistic is reached, i.e. finding the set of parametersthat best match that model to the data. This method has the advantagethat the fit to the model does not depend on how the data are binned.In the case of the flux density variations, we find a Cauchy distri-bution centred at 0 fits the distribution well, i.e.PDF (cid:32) Δ 𝐼 I (cid:33) = 𝜋𝛾 + (cid:32) Δ 𝐼 I 𝛾 (cid:33) − (5)with the best-fit 𝛾 = .
22. The KS test ‘D’-statistic between the sym-metrized (non-symmetrized) data and the fitted Cauchy distributionis 0.04 (0.08) and the 𝑝 -value is 1.00 (0.96) indicating a good fiteven to the non-symmetrized data. We show both the symmetrizedand non-symmetrized data in Fig. 7, along with the best-fit Cauchydistribution.In the case of the 𝛼 variation, we find a Gaussian distribution tobe a better fit than a Cauchy distribution, although clearly it doesnot fit the data particularly well. We find a best-fit standard deviation 𝜎 = .
26 with mean 𝜇 fixed to 0. The KS test ‘D’ values between thesymmetrized (non-symmetrized) data and the fitted distribution are0.04 (0.18) and the 𝑝 -values are 1.00 (0.13). The measurements areclearly limited by the relatively large spectral index calibration error,and more accurate measurements with a longer frequency lever armwould be required to investigate further the spectral index variability.Fig. 7 also shows the correlation between variability in (non-symmetrized) total intensity flux density and spectral index, which isweak with Pearson 𝑅 = . 𝑝 -value=0.33, indicating that there isnot enough information in our two-point estimates to investigate thephysical mechanisms underlying the changes. In general one wouldexpect a change in flux density to be accompanied by a change inspectral index and vice versa. There are, however, some cases wherea highly significant variation in spectral index occurs although theflux density is nearly constant; we investigate one of these cases inSection 5.1. A change in spectral index is therefore also a good in-dicator of variability. In total 59% of the sources have a change ineither flux density or spectral index of > 𝜎 (using the systematicerrors as calibrated in Section 3.1). MNRAS , 1–17 (2021) Figure 7.
Histograms of fractional variation in total intensity flux density andspectral index variation (histograms), and correlation between these quantitiescolour-coded by epoch pair. Along with the histograms of the data (solidblue), we also plot the ‘symmetrized’ histograms (red outlines; normalizedby a factor of two for comparison with the non-symmetrized histograms) asdescribed in the text. The black lines show Cauchy (Gaussian) fits to thesymmetrized flux density (spectral index) variation distributions; see text formore details.
We detect 85% of the sources in polarization, considering all ob-serving epochs together. In Fig. 8 we show the distributions of po-larization fractions, spectral indices, RM and intrinsic polarizationangles fit to the polarization fraction spectra as well as correlationsbetween the four quantities. We see fairly low polarizations in mostcases, with a median of 2.2%; this is consistent with other studiesat similar frequencies (e.g. Sajina et al. 2011, Battye et al. 2011)and indicates that we are seeing mostly the flat-spectrum core ratherthan the steep-spectrum lobes, which are more highly polarized. Wealso show a log-normal fit to the distribution, made by minimizingthe KS ‘D’-statistic between the data and a log-normal distribution.This gives fitted parameters ( 𝜇, 𝜎 ) = ( . , . ) which are in goodagreement with the fits to flat-spectrum source data at various fre-quencies from Puglisi et al. (2018) (e.g. their Fig. 3). We see a fairlybroad distribution of spectral indices in polarization (although thesedo have large errors due to the small change observed over the band),with ≈
26% of sources having rising spectra at > 𝜎 significanceindicating depolarization, and ≈
15% having falling spectra at > 𝜎 significance indicating repolarization (we note that repolarizationcan occur as a consequence of physical models more complicatedthan a single Faraday screen; see e.g. O’Sullivan et al. (2012) for ex-ample models and data). We do not attempt fits to the 𝛼 𝑃 and | RM | distributions due to the large and non-uniform errors. As expected,the intrinsic polarization angles are uniformly distributed.We note that the rotation measures in Fig. 8 are observed RMs withno correction for Galactic RM or conversion to the AGN rest framesince we are interested in the observed properties in the context of the QUIJOTE analysis. We detect an RM at > 𝜎 significance in 65% ofthe polarization detections and find a median | RM | ≈ − .Previous VLBI studies have found similarly high AGN core RMsat centimetre wavelengths (e.g. Zavala & Taylor (2004) find coreRMs up to ≈ − between 8 and 15 GHz;Hovatta et al. (2012) find core RMs up to ≈ > − between 8 and 15 GHz) whilethe jet RMs tend to be lower. This lends additional support to theidea that our 30 – 40 GHz observations are mostly probing the core.Hovatta et al. (2019) find that RM increases as a function of frequencyin 3C 273 with RM ∝ 𝜈 in agreement with models for a sheathsurrounding a conically expanding flow; if this effect occurs for themajority of sources it would explain our somewhat higher median.Indeed Algaba, Gabuzda, & Smith (2011) find a median RM of ≈ − between 15 and 43 GHz, although their values maysuffer from 𝑛𝜋 ambiguities. We see one extreme outlier in RM whichwe describe in more detail in Section 5.2.The only significant correlation between the polarization prop-erties is between the 𝛼 𝑃 and | 𝑅𝑀 | measurements, with Pearson 𝑅 = . 𝑝 -value = . 𝛼 𝑃 . However, the sources with signi-ficantly negative 𝛼 𝑃 also seem to show a high RM, suggesting thesituation is more complicated than a single Faraday screen and thepolarization angle spectra just happen to follow a 𝜆 law over the rel-atively small range in wavelength. Given the uncertainties introducedby averaging over our relatively large beam (compared to VLBI stud-ies which resolve the source structure) we do not investigate thisfurther. In Fig. 9 we show the distributions of and correlations between thevariation in polarization fraction, spectral index, RM and intrinsicpolarization angle for the sample. There are 34 sources with morethan one detection; of these 65% are variable in 𝑝 and all arevariable in one or more of these quantities at 3 𝜎 (where 𝜎 includessystematic and statistical contributions as described in Section 3.2).The notable outlier in polarization fraction is PCCS1 G145.78+43.13which we investigate further in Section 5.3. Aside from this extremeoutlier, the distribution of changes in polarization fraction is relativelyGaussian; a fit to the symmetrized data excluding the outlier givesa best-fit 𝜎 = . 𝑝 -value of 1.00 (0.96) indicating goodagreement.The changes in polarization fraction spectral index and RM areclustered around 0, with the exception of PCCS1 G145.78+43.13.It is more difficult to assess the Gaussianity of these distributionsgiven the large variations in measurement errors. Most of the out-liers do have large errors. Of the two highly significant outliers, oneis PCCS1 G145.78+43.13, as already mentioned and OVRO mon-itoring of the second, PCCS1 G129.09-13.46 shows that the twoobservations similarly took place during a quiescent period and aflare so these large changes are likely genuine.We see a significant correlation between Δ 𝛼 𝑃 and Δ | RM | withPearson 𝑅 = − .
42 and 𝑝 -value=0.01. This is in line with the cor-relation between 𝛼 𝑃 and | RM | discussed in Section 4.2; a sourceundergoing a higher degree of depolarization has a higher RM. Re-moving PCCS1 G145.78+43.13 strengthens the correlation slightlyto 𝑅 = − . 𝑝 -value=0.006. We also see a highly significant correl-ation between Δ 𝑝 and Δ | RM | with 𝑅 = − . 𝑝 -value = × − , MNRAS000
42 and 𝑝 -value=0.01. This is in line with the cor-relation between 𝛼 𝑃 and | RM | discussed in Section 4.2; a sourceundergoing a higher degree of depolarization has a higher RM. Re-moving PCCS1 G145.78+43.13 strengthens the correlation slightlyto 𝑅 = − . 𝑝 -value=0.006. We also see a highly significant correl-ation between Δ 𝑝 and Δ | RM | with 𝑅 = − . 𝑝 -value = × − , MNRAS000 , 1–17 (2021) Y. Perrott et al.
Figure 8.
Histograms showing the distribution of polarization fractions (in percent), polarization spectral indices, RM (in rad m − ) and intrinsic polarizationangles (in degrees) fit across the VLA bands, and correlations between the four quantities. The histogram of polarization fractions also shows a log-normal fitto the distribution and the median RM is indicated with a vertical line on the histogram. Points are colour-coded by observing epoch as in Fig. 2. however this correlation is largely driven by PCCS1 G145.78+43.13and removing this outlier decreases the correlation strength to 𝑅 = − . 𝑝 -value=0.13. This correlation is also consistent withthe idea that a higher RM implies a higher degree of depolarization,so increase in rotation measure correlates with decrease in polariza-tion fraction. In Fig. 10 we show the correlation between total intensity parametersand polarization parameters. All appear uncorrelated. Of particularimportance for predicting the contamination to B-mode analysis isthe correlation between total intensity flux density and polarizationfraction; i.e. if simulating point source contamination, can randomvalues be drawn independently from the total intensity source countand polarization fraction distribution, or is there a correlation? Al-though the Pearson 𝑅 value for 𝑝 vs I is low ( 𝑅 = − . 𝑝 = . I values. We test this idea by dividing the sources into three bins in I and calculating median polarization fractions in each bin. Wesee a slight positive correlation between flux density and medianpolarization fraction. To test the statistical significance of the trend,we calculate KS test statistics between the polarization fractions ineach bin and the overall fitted log-normal distribution. The resultsfor the bins are summarised in Table 5; we find that for the lower andmiddle flux density bins, the 𝑝 -values are ≈ ≈ MNRAS , 1–17 (2021) Figure 9.
Histograms showing the distribution of variation in polarization fraction (in percent), polarization spectral index, RM (in rad m − ) and intrinsicpolarization angle (in degrees) fit across the VLA bands, and correlations between the four quantities. Points are colour-coded by observing epoch pair as inFig. 7. The histogram of polarization fraction (in percent) also shows a normal distribution centred at 0 with standard deviation equal to the standard deviationmeasured from the observed distribution, excluding the large outlier. and polarization fraction. Here we see no evidence for correlation( 𝑅 = − . 𝑝 = .
85) and no evidence for a different distribution ofpolarization fractions when we divide into two bins at 𝛼 = − . 𝑝 -values are also reported in Table 5.In Fig. 11 we show correlations between the variations in total in-tensity flux density and spectral index with variation of the polariza-tion parameters. Here too we see very little correlation, as evidencedby the very low Pearson 𝑅 -coefficients shown on the plots. Thereseems to be a slight correlation between Δ 𝛼 and Δ 𝛼 𝑃 ( 𝑅 = − . 𝑝 = .
09) which may indicate that sources coming down from a flare(becoming more optically thin) are also becoming less depolarized( 𝛼 𝑃 becoming less positive), and vice versa. Here we investigate in more detail some of the interesting results forindividual sources.
Table 5.
Median polarization fraction (in percent) and KS test 𝑝 -values forpolarization percentages in different bins in total intensity flux density (top)and spectral index (bottom), compared to the log-normal fit for the overalldistribution. Flux densities are in Jy and 𝑛 is the number of sources in thebin. The top row is the overall distribution. I min I max 𝑛 𝑝 median 𝑝 -value0.23 18.4 82 2.06 0.950.23 0.60 9 1.44 0.0600.60 2.0 52 2.29 0.252.0 18.4 21 2.03 0.59 𝛼 min 𝛼 max 𝑛 𝑝 median 𝑝 -value-1.2 0.7 82 2.06 0.95-1.2 -0.5 12 2.17 0.80-0.5 0.7 70 2.06 0.94MNRAS000
Median polarization fraction (in percent) and KS test 𝑝 -values forpolarization percentages in different bins in total intensity flux density (top)and spectral index (bottom), compared to the log-normal fit for the overalldistribution. Flux densities are in Jy and 𝑛 is the number of sources in thebin. The top row is the overall distribution. I min I max 𝑛 𝑝 median 𝑝 -value0.23 18.4 82 2.06 0.950.23 0.60 9 1.44 0.0600.60 2.0 52 2.29 0.252.0 18.4 21 2.03 0.59 𝛼 min 𝛼 max 𝑛 𝑝 median 𝑝 -value-1.2 0.7 82 2.06 0.95-1.2 -0.5 12 2.17 0.80-0.5 0.7 70 2.06 0.94MNRAS000 , 1–17 (2021) Y. Perrott et al.
Figure 10.
Correlations between total intensity flux density (left column) andspectral index (right column) with polarization parameters. In the 𝑝 plots,the vertical lines show a division into flux density/spectral index bins, and theblack stars show the median polarization fraction within the bins. PCCS1 G156.86-39.13 is one of the cases where Δ 𝛼 is large yet itsflux density is nearly constant. It has coverage from OVRO and also at43 GHz from the VLBA-BU Blazar Monitoring Project . The VLAdata-points are consistently above the VLBA lightcurve since theVLBA resolves out some of the flux from the source, however we cansee that the two VLA epochs of observing happen to catch the sourceon either side of a peak so that the flux density is approximately thesame. Comparing the OVRO and VLBA light-curves we can see thatthe 43 GHz flux density decreases more quickly than the 15 GHz fluxdensity for both this peak and the earlier peak, indicating a changein optical depth. This is consistent with the VLA spectral indiceschanging from slightly positive (optically thick) to slightly negative(optically thin) between the two epochs. Figure 11.
Correlations between total intensity flux density variation (leftcolumn) and spectral index variation (right column) with variability in thepolarization parameters. Points are coloured by observing epoch pair as inFig. 7.
Figure 12.
Light curves for PCCS1 G156.86-39.13 from OVRO at 15 GHz(black errorbars); the VLBA at 43 GHz (red points) and the VLA observations(cyan lines from 28 GHz to 40 GHz, with triangles at the 40 GHz end). Seetext for more detail.MNRAS , 1–17 (2021) PCCS1 G147.84-44.04 or 4C 15.05 is an extreme outlier in RM,having RM ≈ − − − in the two observa-tion epochs; the polarization fits for the first epoch are shownin Fig. 3. At the redshift 𝑧 = .
833 of the source our observedfrequency 𝜈 obs =
34 GHz corresponds to an emitted frequency 𝜈 em = ( + 𝑧 ) 𝜈 obs =
62 GHz and the observed RMs correspondto an intrinsic RM of ( + 𝑧 ) RM obs ≈ × rad m − in the sourcerest-frame. To our knowledge this source has not been identified ashaving a particularly high RM in any other work. Zavala & Taylor(2004) and Hovatta et al. (2012) both fail to measure RMs for thissource which may be a consequence of the rapid depolarizationdown to their lower frequencies (see Fig. 3). The intrinsic RM islarger than the intrinsic RM measured by Hovatta et al. (2019) for3C 273 of ≈ . × rad m − at ≈
62 GHz (extrapolating the meas-urement of 6 . × rad m − at observed wavelength ≈
234 GHzusing RM int ∝ 𝜈 ) and ≈ × rad m − found for 3C 84 (extrapol-ating the measurement of 8 . × rad m − at observed wavelength ≈
230 GHz) from Plambeck et al. (2014). It may even approach thecurrent largest-known RM observed for the lensed quasar PKS 1830-211 Martí-Vidal et al. (2015), which is 10 rad m − at 𝜈 em = × rad m − at 𝜈 em =
62 GHz if theRM int ∝ 𝜈 law holds over such a wide range in frequency. PCCS1 G145.78+43.13 is a notable outlier in polarization fractionvariation ( > 𝜎 ). The total intensity and polarization fraction (inpercent) light curves for this source are shown in Fig. 13. It under-went an exceptionally high optical flare in 2015 (MJD = 57067),coinciding with the emergence of a new knot detected by the VLBA-BU-Blazar Monitoring Project (MAGIC Collaboration et al. 2018).Our first observing epoch of this source happened to coincide withthe optical flare, and the polarization fraction measured at this epochof 13.6% agrees with the VLBA measurement of ≈ ≈ + . + .
62; OVRO and VLBA monitoring data confirm that the sourceis undergoing a radio flare. The polarization fraction has decreasedto 0.7%, in line with the VLBA polarization values; the polariza-tion fraction spectral index has steepened significantly from + .
06 to − .
92; and the RM has significantly changed from −
440 to − − . These changes could be attributed to the integrated fluxdensity during the flare state containing a significant contributionfrom the emerging knot of plasma. As shown in, e.g. Battye et al. (2011) and Puglisi et al. (2018),polarized extra-galactic sources produce significant contaminationto the B-mode power spectrum. The strong variability of the sourcesthat dominate at the QUIJOTE frequencies which will be used forcosmological analysis means that care must be taken to correctlyaccount for their presence.Following the methodology described in Tucci & Toffolatti (2012),we now estimate the contribution of unresolved polarized radiosources to the angular power spectra at 30 GHz, and discuss the im-plications for that frequency channel of the QUIJOTE experiment.For a Poisson distribution of point sources with flux densities below
Figure 13.
Light curves for PCCS1 G145.78+43.13 from OVRO at 15 GHz(black errorbars); the VLBA at 43 GHz (red points) and the VLA observations(cyan lines from 28 GHz to 40 GHz, with triangles at the 40 GHz end). Thetop shows total intensity while the bottom shows polarization fraction (inpercent). The vertical black line shows the date of the optical flare. See textfor more detail. a certain cut-off value 𝑆 C , the contribution to the B-mode angularpower spectrum can be estimated as: 𝐶 𝐵𝐵ℓ = (cid:18) 𝑑𝐵𝑑𝑇 (cid:19) − (cid:104) 𝛱 (cid:105) ∫ 𝑆 C 𝑛 ( 𝑆 ) 𝑆 𝑑𝑆 (6)where 𝑛 ( 𝑆 ) is the differential number of sources per steradian, 𝑑𝐵 / 𝑑𝑇 is the conversion factor from brightness to temperature, and 𝛱 cor-responds to the fractional polarization.The model for the differential source counts is taken from de Zottiet al. (2005, 2010) . The average value (cid:104) 𝛱 (cid:105) can be computed usingthe fitted log-normal distribution function in Section 4.2, using theequations given in Battye et al. (2011) and Puglisi et al. (2018). Inour case, we have √︁ (cid:104) 𝛱 (cid:105) ≈ . 𝐶 𝐵𝐵ℓ ≈ . × − 𝜇 K for 𝑆 C =
20 Jy,and 𝐶 𝐵𝐵ℓ ≈ . × − 𝜇 K for 𝑆 C = 𝑟 = .
2, 0 . .
05. As an indication, wealso include the prediction for 𝑆 C = .
3, 0 . . ℓ =
80 (i.e. close to the recombinationbump of the expected cosmological signal), we find that the sourcecontribution is reduced from 0 . 𝜇 K (for 𝑆 C =
20 Jy) to 0 . 𝜇 K Available online http://w1.ira.inaf.it/rstools/srccnt/srccnt_tables.html . MNRAS000
20 Jy) to 0 . 𝜇 K Available online http://w1.ira.inaf.it/rstools/srccnt/srccnt_tables.html . MNRAS000 , 1–17 (2021) Y. Perrott et al.
Figure 14.
Predicted primordial E-mode signal (upper black line), lensingB-mode signal (lower black line), primordial B-mode signal depending ontensor-to-scalar ratio value (red and pink lines, labelled by 𝑟 -value), and pointsource contamination depending on total intensity cut-off value (coloureddashed and dotted lines, as labelled in caption), all at 30 GHz. It can beseen that when masking sources with total intensity > 𝑟 = . 𝑟 = . (for 𝑆 C = 𝑟 = .
2, consistently with the results presentedin Rubiño-Martín et al. (2012, see their Fig. 7) and Tucci et al.(2005). This highlights the importance of removing the contributionfrom these sources; since we find that both their total intensity andpolarization vary unpredictably, it is likely that they will need to bemasked in the analysis rather than subtracted directly, in the absenceof simultaneous polarimetric monitoring observations.Finally, the QUIJOTE experiment aims to reach a limit of 𝑟 = . 𝑟 = .
05 when combining the TGI results with thosefrom the Forty-GHz Instrument (Rubiño-Martín et al. 2012). As areference, the best current constraints on tensor-to-scalar ratio are 𝑟 < .
064 (Planck Collaboration X 2020). This study suggests thatwe will need to remove sources down to ≈
500 mJy in total intensity(see Fig. 14) to achieve the 𝑟 = . ≈
300 mJy at 30 GHz in order to reach the 𝑟 = .
05 goal.However more studies will need to be done on the source populationat these lower flux densities to verify this. We also emphasize thatthese results are frequency-dependent and other experiments oper-ating at different frequencies will be differently affected by pointsource contamination.
In order to assess the contamination of the QUIJOTE cosmologicalfields by polarized emission from radio sources, we observed 51sources, selected to be brighter than 1 Jy at 30 GHz, at 28 – 40 GHzwith the VLA. The sample is dominated by flat-spectrum radio quas-ars. For 41 of the sources, we have two epochs of observation which allows us to investigate the variability of the sources both in totalintensity and polarization. We find that:(i) Our in-band spectral indices agree well with simultaneouslymeasured 15 to 34 GHz spectral indices using OVRO monitoringdata, with some indication of spectral steepening at the higher band.(ii) The median polarization fraction of our sample at 34 GHzis 2.2%, with the largest being 14%; the distribution of polariza-tion fractions agrees well with polarization fraction distributions atvarious frequencies summarised in Puglisi et al. (2018).(iii) We find a median rotation measure of |RM| ≈ − ,with one extreme outlier (4C 15.05) having RM ≈ − − − in the two observation epochs. This may be amongstthe highest RMs measured up to now in quasar cores.(iv) We find hints of a correlation between the total intensity fluxdensity and the median polarization fraction, however a larger samplecomplete to a lower flux density level would be required to confirmthis. We find no correlation between the total intensity spectral indexof a source and its polarization fraction.(v) 59% of the sources are variable in total intensity, while all arevariable in polarization at 3 𝜎 level. Changes in polarization fractionare roughly Gaussian-distributed with 𝜎 = . ≈
300 mJy will need to be masked to reachthe QUIJOTE goal of 𝑟 < .
05. For general experiments aiming todetect inflationary B-modes the point source population will need tobe studied at the frequency of the experiment to determine the levelof masking required.
ACKNOWLEDGEMENTS
We thank an anonymous reviewer who substantially improved thepresentation of the paper. The
QUIJOTE
MNRAS , 1–17 (2021) DATA AVAILABILITY
Raw data used in this study are publicly available from the VLAdata archive ( https://archive.nrao.edu ) and can be accessedusing the project codes in Table 3. Calibrated data are available onreasonable request from the authors.
References
Algaba J. C., Gabuzda D. C., Smith P. S., 2012, MNRAS, 420, 542Algaba J. C., Gabuzda D. C., Smith P. S., 2011, MNRAS, 411, 85Battye R. A., Browne I. W. A., Peel M. W., Jackson N. J., Dickinson C., 2011,MNRAS, 413, 132Bonavera L., González-Nuevo J., Argüeso F., Toffolatti L., 2017, MNRAS,469, 2401de Vaucouleurs G., de Vaucouleurs A., Corwin H. G., Buta R. J., Paturel G.,Fouque P., 1991, rc3..bookde Zotti G., Ricci R., Mesa D., Silva L., Mazzotta P., Toffolatti L., González-Nuevo J., 2005, A&A, 431, 893de Zotti G., Massardi M., Negrello M., Wall J., 2010, A&ARv, 18, 1Galluzzi V., Massardi M., 2016, IJMPD, 25, 1640005Galluzzi V., et al., 2017, MNRAS, 465, 4085Galluzzi V., et al., 2018, MNRAS, 475, 1306Grainge K., et al., 2003, MNRAS, 341, L23Hales C. A., 2017, CASA Interferometric Pipeline Polarization Calibration& Imaging Requirement & Design Specifications, EVLA Memo Series,201; also ALMA Memo Series, 603Hales C. A., 2017, AJ, 154, 54Hales, C. A., 2017, Zenodo, http://doi.org/10.5281/zenodo.801336Hovatta T., et al., 2012, AJ, 144, 105Hovatta T., O’Sullivan S., Martí-Vidal I., Savolainen T., Tchekhovskoy A.,2019, A&A, 623, A111Huchra J. P., Vogeley M. S., Geller M. J., 1999, ApJS, 121, 287Huffenberger K. M., et al., 2015, ApJ, 806, 112Jackson N., Browne I. W. A., Battye R. A., Gabuzda D., Taylor A. C., 2010,MNRAS, 401, 1388Jones T. W., Rudnick L., Aller H. D., Aller M. F., Hodge P. E., Fiedler R. L.,1985, ApJ, 290, 627Jorstad S. G., et al., 2007, AJ, 134, 799Kravchenko E. V., Cotton W. D., Kovalev Y. Y., 2015, IAUS, 313, 128Liu H. B., Hasegawa Y., Ching T.-C., Lai S.-P., Hirano N., Rao R., 2018,A&A, 617, A3Lister M. L., 2001, ApJ, 562, 208López-Caniego M., Massardi M., González-Nuevo J., Lanz L., Herranz D.,De Zotti G., Sanz J. L., Argüeso F., 2009, ApJ, 705, 868MAGIC Collaboration, et al., 2018, A&A, 617, A30Martí-Vidal I., Muller S., Vlemmings W., Horellou C., Aalto S., 2015, Sci,348, 311Massaro E., Maselli A., Leto C., Marchegiani P., Perri M., Giommi P., Pir-anomonte S., 2015, Ap&SS, 357, 75Massardi M., et al., 2008, MNRAS, 384, 775Meisner A. M., Romani R. W., 2010, ApJ, 712, 14Orienti M., Dallacasa D., 2008, A&A, 479, 409O’Sullivan S. P., et al., 2012, MNRAS, 421, 3300Park J., et al., 2018, ApJ, 860, 112Perley R. A., Butler B. J., 2013, ApJS, 206, 16 Perley R. A., Butler B. J., 2017, ApJS, 230, 7Plambeck R. L., et al., 2014, ApJ, 797, 66Planck Collaboration XXVIII, 2014, A&A, 571, A28.Planck Collaboration XXVI, 2016, A&A, 594, A26Planck Collaboration XLV, 2016, A&A, 596, A106Planck Collaboration X, 2020, A&A, 641, A10.Popescu C. C., Hopp U., Hagen H. J., Elsaesser H., 1996, A&AS, 116, 43Puglisi G., et al., 2018, ApJ, 858, 85Richards J. L., et al., 2011, ApJS, 194, 29Rubiño-Martín J. A., et al., 2010, ASSP, 14, 127Rubiño-Martín J. A., et al., 2012, SPIE, 84442Y, SPIE.8444Rudnick L., et al., 1985, ApJS, 57, 693Sajina A., Partridge B., Evans T., Stefl S., Vechik N., Myers S., Dicker S.,Korngut P., 2011, ApJ, 732, 45Schmidt M., 1965, ApJ, 141, 1Terasranta H., et al., 1992, A&AS, 94, 121Tucci M., Martínez-González E., Vielva P., Delabrouille J., 2005, MNRAS,360, 935Tucci M., Toffolatti L., 2012, AdAst, 2012, 624987Trager S. C., Faber S. M., Worthey G., González J. J., 2000, AJ, 119, 1645Zavala R. T., Taylor G. B., 2004, ApJ, 612, 749This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000