© ESO 2023

1 Introduction

There are many categories of prominent radio, millimeter-wave, and submillimeter sources, including Galactic sources such as H II regions and supernovae remnants, and extragalactic systems such as quasars and dusty star-forming galaxies. The brightest of these sources at high frequencies, however, are typically either Galactic or active galactic nuclei (AGN). Both classes are considered here, though the bulk of our objects are extragalactic. AGN are the most powerful long-lasting sources in the Universe. Unlike much shorter timescale events such as black hole mergers or gamma-ray bursts, the active phase lasts long enough so that its properties and evolution can be studied in detail. Radio and submillimeter observations are central to understanding these systems. The radio emission of AGN is almost certainly related to a black hole at the center of the active galaxy and to the jets it engenders. These massive black holes in turn appear to be intimately related to the evolution of their host galaxies (Magorrian et al. 1998). AGN feedback in the form of radiation, winds, and relativistic jets is now invoked to explain the galaxy luminosity function and the M–σ relation (e.g., Croton et al. 2006; Fabian 2012; King & Pounds 2015). The total power in the jets is a critical parameter in estimating the magnitude of this feedback. The physics required to estimate the jet power, however, is poorly understood. Important clues are provided by studies of both polarization and variability of emission from the jets. Such observations can illuminate the launching and initial propagation of relativistic jets, as well as the presence of jet “sheaths” and surrounding winds that could contribute to feedback (Blandford et al. 2019). The optical depth of the material surrounding the central black hole increases with wavelength, however. Hence, millimeter and submillimeter observations are particularly useful: they allow us to study the properties of both the black hole and its surrounding accretion disk, as has been done for a few bright, relatively nearby, radio galaxies such as M87 (3C 274) (Junor et al. 1999; Hada et al. 2011) and 3C84 (NGC 1275) (Giovannini et al. 2018). In these objects, as well as more distant radio galaxies and blazars, the millimeter-submillimeter continuum spectrum, time variability, and linear polarization probe the acceleration and collimation zone of the jet.

The observed properties of AGN strongly depend on the orientation between the line of sight and the axis of the jets. Emission from the jet itself is dominant when the line of sight and jet axis align closely. This class of AGN, the blazars, includes BL Lac objects and flat-spectrum radio quasars (FSRQs). Doppler boosting causes the emission of the aligned jet to be strong and the counter-jet to be weak. A substantial fraction of the AGN detected by Planck and other cosmic microwave background (CMB) experiments are blazars (Planck Collaboration XIII 2011). That is not surprising, given the high apparent luminosity of blazars and their continuum spectra, which peak at millimeter to submillimeter wavelengths. The close alignment of jet and line of sight also allows for variability introduced by small changes in the angle between the line of sight and the jet (Begelman et al. 1984; Chen et al. 2013).

To maximize the science that can be extracted from radio and submillimeter observations of AGN, we ideally need high angular resolution, wide frequency coverage (extending to the submillimeter range) and repeated observations over time-scales of days to years including both total power and polarization measurements. Ground based interferometers, including the Atacama Large Millimeter Array (ALMA), meet the first of these needs. We show here that data derived as a byproduct of CMB experiments, properly analysed, can meet the last two. The sensitivity of current CMB experiments, including the Planck mission, is adequate to detect thousands of the brightest AGN. Most current CMB experiments now employ polarization-sensitive detectors and combine observations at several frequencies. Planck, for instance made polarized observations of the entire sky at seven frequencies from 30 to 353 GHz, some of them inaccessible from the ground. Finally, CMB searches generally cover large areas of the sky (the entire sky, in the case of Planck) repeatedly, providing observations of sources with cadences ranging from minutes to years.

To date, however, the potential contributions of CMB experiments and data sets have not been fully realized because the repeated observations of a given sky region are merely added together to produce static sky maps (Partridge et al. 2017). Combining repeated observations such as this makes perfect sense for the CMB and for Galactic foregrounds, both of which are static, but all information in the time domain is lost. We remedy that situation here by working with the time-ordered data (ToD) available from the multiyear Planck mission. The methods we employ can easily be extended to data coming from other CMB surveys. Here we provide a catalog of extragalactic sources measured in both total power and polarization by Planck, with frequency coverage from 30 to 353 GHz, and time-scales ranging from days to years.

We discuss the Planck mission, source selection, calibration and related topics in Sect. 2. In Sect. 3, we describe our technical approach to extracting flux density and polarization measurements from the Planck ToD, source by source and epoch by epoch. Section 4 describes the catalog we construct of these observations. Section 5 then treats the results and some consequences for models of AGN. We conclude in Sect. 6.

2 Planck data

All the data employed here were taken from archived time-ordered data (ToD) from the Planck mission as found in the Planck Legacy Archive¹. We first describe aspects of the mission relevant to our study, then describe how we selected bright sources. While most bright sources are AGN, we include some Galactic and thermal sources in our study. Galactic sources, such as the planetary nebula NGC 6572, are not expected to vary on short time scales, and hence provide a useful test of our pipeline.

2.1 The Planck mission

Planck was a European Space Agency (ESA) mission designed to map fluctuations in the CMB to the limit set by unavoidable cosmic variance down to angular scales of 5′ or so (Planck Collaboration I 2020). Two instruments shared a common focal plane: the low-frequency instrument (LFI), with HEMT radiometers operating at 30, 44 and 70 GHz, and the high-frequency instrument employing lower-noise bolometric detectors in six frequency bands from 100 to 857 GHz. The fractional bandwidth at each frequency ranged from 20% to 30%. Polarization was measured at all but the two highest-frequency bands (545 and 857 GHz). Table 1 lists some of the relevant properties of the Planck mission. While mapping the CMB, Planck detected thousands of compact sources, both Galactic and extra-galactic. Planck thus provides a valuable archive of observations of these sources, which we explore in this paper. These sources are cataloged, frequency by frequency, in several Planck catalogs; we employ only the latest, PCCS2 (Planck Collaboration XXVI 2016). Approximate completeness limits in this catalog are provided in Table 1. These limits give a rough indication of the sensitivity of Planck to compact (unresolved) sources in each frequency band. A multifrequency catalog of Planck non-thermal sources is provided in Planck Collaboration Int. LIV (2018).

