An approach to detection of point sources in very highresolution microwave maps
^{1}
Chip Computers Consulting s.r.l., Viale Don L. Sturzo 82, S. Liberale di
Marcon,
30020
Venice,
Italy
email:
robertovio@tin.it
^{2}
ESO, Karl Schwarzschild strasse 2, 85748
Garching,
Germany
email:
pandrean@eso.org
^{3}
Centro de Astrofísica, Universidade do Porto,
Rua das Estrelas,
4150762
Porto,
Portugal
email: eramos@astro.up.pt; asilva@astro.up.pt
^{4}
Departamento de Física e Astronomia, Faculdade de Ciências,
Universidade do Porto, Rua do Campo
Alegre, 4169007
Porto,
Portugal
Received:
13
June
2012
Accepted:
17
May
2013
This paper deals with the detection problem of extragalactic point sources in multifrequency, microwave sky maps that will be obtainable in future cosmic microwave background (CMB) radiation experiments with instruments capable of very high spatial resolution. With spatial resolutions that can be 0.1−1.0 arcsec or better, the extragalactic point sources will appear isolated. The same also holds for the compact structures due to the SunyaevZeldovich (SZ) effect (both thermal and kinetic). This situation is different from the maps obtainable with instruments such as WMAP or Planck where, because of the lower spatial resolution (≈5−30 arcmin), the point sources and the compact structures due to the SZ effect form a uniform noisy background (confusion noise). The point source detection techniques developed in the past are therefore based on the assumption that all the emissions that contribute to the microwave background can be modeled with homogeneous and isotropic (often Gaussian) random fields and make use of the corresponding spatial power spectra. In the case of very highresolution observations, such an assumption cannot be adopted since it still holds only for the CMB. Here, we propose an approach based on the assumption that the diffuse emissions that contribute to the microwave background can be locally approximated by twodimensional loworder polynomials. In particular, two sets of numerical techniques are presented that contain two different algorithms each. The first set makes use of the a priori information about the spectral properties of CMB and SZ and is suited to detecting an extragalactic point source with a different spectrum for these emissions. In this set, one algorithm is a modification of the internal linear combination method, which is widely used in cosmology to extract the component of interest from a mixture of signals, and it is appropriate for extragalactic point sources with a known spectrum. The other one does not make use of this piece of information. The second set is tailored to detecting of extragalactic point sources with a similar spectrum to that of the CMB or SZ. Also in this set one algorithm is specific for extragalactic point sources with known spectrum whereas the other does not make use of this information. The performance of the algorithms is tested with numerical experiments that mimic the physical scenario expected for high Galactic latitude observations with the Atacama Large Millimeter/Submillimeter Array (ALMA).
Key words: methods: data analysis / methods: statistical / cosmic background radiation
© ESO, 2013
1. Introduction
The detection of extragalactic point sources in experimental microwave maps is a critical step in analyzing of the cosmic microwave background (CMB) maps. Besides the specific interest related to constructing of dedicated catalogs, these sources can, if not properly removed, have adverse effects on the estimation of the power spectrum and/or the test of Gaussianity of the CMB component. Much effort has been dedicated to multiple frequency maps of the same sky area, and many algorithms have been proposed (see Herranz & Sanz 2008; Herranz et al. 2012,and references therein). Apart from a recent Bayesian approach (Carvalho et al. 2009), most of them belong to two broad classes of techniques. The first class, suited to the extragalactic point sources with known spectra, is based on the NeymanPearson (NP) criterion that consists in maximizing the probability of detection P_{D} under the constraint that the probability of false alarm P_{FA} (i.e., the probability of a false detection) does not exceed a fixed value α (Kay 1998). The resulting algorithms are extensions of the classic matched filter (MF; Herranz et al. 2002; Ramos et al. 2011). The second class, appropriated to the extragalactic point sources with unknown spectra, is based on maximizing of the signaltonoise ratio (S/N) of the source intensity with respect to the underlying background (Herranz et al. 2009; Lanz et al. 2010). Both classes require the spatial power spectrum of the emitting components and therefore are based on the assumption that the emissions contributing to the microwave sky can be modeled by means of homogeneous and isotropic random fields.
With a spatial resolution worse than 5 arcmin typical of Planck and WMAP experiments, such an assumption is excellent for the CMB, good for the confusion noise due to the point sources and the SunyaevZeldovich (SZ) effect component, and locally barely acceptable for the diffuse Galactic emission. In the case of observation at very high spatial resolution, of order of 0.1−1.0 arcsec (that will be possible, for instance, with instruments as ALMA), this is no longer true. Indeed, since almost all of the extragalactic point sources but also of the compact structures due to the SZ effect will appear isolated, they can be thought of as the realization of (not necessarily stationary) shotnoise processes. Moreover, on small observing areas, it is realistic to expect that any other additional diffuse component (e.g. due to the Galactic emissions) is almost constant or slowly changing. This experimental scenario is different from the past or ongoing CMB experiments. New techniques of point source detection are therefore necessary.
From the consideration that the detection of extragalactic point sources is typically done on very small areas of the sky, where the contribution of the diffuse components can be approximated well by means of lowdegree, twodimensional polynomials, we propose an approach based on two sets of algorithms. The first set makes use of the a priori information about the spectral properties of CMB and SZ, and the second one does not exploit this piece of information. Each set contains two algorithms that are suited to detecting point sources with known and unknown spectra, respectively. The reason for two different sets of algorithms is that the use of the spectral properties of the CMB and SZ permits removal of the contribution from these components. As a consequence, the knowledge of their spatial power spectrum is no longer necessary, and at the same time it is possible to unambiguously distinguish between true extragalactic point sources and the compact structures due to SZ. The price is a reduced, or even null, detection capability for the point sources with spectra that are similar, or even identical, to that of CMB, SZ, or to a linear combination of these. The algorithms that do no make use of this piece of information do not suffer this limitation, but they cannot distinguish the true extragalactic point sources from the compact structures due to SZ.
The reason for two algorithms within each set is that those specialized for the extragalactic point sources with a specific spectrum are characterized by a greater detection capability. They are obtained by modifying the internal linear combination (ILC) method that in cosmology is used to extract a component of interest from the mixture of signals that contribute to the microwave sky emission (Eriksen et al. 2004; Hinshaw et al. 2007; Vio & Andreani 2008). Their main limitation is the necessity of multiple applications for detecting point sources with different spectra. Such a necessity is avoided by the algorithms that do not exploit this piece of information but at the price of less detection capability. The combined use of a couple of these algorithms belonging to different sets permits the detection of extragalactic point sources independently of their spectral emission and reduced contamination due to the SZ compact structures. The performances of both sets of algorithms is tested via numerical experiments based on simulated maps of high Galactic latitude that might be the area of interest of CMB high spatial resolution observations.
The paper is organized as follows. Section 2 introduces the mathematical framework and explains the algorithms used. Section 3 discusses practical problems related to the choice of the experimental parameters. Numerical experiments are reported in Sect. 4, while conclusions are summarized in Sect. 5.
2. Formalization and solution of the problem
By searching for a single extragalactic point source in a small area of sky, the microwave emission can be modeled with bidimensional discrete patches , each of them containing N_{p} = N_{p1} × N_{p2} pixels, corresponding to N_{f} different observing frequencies (channels), with the form (1)Here, is the contribution of the extragalactic point source at the ith frequency, , , and are the backgrounds due to CMB, extragalactic emission, and some other possible diffuse component (e.g. Galactic emission), respectively, and is the instrumental noise. In this model, the contribution of the extragalactic point sources is assumed in the form (2)with “a_{i}” the intensity of the source at the ith channel. According to Eq. (2), and without loss of generality, all the sources are assumed to have the same profile ℱ independently of the observing frequency. In practical applications, this is not true. However, it is possible to meet this condition by convolving the images with an appropriate kernel (see below).
For computational reasons that soon will become evident, it is useful to convert the twodimensional model (1) into the onedimensional form (3)Here, , with VEC [H] the operator that transforms a matrix H into a vector by stacking its columns one underneath the other. Something similar holds for the other quantities.
2.1. Detection with ILC background removal
One classical solution to deal with many maps of the same sky area taken at different observing frequencies consists in a linear composition by means a set of weights w = [w_{1},w_{2},...,w_{Nf}] ^{T}. In this way it is possible to work with a single map given by (4)where X = (x_{1},x_{2},...,x_{Nf}) is an N_{p} × N_{f} matrix. The obvious question is how to fix such weights.
Before proceeding, it is necessary to take into account that there is an a priori information about the various components in Eq. (1). In particular,

