# A discussion on the anomalous threshold enhancement of $J/\psi$ –
$\psi(3770)$ couplings and $X(6900)$ peak

###### Abstract

The attractive interaction between $J/\psi$ and $\psi(3770)$ has to be strong
enough to form the molecular resonant structure, i.e., $X(6900)$. We argue
that since $\psi(3770)$ decays predominantly into a $D\bar{D}$ pair, the
interactions between $J/\psi$ and $\psi(3770)$ may be significantly enhanced
due to the three point $D\bar{D}$ loop diagram. The enhancement comes from the
anomalous threshold located at $t=-1.288$GeV2, which effect propagates into
the $s$-channel partial wave amplitude in the vicinity of $\sqrt{s}\simeq
6.9$GeV. This effect may be helpful in forming the $X(6900)$ peak.

Ye Lu${}^{\dagger\,}$, Chang Chen${}^{\dagger\dagger\,}$, Guang-you
Qin${}^{\dagger\,}$, Han-Qing Zheng${}^{\heartsuit\,}$
† Institute of Particle Physics and Key Laboratory of Quark and Lepton Physics
(MOE), Central China Normal University, Wuhan, Hubei 430079, China
†† Department of Physics, Peking University, Beijing 100871, China
♡ College of Physics, Sichuan University, Chengdu, Sichuan 610065, China

The $X(6900)$ peak observed by the LHCb Collaboration in the di-$J/\psi$
invariant mass spectrum [1][2], and later in the $J\psi\psi(3686)$ invariant
mass spectrum [3]. has stimulated many discussions in theory aspects (see for
example Ref. [4] for an incomplete list of references). On the other side,
since $X(6900)$ is close to the threshold of $J/\psi\psi(3770)$,
$J/\psi\psi_{2}(3823)$, $J/\psi$$\psi_{3}(3842)$, and $\chi_{c0}\chi_{c1}$
(and $X(7200)$ is close to the threshold of $J/\psi\psi(4160)$ and
$\chi_{c0}$$\chi_{c1}(3872)$), inspired by this, it is studied in Ref. [5] the
properties of $X(6900)$ and $X(7200)$, by assuming the $X(6900)$ coupling to
$J/\psi J/\psi$, $J/\psi\psi(3770)$, $J/\psi\psi_{2}(3823)$,
$J/\psi\psi_{3}(3842)$ and $\chi_{c0}\chi_{c1}$ channels, and $X(7200)$ to
$J/\psi J/\psi$, $J/\psi\psi(4160)$ and $\chi_{c0}\chi_{c1}(3872)$ channels.
For the $S$-wave $J/\psi J/\psi$ coupling, the pole counting rule (PCR) [6],
which has been applied to the studies of “$XYZ$” physics in Refs. [7, 8, 9,
10] is employed to analyze the nature of the two structures. It is found that
the di-$J/\psi$ data alone is not enough to judge the intrinsic properties of
the two states. It is also pointed out that the $X(6900)$ is unlikely a
molecule of $J/\psi\,\psi(3686)$ [5] – a conclusion drawn before the discovery
of Ref. [3]. More recently, Refs. [11, 12, 4] made an analysis on $X(6900)$
using a combined analysis on di-$J/\psi$ and $J/\psi\psi(3686)$ data and
concluded that $X(6900)$ cannot be a $J/\psi\psi(3686)$ molecule.

Nevertheless, as already stressed in Ref. [4], even though $X(6900)$ is very
unlikely a molecule of $J/\psi\psi(3686)$, it does not mean that it has to be
an “elementary state” (i.e., a compact $\bar{c}\bar{c}cc$ tetraquark state).
It is pointed out that, it is possible that $X(6900)$ be a molecular state
composed of other particles, such as $J/\psi\psi(3770)$, which form thresholds
closer to $X(6900)$, if the channel coupling is sufficiently large.

This note will discuss a possible mechanism for the enhancement of the
$J/\psi\psi(3770)$ channel coupling. The $D\bar{D}$ (or $D^{*}\bar{D}^{*}$)
component inside $\psi(3770)$ may play an important role, so far ignored in
the literature, in explaining the $X(6900)$ resonant peak, through the
anomalous threshold emerged from the triangle diagram generated by $D$
($D^{*}$) loop, as depicted in Fig. 1.

Noticing that $\psi(3770)$ or $\psi^{\prime\prime}$ couples dominantly to
$D\bar{D}$, we start from the Feynman diagram as depicted in Fig. 1 by
assuming it contributes $J/\psi\psi^{\prime\prime}$ elastic scatterings near
$J/\psi\psi^{\prime\prime}$ threshold. 111It could be wondered that the
$X(6900)$ be a $\bar{D}DJ/\psi$ three body halo state, see Refs. [14, 15].

Figure 1: The $D\bar{D}$ triangle diagram in $J/\psi\psi^{\prime\prime}$
scattering process.

Assuming an interaction lagrangian222We neglect all the complexities such as
form factors, hence our calculations are only qualitative or at best semi-
quantitative.

$\displaystyle\mathcal{L}=$ $\displaystyle-
ig(D^{0}\partial_{\mu}\bar{D}^{0}-\bar{D}^{0}\partial_{\mu}D^{0})\psi^{\prime\prime\mu}-ig(D^{+}\partial_{\mu}D^{-}-D^{-}\partial_{\mu}D^{+})\psi^{\prime\prime\mu}$
(1)
$\displaystyle+g^{\prime}D^{0}\bar{D}^{0}J/\psi^{\mu}J/\psi_{\mu}+g^{\prime}D^{+}D^{-}J/\psi^{\mu}J/\psi_{\mu}\
,$

after performing the momentum integration the amplitude as depicted by Fig. 1
is

$\displaystyle i\mathcal{M}$
$\displaystyle=(-16g^{2}g^{\prime})~{}(\epsilon_{2}\cdot\epsilon_{4})~{}\epsilon_{1}^{\mu}\epsilon_{3}^{\nu}\int
d^{D}k\frac{k_{\mu}k_{\nu}}{((k-p_{1})^{2}-m_{D}^{2})((k-p_{3})^{2}-m_{D}^{2})(k^{2}-m_{D}^{2})}$
(2)
$\displaystyle\equiv(-16g^{2}g^{\prime})~{}(\epsilon_{2}\cdot\epsilon_{4})~{}\epsilon_{1}^{\mu}\epsilon_{3}^{\nu}~{}\mathcal{A}_{\mu\nu}$

where

$\displaystyle\mathcal{A}^{\mu\nu}$
$\displaystyle=\frac{-i}{16\pi^{2}}\int_{0}^{1}dx\int_{0}^{1-x}dy~{}\left\\{\frac{p_{3}^{\mu}p_{1}^{\nu}~{}xy}{\Delta}-\frac{g^{\mu\nu}}{4}\Gamma(\epsilon)\frac{1}{\Delta^{\epsilon}}\right\\}$
(3)
$\displaystyle\equiv\frac{-i}{16\pi^{2}}\left\\{p_{3}^{\mu}p_{1}^{\nu}~{}\mathcal{B}+\int_{0}^{1}dx\int_{0}^{1-x}dy\left(-\frac{g^{\mu\nu}}{4}\Gamma(\epsilon)\frac{1}{\Delta^{\epsilon}}\right)\right\\}\
,$

and

$\displaystyle\mathcal{B}(t)=$
$\displaystyle\int_{0}^{1}dx\frac{x}{2M^{2}}\left(2(M^{2}(1-2x)+tx)\frac{\mathrm{ArcTan}(\frac{M^{2}-tx}{\Lambda(t,x)})-\mathrm{ArcTan}(\frac{M^{2}(2x-1)-tx}{\Lambda(t,x)})}{\Lambda(t,x)}\right.$
(4) $\displaystyle+\left.\ln\frac{m^{2}+t(x-1)x}{m^{2}+M^{2}(x-1)x}\right)\ .$

In above
$\Lambda(t,x)=\sqrt{4m^{2}M^{2}-M^{4}+4M^{2}tx^{2}-2M^{2}tx-t^{2}x^{2}}$,
$\Delta=M^{2}(x^{2}+y^{2})+(2M^{2}-t)xy-M^{2}(x+y)+m^{2}$,
$\Gamma(\epsilon)\frac{1}{\Delta^{\epsilon}}=\frac{1}{\epsilon}-\ln\Delta-\gamma+\ln
4\pi+\mathcal{O}(\epsilon)$. On the right hand side of Eq. (3), only the
$\mathcal{B}$ term will be considered since the rest will be absorbed by the
contact interactions to be introduced latter. $M$ is the mass of $\psi(3770)$,
$m$ is the mass of $D(\bar{D})$. The parameter $g$ is the coupling strength of
$\psi^{\prime\prime}$ $D\bar{D}$ three point vertex, $g’$ the coupling
strength of $J/\psi J/\psi D\bar{D}$ four point vertex; $t=(p_{2}-p_{4})^{2}$.
The parameter $g$ can be determined by the process $\psi^{\prime\prime}\to
D^{0}\bar{D}^{0}$ decay,

$\displaystyle
i\mathcal{M}_{\psi^{\prime\prime}DD}=i\,g\,\epsilon(\psi^{\prime\prime})\cdot[p(D^{0})-p(\bar{D}^{0})]\
,$ (5)
$\displaystyle\Gamma=\frac{1}{8\pi}\,|i\mathcal{M}_{\psi^{\prime\prime}DD}|^{2}\,\frac{q(DD)}{M_{\psi^{\prime\prime}}^{2}}=\frac{1}{6\pi}\,g^{2}\,\frac{q(DD)^{3}}{M_{\psi^{\prime\prime}}^{2}}\
,$

where $q(DD)$ is the norm of three-dimensional momentum of $D^{0}$ or
$\bar{D}^{0}$ in final state. The PDG value $\Gamma_{\psi^{\prime\prime}\to
D^{0}\bar{D}^{0}}\sim 27.2\times 52\%\times 10^{-3}$ [16] determines $g\sim
12$. 333In Ref. [17], g is estimated to be larger ($\sim 30$). Parameter
$g^{\prime}$ is unknown and left as a free parameter. The amplitude Eq. (4)
contains rather complicated singularity structure, especially the well-known
anomalous threshold, which in history was discovered by Mandelstam who use it
to explain the looseness of deuteron wave function [13]. The anomalous
threshold locates at

$s_{A}=4m^{2}-\frac{(M^{2}-2m^{2})^{2}}{m^{2}}\ .$ (6)