The scan strategy employed by Planck ensured that the full sky was covered twice in one year (Planck Collaboration 12014). Typically, a given source was observed at a given Planck frequency for a few days to one week in each 6-month survey (sources near the ecliptic poles, where Planck scans overlapped, were viewed much more frequently). HFI completed approximately five full sky surveys from 2009 August to 2012 January before its coolant ran out. LFI completed a fraction more than eight six-month surveys stretching from 2009 August to 2013 October. Thus we have five independent samples of flux densities at nine frequencies for each source, and an additional three or four measurements at 30, 44 and 70 GHz. The flux densities recorded in PCCS2 are averages over the entire mission, however, and thus smooth out any variability. In this paper, we work with the independent samples by analysing the time-ordered data (ToD).

Planck calibration at the frequencies employed here is determined from the satellite’s annual motion around the solar system barycenter (Planck Collaboration II 2020 for LFI; Planck Collaboration III 2020 for HFI). The calibration is thus absolute and precise to subpercent accuracy. This method of calibration has two consequences. First, the annual motion of the satellite produces a dipole signal in the CMB, much larger in angular scale than the beam width. In contrast the compact sources we investigate are almost never resolved, even at Planck’s highest frequency. Accurate calculations of flux density therefore require exact knowledge of the beam solid angle. The beam solid angles are known at each Planck band to percent-level accuracy (Planck Collaboration I 2020). Second, the substantial bandwidth of the Planck detectors means that they respond differently to sources with spectra that differ from the blackbody spectrum of the CMB. Thus, for precision work, each measured flux density needs a small, spectrum-dependent, “color correction.” We return to this issue in Sect. 4 below.

2.2 Source selection

Given Planck’s relatively low sensitivity to compact sources and the low polarization fraction expected (a few percent; see, e.g., Datta et al. 2019), we elected to work only with the brightest sources in the Planck PCCS2 catalog (Planck Collaboration XXVI 2016). Specifically, we selected all PCCS2 sources with 100 GHz flux density greater than 1 Jy (roughly 15–20 times the 1σ uncertainty). Recall that the cataloged Planck flux densities are averages over the entire duration of the mission. Since most of these bright sources are variable, some may have fallen below the 1 Jy threshold for part of that period. The PCCS2 catalog yielded 148 sources above the 1 Jy threshold at 100 GHz; 111 were at Galactic latitude >20° and therefore likely to be extragalactic. This list included many well-known AGN such as 3C 454.3 and BL Lac. To this initial list, we added five sources well-studied in other high-frequency monitoring programs. In the end, we examined 153 sources, both Galactic and extragalactic. The majority were extragalactic, and most of these were blazars. The full list of sources with classification is shown in Table A.2. Positions of all of the sources are indicated in Fig. 2.

2.3 NPIPE processing

The Planck NPIPE data release (also known as PR4) differs from earlier data releases in several critical ways (Planck Collaboration Int. LVII 2020). It is the first data release based on a single, stand-alone pipeline that processes both raw LFI and HFI data into calibrated, bandpass-corrected, noise-subtracted frequency maps. The release comprises frequency maps, Monte Carlo simulations, and calibrated timelines which we used as inputs to our processing. The timelines are accompanied by an updated instrument model that includes notable corrections to HFI polarization angles and efficiencies.

The NPIPE release is the first to include detector data acquired during so-called repointing maneuvers. Including these 4-min periods between stable science scans not only increases integration time by about 9% but also improves sky sampling by adding detector samples between the regular scanning paths, which are separated by two arcminutes. NPIPE suppresses degree-scale noise in the maps by modeling the 1/f noise fluctuations in the data with 167 ms noise offsets. In contrast, earlier LFI data were fitted with one-second offsets and earlier HFI data with full pointing period (45-min) offsets.

We list here other systematic effects and artifacts corrected in NPIPE processing:

• electrical interference from housekeeping electronics (for LFI) and the 4-kelvin cooler (for HFI);

• jumps in signal offset;

• cosmic ray glitches (HFI);

• bolometric nonlinearity (HFI);

• thermal common modes (HFI);

• ADC nonlinearity (HFI);

• gain fluctuations;

• bandpass mismatch from both continuum emission and CO lines;

• bolometric transfer function and transfer function residuals (HFI);

• far-sidelobe pickup.

In addition, NPIPE removes time-dependent signals from the orbital dipole and Zodiacal light that would otherwise compromise mapmaking.

Fig. 1 Footprint of the Planck focal plane as seen by an observer at infinity. The size of each colored spot is a rough indication of the relative resolution. The small crosses indicate polarization sensitivity.

3 Technical approach

As noted above, Planck sees most compact sources on the sky only once during each six-month sky survey. We call the roughly week-long period during which the Planck detectors collect time samples in the vicinity of a source an observation. During each observation, the orientation of the detectors with respect to the source remains fixed. Consequently, individual detectors cannot resolve the polarized flux density from the intensity in a single observation. To overcome this restriction, different detectors in a given Planck band are given different polarization angles (see Fig. 1).

Fig. 2 Positions of all 153 sources plotted over the Planck 857 GHz intensity map.

3.1 Individual detectors

To solve for the polarized flux density, we first measured the polarization-modulated flux density seen by each of the Planck detectors. The data model can be cast in the typical linear regression form:

