text:
# Grand Unification signal from ultrahigh energy cosmic rays?

\[

## Abstract

The spectrum of ultrahigh energy (above $`10^9`$ GeV) cosmic rays is consistent with the decay of GUT scale particles. The predicted mass is $`m_X=10^b`$ GeV, where $`b=14.6_{1.7}^{+1.6}`$.

ITP-Budapest 560

\]

The interaction of protons with photons of the cosmic microwave background radiation (CMBR) predicts a sharp drop in the cosmic ray flux above the Greisen-Zatsepin-Kuzmin (GZK) cutoff around $`510^{19}`$ eV . The available data shows no such drop. About 20 events above $`10^{20}`$ eV were observed by experiments such as Akeno Giant Air Shower Array (AGASA) , Fly’s Eye , Haverah Park , Yakutsk  and HiRes . Future experiments, particularly Pierre Auger , will have a much higher statistics.

Usually it is assumed that at these energies the galactic and extragalactic (EG) magnetic fields do not affect the orbit of the cosmic rays, thus they should point back to their origin within a few degrees. Though there are clustered events  the overall distribution is practically isotropic , which usually ought to be interpreted as a signature for EG origin.

Since above the GZK energy the attenuation length of particles is a few tens of megaparsecs  if an ultrahigh energy cosmic ray (UHECR) is observed on Earth it must be produced in our vicinity (except for UHECR scenarios based on weakly interacting particles, e.g. neutrinos ). Sources of EG origin (e.g. AGN , topological defects  or the local supercluster ) should result in a GZK cutoff, which is in disagreement with experiments. It is generally believed  that there is no conventional astrophysical explanation for the observed UHECR spectrum.

An interesting idea suggested by refs. is that superheavy particles (SP) as dark matter could be the source of UHECRs. (Note, that metastable relic SPs were proposed  before the observation of UHECRs beyond the GZK cutoff.) In  EG SPs were studied. Ref.  made a crucial observation and analyzed the decay of SPs concentrated in the halo of our galaxy. They used the modified leading logarithmic approximation (MLLA)  for ordinary QCD and for supersymmetric QCD . A good agreement of the EG spectrum with observations was noticed in . Supersymmetric QCD is treated as the strong regime of the minimal supersymmetric standard model (MSSM). To describe the decay spectrum more accurately HERWIG Monte-Carlo was used in QCD  and discussed in supersymmetric QCD , resulting in $`m_X10^{12}`$ GeV and $`10^{13}`$ GeV for the SP mass in SM and in MSSM, respectively.

SPs are very efficiently produced by the various mechanisms at post inflatory epochs . Note, that our analysis of SP decay covers a much broader class of possible sources. Several non-conventional UHECR sources (e.g. EG long ordinary strings  or galactic vortons , monopole-antimonopole pairs connected by strings ) produce the same UHECR spectra as decaying SPs.

In this letter we study the scenario that the UHECRs are coming from decaying SPs and we determine the mass of this $`X`$ particle $`m_X`$ by a detailed analysis of the observed UHECR spectrum. We discuss both possibilities that the UHECR protons are produced in the halo of our galaxy and that they are of EG origin and their propagation is affected by CMBR. Here we do not investigate how can they be of halo or EG origin, we just analyze their effect on the observed spectrum instead. We assume that the SP decays into two quarks (other decay modes would increase $`m_X`$ in our conclusion). After hadronization these quarks yield protons. The result is characterized by the fragmentation function (FF) $`D(x,Q^2)`$ which gives the number of produced protons with momentum fraction $`x`$ at energy scale $`Q`$. For the proton’s FF at present accelerator energies we use ref. . We evolve the FFs in ordinary  and in supersymmetric  QCD to the energies of the SPs. This result can be combined with the prediction of the MLLA technique , which gives the initial spectrum of UHECRs at the energy $`m_X`$. Altogether we study four different models: halo-SM, halo-MSSM, EG-SM and EG-MSSM.

Ref.  showed that both AGASA and Fly’s Eye data demonstrated a change of composition, a shift from heavy –iron– at $`10^{17}`$ eV to light –proton– at $`10^{19}`$ eV. Thus the UHECRs are most likely to be dominated by protons and in our analysis we use them exclusively.

The FF of the proton can be determined from present experiments . (Note, that QCD event generators e.g. HERWIG  predict the overall proton multiplicity correctly, however they describe the large $`x`$ region of the FF inaccurately.) The FFs at $`Q_0`$ energy scale are $`D_i(x,Q_0^2)`$, where $`i`$ represents the different partons (quark/squark or gluon/gluino). The FFs can not be determined in perturbative QCD; however, their evolution in $`Q^2`$ is governed by the DGLAP equations :

$$\frac{D_i(x,Q^2)}{\mathrm{ln}Q^2}=\frac{\alpha _s(Q^2)}{2\pi }\underset{j}{}_x^1\frac{dz}{z}P_{ji}(z,\alpha _s(Q^2))D_j(\frac{x}{z},Q^2),$$
(1)
One can interpret $`P_{ji}(z)`$, the splitting function, as the probability density that a parton $`i`$ produces a parton $`j`$ with momentum fraction $`z`$. Analogous evolution equations can be obtained by using coherent branching (angular ordering) for the emitted gluons . We use this technique too. Results of direct Monte-Carlo jet simulations are also available (cf. ). We include the uncertainties coming from the different choices of evolution in our final error estimates. We solve the DGLAP equations numerically with the conventional QCD (SM case) splitting functions and with the supersymmetric (MSSM case) ones . For the top and the MSSM partons we used the FFs of ref. . While solving the DGLAP equations each parton is included at its own threshold energy. Table I shows all the FFs we used.

At small values of $`x`$, multiple soft gluon emission can be described by the MLLA . This gives the shape of the total hadronic FF for soft particles (not distinguishing individual hadronic species) $`xF(x,Q^2)\mathrm{exp}[\mathrm{ln}(x/x_m)^2/(2\sigma ^2)],`$ which is peaked at $`x_m=\sqrt{\mathrm{\Lambda }/Q}`$ with $`2\sigma ^2=A\mathrm{ln}^{3/2}(Q/\mathrm{\Lambda })`$. According to  the values of $`A`$ are $`\sqrt{7/3}/6`$ and $`1/6`$ for SM and MSSM, respectively. The MLLA describes the observed hadroproduction quite accurately in the small $`x`$ region . For large values of $`x`$ the MLLA should not be used. We smoothly connect the solution for the FF obtained by the DGLAP equations and the MLLA result at a given $`x_c`$ value. Our final result on $`m_X`$ is rather insensitive to the choice of $`x_c`$, the uncertainty is included in our error estimate. We also determined the FF of the pion. Fig. 1 shows the FF for the proton and pion at $`Q=10^{16}`$ GeV in SM and MSSM.

UHECR protons produced in the halo of our galaxy can propagate practically unaffected and the production spectrum should be compared with the observations.

Particles of EG origin and energies above $`510^{19}`$ eV loose a large fraction of their energies due to interactions with CMBR . This effect can be quantitatively described by the function $`P(r,E,E_c)`$, the probability that a proton created at a distance $`r`$ with energy $`E`$ arrives at Earth above the threshold energy $`E_c`$ . This function has been calculated for a wide range of parameters in , which we use in the present calculation. The original UHECR spectrum is changed at least by two different ways: (a) there should be a steepening due to the GZK effect; (b) particles loosing their energy are accumulated just before the cutoff and produce a bump. We study the observed spectrum by assuming a uniform source distribution for UHECRs.

Our analysis includes the published and the unpublished (from the www pages of the experiments on 17/08/00) UHECR data of . Due to normalization difficulties we did not use the Yakutsk  results. We also performed the analysis using the AGASA data only and found the same value (well within the error bars) for $`m_X`$. Since the decay of SPs results in a non-negligible flux for lower energies $`\mathrm{log}(E_{min}/\text{eV})=18.5`$ is used as a lower end for the UHECR spectrum. Our results are insensitive to the definition of the upper end (the flux is extremely small there) for which we choose $`\mathrm{log}(E_{max}/\text{eV})=26`$. As it is usual we divided each logarithmic unit into ten bins. The integrated flux gives the total number of events in a bin. The uncertainties of the measured energies are about 30% which is one bin. Using a Monte-Carlo method we included this uncertainty in the final error estimates. The predicted number of events in a bin is given by

$$N(i)=_{E_i}^{E^{i+1}}\left[AE^{3.16}+Bj(E,m_X)\right],$$
(2)
where $`E_i`$ is the lower bound of the ith energy bin. The first term describes the data below $`10^{19}`$ eV according to , where the SP decay gives negligible contribution. The second one corresponds to the spectrum of the decaying SPs. A and B are normalization factors.

The expectation value for the number of events in a bin is given by eqn. (2) and it is Poisson distributed. To determine the most probable $`m_X`$ value we used the maximum-likelihood method by minimalizing the $`\chi ^2(A,B,m_X)`$ for Poisson distributed data

$$\chi ^2=\underset{i=18.5}{\overset{26.0}{}}2\left[N(i)N_o(i)+N_o(i)\mathrm{ln}\left(N_o(i)/N(i)\right)\right],$$
(3)
where $`N_o(i)`$ is the total number of observed events in the ith bin. In our fitting procedure we have three parameters: $`A,B`$ and $`m_X`$. The minimum of the $`\chi ^2(A,B,m_X)`$ function is $`\chi _{min}^2`$ at $`m_{Xmin}`$ which is the most probable value for the mass, whereas $`\chi ^2(A^{},B^{},m_X)\chi _o^2(m_X)=\chi _{min}^2+1`$ gives the one-sigma (68%) confidence interval for $`m_X`$. Here $`A^{},B^{}`$ are defined in such a way that the $`\chi ^2(A,B,m_X)`$ function is minimalized in $`A`$ and $`B`$ at fixed $`m_X`$. Fig. 2 shows the measured UHECR spectrum and the best fit in the EG-MSSM scenario. The first bump of the fit represents particles produced at high energies and accumulated just above the GZK cutoff due to their energy losses. The bump at higher energy is a remnant of $`m_X`$. In the halo models there is no GZK bump, so the relatively large $`x`$ part of the FF moves to the bump around $`5\times 10^{19}`$ GeV resulting in a much smaller $`m_X`$ than in the EG case. An interesting feature of the GZK effect is that the shape of the produced GZK bump is rather insensitive to the injected spectrum so the dependence of $`\chi ^2`$ on the choice of the FF is small. The experimental data is far more accurately described by the GZK effect (dominant feature of the EG fit) than by the FF itself (dominant for halo scenarios).

To determine the most probable value for the mass of the SP we studied 4 scenarios. Fig. 3 contains the $`\chi _{min}^2`$ values and the most probable masses with their errors for these scenarios. (The uncertainties coming from the FFs are included in our error estimates on $`m_X`$.)

The UHECR data favors the EG-MSSM scenario. The goodnesses of the fits for the halo models are far worse. The SM and MSSM cases do not differ significantly. The most important message is that the masses of the best fits (EG cases) are compatible within the error bars with the MSSM gauge coupling unification GUT scale .

The SP decay will also produce a huge number of pions which will decay into photons. Our spectrum contains 94% of pions and 6% of protons. This $`\pi /p`$ ratio is in agreement with  which showed that for different classes of models $`m_X<10^{16}`$ GeV , which is the upper boundary of our confidence intervals, the generated gamma spectrum is still consistent with the observational constraints. We performed the whole analysis including the pion produced $`\gamma `$-s in eqn. (3). The results agree with our results of Fig. 3 within errorbars, which is easy to understand. For the EG case high energy $`\gamma `$-s dominate at energies where the observed flux is zero . For the halo case the agreement is resulted by the similarity (except normalization) between $`D_p`$ and $`D_\pi `$ (cf. Fig. 1).

In the near future the UHECR statistics will probably be increased by an order of magnitude . Performing our analysis for such a statistics the uncertainty of $`m_X`$ was found to be reduced by two orders of magnitude.