For each observing frequency i, the spectra of and are known with good accuracy. Moreover, the spatial distribution of these components is independent of the observing frequency;

has a diffuse spatial distribution with a coherence scale of about 10 arcmin, which is much greater than that of the point sources. It is reasonable to expect that something similar could hold for ;

Noises can be reasonably assumed as given by the realization of independent Gaussian whitenoise processes with standard deviation .
The first point implies that, for a given observing frequency, the emission due to the CMB and the SZ components can be obtained from a rescaled linear mixture of two templates (i.e., maps that do not depend on frequency) c and z, respectively. This implies that the cumulative contribution b_{i} = c_{i} + z_{i} of these components is given by the ith column of matrix (5)where M is an N_{e} × N_{f} matrix usually indicated with the term of mixing matrix. In the present context, the number of emission mechanisms is N_{e} = 2 since the kinetic SZ emission has the same spectrum as for the CMB. For this reason, with c_{i} we indicate the CMB plus the kinetic SZ emission from now on. The second point implies that within a small area centered on an extragalactic point source, the CMB and any other diffuse emission vary very little. This suggests that, for any patch with − N_{j} ≤ j ≤ N_{j} and − N_{k} ≤ k ≤ N_{k} (N_{p1} = 2N_{j} + 1, N_{p2} = 2N_{k} + 1), these emissions can be safely approximated by a lowdegree, twodimensional polynomial of degree m(6)where { α_{l} } are real coefficients, whereas q and r are integer numbers permuted accordingly.
2.1.1. Point sources with known spectrum
Starting from these considerations and adopting the criterion of the S/N maximization for a given extragalactic point source with emission spectrum , where “a” is to be estimated and is fixed, the weights w can be computed through the maximization of the quantity (7)with “∥ . ∥” the Euclidean norm, under the constraints Because of the constraint (8), which forces weights w to preserve the intensity a, the numerator of Eq. (7) is a constant. Therefore, the maximization of the SNR is obtained by the minimization of the denominator. The rationale behind this approach is that if each of the maps x_{i} contains the contribution of a point source with shape ℱ and of smooth component g_{i}, approximable by means of a twodimensional polynomial, then the same has to hold for their linear combination x = Xw. The denominator x−Lq thus represents the residuals of the leastsquares fit to the composite map x of a model where a central point spread function (PSF) is superimposed to a bivariate polynomial background. The weights w are computed in such a way as to minimize the standard deviation of these residuals under the constraints (8) and (9). In Eq. (7), q = (a,α^{T})^{T} is an array with size N_{c} = [(m + 1)(m + 2)/2] + 1, α the coefficients of the twodimensional polynomial, whereas L is an N_{p} × N_{c} matrix with the form L = [f,P] with f = VEC [ℱ] and P the N_{p} × (N_{c} − 1) design matrix corresponding to the leastsquares fit of a twodimensional polynomial^{1}. The constraint (9) forces the contribution of the CMB and SZ components to the final map x to zero, whereas the contamination g is removed by means of the twodimensional polynomial. The quantities w and q are unknown and have to be estimated. The maximization of the SNR with the constraints (8)−(9) can be written in the form (10)with , λ an N_{e} + 1 array of Lagrange multipliers, and e_{1} an N_{e} + 1 array of zeros except for the first element, which is “1”.
This method is a modification of the constrained ILC by Remazeilles et al. (2010). We call it modified multiple ILC (MMILC). The basic idea is that, if in the center of the selected patch there is a point source, then the value of “a” should exceed a threshold due to noise. After some algebra, one obtains that the solution of problem (10) is given by the system of equations (11)C_{XX} = X^{T}X, C_{XL} = X^{T}L and C_{LL} = L^{T}L, which provides where (14)with I the identity matrix, and (15)One interesting characteristic of solution (11) is that it does not require knowing the noise level of each map, a quantity that often can only be roughly estimated.
A useful insight into how MMILC works can be obtained if problem (10) is recast in the form^{2}(16)where A = (I−L(L^{T}L)^{1}L^{T})X. Since matrix A is the orthogonal projection of X onto the nullspace of L^{T}, MMILC can be seen to sequentially perform a leastsquares fit of the PSF overlapping a twodimensional polynomial background on the original data for each frequency, followed by a constrained ILC on the residuals.
2.1.2. Point sources with unknown spectrum
The main limitation of MMILC is that it only works optimally for a specific emission spectrum . This assumption can be relaxed by converting the maximization of the SNR with the constraints (8)−(9) in the leastsquares minimization of the quantity (17)with the constraints The constraint (18) is set to avoid the trivial solution w = 0. In this way, problem (10) is converted into (20)After some algebra, it is possible to see that q is again given by Eq. (13) but with w the solution of the eigenvalue problem (21)where (22)and C_{MM} = MM^{T}. The searched for w is given by the eigenvector of H that minimizes quantity S. Presently, the only method that we can suggest is to insert each eigenvector in Eq. (17) and to check numerically which of them provides the smallest S. This is because we have ascertained that there are situations where the eigenvector corresponding to the lowest eigenvalue of H (a criterion typical of the leastsquares problems) does not work. Indeed, matrix H is not symmetric, and it cannot be expected to have any particular property. We call this method nonparametric MMILC (NPMMILC).
Although not specifically optimized for a particular , the results provided by NPMMILC depend on the emission spectrum of the point source. Indeed, if like MILC, this method too is interpreted as a sequential leastsquares fits followed by a constrained ILC, it can be understood that, in the case of a point source with a large amplitude in maps with a low S/N and a small amplitude in maps with a high S/N, this results in a reduced detection capability. A simple procedure for avoiding this problem consists of applying NPMMILC not to all maps but only to those for which the best S/N for the point source is expected. The choice can be based on . This means to use the a priori information on the emission spectrum in a different way from MMILC.
2.1.3. Detection procedure
When searching for extragalactic point sources with MMILC or NPMMILC in a given set of maps, the procedure consists in fixing the size (2N_{j} + 1) × (2N_{k} + 1) of a window that is made to slide, pixel by pixel, across the area of interest. At the end of this procedure a single map is obtained containing the estimated values of “a” for each pixel. Now, the question is to fix the detection threshold below which a given value of “a” is supposed to be only due to noise. In this respect, the direct use of solutions (12)−(14) and (21) is difficult. For this reason, two different procedures are suggested,