$$d=Fa+n,$$(1)

where

• d is the vector of detector samples that fall within some fitting radius from the source coordinates;

• F contains the time-domain templates representing the source and the background. The primary template in F is the detector point spread function (PSF). This is the beam pattern produced when a detector with finite resolution scans a point source.

• a contains the template amplitudes;

• n is the detector noise.

The maximum likelihood solution for the template amplitudes, a, is:

$$a={\left(F^{\text{T}}{N}^{-1}F\right)}^{-1}{F}^{\text{T}}{N}^{-1}d\text{,}$$(2)

where N is the noise covariance matrix N = 〈nn^T〉. The Planck PR4 timelines used as inputs are calibrated and destriped, making a diagonal noise matrix, N = σ²1, an excellent approximation. Furthermore, the detector samples in d are separated in time by 60 s, the time it takes for Planck to spin around, effectively breaking sample-sample noise correlations.

The PSF must be rotated into the appropriate orientation and centered at the cataloged source position. The fitted amplitude of the PSF template is readily translated into a flux density across the detector pass band. Since the detector data are calibrated in thermodynamic units, the fitted amplitudes must be translated into flux density using the unit conversion factors derived from ground-measured band-passes. The full unit conversion is

$$f=T_{\max }{\text{Ω}}_{\det }{u}_{\det },$$(3)

where T_max is the peak temperature read off the PSF fit, Ω_det is the solid angle of the PSF and u_det is the unit conversion factor from kelvin to MJy sr⁻¹, derived by convolving the detector bandpass with a standard source spectrum (slope α = −1).

Additional columns in F model the background [foreground] emission. There is always a constant offset template, but we also considered higher-order templates such as slopes. Keeping the fitting radius close to the beamwidth limits the number of modes needed to fully model the background [foreground].

Small offset between the catalog source position and apparent source location can be modeled by adding spatial derivatives of the PSF in F. Such offsets naturally rise if the source position was originally fitted over a pixelized sky map. There can also be small pointing errors in the Planck data themselves. Finally, the center of the PSF template is slightly arbitrary and can be defined as either the peak or the center of weight. For vast majority of our targets, including position corrections in the fit increases uncertainty without a significant improvement in the quality of the fit. The few cases where the improvement is statistically significant (source is very bright), the change in measured flux density is less than a percent. Therefore the reported flux densities are based on fits directly at the catalog positions.

3.2 Polarized flux density

For a given observation, each detector produces one data point: the polarization-modulated flux density:

$$f_{\text{d}}=I+{\eta }_{\text{d}}\left[Q\cos \left(2\left(\psi +{\psi }_{\text{d}}\right)\right)+U\sin \left(2\left(\psi +{\psi }_{\text{d}}\right)\right)\right],$$(4)

where

• f_d is the polarization-modulated flux density measured with detector d;

• I, Q and U are the Stokes parameters of the source flux density (CMB detectors are not typically sensitive to V, the circular polarization component);

• η is the detector polarization efficiency – an ideal polarization-sensitive detector would have η = 1;

• ψ represents the angle between the source coordinate system and the focal-plane coordinate system;

• ψ_d is the polarization sensitive direction of the detector in the focal-plane coordinates.

A vector of detector fluxes, f, can be modeled as

$$f=Pm+n,$$(5)

where each row of the pointing matrix, P, comprises the pointing weights, [1, η cos(2(ψ + ψ_d)), η sin(2(ψ + ψ_d))], for one detector, and the “map” vector contains the three Stokes parameters: m = [I, Q, U]^T. The noise vector, n, contains the measurement uncertainty for each detector. Analogous to Eq. (2), the maximum likelihood solution reads:

$$m={\left(P^{\text{T}}{N}^{-1}P\right)}^{-1}{P}^{\text{T}}{N}^{-1}f.$$(6)

The theory of linear regression also provides the covariances:

$$N_{\text{m}}=\langle \text{Δ}m\text{Δ}{m}^{\text{T}}\rangle ={\left(P^{\text{T}}{N}^{-1}P\right)}^{-1}=\left[\begin{array}{ccc}{\sigma }_I^2& {\sigma }_{IQ}& {\sigma }_{IU}\\ & {\sigma }_Q^2& {\sigma }_{QU}\\ & & {\sigma }_U^2\end{array}\right],$$(7)

where we have omitted the redundant elements in the symmetric covariance matrix.

In order to derive the polarization fraction and polarization angle, we defined two estimators:

$$\begin{array}{ccc}\widehat{p}=\frac{\sqrt{Q^2+U^2}}{I}& \text{and}& \widehat{\psi }=\frac{1}{2}\arctan \frac{U}{Q}.\end{array}$$(8)

While these estimators are textbook definitions of the polarization parameters, they are biased by uncertainties in I, Q, and U. Finding unbiased or minimally biased estimators for these parameters turns out to be complicated and depends on the signal-to-noise ratio (S/N) of the I, Q, U estimates. For the modest S/N in our case, we have adopted the asymptotic estimator from Montier et al. (2015):

$${\widehat{p}}_{\text{AS}}=\{\begin{array}{ll}\sqrt{\widehat{p}-b^2}\hfill & \text{if\hspace{0.17em}}\widehat{p}>b,\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array},$$(9)

where

$$b^2=\frac{{\tilde{\sigma }}_U^2{\cos }^2\left(2{\psi }_0-\theta \right)+{\tilde{\sigma }}_Q^2{\sin }^2\left(2{\psi }_0-\theta \right)}{I_0^2}$$(10)

$$\theta =\frac{1}{2}\arctan \left(\frac{2\rho {\sigma }_Q{\sigma }_U}{{\sigma }_Q^2-{\sigma }_U^2}\right)$$(11)

