© 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 Neyman-Pearson (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 signal-to-noise 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 Sunyaev-Zeldovich (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) shot-noise 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 low-degree, two-dimensional 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 two-dimensional model (1) into the one-dimensional 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₁,w₂,...,w_Nf] ^T. In this way it is possible to work with a single map given by (4)where X = (x₁,x₂,...,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 white-noise 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 low-degree, two-dimensional 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 two-dimensional polynomial, then the same has to hold for their linear combination x = Xw. The denominator x−Lq thus represents the residuals of the least-squares 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 two-dimensional 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 least-squares fit of a two-dimensional polynomial¹. 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 two-dimensional 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₁ 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^TX, C_XL = X^TL and C_LL = L^TL, 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²(16)where A = (I−L(L^TL)^-1L^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 least-squares fit of the PSF overlapping a two-dimensional 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 least-squares 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 least-squares 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 (NP-MMILC).

Although not specifically optimized for a particular , the results provided by NP-MMILC depend on the emission spectrum of the point source. Indeed, if like MILC, this method too is interpreted as a sequential least-squares 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 NP-MMILC 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 NP-MMILC 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^TL)^-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^TL)^-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 NP-MMILC, 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 NP-MMILC. 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 two-dimensional 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 NP-MMILC 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^Tw 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: (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 point-like 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 (NP-MILC). 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 NP-MMILC, 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, NP-MILC suffers the same dependence on as NP-MMILC.

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 first-degree 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 low-degree 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 two-dimensional one-degree 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 Nj = Nk of the searching patch. Here, a single map is considered with a background given by a two-dimensional, one-degree polynomial. Instrumental noise is Gaussian and white with standard deviation σn. The true value of “a” is 1 in units of σn.

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 least-squares estimate of the coefficients q (i.e., the estimate is unbiased but with greater variance). Something similar is also expected for NP-MILC and NP-MMILC that represent the solution of a linear least-squares problem with a quadratic constraint (i.e., both the quantity to minimize and the constraint are smooth functions). In other words, given the above-mentioned 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, PD, vs. the probability of false alarm, PFA 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 bottom-right panel.

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, a1 and a2 indicate the point source with spectrum, respectively, similar to and unlike that of CMB used in the numerical experiment (see next figures).

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 a1 in Fig. 3). The PSFs are assumed to be Gaussian with a standard deviation of 3 pixels. Noise is Gaussian-white with standard deviation set to 0.12 time the standard deviation of the values in the corresponding noise-free maps. All of the point sources have the same intensity set to 1.7 times the standard deviation of the noise. The two bottom-right panels show the simulated point sources and their position on the 950  GHz map.

Fig. 5 Results provided by MILC, MMILC, NP-MILC, and NP-MMILC when applied to the maps in Fig. 4. The detection threshold has been set to 4σL for all algorithms except for NP-MMILC for which a value of 3.5σL has been adopted (see text) and the background has been approximated by a two-dimensional polynomial of degree one. The top and bottom left panels clearly show that both MMILC and NP-MMILC are not able to retrieve point sources in this case. This happens because the spectrum of the point sources has a frequency-dependence 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 NP-MILC, on the contrary, retrieve all sources and the SZ point-like emissions, because it subtracts the underlying diffuse component with the polynomial approximation. For NP-MILC a greater noise contamination of the detection map is evident. Also notice the detection of two couples of overlapping point sources in the bottom-right and middle-right part of the map.