1.
For MMILC it is set a = 0 if a ≤ kσ_{L}, where k is a constant factor (typically k = 4,5), , σ_{n} = (σ_{n1},σ_{n2},...,σ_{nNf})^{T}, and is the first entry of matrix (L^{T}L)^{1}. This operation corresponds to estimating the standard deviation σ_{a} of “a” for a fixed w. Such an approach has the advantage that matrix (L^{T}L)^{1} can be computed only once since it is the same for all the patches. But, it has the disadvantage that the standard deviations of the noises { n_{i} } have to be known in advance. In the case of NPMMILC, a similar procedure holds, but it is necessary to take the absolute value of a. Indeed, it is readily verified that, with q given by Eq. (13), a change of sign of w does not modify the value of S;

2.
It is set a = 0 if a ≤ kσ_{map}, for MMILC, and  a ≤ kσ_{map}, for NPMMILC. Again, k is a constant factor and σ_{map} is the standard deviation of the entries in the final map. This is an unsophisticated approach, but it does have the advantage of not requiring the standard deviation of the noise in each patch, a quantity usually known only roughly.
Before concluding this section, we underline that the number of rows N_{e} = 2 of the mixing matrix M comes from our interest in exploring the situation in which the extragalactic component consists of secondary anisotropies of the CMB. In particular, we have only considered the SZ effect (both thermal and kinetic), which is the strongest one in galaxy clusters, groups of galaxies, and in protoclusters (i.e., Birkinshaw 1999). However, if the information is available for one or more additional components, then it is sufficient to update M and the same solutions (12)−(14) and (21) and hold for N_{e} = 3 or greater. Similarly, if one decides to remove only one component via ILC, either CMB or SZ, then it is sufficient to eliminate the appropriate row from matrix M and set N_{e} = 1 in the solution. In this way, however, the drawback is that the remaining component has to be removed by means of the twodimensional polynomial. This could be a necessary operation in the case of noisy maps (see below).
2.2. Detection without ILC background removal
The MMILC and NPMMILC detection techniques are potentially quite effective, however they suffer from two main drawbacks:

1.
To remove the CMB and SZ components, one or more of the weights in w have to be negative. As a consequence, since in the final map a = a^{T}w and , “a” is given by the sum of both positive and negative values whereas, σ_{map} is given by the sum of positive values alone. In other words, the background subtraction reduces the S/N with respect to a simple sum of the maps. The situation worsens when the emission of an extragalactic point source has a spectrum similar to that of the CMB or of the SZ since “a” will tend to zero;