Putting the mass of $\psi^{\prime\prime}$ and $D^{0}$ mesons, one gets
$s_{A}=-1.28$GeV2 (for $D^{+}$ loop, -0.98GeV2 ). Numerically function
$\cal{B}$ is plotted in Fig. 2a, where one clearly sees the anomalous
threshold, beside the normal one at $t=4m^{2}$. Notice that if $M^{2}$ were
smaller than $2m^{2}$, the anomalous branch point locates below physical
threshold, but on the second sheet. It touches the physical threshold $4m^{2}$
and turns up to the physical sheet if the value of $M^{2}$ increased to
$2m^{2}$. When further increase $M^{2}$, the anomalous threshold moves towards
left on the real axis and pass the origin when $M^{2}=4m^{2}$, and finally
reach the physical value, i.e., -1.28GeV2. The situation is depicted in Fig.
2b. Notice that $s_{A}$ here is negative, contrary to what happens in
deuteron, because the latter is a bound state with a normalizable wave
function, whereas $\psi^{\prime\prime}$ is an unstable resonance.

Figure 2: Left: triangle diagram contribution (the y axis label is arbitrary).
Right: the trajectory of the anomalous threshold with respect to the variation
of $M^{2}$.

To proceed, one needs further to make the partial wave projection of
$\mathcal{M}$ and gets

$T^{J}_{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}(s)=\frac{1}{32\pi}\frac{1}{2q^{2}(s)}\int^{0}_{-4q^{2}(s)}\,dt\,d^{J}_{\mu\mu^{\prime}}\left(1+\frac{t}{2q^{2}(s)}\right)\mathcal{M}_{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}(t)\
,$ (7)

where the channel momentum square reads:
$q^{2}(s)=(s-(M-M_{J})^{2})(s-(M+M_{J})^{2})/4s$ with $M_{J}$ the mass of
$J/\psi$, $\mu_{i}$ denotes the corresponding helicity configuration and
$\mu=\mu_{1}-\mu_{2}$, $\mu^{\prime}=\mu_{3}-\mu_{4}$. The key observation is
that the integral interval in Eq. (7) will cover $s_{A}$ if $\sqrt{s}>6.96$GeV
(for $D^{+}$ loop, $\sqrt{s}>6.94$GeV). In other words, the partial wave
amplitude will be enhanced in the vicinity of $X(6900)$ peak by the anomalous
threshold enhancement in the $t$ channel, see Fig. 3 for illustration.

Figure 3: The enhancement of the $J/\psi\psi^{\prime\prime}$ scattering
amplitude from the triangle diagram.

Based on the above observation it is suggested that the $X(6900)$ peak may at
least partly be explained by the anomalous threshold generated by the triangle
diagram as depicted in Fig. 1. Furthermore, to get the $L=0$ (s wave)
amplitudes, we need the relation between s wave and helicity amplitudes (see
Refs. [11, 12] for more discussions):

$T_{L=0}^{0,ij}(s)=\frac{1}{3}\left[2T_{++++}^{0,ij}(s)+2T_{++--}^{0,ij}(s)-2T_{++00}^{0,ij}(s)-2T_{00++}^{0,ij}(s)+T_{0000}^{0,ij}(s)\right],$
(8) $\displaystyle T_{L=0}^{2,ij}(s)$
$\displaystyle=\frac{1}{15}[T_{++++}^{2,ij}(s)+T_{++--}^{2,ij}(s)]$ (9)
$\displaystyle+\frac{\sqrt{6}}{15}[T_{+++-}^{2,ij}(s)+T_{+-++}^{2,ij}(s)+T_{++-+}^{2,ij}(s)+T_{-+++}^{2,ij}(s)]$
$\displaystyle+\frac{\sqrt{3}}{15}[T_{+++0}^{2,ij}(s)+T_{+0++}^{2,ij}(s)+T_{++0+}^{2,ij}(s)+T_{0+++}^{2,ij}(s)$
$\displaystyle+T_{++-0}^{2,ij}(s)+T_{-0++}^{2,ij}(s)+T_{++0-}^{2,ij}(s)+T_{0-++}^{2,ij}(s)]$
$\displaystyle+\frac{1}{5}[T_{+00+}^{2,ij}(s)+T_{0++0}^{2,ij}(s)+T_{+00-}^{2,ij}(s)+T_{0-+0}^{2,ij}(s)+T_{+0+0}^{2,ij}(s)$
$\displaystyle+T_{0+0+}^{2,ij}(s)+T_{0+0-}^{2,ij}(s)+T_{+0-0}^{2,ij}(s)]$
$\displaystyle+\frac{\sqrt{2}}{5}[T_{+-+0}^{2,ij}(s)+T_{+0+-}^{2,ij}(s)+T_{+-0+}^{2,ij}(s)+T_{0++-}^{2,ij}(s)$
$\displaystyle+T_{-++0}^{2,ij}(s)+T_{+0-+}^{2,ij}(s)+T_{-+0+}^{2,ij}(s)+T_{0+-+}^{2,ij}(s)]$
$\displaystyle+\frac{2}{15}[T_{++00}^{2,ij}(s)+T_{00++}^{2,ij}(s)+T_{0000}^{2,ij}(s)]+\frac{2\sqrt{6}}{15}[T_{+-00}^{2,ij}(s)+T_{00+-}^{2,ij}(s)]$
$\displaystyle+\frac{2}{5}[T_{+-+-}^{2,ij}(s)+T_{+--+}^{2,ij}(s)]+\frac{2\sqrt{3}}{15}[T_{+000}^{2,ij}(s)+T_{00+0}^{2,ij}(s)+T_{0+00}^{2,ij}(s)+T_{000+}^{2,ij}(s)]\
,$

where $i,j=1,2,3$ denote the corresponding outgoing and incoming channels
(corresponding to $J/\psi J/\psi$, $J/\psi\psi(3685)$ and $J/\psi\psi(3770)$,
respectively). In practice, it is found that the anomalous enhancement gives a
more prominent effect to the $J=2$ amplitude than the $J=0$ amplitude. The
amplitude $T^{J,ij}_{L=0}$ is a $K$\- matrix unitarized amplitude obtained
from the tree level amplitudes generated from the following contact
interaction lagrangian:

$\displaystyle\mathcal{L}_{c}=$ $\displaystyle
c_{1}V_{\mu}V_{\alpha}V^{\mu}V^{\alpha}+c_{2}V_{\mu}V_{\alpha}V^{\mu}V^{\prime\alpha}+c_{3}V_{\mu}V_{\alpha}^{\prime}V^{\mu}V^{\prime\alpha}+c_{4}V_{\mu}V^{\prime\mu}V_{\alpha}V^{\prime\alpha}+c_{5}V_{\mu}V_{\alpha}V^{\mu}V^{\prime\prime\alpha}$
(10)
$\displaystyle+c_{6}V_{\mu}V_{\alpha}^{\prime\prime}V^{\mu}V^{\prime\prime\alpha}+c_{7}V_{\mu}V^{\prime\prime\mu}V_{\alpha}V^{\prime\prime\alpha}+c_{8}V_{\mu}V_{\alpha}^{\prime}V^{\mu}V^{\prime\prime\alpha}+c_{9}V_{\mu}V^{\prime\mu}V_{\alpha}V^{\prime\prime\alpha},$

The production amplitudes $F^{J}_{i}$ are parameterized as,

$F^{J}_{i}(s)=\sum^{3}_{k=1}\alpha_{k}(s)T_{L}^{J,ki}(s)~{},$ (11)

where $\alpha_{k}(s)$ are real polynomial functions in general and are set to
be constants here, and further we set $\alpha_{1}(s)^{2}=1$. Further,

$\frac{dEvents_{i}}{d\sqrt{s}}=N_{i}~{}p_{i}(s)~{}|F_{i}|^{2}~{},$ (12)

where $p_{i}(s)$ refers to the abs of three momentum for corresponding
channel. According to partial wave convention, for $J=0$ and $J=2$ case, they
have a total scale factor that [11, 12]

$|F_{i}|^{2}=|F_{i}^{J=0}|^{2}+5~{}|F_{i}^{J=2}|^{2}\ .$ (13)

The fit is overdone since there are many parameters. One solution is shown in
Fig. 4, and the fit parameters are listed in Tab. 1 for illustration. During
the fit many solutions are examined, nevertheless it is found that the
triangle diagram contributions are all small. It is not totally clear to us
why it behaves like such but one possible reason could be that the peak
position through anomalous threshold contribution, as shown in Fig. 3, is
roughly 60 – 80MeV above the $X(6900)$ peak, and hence the fit does not like
it. One possible way to rescue the problem may be to adopt another
parameterization in which the background contributions are more flexible to be
tuned, hence the interference between the background and the anomalous
enhancement can lead to the shift of the peak position by a few tens of MeV.
We leave this investigation for future studies.

Parameter | $\chi^{2}/d.o.f$ | $g^{\prime}$ | $N_{1}$ | $N_{2}$ | $\alpha_{2}$ | $\alpha_{3}$
---|---|---|---|---|---|---
Fit | $0.93$ | $-18.3\pm 1.6$ | $23\pm 10$ | $0.34\pm 0.23$ | $-4.94\pm 0.02$ | $3.97\pm 0.05$

$c_{5}$ | $c_{6}$ | $c_{7}$ | $c_{8}$ | $c_{9}$
---|---|---|---|---
$11.34\pm 0.05$ | $50.79\pm 0.05$ | $-35.01\pm 0.19$ | $-64.71\pm 0.07$ | $1.529\pm 0.002$

Table 1: The fit parameters of Fig. 4. The $c_{i}$ parameters are defined in
Eq. (10). In this fit, we set $c_{1},\cdots,c_{4}$ to be negligible. Figure 4:
The fit results from table 1.

## Acknowledgements

We would like to thank De-Liang Yao and Ling-Yun Dai for very helpful
discussions. This work is supported in part by National Nature Science
Foundations of China under Contract Number 12335002.

## References

  * [1] R. Aaij et al. (LHCb Collaboration), Sci.Bull. 65 (2020) 23, 1983.
  * [2] CMS collaboration, CMS-PAS-BPH-21-003 (2022).
  * [3] ATLAS collaboration, ATLAS-CONF-2022-040 (2022).
  * [4] Y. Lu et al., Phys. Rev. D107 (2023) 094006.
  * [5] Q. F. Cao $et$ $al.$, Chin. Phys. C 45 (2021) 10, 103102.
  * [6] D. Morgan, Nucl. Phys. A 543 (1992) 632.
  * [7] O. Zhang, C. Meng and H. Q. Zheng, Phys. Lett. B 680 (2009) 453.
  * [8] L. Y. Dai, M. Shi, G. Y. Tang and H. Q. Zheng, Phys. Rev. D 92 (2015) 014020.
  * [9] Q. R. Gong, Z. H. Guo, C. Meng, G. Y. Tang, Y. F. Wang and H. Q. Zheng, Phys. Rev. D 94 (2016) 114019.
  * [10] Q. F. Cao, H. R. Qi, Y. F. Wang and H. Q. Zheng, Phys. Rev. D 100 (2019) 054040.
  * [11] Z. R. Liang, X. Y. Wu, D. L. Yao, Phys. Rev. D104 (2021) 034034.
  * [12] Q. Zhou, D. Guo, S. Q. Kuang, Q. H. Yang and L. Y. Dai, Phys. Rev. D 106 (2022) 11.
  * [13] S. Mandelstam, Phys. Rev. Lett. 4 (1960) 84.
  * [14] E. Braaten, H. W. Hammer, Phys. Rep. 428 (2006) 259.
  * [15] H. W. Hammer, C. Ji and D. R. Phillips, J. Phys. G44 (2017) 103002.
  * [16] M. Tanabashi et al. (Particle Data Group), Phys. Rev. D 98 (2018) 030001.
  * [17] S. Coito and F. Giacosa, Nucl. Phys. A 981 (2019) 38.


