## Abstract

In this work, we deal with the problem of simultaneous multifrequency detection of extragalactic point sources in maps of the Cosmic Microwave Background. We developed a linear filtering technique that takes into account the spatial and the cross-power spectrum information at the same time. By using realistic simulations, we compared this technique with the monofrequency matched filter in terms of completeness, reliability, flux and spectral index accuracy. The multifrequency method outperforms the matched filter in all the studied cases in this work.

## 1 INTRODUCTION

Over the last few years, a big effort has been devoted to the problem of detecting point sources in Cosmic Microwave Background (CMB) experiments. The main reason for this effort is that modern CMB experiments have reached resolution and sensitivity levels such that their capability to estimate the statistics of CMB fluctuations at high multipoles is no longer limited by instrumental noise but by Galactic and extragalactic foreground contamination. Among the extragalactic contaminants, point sources are the most relevant, both in temperature (Toffolatti et al. 1998; De Zotti et al. 1999; Hobson et al. 1999; De Zotti et al. 2005) and in polarization (Tucci et al. 2004, 2005; López-Caniego et al. 2009). Moreover, they are one of the most difficult contaminants to deal with and, at least in the frequency range spanned by CMB experiments, one of the most poorly known. Since the point sources contaminate the CMB radiation, hamper the efforts of cleaning CMB from Galactic foregrounds and alter dramatically the statistics of the CMB fluctuations (e.g. the Gaussianity), it is therefore necessary to detect the maximum possible number of extragalactic point sources (EPS) and to estimate their flux with the lowest possible error before any serious attempt to study the CMB anisotropies is made.

However, EPS are not just a contaminant that should be eliminated. From the point of view of extragalactic astronomy, the study of EPS at microwave frequencies is a very interesting topic. The next generation of CMB experiments will allow one to obtain all-sky EPS catalogues that will fill in the existing observational gap in our knowledge of the Universe in the frequency range from 20 to roughly 1000 GHz. We expect to derive source number counts and spectral indices, to study source variability and to discover rare objects such as inverted spectrum radio sources, extreme gigahertz peaked spectrum sources (GPS) and high-redshift dusty proto-spheroids (see for instance the *Planck Bluebook*; The Planck Collaboration 2006). The first results of this new era of CMB experiments are already here: the *Wilkinson Microwave Anisotropy Probe* (*WMAP*) satellite (Bennett et al. 2003) has allowed us to obtain the first all-sky point source catalogues above ∼0.8–1 Jy in the 23–94 GHz frequency range in temperature (Bennett et al. 2003; Hinshaw et al. 2007; López-Caniego et al. 2007; Chen & Wright 2008; Massardi et al. 2009; Wright et al. 2009) and polarization (López-Caniego et al. 2009), being complete the last three ones. In a recent paper by González-Nuevo et al. (2008), it has been shown how these catalogues can be used to study the statistics of point sources in that range of frequencies. The *Planck* mission (Tauber 2005) will allow to extend these catalogues to lower flux limits and frequencies up to 857 GHz.

We have mentioned that the detection and estimation of the flux of EPS are a difficult task. The main reason for this is that the many different types of EPS distributed in the sky form a very heterogeneous set of objects that do not have a common spectral behaviour. While other foreground contaminants follow a specific emission law that is more or less well known (or can be inferred from observations) and that varies relatively slowly and continuously across the sky, each EPS has an emission law that, in principle, can be totally different to any other and independent from them. From the point of view of statistical signal processing, the problem of detecting EPS is a case of sorely underdetermined component separation problem where the number *M* of components is much larger than the number *N* of frequency channels.

Let us consider the case of another contaminant that is of relevance in CMB images and to which the EPS problem has some similarities: the Sunyaev–Zel'dovich (SZ) effect (Sunyaev & Zeldovich 1970, 1972; Carlstrom, Holder & Reese 2002). As in the case of EPS, the SZ contamination occurs in small and compact regions of the sky and it affects the high-ℓ region of the fluctuations angular power spectrum. However, in the case of the thermal SZ effect the frequency dependence is very well known (ignoring relativistic corrections). This has made possible the development of a number of powerful detection techniques that are specifically suited to the problem of SZ effect detection in CMB images (Diego et al. 2002; Herranz et al. 2002a; Hobson & McLachlan 2003; Pierpaoli et al. 2005; Herranz et al. 2005; Vale & White 2006; Pires et al. 2006; Melin, Bartlett & Delabrouille 2006; Bartlett et al. 2008; Leach et al. 2008) or to use more generic diffuse component separation techniques to obtain SZ maps and catalogues (see for example, among many others, Maino et al. 2002; Stolyarov et al. 2002; Delabrouille, Cardoso & Patanchon 2003; Eriksen et al. 2004; Bedini et al. 2005; Bonaldi et al. 2006).

Unfortunately, these approaches are not directly applicable to EPS detection. The most commonly used approach to this problem consists of working separately in each channel. The key idea is to take advantage of the fact that all the EPS have the same shape (basically, that of the beam). In the field of CMB images, wavelet techniques (Vielva et al. 2001, 2003; González-Nuevo et al. 2006; Sanz et al. 2006; López-Caniego et al. 2007), matched filters (MF; Tegmark & de Oliveira-Costa 1998; Barreiro et al. 2003; López-Caniego et al. 2006) and other related linear filtering techniques (Sanz, Herranz & Martínez-Gónzalez 2001; Chiang et al. 2002; Herranz et al. 2002c; López-Caniego et al. 2004, 2005a,b) have proved to be useful. All these techniques rely on the prior knowledge that the sources have a distinctive spatial behaviour and this fact is used to design some bandpass filter to enhance them with respect to the noise. Detection can be further improved by including prior information about the sources, i.e. some knowledge about their flux distribution, in the frame of a Bayesian formalism (Hobson & McLachlan 2003; Carvalho, Rocha & Hobson 2009).