$${\tilde{\sigma }}_Q^2={\sigma }_Q^2{\cos }^2\theta +{\sigma }_U^2{\sin }^2\theta +\rho {\sigma }_Q{\sigma }_U\sin 2\theta$$(12)

$${\tilde{\sigma }}_U^2={\sigma }_Q^2{\sin }^2\theta +{\sigma }_U^2{\cos }^2\theta -\rho {\sigma }_Q{\sigma }_U\sin 2\theta$$(13)

$$\rho =\frac{{\sigma }_{QU}}{{\sigma }_Q{\sigma }_U}.$$(14)

The estimator from Eq. (8) is typically substituted in place of the true polarization angle, ψ₀.

Confidence regions for the polarization parameters are determined using a Monte Carlo method: we sampled 100 000 triplets of true (I₀, p₀,ψ₀) and recorded the likelihood of the measured (I, p, ψ). From those likelihoods, we inferred the 68% and 95% confidence regions for the true polarization parameters. The Bayesian likelihood, ℒ(p,ψ|I₀, p₀, ψ₀), is defined in Montier et al. (2015), Eqs. (3) and (23). Examples of the sampled likelihood, projected to the p-axis are shown in Fig. 3.

Fig. 3 100 000 sampled values of the likelihood function, projected to the intensity axis (top row), polarization fraction axis (second row) or polarization angle axis (third row). The source, PCCS2 143 G052.40–36.49 also known as 4C–02.81, was chosen because it has both positive detections of and upper limits on polarization fraction. The likelihood threshold that includes 68% of the total likelihood is indicated in green and 95% in orange. Samples with likelihood lower than 95% of the distribution are in blue. When the green part of the polarization fraction distribution excludes zero, we report a 68% confidence interval. When zero is included, we report a 95% upper limit. The bottom row presents the maximum of the likelihood function in the polarization fraction-polarization angle plane. The dashed lines indicate the maximum likelihood values. The likelihoods are based on the first survey shown in Fig. 4.

3.3 Multifrequency fits

In addition to tabulating flux densities at each Planck frequency, we provide multifrequency fits in order to improve the S/N of the measured polarized flux densities. A fit such as this necessarily requires some accommodation of the source spectral energy distribution (SED). We extended the signal model in Eq. (4) to:

$$f_{\text{d}}={\left(\frac{v_{\text{d}}}{v_0}\right)}^{\alpha }\left\{I+{\eta }_{\text{d}}\left[Q\cos \left(2\left(\psi +{\psi }_{\text{d}}\right)\right)+U\sin \left(2\left(\psi +{\psi }_{\text{d}}\right)\right)\right]\right\},$$(15)

where α is the spectral index of the source, ν_d is the central frequency of detector d, and ν₀ is the noise-weighted central frequency of all detectors that are being combined. This model implies strong constraints on the fit:

• the SED is approximated as a power-law across the fitted frequencies;

• the spectral indices of intensity and polarized emission are the same;

• the polarization fraction and angle do not depend on the frequency.

These approximations obviously break down as increasingly wider frequency ranges are considered. Hence we provide multi-frequency fits only up to three adjacent frequency bands. Within these bands, and to the extent the approximations are valid, the additional constraints greatly improve the S/N of our fits.

Fig. 4 Example of a full diagnostic plot showing intensity, polarization fraction, polarization angle, and position for PCCS2 143 G052.40–36.49. The polarization fraction axis is nonlinear between 0.3 and 1.0 to accommodate all values. The lines and points of different colours mark different observing epochs, separated by 6 months. The Galactic-coordinate map (HEALPix formatted, Górski et al. 2005) in the lower right panel shows the position of the source on top of the Planck 857 GHz intensity map. Plots and maps similar to this are provided for each source and each number of jointly fit frequencies.

3.4 Software

Our flux-fitting software is implemented in Python-3 and released as open source software². The package includes Planck-specific specializations but they are modularized and the software can be adapted for other experimental data sets that can be sampled into the small data set format used as inputs to the calculation.

4 The catalog

Our catalog of time- and frequency-dependent, polarized, flux densities is spread over multiple comma-separated value (CSV) files. Each file features a single source and the entries correspond to a fixed number (ranging from 1 to 3) of frequency bands averaged together. The filenames are of the form results_<PCCSNAME>_nfreq=<NFREQ>.csv, where <PCCSNAME> is the name of the source in the 143 GHz Planck PCCS2 catalog and <NFREQ> is the number of frequency bands averaged for each entry (1, 2, or 3). The fields of the files are

1. band(s) – nominal Planck bands comprising the entry.

2. freq [GHz] – noise-weighted average of the nominal frequencies.

3. start – approximate UT calendar date of the start of the observations used in deriving the entry.

4. stop – approximate UT calendar date of the end of the observations.

5. I flux [mJy] – intensity flux density at the effective central frequency, not color-corrected to match freq [GHz].

6. I error [mJy] – 1-σ uncertainty of the I flux [mJy].

7. Q flux [mJy] – linear Q-polarization flux density in IAU convention and Equatorial (J2000) coordinates.

8. Q error [mJy] – 1-σ uncertainty of Q flux [mJy].

9. U flux [mJy] – linear U-polarization flux density in IAU convention and Equatorial (J2000) coordinates.

10. U error [mJy] – 1-σ uncertainty of U flux [mJy].

11. Pol.Frac – debiased polarization fraction in the range [0 … 1].

12. PF low lim – Bayesian 16% quantile lower limit of Pol.Frac or zero.

13. PF high lim – If PF low lim is nonzero the Bayesian 16% quantile upper limit of Pol.Frac. Otherwise the 95% upper limit.

14. Pol.Ang [deg] – Polarization angle in IAU convention and Equatorial (J2000) coordinates.