11institutetext: Instituto de Física de Cantabria, CSIC-UC, Av. de Los Castros
s/n, E-39005 Santander, Spain 22institutetext: Departamento de Matemáticas,
Universidad de Oviedo, C. Federico García Lorca 18, 33007 Oviedo, Spain
33institutetext: Departamento de Física, Universidad de Oviedo, C. Federico
García Lorca 18, 33007 Oviedo, Spain 44institutetext: Instituto Universitario
de Ciencias y Tecnologías Espaciales de Asturias (ICTEA), Escuela de
Ingeniería de Minas, Materiales y Energía de Oviedo, C. Independencia 13,
33004 Oviedo, Spain 55institutetext: Departamento de Física Moderna,
Universidad de Cantabria, 39005-Santander, Spain 66institutetext: ESAC, Camino
Bajo del Castillo s/n, 28692 Villafranca del Castillo, Madrid, Spain

# A Bayesian method for point source polarization estimation

D. Herranz 11 F. Argüeso 2244 L. Toffolatti 3344 A. Manjón-García 1155 M.
López-Caniego 66

The estimation of the polarization $P$ of extragalactic compact sources in
Cosmic Microwave Background images is a very important task in order to clean
these images for cosmological purposes –as, for example, to constrain the
tensor-to-scalar ratio of primordial fluctuations during inflation– and also
to obtain relevant astrophysical information about the compact sources
themselves in a frequency range, $\nu\sim 10$–$200$ GHz, where observations
have only very recently started to be available. In this paper we propose a
Bayesian maximum a posteriori (MAP) approach estimation scheme which
incorporates prior information about the distribution of the polarization
fraction of extragalactic compact sources between 1 and 100 GHz. We apply this
Bayesian scheme to white noise simulations and to more realistic simulations
that include CMB intensity, Galactic foregrounds and instrumental noise with
the characteristics of the QUIJOTE experiment Wide Survey at 11 GHz. Using
these simulations, we also compare our Bayesian method with the frequentist
Filtered Fusion method that has been already used in WMAP data and in the
_Planck_ mission. We find that the Bayesian method allows us to decrease the
threshold for a feasible estimation of $P$ to levels below $\sim 100$ mJy (as
compared to $\sim 500$ mJy that was the equivalent threshold for the
frequentist Filtered Fusion). We compare the bias introduced by the Bayesian
method and find it to be small in absolute terms. Finally, we test the
robustness of the Bayesian estimator against uncertainties in the prior and in
the flux density of the sources. We find that the Bayesian estimator is robust
against moderate changes in the parameters of the prior and almost insensitive
to realistic errors in the estimated photometry of the sources.

###### Key Words.:

methods: data analysis; techniques: image processing; cosmic background
radiation; radio continuum: galaxies; polarization

##  1 Introduction

The polarization properties of extragalactic radio sources (ERS) –i.e., radio
galaxies, radio loud quasars, blazars, etc.– are not well constrained even at
cm wavelengths, given that the total linear polarization of ERS, $P$, in
general constitutes a small fraction of their total flux density, $S$. The
observed value of $P$ being typically a few per cent, with only very few ERS
showing a total polarization fraction, $\Pi=P/S$, as high as $\sim 10$ per
cent of the total flux density (e.g. Sajina et al. 2011; Tucci & Toffolatti
2012). Moreover, at shorter wavelengths, i.e. at $\lambda\leq 1$ cm, these
properties are still poorly known, due to the difficulty to properly calibrate
in the radio to mm regime that afflicted the polarization experiments until a
few years ago. However, the knowledge of the total and polarization fraction
of ERS is rapidly improving at high radio frequencies thanks to large samples
of sources mainly observed by the Australia Telescope Compact Array (ATCA) and
by the Very Large Array (VLA) (Sadler et al. 2006; López-Caniego et al. 2009;
Massardi et al. 2008, 2011, 2013; Murphy et al. 2010; Jackson et al. 2010;
Galluzzi et al. 2017, 2018). More recently, and thanks to the very high
sensitivity of the new detectors of the Atacama Large Millimeter Array (ALMA),
Galluzzi et al. (2019) could extend up to 97.5 GHz the analysis of
polarization properties of ERS performed by Galluzzi et al. (2018), by
polarimetric observations of a complete sample of 32 extragalactic radio
sources. Their findings showed that the distribution of the observed $\Pi$
fractions is, again, well fitted by a log-normal distribution, thus confirming
previous outcomes at lower frequencies (Massardi et al. 2013; Galluzzi et al.
2018) and also the predictions of Tucci & Toffolatti (2012). The analysis of
Galluzzi et al. (2019) also confirmed the absence of any statistically
significant trend of polarization properties of ERS with the frequency or the
flux density.

Recent analyses of the ERS present in the full sky cosmic microwave background
(CMB) anisotropy maps in polarization provided by the European Space Agency
(ESA) Planck mission (Planck Collaboration I 2016; Planck Collaboration XXVI
2016) also indicate typical median polarization fractions of ERS of $2-3$% at
frequencies as high as 300 GHz (Bonavera et al. 2017a, b; Trombetti et al.
2017). Therefore, an accurate characterization of polarization properties of
ERS as well as their efficient detection and subtraction from CMB maps is
especially crucial for measuring the primordial CMB B-mode polarization down
to values of the tensor to scalar ratios $r\sim 0.001$, that could be
achievable by future space probes (i.e. Sekimoto et al. 2018, LiteBird:;
Delabrouille et al. 2018, COrE:). We remind that the simulations by
Remazeilles et al. (2018) have shown that, at these low values of $r$,
unresolved polarized ERS can probably be the dominant foreground at multipoles
$\ell>50$ in the power spectrum of the CMB anisotropy. These results have been
confirmed by Puglisi et al. (2018), by exploiting the state-of-the-art data
sets on polarized point sources over the 1.4–217 GHz frequency range.

In addition to the essential information that the polarization of ERS provides
about the structure and evolution of extragalactic baryonic matter at low to
intermediate redshifts, the study of this polarized radiation is paramount for
cosmology, and in particular for Cosmic Microwave Background (CMB) science.
ERS detection and subtraction is a fundamental part of the component
separation process necessary to achieve the science goals set for the next
generation of CMB experiments. In particular, ERS would significantly affect
the estimation of the CMB polarization angular power spectra and, therefore,
limit the ability of CMB experiments to constrain cosmological parameters such
as the tensor-to-scalar ratio $r$ of primordial perturbations during
inflation. ERS could become an important obstacle for the detection of the
Primordial Gravitational Wave Background (PGWB) for low values of $r$ (Tucci
et al. 2005; Puglisi et al. 2018; Trombetti et al. 2018) due to both the
additional noise they constitute in themselves, and the reduction in delensing
power they cause by degrading lensing potential reconstructions (see, e.g.,
Sailer et al. 2020). Therefore, during recent years the interest in the
development of signal processing techniques specifically tailored for the
detection and characterization of ERS in CMB images has been growing in the
literature.

Signal processing techniques for the detection of polarized ERS must take into
account the spinorial nature of electromagnetic waves. The signal can be
described by not only one but as many as four independent components, one for
the total intensity of the radiation field and three for its polarization
state. It is convenient to use the Stokes’ parameters $S,Q,U,V$ ($S$ for total
intensity111The usual notation for this Stokes parameter is $I$. However, in
this work we have changed the notation in order to avoid confusion between the
intrinsic intensity of a source, that we will call later in this paper
$S_{0}$, and the modified Bessel function of zero order, $I_{0}$, that appears
in several equations in Section 2. in terms of flux density, see, e.g.,
Galluzzi et al. (2019), $Q$ and $U$ for linear polarization and $V$ for
circular polarization), but other representations are also possible. The
Stokes’ $V$ parameter is not usually considered, since Thomson scattering does
not induce circular polarization in the CMB. Circular polarization mechanisms
in active galaxies have been described in the literature (see for example
Rayner et al. 2000), but they are nonetheless considered to be sub-dominant in
comparison to linear polarization mechanisms. Therefore, in this paper we will
consider, as it is customary222Foregrounds can produce circular polarization
under some circumstances, and it has been observed in a few extragalactic
sources. The value of $V$ is typically much lower than the other Stokes’
parameters. As it will be explained in Section 2, the existence of sources
with non-zero circular polarization would not affect our estimations of the
$Q$ and $U$ Stokes parameters. Of course, if there was a significant $V$ term,
neglecting it would lead us to miss a part of the polarization $P$. However,
our method can be easily adapted to work with a third component in the form of
an additional image –corresponding to the $V$ Stokes’ parameter– if
necessary., $V=0$. Then the signal processing of polarized ERS must deal with
three independent quantities, two of them having the mathematical structure of
a spinor field.

The $S$, $Q$, and $U$ signals (or, alternatively, $S$, $E$, and $B$, or any
other set of three quantities obtained from the Stokes’ parameters) can be
treated separately as independent images to which any of the standard compact
component separation techniques could be applied. The main difference with
respect to the classical setting is that, unlike the total flux density $S$,
which is always non-negative, $Q$ and $U$ can be either positive, negative or
zero. From a physical point of view, however, it makes more sense to process
the polarization data jointly (see Herranz et al. 2012, for a review on the
topic). In particular, the total polarization of a source
$P=\sqrt{Q^{2}+U^{2}}$ and its polarization fraction $\Pi=P/S$ are directly
related to the physical processes occurring along the path of photons from the
ERS to Earth, while $Q$ and $U$ are frame-dependent quantities lacking in
physical meaning on themselves.