Since the decay time should be at least the age of the universe it might happen that such SPs overclose the universe. Due to the large mass of the SPs a single decay results in a large number of UHECRs, thus a relatively small number of SPs can describe the observations. We checked that in all of the four scenarios the minimum density required for the best-fit spectrum is more than ten orders of magnitude smaller than the critical one.

Details will be presented in a subsequent paper .

We thank B.A. Kniehl for providing us with the proton’s FF prior to its publication and F. Csikor for useful comments. This work was partially supported by Hungarian Science Foundation grants No. OTKA-T29803/T22929-FKP-0128/1997.


format: text

text:
# Analysis of cosmic microwave background data on an incomplete sky

## 1 Introduction

Since the first measurements of the temperature anisotropy of the cosmic microwave background (CMB) by the Cosmic Background Explorer (COBE) satellite (Smoot et al. 1992), a number of sophisticated experiments have been undertaken to measure the fluctuations at higher resolutions and sensitivities (e.g. Scott et al. 1996; Tanaka et al. 1996; Netterfield et al. 1997; de Oliveira-Costa et al. 1998; Coble et al. 1999; de Bernardis et al. 2000; Wilson et al. 2000; Padin et al. 2001; Halverson et al. 2001; Lee et al. 2001; Netterfield et al. 2001). The primary result of these experiments has been the measurement of the angular power spectrum of the CMB to Legendre multipoles of up to $`l1000`$, which places strong constraints on a number of cosmological parameters (Lineweaver 1998; Efstathiou et al. 1999; de Bernardis et al. 2000; Netterfield et al. 2001; Wang, Tegmark & Zaldarriaga 2001 and references therein). In the future the Microwave Anisotropy Probe (MAP; e.g. Jarosik et al. 1998) and the Planck satellite (e.g. Bersanelli et al. 1996) will produce maps of the microwave sky with resolutions of between 5 and 30 arcmin at a number of frequencies. Such extraordinary data-sets, consisting of millions of independent measurements, will clearly require novel analysis techniques.

One of the many difficulties is the treatment of the non-cosmological contributions to the observed microwave sky. Dust, synchrotron and free-free emission from the Galaxy (e.g. Haslam et al. 1982; Schlegel, Finkbinder & Davies 1998); radio galaxies and other extra-Galactic ‘point’ sources (e.g. Toffolatti et al. 1998); and the Sunyaev-Zel’dovich (1970) effect caused by galaxy clusters (e.g. Birkinshaw 1999) all obscure the CMB at some level (see Hu, Sugiyama & Silk 1997 or Barreiro 2000 for more complete reviews), although these components have quite distinct spectral properties and so can be separated using multi-frequency observations (e.g. Bennett et al. 1992; Tegmark & Efstathiou 1996; Hobson et al. 1998; Bouchet & Gispert 1999; Jones, Hobson & Lasenby 1999; Baccigalupi et al. 2000). However these techniques are not likely to be able to extract the Galactic emissions completely (Stolyarov et al. 2001), leaving the removal of the Galactic plane as the only option. The Galaxy contributes relatively little at high latitudes (e.g. Haslam et al. 1982; Schlegel et al. 1998) so this is an acceptable, if not optimal, solution. For instance, Górski et al. (1994) removed the band within 20 deg of the Galactic plane to estimate the power spectrum of the two-year COBE Differential Microwave Radiometer sky maps, and similar cuts have been proposed by both the MAP and Planck collaborations. An essentially equivalent problem is posed if the survey’s sky coverage is incomplete, although there is less choice about the geometry of the cut in this case.

A number of aspects of the analysis become more difficult on an incomplete sphere, one of the most obvious reasons being that the spherical harmonics are no longer an orthonormal basis set. The most successful component separation techniques to date (e.g. Hobson et al. 1998; Bouchet & Gispert 1999) rely on a mode-by-mode analysis which explicitly utilises the orthogonality of the spherical harmonics, although it may be preferable to remove the Galactic plane only when estimating the CMB power spectrum. Unbiased power spectrum estimation using the spherical harmonics is possible on the cut sky (Wandelt, Górski & Hivon 2001), but the covariance structure of the resulting psuedo-harmonics is far from ideal, so analysis using an orthonormal basis set is preferable. In particular the noise covariance matrix remains diagonal in the case of spatially uniform (and uncorrelated) noise.

It is possible to construct an orthonormal basis set from linear combinations of the spherical harmonics, and an elegant implementation of this, based on Cholesky decomposition of the coupling matrix of the spherical harmonics on the cut sky, was described by Górski (1994). However the coupling matrix becomes ill-conditioned for $`l_{\mathrm{m}ax}50`$, and so this method cannot be used to perform cut-sky orthogonalisation for either the MAP experiment (with $`l_{\mathrm{m}ax}1500`$) or the Planck mission (with $`l_{\mathrm{m}ax}2500`$).

A general formalism for orthogonalisation of the spherical harmonics is presented in Section 2, although implementation to high $`l_{\mathrm{m}ax}`$ is only possible at present in the special case of a constant latitude cut. The relationship between the various harmonic coefficients is discussed in Section 3, and the extension of these results to CMB analysis techniques (specifically component separation and power spectrum estimation) is covered in Section 4. The results are summarised and future possibilities are discussed in Section 5. Finally, the chosen conventions for the spherical harmonics are defined in Appendix A; formulæ for integrals of the products of Legendre functions are given in Appendix B; and the treatment of the non-cosmological monopole and dipole modes are discussed in Appendix C.

## 2 Orthogonalisation of scalar basis
functions

The physics of the CMB is most naturally expressed in Fourier space, and it is standard practice to represent sky maps by their harmonic coefficients. The basis functions chosen here are the real spherical harmonics, $`Y_{l,m}(\widehat{𝒓})`$ (as defined in Appendix A), which form an orthonormal basis on the complete sphere, $`S`$. In general $`l0`$ and $`lml`$, although in practice a finite $`l_{\mathrm{m}ax}`$ must be used, which implies a band-limit. It is convenient to combine the two indices, allowing the basis set to be expressed as a vector, $`𝒀(\widehat{𝒓})=[Y_1(\widehat{𝒓}),Y_2(\widehat{𝒓}),\mathrm{},Y_{i_{\mathrm{m}ax}}(\widehat{𝒓})]^\mathrm{T}`$, where $`i_{\mathrm{m}ax}=(l_{\mathrm{m}ax}+1)^2`$. There are several reasonable choices for the indexing, $`i(l,m)`$, most notably grouping coefficients in $`l`$ or $`m`$, as defined in Appendix A. Grouping in $`l`$ is most natural for power spectrum estimation, but grouping in $`m`$ is more efficient computationally in cases of azimuthal symmetry (Section 2.3).

The spherical harmonics are not orthogonal on the incomplete sphere $`S^{}`$, as can be seen from the structure of their coupling matrix (Section 2.1). A decomposition of the coupling matrix can be used to construct an orthonormal basis set (Section 2.2), but implementation to high-resolution is currently possible only in the special case of constant latitude cuts (Section 2.3).

### 2.1 The coupling matrix

The coupling matrix of a set of functions encodes their orthogonality and normalisation properties over a given range. In the case of the spherical harmonics on the incomplete sphere it is given by

$$\text{C}=_S^{}𝒀(\widehat{𝒓})𝒀^\mathrm{T}(\widehat{𝒓})d\mathrm{\Omega }.$$
(1)
If $`S^{}=S`$ then the harmonics are orthonormal and $`\text{C}=\text{I}`$; otherwise the off-diagonal elements are non-zero, indicating that the basis functions are non-orthogonal.

An alternative formulation is to introduce a window function, $`w(\widehat{𝒓})`$, so that

$$\text{C}=_Sw^2(\widehat{𝒓})𝒀(\widehat{𝒓})𝒀^\mathrm{T}(\widehat{𝒓})d\mathrm{\Omega }.$$
(2)
In some ways this approach is more flexible, as $`w(\widehat{𝒓})`$ can either be a smoothly-varying apodizing function (cf. Tegmark 1997) or take the form