15. PA low lim – Bayesian 16% quantile lower limit of Pol.Ang [deg].

16. PA high lim – Bayesian 16% quantile upper limit of Pol.Ang [deg].

An excerpt of one catalog file is shown in Appendix A. The CSV file is accompanied by two diagnostic plots that visualize its contents:

• flux_fit_<PCCSNAME>_nfreq=<NFREQ>.png;

• results_<PCCSNAME>_nfreq=<NFREQ>.replot.png.

The first of these shows the fit data as a function of frequency (e.g., Fig. 4) and the second as a function of time (e.g., Fig. 5).

For convenience, we also provide a single catalog file, combined.csv, that combines entries from all individual results_*csv files so that it includes data for all 153 sources. This file has three additional columns appended to the beginning of each row:

• target – Name of the target in the Planck 143GHz PCCS2 catalog;

• RA – Right ascension in degrees (J2000);

• Dec – Declination in degrees (J2000).

The flux densities in our catalog are not color-corrected. That is, we made no effort to translate the measured flux density to the nominal observing frequency. The reason is that the spectral indices of the sources can vary as a function of time and across the Planck pass-band, and the small corrections (of order a few percent) depend on a source’s spectral index. Instead, we present in Table 2 the effective central frequencies of the Planck bands for a number of spectral slopes. The formulas for evaluating the central frequency are laid out in Planck Collaboration IX (2014).

Our catalog, PCCS-PV, will be available at the ESA Planck Legacy Archive³. and at the NASA/IPAC Infrared Science Archive⁴.

Fig. 5 Time-series plot of our flux density fits to 3C 454.3 observations. The HFI data for the third pass are missing due to quality flagging.

5 Results and discussion

5.1 Polarized signal-to-noise ratio

The targets in our catalog vary in intensity, polarization fraction and spectral properties. We explore the polarized S/N of our catalog in Fig. 6. The figure shows the scatter in polarized S/N as a function of the intensity flux density and identifies a threshold beyond which half of the observations reach S/N of one. The threshold varies considerably across the Planck frequencies and obviously depends on the spectral and polarization properties of our catalog.

Our most sensitive, single frequency measurements of polarized flux density come from the 143 GHz band where the S/N is higher than unity in a majority of observations when the source is brighter than 1.86 Jy. Fitting for the source spectral index and combining frequencies gives us a combined threshold better than 1.5 Jy in all frequency combinations involving 143 GHz. Combining the Planck LFI frequencies (30–70 GHz) has a combined threshold of 2.20 Jy.

Fig. 6 Polarized S/N: for each polarization-sensitive observation in our catalog, we plot the de-biased polarized flux density S/N against the intensity flux density. Observing bands are indicated in the labels. For cases where de-biasing renders a S/N less than 0.1, we plot a data point at 0.1. The vertical grey bands indicate the lowest and narrowest flux density bin that has a median S/N of 1.0 and at least 20 entries. The center of the bin is spelled out in the label.

5.2 Characterization

We made use of both spectral information and classifications from the literature, mainly NED (Helou et al. 1995) and SIMBAD (Wenger et al. 2000), to divide our sample of 153 sources into several categories. The majority are extragalactic, but nine (6%) are Galactic sources such as H II regions or planetary nebulae, and five more are H II regions located in the Magellanic Clouds. Many of the extragalactic objects are well-studied radio sources such as 3C 84 or M87. We characterize the 139 extragalactic sources as ~25% BL Lac objects, ~60% quasars and ~10% Seyfert-1 galaxies. Some sources exhibit both synchrotron emission with a falling spectrum and dust emission at the highest Planck frequencies (see Fig. 8 for an example).

The few Galactic (and LMC) sources are not beamed and have large physical dimensions, so they are unlikely to be variable on time scales of a few years or less. They therefore offer an opportunity to check the fidelity of our pipeline. Figure 7 provides an example. PCCS2 143 G000.35–19.48 is a Galactic dark cloud, and, as expected, shows no statistically significant variation with epoch at any Planck frequency except for some scatter at 44 GHz, which has the highest noise among the Planck bands. In addition, as shown in Fig. 1, the three 44 GHz beams lie at the periphery of the Planck field of view, and two of the beams have much higher ellipticity than the other one. Consequently, the measured flux density at that frequency depends on how the source crosses the focal plane, and hence may show some instrument-related variation.

Some extragalactic sources also demonstrate no significant variation during the Planck mission. We show two examples in Fig. 8. G021.20+19.62 (3C 353) is an radio galaxy with a spectral upturn at 545 GHz showing no variation in flux density at any frequency. The polarization percentage, while small, also does not appear to vary.

Of more interest are those sources that clearly do vary during the Planck mission. The blazar PCCS2 143G086.11–38.19 (3C 454.3), for instance, is a well-known and well-studied variable source that is known to have flared during the Planck mission (Planck Collaboration XIV 2011). We discuss 3C 454.3 in some detail below (see Fig. 5).

In Fig. 9, we show the SED of another strongly variable source, G304.15–72.18 at RA= 12.32°, Dec= −44.95°. There is clear evidence of spectral change with epoch. Strong variability is also evident in the following sources, among others: G106.95–50.62 (a Seyfert-1 galaxy), G315.80–36.52 (a QSO), and G352.45–08.41 (a blazar).

Fig. 7 SED of a Galactic source, PCCS2143 G000.35–19.48 at RA = 287.56° and Dec = −37.15°. The lines and points of different colors mark different observing epochs, separated by 6 months. It is important to note the larger scatter of the noisier LFI measurements at 30, 44, and 70 GHz. With the possible exception of 44 GHz, there is no evidence for variability, as expected for a large source of this type.