With the previous single-frequency methods and for the case of the *Planck* mission, we expect to reach detection limit fluxes that range from a few hundred mJy to several Jy, depending on the frequency channel, and to obtain catalogues that will contain several hundred to a few thousand objects (López-Caniego et al. 2006; Leach et al. 2008). This is satisfactory, but we feel that we could do better if we were able to use multiwavelength information in some way.

For now, multiwavelength detection of EPS in CMB images remains a largely unexplored field. In recent years, some attempts have been done in this direction (Naselsky, Novikov & Silk 2002; Chen & Wright 2008; Wright et al. 2009). More recently, Herranz & Sanz (2008) have introduced the technique of ‘matched matrix filters’ (MTXF) as the first fully multifrequency, non-parametric, linear filtering technique that is able to find EPS and to do unbiased estimations of their fluxes thanks to the distinctive spatial behaviour of the sources. This method incorporates at the same time some multiwavelength information, without assuming any specific spectral behaviour for the sources. Herranz et al. (2009) have applied the MTXF to realistic simulations of the *Planck* radio channels, showing that it is possible to practically double the number of detections, for a fixed reliability level, for some of the channels with respect to the single-frequency MF approach.

MTXF use multiwavelength information in such a way that it is not necessary to make any assumption about the spectral behaviour of the sources. In fact, in that formalism their spectral behaviour is entirely irrelevant. All the multiwavelength considerations concern only to the noise1 and its correlations. In this sense, the MTXF method deals with only half of the problem. This has its advantages in terms of robustness and reliability, but one could wish to have a technique that uses multiwavelength information in the modelling of *both* the signal (EPS) and the noise. But, as we mentioned before, the spectral behaviour of the EPS is not known.

In this paper, we will show that even if the spectral behaviour of the EPS is unknown a priori, it is still possible to determine it directly from the data by means of an adaptive filtering scheme that incorporates multifrequency information not only through the noise correlations among channels, but also about the sources themselves.

The problem is, in more than one sense, similar to the problem of detecting SZ clusters (that is the reason we discussed that case a few paragraphs above). In the SZ case, the spectral behaviour of the sources is known, but not their size. A way to deal with this problem is to introduce the scale of the source as a free parameter (for example, the cluster core radius *r*_{c}) in the design of a ‘*matched multifilter*’ (MMF) and to optimize the value of this parameter for each source so that a maximum signal to noise ratio is obtained after filtering (see the details in Herranz et al. 2002a,b; Herranz et al. 2005; Melin et al. 2006; Schäfer et al. 2006). As the problem depends on the optimization of one single parameter, the method is easy to implement in codes that are relatively fast.

In this paper, we introduce a modification of the MMF technique in which the mixing coefficients of the frequency dependence vector of the sources are considered as free parameters to be optimized. As we will show, if the number of frequency channels is *N* and we choose wisely a fiducial frequency of reference, the number of free parameters to optimize is *N*− 1. For simplicity, throughout this paper, we will use as an example the case of two channels, *N*= 2, but the method is valid for any number of channels *N* > 1.

The structure of this paper is as follows. In Section 2, we will summarize the formulae of the MMF and we will introduce the modification that allows us to use it for detecting EPS in arbitrarily chosen pairs of frequency channels. We will show that the filter can be in this case factorized in such a way that a particularly fast implementation of the method can be achieved. In Section 3, we will describe the simulations that we use to test the method. The results of the exercise will be described in Section 4. Finally, in Section 5, we will draw some conclusions.

## 2 METHOD

### 2.1 The single-frequency approach

Let us assume a set of images corresponding to the same area of the sky observed simultaneously at *N* different frequencies,

*N*. At each frequency ν,

*y*

_{ν}is the total signal in the pixel

*and*

**x***s*

_{ν}represents the contribution of the point source to the total signal

*y*

_{ν}; for simplicity let us assume that there is only one point source centred at the origin of the image;

*f*

_{ν}is the frequency dependence of the point source; and

*n*

_{ν}is the background or generalized noise (containing not only the instrumental noise, but also the contributions of the rest of the components).

The intrinsic angular size of the point sources is smaller than the angular resolution of the detector. At each observing frequency, each source is convolved with the corresponding antenna beam. For simplicity, we will assume that the antenna beam can be well described by a symmetric 2D Gaussian function. Then, we can write

where*x*= |

**| (since we are considering symmetric beams),**

*x**A*is the amplitude of the source, and τ is the spatial template or profile. The background

*n*

_{ν}(

**) is modelled as a homogeneous and isotropic random field with average value equal to zero and power spectrum**

*x**P*

_{ν}defined by where

*n*

_{ν}(

**) is the Fourier transform of**

*q**n*

_{ν}(

**) and δ**

*x*^{2}

_{D}is the 2D Dirac distribution.

In the single-frequency approach, each channel is processed separately and independently from the other frequency channels. This approach is robust in the sense that it is not necessary to assume anything about the spectral behaviour of the sources. The main drawback of the single-frequency approach, however, is that one misses the potential noise reduction that could be obtained with a wise use of the information present at the other frequencies.