2.
If is an array such that , i.e. belongs to the nullspace of M (i.e., it is given by the linear combination of the column of M), then the system (11) does not have any useful solution since, when , both constraints and Mw = 0 are satisfied, but it happens that the estimate a is such that its expected value is zero.
For this reason, to detect extragalactic point sources with belonging to the nullspace of M, the above procedures have to be adapted to work without the ILC removal of the CMB and the SZ components. Again, two different algorithms are presented.
2.2.1. Point sources with known spectrum
The case of point sources with known spectra can be easily obtained from problem (10) through the substitutions and e_{1} = 1: (23)with solution given by (24)The explicit solution for w and q is given by Eqs. (12)−(14) with . Detection is still carried out as explained in Sect. 2.1.3. With this method, which we call modified ILC (MILC), the CMB and SZ emissions are not removed through the use of the mixing matrix M, but rather by exploiting that the CMB, part of the SZ, and any other component with a diffuse spatial distribution can be removed through the polynomial approximation. As a consequence, the only contribution in the final map beyond that of the extragalactic point sources is the compact component of the SZ (both thermal and kinetic), and this is an unavoidable problem. Without additional information, it is impossible to separate an SZ emission with pointlike spatial distribution from a genuine extragalactic point source. In the case of SZ emission with more extended structures, a possible solution consists of checking if their spatial distribution is compatible with the PSF ℱ. This issue, however, is beyond the scope of the present work.
2.2.2. Point sources with unknown spectrum
If in the minimization of the quantity S as given in Eq. (17) the constraint (19) is relaxed, the nonparametric version of MILC is obtained (NPMILC). It is easily verified that for this problem too, q is given by Eq. (13), but now w is the solution of the eigenvalue problem Hw = γw with (25)As for NPMMILC, the searched w is given by the eigenvector of H that minimizes quantity S, and detection is carried out as explained in Sect. 2.1.3. Also, NPMILC suffers the same dependence on as NPMMILC.
3. Practical uses
In this section we discuss some practical problems and how they can be addressed. The first is related to the degree m of the polynomial used to approximate the background. Obviously, the smaller the sky area of interest the lower the degree of the polynomial. For example, considering that the CMB has a coherence scale of about 10 arcmin, it can be reasonably expected that with a resolution of 0.1−1.0 arcsec a firstdegree polynomial is a good choice. The second question is related to the sizes N_{j} and N_{k} of the patch for testing for the presence of a point source. Two competing requirements arise: on the one hand, N_{j} and N_{k} must be as large as possible to reduce errors in estimating the polynomial parameters, on the other, a small size implies that the approximation of the background with a lowdegree polynomial is a more reliable operation and, at the same time, that the probability two or more sources being in the same patch is low. For illustrative purposes, Fig. 1 shows the standard deviation σ_{a} of the estimated intensity a as provided by MILC in the case of a point source with a Gaussian profile and a dispersion σ_{psf} equal to three pixels. A single map is considered where the background is given by a twodimensional onedegree polynomial, instrumental noise is Gaussian and white with standard deviation σ_{n}, and N_{p1} and N_{p2} are progressively increased. The true value of a is one in unit of σ_{n}. The decrease in σ_{a} is evident. Figure 2 shows the relationship between P_{D} and P_{FA} for different values of the ratio a/σ_{n}. These figures clearly show that N_{j} and N_{k} lying in the range 3σ_{psf}–5σ_{psf} is a reasonable compromise.
Fig. 1 Standard deviation σ_{a} of the estimated intensity a as provided by MILC in the case of a point source, with a Gaussian profile with dispersion σ_{psf} equal to 3 pixels, as a function of the sizes N_{j} = N_{k} of the searching patch. Here, a single map is considered with a background given by a twodimensional, onedegree polynomial. Instrumental noise is Gaussian and white with standard deviation σ_{n}. The true value of “a” is 1 in units of σ_{n}. 

Open with DEXTER 
As shown in the Appendix A, that the exact PSF ℱ could not be known or that some of the point sources could be overlapping and/or could have an extended shape, have no important consequences. However, another issue arises because the shape of the PFSs changes with observing frequency in practical applications. Widespread practice is to convolve maps with a suited kernel function in order to get a common spatial PSF ℱ for all the frequencies. This operation has the beneficial effect of reducing the standard deviation of the instrumental noise, but at the same time it introduces a spurious spatial correlation in it. Actually, even if neglected, this is not expected to be critical since both MILC and MMILC are linear techniques, and the only consequence is a reduction of the efficiency of the leastsquares estimate of the coefficients q (i.e., the estimate is unbiased but with greater variance). Something similar is also expected for NPMILC and NPMMILC that represent the solution of a linear leastsquares problem with a quadratic constraint (i.e., both the quantity to minimize and the constraint are smooth functions). In other words, given the abovementioned reduction in the standard deviation of the noise, this spurious correlation is not expected to have critical consequences. This is especially true if one takes into account that there are other and more important approximations that make the analysis of data less rigorous (e.g., often the level of instrumental noise is only roughly known).
Fig. 2 Relationship between the probability of detection, P_{D}, vs. the probability of false alarm, P_{FA} for the case shown in Fig. 1 but for different values of the ratio a/σ_{n}. Note the different scale used for the abscissa in the bottomright panel. 

Open with DEXTER 
Fig. 3 Comparison of the spectrum of the point sources used in the numerical experiments of Figs. 4−7 with that of CMB and the thermal SZ. All spectra have been normalized in such a way to have absolute intensity equal to 1 at 90 GHz. Here, a_{1} and a_{2} indicate the point source with spectrum, respectively, similar to and unlike that of CMB used in the numerical experiment (see next figures). 

Open with DEXTER 
Fig. 4 Simulations of a sky region at high Galactic declination at the ALMA observing frequencies. 20 randomly distributed point sources with the same intensity have been added. Here, the point sources have a spectrum similar to that of CMB (see text and curve a_{1} in Fig. 3). The PSFs are assumed to be Gaussian with a standard deviation of 3 pixels. Noise is Gaussianwhite with standard deviation set to 0.12 time the standard deviation of the values in the corresponding noisefree maps. All of the point sources have the same intensity set to 1.7 times the standard deviation of the noise. The two bottomright panels show the simulated point sources and their position on the 950 GHz map. 

Open with DEXTER 
Fig. 5 Results provided by MILC, MMILC, NPMILC, and NPMMILC when applied to the maps in Fig. 4. The detection threshold has been set to 4σ_{L} for all algorithms except for NPMMILC for which a value of 3.5σ_{L} has been adopted (see text) and the background has been approximated by a twodimensional polynomial of degree one. The top and bottom left panels clearly show that both MMILC and NPMMILC are not able to retrieve point sources in this case. This happens because the spectrum of the point sources has a frequencydependence similar to the CMB and SZ, and therefore the subtraction process gets rid of all of them. The top and bottom right panels show that both MILC and NPMILC, on the contrary, retrieve all sources and the SZ pointlike emissions, because it subtracts the underlying diffuse component with the polynomial approximation. For NPMILC a greater noise contamination of the detection map is evident. Also notice the detection of two couples of overlapping point sources in the bottomright and middleright part of the map. 