Fig. 8 Same as in Fig. 7, but the lines and points of different colors mark different observing epochs, separated by 6 months. For the polarization fraction (right column) we plot the maximum likelihood value as a solid line and indicate a 68% confidence region or 95% upper limit with error bars. Top: An AGN (3C 353) displaying no significant changes in SED, and small, but apparently constant polarization. Bottom: a BL Lac object, PKS 1514–24, showing only mild variability, and no significant change in polarization.

Fig. 9 Strongly variable source, PCCS2143 G304.15–72.18. It is apparent that the spectral shape, as well as luminosity, changes with epoch.

Fig. 10 Flux density and polarization versus time of the BL Lac object OJ 287. Frequencies and sources of the data are indicated in the figure. The figure does not plot the following upper limits: 44 GHz: MJD 55301-55309 (≤34%), 55497–55505 (≤30%), 55666–55674 (≤24%), 56038–56050 (≤27%), 56236–56247 (≤42%); 70 GHz: MJD 56238–56241 (≤31%), 56410–56413 (≤31.2%); 353 GHz: MJD 55304–55305 (≤77%), 55500–55501 (≤52%), 55669–55670 (≤37%). No polarization was detected at 545 or 857 GHz, since the HFI detectors are unpolarized at those frequencies.

5.3 Concurrent observations with other instruments

In order to check our measurements and to illustrate possible uses of the data presented here, we compare the Planck polarization data of selected objects with published data from other astronomical instruments. Figure 10 presents the Planck data of the BL Lac object OJ287 alongside the 37 GHz light curve from the Metsähovi Radio Observatory (as summarized in Planck Collaboration Int. XLV 2016), linear polarization at 86 and 230 GHz measured with the 30-meter antenna of the Institut de Radioastronomie Millimétrique (IRAM, Agudo et al. 2018), and Very Long Baseline Array (VLBA) images that include linear polarization (Jorstad et al. 2017). Fig. 11 provides a similar plot for the quasar 3C 454.3. For comparison, we present in Figs. 12 and 13 VLBA images at 43 GHz of OJ 287 and 3C 454.3 obtained within one month of an epoch for which we have Planck polarization measurements. The images are from the VLBA-BU Blazar Monitoring Program⁵. We find that the flux densities at 30 and 44 GHz agree with the Metsähovi 37 GHz values at nearby epochs. In addition, there is general correspondence between the Planck and IRAM polarization. These similarities serve to validate our measurements.

Figure 10 demonstrates that the Planck polarization of OJ287 is similar to that of the core seen in Fig. 13. We note that the optical electric-vector position angle (EVPA) of the knot deviates at some epochs from the source-integrated polarization at 30-353 GHz. Sasada et al. (2018) have shown that the optical polarization is related to that of either the core or the knot, depending on the epoch. Our data are consistent with this conclusion.

Figure 14 presents the frequency dependence of the flux density, degree of polarization, and EVPA of 3C 454.3 during the first Planck survey in late 2009. At that epoch, an extremely high-amplitude two-year outburst was in its early, rising stage. The core completely dominated the polarization in a VLBA image obtained 17 days earlier (from the VLBA-BU-BLAZAR website). Figure 14 reveals that the degree of polarization was higher by a factor of ~2 from 44–100 GHz below the spectral turnover than at 217–353 GHz above the turnover. There is a frequency gradient in EVPA from ~5° to ~−30° from 70 to 353 GHz. At the same epoch, the optical R-band EVPA was −37° ± 2.1°, having shifted from +25° nine days earlier (Jorstad et al. 2013). The degree of polarization at R band was 15.4 ± 1.8%. This frequency dependence strongly implies that the outburst involved a component that dominated the emission from ~200 GHz to optical frequencies. We are therefore able to conclude that the optical flare occurred within the same region as the ~200 GHz flare, which is distinct from the regions (probably farther downstream in the jet) that dominate the emission at lower frequencies.

Fig. 11 Flux density and polarization versus time of the quasar 3C 454.3. Frequencies and sources of the data are indicated in the figure. The figure does not plot the following upper limits: 30 GHz: MJD 55179–55182 (≤2.4%), 56464–56465 (≤11%); 44 GHz: MJD 55357–55369 (≤3.6%), 55535–55546 (≤4.6%), 55899–55914 (≤36%), 56089–56097 (≤42%), 56277–56285 (≤41%), 56464–56472 (≤29%); 70 GHz: MJD 55912–55913 (≤47%), 56101–56103 (≤64%), 56277–56279 (≤36%), 56465–56467 (≤18%); 217 GHz: MJD 55902–55902 (≤27%); 353 GHz: MJD 55901–55901 (≤100%). No polarization was detected at 545 and 857 GHz.

Fig. 12 VLBA image of OJ 287 at 43 GHz at epoch 2012 December 21 (MJD 56282). The contours represent total intensity, starting at 0.25% of the maximum of 0.79 Jy beam−1 and increasing by factors of 2. The color scale represents linearly polarized intensity, with maximum value (yellow) of 0.032 Jy beam−1. The yellow line segments indicate the EVPA. The restoring beam FWHM, representing the angular resolution, is displayed in the lower left corner.

Fig. 13 VLBA image of 3C 454.3 at 43 GHz at epoch 2011 May 22 (MJD 55704). The contours represent total intensity, starting at 0.25% of the maximum of 6.00 Jy beam−1 and increasing by factors of 2. The color scale represents linearly polarized intensity, with maximum value (yellow) of 0.27 Jy beam−1. The yellow line segments indicate the EVPA. The restoring beam FWHM, representing the angular resolution, is displayed in the lower left corner.

Fig. 14 Frequency dependence of the flux density, degree of linear polarization, and EVPA of 3C 454.3 in 2009 December. For frequencies ≤353 GHz, the data are from this work. The high-frequency point is at optical R band from Jorstad et al. (2013).