The standard single-frequency point-source detection method in the literature is based on the MF (Tegmark & de Oliveira-Costa 1998; Barreiro et al. 2003; López-Caniego et al. 2006). The MF is the optimal linear detector for a single map in the sense that it gives the maximum signal-to-noise amplification. The MF can be expressed in Fourier space in the following way:

Here,*a*is a normalization factor that preserves the source amplitude after filtering.

### 2.2 The matched multifilter

In the multifrequency approach, we take into account the statistical correlation of the noise between different frequency channels and the frequency dependence of the sources. Now let us model the background *n*_{ν}(** x**) as a homogeneous and isotropic random field with average value equal to zero, and cross-power spectrum defined by

*n*

_{ν}(

**) is the Fourier transform of**

*q**n*

_{ν}(

**) and δ**

*x*^{2}

_{D}is the 2D Dirac distribution. Let us define a set of

*N*linear filters ψ

_{ν}that are applied to the data here,

*defines a translation. The right-hand part of equation (6) shows the filtering in the Fourier space, where*

**b***y*

_{ν}(

**) and ψ**

*q*_{ν}(

*q*) are the Fourier transforms of

*y*

_{ν}(

**) and ψ**

*x*_{ν}(

**), respectively. The quantity**

*x**w*

_{ν}(

**) is the filtered map ν at the position**

*b**. The total filtered map is the sum Therefore, the total filtered field is the result of two steps: (a) filtering and (b) fusion. During the first step, each map*

**b***y*

_{ν}is filtered with a linear filter ψ

_{ν}, during the second step the resulting filtered maps

*w*

_{ν}are combined so that the signal

*s*is boosted while the noise tends to cancel out. Note that the combination in equation (7) is completely general, since any summation coefficients different from the other can be absorbed in the definition of the filters ψ

_{ν}. Then the problem consists of how to find the filters ψ

_{ν}so that the total filtered field is optimal for the detection of point sources.

The total filtered field *w* is optimal for the detection of the sources if:

*w*() is an unbiased estimator of the amplitude of the source, so 〈*0**w*()〉=*0**A*;the variance of

*w*() is minimum, that is, it is an efficient estimator of the amplitude of the source.*b*

If the profiles τ_{ν} and the frequency dependence *f*_{ν} are known, and if the cross-power spectrum is known or can be estimated from the data, the solution of the problem is already known: the MMF (Herranz et al. 2002a),

**Ψ**(

*q*) is the column vector

**Ψ**(

*q*) =[ψ

_{ν}(

*q*)],

**is the column vector**

*F***=[**

*F**f*

_{ν}τ

_{ν}] and

^{−1}is the inverse matrix of the cross-power spectrum . Finally, we can obtain the variance of the total filtered field, given by the following expression:

### 2.3 MMF with unknown source frequency dependence

As it was previously discussed, the problem is that the frequency dependence *f*_{ν} of the sources is not known a priori. Then, the possibilities are either (a) to admit defeat, returning to the single-frequency approach, (b) to devise a filtering method that does not use the frequency dependence of the sources altogether or (c) to model somehow the unknown frequency dependence in the framework of some optimization scheme. The second approach was explored in Herranz & Sanz (2008) and Herranz et al. (2009). In this work, we will study the third approach to the problem.

Before addressing this problem, it will be useful to rewrite equation (8) in a slightly different way. Let us write the vector ** F**=[

*f*

_{ν}τ

_{ν}] in matrix form as

*q*) = diag [τ

_{1}(

*q*), …, τ

_{N}(

*q*)] and

**=[**

*f**f*

_{ν}] the vector of frequency dependence. Note that all the dependence in

*q*is included in the matrix ; this fact will be useful later on.

Now imagine that * f* describes the true (unknown) frequency dependence of the sources and that

**=[**

*g**g*

_{ν}], ν= 1, …,

*N*is a new vector of equal size as

*, but whose elements can take any possible value. We can define the MMF for vector*

**f***,*

**g**^{t}= and that vector

*does not depend on*

**g***and can therefore go out of the integral. When applied to a set of images where there is present a point source with true amplitude*

**q***A*and true frequency dependence

*, the filters*

**f****Ψ**

_{g}will lead to an estimation of the amplitude, Note that if

**≠**

*g***, then**

*f**A*

_{g}≠

*A*. On the other hand, the variance of the filtered field would be, in analogy with equation (9), σ

^{2}

_{g}=α

_{g}. Let us finally define the signal-to-noise ratio (SNR) of the source in the total filtered map as Then, we can ask what the vector

*is that maximizes the signal-to-noise ratio*

**g***SNR*

_{g}. Intuition alone indicates that

*SNR*

_{g}is maximum if and only if

**=**

*g***. This can be formally proved with little effort by taking variations of**

*f**.*

**g**Then the problem can be solved via a maximization algorithm. In the case of a non-blind search, where the position of a given point source is known, one can focus on that point source and iteratively try values of the elements of * g* until a maximum signal to noise is reached. In the case of a blind search, the situation is a little bit more difficult because in a given image there may be many different objects

*s*with a different solution

_{i}

*g*_{i}. A way to proceed is to filter many times the image, using each time a different set of values of the elements of

*so that the appropriate range of frequency dependence is sufficiently well sampled, and then to proceed counting one by one all the possible detections and associating to each one the values of*