Open with DEXTER 
A final question regards whether in practical applications it is more convenient to use the MILC and MMILC algorithms or the NPMILC and NPMMILC ones. Indeed, MILC and MMILC work optimally only for a specific emission spectrum , a feature common to other detection techniques such as the matched multifilter (Herranz et al. 2012). In principle, this is not a critical question. It is sufficient to apply the detection algorithm to a set of prefixed obtained by grouping sources in broad families – radio flat, radio steep, dusty galaxies of a certain type, etc. – and defining average spectral laws for each family (Ramos et al. 2011; Herranz et al. 2012). Such an approach is viable since MILC and MMILC are fast algorithms, and they require the numerical solution of linear systems containing no more than a few of tens of linear equations. Moreover, as shown in Ramos et al. (2011), where a version of MMILC without background subtraction is applied to highGalactic latitude WMAP maps, strong degradation of the detection capability has to be expected only if the spectrum of the point sources is quite different from the one for which the MMILC algorithm has been optimized. Of course, this kind of problem can be avoided using NPMILC and NPMMILC. However, since they are not optimized for specific emission characteristics, the price is a lower detection capability for specific spectra. Given the inexpensive computational cost of the four algorithms, the best choice is to try all of them and check the results.
Fig. 6 Simulations of a sky region at high Galactic declination at the ALMA observing frequencies, with 20 randomly distributed point sources with the same intensity added. In this case the point sources have a spectrum given by curve a_{2} in Fig. 3. The PSFs are assumed to be Gaussian with a standard deviation of 3 pixels. Noise is Gaussianwhite with standard deviation set to 0.12 times the standard deviation of the values in the corresponding noisefree maps. The PSFs are assumed to be Gaussian with a standard deviation of 3 pixels. The two bottomright panels show the simulated point sources and their position on the 950 GHz map. 

Open with DEXTER 
Fig. 7 Results provided by MILC, MMILC, NPMILC and NPMMILC when applied to the maps in Fig. 6. The detection threshold has been set to 4σ_{L} for all algorithms except for NPMMILC for which a value of 3.5σ_{L} has been adopted (see text), and the background has been approximated by a twodimensional polynomial of degree one. The top and bottom left panels clearly show that in this case both MMILC and NPMMILC miss only one point source and get rid of the SZ pointlike emissions. A greater noise contamination in the detection map of NPMMILC is also evident. The top and bottom right panels show that MILC and NPMILC retrieve all point sources but also the SZ pointlike emissions. This is a consequence of subtracting the underlying diffuse component with the polynomial approximation. Also for NPMILC a greater noise contamination of the detection map is evident, and notice the detection of two couples of overlapping point sources in the bottomright and middleright part of the map. 

Open with DEXTER 
Fig. 8 As in Fig. 7 with the difference that maps has not been thresholded. This is only to show that the CMB and the extended SZ components have been effectively removed by the polynomial approximation of the background. 