A final question regards whether in practical applications it is more convenient to use the MILC and MMILC algorithms or the NP-MILC and NP-MMILC 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 high-Galactic 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 NP-MILC and NP-MMILC. 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 a2 in Fig. 3. The PSFs are assumed to be Gaussian with a standard deviation of 3 pixels. Noise is Gaussian-white with standard deviation set to 0.12 times the standard deviation of the values in the corresponding noise-free maps. The PSFs are assumed to be Gaussian with a standard deviation of 3 pixels. The two bottom-right panels show the simulated point sources and their position on the 950  GHz map.

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

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.

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² at 3′′ angular resolution with several components, namely, the CMB and the Sunyaev-Zel’dovich effects (SZ), both kinetic and thermal. To produce these maps we used hydrodynamic/N-body 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 σ₈ = 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 full-sky 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² 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 co-added, 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 white-noise 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₁ 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 point-like emission. Figure 5 shows the results obtained with the four algorithms presented above. For MILC, MMILC, and NP-MILC the detection threshold has been set to 4σ_L whereas a value of 3.5σ_L has been used for NP-MMILC. Background has been approximated by a two-dimensional first-degree polynomial. A sliding square window of 19 × 19 pixels has been adopted for the local search of point sources. As expected, the MMILC and NP-MMILC do not work. On the other hand, MILC and NP-MILC 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 point-like 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₂ 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 NP-MMILC 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 NP-MMILC 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 NP-MMILC 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 nonparametric-MILC (NP-MILC), and the nonparametric-MMILC (NP-MMILC), have been presented for detection of extragalactic point sources in multifrequency, very high-resolution CMB maps. In particular, MMILC and NP-MMILC 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 low-degree, two-dimensional 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.

Acknowledgments

E. P. Ramos is supported by grant POPH-QREN-SFRH/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

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.

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 point-like 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 one-dimensional signal x = af + n. Here, f is the one-dimensional 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 least-squares 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^Tx (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 least-squares fit using the body of the theory of the linear filters. Therefore, f can be interpreted as a linear, typically low-pass, 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 point-wise 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 two-dimensional 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 white-noise 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.

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.

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.

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.

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 first-degree, two-dimensional polynomial background and a white noise with standard deviation equal to half the peak value of the sources are added.

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 two-dimensional, 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 least-squares 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 first-degree, two-dimensional 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 least-squares fit of a point spread function (PSF) overlapped to a two-dimensional 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, NP-MILC, and NP-MMILC. 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 least-squares fit is done using essentially the same data. The experiments concerning MILC and NP-MILC 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 Nj = Nk of the searching patch. Here, a single map is considered with a background given by a two-dimensional, one-degree polynomial. Instrumental noise is Gaussian and white with standard deviation σn. The true value of “a” is 1 in units of σn.

In the text

	Fig. 2 Relationship between the probability of detection, PD, vs. the probability of false alarm, PFA 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 bottom-right panel.
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, a1 and a2 indicate the point source with spectrum, respectively, similar to and unlike that of CMB used in the numerical experiment (see next figures).
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 a1 in Fig. 3). The PSFs are assumed to be Gaussian with a standard deviation of 3 pixels. Noise is Gaussian-white with standard deviation set to 0.12 time the standard deviation of the values in the corresponding noise-free maps. All of the point sources have the same intensity set to 1.7 times the standard deviation of the noise. The two bottom-right panels show the simulated point sources and their position on the 950  GHz map.

In the text

Fig. 5 Results provided by MILC, MMILC, NP-MILC, and NP-MMILC when applied to the maps in Fig. 4. The detection threshold has been set to 4σL for all algorithms except for NP-MMILC for which a value of 3.5σL has been adopted (see text) and the background has been approximated by a two-dimensional polynomial of degree one. The top and bottom left panels clearly show that both MMILC and NP-MMILC are not able to retrieve point sources in this case. This happens because the spectrum of the point sources has a frequency-dependence 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 NP-MILC, on the contrary, retrieve all sources and the SZ point-like emissions, because it subtracts the underlying diffuse component with the polynomial approximation. For NP-MILC a greater noise contamination of the detection map is evident. Also notice the detection of two couples of overlapping point sources in the bottom-right and middle-right part of the map.

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 a2 in Fig. 3. The PSFs are assumed to be Gaussian with a standard deviation of 3 pixels. Noise is Gaussian-white with standard deviation set to 0.12 times the standard deviation of the values in the corresponding noise-free maps. The PSFs are assumed to be Gaussian with a standard deviation of 3 pixels. The two bottom-right panels show the simulated point sources and their position on the 950  GHz map.

In the text

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

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.
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.
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 white-noise 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.

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.

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.

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.
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 first-degree, two-dimensional polynomial background and a white noise with standard deviation equal to half the peak value of the sources are added.
In the text