5.4 Planck’s contribution to the overall spectral energy distribution

In Figs. 15 and 16, we present the Planck data from PCCS2 (Planck Collaboration XXVI 2016). For the remaining three sources, Centaurus A, 3C 84, and PictorA, the Planck data help to fill in some missing spectral coverage. In particular, in Centaurus A, the Planck data help to define the rise into the spectral region dominated by the infrared dust component.

6 Summary

In this paper, we have introduced the new Planck Catalog of Polarized and Variable Compact Sources, PCCS-PV, comprising 153 sources, the majority of which are extragalactic, measured in both total power and polarization by Planck, with frequency coverage from 30 to 353 GHz, and time-scales ranging from days to years. We classify 135 of the 153 sources as beamed, extragalactic radio sources (blazars), four as well-studied radio galaxies such as 3C 84 and M87, and 14 as Galactic or Magellanic Cloud sources, including H II regions, planetary nebulae, etc. To allow an assessment of variability of polarized sources, the catalog was generated using an advanced extraction method, tod2flux, applied directly to the multifrequency Planck time-ordered data, rather than the mission sky maps. We used the calibrated timelines from the Planck NPIPE data release, PR4, as input to our processing.

To check our measurements and to illustrate possible uses of our catalog, we have compared the Planck polarization data of selected objects with published data from other astronomical instruments. Our preliminary findings show that for the selected sources OJ 287 and 3C 454.3, the flux densities at 30 and 44 GHz agree with the Metsähovi 37 GHz values at nearby epochs, while there is general correspondence between the Planck and IRAM polarization. This general agreement confirms our measurements. Furthermore we have found that the Planck polarization of OJ 287 is similar to that of the core, supporting the conclusion of Sasada et al. (2018) that the optical polarization is related to that of the core at some epochs and the knot at others. From joint analysis of our data and optical R-band data, we have shown a strong frequency dependency of the degree of polarization and EVPA of 3C 454.3, with the EVPA at the highest Planck frequencies becoming similar to the optical EVPA during the early stages of a major outburst. This has led us to conclude that the optical flare during the first Planck survey in late 2009 occurred within the same region as the 200 GHz flare, which is distinct from the regions that dominate the emission at lower frequencies.

We have also presented Planck data along with observations compiled in NED for five non-blazar extragalactic sources to demonstrate Planck’s contribution to the overall SED of the sources. Planck observations of M87 and the quasar 2134+004 are in good agreement with the observations in NED, confirming the flux density calibration of the Planck data. For the remaining three sources, Centaurus A, 3C 84, and Pictor A, the Planck data help to fill in missing spectral coverage. In particular, in Centaurus A, the Planck data help to define the rise as the infrared dust component starts to dominate the spectrum.

These preliminary results illustrate the usefulness of our Planck catalog, PCCS-PV, for application with follow-up studies. Furthermore, our versatile ToD-based flux-fitting software, tod2flux, is applicable to other experimental data sets.

PCCS-PV will be available at the ESA Planck Legacy Archive⁶. and at the NASA/IPAC Infrared Science Archive⁷, while the flux-fitting algorithm, tod2flux, is released as open source software⁸.

Fig. 15 Spectral energy distribution (SED) for a few sources: red points are photometry from NED and the blue points are our Planck observations. The Planck data for Cen A add some new spectral coverage and help to define the rise of the dust component.

Fig. 16 Spectral energy distribution (SED) for a few sources: red points are photometry from NED and the blue points are our Planck observations. The Planck observations of 2134+004 and M87 show good agreement with previous measurements.

Acknowledgements

The authors would like to acknowledge Charles Lawrence for insightful comments that helped improve this paper and Mainak Singha for his help with Figs. 15 and 16 and Mike Peel, Clive Dickinson and Ricardo Tanausú Génova Santos for pointing out diculties in interpreting a previous version of Table 2 and an associated error in LFI bandpass handling. This research was conducted under the auspices of a NASA Astrophysics and Data Analysis Program (ADAP) award NNH18ZDA001N-ADAP. This work was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration⁹. This study made use of Very Long Baseline Array (VLBA) data from the VLBA-BU-BLAZAR project, funded by NASA through the Fermi Guest Investigator Program. The VLBA is an instrument of the National Radio Astronomy Observatory. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated by Associated Universities, Inc. The work at Boston University was supported in part by NASA through Fermi Guest Investigator Program grants 80NSSC17K0649 and 80NSSC20K1567. C. O’Dea received support from the Natural Sciences and Engineering Research Council (NSERC) of Canada. This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11231 This research has made use of the SIMBAD database (Wenger et al. 2000), operated at CDS, Strasbourg, France. Softwares used in this work: Astropy (Astropy Collaboration 2013, 2018), HEALPix (Górski et al. 2005), Matplotlib (Hunter 2007; Caswell et al. 2021), Numpy (Harris et al. 2020), TOAST¹⁰, and Scipy (Virtanen et al. 2020).

Appendix A Sample catalog file

We show an example of a catalog file in Table A.1. This one is for These results combine two adjacent frequency bands to boost PCCS2 143 G009.33-19.61, better known as QSO B1921-293. the signal−to−noise ratio of each entry.

References

All Tables

All Figures

	Fig. 1 Footprint of the Planck focal plane as seen by an observer at infinity. The size of each colored spot is a rough indication of the relative resolution. The small crosses indicate polarization sensitivity.
In the text

	Fig. 2 Positions of all 153 sources plotted over the Planck 857 GHz intensity map.
In the text