**g***that maximize the SNR of that source in particular.*

**g**We would like to remark that this situation is very similar to the case of the detection of galaxy clusters with unknown angular size described in Herranz et al. (2002a,b). In that case, the frequency dependence * f* was known but the size of the clusters (their source profile) was not. The cluster profile can be parametrized as a modified beta profile with a free scale parameter

*r*

_{c}(typically, the cluster core radius). In Herranz et al. (2002a,b), it was shown that the true scale of the clusters can be determined by maximizing the SNR of the detected clusters as a function scale

*r*

_{c}of the filter.

In our case, factorization (10) leads to equations (11) and (12); this is very convenient for implementation of the MMF when many filtering steps are necessary. The most time-consuming part of the filtering is the calculation of matrices and because they must be calculated for all values of * q*. In the case of clusters with unknown size, had to be calculated for every value of

*r*

_{c}. However, in the case we are considering in this paper the only quantities that vary during the maximization process are the elements of vector

*. This allows us to compute the integrals of matrix only once for each set of images. As a result, applying the MMF to large numbers of point sources with unknown frequency dependence is, in general, much faster than applying the MMF to the same number of clusters with unknown source profile.*

**g**The main difference is that while in the case of galaxy clusters it was necessary to maximize with respect to only one single parameter (the core radius), in the case of the unknown frequency dependence it is necessary to maximize with respect to the *N* components of vector * g*. This procedure can require a very large number of computations if

*N*is big. Although we have just seen that each free parameter of the frequency dependence can be mapped much faster than each free parameter of the source profile, we are still interested in reducing the number of computations as much as possible.

### 2.4 Number of degrees of freedom of vector *g*

For *N* images, vector ** g**=[

*g*

_{1}, …,

*g*] has

_{N}*N*degrees of freedom. This makes the optimization procedure more complex and computationally expensive. The situation can be lightened if we choose one of the frequencies under consideration as our fiducial frequency of reference. Let us choose for example a concrete frequency

*j*∈{1, …,

*N*} to be our fiducial frequency of reference, then

_{j}is normalized to unity,

*f*must be equal to one. Hence, we must look for vectors

_{j}**=**

*g**g*

_{1}, …,

*g*

_{j−1}, 1,

*g*

_{j+1}, …,

*N*and the number of independent degrees of freedom is

*N*− 1. If the number of channels is

*N*= 2, there is only one degree of freedom for the optimization problem. The case

*N*> 2 does not introduce any modification in the strategy. As mentioned in Section 2.3, the SNR of the sources has a global maximum if and only if

**=**

*g***. This implies that the problem is not degenerated. The only complication is computational: as the number of free parameters to determine increases, more operations are needed in order to sample the parameter space. Fortunately, there are many efficient maximization algorithms available in the literature (amoeba, Metropolis algorithms, simulated annealing, etc.). We will explore this in a future work. In the next sections of this paper, we will consider, for simplicity, the case of two channels, but we would like to remark that the extension to**

*f**N*> 2 frequencies is straightforward.

### 2.5 Optional parametrization of vector *g*

Another way to reduce the number of degrees of freedom of vector * g* is to find a suitable parametrization for it. For example, the power-law relationship

*I*(ν) is the flux at frequency ν, ν

_{0}is a frequency of reference,

*I*

_{0}is the flux at that frequency of reference and γ is the spectral index. This equation is widely used in the literature. By using equation (15), the reference flux

*I*

_{0}can easily be related to the reference amplitude

*A*of the sources and the number of degrees of freedom is just one, the spectral index γ.

However, the parametrization of vector * g* has its own risks. For instance, it is known that equation (15) is valid only as an approximation for any given frequency interval and that its validity decreases as the size of the interval grows. If we choose to follow the parametric approach we know for sure that the results will be less and less accurate as the number

*N*of channels grows, especially if the separation between frequency channels is large. On the other hand, equation (15) is always exact if

*N*= 2. Therefore, we can safely use the parametrization (15) for

*N*= 2 channels, without loss of generality. Since the number of degrees of freedom is one either if (15) is used or not, the use of the parametrization is irrelevant in this case. However, we may be interested in using it for historical, didactic and practical motivations. For example, equation (15) is useful to express the physical properties of the sources in terms of their (steep, flat, inverted, etc.) spectral index.

We would like to remark that frequency dependence parametrization, in the form of equation (15) or any other way, may or may not be useful in some cases, but it is not essential for the method we propose in this paper.

## 3 SIMULATIONS

In order to illustrate the MMF method described above and to compare the MMF multifrequency approach with the single-frequency approach, we have performed a set of basic, yet realistic, simulations. We consider the case of two frequency channels (*N*= 2). Generalization to more frequency channels is possible, as discussed in Section 2.3, but we choose to keep things simple in this paper.

For this example, we take the case of the *Planck* mission (Tauber 2005). We will consider the 44 and 100 GHz *Planck* channels. The choice of the pair channels is not essential. Any other pair of channels would have served the same for this exercise. This particular choice allows to study the case of radio sources in two, not adjoining, channels with different instrumental settings: the 44 GHz channel belongs to the Low-Frequency Instrument of *Planck* and the 100 GHz channel belongs to the High-Frequency Instrument.