The main two problems arising when dealing with $P$ are the typically low
signal-to-noise ratio of the polarization signal coming from ERS and the non-
Gaussian distribution of its noise statistics. Regarding the former, as
mentioned above the typical polarization fractions of ERS at frequencies below
$\sim 10$ GHz are at most $10\%$. This means that only a few ERS are bright
enough to be detected in polarization with present-day technology. A standard
procedure to avoid false detections in polarization is to _detect_ sources in
total intensity and then to try to estimate their polarization properties in a
non-blind way333That is, focusing efforts on the precise position of the
source once it has been detected in intensity, i.e. the non-blindness is only
related to the positions of targets, not to any other quantity.. We will
follow this approach in this paper. Regarding the latter problem, assuming
that the $Q$ and $U$ noises are Gaussian-distributed, $P$ will have a non-
Gaussian Rice distribution (Rice 1945). Rician distribution has strictly non-
negative support and heavy tails, which a) biases the estimation of the
polarization of the sources and b) disrupts the intuitive interpretation of
signal-to-noise in terms of $\sigma$ thresholds that is used virtually
everywhere else in radio Astronomy. Simmons & Stewart (1985) discussed four
estimators which attempted to correct for biasing in the degree of linear
polarization in the presence of low signal-to-noise ratios. More recently,
Argüeso et al. (2009) studied the problem in the context of CMB astronomy and
developed two methods for the detection/estimation of ERS in polarization
data: one that applies the Neyman-Pearson lemma to the Rice distribution, the
Neyman-Pearson filter (NPF), and another based on pre-filtering before fusion
of $Q$ and $U$ to obtain $P$, the filtered fusion (FF) method. That work found
that under typical CMB-experiment settings the FF outperforms the NPF both in
terms of computational simplicity and accuracy, especially for low fluxes.
López-Caniego et al. (2009) applied the FF to the WMAP five-year data. The
same method has been used to study the polarization of the _Planck_ Second
Catalogue of Compact Sources (PCCS2, Planck Collaboration XXVI 2016) and of
the QUIJOTE experiment Wide Survey Source Catalogue (Herranz et al. 2020).
Alternatively, a novel method for the estimation of the polarization intensity
and angle of compact sources in the $E$ and $B$ modes of polarization based on
steerable wavelets has been recently proposed by Diego-Palazuelos et al.
(2020).

All the previously mentioned methods attempt to estimate the ERS polarization
by minimizing as much as possible the impact of noise and Galactic and
extragalactic foregrounds on the observed signal. The _expected_ value of the
polarization does not intervene in the estimation process. In other words, no
_a priori_ information is used in the estimation. Until very recently, this
has been the most sensible choice, as the polarization properties of
extragalactic sources were virtually unknown at microwave frequencies.
However, as recent experiments and facilities such as the ALMA, _Planck_ and
the upgraded versions of ATCA and VLA start shedding light on the $\lambda\leq
1$ cm polarized sky, the possibility of adding physical priors to our signal
processing techniques is gradually opening. In this paper, we propose a
Bayesian maximum _a posteriori_ (MAP) method for the estimation of the
polarization properties of point sources.

The structure of this paper is as follows. In section 2 we review the current
observational evidence to construct physical priors on the polarization
fraction of ERS and incorporate that information into two possible MAP
estimators of the polarization of a compact source of known flux density $S$.
These MAP estimators take a form analogous to the Neyman-Pearson filter and
Filtered Fusion by Argüeso et al. (2009), respectively, plus additional terms
that contain the _a priori_ physical information of the probability
distribution function of $P$ for ERS. We call these two methods Bayesian Rice
and Bayesian Filtered Fusion, respectively. The Bayesian Filtered Fusion is
easily applicable for both white and colour noise. For this reason, and
because in Argüeso et al. (2009) it was shown that the FF outperforms the NPF,
we focus the rest of the paper on the Bayesian Filtered Fusion. In section 3
we describe the simulations we have used to test the Bayesian Filtered Fusion
method. We first make simplistic simulations containing just white noise in
$Q$ and $U$ and then we upgrade to realistic simulations with polarized
Galactic foregrounds and CMB emissions. In both cases, we use angular
resolution, pixel scale and noise levels similar to the upcoming QUIJOTE
experiment Wide Survey data at 11 GHz (Rubiño-Martín et al. 2020). The results
of applying the Bayesian Filtered Fusion to our simulations are discussed in
section 4, where we also make a brief discussion about the robustness of the
method against uncertainties on the priors and the determination of the total
flux density of the sources. Finally, we draw our conclusions in section 5.

##  2 Method

When we try to detect or estimate the polarization $P_{0}$ of a compact source
embedded in Gaussian noise –what can be a good approximation when dealing with
sources present in CMB maps–, and we consider the measured polarization
$P=\sqrt{Q^{2}+U^{2}}$ with similar Gaussian noise dispersions in $Q$ and $U$,
i.e. $\sigma_{Q}=\sigma_{U}=\sigma$, the distribution of $P$ given $P_{0}$
follows the Rice distribution

$f\left(P|P_{0}\right)=P/\sigma^{2}\exp\left[-\left(P^{2}+P_{0}^{2}\right)/2\sigma^{2}\right]I_{0}\left(PP_{0}/\sigma^{2}\right)$
(1)

This is the conditional probability distribution of $P$, measured
polarization, given $P_{0}$, source polarization, with $I_{0}$ the modified
Bessel function of zero order (Rice 1945). This distribution has been used to
obtain suitable estimators of $P_{0}$ (Simmons & Stewart 1985; Argüeso et al.
2009; López-Caniego et al. 2009; Herranz et al. 2012). However, no previous
knowledge about $P_{0}$ is assumed in these papers and this information, if
available, could be very useful after being incorporated in a Bayesian scheme.
Recently, relevant data about the distribution of the polarization fraction at
different frequencies, $\Pi=P_{0}/S_{0}$ , with $S_{0}$ the source flux
density, have been presented (Massardi et al. 2013; Galluzzi et al. 2017,
2019). According to these authors, this distribution can be represented by a
log-normal probability density function (pdf)

$g(\Pi)=\frac{1}{\sqrt{2\pi}\sigma_{\Pi}\Pi}\exp\left[-(\log(\Pi)-\mu)^{2}/2\sigma_{\Pi}^{2}\right]$
(2)

where $\mu=\log\Pi_{med}$ –with $\Pi_{med}$ the median polarization fraction–
and $\sigma_{\Pi}$ are easily obtained from $\langle\Pi\rangle$ and
$\langle\Pi^{2}\rangle$ (Crow & Shimizu 1988). Since $P_{0}=\Pi S_{0}$, if we
assume that the value of $S_{0}$ is known, then the distribution of $P_{0}$ is
readily calculated

$g\left(P_{0}\right)=\frac{1}{\sqrt{2\pi}\sigma_{\Pi}P_{0}}\exp\left[-\left(\log
P_{0}-\mu_{1}\right)^{2}/2\sigma_{\Pi}^{2}\right]$ (3)

with $\mu_{1}=\mu+\log\left(S_{0}\right)=\log\left(S_{0}\Pi_{med}\right)$.
This assumption is safe, since estimating $S_{0}$ is much simpler than
estimating $P_{0}$ and point sources for which the polarized emission is
detectable tend to have high flux densities. Therefore, we will work, in
general, with non-blind detection in polarization. The knowledge of $S_{0}$
allows us to write the joint probability distribution of $P$ and $P_{0}$,
$h(P,P_{0})=f(P|P_{0})g(P_{0})$, and by using Bayes’ theorem

$f\left(P_{0}|P\right)=\frac{h\left(P,P_{0}\right)}{g(P)}$ (4)

with $g(P)=\int h\left(P,P_{0}\right)dP_{0}$. Finally, we obtain, by
substituting (1) and (3) in (4),

$f\left(P_{0}|P\right)=\frac{e^{-P_{0}^{2}/2\sigma^{2}}I_{0}\left[\frac{PP_{0}}{\sigma^{2}}\right]\exp\left[-\frac{\left(\log
P_{0}-\mu_{1}\right)^{2}}{2\sigma_{\Pi}^{2}}\right]\frac{1}{P_{0}}}{\int{e^{-P_{0}^{2}/2\sigma^{2}}I_{0}\left[\frac{PP_{0}}{\sigma^{2}}\right]\exp\left[-\frac{\left(\log
P_{0}-\mu_{1}\right)^{2}}{2\sigma_{\Pi}^{2}}\right]\frac{dP_{0}}{P_{0}}}}.$
(5)

The integral in the denominator is just a normalization. We have found the
distribution of $P_{0}$ given $P$, the posterior distribution, simply by
assuming Gaussian noise in $Q$ and $U$ with the same dispersion and a log-
normal pdf for $\Pi$ (prior distribution). Everything has been calculated for
a source located at a central pixel and without taking into account any
information about the beam and the data in a certain patch around the source.
If we consider a polarized source at the central pixel of an $n$-pixel patch,
a beam with profile $\tau(\mathbf{x})$ and values $P_{i}$ for the polarization
measured at each pixel $i=1,\ldots,n$, we can write the following expression
for the conditional pdf of $P_{0}$ given the values $P_{i}$ at the different
pixels:

$\displaystyle f\left(P_{0}|P_{i}\right)$ $\displaystyle\propto$
$\displaystyle\exp\left[-\frac{\left(\log
P_{0}-\mu_{1}\right)^{2}}{2\sigma_{\Pi}^{2}}\right]\frac{1}{P_{0}}\times$ (6)
$\displaystyle\prod_{i}\exp\left[-P_{0}^{2}\frac{\tau_{i}^{2}}{2\sigma_{i}^{2}}\right]I_{0}\left[P_{0}\frac{P_{i}\tau_{i}}{\sigma_{i}^{2}}\right]$

Here $\Pi_{i}$ is the product symbol, $\tau_{i}$ is the profile at each pixel
and $\sigma_{i}$ the noise dispersion (it could be different from pixel to
pixel). We take the natural logarithm of the right-hand side and change sign,
this is minus the log-posterior of the distribution, save constant terms. In
this way, we obtain a simplified expression that we will later minimize to
find the estimator that makes the posterior distribution maximum

$\displaystyle-\log f\left(P_{0}|P_{i}\right)$ $\displaystyle=$
$\displaystyle\frac{\left(\log
P_{0}-\mu_{1}\right)^{2}}{2\sigma_{\Pi}^{2}}+\log\left(P_{0}\right)+P_{0}^{2}z$
(7) $\displaystyle-$ $\displaystyle\sum_{i}\log
I_{0}\left[P_{0}\,y_{i}\right]+K$

with

$z=\sum_{i}\frac{\tau_{i}^{2}}{2\sigma_{i}^{2}}$ (8)

and

$y_{i}=\frac{P_{i}\tau_{i}}{\sigma_{i}^{2}},$ (9)

and where $K$ is a constant term that encloses the proportionality terms not
directly included in (7). If we differenciate with respect to $P$ and equate
to zero, the estimator $\hat{P}_{0}$ will satisfy

$\frac{\log\hat{P}_{0}-\mu_{1}}{\hat{P}_{0}\sigma_{\Pi}^{2}}+\frac{1}{\hat{P}_{0}}+2\hat{P}_{0}z-\sum_{i}\frac{I_{1}\left[\hat{P}_{0}y_{i}\right]}{I_{0}\left[\hat{P}_{0}y_{i}\right]}\,y_{i}=0$
(10)