Fig. 3 100 000 sampled values of the likelihood function, projected to the intensity axis (top row), polarization fraction axis (second row) or polarization angle axis (third row). The source, PCCS2 143 G052.40–36.49 also known as 4C–02.81, was chosen because it has both positive detections of and upper limits on polarization fraction. The likelihood threshold that includes 68% of the total likelihood is indicated in green and 95% in orange. Samples with likelihood lower than 95% of the distribution are in blue. When the green part of the polarization fraction distribution excludes zero, we report a 68% confidence interval. When zero is included, we report a 95% upper limit. The bottom row presents the maximum of the likelihood function in the polarization fraction-polarization angle plane. The dashed lines indicate the maximum likelihood values. The likelihoods are based on the first survey shown in Fig. 4.

In the text

Fig. 4 Example of a full diagnostic plot showing intensity, polarization fraction, polarization angle, and position for PCCS2 143 G052.40–36.49. The polarization fraction axis is nonlinear between 0.3 and 1.0 to accommodate all values. The lines and points of different colours mark different observing epochs, separated by 6 months. The Galactic-coordinate map (HEALPix formatted, Górski et al. 2005) in the lower right panel shows the position of the source on top of the Planck 857 GHz intensity map. Plots and maps similar to this are provided for each source and each number of jointly fit frequencies.

In the text

	Fig. 5 Time-series plot of our flux density fits to 3C 454.3 observations. The HFI data for the third pass are missing due to quality flagging.
In the text

Fig. 6 Polarized S/N: for each polarization-sensitive observation in our catalog, we plot the de-biased polarized flux density S/N against the intensity flux density. Observing bands are indicated in the labels. For cases where de-biasing renders a S/N less than 0.1, we plot a data point at 0.1. The vertical grey bands indicate the lowest and narrowest flux density bin that has a median S/N of 1.0 and at least 20 entries. The center of the bin is spelled out in the label.

In the text

Fig. 7 SED of a Galactic source, PCCS2143 G000.35–19.48 at RA = 287.56° and Dec = −37.15°. The lines and points of different colors mark different observing epochs, separated by 6 months. It is important to note the larger scatter of the noisier LFI measurements at 30, 44, and 70 GHz. With the possible exception of 44 GHz, there is no evidence for variability, as expected for a large source of this type.

In the text

Fig. 8 Same as in Fig. 7, but the lines and points of different colors mark different observing epochs, separated by 6 months. For the polarization fraction (right column) we plot the maximum likelihood value as a solid line and indicate a 68% confidence region or 95% upper limit with error bars. Top: An AGN (3C 353) displaying no significant changes in SED, and small, but apparently constant polarization. Bottom: a BL Lac object, PKS 1514–24, showing only mild variability, and no significant change in polarization.

In the text

	Fig. 9 Strongly variable source, PCCS2143 G304.15–72.18. It is apparent that the spectral shape, as well as luminosity, changes with epoch.
In the text

Fig. 10 Flux density and polarization versus time of the BL Lac object OJ 287. Frequencies and sources of the data are indicated in the figure. The figure does not plot the following upper limits: 44 GHz: MJD 55301-55309 (≤34%), 55497–55505 (≤30%), 55666–55674 (≤24%), 56038–56050 (≤27%), 56236–56247 (≤42%); 70 GHz: MJD 56238–56241 (≤31%), 56410–56413 (≤31.2%); 353 GHz: MJD 55304–55305 (≤77%), 55500–55501 (≤52%), 55669–55670 (≤37%). No polarization was detected at 545 or 857 GHz, since the HFI detectors are unpolarized at those frequencies.

In the text

Fig. 11 Flux density and polarization versus time of the quasar 3C 454.3. Frequencies and sources of the data are indicated in the figure. The figure does not plot the following upper limits: 30 GHz: MJD 55179–55182 (≤2.4%), 56464–56465 (≤11%); 44 GHz: MJD 55357–55369 (≤3.6%), 55535–55546 (≤4.6%), 55899–55914 (≤36%), 56089–56097 (≤42%), 56277–56285 (≤41%), 56464–56472 (≤29%); 70 GHz: MJD 55912–55913 (≤47%), 56101–56103 (≤64%), 56277–56279 (≤36%), 56465–56467 (≤18%); 217 GHz: MJD 55902–55902 (≤27%); 353 GHz: MJD 55901–55901 (≤100%). No polarization was detected at 545 and 857 GHz.

In the text

Fig. 12 VLBA image of OJ 287 at 43 GHz at epoch 2012 December 21 (MJD 56282). The contours represent total intensity, starting at 0.25% of the maximum of 0.79 Jy beam−1 and increasing by factors of 2. The color scale represents linearly polarized intensity, with maximum value (yellow) of 0.032 Jy beam−1. The yellow line segments indicate the EVPA. The restoring beam FWHM, representing the angular resolution, is displayed in the lower left corner.

In the text

Fig. 13 VLBA image of 3C 454.3 at 43 GHz at epoch 2011 May 22 (MJD 55704). The contours represent total intensity, starting at 0.25% of the maximum of 6.00 Jy beam−1 and increasing by factors of 2. The color scale represents linearly polarized intensity, with maximum value (yellow) of 0.27 Jy beam−1. The yellow line segments indicate the EVPA. The restoring beam FWHM, representing the angular resolution, is displayed in the lower left corner.

In the text

	Fig. 14 Frequency dependence of the flux density, degree of linear polarization, and EVPA of 3C 454.3 in 2009 December. For frequencies ≤353 GHz, the data are from this work. The high-frequency point is at optical R band from Jorstad et al. (2013).
In the text

	Fig. 15 Spectral energy distribution (SED) for a few sources: red points are photometry from NED and the blue points are our Planck observations. The Planck data for Cen A add some new spectral coverage and help to define the rise of the dust component.
In the text

	Fig. 16 Spectral energy distribution (SED) for a few sources: red points are photometry from NED and the blue points are our Planck observations. The Planck observations of 2134+004 and M87 show good agreement with previous measurements.
In the text