For the simulations, we have used the *Planck* Sky Model2 (PSM; Delabrouille et al., in preparation), a flexible software package developed by *Planck* Working Group 2 (WG2) for making predictions, simulations and constrained realizations of the microwave sky. The simulated data used here are the same as in Leach et al. (2008), where the characteristics of the simulations are explained in more detail. Maps are expressed in (Δ*T*/*T*), thermodynamic units. Simulations include all the relevant astrophysical components: the CMB sky is based on a Gaussian realization assuming the *WMAP* best-fitting *C*_{ℓ} at higher multipole; Galactic emission is described by a three component model of the interstellar medium comprising free–free, synchrotron and dust emissions. Free–free emission is based on the model of Dickinson, Davies & Davis (2003) assuming an electronic temperature of 7000 K. The spatial structure of the emission is estimated using an Hα template corrected for dust extinction. Synchrotron emission is based on an extrapolation of the 408 MHz map of Haslam et al. (1982) from which an estimate of the free–free emission was extracted. A limitation of this approach is that this synchrotron model also contains any anomalous dust emission seen by *WMAP* at 23 GHz. The thermal emission from interstellar dust is estimated using model 7 of Finkbeiner, Davis & Schlegel (1999).

For the purely descriptive purposes of this example, we take eight different regions of the sky located at intermediate Galactic latitudes (four of them uniformly distributed across the 40° North Galactic latitude parallel and four of them distributed in the same way 40° South of the Galactic plane). For each region, we select a 512 × 512 pixel flat and square patch (at 44 and 100 GHz). Pixel size is 1.72 arcmin for the two frequencies. Therefore, each patch covers an area of 14.656 deg^{2} of the sky. Once the region has been selected, we add simulated EPS with a spectral behaviour described by equation (15). We take as frequency of reference ν_{0}= 44 GHz. The antenna beam is also taken into account: the full width at half maximum (*FWHM*) = 24 arcmin for the 44 GHz channel and *FWHM*= 9.5 arcmin at 100 GHz. Finally, we have added to each patch uniform white noise with the nominal levels specified for *Planck* (The Planck Collaboration 2006) and this pixel size.

We are interested in comparing the performance of the multifrequency approach with that of the single-frequency MF. In particular, we expect to be able to detect fainter sources with the MMF than with the MF. From recent works (López-Caniego et al. 2006; Leach et al. 2008), we know that in this kind of *Planck* simulations, the MF can detect sources down to fluxes ∼0.3 Jy (the particular value depends on the channel and the region of the sky). Here, we will simulate sources in the interval [0.1, 1.0] Jy plus a few cases, that will be described below, where even lower fluxes are necessary. Regarding the spectral index of the sources, according to González-Nuevo et al. (2008), most radio galaxies observed by *WMAP* at fluxes ∼1 Jy show spectral indices that lie in the range (−1.0, 1.4).

We sample the interesting intervals of flux and spectral index by simulating sources with fluxes at 44 GHz *I*_{0}={0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0} Jy and spectral indices γ={−1.0, −0.7, −0.4, −0.1, 0.2, 0.5, 0.8, 1.1, 1.4}. For each pair of values (*I*_{0}, γ), we have simulated 100 point sources. The point sources are randomly distributed in the maps (the same source is placed in the same pixel in both frequencies), with only one constraint: it is forbidden to place a source closer than *FWHM*_{44}/2 pixels from any other. In this way, we avoid source overlapping. Image borders are also avoided. For each set of 100 sources, we proceed in the following way: we randomly choose one among the eight patches we have and place 10 sources in it. Then, we randomly choose another patch (allowing repetition) and place the next 10 sources and so on. In total, we have simulated 9000 sources for this exercise (the additional simulations at fluxes below 0.1 Jy are not included).

## 4 RESULTS AND DISCUSSION

In order to compare the MF and the MMF, we use the same maps with both methods. It means that not only the maps but also the sources (their intrinsic fluxes and positions) are identical for the two filters. Each simulation is filtered separately with the MF (4) and the MMF (8). In order to have a better estimation of the power spectra, avoiding as much as possible aliasing effects, we implement a power spectrum estimator that uses the 2D Hann window (Jiang, Chen & Liu 2002).

We will compare the performance of the two methods in terms of the following aspects: spectral index estimation, source detection and flux estimation.

### 4.1 Source detection

Direct comparison of source detection between the MF and the MMF is not an obvious task because for *N*= 2 input images the MF produces two filtered images, whereas the MMF produces only one combined filtered map. The meaning of ‘detection’ in the MMF case is straightforward (once a detection criterion is chosen, a source is detected or not), but in the MF case the situation is not so clear. Imagine we decide to apply the same detection criterion to the two MF filtered images (which is not an obvious option), then for a given source we can have three different outcomes of the detection.

The source may be detected in both maps.

The source may be detected in only one of the maps.

The source may not be detected in any of the maps.

In the first two cases, we can obtain point source catalogues (with different number of objects, in principle, for the two frequencies), but only in the first case are we able to estimate the flux at the two frequencies and therefore the spectral index. In the case of the MMF, if the source is detected we automatically know the spectral index and we can use equation (15) to give the fluxes at the two frequencies.

Therefore, in the following lines we will distinguish between two different cases when we speak about detections with the MF. On the one hand, the intersection of the detections in the two channels gives us the objects that can be used for studying the spectral index distribution; on the other hand, the union of the two sets gives us the total number of objects that can be detected in at least one of the channels. For the MMF, both sets are the same by definition.

Regarding the detection criterion, for simplicity we will apply the same criterion to all the filtered maps: the widespread 5σ threshold. Note that the 5σ threshold corresponds to different flux values for different filters. However, in this paper, we will follow the standard 5σ criterion for simplicity.

Fig. 1 shows the real sources (in per cent) that we detect above a 5σ level detection whose intrinsic fluxes (values introduced by us in the simulations) in the reference frequency (*I*_{0}) are the corresponding values in the horizontal axis. Table 1 shows the flux at which we are able to detect at least the 95 per cent of the sources. We can observe several interesting aspects. The first one is the fact that the MMF improves the level of detection with respect to the MF level for all the values of γ we have inserted.

γ |
I
_{0(MMF)95 per cent}(Jy) |
I
_{0(MFi)95 per cent}(Jy) |
I
_{0(MFu)95 per cent}(Jy) |
||

1.4 | 0.3 | 0.9 | >1.0 | >1.0 | 0.9 |

1.1 | 0.2 | 0.8 | 1.0 | 1.0 | 0.8 |

0.8 | 0.2 | 0.9 | 1.0 | 0.9 | 0.7 |

0.5 | 0.2 | 0.8 | 0.6 | 0.8 | 0.6 |

0.2 | 0.2 | 0.9 | 0.5 | 0.9 | 0.5 |

−0.1 | 0.1 | 0.9 | 0.4 | 0.9 | 0.4 |

−0.4 | 0.1 | 0.9 | 0.3 | 0.9 | 0.3 |

−0.7 | 0.1 | 0.8 | 0.3 | 0.8 | 0.3 |

−1.0 | 0.05 | 0.8 | 0.2 | 0.8 | 0.2 |

γ |
I
_{0(MMF)95 per cent}(Jy) |
I
_{0(MFi)95 per cent}(Jy) |
I
_{0(MFu)95 per cent}(Jy) |
||

1.4 | 0.3 | 0.9 | >1.0 | >1.0 | 0.9 |

1.1 | 0.2 | 0.8 | 1.0 | 1.0 | 0.8 |

0.8 | 0.2 | 0.9 | 1.0 | 0.9 | 0.7 |

0.5 | 0.2 | 0.8 | 0.6 | 0.8 | 0.6 |

0.2 | 0.2 | 0.9 | 0.5 | 0.9 | 0.5 |

−0.1 | 0.1 | 0.9 | 0.4 | 0.9 | 0.4 |

−0.4 | 0.1 | 0.9 | 0.3 | 0.9 | 0.3 |

−0.7 | 0.1 | 0.8 | 0.3 | 0.8 | 0.3 |

−1.0 | 0.05 | 0.8 | 0.2 | 0.8 | 0.2 |

The second one is a natural selection effect. We detect more flat/inverted sources (γ∼ 0/negative values of γ) at low fluxes than steep ones (positive values of γ). Keeping in mind that the reference frequency ν_{0} is equal to 44 GHz, and according to equation (15), it can be seen that for γ > 0, the simulated sources appear less bright at 100 GHz (*I*_{100} < *I*_{44}). In these conditions, it is quite difficult to detect sources at 100 GHz, and for this reason we are not able to give the spectral indices for most of these point sources by means of the MF method when γ is strongly positive (for instance, γ≳ 1). On the other hand, for the opposite reason, we have added two additional bins (*I*_{0}= 25 and 50 mJy) for the cases γ≤−0.1 because the MF (and the MMF too) seems to perform better at 100 GHz for low values of γ.

Another aspect we have to remark is the similar trend of the detection curve for the MF at 44 GHz for all the γ values (see Fig. 1). The reason of this similarity is that the reference frequency is 44 GHz and the maps we have simulated at this frequency are the same, independently of γ, with only one exception: the position of the sources. This means that statistically they are equivalent (with the inherent fluctuations due to the variation in the positions of the sources and the patches randomly selected). We find a similar number of detections by means of the MF at 44 GHz, independently of γ. For this reason, it is difficult to characterize the spectral behaviour of the sources for *I*_{0}≲ 0.5 Jy because we do not have a high percentage of detected sources at 44 GHz below that value of *I*_{0}.

Additionally, we observe that the MMF is capable to detect sources with *I*_{0} < 0.1 Jy for γ≲−0.1. It is interesting to compare this with the MF which does not detect sources below 0.1 Jy in the conditions of this work. The method presented here, allows us to detect point sources whose *I*_{0} is too low to be detected with the traditional MF.

To summarize, we can say that the MMF improves the detection level. Specially, remarkable are the cases where the sources are near to be flat (central row of Fig. 1). At 100 GHz, the MF recovers the 100 per cent of the sources for *I*_{0} 0.4–0.6 Jy. Meanwhile, the MMF reaches this level for *I*_{0} 0.1 Jy. This particular case is really interesting because high-frequency surveys show that most of the sources have this spectral behaviour.

### 4.2 Spectral index estimation

One interesting quantity commonly used to describe the different subpopulations of sources as a function of its frequency behaviour is the spectral index. Normally, in monofrequency detection methods, the spectral indices are estimated using the detected fluxes in two different frequencies (equation 15). On the contrary, in our proposed method, we detect simultaneously the spectral indices (the one that maximize the SNR) and the fluxes only at the reference frequency. Then, the fluxes at the other frequencies are estimated using the usual relation (equation 15). In Fig. 2, we see how the spectral indices are recovered by means of the MMF and by the MF. In general, we observe that the MMF is able to recover the value of γ with more accuracy and less uncertainty than the traditional MF.

Another aspect is that the error bars increase when *I*_{0} is smaller. It seems logical because we have fainter sources and a smaller number of detections (see Fig. 1). Then, at *I*_{0}= 0.1 Jy, we can see that the estimation of γ is not as good as we wish because it has a great uncertainty. The main reason is that the SNR is close to the threshold level we have imposed.

In the case of the MF, the spectral indices are correctly estimated for *I*_{0}≳ 0.7 Jy at γ≤ 0.8. At higher values of γ, we find the same problem that we have mentioned before: there are a few detections at 100 GHz below 0.7 Jy (Fig. 1). Since the detected sources are close to the noise level, the fluxes recovered present an overestimation with respect to the input value due to the Eddington bias (Eddington 1913).

Finally, we can observe an interesting aspect of the MF. When we do not have the sufficient detections in at least one channel (the sources detected are below the ∼40 per cent of the total number of sources), the estimation of the spectral index is not good. For γ≥ 1.1, the sources at 100 GHz are much fainter and the number of detections at this channel is really small. Then, because of the Eddington bias, the flux at this frequency is overestimated, and consequently, the value of γ is underestimated. The Eddington bias explains as well the overestimation of γ in the other cases. The only difference is that now it is at 44 GHz where we have a smaller number of detections. If we also see Fig. 1, we observe that for values of *I*_{0}≲ 0.6 Jy, we are pretty close to the noise level. It means that the noise fluctuations in the maps produce an overestimation in the flux at 44 GHz (*I*_{0}) and, as a consequence, an overestimation in γ too. In summary, for γ= 1.4 and 1.1, the Eddington bias appears at 100 GHz (underestimation of the spectral index). For the rest of γ values, this bias appears at 44 GHz (overestimation of the spectral index).

### 4.3 Flux estimation

Fig. 3 shows the recovered flux at the reference frequency (44 GHz) for a given value of the spectral index. The error bars recovered with the MF are, in general, larger than the ones we obtain with the MMF. It is particularly notorious at small values of *I*_{0}, where the recovered values of the flux have a good agreement with respect to the input values, with small error bars.

In general, for all the values of γ that we have studied, the MMF is a suitable and effective tool to estimate the *I*_{0} of the sources. For the MF, we observe a good determination of *I*_{0} for input values above 0.7 Jy. For smaller values, *I*_{0} has a higher value than its real one. That is due to the Eddington bias at 44 GHz. At this frequency, in Fig. 1, we observe that for values smaller than 0.6 Jy, we only detect a ∼40 per cent of the total sources. That means that many of these sources are close to this noise level. For the correct estimation of *I*_{0} and the spectral index, we need a good detection of the sources at the two channels. For low values of *I*_{0}, the number of detected objects is small and we have few statistics.

The same conclusions are applicable to the estimation of the fluxes at the second frequency (100 GHz in our example), even if in the MMF case these fluxes are estimated using equation (15) with the consequent propagation of the detection errors.

### 4.4 Reliability

For academic purposes, in the previous sections, we have produced the simulations introducing 100 sources for each of the pair values of flux at reference frequency and spectral index (see Section 3), which simplifies the comparison between both filters for all the cases under study. However, it is well known that the number of sources per flux interval, the source number counts, is not constant (De Zotti et al. 2005; González-Nuevo et al. 2008) nor is the spectral index distribution (González-Nuevo et al. 2008; Sadler et al. 2008; Massardi et al. 2009). In order to study the performance of this multifrequency method under more realistic simulations, we produced a new set of simulations (100) with the following characteristics.

We used as a background the same eight regions described in the previous sections.

The sources were simulated with an almost Poissonian distribution (see González-Nuevo, Toffolatti & Argüeso 2005 for more details about the method) at 44 GHz, with fluxes that follow the source number counts model of De Zotti et al. (2005).

The fluxes at 100 GHz were estimated assuming random spectral indices from the González-Nuevo et al. (2008) distribution.

The point source maps were filtered with the same resolution as the background maps and randomly added to them.

There is also another interesting quantity commonly used in the study of the performance of a source detector: the number of spurious sources. Spurious sources are fluctuations of the background that satisfy the criteria of the detection method and therefore are considered as detected sources. It is clear that the best method will be the one that has the best detections versus spurious ratio. The maps are filtered using the MMF and the MF at both frequencies. We estimate the position and flux of the sources above 3σ level and, by comparing them with the input source simulations, we count the number of real and spurious sources that we are able to detect. We have also changed the detection level from 5σ to 3σ in order to increase the number of spurious sources to make the analysis.

In Fig. 4, we observe the number of real sources that both methods are capable to detect, whose intrinsic fluxes are higher than the corresponding value in the horizontal axis. As we can see at 44 GHz, MMF detects a higher number of real sources for fluxes below ∼0.4–0.5 Jy, being this difference very important at lower fluxes. At 100 GHz, we observe a similar behaviour, but in this case the differences between the MF and the MMF start at ∼0.2 Jy. If we observe Fig. 1, we note that the number of sources detected with the MF is higher at 100 GHz than at 44 GHz for values of the spectral index between 0 and 0.5. These values of γ, according to the model used to simulate the point sources in this section, are the most frequent ones. This gives us an idea about why the detection level of the MF is higher at 100 GHz.

In Fig. 5, the reliability of both methods at 44 and 100 GHz is compared. Reliability above a certain recovered flux is defined as *r*=*N*_{d}/(*N*_{d}+*N*_{s}), where *N*_{d} is the number of real sources above that flux, and *N*_{s} is the number of spurious sources above the same flux. At 44 GHz, we reach a ∼100 per cent of reliability at fluxes of ∼0.3 Jy with the MMF. However, the MF at this frequency reaches this level of reliability only for ∼1 Jy. On the other hand, at 100 GHz we obtain better levels of reliability. For instance, with the MMF we have at 0.1 Jy more than 95 per cent of reliability (∼0.3 Jy for the MF). Therefore, we can say that the MMF is more reliable than the MF, specially at lower fluxes.

Additionally, we can establish the flux for which we have the 5 per cent of spurious sources. Using the MF, we have found ∼0.5–0.6 and 0.25 Jy at 44 and 100 GHz, respectively, to be compared with 0.15 and <0.1 Jy in the case of the MMF.

Moreover, we make an additional plot where we represent, for both frequencies, the number of real sources detected versus the number of the spurious sources (Fig. 6). In this way, what we represent is the number of sources that a method detects, given a number of spurious sources. If we compare both plots, we can see that the curve of the MMF is always above the MF. It means that, when we have a fixed number of spurious detections, the MMF method detects more real sources.

Finally, we have to point out that the plots that we have introduced here are not directly comparable to Fig. 1, due to three basic but important differences.

A different way to simulate the point sources.

A different level of detection (in this case, a 3σ level).

In Fig. 1, we represent the number of sources with the corresponding flux in the horizontal axis. In the plots of this section, what we represent is the number of sources whose fluxes are higher than the corresponding value in the horizontal axis.

## 5 CONCLUSIONS

The detection of EPS in CMB maps is a challenge taking into account that they have to be removed to do a proper study of the cosmic radiation. But at the same time, it is of great interest to study their properties, spatial and spectral distributions, etc. For these reasons, we need suitable tools to detect and extract these sources. There are many filtering techniques that have been used in this context. In this work, we have compared the MF, one of the most studied techniques, with a multifrequency one based on the MMF. The great difference is that the latter takes into account information from all the channels of the same sky region in a simultaneous way. In particular, we show an example for *N*= 2.

The different tests that we have used have shown an improvement in the results obtained by the MMF with respect to the traditional MF. The number of detections is always higher when the MMF is used.

Another important aspect is to give a good estimation of the quantities that we have chosen to determine the sources, basically the spectral index and the flux at the reference frequency. In both cases, we can see that the MMF improves the results obtained with the MF: the values are close to the input values with smaller error bars (with one exception, the determination of the spectral index for *I*_{0}≲ 0.5 Jy at positive values of the input γ). This is a significant fact in order to be able to detect and study properly these kind of sources.

Additionally, we have made a set of more realistic simulations in order to study and compare both filters in terms of the spurious sources. We have also changed the threshold detection from 5σ to 3σ to find more spurious sources and make a more complete statistical analysis. First of all, we compare the number of real detections that we obtain with both techniques at 44 and 100 GHz. By comparing the plots of Fig. 4, we appreciate that, at lower fluxes, we detect more real sources with the MMF than with the MF. This aspect is more notorious at 44 GHz.

Moreover, the reliability of the MMF is much higher than the reliability of the MF for low fluxes (see Fig. 5). This difference is particularly important at 44 GHz, where the MF obtains similar values to the reliability of the MMF only for fluxes close to 1 Jy. At 100 GHz, the MF has the same reliability of the MMF only for fluxes greater than 0.3 Jy.

The last aspect that we have used to compare both methods is to look at the number of real sources that we have for a fixed number of spurious detections. The most efficient method is the one that has higher number of real detections for the same value of spurious detections. If we see Fig. 6, the best method is the MMF because its curves are always above the MF, i.e. if we take a number of spurious sources, the MMF recovers a larger number of real objects.

Finally, we have commented at Section (1) and the beginning of Section (2.3) the possibility of devising a filtering method (the MTXF) that does not use the frequency dependence of the sources altogether, totally independent of the frequency behaviour of the sources (flat, steep or inverted). This fact is significant in the sense that this filtering method is a robust technique for changes of *f*_{ν}. By contrast, it is necessary to impose the condition of orthonormalization of the matrix of the filters (see Herranz & Sanz 2008; Herranz et al. (2009) for more details), reducing the power of the method. On the contrary, the MMF is optimal in the sense of the SNR (see Section 2.3), but more complicated because we have to optimize another set of parameters (*f*_{ν}). For these reasons, the MTXF and the MMF are complementary.

The authors acknowledge partial financial support from the Spanish Ministry of Education (MEC) under project ESP2004-07067-C03-01 and the joint CNR-CSIC research projects 2006-IT-0037 and 2008IT0059. LFL acknowledges the Spanish CSIC for a JAE-Predoc fellowship and the hospitality of the Osservatorio Astronomico di Padova (INAF) during a research stay. Partial financial support for this research has been provided to JLS by the Spanish MEC and to JG-N by the Italian ASI (contracts Planck LFI Activity of Phase E2 and I/016/07/0 COFIS) and MUR. JG-N also acknowledges a researcher position grant at the SISSA (Trieste). ML-C acknowledges a post-doctoral fellowship from EGEE-III (FP7 INFSO-RI 222667). The authors are grateful to Raquel Fraga Encinas for reading carefully this paper and her style and grammatical suggestions and acknowledge the use of the PSM, developed by the Component Separation Working Group (WG2) of the Planck Collaboration.