with $I_{1}$ the modified Bessel function of order one. On the other hand, in
Argüeso et al. (2009) a method called filtered fusion (FF) was shown to
perform better than the one derived from the Rice distribution. The FF
calculates the square root of the sum of the squares of the maps in Q and U to
which a matched filter has been previously applied. This method is just a
maximization of the conditional probability of the data $Q_{i}$, $U_{i}$ given
the source polarization $Q_{0}$, $U_{0}$, assuming that the noise is Gaussian
with zero mean and independent for each pixel

$f(Q_{i},U_{i}|Q_{0},U_{0})=\prod_{i}\exp\left[\frac{\left(Q_{i}-Q_{0}\tau_{i}\right)^{2}}{2\sigma_{Q_{i}}^{2}}+\,\frac{\left(U_{i}-U_{0}\tau_{i}\right)^{2}}{2\sigma_{U_{i}}^{2}}\right]$
(11)

If we also assume that $P_{0}$ follows a log-normal pdf and the polarization
angle distribution is uniform, we can combine the previous formula with the
prior distribution of $Q_{0},U_{0}$ and write, by applying Bayes’ theorem,
minus the log-posterior of $Q_{0}$, $U_{0}$ given the data $Q_{i}$, $U_{i}$

$\displaystyle-\log f\left(Q_{0},U_{0}|Q_{i},U_{i}\right)$ $\displaystyle=$
$\displaystyle\frac{\left(\log\sqrt{Q_{0}^{2}+U_{0}^{2}}-\mu_{1}\right)^{2}}{2\sigma_{\Pi}^{2}}$
(12) $\displaystyle+$ $\displaystyle\log\left(Q_{0}^{2}+U_{0}^{2}\right)$
$\displaystyle+$
$\displaystyle\sum_{i}\frac{\left(Q_{i}-Q_{0}\tau_{i}\right)^{2}}{2\sigma_{Q_{i}}^{2}}$
$\displaystyle+$
$\displaystyle\sum_{i}\frac{\left(U_{i}-U_{0}\tau_{i}\right)^{2}}{2\sigma_{U_{i}}^{2}}+K.$

This expression could be very easily generalized, allowing even the treatment
of correlations between the noise in different pixels. In that case, (12) can
be written

$\displaystyle-$ $\displaystyle\log f\left(Q_{0},U_{0}|Q_{i},U_{i}\right)=$
(13)
$\displaystyle\frac{\left(\log\sqrt{Q_{0}^{2}+U_{0}^{2}}-\mu_{1}\right)^{2}}{2\sigma_{\Pi}^{2}}$
$\displaystyle+\log\left(Q_{0}^{2}+U_{0}^{2}\right)$
$\displaystyle+\frac{1}{2}\sum_{\mathbf{k}}\left(Q_{\mathbf{k}}-Q_{0}\tau_{\mathbf{k}}\right)^{t}\mathbf{P}^{-1}_{Q,\mathbf{k}}\left(Q_{\mathbf{k}}-Q_{0}\tau_{\mathbf{k}}\right)$
$\displaystyle+\frac{1}{2}\sum_{\mathbf{k}}\left(U_{\mathbf{k}}-U_{0}\tau_{\mathbf{k}}\right)^{t}\mathbf{P}^{-1}_{U,\mathbf{k}}\left(U_{\mathbf{k}}-U_{0}\tau_{\mathbf{k}}\right)$
$\displaystyle+K,$

where the subindex $\mathbf{k}$ refers to the Fourier wave vector (or, for the
case of spherical data, the appropriate spherical harmonic) and
$P_{x,\mathbf{k}}$ is the power spectrum (or angular power spectrum) of the
noise for the Stokes parameter $x$. This formula is expressed in Fourier space
for the sake of computational efficiency, but could be also expressed in real
space by means of the correlation matrix of the noise. In (13) it is immediate
to recognize that the third and fourth terms in the right side of the equation
are analogous to the matched filter on the $Q$ and $U$ maps, that in turn are
the solution of a Maximum Likelihood estimator (MLE). The first term adds the
prior information, whereas the second term acts as a penalty for large values
of the estimated polarization, and finally the last term is a constant that is
irrelevant for the solution.

Now we can obtain the estimators $\hat{Q}_{0}$ and $\hat{U}_{0}$ that minimize
the previous expression. Finally, we find

$\hat{P}_{0}=\sqrt{\hat{Q}^{2}_{0}+\hat{U}^{2}_{0}}$ (14)

as our estimator of $P_{0}$. In all these formulas we have not considered the
effect of the circular polarization, $V$. As commented in the introduction,
this effect is very small and can be, in general, neglected. At any rate, the
generalization of the Rice and FF methods to include circular polarization has
been presented in Argüeso et al. (2011). However, in order to extend our
Bayesian methods to this case, we would have to use a prior distribution for
$V$ which is not known yet.

To sum up the previous paragraphs: we have presented four possible estimators
of $P_{0}$, the old ones are the Rice method and the FF (Argüeso et al. 2009;
López-Caniego et al. 2009; Herranz et al. 2012; Planck Collaboration XXVI
2016), and the new ones, obtained by minimizing the right-hand sides of (7)
and (12) or (13), that we will call Bayesian Rice method and Bayesian FF,
respectively. These new methods incorporate in a natural way our information
about the source polarization distribution. The FF and Rice methods are
implemented by minimizing (7) and (12)-(13) without the first two terms, which
come from the Bayesian prior.

##  3 Simulations

###  3.1 White noise

As a first test bed to assess the performance of the Bayesian techniques, we
have run 10000 simulations using only white noise as background for simulated
sources. Table 1 shows the simulation parameters for these simulations; the
pixel size, beam FWHM and white noise rms emulate those of the QUIJOTE (Q-U-I
JOint TEnerife, Rubiño-Martín et al. 2010, 2012; Génova-Santos et al. 2015;
López-Caniego 2016) experiment Wide Survey at 11 GHz (Rubiño-Martín et al.
2020). We will use the same simulation parameters for the full sky simulations
to be discussed in Section 3.2.

For these white noise simulations we directly create flat images with
uncorrelated Gaussian noise and inject at the center a point source with he
FWHM listed in Table 1, a given flux density $S_{0}$ and polarization fraction
$\Pi$ randomly drawn from the log-normal distribution (2) with the mean and
standard deviation values $\langle\Pi\rangle$ and $\sigma_{\Pi}$ described in
Table 1. We simulate intensities $S_{0}$ in ten logarithmically spaced values
between 0.1 and 100 Jy. In this way we get a sample of sources from moderately
faint to extra bright (and, since the polarization fraction follows
distribution (2) with average $\langle\Pi\rangle=0.02$, values of $P$ from
below 1 mJy to a few tens of Janskys).

###  3.2 Full sky simulations

In order to assess the performance of our Bayesian techniques under realistic
conditions, we use realistic simulations of the QUIJOTE experiment Wide Survey
(Rubiño-Martín et al. 2020). The QUIJOTE Wide Survey is observing
approximately half the sky at 11, 13, 17 and 19 GHz. These simulations have
been produced thanks to the EU RADIOFOREGROUNDS project444The RADIOFOREGROUNDS
project aims to combine two unique datasets, the nine Planck all-sky (30-857
GHz) maps and the four QUIJOTE Northern sky (10-20 GHz) maps, to provide the
best possible characterization of the physical properties of polarized
emissions in the microwave domain, together with an unprecedentedly thorough
description of the intensity signal. This legacy information will be essential
for future sub-orbital or satellite experiments. See more information in
http://www.radioforegrounds.eu/, but are not yet public. The simulations make
use of the Planck Sky Model (PSM, Delabrouille et al. 2013), a global
representation of the multi-component sky at frequencies ranging from 1 GHz to
1 THz that summarizes in a synthetic way as much as possible of our present
knowledge of the GHz sky. The PSM is a public code555http://www.apc.univ-
paris7.fr/$\sim$delabrou/PSM/psm.html developed mostly by members of the
Planck collaboration as a simulation tool for Planck activities and it makes
possible to simulate the sky in total intensity and the $Q$, $U$ Stokes
parameters for any experimental configuration in the GHz range. For this work,
we choose to simulate the QUIJOTE wide Survey at 11 GHz. Figures 1 and 2 show
the full $Q$ Stokes Wide Survey simulated at 11 GHz and the simulated QUIJOTE
instrumental noise for the same Stokes parameter, area and frequency. Table 1
indicates the main parameters used for this simulation.

Frequency (GHz) | 11
---|---
FWHM (degrees) | 0.85
Nside parameter | 246
Pixel resolution (arcmin) | 13.74
White noise rms (Jy) | 0.386
$\langle\Pi\rangle$ | 0.02
$\sigma_{\Pi}$ | 1.0
Table 1: Simulation parameters used for this work.

Formula (13) could be applied to the whole sky, but since the statistical
properties of the foregrounds vary strongly with Galactic latitude we prefer
to apply the Bayesian estimator locally. In order to test the method, we
compute (13) on flat sky patches, projecting the HEALPix666Hierarchical Equal
Area isoLatitude Pixelation of a sphere, http://healpix.sf.net. simulations
described above on $64\times 64$ pixel (that is, a $14.658\times 14.658$
square degrees area) planar images.

We have run the estimator on 2000 flat patches as described above. In order to
study the effect of the level Galactic contamination, we divide the sky in two
areas: 10000 simulations within a ‘Galactic’ band with Galactic latitude
$|b|\leq 10^{\circ}$ and 10000 within an external region with
$|b|>10^{\circ}$. The center sky coordinate of each patch is chosen randomly,
according to these latitude intervals and inside the simulated Wide Survey
observed area (see Figures 1 and 2). For each patch, we inject at the center a
point source with the FWHM listed in Table 1, a given flux density $S_{0}$ and
polarization fraction $\Pi$ randomly drawn from the log-normal distribution
(2) with the mean and standard deviation values $\langle\Pi\rangle$ and
$\sigma_{\Pi}$ described in Table 1. Please note that the PSM simulations
already contain resolved and unresolved polarized point sources apart from the
synthetic test sources we are injecting at the central position of each
simulated patch. Figure 3 shows the $Q$ and $U$ Stokes parameters for one of
our simulations. We simulate intensities $S_{0}$ in ten logarithmically spaced
values between 0.1 and 100 Jy (that is, 200 sources –100 of them in the band,
100 of them outside it– with $S_{0}=0.1$ Jy, 200 with $S_{0}=0.2154$ Jy, and
so on). In this way we get a sample of sources from moderately faint to extra
bright (and, since the polarization fraction follows distribution (2) with
average $\langle\Pi\rangle=0.02$, values of $P$ from below 1 mJy to a few tens
of Janskys).

Figure 1: $Q$-Stokes simulated QUIJOTE Wide Survey sky at 11 GHz. The false
colour bar is expressed in Jy. Figure 2: $Q$-Stokes simulated QUIJOTE Wide
Survey instrumental noise at 11 GHz. The false colour bar is expressed in Jy.
The non uniformity of the noise reflects the non uniform sky scanning strategy
of the telescope.