Open with DEXTER 
4. Numerical experiments
To support the arguments presented above, we present some numerical experiments here with simulated maps at high Galactic latitude (where the Galactic contamination is negligible) that is the region of interest for future CMB experiments. Since realistic experimental conditions are not yet available, such simulations are only presented for illustration.
We produced small sky patches of 0.86 deg^{2} at 3′′ angular resolution with several components, namely, the CMB and the SunyaevZel’dovich effects (SZ), both kinetic and thermal. To produce these maps we used hydrodynamic/Nbody simulations with cosmological parameters that are consistent with WMAP parameters for a flat Universe and standard ΛCDM model, with an equation of state for the dark energy component of w = −1. The adopted present time density parameters expressed in terms of the critical density are (Ω_{cdm},Ω_{Λ},Ω_{b}) = (0.256,0.7,0.044), a dimensionless Hubble constant of h = 0.71, and a mean CMB temperature of T = 2.725 K. Adiabatic initial conditions are assumed, a spectral index of n_{s} = 1, and full reionization at redshift 7.
For the present epoch, we considered a normalization power spectrum of σ_{8} = 0.9 and a shape parameter of Γ = 0.17. The CMB component is produced with the CAMB code (Lewis et al. 2000) to obtain the linear CMB power spectrum. The fullsky CMB temperature anisotropy map was generated with the HEALPix software (Górski et al. 2005) with Nside = 8192. From this map a small sky region was extracted with an area of about 0.86 deg^{2} around the equator, projected on a squared map. Details about the simulations of the SZ effect components can be found in da Silva et al. (2001) and Ramos et al. (2012). The frequencies chosen were 90, 150, 250, 330, 440, 675, and 950 GHz, which correspond to the ALMA receiver bands.
All components were coadded, resulting in a final map, ΔI_{CMB + SZ}/I, with a pixel size of 3 arcsec. We use the central part of the maps (300 × 300 pixels) and convolve them for each frequency with a Gaussian PSF with a dispersion of three pixels. To each map a whitenoise process has also been added with standard deviations σ_{ni} set to 0.12 time the standard deviation of the values of map itself. Finally, 20 randomly distributed point sources were included with a_{i} = 1.7σ_{ni}. In this way, maps with the same S/N are obtained. The values of σ_{ni} and a_{i} have been arbitrarily chosen to test algorithms under very bad operational conditions, but at the same time to obtain stable results (i.e. with different realizations of the noise process almost all the sources are correctly detected with no false detections).
The simulated experimental scenario corresponds to an adverse situation of rather low S/N and, since σ_{ni} increases with frequency, with a spectrum (see curve a_{1} in Fig. 3) that mimics that of the CMB plus SZ background (i.e. is close to the nullspace of M, or ). Figure 4 displays the simulated maps. We note that the point sources are not even visible, and they are by far exceeded by the SZ pointlike emission. Figure 5 shows the results obtained with the four algorithms presented above. For MILC, MMILC, and NPMILC the detection threshold has been set to 4σ_{L} whereas a value of 3.5σ_{L} has been used for NPMMILC. Background has been approximated by a twodimensional firstdegree polynomial. A sliding square window of 19 × 19 pixels has been adopted for the local search of point sources. As expected, the MMILC and NPMMILC do not work. On the other hand, MILC and NPMILC have effectively removed the CMB, as well as the SZ diffuse components, and correctly detected all the point sources in the map (see Fig. 8). However, many of the SZ pointlike emissions are also present.
The situation greatly improves if { a_{i} } = [2.00,1.00,0.60,0.40,0.20,0.06,0.02] ^{T} ⊙ σ_{n}, which simulates an emission spectrum decreasing with frequency (see curve a_{2} in Fig. 3). In this way . Figure 6 shows the maps corresponding to this case. As is visible in Fig. 7, now both MMILC and NPMMILC are able to correctly detect all the point sources except one, as well to completely remove the CMB and SZ contamination. We note the detection of two couples of overlapping point sources in the bottom right hand and middle right hand parts of the map. However, as expected, the NPMMILC map presents various false detections due to noise (the detection threshold has been set to 3.5σ_{L} in order to obtain the same number of true detections as MMILC). A point to stress is that, in the present high Galactic latitude experiment, the diffuse components { g_{i} } can be set to zero. As a result, in models (10) and (20) it should have been possible to set L = f (i.e., to avoid the polynomial approximation of the background). However, the use of the general model permits the stability of MMILC and NPMMILC to be tested for high level of noise.
5. Summary and conclusions
In this paper four algorithms, the modified ILC (MILC), the modified multiple ILC (MMILC), the nonparametricMILC (NPMILC), and the nonparametricMMILC (NPMMILC), have been presented for detection of extragalactic point sources in multifrequency, very highresolution CMB maps. In particular, MMILC and NPMMILC make use of the a priori information about the spectral properties of the CMB and SZ with MMILC tailored to detecting extragalactic point sources with a given spectrum that has to be different from that of these emissions. The other two algorithms are suited to detecting the extragalactic point sources with similar spectra to the CMB or SZ with MILC tailored to the specific spectra. The main property of the proposed algorithms is that they do not require any a priori knowledge of the statistical characteristics of spatial distribution of the diffuse emissions that contribute to the microwave background. Indeed, these can be locally approximated with a lowdegree, twodimensional polynomial. The two proposed sets of algorithms used in conjunction are effective in detecting extragalactic point sources independently of the spectral characteristics of their emission. Their potential performance has been illustrated with some numerical experiments.
If the degree is one, P = [δ_{1},δ_{2},1], whereas for a degree two P = [δ_{1} ⊙ δ_{1},δ_{2} ⊙ δ_{2},δ_{1} ⊙ δ_{2},δ_{1},δ_{2},1], where “⊙” represents the elementwise matrix multiplication (Hadamard product), 1 is a vector of ones, and δ_{1} = VEC [Δ_{1}], δ_{2} = VEC [Δ_{2}], where Δ_{1} is a matrix with 2N_{j} + 1 identical columns [−N_{k},−N_{k} + 1,...,0,...,N_{k} −1,N_{k}] ^{T}, whereas Δ_{2} is a matrix with 2N_{k} + 1 identical rows [−N_{j}, −N_{j} + 1,...,0,...,N_{j} −1,N_{j}].
Acknowledgments
E. P. Ramos is supported by grant POPHQRENSFRH/BD/45613/2008, from FCT (Portugal). E. P. Ramos and R. Vio thank ESO for its hospitality and support through the DGDF funding program.
References
 Birkinshaw, M. 1999, Phys. Rep., 310, 97 [NASA ADS] [CrossRef] [Google Scholar]
 Carvalho, P., Rocha, G., & Hobson, M. P. 2009, MNRAS, 393, 681 [NASA ADS] [CrossRef] [Google Scholar]
 da Silva, A. J. C., Barbosa, D., Liddle, A. R., & Thomas, P. A. 2001, MNRAS, 326, 155 [NASA ADS] [CrossRef] [Google Scholar]
 Eriksen, H. K., Banday, A. J., Górski, K. M., & Lilje, P. B. 2004, ApJ, 612, 633 [NASA ADS] [CrossRef] [Google Scholar]
 Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [NASA ADS] [CrossRef] [Google Scholar]
 Herranz, D., & Sanz, J. L. 2008, IEEE J. Selected Topics in Signal Processing, 5, 727 [NASA ADS] [CrossRef] [Google Scholar]
 Herranz, D., Sanz, J. L., Hobson, M. P., et al. 2002, MNRAS, 336, 1057 [NASA ADS] [CrossRef] [Google Scholar]
 Herranz, D., LópezCaniego, M., Sanz, J. L., & GonzálezNuevo, J. 2009, MNRAS, 394, 510 [NASA ADS] [CrossRef] [Google Scholar]
 Herranz, D., Argüeso, F., & Carvalho, P. 2012, Adv. Astron., 2012, 410965 [NASA ADS] [CrossRef] [Google Scholar]
 Hinshaw, G., Nolta, M. R., Bennett, C. L., et al. 2007, ApJS, 170, 288 [NASA ADS] [CrossRef] [Google Scholar]
 Kay, S. M. 1998, Fundamentals of Statistical Signal Processing: Detection Theory (London: Prentice Hall) [Google Scholar]
 Lanz, L. F., Herranz, D., Sanz, J. L., GonzalezNuevo, J., & LopezCaniego, M. 2010, MNRAS, 403, 212 [CrossRef] [Google Scholar]
 Lewis, A., Challinor, A., & Lasenby, A. 2000, ApJ, 538, 473 [NASA ADS] [CrossRef] [Google Scholar]
 Ramos, E. P. R. G., Vio, R., & Andreani, P. 2011, A&A, 528, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ramos, E. P. R. G., da Silva, A., J. C., & Liu, G. C. 2012, ApJ, submitted [Google Scholar]
 Remazeilles, M., Delabrouille, J., & Cardoso, J. F. 2011, MNRAS, 418, 467 [NASA ADS] [CrossRef] [Google Scholar]
 Vio, R., & Andreani, P. 2008, A&A, 487, 775 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Some additional questions
Fig. A.1 Central slice of the spectrum P(ν) of three bivariate, circularly symmetric, Gaussian PSFs with dispersion σ_{G} = 3,4,10, respectively. Frequency ν is in Nyquist units. 