$$w_S^{}(\widehat{𝒓})=\{\begin{array}{ccc}1,\hfill & \mathrm{i}f\hfill & \widehat{𝒓}S^{},\hfill \\ & & \\ 0,\hfill & \mathrm{i}f\hfill & \widehat{𝒓}SS^{},\hfill \end{array}$$
(3)
mimicing the effect of the sharp cut defined above. However this definition of the window function can lead to inconsistencies if a band-limited analysis is carried out, as $`w_S^{}(\widehat{𝒓})`$ cannot be properly represented by a finite analysis (see Section 3.2). It is for this reason that the first formalism is used here, although most of the subsequent results can also be derived using window functions.

For a pixel-based analysis, the coupling matrix can be defined by replacing the integral in equation (1) by a sum over points on the sphere (i.e. pixel centres), $`\widehat{𝒓}_p`$, where $`1pN_\mathrm{p}`$ and $`N_\mathrm{p}`$ is the number of pixels. In this case

$$\text{C}=\underset{p=1;\widehat{𝒓}_pS^{}}{\overset{N_\mathrm{p}}{}}𝒀(\widehat{𝒓}_p)𝒀^\mathrm{T}(\widehat{𝒓}_p)\mathrm{\Omega }_p,$$
(4)
where $`\mathrm{\Omega }_p`$ is the area of the $`p`$th pixel and there is no pixel-smoothing (cf. Górski 1994). In the limit $`N_\mathrm{p}\mathrm{}`$ equations (1) and (4) become equivalent and pixelisation issues become irrelevant. If the points are uniformly distributed over the sphere C should be close to the identity, the small discrepancies merely reflecting the approximation of the integral as a sum; otherwise C reflects the spatial distribution of the points much as in the continuum case, but there is freedom to represent apodizing filters as well as discrete cuts.

Assuming $`\mathrm{\Omega }_S^{}>0`$, the coupling matrix is formally symmetric, positive definite and invertible, irrespective of which of the above definitions is used. However C rapidly becomes numerically singular: e.g. if $`l_{\mathrm{m}ax}=50`$, the condition number111The condition number of a matrix is the (absolute value of) the ratio of its greatest and smallest eigenvalues; it is large for ill-conditioned matrices, and infinite for singular matrices. of C is $`5\times 10^9`$ for a constant latitude cut of $`b_{\mathrm{c}ut}=\pm 20`$ deg. This can be further understood in terms of the eigenstructure of the coupling matrix.

#### 2.1.1 Eigenstructure

The coupling matrix has $`i_{\mathrm{m}ax}=(l_{\mathrm{m}ax}+1)^2`$ eigenvectors, $`𝒗_i`$, and eigenvalues, $`\lambda _i`$, which satisfy

$$\text{C}𝒗_i=\lambda _i𝒗_i.$$
(5)
Premultiplying by $`𝒀^\mathrm{T}(\widehat{𝒓})`$ and expanding out the implicit summations gives

$$_S^{}\underset{k=1}{\overset{i_{\mathrm{m}ax}}{}}Y_k(\widehat{𝒓}^{})(v_i)_k\underset{j=1}{\overset{i_{\mathrm{m}ax}}{}}Y_j(\widehat{𝒓})Y_j(\widehat{𝒓}^{})\mathrm{d}\mathrm{\Omega }^{}=\lambda _iv_i(\widehat{𝒓}),$$
(6)
where $`v_i(\widehat{𝒓})=𝒀^\mathrm{T}(\widehat{𝒓})𝒗_i`$ is the $`i`$th eigenfunction of the coupling matrix. The completeness of the spherical harmonics in the $`l_{\mathrm{m}ax}\mathrm{}`$ limit implies that $`𝒀^\mathrm{T}(\widehat{𝒓})𝒀(\widehat{𝒓}^{})=\delta (\widehat{𝒓}\widehat{𝒓}^{})`$ \[where $`\delta (x)`$ is the Dirac delta function\], so that equation (6) reduces to

$$_S^{}v_i(\widehat{𝒓}^{})\delta (\widehat{𝒓}\widehat{𝒓}^{})d\mathrm{\Omega }^{}=\lambda _iv_i(\widehat{𝒓}).$$
(7)
For a given $`i`$ this must be true at all $`\widehat{𝒓}`$, which implies that either: $`v_i(\widehat{𝒓})=0`$ in the cut region, $`SS^{}`$, in which case $`\lambda _i=1`$; or $`v_i(\widehat{𝒓})=0`$ in $`S^{}`$, in which case $`\lambda _i=0`$. In other words these eigenfunctions are completely localised in either the cut sphere or the removed region. This bimodality is only strictly true in the $`l_{\mathrm{m}ax}\mathrm{}`$ limit, but, as shown in Fig. 1, is a good approximation for $`l_{\mathrm{m}ax}500`$.

As the coupling matrix is symmetric, those eigenvectors with different eigenvalues are orthogonal, and those with the same eigenvalues can be made orthogonal by a rotation in the subspace defined by the eigenvalue in question (e.g. Arfken 1985). Thus the eigenfunctions with $`\lambda _i=1`$ represent an orthogonal basis set on $`S^{}`$, whereas those with $`\lambda _i=0`$ have no support in this region and so cannot be orthogonal (or normalised) on the cut sky. The freedom in choice of basis does not extend to mixing the $`\lambda _i=0`$ modes (i.e. those corresponding to the null-space of C) with the $`\lambda _i=1`$ modes (i.e. those in the range of C), and so the number of supported modes is determined by a combination of the band-limit and the cut.

The number of orthonormal basis functions (i.e. the rank of C) is proportional to the area of the sphere retained, $`\mathrm{\Omega }_S^{}`$. Hence it is possible to define only

$$i_{\mathrm{m}ax}^{}\frac{\mathrm{\Omega }_S^{}}{\mathrm{\Omega }_S}i_{\mathrm{m}ax}=\frac{\mathrm{\Omega }_S^{}}{4\pi }i_{\mathrm{m}ax}$$
(8)
orthonormal functions on the cut sphere for a given (large) band-limit. The relative reduction in the basis set is the same as would occur in the equivalent pixel analysis: the number of pixels retained is also given by $`N_\mathrm{p}^{}\mathrm{\Omega }_S^{}/\mathrm{\Omega }_SN_\mathrm{p}`$. For low $`l_{\mathrm{m}ax}`$ these arguments do not hold, and it is possible to create a basis set with more than $`\mathrm{\Omega }_S^{}/\mathrm{\Omega }_Si_{\mathrm{m}ax}`$ elements. Moreover, all these functions are required to ensure that the cut-sphere basis set is complete (as well as orthonormal) in the case of a low band-limit.

### 2.2 Construction of an orthonormal basis

The construction of an orthonormal basis set from a set of linearly-independent functions is a well-established mathematical technique, and a number of orthogonalisation methods are possible. The most basic is Gram-Schmit orthogonalisation (e.g. Arfken 1985), in which the new basis functions are built-up sequentially, but this algorithm is numerically unstable. The modified Gram-Schmidt algorithm (e.g. Golub & van Loan 1996) is stable, but it is generally preferable to use matrix techniques to create all the new basis functions simultaneously.

Starting with the spherical harmonics, $`𝒀(\widehat{𝒓})`$, the task is to find a set of functions $`𝒀^{}(\widehat{𝒓})`$ which are orthonormal on the incomplete sphere. In terms of a conversion matrix, B, the two sets of functions are related by

$$𝒀^{}(\widehat{𝒓})=\text{B}𝒀(\widehat{𝒓}).$$
(9)
Note that B has dimensions $`i_{\mathrm{m}ax}^{}\times i_{\mathrm{m}ax}`$, where $`i_{\mathrm{m}ax}=(l_{\mathrm{m}ax}+1)^2`$ and $`i_{\mathrm{m}ax}^{}i_{\mathrm{m}ax}`$ is determined by the band-limit and the cut, as described in Section 2.1. It is also important to note that the indexing of the $`Y_i^{}(\widehat{𝒓})`$ is qualitatively different from the $`Y_i(\widehat{𝒓})`$. The latter are really two-index objects, with their characteristic scale given by $`\pi /l`$ and $`m`$ relating to ‘orientation’. However the new basis functions include contributions from spherical harmonics with different $`l`$-values, and thus do not have a well-defined angular scale. Hence their single index contains no physical information, and the ordering or grouping of the new basis functions is arbitrary.

From equation (1), the coupling matrix of these new functions is

$`\text{C}^{}`$ $`=`$ $`{\displaystyle _S^{}}𝒀^{}(\widehat{𝒓})𝒀_{}^{}{}_{}{}^{\mathrm{T}}(\widehat{𝒓})d\mathrm{\Omega }`$ (10)
$`=`$ $`\text{B}\text{C}\text{B}^\mathrm{T}.`$
Hence any conversion matrix which satisfies222Here $`\stackrel{~}{\text{I}}`$ is the $`i_{\mathrm{m}ax}^{}\times i_{\mathrm{m}ax}^{}`$ identity matrix, as distinct from the (potentially larger) $`i_{\mathrm{m}ax}\times i_{\mathrm{m}ax}`$ identity matrix, I.

$$\text{B}\text{C}\text{B}^\mathrm{T}=\stackrel{~}{\text{I}}$$
(11)
yields basis functions which are orthonormal on the cut sphere, and the task of orthogonalisation is reduced to finding a solution for B given C. Whilst such a solution does not exist for arbitrary C, in all cases of practical interest a suitable conversion matrix can be constructed from the coupling matrix. One possible method is direct calculation of the eigenstructure of C, which yields a conversion matrix with elements given by $`B_{i^{},i}=(v_i^{})_i\lambda _i^{}^{1/2}`$, where the $`𝒗_i^{}`$ are the eigenvectors of C and the $`\lambda _i^{}`$ its positive eigenvalues. However it is advantageous to include the symmetry of the coupling matrix explicitly, which leads to a factorization of the form

$$\text{C}=\text{A}\text{A}^\mathrm{T},$$
(12)
where A is an $`i_{\mathrm{m}ax}\times i_{\mathrm{m}ax}^{}`$ matrix, the form of which is determined by the decomposition method. Combining equations (11) and (12), the task of orthogonalisation is reduced to finding B such that

$$\text{B}\text{A}=\stackrel{~}{\text{O}},$$
(13)
where $`\stackrel{~}{\text{O}}`$ is an $`i_{\mathrm{m}ax}^{}\times i_{\mathrm{m}ax}^{}`$ orthogonal matrix (i.e. $`\stackrel{~}{\text{O}}\stackrel{~}{\text{O}}^\mathrm{T}=\stackrel{~}{\text{O}}^\mathrm{T}\stackrel{~}{\text{O}}=\stackrel{~}{\text{I}}`$).

Whilst equations (12) and (13) are general expressions which must be satisfied by the conversion matrix, they do not define a definite algorithm for the orthogonalisation. In practice it is simplest to choose $`\stackrel{~}{\text{O}}=\stackrel{~}{\text{I}}`$, leading to the requirement that

$$\text{B}\text{A}=\stackrel{~}{\text{I}}.$$
(14)
However the optimal choice of decomposition method used to generate A depends on whether the coupling matrix is (numerically) invertible, and hence on the band-limit of the analysis.

#### 2.2.1 Low-resolution analysis

If $`l_{\mathrm{m}ax}50`$ and most of the sphere is retained (i.e. $`\mathrm{\Omega }_S^{}\mathrm{\Omega }_S/2=2\pi `$) then the coupling matrix is numerically invertible and can be treated as positive definite in practice. Consequently A \[defined in equation (12)\] is also invertible, and $`\text{B}=\text{A}^1`$, so that, from equation (9), the orthonormal basis set is given by

$$𝒀^{}(\widehat{𝒓})=\text{A}^1𝒀(\widehat{𝒓}).$$
(15)
The form of A depends on the factorization method; of the wide variety available (e.g. Golub & van Loan 1996), the two most useful here are singular value decomposition (SVD) and Cholesky decomposition.

The SVD of the covariance matrix is defined in terms of equation (12) by $`\text{A}=\text{V}\text{W}^{1/2}`$ (i.e. $`\text{C}=\text{V}\text{W}\text{V}^\mathrm{T}`$), where V is orthogonal and W is diagonal333If M is diagonal then the notation $`\text{M}^{\pm 1/2}`$ is used here to denote the matrix defined by $`(M^{\pm 1/2})_{i,j}=\delta _{i,j}M_{i,i}^{\pm 1/2}`$, where $`\delta _{i,j}`$ is the Kronecker delta function. Thus $`\text{M}^{1/2}`$ only exists if the diagonal elements of M are non-negative and $`\text{M}^{1/2}`$ only exists if the diagonal elements of M are strictly positive.. The diagonal elements of W are the eigenvalues of C and, as their ordering is arbitrary, W can be defined such that $`W_{i,i}W_{i+1,i+1}`$, provided the columns of V are permuted in the same way. The columns of V, in turn, are the eigenvectors of C, and V is an orthogonal matrix (i.e. $`\text{V}^1=\text{V}^\mathrm{T}`$). Hence the conversion matrix is given by $`\text{B}=\text{A}^1=\text{W}^{1/2}\text{V}^\mathrm{T}`$, which is trivially computed once the SVD has been performed. Note that this approach is effectively the same as the direct calculation of the eigenstructure of C mentioned above in Section 2.2.

Whilst SVD is a powerful technique, it is computationally expensive – a Cholesky decomposition is approximately ten times faster, although it can only be performed on symmetric matrices which are numerically positive definite. The Cholesky decomposition of the covariance matrix takes the form $`\text{C}=\text{L}\text{L}^\mathrm{T}`$ (i.e. $`\text{A}=\text{L}`$), where L is lower triangular. Hence the conversion matrix, $`\text{B}=\text{L}^1`$ can be computed quickly from the initial factorization in this case as well. The triangular structure of the conversion matrix also ensures that the new basis functions are the same as those formed by a numerically-stable Gram-Schmidt orthogonalisation (Górski 1994).

Despite the fact that the SVD and the Cholesky decomposition result in quite different sets of basis functions, there is no reason to prefer one over the other in general. In the case of CMB analysis, however, the triangular structure of A and B as generated by the Cholesky decomposition is preferable as it ensures that the non-cosmological monopole and dipole modes are kept separate from the $`l2`$ modes, assuming $`l`$-ordering is used (Górski 1994). If the SVD route is taken (or another indexing scheme used) the separation of the $`l=0`$ and $`l=1`$ modes can be ensured using the partial Householder transform described in Appendix C. Nonetheless, if the coupling matrix is sufficiently non-singular, a Cholesky decomposition should be used to create the orthonormal basis set, due to both its computational efficiency and the simplicity with which the non-cosmological modes are handled.

#### 2.2.2 High-resolution analysis

If $`l_{\mathrm{m}ax}50`$ the coupling matrix is numerically singular, and thus A \[defined in equation (12)\] is non-invertible. Cholesky decomposition of C is thus impractical and, whilst an SVD is possible, the conversion matrix as defined in Section 2.2.1 cannot be computed, as the smallest elements of W (i.e. the smallest eigenvalues of C) are so close to zero. This implies that the corresponding columns of V do not contribute to the reconstruction of C and can be ignored. Hence it is possible to perform an approximate SVD of the coupling matrix, defined by $`\text{C}\stackrel{~}{\text{V}}\stackrel{~}{\text{W}}\stackrel{~}{\text{V}}^\mathrm{T}`$ (i.e. $`\text{A}=\stackrel{~}{\text{V}}\stackrel{~}{\text{W}}^{1/2}`$), where $`\stackrel{~}{\text{W}}`$ is an $`i_{\mathrm{m}ax}^{}\times i_{\mathrm{m}ax}^{}`$ diagonal matrix containing the largest elements of W and $`\stackrel{~}{\text{V}}`$ is an $`i_{\mathrm{m}ax}^{}\times i_{\mathrm{m}ax}`$ matrix consisting of the corresponding columns of V. The value of $`i_{\mathrm{m}ax}`$ is determined by the choice of $`W_{\mathrm{m}in}`$ used to truncate W, but the bimodality of the eigenvalue distribution means that any value between $`10^5`$ and $`0.1`$ is acceptable. The resultant conversion matrix is $`\text{B}=\stackrel{~}{\text{W}}^{1/2}\stackrel{~}{\text{V}}^\mathrm{T}`$ (satisfying $`\text{B}\text{A}\stackrel{~}{\text{I}}`$) and the $`i_{\mathrm{m}ax}^{}`$ new basis functions are given by

$$𝒀^{}(\widehat{𝒓})=\stackrel{~}{\text{W}}^{1/2}\stackrel{~}{\text{V}}^\mathrm{T}𝒀(\widehat{𝒓}).$$
(16)
These basis functions represent an orthonormal basis set on the incomplete sphere, but they are not formally complete to the nominal band-limit due to the slightly approximate nature of the reduced SVD. The decomposition becomes exact in the limit $`l_{\mathrm{m}ax}\mathrm{}`$ as $`\stackrel{~}{\text{W}}\stackrel{~}{\text{I}}`$ (from the eigenvalue arguments described in Section 2.1.1) and the reduced SVD becomes $`\text{C}\stackrel{~}{\text{V}}\stackrel{~}{\text{V}}^\mathrm{T}`$ (i.e. $`\text{A}\stackrel{~}{\text{V}}`$ and $`\text{B}\stackrel{~}{\text{V}}^\mathrm{T}`$) in this limit. Note also that these basis functions are still orthogonal on the full sphere, although they are no longer normalised to unity; this is another potential advantage of the SVD-based method.

Several examples of orthonormal cut-sphere basis functions are shown in Fig. 2, for both symmetric and asymmetric constant latitude cuts. The link between these functions and the spherical harmonics is apparent – they have the same cellular structure, but, for the most part, are combined in such a way that their support is localised to $`S^{}`$. However the functions shown in Fig. 2 (g) and (h) are from the small subset with intermediate values of $`W_{i,i}`$ and, as such, have considerable support in the removed region.

In the case of CMB analysis, one short-coming of this approach is that the non-cosmological monopole and dipole modes are not distinguished from the higher moments (see Section 3.2), although this separation can be achieved post facto by using a partial Householder transform, as described in Appendix C. In doing this some of the useful properties of the SVD are lost, but this operation need only be performed as the last step in generating the orthonormal basis set, by which stage all the computationally intensive matrix operations have already been performed.

Despite the inconvenience caused by the mixing of the monopole and dipole modes, SVD is clearly the most flexible and general method of orthogonalisation. In part, this stems from the fact that it can be applied to the coupling matrix without any prior knowledge of its singular properties. However, in the high-$`l_{\mathrm{m}ax}`$ limit, the number of supported basis functions is determined by the area of the cut and given, to a good approximation, by equation (8). Hence faster, if less powerful, techniques can be used to orthogonalise the spherical harmonics for $`l_{\mathrm{m}ax}1000`$ as the coupling matrix is guaranteed to have at least $`(l_{\mathrm{m}ax}+1)^2\mathrm{\Omega }_S^{}/(4\pi )`$ positive eigenvalues. An example of this idea would be to modify the pivoting algorithm in the psuedo-Cholesky decomposition described in Section 4.2.9 of Golub & van Loan (1996) so that the decomposition halts when the predetermined number of basis functions have been generated, as opposed to using the less robust threshold based on the values of the diagonal elements of C.

### 2.3 Constant latitude cuts

In principle the method presented above is a complete solution to the problem of constructing orthonormal bases on the cut sphere, but the coupling matrix requires $`O(l_{\mathrm{m}ax}^4)`$ storage, limiting a general implementation to $`l_{\mathrm{m}ax}200`$ on most current computers. Furthermore, the SVD of an $`n\times n`$ matrix requires $`O(n^3)`$ operations, and so the orthogonalisation operation count scales as $`O(l_{\mathrm{m}ax}^6)`$. Similar difficulties are encountered in merely evaluating C, regardless of whether numerical integration or recursive techniques are used (e.g. Hivon et al. 2001).

Fortunately all these difficulties are significantly reduced in the case of a constant latitude cut (cf. Oh, Spergel & Hinshaw 1999; Wandelt et al. 2001), defined by ignoring all $`\theta `$ for which $`\theta _1\theta \theta _2`$. This could be the symmetric removal of the Galactic plane (i.e. $`\theta _1=\pi /2b_{\mathrm{c}ut}`$ and $`\theta _2=\pi /2+b_{\mathrm{c}ut}`$, where $`b_{\mathrm{c}ut}`$ is the latitude of the cut) or the absence of data round one pole (i.e. $`\theta _1=0`$ and $`\theta _2=\theta _{\mathrm{c}ut}`$). The formalism derived below can also be trivially extended to include multiple cuts, as would be required for a CMB experiment which did not observe either ecliptic pole.

Explicitly including the constant latitude cut in equation (1), the elements of the coupling matrix are given by

$`C_{i(l,m),i(l^{},m^{})}`$ $`=`$ $`{\displaystyle _0^{2\pi }}s_m(\varphi )s_m^{}(\varphi )d\varphi `$
$`\times `$ $`[{\displaystyle _1^{\mathrm{cos}(\theta _2)}}\lambda _{l,|m|}(x)\lambda _{l^{},|m^{}|}(x)\mathrm{d}x`$
$`+{\displaystyle _{\mathrm{cos}(\theta _1)}^1}\lambda _{l,|m|}(x)\lambda _{l^{},|m^{}|}(x)\mathrm{d}x],`$
where $`s_m(\varphi )`$ is defined in equation (44), and the $`\lambda _{l,m}(x)`$ are normalised associated Legendre functions, given in equation (45). From Appendix A, the first integral in equation (2.3) reduces to $`2\pi \delta _{m,m^{}}`$, and so

$`C_{i(l,m),i(l^{},m^{})}`$ $`=`$ $`\delta _{m,m^{}}`$
$`\times `$ $`[\delta _{l,l^{}}2\pi {\displaystyle _{\mathrm{cos}(\theta _2)}^{\mathrm{cos}(\theta _1)}}\lambda _{l,|m|}(x)\lambda _{l^{},|m|}(x)dx].`$
The remaining integral can be evaluated using a combination of analytical formulæ and recursion relations, as described in Appendix B.

The most important aspect of equation (2.3) is that the coupling matrix C is extremely sparse (only one element in $`l_{\mathrm{m}ax}`$ is non-zero) and, if stored using the indexing scheme defined in equation (54) (i.e. grouped into sub-matrices of fixed $`m`$), is block diagonal. C can thus be stored in the form of $`2l_{\mathrm{m}ax}+1`$ sub-matrices, the $`m`$th of which has $`(l_{\mathrm{m}ax}+1|m|)^2`$ elements, and the storage requirements thus scale as $`O(l_{\mathrm{m}ax}^3)`$ rather than $`O(l_{\mathrm{m}ax}^4)`$. Whilst it is convenient to store all the blocks simultaneously, there is no need to do so, which can further reduce the storage requirements to $`O(l_{\mathrm{m}ax}^2)`$. It is also clear from equation (2.3) that only the $`m0`$ terms need be treated explicitly and that $`l`$ and $`l^{}`$ are interchangeable, decreasing the storage requirements by a further factor of four. Finally, in the case of a symmetric cut (i.e. $`\theta _2=\pi \theta _1`$) the parity of $`\lambda _{l,m}(x)`$ is such that all terms for which $`l+l^{}`$ is odd vanish, resulting in an additional halving of the memory requirements.

The orthogonalisation can be performed by decomposing each sub-matrix separately, reducing the operation count from $`O(l_{\mathrm{m}ax}^6)`$ to $`O(l_{\mathrm{m}ax}^4)`$. The removal of the poorly-supported basis functions is achieved in the same manner as described in Section 2.2, although the book-keeping is more complicated. Similarly the partial Householder transform required to separate the $`l=0`$ and $`l=1`$ modes need only be applied to the $`m=0`$ and $`m=\pm 1`$ blocks of the resultant conversion matrix (Appendix C). An important side-effect of the separation in $`m`$ is that the $`Y_i^{}(\widehat{𝒓})`$ have the same trigonometric $`\varphi `$-dependence as the full-sky spherical harmonics (Appendix A). This also implies that the $`Y_i^{}(\widehat{𝒓})`$ can be treated as two-index quantities, defined by $`m`$ and a second, arbitrary index in place of $`l`$.

The algorithms described here were implemented on the Cambridge Centre for Mathematical Science’s COSMOS 64-processor Silicon Graphics Origin 2000 and the evaluation and decomposition of the coupling matrix at the highest Planck resolution of $`l_{\mathrm{m}ax}`$ 2500 required about an hour. The majority of the time was spent factorizing the sub-matrices, and thus significant accelerations are unlikely, the highly-optimised linear algebra package (lapack; Anderson 1992) routines having been used for all the decompositions. For a given choice of $`\theta _1`$ and $`\theta _2`$, the decomposition of C need only be performed once, so orthogonalisation of the spherical harmonics on an incomplete sky should comprise only a small fraction of the analysis required for the forthcoming MAP and Planck missions.

## 3 Harmonic analysis

Methods for constructing an orthonormal basis set on the incomplete sphere from the spherical harmonics were discussed in Section 2, but in most cases of data analysis it is the harmonic coefficients, representing functions on the sphere, that are of interest. There are at least three useful harmonic expansions of a general function on the sphere, and the relationships between these coefficients, which are summarised in Table LABEL:table:harm, are derived here.

A band-limited function, $`a(\widehat{𝒓})`$, can be completely specified by a finite number of harmonic coefficients as (cf. Appendix A)

$$a(\widehat{𝒓})=𝒀^\mathrm{T}(\widehat{𝒓})𝒂,$$
(19)
where it is assumed that $`l_{\mathrm{m}ax}`$ is greater than or equal to the band-limit of $`a(\widehat{𝒓})`$ and the harmonic coefficients are defined by

$$𝒂=_S𝒀(\widehat{𝒓})a(\widehat{𝒓})d\mathrm{\Omega }.$$
(20)
The invertibility of these transformations is due to the orthonormality of the spherical harmonics on $`S`$ and the fact that they represent a complete basis set given the band-limit.

If $`a(\widehat{𝒓})`$ is only known over some fraction of the sphere $`S^{}S`$, then $`𝒂`$ cannot be determined as above, as the integral in equation (20) is incomplete. In this case the psuedo-harmonics

$$\stackrel{~}{𝒂}=_S^{}𝒀(\widehat{𝒓})a(\widehat{𝒓})d\mathrm{\Omega }$$
(21)
fully specify $`a(\widehat{𝒓})`$ in $`S^{}`$ due to the band-limit. From equation (19) they are related to the full harmonic coefficients by

$$\stackrel{~}{𝒂}=\text{C}𝒂,$$
(22)
where C is the coupling matrix, defined in equation (1).

The psuedo-harmonics are useful quantities, but it is preferable to work with basis functions that are orthonormal on $`S^{}`$. Denoted $`𝒀^{}(\widehat{𝒓})`$ in Section 2, their harmonic coefficients are given by

$$𝒂^{}=_S^{}𝒀^{}(\widehat{𝒓})a(\widehat{𝒓})d\mathrm{\Omega }.$$
(23)
The relationship between $`𝒀(\widehat{𝒓})`$ and $`𝒀^{}(\widehat{𝒓})`$ given in equation (9) flows through to the harmonic coefficients and applying equations (9), (12) and (19) to equation (23) gives

$`𝒂^{}`$ $`=`$ $`\text{B}{\displaystyle _S^{}}𝒀(\widehat{𝒓})\left[𝒀^\mathrm{T}(\widehat{𝒓})𝒂\right]d\mathrm{\Omega }`$ (24)
$`=`$ BC$`𝒂`$
$`=`$ $`\text{A}^\mathrm{T}𝒂.`$
The form of the above transformation depends on the decomposition used to generate A (cf. Section 2.2), but, for CMB analysis, it is desirable to separate the non-cosmological modes. This amounts to demanding that the only four of the $`a_i^{}`$ have any contribution from the $`l=0`$ and $`l=1`$ spherical harmonic coefficients. The uppper triangular structure of $`\text{A}^\mathrm{T}=\text{L}`$ as generated by a Cholesky decomposition inherently satisfies this requirement, but in general the conversion matrix needs to be transformed explicitly. One option is to use successive partial Householder transforms, as described in detail in Appendix C. As can be seen from Fig. 3, the index-ordering and decomposition method combine to give a wide variety of conversion matrices; which of these is most suitable depends on the application.

It is also possible to convert between the psuedo-harmonics and the cut-sky harmonics, as they both contain information about $`a(\widehat{𝒓})`$ in $`S^{}`$ alone. Combining equations (12), (21) and (23) implies that $`\stackrel{~}{𝒂}=\text{A}𝒂^{}`$ and $`𝒂^{}=\text{B}\stackrel{~}{𝒂}`$. However it is not always possible to determine $`𝒂`$ from either $`𝒂^{}`$ or $`\stackrel{~}{𝒂}`$. In the low-$`l_{\mathrm{m}ax}`$ limit these inversions are defined (Section 3.1), but for appreciable band-limits only a projection onto the cut sphere is possible (Section 3.2).

### 3.1 Low-resolution analysis

If the coupling matrix is numerically non-singular (i.e. $`l_{\mathrm{m}ax}50`$) then equation (24) can be inverted to give

$$𝒂=\left(\text{A}^\mathrm{T}\right)^1𝒂^{}=\text{B}^\mathrm{T}𝒂^{},$$
(25)
and equation (21) implies that

$$𝒂=\text{C}^1\stackrel{~}{𝒂}.$$
(26)
These are specific examples of the fact that a band-limited function is completely defined if it is known over any finite portion of the sphere, and a cut-sky analysis serves no purpose – any apparently localised contaminants infect the entire sky. However any measurement of a field on the sphere is subject to noise which is not band-limited, in which case the application of a cut has the effect of greatly amplifying the noise in the removed region (see Section 4.4), justifying the use of a cut-sky analysis in the low-resolution limit.

### 3.2 High-resolution analysis

If $`l_{\mathrm{m}ax}50`$ and the coupling matrix is numerically singular, it is impossible to reconstruct (even) a band-limited function that is known only on $`S^{}`$. The loss of information about modes constrained to the cut makes it clear that the analysis has the desired effect of removing contaminated (or otherwise problematic) regions, but the most appropriate transformation from the cut-sphere basis to conventional harmonics is less obvious.

A least squares-like approach leads to a definition of the reconstructed full sky coefficients as

$$\widehat{𝒂}\text{B}^\mathrm{T}𝒂^{}=\left(\text{A}\text{B}\right)^\mathrm{T}𝒂.$$
(27)
Similarly, equation (22) implies that

$$\widehat{𝒂}\text{B}^\mathrm{T}\text{B}\stackrel{~}{𝒂}=\left(\text{A}\text{B}\right)^\mathrm{T}𝒂,$$
(28)
where $`(\text{A}\text{B})^\mathrm{T}`$ is a projection operator444The definition of B given in Section 2.2 implies that $`[(\text{A}\text{B})^\mathrm{T}]^2=(\text{A}\text{B})^\mathrm{T}`$ and it is hence a projection operator if $`l_{\mathrm{m}ax}\mathrm{}`$. onto the range of C, which in real space is $`S^{}`$. If $`l_{\mathrm{m}ax}\mathrm{}`$ it is possible to write $`\widehat{a}(\widehat{𝒓})=w_S^{}(\widehat{𝒓})a(\widehat{𝒓})`$, where $`w_S^{}(\widehat{𝒓})`$ is the sharp window function defined in equation (3).

It is at this point that the subtle distinctions between the use of a discrete cut and a window function become apparent. These results only hold for band-limited functions, but, as defined above, $`\widehat{a}(\widehat{𝒓})`$ is not band-limited, and so cannot be analysed self-consistently. Whether a discrete cut or an apodizing function is to be preferred depends on the situation in which the incomplete sky analysis is required.

## 4 CMB data analysis

In order to determine the properties of the CMB from noisy observations of the microwave sky (Section 4.1) a number of non-trivial analysis steps are required, including: map-making (Section 4.2); component separation (Section 4.3) and power spectrum estimation (Section 4.4). Several algorithms have been suggested for all these steps, and these are discussed briefly below, but the main focus is on when and how to apply a sky cut. Further, whilst it is possible to analyse the data in either real or Fourier space, the latter approach is emphasised here as it is more directly related to the formalism described in Sections 2 and 3, as well as being the focus of a related series of papers (van Leeuwen et al. 2001; Challinor et al. 2001; Stolyarov et al. 2001). Note that the term ‘map’ is used here to denote any representation of a field on the sky and can imply either a set of real space pixel values or a vector of spherical harmonic coefficients.

### 4.1 Observations

Observations of the CMB can be made using a number of quite distinct techniques. Data have been obtained from the ground, high altitude balloons, and satellites, but the more important distinction is the type of telescope. The experiments listed in Section 1 include: straightforward single dish telescopes, such as BOOMERanG (e.g. Netterfield et al. 2001) and Planck (e.g. Bersanelli et al. 1996); differencing experiments, such as COBE (e.g. Smoot et al. 1992) and MAP (e.g. Jarosik et al. 1998); and interferometers, such as the Cambridge Anisotropy Telescope (CAT; Scott et al. 1996), the Cosmic Backround Imager (CBI; Padin et al. 2001) and the Degree Angular Scale Interferometer (DASI; Halverson et al. 2001). The interferometry surveys inevitably cover only a small fraction of the sky, and so a flat-sky Fourier analysis becomes possible. However the both the differencing and single dish surveys can, in principle, cover most of the celestial sphere, and should yield maps of the microwave sky that are limited only by (the combined effects of) instrument noise and the finite telescope beam.

#### 4.1.1 Noise

The typical receivers used in the above experiments have two main noise contributions: random white noise and a correlated low-frequency (i.e. ‘$`1/f`$’) component. The latter is potentially troublesome, leading to ‘stripy’ maps with correlated errors, and is the main reason for the popularity of differencing experiments which remove low frequency noise at the moment of observation. However data from single dish surveys can be ‘de-striped’ if the scan strategy includes sufficiently many multiply-observed points (e.g. Tegmark 1997; Delabrouille 1998; Maino et al. 1999) or the time-time noise covariance matrix can be fully included in the map-making process (Wright, Hinshaw & Bennett 1996; Natoli et al. 2001; Challinor et al. 2001). Hence correlated errors are ignored in the simple analysis presented here.

This leaves only the white component, which can be analysed most simply in the case of a single beam experiment. Following Knox (1995), a receiver is characterised by its sensitivity, $`s`$ (generally chosen to have units of temperature time1/2). Assuming the noise is Gaussian it has expectation values $`n=0`$ and $`n^2=s^2t`$ over an integration time $`t`$. The manner in which this noise projects onto a sky map depends on the map-making algorithm, the scan strategy, and the beam.

#### 4.1.2 Beam convolution

All telescopes necessarily have a finite point-spread function or beam, which, for a given detector, can be characterised by $`b(\widehat{𝒓})`$, the fraction of photons from direction $`\widehat{𝒓}`$ that are registered, given a nominal orientation towards the north pole (i.e. $`\theta =0`$). The harmonic expansion of the beam in this orientation is denoted $`𝒃=b_{i(l,m)}`$, with the band-limit being related to the nominal resolution of the detector. For a given type of telescope the resolution improves with frequency due to diffraction effects; this places limitations on the component separation algorithms that are used on the incomplete sky (Section 4.3).

Most experiments have beams that are manifestly asymmetric, a fact which must be accounted for explicitly by the data analysis algorithms, but the cut-sky issues of interest here can be explored more clearly if the beam is approximated by its azimuthally averaged counterpart (e.g. Challinor et al. 2001). Defined by

$$\overline{b}(\theta )=\frac{1}{2\pi }_0^{2\pi }b(\theta ,\varphi )d\varphi ,$$
(29)
its harmonic coefficients are simply

$`\overline{b}_{l,m}`$ $`=`$ $`\delta _{m,0}b_{l,0}`$ (30)
$`=`$ $`\delta _{m,0}\sqrt{(2l+1)\pi }{\displaystyle _1^1}\overline{b}\left[\mathrm{arccos}(x)\right]P_l(x)dx,`$
where $`P_l(x)`$ is a Legendre polynomial (Appendix A). The use of $`\overline{b}(\widehat{𝒓})`$ allows the definition of a beam-smoothed sky, $`\overline{s}(\widehat{𝒓})`$, given in terms of the true sky, $`s(\widehat{𝒓})`$, by

$$\overline{s}(\widehat{𝒓})=_S\overline{b}\left[\mathrm{arccos}\left(\widehat{𝒓}\widehat{𝒓}^{}\right)\right]s(\widehat{𝒓}^{})d\mathrm{\Omega }^{}.$$
(31)
This convolution is much simpler in harmonic space, and applying equation (48) to equation (31) yields

$$\overline{𝒔}=\overline{\text{B}}𝒔,$$
(32)
where $`\overline{\text{B}}`$ (as distinct from the conversion matrix, B) is a diagonal ‘convolution matrix’ with

$`\overline{B}_{i,i}`$ $`=`$ $`\sqrt{{\displaystyle \frac{4\pi }{2l(i)+1}}}\overline{b}_{l(i),0}`$ (33)
$`=`$ $`2\pi {\displaystyle _1^1}\overline{b}\left[\mathrm{arccos}(x)\right]P_{l(i)}(x)dx,`$
$`l(i)`$ being defined in Appendix A. The simple form of equation (32) is often utilised explicitly in CMB analysis algorithms (e.g. Knox 1995; Hobson et al. 1998; Oh et al. 1999; Stolyarov 2001), but in all cases full sky coverage is – sometimes implicitly – assumed.

Turning to convolution on the incomplete sphere, $`S^{}`$, application of equation (24) to equation (32) yields

$$\overline{𝒔}^{}=\text{A}^\mathrm{T}\overline{\text{B}}𝒔,$$
(34)
where A is defined implicitly in equation (12). In the low-resolution limit equation (25) then gives the cut-sky analogue of equation (32) as

$$\overline{𝒔}^{}=\overline{\text{B}}^{}𝒔^{},$$
(35)
with the new convolution matrix defined by $`\overline{\text{B}}^{}=\text{A}^\mathrm{T}\overline{\text{B}}\text{B}^\mathrm{T}`$.

For higher band-limits no such relation exists, the loss of modes in the cut region rendering the convolution ill-defined. This is simply understood in real space, as the value of $`\overline{s}(\widehat{𝒓})`$ near the edge of $`S^{}`$ is given by an integral that extends several beam widths into the removed region. Thus it is impossible to relate $`\overline{s}(\widehat{𝒓})`$ to $`s(\widehat{𝒓})`$ with $`\widehat{𝒓}S^{}`$. These arguments are true independent of the representation chosen, but in harmonic space they mean that it is impossible to relate $`\overline{𝒔}^{}`$ to $`𝒔^{}`$.

Whilst equation (35) is formally incorrect in high-$`l_{\mathrm{m}ax}`$ cases, it is potentially useful as a practical approximation. It is equivalent to assuming that the signal is given by $`𝒀^{}(\widehat{𝒓})𝒔^{}`$, which implies that $`s(\widehat{𝒓})0`$ in $`SS^{}`$. This is particularly inaccurate if the removed region contains anomalously strong sources, such as the Galactic plane. Nonetheless, equation (35) gives $`s^{}(\widehat{𝒓})`$ correctly for all $`\widehat{𝒓}`$ more than a few beam widths away from the edge of the cut region. However, even if this is an acceptable approximation, there is the further inconvenience that the effective cut-sky beam, $`\overline{\text{B}}^{}`$, is not diagonal, introducing couplings between all the modes.

In short, it is preferable to avoid performing any sort of convolution (or deconvolution) on the cut sky, although it is clear that this situation is encountered in any survey with incomplete sky coverage. The one, albeit trivial, exception to this rule is if the beam is a delta function, or at least the closest approximation to a delta function possible given the band-limit under consideration. In this case $`\overline{\text{B}}=\text{I}`$ and hence, from equation (14), $`\overline{\text{B}}^{}=\stackrel{~}{\text{I}}`$ as well. Equation (35) then implies that $`\overline{𝒔}^{}`$ ($`=𝒔^{}`$) is the true, unconvolved sky map, estimation of which is addressed next.

### 4.2 Map-making

Some of the most important products of the next generation of CMB survey will be high-resolution maps of the sky at each of several frequencies. Such maps can be created in a number of ways, but care must be taken to account for a huge variety of systematics whilst retaining as much information as possible. Both real space (e.g. Wright et al. 1996; Bond et al. 1998; Natoli et al. 2001) and Fourier space (van Leeuwen et al. 2001; Challinor et al. 2001) algorithms have been proposed as being suited to particular apsects of the map-making problem. The resultant uncertainties in the pixel values or harmonic coefficients depend on both the data itself (i.e. the scan strategy, noise properties, etc.) and the map-making algorithms used, and can vary quite markedly from experiment to experiment.

Here only the idealised case of uniform sky coverage is considered as the discussion which follows is not significantly changed by this useful simplification. Under this assumption the optimal estimator for the unsmoothed sky, $`\widehat{𝒔}`$ ($`=𝒔+𝒏`$), would be unbiased (i.e. $`\widehat{𝒔}=𝒔`$) and have covariance given by

$$\text{N}=𝒏𝒏^\mathrm{T}=\left(\widehat{𝒔}𝒔\right)\left(\widehat{𝒔}𝒔\right)^\mathrm{T}=\sigma ^2\overline{\text{B}}^2,$$
(36)
where (cf. Knox 1995) $`\sigma ^2=4\pi s^2/(N_\mathrm{d}t_{\mathrm{o}bs})`$, $`t_{\mathrm{o}bs}`$ is the total observation time of the survey, and $`N_\mathrm{d}`$ is the number of detectors at the frequency in question (all of which are assumed to have the same beam). An important issue at this point is the band-limit chosen. Clearly only a finite analysis is possible in practice, and $`\overline{\text{B}}`$ becomes increasingly singular as $`l_{\mathrm{m}ax}\mathrm{}`$; these two points are related in so far as the sky can never be reconstructed with infinite resolution. The choice of $`l_{\mathrm{m}ax}`$ is somewhat arbitrary, although any value a factor of a few greater than the effective beam width will ensure that $`\overline{\text{B}}^1`$ exists whilst discarding only multipoles that are noise-dominated. The fact that, unlike the useful signal, the noise is not subject to any band-limit is critical to the understanding of the low-$`l_{\mathrm{m}ax}`$ cut-sky power spectrum estimation discussed in Section 4.4.

Note also that, due to the assumption of uniform sky coverage, the covariance matrix is diagonal. Transforming this estimator into real space yields maps with covariance given by

$$n(\widehat{𝒓}_1)n(\widehat{𝒓}_2)=\sigma ^2𝒀^\mathrm{T}(\widehat{𝒓}_1)\overline{\text{B}}^2𝒀(\widehat{𝒓}_2).$$
(37)
As the noise term in the data is not beam-convolved the removal of the beam results in spatial correlations of the noise (as encoded in $`\overline{\text{B}}`$), as well as correlations due to the finite resolution analysis (the sums over spherical harmonics), which are essentially equivalent to pixel smoothing. In the more realistic case of non-uniform sky coverage, the covariance matrix is non-diagonal in both bases, a point discussed further by Oh et al. (1999).

The above estimator for the true sky is closely linked to the more commonly used estimator for the smoothed sky, $`\widehat{\overline{𝒔}}`$. Being related by $`\widehat{\overline{𝒔}}=\overline{\text{B}}^1\widehat{𝒔}`$, it is clear they contain the same information (under the assumption of a symmetric beam). The covariance structure of $`\widehat{\overline{𝒔}}`$ is simpler as the correlations discussed above are not introduced, but $`\widehat{𝒔}`$ is a more natural data object in the context of this discussion as it is the true sky that is of interest. In particular, unsmoothed maps allow more flexibility in applying a sky cut, as the problems with convolution on the incomplete sphere described in Section 4.1.2 do not arise. In practice the best compromise may be to reconstruct the sky convolved with the azimuthally averaged beam, thus creating maps with the simplest covariance structure possible without information loss. This can be done in either real space (e.g. Bond et al. 1998) or harmonic space (Challinor et al. 2001), although the real space pointing matrix is more complicated if beam asymmetry information is included.

In summary, if a survey covers the entire celestial sphere it is preferable to use full-sky frequency maps. However it is possible that small parts of the sky will be missed due to either the scan strategy (cf. Maino et al. 1999) or hardware problems during the survey itself. If this is the case the best unsmoothed map that could be constructed would be larger than the actual observed region, but the errors around the boundary of this area would be very high. An inferential approach is possible, but significant difficulties are encountered, especially in Fourier space (Challinor et al. 2001). Fortunately, it is probable that both MAP and Planck will produce full sky maps at several frequencies, which can then be used to construct maps of the various astrophysical components.

### 4.3 Component separation

The microwave sky consists of several distinct astrophysical components, as listed in Section 1. Fortunately they have sufficiently distinct spectra that they can be separated using multi-frequency data. Given that MAP and Planck will produce maps in five and ten bands, respectively, it should be possible to produce maps of the various components (particularly the CMB) that are relatively free of contamination. As with map-making, a number of algorithms have been put forward for this stage of the data processing, although the main focus has been on Fourier space methods (e.g. Tegmark & Efstathiou 1996; Hobson et al. 1998; Bouchet & Gispert 1999; Prunet et al. 2001; Stolyarov et al. 2001). Aside from the expected statistical isotropy of the CMB signal, one of the reasons for this emphasis has been the simplicity of beam convolution in harmonic space (Section 4.1.2). This is critical if smoothed frequency maps are used as the effective smoothing scale will vary with frequency if the telescope is (close to) diffraction limited. However if unsmoothed maps are used real space component separation methods (e.g. Baccigalupi et al. 2000) must also come into consideration, the optimal choice of basis being less clear.

One common aspect of all the separation techniques referenced above is that much of the (prior) information about both signal and noise correlations is disregarded in order to render the problem computationally feasible. In real space the correlations between nearby pixels are ignored, and in Fourier space it is the mode-mode couplings that are neglected. Surprisingly, these approximations appear to be unimportant in practice – even the Galactic components have been recovered with striking accuracy. The most relevant result to this discussion is the all-sky component separation to $`l_{\mathrm{m}ax}2500`$ presented by Stolyarov et al. (2001), as it provides clear evidence that whatever correlations are present in the full-sky harmonic basis are unimportant – there are some errors close to the Galactic centre, but they are localised, and there is no sign of this affecting the reconstruction globally.

If the Galactic plane is removed prior to component separation this one troublesome region is no longer present in the analysis, but new problems arise. Firstly, smoothed maps (with frequency-dependent beam-widths) cannot be used as input data without inducing errors around the edges of $`S^{}`$ due to the ill-defined nature of convolution on the cut sky (Section 4.1.2). Even if such errors are deemed acceptable (e.g. Prunet et al. 2001) or beam-deconvolved maps are used, the transformation described in equation (24) completely changes the correlation structure of the harmonics. In particular, the signal-signal correlation matrices are non-diagonal for all components, including random fields like the CMB (see Section 4.4). Whereas the couplings between the spherical harmonic coefficients can apparently be disregarded, this has not been demonstrated for these induced correlations in the orthonormal basis. Prunet et al. (2001) performed cut-sky component separation including them in full, but were thus limited to $`l_{\mathrm{m}ax}500`$, the computational task being made considerably larger.

In real space the application of a cut is trivial, provided that beam-deconvolved maps are used, as it simply requires that pixels in the removed region be ignored. Thus the Baccigalupi et al. (2000) method should be well suited to a cut-sky analysis.

Given that realistic component separation simulations have only recently become available, it is likely that important developments in this field will be made in the near future. For the moment, however, it appears that separation can usefully be performed on either the full or cut sky without introducing catastrophic errors. Provided the models of microwave emissions from the Galactic plane used in the above simulations are sufficiently realistic, it may thus be preferable to generate the full-sky maps of the various astrophysical components, retaining the option of masking unwanted regions at a later stage.

### 4.4 Power spectrum estimation

If the fluctuations in the early universe were the result of inflation (e.g. Linde 1990) then the CMB is expected to be a Gaussian random field, the statistical properties of which can be specified completely by its angular power spectrum, $`C_l`$. Even if this is not the case, the power spectrum should encode much of the cosmological information present. It is thus unsurprising that, as with map-making and component separation (Sections 4.2 and 4.3, respectively), many different methods of power spectrum estimation have been developed (e.g. Tegmark 1997; Górski 1994; Bond et al. 1998; Oh et al. 1999; Szapudi et al. 2001; Wandelt et al. 2001; Hivon et al. 2001). Further, sky cuts have been incorporated into many of these algorithms as it seems certain that the strength of the Galactic microwave emissions will prevent the CMB from ever being accurately measured in this region. Due to the proliferation of papers on this subject, this discussion of power spectrum estimation is limited to a description of a maximum likelihood formalism using the orthonormal basis functions described in Section 2, with reference to how their behaviour differs in the low- and high-resolution regimes.

#### 4.4.1 Maximum likelihood formalism

The most powerful method of power spectrum estimation is maximum likelihood (e.g. Press et al. 1992), although this has only been implemented to MAP resolution to date (Oh et al. 1999). By invoking Gaussian statistics for both the CMB and the noise, it is possible to write down the exact likelihood for the observed map (e.g. Górski 1994; Borrill 1999). On the full sky the effective data-vector is $`\widehat{𝒔}=𝒔+𝒏`$ (i.e. the estimator for the true sky, of the form discussed in Section 4.2 and not the quantity being estimated here) with the non-cosmological $`l=0`$ and $`l=1`$ modes removed. The full likelihood is given by

$$\frac{\mathrm{d}p}{\mathrm{d}\widehat{𝒔}}=\frac{1}{\sqrt{(2\pi )^{i_{\mathrm{m}ax}4}\mathrm{d}et(\text{S}+\text{N})}}\mathrm{exp}\left[\frac{1}{2}\widehat{𝒔}^\mathrm{T}(\text{S}+\text{N})^1\widehat{𝒔}\right],$$

(38)
where S and N are the signal and noise covariance matrices, respectively. The assumption of Gaussianity implies that $`S_{i,i^{}}=\delta _{i,i^{}}C_{l(i)}`$, where $`C_l`$ is the CMB power spectrum and $`l(i)`$ is defined in Appendix A. The form of N is determined by a combination of the survey method and the data processing up to this point, but is unlikely to have the simple form of equation (36) due to the imperfect component separation. The maximum likelihood calculation consists of finding an estimator for the underlying power spectrum, $`\widehat{C_l}`$, such that equation (38) is maximised, and there are a number of algorithms for finding this quantity (e.g. Bond et al. 1998; Oh et al. 1999).

The maximum likelihood formalism on the cut sky takes the same form as on the full sky, but with the data-vector $`\widehat{𝒔}^{}=\text{A}^\mathrm{T}\widehat{𝒔}`$ and the covariance matrices suitably transformed to give (cf. Górski 1994)

$`{\displaystyle \frac{\mathrm{d}p}{\mathrm{d}\widehat{𝒔}^{}}}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{(2\pi )^{i_{\mathrm{m}ax}^{}4}\mathrm{d}et(\text{S}^{}+\text{N}^{})}}}`$ (39)
$`\times `$ $`\mathrm{exp}\left[{\displaystyle \frac{1}{2}}\widehat{𝒔}^\mathrm{T}(\text{S}^{}+\text{N}^{})^1\widehat{𝒔}^{}\right],`$
with the four modes containing information on the monopole and dipole (Section 2.2) again excluded. The signal covariance matrix is subject to a simple similarity transform, $`\text{S}^{}=\text{A}^\mathrm{T}\text{S}\text{A}`$, but the same is not true for the noise covariance matrix as the noise field is not band-limited (a fact critical to the use use of a cut-sky analysis at low resolution, as discussed below). The coupling of the cut-sky modes makes the maximisation of equation (39) non-trivial (cf. Oh et al. 1999), even if $`\text{S}+\text{N}`$ is diagonal on the full sky. Nonetheless it is useful to work under this idealised assumption in order to see how the application of the cut ensures that the maximum likelihood solution is independent of the Galactic signal; the manner in which this is achieved is quite different in the low- and high-resolution cases.

#### 4.4.2 Low-resolution analysis

The effect of a sky cut on power spectrum estimation is not entirely obvious in the low-$`l_{\mathrm{m}ax}`$ case in which the coupling matrix (Section 2.1) is invertible. The effective band-limit, produced by the combined effects of the beam and noise (Sections 4 and 4.2), means the signal over the whole sky (including e.g. the Galactic plane) is encoded in the cut-sky coefficients. Thus the application of a cut would be redundant were it not for the presence of non-band-limited noise which cannot be characterised properly a finite harmonic analysis.

Applying the incomplete spherical transform defined in equation (23) to a purely white noise field $`n(\widehat{𝒓})`$ \[i.e. $`n(\widehat{𝒓})=0`$ and $`n(\widehat{𝒓}_1)n(\widehat{𝒓}_2)=\delta (\widehat{𝒓}_1\widehat{𝒓}_2)\sigma ^2`$; cf. equation (37)\] gives cut-sky harmonic coefficients with $`𝒏^{}=\mathrm{𝟎}`$ and

$$\text{N}^{}=𝒏^{}𝒏_{}^{}{}_{}{}^{\mathrm{T}}=\sigma ^2\stackrel{~}{\text{I}}.$$
(40)
Projecting back into real space gives a field $`n^{}(\widehat{𝒓})`$ which satisfies $`n^{}(\widehat{𝒓})=0`$ and

$$n^{}(\widehat{𝒓}_1)n^{}(\widehat{𝒓}_2)=\sigma ^2𝒀^\mathrm{T}(\widehat{𝒓}_1)\text{C}^1𝒀(\widehat{𝒓}_2).$$
(41)
On the full sky the same procedure (i.e. a finite spherical harmonic analysis followed by a transformation back into real space) would yield a noise field with covariance structure given by

$`\widehat{n}(\widehat{𝒓}_1)\widehat{n}(\widehat{𝒓}_2)`$ $`=`$ $`\sigma ^2𝒀^\mathrm{T}(\widehat{𝒓}_1)𝒀(\widehat{𝒓}_2)`$ (42)
$`=`$ $`{\displaystyle \frac{\sigma ^2}{4\pi }}{\displaystyle \underset{l=0}{\overset{l_{\mathrm{m}ax}}{}}}(2l+1)P_l(\widehat{𝒓}_1\widehat{𝒓}_2),`$
where $`P_l(x)`$ is a Legendre polynomial (Appendix A). Taking the limit $`\widehat{𝒓}_2\widehat{𝒓}_1`$ this implies that $`\widehat{n}^2(\widehat{𝒓})=i_{\mathrm{m}ax}\sigma ^2/(4\pi )`$, which represents smoothing relative to the original noise field caused by the use of a finite analysis.

Whilst this smoothing occurs in both the cut- and full-sky formalisms, the presence of $`\text{C}^1`$ in the former case \[equation (41)\] implies a spatial dependence. As can be seen from Fig. 4, the noise in the cut region is greatly increased, which is a natural way of formally encoding the qualitative fact that, for whatever reason, the data in the cut is contaminated by more than just the white noise field. Thus, despite the invertibility of the coupling matrix (and the band-limited cut-sky analysis), the application of a cut has the desired effect of greatly reducing the impact of any spurious signal, such as the Galaxy. However Fig. 4 also implies that a similar effect could be achieved without performing a cut (and hence leaving the signal unchanged), instead adding a high level of artificial noise in the offending region(s). Finally, it is important to note that the dual assumptions used in the derivation of equation (41) – uniform noise and no beam – are unrealistic, but the manner in which a low-resolution cut-sky analysis works is the same in less idealised scenarios.

#### 4.4.3 High-resolution analysis

The high-resolution case is more straightforward, as the application of the cut results in a data-vector, $`\widehat{𝒔}^{}`$, which contains little information about the removed region. This is quite distinct from the low-resolution case discussed above, in that here it is the predominantly the signal that is changed, rather than the noise. That said, the noise close to the boundary of the cut is increased in the same manner as explained above. This has the same effect as the apodizing function formalism described by Tegmark (1997), downweighting points around which there is not full correlation information.

Another difference between the low- and high-resolution analyses is that $`\widehat{𝒔}^{}`$ is smaller than $`\widehat{𝒔}`$, from Section 2.1.1. Although this does not result in any significant computational saving, it serves to emphasise the information loss associated with removing part of the sky, and is an independent derivation of the fact that the uncertainties in the estimated power spectrum increase as $`4\pi /\mathrm{\Omega }_S^{}`$ (cf. Hobson & Magueijo 1996; Tegmark 1997).

## 5 Conclusions

The upcoming microwave surveys will require a cut-sky analysis to prevent the strong Galactic emissions from contaminating the CMB signal. The spherical harmonics are non-orthogonal on the cut-sphere, but an orthonormal basis set can be constructed from them using SVD-based techniques (Section 2). The application of the resultant conversion matrix to the conventional multipoles results in cut-sphere harmonics that contain only the desired information. In the low-resolution case the influence of the Galaxy is reduced by increasing the effective noise in the cut; in the high-resolution limit the orthonormal basis functions can model the infinitely sharp cut sufficiently well that they have no support in the removed region. It is also important to note that the cut should probably only be applied after beam-deconvolution has been attempted, as convolution is ill-defined on the incomplete sphere.

The algorithms described here were implemented to Legendre multipoles of $`l_{\mathrm{m}ax}2500`$ for a constant latitude cut, in which case the coupling matrix of the spherical harmonics is block-diagonal. At present computational limitations make a general orthogonalisation impractical for $`l_{\mathrm{m}ax}200`$, although there are some possibilities to extend this. For instance only $`1`$ per cent of the coupling matrix contains significant information if the cut is well-chosen (e.g. rectangular in $`\theta `$ and $`\varphi `$) and so sparse matrix techniques should thus allow orthogonalisation to $`l_{\mathrm{m}ax}1000`$ in this case.

Another requirement is orthogonalisation of tensor basis functions on the incomplete sphere, as both the MAP and Planck satellites will measure polarization. The resultant formalism is more complicated, but the same general principles hold; this issue is explored further by Lewis, Challinor & Turok (2001).

## Acknowledgments

This paper benefited from useful discussions with several members of the Planck collaboration, in particular Mark Ashdown, François Bouchet, Martin Bucher, Rob Crittenden, Jacques Delabrouille, George Efstathiou, Krzysztof Górski, Floor van Leeuwen and Ben Wandelt. DJM was funded by PPARC. ADC acknowledges a PPARC Postdoctoral Research Fellowship. MPH acknowledges a PPARC Advanced Fellowship.

## Appendix A Spherical harmonics

The spherical harmonics form a complete set of orthonormal basis functions over the entire sphere. They are most commonly defined as complex functions (e.g. Landau & Lifshitz 1976; Brink & Satchler 1993), but it is more convenient to use real harmonics in this application. Adapting the notation of Górski (1994), the real spherical harmonics are given by

$$Y_{l,m}(\widehat{𝒓})=Y_{l,m}(\theta ,\varphi )=\lambda _{l,|m|}[\mathrm{cos}(\theta )]s_m(\varphi ),$$
(43)
where $`l0`$ and $`|m|l`$ and

$$s_m(\varphi )=\{\begin{array}{ccc}\sqrt{2}\mathrm{sin}(|m|\varphi ),\hfill & \mathrm{i}f\hfill & m<0,\hfill \\ & & \\ 1,\hfill & \mathrm{i}f\hfill & m=0,\hfill \\ & & \\ \sqrt{2}\mathrm{cos}(|m|\varphi ),\hfill & \mathrm{i}f\hfill & m>0,\hfill \end{array}$$
(44)
implying that $`_0^{2\pi }s_m(\varphi )s_m^{}(\varphi )d\varphi =2\pi \delta _{m,m^{}}`$. For $`0ml`$ and $`1x1`$ the normalised associated Legendre functions are defined by

$$\lambda _{l,m}(x)=\sqrt{\frac{2l+1}{4\pi }\frac{(lm)!}{(l+m)!}}P_{l,m}(x).$$
(45)
Hence $`_1^1\lambda _{l,m}(x)\lambda _{l^{},m}(x)dx=\delta _{l,l^{}}/(2\pi )`$. Under the same conditions the (unnormalised) associated Legendre functions are given by555This definition, with the $`(1)^m`$ term, correpsonds to that given by Abramowitz & Stegun (1971) and Gradshteyn & Ryzhik (2000) but differs from that used by Arfken (1985) and Brink & Satchler (1993).

$$P_{l,m}(x)=(1)^m\left(1x^2\right)^{m/2}\frac{\mathrm{d}^m}{\mathrm{d}x^m}P_l(x),$$
(46)
with the Legendre polynomials given by

$$P_l(x)=\frac{(1)^l}{2^ll!}\frac{\mathrm{d}^l}{\mathrm{d}x^l}\left(1x^2\right)^l.$$
(47)
A real field on the sphere, $`a(\widehat{𝒓})`$, can be expanded in terms of spherical harmonic coefficients, given by

$$a_{l,m}=_SY_{l,m}(\widehat{𝒓})a(\widehat{𝒓})d\mathrm{\Omega }.$$
(48)
This can be inverted to give

$$a(\widehat{𝒓})=\underset{l=0}{\overset{l_{\mathrm{m}ax}}{}}\underset{m=l}{\overset{l}{}}a_{l,m}Y_{l,m}(\widehat{𝒓}),$$
(49)
provided that $`l_{\mathrm{m}ax}\mathrm{}`$, due to the orthonormality of the spherical harmonics on the full sphere:

$$_SY_{l,m}(\widehat{𝒓})Y_{l^{},m^{}}(\widehat{𝒓})d\mathrm{\Omega }=\delta _{l,l^{}}\delta _{m,m^{}}.$$
(50)
If a finite $`l_{\mathrm{m}ax}`$ is used this inversion is no longer possible for general $`a(\widehat{𝒓})`$, although it does still hold for band-limited functions666A band-limited function can, by (somewhat circular) definition, be constructed from a finite sum over spherical harmonics..

Whilst the two indices $`l`$ and $`m`$ have quite distinct interpretations it is convenient to combine them into a single index, $`i`$, which allows the definition of vectors $`𝒀(\widehat{𝒓})=Y_{i(l,m)}(\widehat{𝒓})`$ and $`𝒂=a_{i(l,m)}`$. Two obvious indexing schemes present themselves: grouping in $`l`$ and $`m`$. The first, as introduced by Górski (1994), is natural for power spectrum estimation and very simple:

$$i(l,m)=l^2+l+m+1.$$
(51)
The two ‘inverses’ of this relationship are

$$l=\mathrm{i}nt\left[(i1)^{1/2}\right]$$
(52)
and

$$m=i(l^2+l+1).$$
(53)
The second choice of ordering is useful in cases of azimuthal symmetry in which the orthogonality expressed in equation (44) is maintained, and grouping in $`m`$ is achieved by defining

$$i(l,m)=\{\begin{array}{ccc}l+m+1\hfill & & \\ +(l_{\mathrm{m}ax}+m)(l_{\mathrm{m}ax}+m+1)/2,\hfill & \mathrm{i}f\hfill & m0,\hfill \\ & & \\ ll_{\mathrm{m}ax}+(l_{\mathrm{m}ax}+1)^2\hfill & & \\ (l_{\mathrm{m}ax}m)(l_{\mathrm{m}ax}m+1)/2,\hfill & \mathrm{i}f\hfill & m>0.\hfill \end{array}$$

(54)
The ‘inverses’ in this case are given by

$$m=\{\begin{array}{ccc}\multicolumn{3}{c}{\mathrm{i}nt\left[l_{\mathrm{m}ax}+\frac{1}{2}\left(\sqrt{8i+1}3\right)\right],}\\ \multicolumn{3}{c}{\mathrm{if}i(l_{\mathrm{m}ax}+1)(l_{\mathrm{m}ax}+2)/2,}\\ & & \\ \multicolumn{3}{c}{\mathrm{i}nt[[l_{\mathrm{m}ax}\frac{1}{2}\left\{\sqrt{8[(l_{\mathrm{m}ax}+1)^2i+1]+1}3\right\}]],}\\ \multicolumn{3}{c}{\mathrm{if}i>(l_{\mathrm{m}ax}+1)(l_{\mathrm{m}ax}+2)/2}\end{array}$$

(55)
and

$$l=\{\begin{array}{ccc}i[m+1\hfill & & \\ +(l_{\mathrm{m}ax}+m)(l_{\mathrm{m}ax}+m+1)/2],\hfill & \mathrm{i}f\hfill & m0,\hfill \\ & & \\ i[l_{\mathrm{m}ax}+(l_{\mathrm{m}ax}+1)^2\hfill & & \\ (l_{\mathrm{m}ax}m)(l_{\mathrm{m}ax}m+1)/2],\hfill & \mathrm{i}f\hfill & m>0.\hfill \end{array}$$
(56)
Other indexing schemes have been used in the more specific case of simulated Planck data-sets in which the sky coverage is periodic in azimuth (van Leeuwen, private communication), but are beyond the scope of this paper.

## Appendix B Integration of products of associated Legendre functions

In Section 2.3 integrals of the form

$$I_{l,l^{},m}(x_1,x_2)=_{x_1}^{x_2}\lambda _{l,m}(x)\lambda _{l^{},m}(x)dx$$
(57)
arose; here the $`\lambda _{l,m}(x)`$ are the normalised associated Ledengre functions, defined in equation (45), and $`m`$ is assumed to be non-negative. These integrals can be evaluated quickly and accurately using a combination of closed formulæ and recursion relations.

The associated Legendre functions, $`P_{l,m}(x)`$ \[defined in equation (46)\], are solutions of the ordinary differential equation (e.g. Arfken 1985)

$$\frac{\mathrm{d}}{\mathrm{d}x}\left[\left(1x^2\right)\frac{\mathrm{d}P_{l,m}}{\mathrm{d}x}\right]+\left[l(l+1)\frac{m^2}{1x^2}\right]P_{l,m}(x)=0.$$
(58)
Multiplying this equation by $`P_{l^{},m}(x)`$ and integrating (from $`x_1`$ to $`x_2`$) by parts twice yields

$$(ll^{})(l+l^{}+1)_{x_1}^{x_2}P_{l,m}(x)P_{l^{},x}(x)dx$$
(59)
$$=[\left(1x^2\right)P_{l,m}(x)\frac{\mathrm{d}P_{l^{},m}}{\mathrm{d}x}\left(1x^2\right)P_{l^{},m}(x)\frac{\mathrm{d}P_{l,m}}{\mathrm{d}x}]|_{x=x_1}^{x=x_2}.$$

This is a reflection of the standard result that integrals of solutions of a self-adjoint differential equation \[as equation (58) is\] can be expressed as boundary terms (e.g. Arfken 1985). The derivatives in equation (59) can be removed by using the standard recursion relationship (e.g. Gradshteyn & Ryzhik 2000)

$$\left(1x^2\right)\frac{\mathrm{d}P_{l,m}}{\mathrm{d}x}=(l+m)P_{l1,m}(x)lxP_{l,m}(x)$$
(60)
to yield, for $`ll^{}`$,

$$_{x_1}^{x_2}P_{l,m}(x)P_{l^{},x}(x)dx=\frac{1}{(ll^{})(l+l^{}+1)}$$

$`\times `$ $`[(l^{}+m)P_{l,m}(x)P_{l^{}1,m}(x)`$ (61)
$`+(ll^{})xP_{l,m}(x)P_{l^{},m}(x)`$
$`(l+m)P_{l1,m}(x)P_{l^{},m}(x)\left]\right|_{x=x_1}^{x=x_2}.`$
Note that the first term must be omitted if $`l^{}=m`$ and that the third term must be omitted if $`l=m`$; these Legendre functions are implicitly zero from equation (46). Finally, this can be normalised according to equation (45), giving

$$I_{l,l^{},m}(x_1,x_2)=\frac{1}{(ll^{})(l+l^{}+1)}$$
(62)
$`\times `$ $`[\sqrt{{\displaystyle \frac{2l^{}+1}{2l^{}1}}(l^2m^2)}\lambda _{l,m}(x)\lambda _{l^{}1,m}(x)`$
$`+(ll^{})x\lambda _{l,m}(x)\lambda _{l^{},m}(x)`$
$`\sqrt{{\displaystyle \frac{2l+1}{2l1}}(l^2m^2)}\lambda _{l1,m}(x)\lambda _{l^{},m}(x)\left]\right|_{x=x_1}^{x=x_2}.`$
An alternative derivation of this result was presented by Wandelt et al. (2001); it is also in principle equivalent to Eq. 5.9(13) of Varshalovich, Moskalev & Khersonskii (1988), but their application of equation (60) is in error.

For the case $`l=l^{}`$ a recursion relation is required, starting with $`l=m`$. Combining equations (45) and (46),

$$\lambda _{m,m}(x)=(1)^m\sqrt{\frac{2m+1}{4\pi }}(2m1)!!\left(1x^2\right)^{m/2},$$
(63)
where $`n!!=1\times 3\times \mathrm{}\times (n2)\times n`$ for odd $`n`$. Integrating by parts and using equation (63) again gives

$$I_{m,m,m}(x_1,x_2)$$
(64)
$$=\{\begin{array}{ccc}\frac{x_2x_1}{4\pi },\hfill & \mathrm{i}f\hfill & m=0,\hfill \\ & & \\ I_{m1,m1,m1}(x_1,x_2)+\frac{x\lambda _{m,m}^2(x)|_{x_1}^{x_2}}{2m+1},\hfill & \mathrm{i}f\hfill & m>0.\hfill \end{array}$$

Moving to $`l=m+1`$, the standard relationship (e.g. Gradshteyn & Ryzhik 2000) that

$$\lambda _{m+1,m}(x)=(2m+3)x\lambda _{m,m}(x)$$
(65)
combines with equation (63) to give

$`I_{m+1,m+1,m}(x_1,x_2)`$ $`=`$ $`I_{m,m,m}(x_1,x_2)`$ (66)
$``$ $`{\displaystyle \frac{2m+2}{2m+3}}x\lambda _{m+1,m+1}^2(x)|_{x_1}^{x_2},`$
where the first term is given in equation (64).

The last step is to derive a recursion relation relating $`I_{l,l,m}(x_1,x_2)`$ to $`I_{l1,l1,m}(x_1,x_2)`$ and $`I_{l2,l2,m}(x_1,x_2)`$. Eq. (C22) of Wandelt et al. (2001) gives a four-term recursion to obtain $`I_{l,l^{},m}(x_1,x_2)`$; it can be applied successively (once swapping $`l`$ and $`l^{}`$) to obtain

$$I_{l,l,m}(x_1,x_2)=\frac{2m^22l^2+2l1}{(lm)(l+m)}I_{l1,l1,m}(x_1,x_2)$$

$``$ $`{\displaystyle \frac{(l1m)(l1+m)}{(lm)(l+m)}}I_{l2,l2,m}(x_1,x_2)`$
$``$ $`\sqrt{{\displaystyle \frac{2l1}{2l+1}}{\displaystyle \frac{1}{(lm)(l+m)}}}`$
$`\times \left(1x^2\right)\lambda _{l,m}(x)\lambda _{l1,m}(x)|_{x=x_1}^{x=x_2}`$
$`+`$ $`\sqrt{{\displaystyle \frac{2l1}{2l3}}{\displaystyle \frac{(l1m)(l1+m)}{(lm)^2(l+m)^2}}}`$
$`\times \left(1x^2\right)\lambda _{l1,m}(x)\lambda _{l2,m}(x)|_{x=x_1}^{x=x_2}.`$
In summary, equation (62) can be used to evaluate all $`I_{l,l^{},m}(x_1,x_2)`$ for which $`ll^{}`$, and equations (64), (66) and (B) combine to give all $`I_{l,l,m}(x_1,x_2)`$ recursively.

## Appendix C Treatment of non-cosmological modes

All the cosmological information encoded in the CMB is expected to be contained in the $`l2`$ modes; the $`l=0`$ mode in an isotropic universe can be normalised arbitrarily and the $`l=1`$ modes can be set to zero by adopting an appropriate reference frame. Nonetheless, observations of the microwave sky will yield non-zero monopole and dipole values for a number of reasons (e.g. the observer’s motion; Galactic emission; extra-Galactic point-sources). Hence these low-order modes must be included in the analysis of CMB data, but should be kept separate from the cosmological modes, as is naturally the case if spherical harmonic coefficients are used to describe the data. It is also important to note that the properties of the basis functions themselves are unimportant – the essential requirement is that only four of the cut-sky harmonic coefficients contain information on the unwanted modes.

The method of orthogonalisation summarised in equations (11), (12) and (14) does not explicitly impose any particular structure on the conversion matrix, A \[which relates harmonic coefficients on the incomplete sphere to those on the full sphere by $`𝒂^{}=\text{A}^\mathrm{T}𝒂`$; equation (24)\]. The non-cosmological modes are kept separate from the cosmological modes if the first four columns of $`\text{A}^\mathrm{T}`$ have only zeros from the fifth row on, assuming the full-sky harmonic coefficients are indexed using $`l`$-ordering (Appendix A). This is achieved naturally if $`\text{A}^\mathrm{T}`$ is constructed to be upper triangular, as in the case of the Cholesky decomposition described in Section 2.2.1. The other decomposition methods discussed in Section 2.2 do not share this property, and so the resultant conversion matrices must be adjusted explicitly.

One way of achieving this is to use a partial Householder transform (e.g. Press et al. 1992). The last $`i_{\mathrm{m}ax}^{}i`$ elements of the $`i`$th column of a general $`i_{\mathrm{m}ax}^{}\times i_{\mathrm{m}ax}`$ matrix can be set to zero by the transformation $`\text{M}^{}=\text{P}_i\text{M}`$, with the orthogonal Householder matrix defined by

$$\text{P}_i=\stackrel{~}{\text{I}}2\frac{𝒎_i𝒎_i^\mathrm{T}}{𝒎_i^\mathrm{T}𝒎_i},$$
(67)
where $`𝒎_i`$ is given by

$$(m_i)_j=\{\begin{array}{ccc}0,\hfill & \mathrm{i}f\hfill & ji,\hfill \end{array}$$
(68)
and $`i\mathrm{min}(i_{\mathrm{m}ax}^{},i_{\mathrm{m}ax}`$) is assumed. Provided that the Householder matrix applied to B is that generated from $`\text{A}^\mathrm{T}`$, the transformations $`\text{A}_{}^{}{}_{}{}^{\mathrm{T}}=\text{P}_i\text{A}^\mathrm{T}`$ and $`\text{B}^{}=\text{P}_i\text{B}`$ leave equations (12) and (14) unaffected as $`\text{P}_i`$ is orthogonal by construction. Applying $`\text{P}_1`$, $`\text{P}_2`$, $`\text{P}_3`$ and then $`\text{P}_4`$ to the successively updated $`\text{A}^\mathrm{T}`$ ensures that the $`l=0`$ and $`l=1`$ modes influence only the first four cut-sky harmonic coefficients, as required. This procedure could be continued, moving $`\text{A}^\mathrm{T}`$ successively closer to upper triangular form, although this cannot be achieved in full as $`\text{A}^\mathrm{T}`$ has more columns than rows.

Special mention must be made of the constant latitude cut case, the symmetry of which can only be utilised if the spherical harmonics are indexed using $`m`$-ordering (Appendix A). In this case only the $`m=0`$ and $`m=\pm 1`$ blocks have any contribution from the monopole or dipole, and each can be treated separately. Further, the ordering within these blocks is such that the non-cosmological modes are in the first rows, and so the above algorithm can be applied to each of three blocks as is. The only slight inconvenience is that it is no longer the first four cut-sky modes that contain the non-cosmological information, and the relevant modes must be flagged explicitly.