(a) Stokes Q (b) Stokes U

Figure 3: Projected sky patches for one of our simulations. Map units are
expressed in Jy.

##  4 Results

###  4.1 White noise

In all the cases analyzed with the white noise simulations, the performance of
the Bayesian FF at estimating the source polarization is better than that of
the Bayesian Rice method. This was also the case for their non-Bayesian
counterparts (Argüeso et al. 2009, etc.). Due to this and taking also into
account that the generalization of the Bayesian Rice method to full-sky
simulations is far from trivial777For spatially correlated noise, such as the
polarization produced by Galactic foregrounds, the distribution of $P$ is not
Ricean any more., from now on, we will only compare the Bayesian FF and the FF
techniques.

The left panel of Figure 4 shows the estimation of the polarized flux density
for the 10000 white noise simulations. The results have been binned into
eleven logarithmically spaced intervals in input $P_{0}$. We show in blue the
results from the Bayesian estimator (13). The error bars show the $68.27\%$
intervals of the corresponding empirical distributions. The filled circles
indicate the median value of the distribution of results; the dots indicate
the average value of the distribution. For comparison, the results of a MLE,
which correspond to the third and fourth terms of (13), are shown in orange.
The red dotted line shows the $P=P_{0}$ line. The MLE is equivalent to the FF
technique (Argüeso et al. 2009; López-Caniego et al. 2009; Herranz et al.
2012). Both the Bayesian FF and the MLE estimator work well for highly
polarized sources ($P_{0}\gtrsim 0.5$ Jy). For input polarization levels below
$\sim 0.2$ Jy (which is approximately the rms of the filtered noise of the
simulations), however, the MLE reaches a plateau: it is naturally limited by
the level filtered noise. The Bayesian estimator, on the contrary, uses the a
priori information on the $\Pi$ distribution and the knowledge of the source
flux density $S$ to predict lower $P$ values. As a matter of fact, the
Bayesian FF tends to overcompensate and predict, for the lower end of values
of the input $P_{0}$, polarized fluxes $P_{est}\sim 0.2P_{0}$. To see why this
happens, we will carry out a short theoretical calculation based on (12) . We
write the first part of (12) as a function of $P_{0}$.

$\displaystyle-$ $\displaystyle\log{f(Q_{0},U_{0}/Q_{i},U_{i})}=(\log
P_{0}-\log(\Pi_{med}S_{0}))^{2}/2\sigma_{\Pi}^{2}+$ $\displaystyle
2\,\log{(P_{0})}+\Sigma_{i}\displaystyle\frac{(Q_{i}-Q_{0}\tau_{i})^{2}}{2\sigma_{i}^{2}}+\Sigma_{i}\displaystyle\frac{(U_{i}-U_{0}\tau_{i})^{2}}{2\sigma_{i}^{2}}+K$

If we define

$\Sigma_{i}\displaystyle\frac{Q_{i}\tau_{i}}{\sigma_{i}^{2}}=P_{1}\cos{\phi_{1}}\,,Q_{0}=P_{0}\cos{\phi}$
(16)

and

$\Sigma_{i}\displaystyle\frac{U_{i}\tau_{i}}{\sigma_{i}^{2}}=P_{1}\sin{\phi_{1}}\,,U_{0}=P_{0}\sin{\phi}$
(17)

the previous formula can be expressed in terms of $P_{0}$, $P_{1}$ and the
polarization angles

$\displaystyle-$ $\displaystyle\log{f(Q_{0},U_{0}/Q_{i},U_{i})}=(\log
P_{0}-\log(\Pi_{med}S_{0}))^{2}/2\sigma_{\Pi}^{2}+$ $\displaystyle
2\,\log{(P_{0})}-P_{0}P_{1}\cos{(\phi-\phi_{1})}+P_{0}^{2}\,\Sigma_{i}\displaystyle\frac{\tau_{i}^{2}}{2\sigma_{i}^{2}}+K_{2}$

with

$\Sigma_{i}\,\displaystyle\frac{Q_{i}^{2}+U_{i}^{2}}{2\sigma_{i}^{2}}+K=K_{2}$
(19)

Taking the partial derivatives of (18) with respect to $P_{0}$ and $\phi$ and
equating them to zero, we obtain the estimators $\hat{P_{0}}$ and $\hat{\phi}$
that minimize minus the log-posterior. It can be easily seen that
$\hat{\phi}=\phi_{1}$ and $\hat{P_{0}}$ is the solution of the following
equation

$(\log\hat{P_{0}}-\log(\Pi_{med}S_{0}))/\sigma_{\Pi}^{2}+2-\hat{P_{0}}\,P_{1}+\hat{P_{0}}^{2}\,\Sigma_{i}\,\displaystyle\frac{\tau_{i}^{2}}{\sigma_{i}^{2}}=0$
(20)

This equation is very interesting: if we assume that $\hat{P_{0}}\,P_{1}\ll 1$
and
$\hat{P_{0}}^{2}\,\Sigma_{i}\,\displaystyle\frac{\tau_{i}^{2}}{\sigma_{i}^{2}}\ll
1$, that is, the source polarization is much lower than the noise, the
estimation will be dominated by the Bayesian prior and neglecting the last two
terms in (20), we find a constant value for the estimator, independently of
the data,

$\hat{P_{0}}=\Pi_{med}\,S_{0}\,e^{-2\sigma_{\pi}^{2}}$ (21)

For a lognormal distribution

$\Pi_{med}=\langle\Pi\rangle\exp\left(-\sigma_{\Pi}^{2}/2\right),$ (22)

Taking into account the values of $\langle\Pi\rangle$ and $\sigma_{\Pi}$ given
in Table 1, $\Pi_{med}=0.012$ and we obtain

$\hat{P_{0}}=0.00164\,S_{0}.$ (23)

In order to check these theoretical results, we have carried out $10000$
simulations with $S_{0}=1$ Jy. In this case, the noise, $\sigma=0.386$ Jy, is
much higher than the source polarization. We find, for all our simulations,
the estimated value $\langle\hat{P}_{0}\rangle=0.00166\pm 0.00002$, compatible
with the calculation above.

Though the estimator is constant in this case, this value is closer to the
real value than that obtained by using the matched filter, which is completely
dominated by the noise.

For higher values of $S_{0}$, e. g $S_{0}=10$ Jy, there are around $5000$
simulations, corresponding to the lower polarizations, that produce an
estimator close to the default value $0.016-0.020$ Jy. For higher values of
the real polarization, there is a combination of the prior and the matched
filter terms in the solution of (20). At any rate, the performance of the
Bayesian FF is better than that of the plain FF.

The right panel of Figure 4 shows the polarized flux estimation
error888Defined as $P_{0}-\hat{P}_{0}$, where $P_{0}$ is the input value and
$\hat{P}_{0}$ the estimated value of the polarization of the source, either
through the Bayesian method or through the MLE. as a function of the input
flux density (temperature) of the source $S_{0}$. For low flux densities, the
figure shows both the systematic overestimation of the MLE, due to the noise
limit, and the underestimation of the Bayesian FF estimator, due to the
reasons discussed above. In absolute terms the statistical error of the MLE is
much larger than that of the Bayesian FF estimator in the low flux density
regime. There is an interesting interval at intermediate flux densities ($\sim
10$ Jy) at which the errors of the MLE and Bayesian FF estimator are of the
same order, but in opposite directions. The Bayesian estimator seems to reach
a plateau (i.e., is noise-limited) around $P_{0}\sim 10$ mJy, a order of
magnitude below in polarized flux than the MLE.

(a)  (b)

Figure 4: Left: Binned estimates of the polarized flux density as a function
of binned input $P_{0}$, in Janskys, for the set of 10000 simulated white
noise patches. Right: Error in the estimation $P_{0}-\hat{P}_{0}$, as a
function of binned input $S_{0}$ (Jy), for the same set of simulations.
Bayesian estimations appear in blue, whereas maximum likelihood estimations
are shown with orange colour. Error bars show the $68.27\%$ interval of the
distribution of results in each case. Filled circles indicate the median of
the distribution; the diamonds (for the MLE) and crosses (for the Bayesian
estimator) indicate the average value of the distribution.

Figure 5 shows the absolute polarization angle error,

$|\Delta\phi|=|\phi_{0}-\hat{\phi}_{0}|,$ (24)

where $\phi_{0}$ and $\hat{\phi}_{0}$ are the input and estimated polarization
angles (in degrees), as a function of the input polarization of the source
$P_{0}$. The figure shows that there is little difference between the Galactic
and extragalactic areas, and between the Bayesian estimator and the MLE. This
is not a surprise, since the priors in (13) are constant with respect to
$\phi_{0}$, that is, the Bayesian estimator and the MLE should perform
similarly, as it is the case.

Figure 5: Binned estimates of $|\Delta\phi|=|\phi_{0}-\hat{\phi}_{0}|$ as a
function of binned input $P_{0}$, in Janskys, for the set of 10000 simulated
white noise images. Bayesian estimations are marked with blue crosses and
error bars, whereas maximum likelihood estimations are shown with orange
diamonds and error bars. The median values are indicated by means of large
blue filled circles (Bayesian estimator) and small organge squares (MLE).

###  4.2 Full sky simulations

Figure 6 shows the estimation of the polarized flux density for a) the 10000
‘Galactic’ ($|b|\leq 10^{\circ}$) and b) the 10000 ‘extragalactic’
($|b|>10^{\circ}$) simulated QUIJOTE sky patches. The results have been binned
into eleven logarithmically spaced intervals in input $P_{0}$. We show in blue
the results from the Bayesian estimator (13). The error bars show the
$68.27\%$ intervals of the corresponding empirical distributions. For
comparison, the results of a MLE, which correspond to the third and fourths
terms of (13), are shown in orange. The MLE is equivalent to the FF technique
(Argüeso et al. 2009; López-Caniego et al. 2009; Herranz et al. 2012). As it
happened in the case of the white noise simulations (Section 4.1), for low
flux density sources the MLE estimator reaches a plateau dominated by the
noise level (higher for Galactic than for extragalactic sources). The Bayesian
estimator reaches a similar plateau at much lower polarized fluxes, again in
the $\sim 10$ mJy regime instead of the $\sim 100$ mJy regime of the MLE.
Please note, however, that the distribution of the estimated $\hat{P}_{0}$ by
means of the Bayesian estimator becomes more and more skewed as $P_{0}$
decreases999This can be quickly seen by the growing differences between the
median and the average values of the distribution, as shown in the Figure..

(a)  (b)

Figure 6: Binned estimates of the polarized flux density as a function of
binned input $P_{0}$, in Janskys, for the set of 10000 simulated QUIJOTE sky
patches with $|b|\leq 10^{\circ}$ (left) and the set of 10000 simulated
QUIJOTE sky patches with $|b|>10^{\circ}$ (right) at 11 GHz. Bayesian
estimations appear in blue, whereas maximum likelihood estimations are shown
with orange colour. Error bars show the $68.27\%$ interval of the distribution
of results in each case. Filled circles indicate the median of the
distribution; the diamonds (for the MLE) and crosses (for the Bayesian
estimator) indicate the average value of the distribution.