Open with DEXTER 
The detection techniques presented in Sect. 2 are based on the assumptions that a) the true PSF ℱ is known; b) all the sources have a pointlike shape (i.e. the shape of the PSF); c) all the sources are isolated. In practical applications, however, the PSF has to be estimated, some sources can have an extended shape, and overlapping is possible. Here, we show that all this can be expected to have only secondary consequences. In this respect, for sake of ease in the formalism, we start by considering an experimental onedimensional signal x = af + n. Here, f is the onedimensional PSF and n an additional white noise. In practice, x is due to a point source with amplitude a embedded in noise. As is well known, the leastsquares fit estimate of a provides the solution (A.1)At this point, we point out that in detection problems the test statistic T(x) = f^{T}x (e.g. see Kay 1998) is used where f is the well known matched filter (MF). In other words, unless of a constant factor, the estimate is identical to the test statistic T(x). Indeed, the rhs of Eq. (A.1) is nothing else that the normalized correlation of x with f. This means that in the present context, the least squares fit and the matched filter are two equivalent techniques. This allows us to analyze the characteristic of the leastsquares fit using the body of the theory of the linear filters. Therefore, f can be interpreted as a linear, typically lowpass, discrete filter with N_{p} entries that is made to slide across an array with N_{x} ≫ N_{p}. In this way, for each element x_{i} one can obtain the quantity , where . By means of the discrete Fourier transform (DFT) it is possible to compute via (A.2)with and the DFT of x and the zero padded version of f, IDFT [.] the inverse DFT operator and symbols “”, “⊙” indicating the complex conjugate and the pointwise multiplication, respectively. This implies that the statistical characteristics of can be analyzed by means of the spectral properties of f. These considerations can be trivially extended to the twodimensional case.
Fig. A.2 Upper left panel: image of the original point source when the PSF is a bivariate, circularly symmetric, Gaussian PSF with dispersion σ_{G} = 4 pixels. Upper right panel: original point source added with a whitenoise with standard deviation equal to half the peak value of the source itself; noisy image filtered with an improper matched filter (MF) which has a Gaussian shape and σ_{G} = 3 pixels; bottom right panel: noisy image filtered with the correct MF. 

Open with DEXTER 
Fig. A.3 Upper left panel: image of the original point source when the PSF is a bivariate, circularly symmetric, Gaussian PSF with dispersion σ_{G} = 3 pixels. Upper right panel: original point source added with a white noise with standard deviation equal to half the peak value of the source itself. Bottom left panel: noisy image filtered with an improper matched filter (MF) that has a Gaussian shape and σ_{G} = 4 pixels. Bottom right panel: noisy image filtered with the correct MF. 

Open with DEXTER 
Fig. A.4 Upper left panel: image of an extended point source with a circularly symmetric bivariate Gaussian shape with dispersion σ_{g} = 10 pixels. Upper right panel: original point source added with a white noise with standard deviation equal to half the peak value of the source itself. Bottom left panel: noisy image filtered with an improper matched filter (MF), which has a Gaussian shape and σ_{G} = 3 pixels. Bottom right panel: noisy image filtered with the correct MF. 

Open with DEXTER 
Fig. A.5 Results obtainable with the matched filter (MF) in the case of two identical overlapping point sources with shape given by a bivariate, circularly symmetric, Gaussian PSF with σ_{G} = 3 pixels, when their peaks are 6 and 9 pixels apart. A white noise is added with standard deviation equal to half the peak value of the sources. 

Open with DEXTER 
Fig. A.6 Results obtainable with the matched filter (MF) in the case of two identical overlapping point sources with shape given by a bivariate, circularly symmetric, Gaussian PSF with σ_{G} = 3 pixels, when their peaks are 6 and 9 pixels apart. A firstdegree, twodimensional polynomial background and a white noise with standard deviation equal to half the peak value of the sources are added. 

Open with DEXTER 
Here, it is useful to present some simple examples to see that the three facts above are actually not important. Concerning point a), Fig. A.1 shows a central slice of the spectrum of three twodimensional, circularly symmetric, Gaussian PSFs with dispersion σ_{G} = 3 and 4 pixels, respectively. It is clearly visible that the value of σ_{G} determines the bandwidth of in the sense that the larger σ_{G}, the stronger the filtering action. This means that, if a point source with σ_{G} = 4 pixels is filtered assuming a PSF with σ_{G} = 3 (i.e. with an error of 25%), the only consequence is a slighter reduction of the noise with respect to what is obtainable with the correct PSF. In contrast, if the point source has σ_{G} = 3 and a PSF with σ_{G} = 4 is assumed (i.e. with an error of 33%), a stronger reduction of the noise is obtained. In this case, however, part of the signal of interest is also filtered out, again with a slighter reduction of the detection capability with respect to what is obtainable with the correct PSF. These two examples are shown in Figs. A.2, A.3.
Similar arguments hold for the point b) when the source has an extended shape. In this case, however, since the PSF is much “narrower” than the source, the noise filtering will be much less effective. However, as shown by the example in Fig. A.4, where an extended object with Gaussian shape, σ_{G} = 10, and a Gaussian PSF with σ_{G} = 3 pixels are considered, this does not mean that detection is not possible, but rather that it is less effective.
Finally, concerning point c), i.e. the case of two close point sources, the situation does not change very much. The only consequence is that the filtered point sources will appear more overlapped than the unfiltered ones (see Fig. A.5).
In the cases examined above, the presence of a smooth background p has not been considered. However, again, things can be expected not to change very much. Indeed, Eq. (A.1) becomes (A.3)Now, if instead of p the result of a leastsquares polynomial fit is used, then one obtains (A.4)With a smooth function. In other words, a map is produced where a smoothed point source is superimposed to a smoothed (presumably weak) background and noise is reduced. This is visible in Fig. A.6, which shows the results obtained for a situation similar to that of Fig. A.5 when a firstdegree, twodimensional polynomial background is present.
To understand why it is reasonable to assume that these facts are not critical for MMILC, it is useful to interpret this method as done in Sect. 2.1.1; i.e., a sequential leastsquares fit of a point spread function (PSF) overlapped to a twodimensional polynomial background on the original data for each frequency, followed by a constrained ILC on the residuals. Since, each of the fits is not very sensitive to the issues mentioned above, it is reasonable to assume that the same is valid for their linear combination. Something similar holds for MILC, NPMILC, and NPMMILC. These arguments can also be extended with very good approximation to the case that, as in Sect. 2, the background is computed for each position of the sliding detection window. Indeed, since the window is made to slide one pixel at a time, for close pixels the leastsquares fit is done using essentially the same data. The experiments concerning MILC and NPMILC in Figs. 4−7 confirm this.
All Figures
Fig. 1 Standard deviation σ_{a} of the estimated intensity a as provided by MILC in the case of a point source, with a Gaussian profile with dispersion σ_{psf} equal to 3 pixels, as a function of the sizes N_{j} = N_{k} of the searching patch. Here, a single map is considered with a background given by a twodimensional, onedegree polynomial. Instrumental noise is Gaussian and white with standard deviation σ_{n}. The true value of “a” is 1 in units of σ_{n}. 