Figure 7 shows the error of the estimation of $P$ as a function of the input
flux density of the sources (in Janskys). As in the case of white noise, the
MLE estimator tends to overestimate the polarized flux of faint sources
whereas the Bayesian FF tends to underestimate it. This error is a systematic
bias that tends to a constant value in relative terms, but decreases to zero
Janskys in absolute terms for $S_{0}\rightarrow 0$. Error bars are smaller for
extragalactic sources than for Galactic sources, which are embedded in more
intense foreground emission. Figures 8 and 9 show the normalized histograms of
the difference $\Delta P$ between the input polarization $P_{0}$ and the
estimated polarization $P$,

$\Delta P=P_{0}-\hat{P}_{0},$ (25)

for eleven different values of the total (Stokes I) flux density $S_{0}$.
‘Galactic’ sources are shown in Figure 8 and ‘extragalactic’ sources are shown
in Figure 9. The estimation $\hat{P}_{0}$ has been obtained with the Bayesian
Filtered Fusion method introduced in this paper (in blue) and the MLE (in red
color). For bright sources ($S_{0}>10$ Jy) the histograms are approximately
symmetric and centered around $\Delta P=0$, but for fainter sources the MLE
histograms are skewed to the left, showing the same kind of overestimation
already observed in Figure 6. The histograms for the Bayesian estimator,
however, are skewed to the right but much narrower than the MLE histograms,
which indicates that the Bayesian estimator predicts the polarization of a
source with a smaller margin of error. Both types of error, MLE-overestimation
and Bayesian FF-underestimation, must be dealt with in CMB polarization
experiments, but the amount of bias is signficantly smaller for the Bayesian
FF estimator.

Finally, Figure 10 shows the absolute polarization angle error, as a funcion
of the input polarization $P_{0}$. The figure shows that there is little
difference between the Galactic and extragalactic areas, and between the
Bayesian estimator and the MLE. This is not a surprise, since the priors in
(13) are constant with respect to $\phi_{0}$, that is, the Bayesian estimator
and the MLE should perform similarly, as it is the case.

(a)  (b)

Figure 7: Error in the estimation $P_{0}-\hat{P}_{0}$ (Jy), as a function of
binned input $S_{0}$ (Jy), for the set of 10000 simulated QUIJOTE sky patches
with $|b|\leq 10^{\circ}$ (left) and the set of 10000 simulated QUIJOTE sky
patches with $|b|>10^{\circ}$ (right) at 11 GHz. Bayesian estimations appear
in blue, whereas maximum likelihood estimations are shown with orange colour.
Error bars show the $68.27\%$ interval of the distribution of results in each
case. Filled circles indicate the median of the distribution; the diamonds
(for the MLE) and crosses (for the Bayesian estimator) indicate the average
value of the distribution. Figure 8: Normalized histogram of the difference
$\Delta P$ between the input polarization $P_{0}$ and the estimated
polarization for the Bayesian estimator (blue, /) and the MLE (red, \\), for
sources located within the Galactic band $|b|\leq 10^{\circ}$, and for nine
different values of the input total flux density $S_{0}$. Figure 9: Normalized
histogram of the difference $\Delta P$ between the input polarization $P_{0}$
and the estimated polarization for the Bayesian estimator (blue, /) and the
MLE (red, \\), for sources located outside the Galactic band $|b|>10^{\circ}$,
and for nine different values of the input total flux density $S_{0}$.

(a)  (b)

Figure 10: Binned estimates of $|\Delta\phi|=|\phi_{0}-\hat{\phi}_{0}|$ as a
function of binned input $P$, in Janskys, for the set of 500 simulated QUIJOTE
sky patches with $|b|\leq 10^{\circ}$ (left) and the set of 500 simulated
QUIJOTE sky patches with $|b|>10^{\circ}$ (right) at 11 GHz. Bayesian
estimations are marked with blue crosses, whereas maximum likelihood
estimations are marked with orange diamonds. Median values are indicated with
large blue filled circles and small orange filled squares, respectively.

###  4.3 A note on the robustness of the Bayesian estimator

Every time some prior information is used in the Bayesian framework the
inevitable question arises: what is the effect of a wrong guess of the prior
in the estimation? In order to shed some light on this we have re-analyzed the
one hundred simulations of ‘extragalactic’ sources with flux density
$S_{0}=10$ Jy101010We have chosen this particular flux density value because
according to Figure 9 it marks the flux density for which the Bayesian
estimator begins to outperform the MLE.. Instead of using the correct value of
the median polarization fraction $\Pi_{med}$ in equation (13) we use a biased
parameter $\Pi_{med}^{b}=b\,\Pi_{med}$ with $b=[0.5,0.6,\ldots,2.0]$, that is,
we have tested what happens if our guess of the median polarization fraction
is wrong by a factor from $50\%$ to $200\%$. Figure 11 shows the average
estimation error $P_{0}-\hat{P}_{0}$ as a function of the bias factor $b$.
Error bars show the $68.27\%$ intervals of the corresponding empirical
distributions. The figure shows that the average error of the Bayesian
estimator varies smoothly with the bias in the prior. For comparison, for the
same simulations the MLE produces a (bias independent, as the maximum
likelihood estimator does not use prior information) value for the error
$P_{0}-\hat{P}_{0}=-0.22\pm 0.12$, larger than the Bayesian estimation (for
this particular value of $S_{0}$) even when the prior is wrong by a factor of
two.

Figure 11: Error in the estimation of the total polarization of a test source
with $S_{0}=10$ Jy as a function of the bias factor $b$ (defined as
$\Pi_{med}^{b}=b\,\Pi_{med}$) affecting the Bayesian prior on $\Pi_{med}$.

Another potential source of bias is the uncertainty on the true flux density
of the source. The estimators in equations (12) and (13) depend implicitly on
an a priori knowledge of the source flux density $S_{0}$ through the factor
$\mu_{1}=\mu+\log(S_{0})$. In the previous tests we have assumed that $S_{0}$
is known with arbitrary precision, but in practice this will not be the case.
In a real experiment one expects to know some reasonable estimation
$\hat{S}_{0}$ of the true flux density of the source. In a typical CMB
experiment setting the difference between $S_{0}$ and $\hat{S}_{0}$ will be
relatively small (at least in comparison with the relative difference between
$P_{0}$ and $\hat{P}_{0}$), but not zero. The uncertainty on the source flux
density can bias the estimators (12) and (13) even if the distribution of
$\hat{S}_{0}$ is symmetric around $S_{0}$, as $S_{0}$ enters the estimators in
a non-linear fashion. Moreover, one expects the uncertainty in $S$ to increase
the statistical error of the estimators.

In order to test the effect of the uncertainty on $S$ on our Bayesian
estimator, we have conducted a new batch of 10000 simulations in the same
fashion as described in Section 4.2. The analysis follows the same pipeline as
described above, but every time we compute the estimator (13) we introduce a
random photometric error in $S_{0}$. These photometric errors follow a
Gaussian distribution of standard deviation $\sigma=0.3$ Jy, a little smaller
than the QUIJOTE simulation noise rms level111111We assume that the rms of the
photometric errors has been lowered by means of some filtering scheme, such as
a matched filter or a Mexican Hat Wavelet, or any other suitable signal
processing technique. Then the $\sim 0.3$ Jy uncertainty becomes a more
realistic approximation of error in the determination the flux density of
compact sources in the QUIJOTE Wide Survey..

(a)  (b)

Figure 12: Error in the estimation $P-P_{0}$ (Jy), as a function of binned
input $S_{0}$ (Jy), for the set of 10000 simulated QUIJOTE sky patches with
$|b|\leq 10^{\circ}$ (left) and the set of 10000 simulated QUIJOTE sky patches
with $|b|>10^{\circ}$ (right) at 11 GHz. Bayesian estimations for a perfect
photometry of the source total flux density $S_{0}$ appear in blue, whereas
Bayesian estimations including a $0.3$ Jy uncertainty in $S_{0}$ are shown
with orange. Error bars show the $68.27\%$ interval of the distribution of
results in each case. Filled circles indicate the median of the distribution;
dots indicate the average value of the distribution. Orange points and lines
have been slightly displaced to the right in order to make the figure more
readable.

Figure 12 shows the average error of the estimation of the polarization of our
simulated sources comparing the two cases: if the source flux density $S_{0}$
is perfectly known in advance (blue dots and error bars) or if a $0.3$ Jy
photometric uncertainty is present in the analysis (orange dots and bars,
slightly displaced to the right for the sake of clarity). Galactic and
extragalactic cases (as defined above) are shown in the left and right panels,
respectively. The effect of a $\sim 0.3$ Jy photometric error on the flux
density of the sources is negligible in our simulated experimental setting.
This comes not as a suprise, since a $\sim 0.3$ Jy variation in $S_{0}$
produces only a $\sim 10\%$ change in the $\mu_{1}$ term that appears in (12)
and (13) in the worst case (1 Jy sources)121212Moreover, the estimation of the
polarization is not given directly by (12) and (13), but by the minimization
of these functions. A small variation in one of the terms of the functions
does not necessarily mean that the _position_ of the minimum of the function
changes in a noticeable way. The non linear way in which $S_{0}$ appears in
these equations makes it difficult to find an analytical expression of how an
uncertainty in $S_{0}$ affects the minimization. This question is better
answered by simulations, just as we have done in this section.. This
discrepancy quickly decreases as $S_{0}$ grows. Moreover, the rms around the
mean $\langle\mu_{1}\rangle$ also decreases very quickly with $S_{0}$.
Therefore, we conclude that our Bayesian estimator is robust against moderate
uncertainties on the prior and the flux density of the sources.

##  5 Conclusions

The estimation of the polarimetric properties of extragalactic compact sources
at microwave wavelengths will be very relevant in the upcoming years. In this
work, we have introduced a Bayesian approach for the estimation of the
polarized flux density $P$ of this kind of sources. Following recent works by
Massardi et al. (2013); Galluzzi et al. (2017, 2019) among others, we have
proposed an analytical prior for the polarization fraction of extragalactic
radio sources which takes the form of a log-normal distribution whose
parameters (median, average and variance values of the polarization fraction)
can be constrained by the latest observational data. Using this prior, we have
proposed two maximum _a posteriori_ (MAP) estimators of the polarization of a
given source given observations of its $Q$ and $U$ Stokes parameters. The
first method works directly on the quadratic combination $P^{2}=Q^{2}+U^{2}$
whereas the second method produces individual estimators of the ground-truth
values $Q_{0}$ and $U_{0}$ that are then quadratically added to give an
estimator of the ground-truth polarization $P_{0}$ of the source. We have
called these methods Bayesian Rice and Bayesian Filtered Fusion (BFF),
respectively. Both can be considered as natural Bayesian extensions of the
frequentist Neyman-Pearson and standard Filtered Fusion (FF) methods
introduced by Argüeso et al. (2009). The standard FF is shown to be equal to
the Maximum Likelihood Estimator (MLE) for $P$, whereas the BFF adds to the
MLE a number of additional terms that include the _a priori_ information on
the distribution of the polarization fraction. The BFF method can be easily
accommodated to non-white noise and foregrounds. For this reason we have
focused on this method in most of our paper.