Open with DEXTER  
In the text 
Fig. 2 Relationship between the probability of detection, P_{D}, vs. the probability of false alarm, P_{FA} for the case shown in Fig. 1 but for different values of the ratio a/σ_{n}. Note the different scale used for the abscissa in the bottomright panel. 

Open with DEXTER  
In the text 
Fig. 3 Comparison of the spectrum of the point sources used in the numerical experiments of Figs. 4−7 with that of CMB and the thermal SZ. All spectra have been normalized in such a way to have absolute intensity equal to 1 at 90 GHz. Here, a_{1} and a_{2} indicate the point source with spectrum, respectively, similar to and unlike that of CMB used in the numerical experiment (see next figures). 

Open with DEXTER  
In the text 
Fig. 4 Simulations of a sky region at high Galactic declination at the ALMA observing frequencies. 20 randomly distributed point sources with the same intensity have been added. Here, the point sources have a spectrum similar to that of CMB (see text and curve a_{1} in Fig. 3). The PSFs are assumed to be Gaussian with a standard deviation of 3 pixels. Noise is Gaussianwhite with standard deviation set to 0.12 time the standard deviation of the values in the corresponding noisefree maps. All of the point sources have the same intensity set to 1.7 times the standard deviation of the noise. The two bottomright panels show the simulated point sources and their position on the 950 GHz map. 

Open with DEXTER  
In the text 
Fig. 5 Results provided by MILC, MMILC, NPMILC, and NPMMILC when applied to the maps in Fig. 4. The detection threshold has been set to 4σ_{L} for all algorithms except for NPMMILC for which a value of 3.5σ_{L} has been adopted (see text) and the background has been approximated by a twodimensional polynomial of degree one. The top and bottom left panels clearly show that both MMILC and NPMMILC are not able to retrieve point sources in this case. This happens because the spectrum of the point sources has a frequencydependence similar to the CMB and SZ, and therefore the subtraction process gets rid of all of them. The top and bottom right panels show that both MILC and NPMILC, on the contrary, retrieve all sources and the SZ pointlike emissions, because it subtracts the underlying diffuse component with the polynomial approximation. For NPMILC a greater noise contamination of the detection map is evident. Also notice the detection of two couples of overlapping point sources in the bottomright and middleright part of the map. 

Open with DEXTER  
In the text 
Fig. 6 Simulations of a sky region at high Galactic declination at the ALMA observing frequencies, with 20 randomly distributed point sources with the same intensity added. In this case the point sources have a spectrum given by curve a_{2} in Fig. 3. The PSFs are assumed to be Gaussian with a standard deviation of 3 pixels. Noise is Gaussianwhite with standard deviation set to 0.12 times the standard deviation of the values in the corresponding noisefree maps. The PSFs are assumed to be Gaussian with a standard deviation of 3 pixels. The two bottomright panels show the simulated point sources and their position on the 950 GHz map. 

Open with DEXTER  
In the text 
Fig. 7 Results provided by MILC, MMILC, NPMILC and NPMMILC when applied to the maps in Fig. 6. The detection threshold has been set to 4σ_{L} for all algorithms except for NPMMILC for which a value of 3.5σ_{L} has been adopted (see text), and the background has been approximated by a twodimensional polynomial of degree one. The top and bottom left panels clearly show that in this case both MMILC and NPMMILC miss only one point source and get rid of the SZ pointlike emissions. A greater noise contamination in the detection map of NPMMILC is also evident. The top and bottom right panels show that MILC and NPMILC retrieve all point sources but also the SZ pointlike emissions. This is a consequence of subtracting the underlying diffuse component with the polynomial approximation. Also for NPMILC a greater noise contamination of the detection map is evident, and notice the detection of two couples of overlapping point sources in the bottomright and middleright part of the map. 

Open with DEXTER  
In the text 
Fig. 8 As in Fig. 7 with the difference that maps has not been thresholded. This is only to show that the CMB and the extended SZ components have been effectively removed by the polynomial approximation of the background. 

Open with DEXTER  
In the text 
Fig. A.1 Central slice of the spectrum P(ν) of three bivariate, circularly symmetric, Gaussian PSFs with dispersion σ_{G} = 3,4,10, respectively. Frequency ν is in Nyquist units. 

Open with DEXTER  
In the text 
Fig. A.2 Upper left panel: image of the original point source when the PSF is a bivariate, circularly symmetric, Gaussian PSF with dispersion σ_{G} = 4 pixels. Upper right panel: original point source added with a whitenoise with standard deviation equal to half the peak value of the source itself; noisy image filtered with an improper matched filter (MF) which has a Gaussian shape and σ_{G} = 3 pixels; bottom right panel: noisy image filtered with the correct MF. 

Open with DEXTER  
In the text 
Fig. A.3 Upper left panel: image of the original point source when the PSF is a bivariate, circularly symmetric, Gaussian PSF with dispersion σ_{G} = 3 pixels. Upper right panel: original point source added with a white noise with standard deviation equal to half the peak value of the source itself. Bottom left panel: noisy image filtered with an improper matched filter (MF) that has a Gaussian shape and σ_{G} = 4 pixels. Bottom right panel: noisy image filtered with the correct MF. 

Open with DEXTER  
In the text 
Fig. A.4 Upper left panel: image of an extended point source with a circularly symmetric bivariate Gaussian shape with dispersion σ_{g} = 10 pixels. Upper right panel: original point source added with a white noise with standard deviation equal to half the peak value of the source itself. Bottom left panel: noisy image filtered with an improper matched filter (MF), which has a Gaussian shape and σ_{G} = 3 pixels. Bottom right panel: noisy image filtered with the correct MF. 

Open with DEXTER  
In the text 
Fig. A.5 Results obtainable with the matched filter (MF) in the case of two identical overlapping point sources with shape given by a bivariate, circularly symmetric, Gaussian PSF with σ_{G} = 3 pixels, when their peaks are 6 and 9 pixels apart. A white noise is added with standard deviation equal to half the peak value of the sources. 

Open with DEXTER  
In the text 
Fig. A.6 Results obtainable with the matched filter (MF) in the case of two identical overlapping point sources with shape given by a bivariate, circularly symmetric, Gaussian PSF with σ_{G} = 3 pixels, when their peaks are 6 and 9 pixels apart. A firstdegree, twodimensional polynomial background and a white noise with standard deviation equal to half the peak value of the sources are added. 

Open with DEXTER  
In the text 