We have tested the performance of the BFF method and compared it to that of FF
using two sets of simulations: polarized sources embedded in $Q$ and $U$ white
noise, and more realistic simulations that include also polarized CMB and
Galactic foreground emission. In both cases we have used the pixel and beam
scales plus the noise levels and sky coverage of the QUIJOTE experiment Wide
Survey (Rubiño-Martín et al. 2020; Herranz et al. 2020) at 11 GHz. For the BFF
we assumed that the flux density $S_{0}$ of the sources is perfectly known.
For highly polarized sources the two methods yield the same results, but for
medium to low polarizations ($P_{0}\leq 400$ mJy in our simulations) the BFF
gets more accurate estimations of the polarization of the sources. The FF gets
noise-limited around a polarization flux $P_{0}\sim 500$ mJy, whereas the BFF
allows us to reach polarized fluxes well below $P_{0}\sim 100$ mJy before
becoming noise-limited itself. Both estimators are biased for low polarization
(i.e. $P_{0}\lesssim 500$ mJy) sources: the BFF tends to underestimate the
polarization, whereas the standard FF overestimates the polarization of these
sources. In the case of the FF the bias is due to noise boosting of the signal
(akin to Eddington bias). In the case of the BFF, the bias is originated by
the extra terms in the estimator formula that come from the physical prior.
However, the absolute value of bias is significantly smaller for the BFF than
for the FF, specially for faint sources.

In the above discussion, we have assumed that the prior describes the real
distribution of polarization of the sources adequately and that the total flux
density $S_{0}$ of each source is perfectly known. However, information about
the polarization properties of extragalactic sources at microwave frequencies
above $\simeq 10$ GHz is still scarce. Moreover, for any given source $S_{0}$
is known with a certain degree of uncertainty (due to instrumental noise,
less-than-perfect modelling of the spectral energy distribution of the source
and variability, among other possible causes). In the last part of this work,
we have tested the robustness of the BFF estimator against moderate changes in
the prior parameters and realistic uncertainties in the flux density of the
sources. Our simulations indicate that assuming the wrong prior has a mild
effect on the Bayesian estimator. For example, for a $S_{0}=10$ Jy source, a
change by a factor of two in the assumed median polarization fraction of the
sources introduces errors or the order $\lesssim 100$ mJy in the estimation of
$P$. Regarding uncertainties in the flux density of the sources, we find that
non-catastrophic photometric error bars have a minimal impact on the
estimation of $P$.

We therefore conclude that the Bayesian approach can significantly improve the
estimation of the polarization of extragalactic radio sources in current and
upcoming CMB polarization experiments. In an upcoming work, we will explore
the extension of the Bayesian framework to the multi-frequency case.

###### Acknowledgements.

We thank the Spanish MINECO and the Spanish Ministerio de Ciencia, Innovación
y Universidades for for partial financial support under projects
AYA2015-64508-P and PGC2018-101814-B-I00, respectively. DH also acknowledges
funding from the European Union’s Horizon 2020 research and innovation
programme (COMPET-05-2015) under grant agreement number 687312
(RADIOFOREGROUNDS). Some of the results in this paper have been derived using
the HEALPix (Górski et al. 2005) and healpy (Zonca et al. 2019) packages. This
research made use of astropy,131313http://www.astropy.org a community-
developed core Python package for Astronomy (Astropy Collaboration et al.
2013, 2018), matplotlib, a Python library for publication quality graphics
(Hunter 2007), and SciPy, a Python-based ecosystem of open-source software for
mathematics, science, and engineering (Virtanen et al. 2020). We acknowledge
Santander Supercomputacion support group at the University of Cantabria (UC)
who provided access to the supercomputer Altamira Supercomputer at the
Institute of Physics of Cantabria (IFCA-UC-CSIC), member of the Spanish
Supercomputing Network141414https://www.res.es/en/about, for performing
simulations/analyses.

## References

  * Argüeso et al. (2009) Argüeso, F., Sanz, J. L., Herranz, D., López-Caniego, M., & González-Nuevo, J. 2009, MNRAS, 395, 649
  * Argüeso et al. (2011) Argüeso, F., Luis Sanz, J., & Herranz, D. 2011, Signal Processing, 91, 1527
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Bonavera et al. (2017a) Bonavera, L., González-Nuevo, J., Argüeso, F., & Toffolatti, L. 2017a, MNRAS, 469, 2401
  * Bonavera et al. (2017b) Bonavera, L., González-Nuevo, J., De Marco, B., Argüeso, F., & Toffolatti, L. 2017b, MNRAS, 472, 628
  * Crow & Shimizu (1988) Crow, E. & Shimizu, K. 1988, Lognormal distributions : theory and applications (New York: M. Dekker)
  * Delabrouille et al. (2013) Delabrouille, J., Betoule, M., Melin, J.-B., et al. 2013, A&A, 553, A96
  * Delabrouille et al. (2018) Delabrouille, J., de Bernardis, P., Bouchet, F. R., et al. 2018, Journal of Cosmology and Astroparticle Physics, 2018, 014
  * Diego-Palazuelos et al. (2020) Diego-Palazuelos, P., Vielva, P., & Herranz, D. 2020, submitted to Journal of Cosmology and Astroparticle Physics
  * Galluzzi et al. (2018) Galluzzi, V., Massardi, M., Bonaldi, A., et al. 2018, MNRAS, 475, 1306
  * Galluzzi et al. (2017) Galluzzi, V., Massardi, M., Bonaldi, A., et al. 2017, MNRAS, 465, 4085
  * Galluzzi et al. (2019) Galluzzi, V., Puglisi, G., Burkutean, S., et al. 2019, MNRAS, 489, 470
  * Génova-Santos et al. (2015) Génova-Santos, R., Rubiño-Martín, J. A., Rebolo, R., et al. 2015, in Highlights of Spanish Astrophysics VIII, ed. A. J. Cenarro, F. Figueras, C. Hernández-Monteagudo, J. Trujillo Bueno, & L. Valdivielso, 207–212
  * Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759
  * Herranz et al. (2012) Herranz, D., Argüeso, F., & Carvalho, P. 2012, Advances in Astronomy, 2012, 410965
  * Herranz et al. (2020) Herranz, D., López-Caniego, M., Génova-Santos, R., et al. 2020, A&A, in preparation
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90
  * Jackson et al. (2010) Jackson, N., Browne, I. W. A., Battye, R. A., Gabuzda, D., & Taylor, A. C. 2010, MNRAS, 401, 1388
  * López-Caniego (2016) López-Caniego, M. 2016, IAU Focus Meeting, 29, 54
  * López-Caniego et al. (2009) López-Caniego, M., Massardi, M., González-Nuevo, J., et al. 2009, ApJ, 705, 868
  * Massardi et al. (2011) Massardi, M., Bonaldi, A., Bonavera, L., et al. 2011, MNRAS, 415, 1597
  * Massardi et al. (2013) Massardi, M., Burke-Spolaor, S. G., Murphy, T., et al. 2013, MNRAS, 436, 2915
  * Massardi et al. (2008) Massardi, M., Ekers, R. D., Murphy, T., et al. 2008, MNRAS, 384, 775
  * Murphy et al. (2010) Murphy, T., Sadler, E. M., Ekers, R. D., et al. 2010, MNRAS, 402, 2403
  * Planck Collaboration I (2016) Planck Collaboration I. 2016, A&A, 594, A1
  * Planck Collaboration XXVI (2016) Planck Collaboration XXVI. 2016, A&A, 594, A26
  * Puglisi et al. (2018) Puglisi, G., Galluzzi, V., Bonavera, L., et al. 2018, ApJ, 858, 85
  * Rayner et al. (2000) Rayner, D. P., Norris, R. P., & Sault, R. J. 2000, Monthly Notices of the Royal Astronomical Society, 319, 484
  * Remazeilles et al. (2018) Remazeilles, M., Banday, A. J., Baccigalupi, C., et al. 2018, Journal of Cosmology and Astroparticle Physics, 2018, 023
  * Rice (1945) Rice, S. O. 1945, Bell System Technical Journal, 24, 46
  * Rubiño-Martín et al. (2012) Rubiño-Martín, J. A., Rebolo, R., Aguiar, M., et al. 2012, in Proc. SPIE, Vol. 8444, Ground-based and Airborne Telescopes IV, 84442Y
  * Rubiño-Martín et al. (2010) Rubiño-Martín, J. A., Rebolo, R., Tucci, M., et al. 2010, Astrophysics and Space Science Proceedings, 14, 127
  * Rubiño-Martín et al. (2020) Rubiño-Martín, J. A. et al. 2020, A&A, in preparation
  * Sadler et al. (2006) Sadler, E. M., Ricci, R., Ekers, R. D., et al. 2006, MNRAS, 371, 898
  * Sailer et al. (2020) Sailer, N., Schaan, E., & Ferraro, S. 2020, arXiv e-prints, arXiv:2007.04325
  * Sajina et al. (2011) Sajina, A., Partridge, B., Evans, T., et al. 2011, ApJ, 732, 45
  * Sekimoto et al. (2018) Sekimoto, Y., Ade, P., Arnold, K., et al. 2018, in Space Telescopes and Instrumentation 2018: Optical, Infrared, and Millimeter Wave, ed. M. Lystrup, H. A. MacEwen, G. G. Fazio, N. Batalha, N. Siegler, & E. C. Tong, Vol. 10698, International Society for Optics and Photonics (SPIE), 613 – 629
  * Simmons & Stewart (1985) Simmons, J. F. L. & Stewart, B. G. 1985, A&A, 142, 100
  * Trombetti et al. (2017) Trombetti, T., Burigana, C., De Zotti, G., Galluzzi, V., & Massardi, M. 2017, ArXiv e-prints, arXiv:1712.08412
  * Trombetti et al. (2018) Trombetti, T., Burigana, C., De Zotti, G., Galluzzi, V., & Massardi, M. 2018, A&A, 618, A29
  * Tucci et al. (2005) Tucci, M., Martínez-González, E., Vielva, P., & Delabrouille, J. 2005, MNRAS, 360, 935
  * Tucci & Toffolatti (2012) Tucci, M. & Toffolatti, L. 2012, Advances in Astronomy, 2012, 624987
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Zonca et al. (2019) Zonca, A., Singer, L., Lenz, D., et al. 2019, Journal of Open Source Software, 4, 1298