text:
# Damped Lyman 𝜶 systems at 𝒛=0.68 and 𝒛=1.15 towards HE 1122-1649 and HE 0515-4414 Based on observations performed at ESO/La Silla, on observations with the NASA/ESA Hubble Space Telescope and on observations at Keck Observatory.

## 1 Introduction

Damped Lyman $`\alpha `$ absorbers (DLA) are the population of quasar absorption systems showing the highest neutral hydrogen column density ($`N`$(H i )$`2\times 10^{20}`$ cm-2). At high redshift they contain the bulk of the cosmological mass of neutral gas in the Universe and are believed to be progenitors of current galaxies (e.g. Wolfe et al. 1986), although their nature and properties still remain unclear. They provide information about the chemical enrichment history of galaxies which can be contrasted with the studies of stellar populations in the local galaxies. Several samples of DLA at $`z>1.7`$ have been studied thoroughly (e.g Pettini et al. 1990, Lu et al. 1996, Pettini et al. 1997a,b) but very few systems are known at intermediate or low redshift (Pettini & Bowen 1997, and references therein; Boissé. et al. 1998, de la Varga & Reimers 1998, Pettini et al. 1999a,b). To interpret the high-$`z`$ data correctly, it is essential to follow the evolution of the DLA population to the present time. Moreover, at low $`z`$ it is possible to search for luminous counterparts via imaging (e.g. Le Brun et al. 1997, Steidel et al. 1994, 1995, 1997, Rao & Turnshek 1998).

We herein present observations of two DLA, at $`z=0.68`$ and $`z=1.15`$, in the line of sight to HE 1122$``$1649 and HE 0515$``$4414 , respectively. The combination of HST spectra with high resolution spectra (CASPEC and HIRES) allows a detailed study of metal abundances and kinematics at low redshift.

Total metal abundances provide information about the population of absorbers which constitute DLA. Measurements of Zn and Cr show abundances of 1/15 solar, but with a considerable spread (Pettini et al. 1997b). The metallicity is not significantly higher for the few systems found at $`z<1.5`$ (Pettini et al. 1999a,b). The underabundance of Cr relative to Zn (Pettini et al. 1997a, 1999a,b) suggests that dust is also present, though the depletion level seems to be lower than in the Galaxy. The same conclusion was drawn from reddening studies of QSOs with and without DLA in their line of sight (e.g. Pei et al. 1991). The interplay between nucleosynthesis and dust depletion from the study of relative abundance patterns in DLA is still unclear (cf. Lu et al. 1996, Prochaska & Wolfe 1996, Pettini et al. 1997b, Kulkarny et al. 1997) and, however, crucial to understand the evolution and chemical enrichment of galaxies.

Studying the kinematics of the absorption line profiles should clarify the nature of DLA and could help to discern among theories of galaxy formation: are DLA related to rotating disks or rather to irregular protogalactic clumps merging, as expected in scenarios of hierarchical cosmogony? (cf. Prochaska & Wolfe 1997a,b,1998, Haehnelt et al. 1998, McDonald & Miralda-Escudé submitted). Comparing profiles at low and high $`z`$ could help to understand whether the asymmetries arise from a mixture of kinematical structures (random motions, rotation, infall, and merging…), or whether a trend can be recognized.

Our paper is organized as follows: in Sect. 2 we describe the observations and the data reduction. The analysis of the absorption lines is discussed in Sect. 3, the metal abundances in Sect. 4, and the relative abundances in Sect. 5. The neutral carbon and molecular content are studied in Sects. 6 and 7, respectively. In Sect. 8 the kinematics of the absorption profiles is considered, and in Sect. 9 we study the highly-ionized gas related to the systems. Finally, concluding remarks are summarized in Sect. 10.

## 2 Observations and Data Reduction

### 2.1 HST Spectra

Ultraviolet spectra of HE 1122$``$1649 ($`z_{em}=2.4`$, V=16.5, Reimers et al. 1995) were obtained with the Faint Object Spectrograph (FOS) onboard the Hubble Space Telescope (HST) in 1995. With gratings G270H and G190H (red digicon and entrance aperture of $`3.7^{\prime \prime }\times 3.7^{\prime \prime }`$), the wavelength range from 1600 to 3280 Å was covered with a maximum signal-to-noise ratio of 28. The spectral resolution is FWHM = 1.4 and 2.0 Å for G190H and G270H, respectively (Schneider et al. 1993). The wavelength calibration of the pipeline reduced data was examined by comparing the overlap region between the G190H and the G270H spectra, and also with the Keck spectra, which have an overlap of 100 Å with the FOS spectra.

The spectra were corrected for interstellar reddening using the Seaton law (Seaton 1979) with $`E(BV)=0.08`$ , corresponding to $`N`$(H i ) =4.16$`\times 10^{20}`$ cm-2. The high density of absorption lines complicates the measurement of a local continuum. Regions apparently free of absorption were selected by checking the error in the mean flux for consistency with the noise, and cubic splines were fitted to the mean flux values in those regions.

### 2.2 Keck HIRES Spectra

Optical data of HE 1122$``$1649 were taken with the HIRES spectrograph (Vogt et al. vog94 (1994)) on the W.M. Keck1 10-m telescope with a total of 3 hours exposure in April 1997. The spectral range was from 3150 to 4710 Å with a S/N of 19/35 at 3500/4500 Å. We used a 1.14” x 7.0” slit, which gave a resolution of 8 km s-1 FWHM (R=37500) and adequate sky coverage, and used the UV cross-disperser to increase efficiency at wavelengths below 4200 Å. The images were processed and the spectra were optimally extracted using an automated package specifically designed by T. Barlow for HIRES spectra. Thorium-Argon lamp images were obtained immediately after the observations to provide wavelength calibrations in each echelle order. The root-mean-square residuals in the wavelength calibration for each echelle order were less than 0.3 km s-1 . All wavelengths are vacuum values in the heliocentric frame. Each echelle order was continuum fit with a Legendre polynomial to normalize the intrinsic QSO flux to unity.

### 2.3 CASPEC Spectra

Further optical spectra were obtained with CASPEC (Long Camera) at the 3.6 m telescope at ESO/La Silla. The data reduction was done with the ECHELLE software program available in MIDAS, supplemented by software programmes for an optimal extraction of the echelle orders developed by S. Lopez. For wavelength calibration purposes, Th-Ar spectra were taken immediately before and after each quasar observation. The uncertainty in the calibration, $`\mathrm{\Delta }\lambda `$, is given by the rms residuals after fitting third order polynomials to the lines identified in the arc lamp observations. For all the CASPEC observations we obtain $`\mathrm{\Delta }\lambda \mathrm{\Gamma }<0.2`$ of a resolution element. The wavelength scales were corrected to vacuum-heliocentric values and individual observations rebinned to the same wavelength scale and combined by weighting by the inverse variance. The continuum for each echelle order was determined by iterative cubic spline fitting to regions free of absorption lines.

HE 1122$``$1649 :
In the observations of HE 1122$``$1649 in 1994, 43 echelle orders covered the wavelength range from 3750 to 5250 Å. The achieved resolving power, measured from the emission lines of the Th-Ar exposures, was R=17 000, and the signal-to-noise ratio for the combined spectrum was 34, 43 and 53 at $`\lambda _{\mathrm{obs}}`$= 4000, 4500 and 5000 Å, respectively. In 1996, 48 orders covered from 3320 to 4600 Å with S/N of 8/29 at 3500/4500 Å and R=20 000.

HE 0515$``$4414 :
We observed HE 0515$``$4414 using CASPEC in October 1996 from 3530 to 4830 Å with R=21 500, and 4865 to 6200 Å with R=24 000. Further observations in November 1997 provided us with the ranges from 3250 to 4560 Å (R=27 000), 5430 to 7870 Å (R=28 600) and 6140 to 8570 Å (R=29 000). The combined spectrum covering the spectral range from 3250 to 8570 (with a small gap between 4830 and 4865 Å) has a maximum signal-to-noise ratio of $`55`$.

## 3 Analysis of the Absorption Lines

Using the FITLYMAN program in the MIDAS package (Fontana & Ballester 1995), the absorption lines were fitted with Voigt profiles convolved with the instrumental profile, simultaneously determining the redshift, column density, $`b`$ values and associated errors, while minimizing a $`\chi ^2`$ parameter matrix. Line parameters were taken from Morton (1991) and from Verner et al. (1994). Different transitions from the same ion were fitted simultaneously, and different ions and elements were fitted independently. We chose the minimum number of components which gave an adequate fit. Several fits were performed for the same complex by shifting slightly the initial $`z`$ values and testing new sets. The considerable arbitrariness in the choice of components does not seriously affect the total column densities derived for most ions, as the presence of several transitions allow us to discern blends and hidden components –at least in the high resolution spectra– by means of the apparent optical depth method (Savage & Sembach 1991).

### 3.1 DLA at $`z=0.68`$ towards HE 1122$``$1649

The profile of the damped Lyman $`\alpha `$is well fitted in the FOS spectrum by a column density of neutral hydrogen of log$`N`$(H i )$`=20.45\pm 0.15`$ at $`z=0.6819`$ (Fig. 1). A small correction for scattered light ($`1.5`$ % of the continuum) was applied to bring the core of the absorption line to zero flux.

We detect absorption at the expected position of numerous metal lines associated with this system. Table 1 lists the line parameters resulting from the fits, and Fig. 2 and 3 show some of the absorption features in velocity space.

We adopt $`3\sigma `$ upper limits for lines which were not detected.

Comparable $`b`$ values resulting from fits of individual components for different ions suggest that turbulent broadening is dominant. Under this assumption, for the cases where blending complicated the fitting, a typical $`b`$ value was adopted as derived from the high resolution spectra for resolved lines. Keck spectra were filtered at the FOS instrumental resolution and new fits of the Fe ii transitions longwards of Ly$`\alpha `$ were performed. We obtain $`b`$ values of the order of $`30`$ km s-1 by fitting the whole complex as a unique absorption feature at the resolution of the G270H spectra, and $`b25`$ and 16 km s-1 fitting with two components at the resolution of G190H. Because of the severe blending in the Ly$`\alpha `$ forest, the line fits performed at the limited FOS resolution should be considered upper limits. However, the study on metal abundances will be based on column densities derived from the high resolution spectra.

The numerous Fe ii transitions observed in the Keck spectra allow a reliable estimation of N(Fe) by fitting the strongest components in non-saturated transitions (Fe ii 2260,2374). The best fit for Fe ii  requires at least 14 components (13 in the main feature plus the satellite component at $`v100`$ km s-1 ).

The non-detection of Zn ii 2026 indicates that logN(Zn) $`<11.76`$. Several tests that modified the continuum level locally by means of fitting adjacent Ly$`\alpha `$ lines did not alter significantly the quoted value for the upper limit. Cr ii transitions are blended with Ly$`\alpha `$ lines and with C iv at $`z=1.234`$, so that an estimate of the column density for Cr ii could not be derived. For the absorption feature between $`v100`$ and 170 km s-1 , Cr ii 2062 yields an upper limit of logN(CrII)$``$ 12.2 while logN(ZnII)$``$ 11.8. Ti ii 1910.6,1910.9 are in a blended and noisy region of the spectrum, but recent Keck spectra (kindly provided by A. Songaila) covering up to 6060 Å allowed the detection of Ti ii 3384 (at $`4\sigma `$ level) at $`z=0.6818927`$ with logN(Ti ii )$`=11.79\pm 0.28`$, $`b=1.22\pm 0.77`$.

### 3.2 DLA at $`z=1.15`$ towards HE 0515$``$4414

According to IUE data, the extremely bright QSO HE 0515$``$4414 ($`z`$=1.7, V=14.9, Reimers et al. 1998) has a DLA at $`z`$=1.15 with an estimated column density for neutral hydrogen of logN(H i ) $`20.4`$. We observe absorption associated with this system at the expected positions of Al ii 1670, Al iii 1854,1862, C i 1656, C iv 1548,1550, Ca ii 3934, Fe ii 1608,2344,2374,2382,2586,2600, Mg i 1827,2026,2852, Mg ii 2796,2803, Mn ii 2576,2594,2606, Ni ii 1741, Si ii 1526,1808, Zn ii 2026,2062. Upper limits for the column density of ions not detected over $`3\sigma `$ are determined for Cr ii 2056,2062, Fe i 2484, Na i 3303, S i 1807 and Ti ii 3242. Figs. 4 and 5 show some of these ions in velocity space, and Table 3 the line parameters resulting from the fits. The line profiles of transitions of Mg ii , Fe ii and Al ii reveal that this remarkably complex system extends over $`700`$ km s-1 . At least three substructures can be distinguished: from approximately $`575`$ to $`270`$ km s-1 , from $`260`$ to $`100`$ km s-1 and the main absorption feature from $`80`$ to 140 km s-1 . Table 4 shows the derived column densities and metal abundances. The number of components is possibly underestimated because of the blending and blanketing problems. The $`b`$ values resulting from the fit of a blend cannot be well determined, but the column density of a blend should be a good approximation of the column densities of individual lines, when the system is not heavily saturated (Jenkins 1986). The best fit for Mg ii and Fe ii requires at least 20 components for the whole system.

The main features in the Fe ii  complex were fitted on the weakest transitions (unsaturated), and the secondary ones on the strongest transitions (Fe ii 2382,2600 were used for the absorption from $`575`$ to $`270`$ km s-1 and Fe ii 2382,2600,2344 for the range $`260`$ to $`100`$ km s-1 ). Al ii 1670 is saturated in the strongest substructure, and blended with Mg ii at $`z=0.2813`$ in the secondary one. Al iii 1862 is blended with Mg ii 2803 at $`z=0.4291`$.

A marginal $`2\sigma `$ Ni ii detection was obtained summing the spectral regions of Ni ii 1741, Ni ii 1709 and Ni ii 1751 to improve the signal to noise ratio (stacking technique, e.g. Lu 1991). The upper limit for N(Cr) was derived in the same way.

### 3.3 Mg ii System at $`z=0.28`$ towards HE 0515$``$4414

The strong Mg ii 2796 absorption features at $`z=0.281`$ reveal the presence of a multicomponent system spanning over $``$200 km s-1 from $`40`$ to 175 km s-1 , if $`v=0`$ km s-1 corresponds to $`z=0.28135`$ (Fig. 6). The corresponding Mg ii 2803 lines could only be seen for the main structure (from $`50`$ to 40 km s-1 ), since the rest is blended with Al ii 1670 at $`z=1.15`$. The low redshift does not allow to identify many associated features in our optical data. Absorption lines at the expected positions of Fe ii 2600, 2586, Fe i 3021, Mn ii 2576,2606, Ca ii 3969 and Al i 3803 have been looked for and the results of the fits are shown in Table 6. We obtain a Mn ii /Fe ii ratio about 15 times larger than solar, an improbable high value for a possible DLA. This suggests that the ionization is more complex, and/or that the continuum in the noisy region where Fe ii is expected could be underestimated. Several tests considering higher local continua yield an upper limit of N(Fe)$`13.7`$ which would mean an abundance of Mn relative to Fe five times larger than solar, if the Mn ii 2576 feature is real. It remains unclear whether this system is a new DLA candidate at very low redshift, maybe with an important part of Fe ii depleted onto dust. Based on N(Mn ii ) and assuming a Mn abundance relative to H of $`1/10`$ solar, we obtain log$`N`$(H i )$`19.9`$.

## 4 Total Metal Abundances

In Table 2, 4 and 5 we show the total metal abundances derived assuming that H i dominates and, therefore, that ionization corrections are not important, as shown by ionization models (Lu et al. 1995, Viegas 1995).

Both DLA are metal-poor with values for the metallicity of 1/10 to 1/20 solar. These values are comparable to the abundances found in other DLA and confirm the results from previous studies (Lu et al. 1996, Pettini et al. 1997a,b, 1999a,b and references therein.):
(i) The large scatter in the metallicity of DLA at any redshift (1.6 in \[Zn/H\] over the whole range of $`z`$) is consistent with the heterogeneous types of galaxies seen in the cases where the associated galaxy could be identified (e.g. Le Brun et al. 1997, Steidel et al. 1994, 1995, Rao & Turnshek 1998). Different star formation rates in spiral galaxies (Lindner et al. 1999), different impact parameters, and possible local inhomogeneities, would contribute to the scatter.
(ii) The enlarged sample shows no metallicity increase with decreasing $`z`$, contrary to expectations from models of cosmic chemical enrichment (e.g. Madau et al. 1996). Including the new values in the sample compiled from the literature, we obtain a column density-weighted mean abundance of Zn (Pettini et al. 1997a) of $``$\[Zn/H\]$`1.05`$ at $`z=1`$, or $``$\[Zn/H\]$`_{\mathrm{corr}}0.9`$ if Zn is weakly depleted onto dust as in the Galaxy.

A bias in the observations against dust-rich absorbers could mean that the DLA found at low and moderate $`z`$, where the dust presence should be more important, are not representative galaxies (e.g. Boissé et al. 1998 and ref. therein). However, even considerable amounts of dust would not have been able to obscure the line of sight of the bright HE 0515$``$4414 , and the intervening DLA still shows a rather low metallicity \[Zn/H\]$`1`$. The upper limit of \[Zn/H\]$`<1.3`$ at $`z`$=0.68, obtained at the Keck resolution, is also noteworthy.

## 5 Relative Metal Abundances

The ratio of two elements measured in DLA can vary from the solar value due to a different intrinsic pattern and/or to the condensation of refractory elements into dust depletion. We expect that the DLA might have abundance patterns resembling those in Halo stars, which show an overabundance of $`\alpha `$-group elements relative to Fe, and the odd-even effect (e.g. Mn underabundant relative to Fe.)

Assuming that elements in DLA are depleted as in the Galactic interstellar medium, Vladilo (1998) corrected the observed relative abundances from dust effects. He found the 17 systems of his sample chemically homogeneous over a range of metallicities from nearly solar down to $`1/32`$ solar and the dust-corrected abundances remarkably close to the solar values. Applying the same correction to our systems, we obtain the dust-corrected abundances shown in Table 7.

Due to the probable heterogeneity of the DLA population, it is interesting to discuss the abundance pattern in the systems individually.

### 5.1 DLA at $`z=1.15`$ towards HE 0515$``$4414: evidence for dust depletion

The refractory elements Cr, Ti, Ni are strongly depleted in this system: \[Cr/Zn\]$`<1.04`$, \[Ni/Zn\]$`<1.19`$, \[Fe,Ti/Zn\]$`<0.66`$. These abundances are comparable with those observed in the warm disk and halo of our Galaxy (Savage & Sembach 1996) . The corrected ratios are, in general, in line with those in Vladilo’s sample and consistent with solar. In particular, $`\alpha `$ elements (e.g. Si, Ti) are not found to be overabundant relative to Fe, and Mn is not significantly underabundant relative to Fe-peak elements. 111The two Mn/Fe measurements from Vladilo (1998) yield \[Mn/Fe\]$`{}_{\mathrm{corrected}}{}^{}0.2`$ when adopting log(Mn/H)$`{}_{}{}^{}=6.47`$, from Anders & Grevese (1989).

The observed ratio N(Ca ii )/N(Mg i ) in this DLA, like in two DLA at $`z=0.558,0.5240`$ (Vladilo et al. 1997 and references therein), is consistent with the predictions for warm gas with the base depletion representative of low density lines of sight in our Galaxy (Jenkins 1987). A selection bias would favour lines of sight with a base depletion pattern rather than a dense depletion, which would be present in gas with a smaller cross-section.

Recent observations of diffuse ISM along several lines of sight (Meyer et al. 1998 and ref. therein) show evidence for an ISM abundance of 2/3 solar, rather than solar, mainly based on measurements of O and Kr but also consistent with N, Zn, S and maybe C. We revised the correction for dust depletion in the abundance pattern of HE 0515$``$4414 , assuming intrinsic abundances of 2/3 solar for all elements. The pattern remains almost unchanged, with deviations lower than 0.05 dex (mostly of $``$ 0.02 dex.)

The degree to which an element is depleted by dust grains in the Galactic ISM is known to anticorrelate roughly with the temperature at which half the atoms in the gas phase condense into solid form (e.g. Jenkins 1987). A trend suggesting a depletion pattern can be recognized for this DLA (Fig. 7), but it is mainly based on Zn, the only non-refractory element for which abundances could be determined.

### 5.2 DLA at $`z=0.68`$ towards HE 1122$``$1649: no depletion pattern

The relative abundances observed in the DLA at $`z=0.68`$ (Table 2) differ considerably from solar, but follow a different pattern than the system at $`z=1.15`$:

* The comparable Fe, Ti and Zn content in this system indicates that the effect of depletion onto dust is negligible, as expected at this low metallicity. Corrections for possible dust depletion are negligible.
* Mn is underabundant relative to Fe, as in the Galactic Halo, and \[Al/Si\]$`=0.15\pm 0.15`$ is also compatible with the odd-even effect.
* The overabundance of $`\alpha `$-elements relative to Fe is observed for Si/Fe and compatible with the upper limit of O/Fe; however, \[Ti/Fe\]$`=0.3\pm 0.4`$.

Therefore, the abundance pattern shows similarities with metal-poor Halo stars. The low depletion degree could indicate that grain destruction is more efficient than in the Galaxy, maybe due to increased interstellar shocks from supernova explosions, as expected in a galaxy with intense star formation. Further abundance measurements of other weakly depleted elements (like S and O) would help to discern the intrinsic abundance pattern in DLA galaxies.

Plotting the measured abundances versus condensation temperature is consistent with the absence of a depletion pattern (Fig. 7). However, as argued by Pettini et al. (1997a), this lack of correlation should not be overinterpreted at such a low degree of chemical enrichment.

Other lines of sight in the enlarged $`z<1.5`$ DLA sample where the observed abundances could be corrected for depletion effects (Table 7), also show diverse patterns: e.g., low dust content is found at $`z=0.8598`$, $`z=1.3726`$, and a pattern compatible with a solar one (or 2/3 of solar) for the corrected abundances at $`z=1.223,1.149,1.42,0.692`$. The different behaviours confirm that DLA arise in a heterogeneous group, as also inferred from the study of total abundances.

### 5.3 DLA Abundances and the Nucleosynthesis of Nitrogen

The Nitrogen abundance help us to decipher the Galactic chemical enrichment history. At low metallicities, N primary production, which takes place during the asymptotic giant branch phase of intermediate mass stars, is expected to dominate the secondary N production (e.g. Lu et al. 1998 and references therein). The time-delay between the production of primary N from intermediate mass stars and elements like O, released from massive stars, would imply a large scatter in N/O at low metal abundances (e.g. Marconi et al. 1994). However, observations of nearby metal-poor galaxies do not unambiguously show this pattern.

DLA offer the possibility of studying N abundances for metal-poor systems tracing early stages of the galactic chemical evolution and, therefore, of clarifying the N nucleosynthesis. DLA for which N abundances –or at least reliable upper limits– are available, are at high $`z`$ (Pettini et al. 1995, Lu et al. 1998, Centurión et al. 1998 and references therein). The system at $`z=0.68`$ allows the first study of the N content at low redshift. In spite of the numerous N i transitions expected in our spectral range, only upper limits could be derived for N(N i ) due to the severe blending in the Ly$`\alpha `$ forest, yielding \[N/H\]$`<2,`$\[N/Si\]$`<1.36`$ and \[N/Fe\]$`<0.64`$.

To compare with other systems, we assume that O and the $`\alpha `$ elements (like S, Si) trace each other in DLA, as is the case in the Galaxy and in extragalactic H ii regions. On the other hand, the lack of depletion onto galactic-like dust grains in this system implies that Fe and Si measurements correspond to the total elemental values, and therefore \[N/O\]$``$ \[N/S\] $``$ \[N/Si\] is assumed.

Lu et al. (1998) and Centurión et al. (1998) found that the N/Si, N/Fe and N/S ratios in a sample of DLA at high $`z`$ show considerable scatter at low metallicities, as predicted by the standard model of delayed N primary production. The observed N/Fe for our DLA at $`z=0.68`$ is compatible with the range of measurements in Halo stars (cf. Figs. 6 and 7 in Centurión et al.). The low N/Si ($``$ N/O) limit agrees with the large dispersion and hints at a galaxy with recent star formation: O and Si have recently been released from massive stars and are overabundant relative to N, which has not had time to be primarily delivered from intermediate mass stars in significant amount and has not been secondarily produced in an efficient way at such a low metallicity level. The $`\alpha `$/Fe enhancement expected in a galaxy with a recent starburst is also observed in this system. On the contrary, Centurión et al. found no $`\alpha `$ overabundance in the few DLA with low N/O values, and suggested that alternative models are required to explain those abundance patterns (e.g. primary N production in massive stars).

## 6 Neutral Carbon

Neutral species should be detectable in gas with a large H i column density. The detection of C i associated with DLA would be particularly interesting, since the ground-state fine-structure excitation depends on the cosmic microwave background radiation, $`T_{CMB}`$, and measuring the relative ground-state populations one can estimate the excitation temperature and an upper limit for $`T_{CMB}`$ at that epoch (e.g. Ge et al. 1997 and references therein).

In Table 8 we compare the N(C i ) measured in our systems with other values from the literature for a large range of $`z`$, and with those observed in Galactic lines of sight towards $`\zeta `$ Oph (diffuses interstellar medium representative of the Galactic disk) and $`\zeta `$ Pup (more tenuous intercloud medium and H ii region). For the DLA at $`z=1.15`$, the absorption features at the expected wavelength of C i 1656,1560 yield log$`N`$(C i )$`=13.68\pm 0.09`$ and log$`N`$(C i )$`=13.44\pm 0.07`$. The relative population ratios imply an excitation temperature of 14.2 K, which is an upper limit to the $`T_{CMB}`$, since other mechanisms may contribute to the C i excitation (Lu et al. 1996). On the other hand, the complex structure of the $`z=1.15`$ DLA and the considerable noise in the C i 1560 region makes it difficult to deblend C i  from further C i components. For the system at $`z=0.68`$, we obtained an upper limit of log$`N`$(C i )$`<13.45`$.

Comparing these values with the C i /H i ratios known in Galactic gas (Jenkins & Shaya 1979) we find that the $`z=1.15`$ system reproduces the conditions of the Galaxy, as also found for a DLA at $`z=0.8597`$ towards PKS 0454+039 (Boissé et al. 1998) and $`z=1.97`$ towards Q 0013-004 (Ge et al. 1997). The upper limit for the system at $`z=0.68`$ is still compatible with the conditions expected for Galactic gas but could also hint at a lower C i content as observed for a $`z=0.692`$ DLA towards 3C 286 (Boissé et al. 1998) and for the DLA at higher $`z`$ in Table 8.

Plotting the C i /H i  ratio versus redshift (Fig. 8), a trend of lower C i /H i  with increasing $`z`$ can be recognized. The overplotted linear regression is computed including upper limits as exact values. Considering that most of high $`z`$ measurements are actually upper limits, the observed trend would be more pronounced. The evolution of the metagalactic ionizing background (more intense at high $`z`$ than at $`z\mathrm{\Gamma }<1`$, e.g. Haardt & Madau 1996) is probably the main reason of the C i /H i evolution with $`z`$. On the other hand, the local presence of dust could also be playing a role, since dust grains could shield more efficiently the gas from being photoionized. In fact, the DLA in Table 8 with a low dust content show a low C i /H i ratio, regardless of the redshift.

## 7 Molecules in DLA

We do not see (H2) or CO absorption, consistently with the low metal abundances.

Molecular hydrogen (H2) has numerous transitions at $`\lambda <1100`$ Å and carbon monoxide (CO) at $`\lambda <1550`$ Å . Because of the blending problem it is better to look for absorption-free windows where strong H2 is expected than looking for strong lines themselves (Ge & Bechtold 1997). We do not find absorption at the expected positions of some strong H2 lines (e.g. around $`\lambda _{rest}=1094,1079`$) associated with the DLA at $`z=0.68`$. The H2/H ratio is small compared to Galactic interstellar clouds, a common result for DLA (Ge & Bechtold, 1997). For the DLA at $`z=1.15`$, the spectra do not cover the region where H2 absorption is expected.

We searched for strong CO transitions ($`\lambda _{rest}`$ = 1477.568, 1477.355, 1509.75, 1419.046, 1392.525, 1544.451 Å) and obtained $`2\sigma `$ upper limits of logN(CO)$`<14.4`$ for the DLA at $`z=0.68`$ and logN(CO)$`<13.8`$ at $`z=1.15`$ (where only $`\lambda _{rest}`$ = 1544 Å falls in the available wavelength range). Contrasting these values with the upper limits derived by Lu et al. (1998) at the Keck HIRES resolution for several DLA at high $`z`$, we find that our CO limits are roughly comparable to the values deduced for diffuse interstellar clouds in the Milky Way.

## 8 Kinematics: Complex Structure at Moderate $`𝒛`$

The DLA at $`z=0.68`$ towards HE 1122$``$1649 extends over $`230`$ km s-1 in velocity space plus a satellite component at $`100`$ km s-1 . As we can see in Figs. 2 and 3, the profiles of low ions observed at high resolution (cf. Fe ii , Mg ii , Mn ii ) track closely one another and show an edge-leading asymmetry, especially if we consider the feature from $`20`$ to 100 km s-1 .

Remarkably, the Mg i 2852 profile shows a different asymmetry: the main feature lies between 100 and 150 km s-1 , while Mg ii – along with other singly ionized ions – shows the strongest absorption between $`20`$ and 100 km s-1 . This suggests that the region between 100 and 150 km s-1 is less ionized, since the feature could not be identified as a metal line associated with another absorber. The difference between Mg i  and Mg ii  profiles could hint at a non-negligible contribution of an intervening H ii  region for this system. The ionized gas would produce an extra contribution to the column densities of species with an ionization potential larger than 13.6 eV, so that the measured abundances of ions like Fe ii , Mn ii , Si ii and Zn ii  would correspond to both regions. The metal abundances associated with the H i  system would therefore be lower.

The low-ion profiles of the system at $`z=1.15`$ towards HE 0515$``$4414 (Figs. 4 and 5) span more than 700 km s-1 , which is the largest extent in the velocity space seen for a DLA. As pointed out before, three main structures can be recognized: from approximately $`575`$ to $`270`$ km s-1 , from $`260`$ to $`100`$ km s-1 and the main absorber from $`80`$ to 140 km s-1 . This main complex shows no ordered profile and is consistent with absorber components distributed randomly in $`v`$.

The remarkably complex structure of low and intermediate $`z`$ DLA provide new hints about their nature. Two kinds of models have been discussed to explain kinematic characteristics of DLA at high $`z`$: (1) thick rotating disks evolving without merging to form present-day galaxies (Prochaska & Wolfe 1997a,b,1998), and (2) merging protogalactic clumps in the standard Cold Dark Matter cosmology (Haehnelt et al. 1998, McDonald & Miralda-Escudé submitted). The model of rotating thick galactic disks can explain the edge-leading asymmetry of profile structures extending up to $`250`$ km s-1 . An asymmetric structure where the absorption decreases monotonously is found in the system at $`z=0.68`$, but the one at $`z=1.5`$ shows a non-ordered profile. The large extent of the systems suggests that they cannot be single rotating disks.

Alternatively, several authors (e.g. Matteuci et al. 1997, Vladilo et al. 1997) have suggested that dwarf galaxies rather than spiral galaxies could be responsible for DLA. The extent of the kinematic profiles of our two new DLA are again too large to be explained by individual dwarf galaxies. On the other hand, Khersonsky & Turnshek (1996) showed that DLA can be associated with neutral gas in giant hydrogenous clouds that could be associated with any type of galaxy or protogalaxy.

Considering only the deepest part of the absorption ($`\tau 0.7`$) most profiles in a sample of high $`z`$ DLA studied by Ledoux et al. (1998) show $`\mathrm{\Delta }v<150`$ km s-1 and a trend of decreasing velocity with increasing redshift. Under the same assumptions, our DLA at $`z=0.68`$ and $`1.15`$ match fairly well this behaviour, which is expected in CDM structure formation scenario (e.g. Kauffmann 1996), but not if DLA are disks which do not evolve kinematically. Similarly to the models described by Haehnelt et al. 1998 for high $`z`$, the velocity broadening in the clumps at intermediate $`z`$ would be due to a mixture of random motions, infall, merging and also rotation, which would still be important at smaller scales and could be responsible for the edge-leading asymmetry observed in some of the clumps.

Wolfe & Prochaska (1998) found a correlation between kinematics and metallicity for a sample of 17 DLA. At $`z=1.63`$ they detected no systems like our two, with a large extent in the velocity space ($`\mathrm{\Delta }v>120`$ km s-1 ) and low metallicity (\[Zn/H\] $`1`$) at the same time. On the other hand, the new systems do show moderate $`N`$(H i ) of $`10^{20.4}`$ cm-2, following the trend of larger $`\mathrm{\Delta }v`$ to lower $`N`$(H i ) observed for high redshift DLA (Wolfe & Prochaska 1998).

The ions associated with the DLA at $`z=1.15`$ towards HE 0515$``$4414 seem to track each other fairly well, but the lack of unsaturated lines from several ions at some velocities makes it difficult to establish the variation of the relative abundances. However, the strongest feature in the profiles of unsaturated Fe ii and Mn ii show two absorption features of comparable intensity, while for Zn ii and Si ii , the first one is clearly larger. Species weakly related to dust in the Galactic ISM seem to show a different profile than those corresponding to refractory elements. This could mean that there is more dust present in the components between $`50`$ and $`+50`$ km s-1 , possibly associated with a higher metallicity level, and this dust could be partly concealing the edge-leading asymmetry.

## 9 Associated Highly-Ionized Gas.

Damped Ly$`\alpha `$ systems studied at high $`z`$ often show considerable Si iv and C iv absorption with profiles resembling each other, but usually different from low ions (e.g. Lu et al. 1996). This suggests that highly-ionized gas is produced in distinct physical regions, which could be surrounding and shielding the low-ionized gas. In our two DLA at moderate redshift, we do not find strong lines from highly-ionized species covering a much larger extent than the weakly-ionized species.

* DLA at $`z=0.68`$ towards HE 1122$``$1649: Fig. 3 (right panel) shows absorption near the expected positions of C iv , O vi , Si iii and Si iv , none of which is unambiguously present.
* DLA at $`z=1.15`$ towards HE 0515$``$4414: Only weak C iv absorption with $`N`$(C iv ) $`=13.18\pm 0.12`$ was found associated with the main complex at $`z=1.15`$ (from $`80`$ to 140 km s-1 ), contrary to most of DLA, but like $`z=0.692`$ towards 3C286 (Boissé et al. 1998), where no associated high ions could be found except for weak C iv . This suggests that the bulk of gas is weakly ionized and shielded from the ionizing radiation field. Stronger absorption features for the C iv doublet yielding N(C iv )= 14.38 $`\pm 0.17`$ could be fitted in the noisy region extending from $`550`$ to $`200`$ km s-1 (Fig. 4), which is less prominent in the low-ion profiles.

Al iii  measured in the DLA follows the low ionization species (Al iii  is observed for the three main substructures at $`z=1.15`$), but it is present in a much lower amount than Al ii . These results confirm the very low ionization degree in these systems and suggest that Al iii  is probably produced in the same physical region as the low ionization species.

The absence of highly ionized gas associated with the large extents of low ionized regions suggests that the ionizing background is weaker than at high $`z`$, as expected from the cosmological expansion and from models of the ionizing background (e.g. Haardt & Madau 1996).

## 10 Summary

We have studied the properties of the DLA at $`z=0.68`$ towards HE 1122$``$1649 and at $`z=1.15`$ towards HE 0515$``$4414 . A further DLA candidate at low redshift ($`z=0.28`$) has been found in the line of sight of HE 0515$``$4414 . Very few DLA at $`z<1.1`$ had been analyzed up to now and enlarging the sample is essential to understand the nature of DLA, cosmic chemical evolution and possible selection effects which could be biasing our view of the Universe. The main results of our analysis, interpreted in comparison with DLA at other $`z`$, are summarized in this section:

* A large number of associated metal lines have been detected at high resolution. Metal abundances are found to be rather low (of 1/10 and $`<1/20`$ of solar value) and comparable to high-$`z`$ DLA. Therefore, the new values confirm the large scatter at any $`z`$, rather than a clear increase in the metallicity with decreasing $`z`$ for the DLA population. The new systems also confirm the heterogeneity of the DLA population (in agreement with the results from imaging studies, e.g. by Le Brun et al. 1997).
* The relative abundance pattern for the DLA at $`z=1.15`$ in the line of sight of the very bright HE 0515$``$4414 indicates a significant degree of depletion, with a dust-to-metal ratio comparable to the Galactic value. The intrinsic pattern, after correcting for dust effects, is comparable to solar.
* The system at $`z=0.68`$ towards HE 1122$``$1649 presents a negligible dust content and its abundance pattern resembles that in metal-poor Halo stars. Relative abundances suggest that most metals in the $`z=0.68`$ DLA where made recently. Primary N production in massive stars is not required to explain the abundance pattern.
* The search for neutral carbon shows that systems without significant amounts of dust are not as shielded from ionizing sources as systems where dust is present. C i /H i  increases with decreasing $`z`$, as expected from the evolution of the metagalactic ionizing background.
* The non-detection of H2 or CO molecules is congruent with the low metal abundances.
* The high-resolution data reveal a very complex structure for the low-ion profiles. The system at $`z=0.68`$ spans over $``$ 300 km s-1 and shows an edge-leading asymmetry. The $`z=1.15`$ system shows at least three main substructures over more than 700 km s-1 . This suggests that the line of sight intersects a group of galaxies or subgalactic units. The velocity broadening could be due to a mixture of random motions, infall, merging and rotation.
* Only weak absorption from highly-ionized species were detected and it probably originates in regions distinct from the low-ionized gas.

Further surveys will provide new DLA candidates, which should be observed with the new echelle spectrographs and, at moderate redshift, also with imaging techniques in order to detect their galactic counterparts. Enlarging the DLA sample and clarifying their nature will help us to understand the cosmic chemical evolution and to discern among models of galaxy formation.

(1) This work, (2) Boissé et al. 1998, and ref. therein, (3) Steidel et al. 1997, (4) Meyer & York 1992, (5) Lu et al. 1996, (6) Pettini & Bowen 1997, (7) Meyer, Lanzetta & Wolfe 1995, (8) Pettini et al. 1997a, and (9) Pettini et al. 1999.

(1)Meyer & York 1987a, (2)Meyer & York 1987b and ref. therein, (3) Prochaska & Wolfe 1996, (4)Meyer et al. 1986, Chaffee et al. 1986, (5) Ge et al. 1997, (6) Chaffee et al. 1988 and ref, therein, (7) Blades et al. 1985, (8) de la Varga 1999, (9) Boissé et al. 1998 (10) Morton 1978, (11) This paper.


format: text

text:
# The Pseudo-𝐶_𝑙 method: cosmic microwave background anisotropy power spectrum statistics for high precision cosmology

## I Introduction

During the next few years ground based observations and balloon missions  as well as satellite observation  promise exquisite determinations of the angular distribution of cosmic microwave background (CMB) anisotropies.

Inflationary cosmogonies, the most popular theories of structure formation, predict a statistically isotropic CMB with Gaussian fluctuations. To linear order in perturbation theory this means that within this class of models all cosmologically relevant information contained in the map of anisotropies is distilled in the angular power spectrum $`C_l`$.

An accurate estimate of these quantities therefore promises to have tremendous impact on our knowledge of the global properties of the universe and the physical processes which lead to the formation of structure. This knowledge is quantified in terms of cosmological parameters such as $`\mathrm{\Omega }_{tot}`$, $`\mathrm{\Omega }_b`$, $`H_0`$, etc.

The $`C_l`$ are uniquely defined as rotationally invariant quadratic combinations of the coefficients of an expansion of the anisotropy $`\frac{\mathrm{\Delta }T}{T}`$ in spherical harmonics. Their useful statistical properties rest on the the fact that the spherical harmonics are a complete and orthogonal basis set on the sphere. However, in realistic experiments, only part of the celestial sphere is covered, if only because our CMB sky is obscured by the Milky Way. Also, due the particulars of the observational strategy, non–uniform and possibly correlated noise contaminates the observed signal. In order to interpret past observations and forecast the impact of future observations we need a framework for studying the sampling statistics of the $`C_l`$ for incomplete sky coverage in the presence of non–uniform noise.

It is worth noting that there are several different kinds of $`C_l`$ which we will be referring to in this paper.

The stochastic nature of the observed $`C_l`$ is illustrated in Figure 1, where the straightforward computation of $`C_l`$ from a subset of our Monte Carlo simulations of incomplete and noisy skies are overlaid on the theoretical predictions. One or more of the following simplifying assumptions can be invoked in cases where one would like to avoid the full complexity of the problem (e.g.  ):

* (Independence) The observed $`C_l`$ are approximately independent due to nearly full and uniform sky coverage,
* (Scaled $`\chi ^2`$) Their sampling distributions do not change appreciably from the $`\chi ^2`$ distributions which apply for the full sky, apart from a rescaling of the number of degrees of freedom by the sky fraction $`f`$, and
* (Gaussianity) Their sampling distributions are well–approximated by Gaussian distributions which have the same first and second moments as these rescaled $`\chi ^2`$ distributions.

There are several reasons to carefully assess these approximations and, if necessary, go beyond them. One reason is that in the short term balloon and ground based experiments will provide the leading edge science results to the field. Due to practical limitations we cannot hope to even come close to full sky coverage with these types of experiments. The observation regions are often ring–shaped or cover a circular region of the sky which subtends a small solid angle. It is an urgent matter to assess statistically how these imminent experiments can constrain the power spectrum and cosmological parameters.

Further, the planned satellite missions are aiming to determine the $`C_l`$ to sub–percentage accuracy. In the case of the Planck Surveyor mission, this has given us the hope of detecting small effects such as secondary anisotropies which are due to nonlinear gravitational effects on the CMB photons during the free–streaming epoch. An example of such an effect is anisotropies due to gravitational lensing. This is an important issue because a detection would break otherwise present parameter degeneracies and allow a consistent parameter estimation from CMB data alone .

More generally, the impact CMB observations will have on cosmology makes it important to use approximations in a controlled way. For example, in the analysis of COBE–DMR data it was realized that using Gaussian approximations for quadratic quantities introduced systematic biases . At the same time, we need approximations to make feasible the analysis of the huge CMB data sets we expect in the coming years.

The Maximum Likelihood Inversion (MLI) method (cf.  for applications in a cosmological context) is the standard framework for solving this type of problem. To illustrate, we write down the likelihood functional for a Gaussian sky. In order to do so we need to define the vector of observed pixel temperatures on the sky $`𝐝`$ with $`N_{pix}`$ components. A given theory and observational strategy on the sky imply the covariance structure of the elements $`d_i`$ in terms of signal and noise covariance matrices $`𝐒`$ and $`𝐍`$ such that

$$d_id_j=S_{ij}+N_{ij}.$$
(1)
The signal covariance matrix $`S`$ is simply a discretised version of the theoretical angular correlation function

$$\xi \left(\alpha _{ij}\right)=\underset{l=2}{\overset{\mathrm{}}{}}\frac{2l+1}{4\pi }C_lP_l(\alpha _{ij}),$$
(2)
evaluated at the pixel locations. The noise covariance $`𝐍`$ reflects those imperfections of the experiment and observational strategy which lead to spatial correlations in the projected noise.

This leads to an expression for the likelihood

$$L(d|\mathrm{theory})=\frac{\mathrm{exp}\left(\frac{1}{2}𝐝^𝐓(𝐒+𝐍)^1𝐝\right)}{\sqrt{(2\pi )^{N_{pix}}\left|𝐒+𝐍\right|}}$$
(3)
The matrix $`𝐒`$ is a function of predicted $`C_l`$, which are in turn determined by the parameters of some theory. The matrix $`𝐍`$ includes a model of the expected noise covariance given the experimental setup. The structures of both $`𝐒`$ and $`𝐍`$ also depend on the observed region on the sky.

For a given set of data the likelihood is then considered as a function of the parameters of a given theory. The maximum likelihood parameter estimates are determined as those which maximise Eq. (3). While these manipulations are conceptually simple, the practical difficulties are apparent. The number of operations needed to evaluate Eq. (3) (or equivalently an estimator derived from it) scales as the most costly operations contained in it. These are the determinant evaluation in the denominator and the matrix inverse in the exponent. For a general situation both of these require $`N_{pix}^3`$ operations.

In the near future $`N_{pix}`$ will be $`10^5`$ from balloon experiments (such as Maxima), then $`10^6`$ for MAP and in less than a decade $`10^7`$$`10^8`$ for Planck. These numbers imply $`10^{15}`$$`10^{24}`$ computations. Even when considering future growth in computing performance, estimates of the computational time required to evaluate Eq. (3) a single time reach up to many years . To compute the maximum likelihood estimator requires many likelihood evaluations and would thus require an unfeasibly long period of time. For practical purposes, straightforward applications of the MLI method therefore fail to give answers in finite computational time.

Previous studies of this subject have usually directly addressed the problem of ”solving for the $`C_l`$”. Tegmark  rejects the maximum likelihood estimator on the grounds that it is too time–consuming to compute and suggests a minimum variance (and hence “optimal”) method based on a quadratic estimator which makes a form of the Gaussianity assumption above. While this method is powerful in the regime where it is applicable, it does not provide a means of assessing the validity of this assumption. Within our analysis we can make the connection between the minimum variance and the maximum likelihood methods (cf. subsection III A).

The approach of Bond, Jaffe and Knox  is again very different from ours. It consists in finding an approximate form for the likelihood and parametrising the shape of the power spectrum in terms of a set of top hat functions, whose height gives the average value of the $`C_l`$ over its width. The width of these top hats is chosen wide enough that correlations between them become small. Their method is of practical significance because their approximate ansatz for the likelihood is easy to evaluate.

Our perspective is the following: we endeavour to provide an ab initio formalism which lays bare the statistical properties of the power spectrum coefficients. It then turns out that the insight gained from this study can be used to formulate unbiased and well–controlled methods for parameter estimation.

To provide a quantitative basis for this discussion we present in this paper an efficient semi–analytic calculational framework within which the statistical properties of the $`C_l`$ power spectrum coefficients can be studied for theories which generate a Gaussian CMB sky. This framework is exact for the important class of survey geometries and noise patterns which obey rotational symmetry about one arbitrary axis and is not limited to large sky coverage. This class contains both skies with symmetric and asymmetric cuts north and south of the Galactic equator, rings, annuli, and polar cap shaped regions of any size as indicated in Figure 2. For this class we obtain exact analytical solutions for the mean $`C_l`$, the variance $`\left(\mathrm{\Delta }C_l\right)^2`$, the skewness, kurtosis and all higher moments of the sampling distribution for all $`l`$. We derive exact and conveniently computable analytic expressions for the marginalised probability distributions of the cut sky $`C_l`$ for all $`l`$ as well as integral representations for the joint distribution of all $`C_l`$. As an added bonus we find that our methods also allow dealing with arbitrary noise patterns to very good accuracy and conclude by illustrating their use in some applications.

To test our results we perform state–of–the–art Monte Carlo simulations for a high resolution CMB satellite with a Galactic cut of size $`\pm 20^{}`$ (resulting in a sky coverage comparable to COBE). We use a sine modulated noise template as a simple model of the noise pattern which would result from scanning along great circles through the ecliptic poles. The numerical results we obtain provide a solid check on our analytic calculations.

The plan of this paper is as follows: in section II we lay out the general framework. Section III develops the statistical theory of power spectrum coefficients for observation regions which are azimuthally symmetric about some axis. We arrive at expressions for the cut–sky moments of any order and solve the moment problem to obtain the exact sampling distributions. We compare our analytical results to Monte–Carlo simulations in section IV. Several applications of our methods are presented in section V. Section VI contains further discussion and our conclusions. Appendix A discusses the effects of non-cosmological foregrounds on the results presented in this paper. Details of our analytical and numerical techniques are given in three further appendices.

## II The framework

The full sky of CMB temperature fluctuations can be expanded in spherical harmonics, $`Y_{lm}`$, as

$$T(\gamma )=\underset{l=0}{\overset{l_{max}}{}}\underset{m}{}a_{lm}Y_{lm}(\gamma )$$
(4)
where $`\gamma `$ denotes a unit vector pointing at polar angle $`\theta `$ and azimuth $`\varphi `$. Here we have assumed that there is insignificant signal power in modes with $`l>l_{max}`$ and introduce the convenient notation that all sums over $`m`$ run from $`l_{max}`$ to $`l_{max}`$ but all quantities with index $`lm`$ vanish for $`m>l`$.

In a universe which is globally isotropic and contains Gaussian fluctuations, the $`a_{lm}`$ are independent, Gaussian distributed with zero mean and specified variance $`C_l^{theory}\left|a_{lm}\right|^2`$. Hence, for noiseless, full sky measurements, each measured $`C_l`$ independently follows a $`\chi ^2`$–distribution with $`2l+1`$ degrees of freedom and mean $`C^{theory}`$. Owing to Galactic foregrounds, limited surveying time or other constraints inherent in the experimental setup, the temperature map that comes out of an actual measurement will be incomplete. In addition, a given scanning strategy will produce a noise template: high noise per pixel results in regions where the scanning lines are less densely packed compared to regions where they are denser. We model the noise as a Gaussian field with zero mean which is independent from pixel to pixel and modulated by a spatially varying root mean squared amplitude $`W_N(\gamma )W(\gamma )`$ (the factor of $`W`$ appearing because the noise too is zero in unobserved regions). Therefore the observed temperature anisotropy map is in fact

$$\stackrel{~}{T}(\gamma )=W(\gamma )\left[T(\gamma )+W_N(\gamma )T_N(\gamma )\right]$$
(5)
where $`W`$ is unity in the observed region and zero elsewhere.

The coefficients in the spherical harmonic expansion recovered from such an incomplete observing region are

$`\stackrel{~}{a}_{l^{}m^{}}`$ $`=`$ $`{\displaystyle _𝒪}𝑑\mathrm{\Omega }Y_{l^{}m^{}}^{}(\gamma )T(\gamma )`$ (6)
$`=`$ $`{\displaystyle \underset{lm}{}}a_{lm}{\displaystyle _𝒪}𝑑\mathrm{\Omega }Y_{l^{}m^{}}^{}(\gamma )Y_{lm}(\gamma )`$ (7)
The notation “$`_𝒪`$” denotes integration over the observed region. Note that the usual orthogonality property of the $`Y_{lm}(\gamma )`$ does not hold any longer because we do not integrate over all solid angle. This becomes clearer if we define the geometric coupling matrix

$$W_{l^{}m^{}lm}_𝒪𝑑\mathrm{\Omega }Y_{l^{}m^{}}^{}(\gamma )Y_{lm}(\gamma )$$
(8)
and write

$$\stackrel{~}{a}_{l^{}m^{}}=\underset{lm}{}W_{l^{}m^{}lm}a_{lm}.$$
(9)
By definition, the $`W_{l^{}m^{}lm}`$ are just the matrix elements of the $`W(\gamma )`$ window in a spherical harmonic basis. These matrices have been previously discussed in a similar context in the pioneering papers by Peebles and Hauser .

This means that expanding $`\stackrel{~}{T}`$ as in Eq. (7) produces a set of correlated Gaussian variates $`\stackrel{~}{a}_{lm}`$ for the signal and $`\stackrel{~}{a}_{Nlm}`$ for the noise. These combine into power spectrum coefficients

$$\stackrel{~}{C}_l=\frac{1}{2l+1}\underset{m}{}\left|\stackrel{~}{a}_{lm}+\stackrel{~}{a}_{Nlm}\right|^2$$
(10)
whose statistical properties differ from the ones of the $`C_l`$. We therefore refer to these quantities as pseudo-$`C_l`$. In what follows we will discuss the statistical properties of these quantities.

## III Exact pseudo–$`C_ł`$ statistics

Let us focus on the terms of the sum Eq. (10). They are squares of Gaussian distributed variates with zero mean. It follows that the $`\stackrel{~}{C}_l`$ are sums of $`2l+1`$ $`\chi ^2`$ variates with one degree of freedom. However, each term in the sum has a different expectation value and may in general be correlated, so the $`\stackrel{~}{C}_l`$ are not $`\chi ^2`$ distributed with 2l+1 degrees of freedom.

If we assume white noise, $`C_l^NC^N`$, then each term in the sum has expectation value

$$\sigma _{lm}^2=\frac{_{l^{}m^{}}C_l^{}\left|W_{l^{}m^{}lm}\right|^2+C^N_O𝑑\gamma W_N(\gamma )^2\lambda _{lm}^2(\theta )}{2l+1}$$
(11)
A few remarks are in order:

1. Azimuthal symmetry has not been invoked to derive this relation.
2. The cross–term which comes from expanding the square in Eq. (10) has vanishing expectation value because signal and noise are assumed to be uncorrelated.
3. The sum over $`C_l`$ includes the non–cosmological $`l=0`$ and $`l=1`$ terms. The appearance of these terms merits some further discussion, which we relegate to Appendix A.

### A Digression: connecting to the least squares estimator

Summing Eq. (11) over $`m`$ gives a set of linear equations which also results in a different context, namely after minimising the variance of the $`C_l`$ estimator in Tegmark . There it appears as a set of linear equations for the estimator, which the author assumes to be Gaussian distributed.

We see from our treatment that this equation is only valid after the expectation values are taken. The quantities involved are quadratic in Gaussian variates which leads to the presence of a cross term which only disappears after taking the expectation value.

This is the reason why some $`C_l`$ estimates in  turn out negative for values of $`l`$ where the signal to noise ratio becomes of order unity. This is of course unphysical because the $`C_l`$ are properly quadratic and hence positive semi–definite quantities. We observe that it is precisely neglecting this fact which leads to unphysical results.

On the other hand, our treatment provides alternative proof of the fact that the least–squares estimator is unbiased. It also explains the significance of the Fisher matrix which appears in this context: is simply a measure of the geometric couplings between eigenmodes of total angular momentum.

### B Azimuthal symmetry

To elucidate the correlation structure between the $`\stackrel{~}{a}_{lm}`$, Eq. (9), we write Eq. (10) as a quadratic form

$$\stackrel{~}{C}_l=\frac{1}{2l+1}\underset{m}{}\underset{l^{},m^{},l^{\prime \prime },m^{\prime \prime }}{}\xi _{l^{}m^{}}^{}_{l^{}m^{}l^{\prime \prime }m^{\prime \prime }}^{(lm)}\xi _{l^{\prime \prime }m^{\prime \prime }}$$
(12)
in independent normal variates $`\xi _{lm}`$ with zero mean and unit variance. The object $``$ describes the combined effect of statistical correlation and the geometrical couplings between the various expansion coefficients.

For an isotropic Gaussian field on the full sky

$$_{l_1m_1l_2m_2}^{(lm)}=\delta _{ll_1}\delta _{mm_1}(C_l^{theory}+C_l^{noise})\delta _{ll_2}\delta _{mm_2}.\text{(isotropic, full sky)}$$
(13)
We show in Appendix B that for a Gaussian field with an anisotropic but azimuthally symmetric weight function $`W(\theta )`$, $``$ is of the form

$$_{l_1m_1l_2m_2}^{(lm)}=\delta _{mm_1}W_{l_1m_1}^{(l)}\left(C_l^{theory}+C_l^{noise}\right)W_{l_2m_2}^{(l)}\delta _{mm_2},\text{(anisotropic, cut sky)}$$
(14)
where

$$𝒲_{l^{}m}^{(l)}W_{l^{}mlm}=𝑑\mu W(\mu )\lambda _{lm}(\mu )\lambda _{l^{}m}(\mu ),$$
(15)
and we assume identical windows for signal and noise for notational simplicity in Eq. (14). The generalisation for arbitrary azimuthally symmetric noise windows is straightforward.

The structure exhibited in Eqs. (14) and (15) is the key fact which allows the derivation and cheap evaluation of the exact results we obtain. An efficient algorithm for the evaluation of the matrix elements of $`W`$ is given in appendix C. For a maximum $`l`$ of 1024 this takes just over one minute on a single R10000 CPU.

It follows that while the $`\stackrel{~}{a}_{lm}+\stackrel{~}{a}_{Nlm}`$ terms in Eq. (10) are correlated for different $`l`$, they will be independent for different $`m`$. As a consequence, under the assumption of azimuthal symmetry the $`\stackrel{~}{C}_l`$ are sums of independent Chi–squared ($`\chi ^2`$) variates, each with one degree of freedom but different expectations. Below we succeed in finding analytical expressions for the probability density functions of such generalised $`\chi ^2`$ variates and we hence solve for the statistical properties of the pseudo–$`C_l`$ exactly.

### C Pseudo–$`C_ł`$ moments

The problem is to compute the statistical properties of a sum of independent variates, each of which has a known distribution. A useful tool for attacking this kind of problem is the method of characteristic functions . Appendix D provides a heuristic introduction to this method.

The characteristic function of a $`\chi ^2`$ variate with one degree of freedom is given by

$$\frac{1}{(12ix\sigma _{lm}^2)^{\frac{1}{2}}}.$$
(16)
We can multiply them together for each term in the sum Eq. (10) and inverse Fourier transform to obtain $`p(\stackrel{~}{C}_l)`$ as

$$p(\stackrel{~}{C}_l)=\frac{1}{2\pi }_{\mathrm{}}^+\mathrm{}𝑑x\frac{e^{i\stackrel{~}{C}_lx}}{_{m=l}^l(12ix\sigma _{lm}^2)^{\frac{1}{2}}}$$
(17)
The log of the integrand in this expression (not including the complex exponential) is the cumulant generating function, which we can differentiate to obtain all cumulants of the pseudo–$`C_l`$ as

$$\kappa _n=2^{n1}(n1)!\underset{m}{}(\sigma _{lm}^2)^n.$$
(18)
We give the following expressions for the mean, variance, skewness $`\beta _1=\frac{(\mathrm{\Delta }\stackrel{~}{C}_l)^3}{(\mathrm{\Delta }\stackrel{~}{C}_l)^2^{\frac{3}{2}}}`$ and kurtosis $`\beta _2=\frac{(\mathrm{\Delta }\stackrel{~}{C}_l)^4}{(\mathrm{\Delta }\stackrel{~}{C}_l)^2^2}`$ as examples:

| $`\stackrel{~}{C}_l=_m\sigma _{lm}^2`$, | $`(\mathrm{\Delta }\stackrel{~}{C}_l)^2=2_m\sigma _{lm}^4`$ |
| --- | --- |
| $`\beta _1=2^{\frac{3}{2}}\frac{_m\sigma _{lm}^6}{\left(_m\sigma _{lm}^4\right)^{\frac{3}{2}}},`$ | $`\beta _2=12\frac{_m\sigma _{lm}^8}{\left(_m\sigma _{lm}^4\right)^2}`$. |

(19)
### D Pseudo–$`C_ł`$ distributions

It turns out that we do not have to contend with just the cumulants. We will now derive the analytical form for the pseudo–$`C_l`$ distributions. We define $`s_m^{(l)}=(2/\sigma _{lm}^2)`$ and note that $`s_m^{(l)}=s_m^{(l)}`$. Then all terms with equal $`\left|m\right|1`$ pair up and we obtain

$$p(\stackrel{~}{C}_l)=\frac{A^{(l)}}{2\pi }_{\mathrm{}}^+\mathrm{}𝑑x\frac{e^{i\stackrel{~}{C}_lx}}{(s_0^{(l)}ix)^{\frac{1}{2}}_{m=1}^l(s_m^{(l)}ix)}$$
(20)
where $`A^{(l)}\sqrt{s_0^{(l)}}_{m=1}^ls_m^{(l)}`$. We expand in partial fractions and integrate each term analytically (using a standard integral from ) to produce a closed form solution in terms of incomplete gamma functions $`\gamma (\alpha ,x)`$:

$$p(\stackrel{~}{C}_l)=A_{}^{}{}_{}{}^{(l)}\underset{m=1}{\overset{l}{}}\frac{\mathrm{exp}\left(s_m^{(l)}\stackrel{~}{C}_l\right)\gamma (\frac{1}{2},(s_0^{(l)}s_m^{(l)})\stackrel{~}{C}_l)}{\sqrt{(s_0^{(l)}s_m^{(l)})}_{m=1}^{}_{}{}^{}l(s_m^{}^{(l)}s_m^{(l)})}$$
(21)
where the primed product symbol $`^{^{}}`$ only multiplies factors which have $`mm^{}`$ and the normalization constant is

$$A_{}^{}{}_{}{}^{(l)}\frac{\sqrt{s_0^{(l)}}_{m=1}^ls_m^{(l)}}{2\mathrm{\Gamma }(\frac{3}{2})}.$$
(22)
We have, therefore, reduced the problem of determining the sample probability density of the pseudo–$`C_l`$ to computing the $`\sigma _{lm}^2`$. The computation of these coefficients has to be done only once for a given survey geometry and model theory. In fact it is clear from their definition and Eq. (11) that once the $`𝒲_{l^{}m}^{(l)}`$ are computed for a given survey geometry, the $`\sigma _{lm}^2`$ can be generated very easily given a model $`C_l`$ spectrum.

### E Correlations between different $`\stackrel{~}{C}_ł`$

The distributions we have just derived are 1-point or marginal distributions. They contain complete statistical information about each $`\stackrel{~}{C}_l`$ taken for itself, but they only account for the correlations between $`C_l`$ implicitly, giving an effective distribution for each $`\stackrel{~}{C}_l`$ which contains these correlations. To obtain explicit information about the full joint distribution of all $`\stackrel{~}{C}_l`$ we must go a step further.

While we wish to relegate further study of the properties of the joint distribution to future work we will now mention some preliminary results.

The characteristic function of the joint distribution of quadratic forms such as Eq. (12) has been studied in the literature . Fourier transforming we can obtain an integral representation of the joint distribution (analogously to Eq. (17) for the 1-point distributions)

$$p\left(\left\{\stackrel{~}{C}_l\right\}\right)=\frac{1}{(2\pi )^{l_{max}}}\frac{_l\left(e^{iS_l\stackrel{~}{C}_l}dS_l\right)}{\left|𝟙\mathrm{𝟚}𝕚_𝕝^{}𝕊_𝕝^{}^𝕝^{}\right|^{\frac{1}{2}}}$$
(23)
Even though we cannot proceed further in the Fourier inversion or must take recourse to approximate methods, we can still gain useful information about the joint statistics of the set of $`C_l`$. Similarly to subsection III C we can derive the cumulants of the variates from the logarithm of the integrand (disregarding the complex exponential). It is easy to show that for the covariance between the pseudo–$`C_l`$ we can write

$$\mathrm{\Delta }\stackrel{~}{C}_l\mathrm{\Delta }\stackrel{~}{C}_l^{}=\mathrm{tr}^{(l)}^{(l^{})}.$$
(24)
It turns out that higher order joint moments are computed similarly in terms of traces of products of $`^{(l)}`$ of various $`l`$. Note that these expressions are valid in general, but the special form of $`^{(l)}`$ in the case of azimuthal symmetry will reduce the scaling of the computational cost for their evaluation.

### F Approximations for the non–symmetric case

Returning to the case which we have fully solved, let us perform a reality check on the assumptions we have made. There are some important situations where the noise pattern does not follow the azimuthal symmetry of the survey geometry. In the case of the Planck satellite the scanning strategy is approximately centered on the ecliptic poles, while the Galactic cut is tilted through $`60^{}`$ with respect to this. In this case Eq. (21) becomes an approximation. We found it to be very accurate indeed to continue using these distributions with the $`\sigma _{lm}^2`$ computed for an asymmetric $`W_N`$, even for a strongly asymmetric noise pattern. This approximation will be worst in the least interesting, noise dominated regime at very high $`l`$. Note that the $`\sigma _{lm}^2`$ and hence the $`\stackrel{~}{C}_l`$ remain exact (because the $`\lambda _{lm}`$ are independent of the azimuth) but the remarks leading to Eq. (21) are no longer exactly true. For applications the final justification comes from the excellent agreement we find when we check against our Monte Carlo simulations which is illustrated in Figure 3. In fact all the applications we present in the following were computed using strongly asymmetric noise patterns.

## IV Monte–Carlo simulations

To test our formalism we performed 3328 Monte Carlo (MC) simulations for a high resolution CMB satellite, such as MAP or Planck (resulting in a sky coverage comparable to COBE). We simulated realizations of the CMB sky in the standard cold dark matter model($`\mathrm{\Omega }_m=1`$, $`\mathrm{\Omega }_bh^2=0.015`$, $`H_0=70`$km/s/Mpc). From these maps we carved out a $`\pm 20^{}`$ Galactic cut and contaminated the remaining area with spatially modulated Gaussian white noise of maximum root mean squared temperature $`124\mu K`$ per pixel of characteristic size 3.4 arcminutes. We use a tilted noise template $`W_N=\sqrt{\mathrm{sin}\theta _E}`$, where $`\theta _E`$ is the ecliptic latitude, as a simple model of the noise pattern which would result from scanning along meridians through the ecliptic poles. We then Fourier analysed these maps and stored the resulting $`\stackrel{~}{C}_l`$.

To give a visual impression of the resulting probability densities we show four cases in Figure 3. These cases were chosen to represent different regimes. The first case, $`l=2`$ is especially sensitive to effects due to the cut sky and an important region for the overall normalisation of the power spectrum. Two cases for intermediate $`l`$ ($`l=200`$, $`l=310`$) probe cosmologically interesting scales (on top of the first acoustic peak and in the first trough of the standard cold dark matter spectrum we use). At high $`l`$, $`l=1000`$, there is a significant noise contribution.

Next to our analytical results we show the histogrammed results from our MC simulations. A Kolmogorov–Smirnov test failed to detect deviations between the distributions of this MC population and Eq. (21) at 99% confidence, which validates our semi–analytic expressions.

Also shown in Figure 3 are the $`\chi ^2`$ distributions which the $`C_l`$ would follow in the full sky case as well as the commonly used Gaussian approximation . These are mean adjusted to account for the lost solid angle due to the Galactic cut. At $`l30`$ the difference is striking. For higher $`l`$ the Gaussian approximation becomes better as higher moments die away by dint of the Central Limit Theorem, but there remain visible systematic differences to the true distributions. In particular, there is a residual shift in the mean and the approximations tend to be slightly narrower than the histograms for very high $`l1000`$.

This becomes a more quantitative observation when looking at the percentage discrepancies between the mean and variances as a function of $`l`$ in Figure 4. The discrepancy is of the order of 1 % in the mean and 5% in the standard deviation on most scales, except for $`l<70`$ where the effect is larger. These discrepancies are important at the level of precision of future almost full sky missions. For medium and small sky coverage the mode couplings are stronger and we expect this to have an even larger effect on the probability distributions. We also compare the skewness $`\beta _1`$ and kurtosis $`\beta _2`$ of the $`\chi ^2`$ distributions to our distributions. The percentage difference is larger than for the first two moments but arguably less important at large $`l`$, since $`\beta _1`$ and $`\beta _2`$ decay as $`(2l+1)^{\frac{1}{2}}`$ and $`(2l+1)^1`$, respectively.

## V Applications

### A One–dimensional parameter estimation

As a first application we study the effect of approximating the likelihood for parameter estimation. Since our distributions have the correct means, we simply multiply them together for a simple, unbiased approximation to the likelihood

$$\stackrel{~}{}=\underset{l}{}p(\stackrel{~}{C_l})$$
(25)
This is a conservative approximation in the sense that we will not overestimate the estimation accuracy since the marginal distributions have all correlations integrated out and we will therefore overestimate the error bars on the $`C_l`$. Using this likelihood, as well as the Gaussian and $`\chi ^2`$ approximations, we attempt to estimate the baryon parameter $`\mathrm{\Omega }_b`$ (holding all other parameters constant) from several randomly selected realizations in our MC pool. The results are shown in Figure 5. As expected, Gauss and $`\chi ^2`$ consistently find estimates which are biased about 1.6% high, 3 standard errors of the mean away from the true value, while our likelihood gives a perfect fit.

Since the moment discrepancies depend on the underlying cosmological theory, this level of bias will also depend on the true theory. Therefore, this number should only be taken as indicative of the general level of the error introduced by using the naive Gaussian or $`\chi ^2`$ approximations.

### B Multi–parameter estimation

We demonstrated in the previous subsection that we can use our formalism to construct approximate forms for the likelihood which lead to an unbiased parameter estimate. The idea proposed in this subsection is motivated by the fact that higher moments of the $`\stackrel{~}{C}_l`$ distributions die away quickly with $`l`$ (cf. Eq. (19)) especially for observation geometries with large sky coverage. In particular we found in the case we studied that the distributions were visually indistinguishable from Gaussians for $`l100`$. For parameters which are sensitive to high $`l`$ modes, we therefore further approximate Eq. (25) by replacing the pseudo–$`C_l`$ pdfs with Gaussians with the same mean and the same variance,

$$\widehat{}(C_l)=\underset{l>100}{}\mathrm{exp}\left[\frac{\left(C_l\stackrel{~}{C_l}\right)^2}{\mathrm{\Delta }\stackrel{~}{C}_l^2}\right]$$
(26)
In this approximation, maximum likelihood estimation has reduced to simple $`\chi ^2`$ fitting, however with the correct means and variances which implicitly account for inter-$`C_l`$ correlations.

The operation count for each likelihood evaluation has now dropped to $`N_{pix}^{\frac{3}{2}}`$. If the computation of the first and second pseudo-$`C_l`$ moments is organised in a memory–efficient manner the elements of the coupling matrix need to be computed only once. We find that we can evaluate Eq. (26) at a rate of several hundred times per hour. This surprisingly slow scaling of the computational cost with number of theories is presumably due to more efficient use of on–processor caching.

To illustrate, we solve the problem of estimating 3 parameters ($`\mathrm{\Omega }_c`$ $`\mathrm{\Omega }_b`$ and $`H_0`$) simultaneously from a sky with $`12\times 10^6`$ pixels. To compare with the naive Gaussian approach and to show that our method is unbiased, we compute maximum (approximate) likelihood estimates from 100 realisations of the sky and plot a representation of the the empirical distribution of parameter estimates in three dimensions in Figures 6 (naive $`\chi ^2`$) and 7 (our approach).

We find again that our approach is unbiased. The distributions of the estimates are clearly centered on the true values.

We stress that our approach avoids the usual difficulties of the Gaussian approximation as discussed in section III A. For example, even though we use the Gaussian approximation, which of course does not exclude negative $`C_l`$, they are assigned an exceedingly small probability. This is because no attempt is made to subtract out the noise contribution from the pseudo–$`C_l`$ — instead it is modeled consistently and the (signal$`\times `$noise) cross term which is present in each realisation is not allowed to dominate.

### C Test for non–Gaussianity

We now mention an application which reverses our way of thinking about the results in this paper and takes advantage of the exact analytical expressions we compute for the pseudo–$`C_l`$ distributions. Imagine that we have gained knowledge that a certain theory is correct, for example through measurements of the large scale structure of the galaxy distribution and other probes of cosmology such as supernovae etc. Then we can use the formalism shown here as a test of Gaussianity, simply by computing the pseudo–$`C_l`$ and then substituting them into Eq. (21). An interesting feature is that this probes Gaussianity scale by scale which might make it feasible to disentangle confusion effects if the non–Gaussian signal was only apparent on a small range of scales and masked by a Gaussian component on other scales.

### D De–biasing

In the introduction we mentioned the use of iterative methods to attack the difficult problem of power spectrum estimation. A good starting guess can drastically speed up the convergence of such methods. In the case of large sky coverage, where the biases which occur in the $`\stackrel{~}{C}_l`$ are small, we suggest the following recipe as a “poor man’s power spectrum estimator” for de–biasing the power spectrum in only $`𝒪\left(N_{pix}^{\frac{3}{2}}\right)`$ computational time:

1. compute the $`\stackrel{~}{C}_l`$ from the data,
2. fit a smooth curve $`C_l^{smooth}`$ through them,
3. compute $`C_l^s=(\stackrel{~}{C}_lC_l^{noise})`$ and set all negative $`C_l^s`$ to zero,
4. use the $`C_l^s`$ instead of $`C_l^{theory}`$ in Eq. (19) to compute $`\stackrel{~}{C}_l`$.

Then

$$\widehat{C}_l=\stackrel{~}{C}_l+(C_l^{smooth}\stackrel{~}{C}_l)$$
(27)
is an estimator of the underlying theory $`C_l`$. There is a tiny residual bias which is second order in the discrepancy shown in Figure 4. This is utterly negligible for upcoming satellite missions.

## VI Further discussion and conclusions

We have computed the sampling properties of a set of quantities which we call pseudo–$`C_l`$. We have shown the usefulness of these quantities for the easy compression of large Cosmic Microwave Background data sets as well as for the forecasting of the discriminatory power of planned experiments. For example, in situation where percentage accuracy is required, the pseudo–$`C_l`$ sample variances shown in Figure 4 could immediately replace the usual approximation which simply rescales the cosmic variances for the full sky by the sky fraction.

On the other hand this exact framework can be used to justify or derive more approximate methods. We have given an example for the construction of an approximate but unbiased likelihood for cosmological parameters on the basis of a Gaussian approximation to the marginalised pseudo–$`C_l`$ probability distributions. Using this approximation we can evaluate the likelihood for hundreds of theories per hour on a single CPU and find unbiased estimates.

If only a rough treatment of the $`C_l`$ statistics is sufficient, for example in the case of large sky coverage, the results derived here can be used as a justification for other simplifying assumptions which are already in use in the literature.

Apart from the obvious applications to experiments with small and medium sky coverage such as balloons or ground–based missions, many further uses of this framework are conceivable. For example, one could use the statistical framework to design ’optimal’ scanning strategies, encoded in $`W_N`$, and assess more realistically if secondary anisotropies will be detectable with future CMB missions. Finally, all ingredients are there to refine the approximation to the joint likelihood we used in this paper by taking into account the covariances Eq. (24), for example by using a multivariate Edgeworth expansion around the peak of the likelihood.

We note that numerical techniques have recently been developed  which allow the computational solution of the power spectrum estimation for high $`l`$ applications under similar assumptions to the ones which lead to the analytical results derived in this paper. Their methods are efficient for the case of almost full sky observations. A purely numerical approach still requires significant computational resources, especially if there is signal in modes with $`l>1000`$. The results presented in this paper should be seen as complementary to such calculations. An analytical framework admits a more fundamental approach to understanding and is a useful yardstick against which numerical work can be tested or from which approximate methods can be derived as we have demonstrated in the applications presented in section V.

Our methods do not allow for significantly correlated noise. Due to the current state of detector technology, all CMB missions suffer from correlated noise, with the possible exception of the MAP satellite . Apart from brute force numerical calculation with $`𝒪\left(N_{pix}^3\right)`$ operations, which is no longer feasible even with current experiments, there are currently no techniques for computing sampling statistics in the presence of correlated noise. This clearly needs to be addressed as an extremely urgent and important issue.

To summarize, we have presented a theoretical framework for the study of power spectrum statistics which is applicable regardless of the size of the sky area covered. We go beyond current approximations and present a semi–analytic formalism for the computation of sampling distributions of the $`C_l`$ for any Gaussian cosmological model and a large and important class of surveying strategies. We show that these results can be usefully applied to the estimation of cosmological parameters from simulated data. A range of further applications is suggested which demonstrate the power of our formalism.

###### Acknowledgements.

We wish to thank A. J. Banday for stimulating discussions. This research was supported by Dansk Grundforskningsfonden through its funding for TAC.

## A Low order non–cosmological modes

We assume in this paper that we have access to a partial sky map of CMB anisotropies free from non–cosmological signal. Even in the event that foreground contributions can be filtered out using the frequency information provided by multi–channel experiments, this assumption merits further discussion — particularly for large amplitude low order modes such as the monopole and the dipole which are thought to consist exclusively of non-cosmological signal.

An ideal full sky map of the cosmic microwave background anisotropies has zero mean, by definition. Also, in Friedmann–Robertson–Walker models, there is no cosmological dipole. But we expect realistic maps resulting from observations of the microwave sky to contain both monopole and dipole components; Galactic emission is positive definite and will therefore give a monopole contribution. The motion of our solar system with respect to the frame in which the CMB is at rest produces a dipole. The question arises how to deal with the presence of such non–cosmological low order modes. <sup>*</sup><sup>*</sup>*We should mention that there is of course a natural occurrence of spurious monopole and dipole components even in a completely uncontaminated map as soon as a part of the sky is cut away. This is easy to understand: the monopole is only constrained to be zero if the anisotropy is integrated over the full sky. Similarly, a cut sky has a spurious dipole, even if the dipole was zero on the full sky. This statistical effect is fully taken into account in the formalism as it is presented here. All our Monte Carlo simulations of cut skies have monopole and dipole components even though they vanish exactly on the full sky. The sampling distributions of $`\stackrel{~}{C}_0`$ and $`\stackrel{~}{C}_1`$ containing this effect are quoted in Eq. (21).

If the experiment measures total power, we may hope to remove a sizeable part of the monopole using the frequency information of a multi–channel experiment. In the case of a differential experiment such as MAP which is blind to the monopole component, the problem depends on the particulars of the map–making algorithm employed. Frequency information alone cannot distinguish between a non-cosmological and a spurious dipole since both are due to the microwave background. Studies of the peculiar velocity field may provide some information.

We see the three following different approaches which can be used to deal with the problem of mode to mode coupling when non–cosmological components are removed from the sky.

One is orthogonalisation . This solves the problem, but is not feasible computationally when dealing with millions of degrees of freedom.

The other is projecting out the undesired modes either by fitting and subtracting or more formally by operating with projection operators. Those approaches lead to a correlation structure for the remaining modes which is not analytically tractable in the manner presented in this paper.

The third is allowing for the presence of a monopole and a dipole in a Bayesian fashion, by including a rough estimate of their size in the theory $`C_l`$. In particular, a rough a priori measure of the monopole and dipole present in the data can be input into $`C_0`$ and $`C_1`$. Our statistical framework allows this and yields the correct predicted sampling statistics for all $`C_l`$ taking into account all the $`l`$$`l^{}`$ couplings due to the cut. In effect, $`C_0`$ and $`C_1`$ are treated as nuisance parameters which are marginalised over in our formalism.

This will change the sampling statistics of low order modes predicted by our theory in a self-consistent way. Cosmological parameters can then be estimated in an unbiased fashion without monopole or dipole removal. Whether this is feasible in practice would be an interesting avenue for further study. We are not aware of any existing fast method which offers a similarly consistent treatment.

In practice we observe that for large sky coverage the entries of the coupling matrices diminish rapidly with increasing $`l`$. To a good approximation, modes with $`l100`$ are therefore unaffected by the presence or absence of low order modes, $`l<2`$. One of the main applications of our formalism has been to the estimation of parameters which depend mostly on the shape of the $`C_l`$ spectrum around and beyond the first acoustic peak. These analyses are unaffected by the details of how low order modes of non–cosmological origin are dealt with.

## B Factorising the correlation matrix

For a general azimuthally symmetric window $`W(\theta )`$ we can write the $`\stackrel{~}{C}_l`$ as

$$\begin{array}{cc}\hfill \stackrel{~}{C}_l=\underset{m}{}𝑑\mathrm{\Omega }_1W(\theta _1)& d\mathrm{\Omega }_2W(\theta _2)Y_{lm}(\gamma _1)Y_{lm}^{}(\gamma _2)\hfill \\ & \underset{l^{}m^{}}{}\underset{l^{\prime \prime }m^{\prime \prime }}{}a_{l^{}m^{}}a_{l^{\prime \prime }m^{\prime \prime }}^{}Y_{lm}(\gamma _1)Y_{l^{\prime \prime }m^{\prime \prime }}^{}(\gamma _2)\hfill \end{array}$$
(B1)
We define normalised associated Legendre polynomials $`\lambda _{lm}(x)`$ in Appendix C and use the notation $`\mu =\mathrm{cos}\theta `$. Then $`Y_{lm}(\theta ,\varphi )=\lambda _{lm}(\mu )e^{im\varphi }`$ and we notice that the Kronecker–$`\delta `$ from the azimuthal integration select only one term out of the sums over $`m^{}`$ and $`m^{\prime \prime }`$, namely

$$\stackrel{~}{C}_l=\underset{m}{}\left(\underset{l^{}=0}{\overset{\mathrm{}}{}}a_{l^{}m}^{}𝒲_{l^{}m}^{(l)}\right)\left(\underset{l^{\prime \prime }=0}{\overset{\mathrm{}}{}}a_{l^{\prime \prime }m}𝒲_{l^{\prime \prime }m}^{(l)}\right)$$
(B2)
where the

$$𝒲_{l^{}m}^{(l)}W_{l^{}mlm}=𝑑\mu W(\mu )\lambda _{lm}(\mu )\lambda _{l^{}m}(\mu )$$
(B3)
are real matrices. Because the $`W`$ matrices are real, the terms in the brackets are complex conjugates of each other and we obtain Eq. (10) as desired.

## C  Algorithm for computing the coupling matrix

Subsection C 1 reviews some useful properties of spherical harmonics. In subsection C 2 we derive an efficient algorithm for computing the coupling matrix Eq. (B3) due to an azimuthally symmetric mask, Eq. (C16). More general windows can be piecewise defined in terms of masks of varying height. The corresponding coupling matrices are then suitably weighted sums over the coupling matrices due to the individual masks.

### 1 Some properties of spherical harmonics

The Spherical Harmonics are

$$Y_{lm}(\theta ,\varphi )=\lambda _{lm}(\mathrm{cos}\theta )\mathrm{e}^{im\varphi }$$
(C1)
where the

$`\lambda _{lm}(x)`$ $`=`$ $`\sqrt{{\displaystyle \frac{2l+1}{4\pi }}{\displaystyle \frac{(lm)!}{(l+m)!}}}P_{lm}(x),\mathrm{for}m0`$ (C2)
$`\lambda _{lm}`$ $`=`$ $`(1)^m\lambda _{l|m|},\mathrm{for}m0,`$ (C3)
$`\lambda _{lm}=0,\mathrm{for}|m|>l.`$ (C4)
In the following we will only consider positive $`m`$. The associated Legendre Polynomials $`P_{lm}`$ are defined as

$$P_{lm}(x)=(1)^m(1x^2)^{m/2}\frac{d^m}{dx^m}P_l(x),$$
(C5)
where $`P_l`$ is the Legendre function of the first kind defined as

$$P_l(x)=\frac{1}{2^ll!}\frac{d^l}{dx^l}(x^21)^l.$$
(C6)
They are solution of the differential equation

$$(1x^2)\frac{d^2}{dx^2}P_{lm}2x\frac{d}{dx}P_{lm}+\left(l(l+1)\frac{m^2}{1x^2}\right)P_{lm}=0$$
(C7)
Definition Eq. (C5) leads to a large number of relations between polynomials of different order $`l`$ and $`m`$ and their derivatives. Among them the following two will be useful (see, e.g.   §8.700)

$$P_{lm}=\frac{x(m1)}{\sqrt{1x^2}}P_{l(m1)}\sqrt{1x^2}\frac{d}{dx}P_{l(m1)}$$
(C8)
and

$$(1x^2)\frac{d}{dx}P_{lm}=lxP_{lm}+(l+m)P_{(l1)m}$$
(C9)
It is possible to obtain the $`P_{lm}`$ or alternatively the $`\lambda _{lm}`$ through relations that raise from the fact that they are orthogonal polynomials. The following recurrence is stable and very convenient for numerical uses :

starting with $`\lambda _{00}=1/\sqrt{4\pi }`$

$`\lambda _{mm}`$ $`=`$ $`\sqrt{{\displaystyle \frac{2m+1}{2m}}}\sqrt{1x^2}\lambda _{(m1)(m1)}`$ (C10)
$`\lambda _{(m+1)m}`$ $`=`$ $`x\sqrt{2m+3}\lambda _{mm}`$ (C11)
$`\lambda _{lm}`$ $`=`$ $`\left[x\lambda _{(l1)m}{\displaystyle \frac{\lambda _{(l2)m}}{A(l1,m)}}\right]A(l,m)\mathrm{with}A(l,m)=\sqrt{{\displaystyle \frac{4l^21}{l^2m^2}}}`$ (C12)
### 2 Scalar product on the cut sphere

The spherical harmonics have the property that they are an orthonormal basis spanning the Hilbert space of square integrable functions on the sphere $`S^2`$, i.e.  their scalar product is the identity matrix:

$$_0^{2\pi }𝑑\varphi _1^1d\mathrm{cos}\theta Y_{l^{}m^{}}(\theta ,\varphi )Y_{lm}(\theta ,\varphi )=\delta _{l^{}l}\delta _{m^{}m}$$
(C13)
Let us now define the cut sphere scalar product of spherical harmonics under a mask or window $`W`$

$$A_{l^{}l}^{mm^{}}[W]_0^{2\pi }𝑑\varphi _1^1d\mathrm{cos}\theta W(\theta ,\varphi )Y_{l^{}m^{}}(\theta ,\varphi )Y_{lm}(\theta ,\varphi )$$
(C14)
In the following we will consider only a sharp ’strip’ filter whose boundaries are parallel to the equator (but not necessarily symmetric about the equator).

$`W(\mathrm{cos}\theta )=W(x)`$ $`=`$ $`1\mathrm{if}axb`$ (C15)
$`=`$ $`0\mathrm{otherwise}`$ (C16)
This is simply related to physically relevant windows. For example a straight cut removing foregrounds concentrated around the Galactic disk is modeled by $`\stackrel{~}{W}=1W`$ with $`a=b`$ setting the width of the cut. Because $`A_{l^{}l}^{mm^{}}[\stackrel{~}{W}]=\mathrm{𝟏}A_{l^{}l}^{mm^{}}[W]`$, where $`\mathrm{𝟏}=\delta _{l^{}l}\delta _{m^{}m}`$, everything that follows can be translated to this window with the replacements

$$ba,\mathrm{and}A_1^1[\stackrel{~}{W}]=1.$$
(C17)
Because the window Eq. (C16) is independent of $`\varphi `$ it preserves the azimuthal symmetry and only spherical harmonics with the same $`m`$ are coupled. The non–trivial terms of the coupling matrix $`A`$ are then

$`A_{l^{}l}^m`$ $`=`$ $`2\pi {\displaystyle _a^b}𝑑x\lambda _{l^{}m}(x)\lambda _{lm}(x)`$ (C18)
$`=`$ $`{\displaystyle \frac{\sqrt{(2l^{}+1)(2l+1)}}{2}}\sqrt{{\displaystyle \frac{(l^{}m)!(lm)!}{(l^{}+m)!(l+m)!}}}I_{l^{}l}^m(a,b)`$ (C19)
where

$$I_{l^{}l}^m(a,b)_a^b𝑑xP_{l^{}m}(x)P_{lm}(x)$$
(C20)
This expression is symmetric under interchange of $`l^{}`$ and $`l`$ and under $`mm`$. If $`0ll_{max}`$ this implies that, asymptotically for large $`l_{max}`$ there are $`\frac{l_{max}^3}{6}`$ independent components of $`I_{l^{}l}^m(a,b)`$. If the window $`W`$ is symmetric around the equator and therefore preserves reflection symmetry, only multipoles with the same parity of $`l`$ couple to one another, and the number of non–trivial elements reduces to $`l_{\mathrm{max}}^3/12`$.

The problem of evaluating all independent components of Eq. (C20) is to find a way of computing the integrals Eq. (C20) without having to take recourse to time–consuming numerical quadratures.

If one substitutes Eq. (C8) for $`P_{lm}`$ and $`P_{l^{}m}`$ in Eq. (C20), integrates by parts to exhibit the second derivative of $`P_{l(m1)}`$, and makes use of the differential equation Eq. (C7) one obtains

$`I_{l^{}l}^m(a,b)`$ $`=`$ $`(m+l)(ml1)I_{l^{}l}^{m1}(a,b)`$ (C21)
$`+`$ $`(ml1)\left[xP_{l^{}(m1)}P_{l(m1)}\right]_a^b`$ (C22)
$`+`$ $`(m+l1)\left[P_{l^{}(m1)}P_{(l1)(m1)}\right]_a^b`$ (C23)
where we used Eq. (C9) to simplify the last two terms.

This can be simplified further in the case $`ll^{}`$, by noting that the left hand side of Eq. (C23) is by definition symmetric in $`l^{},l`$ whereas the right hand side is not explicitly symmetric. So, by subtracting Eq. (C23) from itself after swapping $`l`$ and $`l^{}`$ one obtains

$`I_{l^{}l}^m(a,b)`$ $`=`$ $`\left((l^{}l)\left[xP_{l^{}m}P_{lm}\right]_a^b+(l+m)\left[P_{l^{}m}P_{(l1)m}\right]_a^b(l^{}+m)\left[P_{(l^{}1)m}P_{lm}\right]_a^b\right)`$ (C24)
$`/`$ $`\left(l_{}^{}{}_{}{}^{2}+l^{}l^2l\right)\mathrm{for}l^{}l.`$ (C25)
The diagonal terms ($`I_l^m(a,b)=_a^bP_{lm}P_{lm}`$) can be obtained thanks to the recurrence relation

$$(lm)P_{lm}=x(2l1)P_{(l1)m}(l+m1)P_{(l2)m}.$$
(C26)
One obtains a recurrence on the diagonal terms involving also the off–diagonal terms computed by Eq. (C25).

$`I_l^m(a,b)`$ $`=`$ $`\left({\displaystyle \frac{2l1}{2l+1}}\right)\left({\displaystyle \frac{l+m}{lm}}\right)I_{l1}^m(a,b)`$ (C27)
$`+`$ $`\left({\displaystyle \frac{2l1}{2l+1}}\right)\left({\displaystyle \frac{lm+1}{lm}}\right)I_{(l+1)(l1)}^m(a,b){\displaystyle \frac{l+m1}{lm}}I_{l(l2)}^m(a,b)`$ (C28)
To conclude, the elements of the coupling matrix $`A`$ are given in terms of recurrence relations as follows. For $`mm^{}`$ we have

$$A_{l^{}l}^{m^{}m}=0,$$

while for $`ll^{}`$

$$\begin{array}{cc}\hfill A_{l^{}l}^m& =((l^{}l)\left[x\lambda _{l^{}m}\lambda _{lm}\right]_a^b+\sqrt{\frac{2l+1}{2l1}(l^2m^2)}\left[\lambda _{l^{}m}\lambda _{(l1)m}\right]_a^b\hfill \\ & \sqrt{\frac{2l^{}+1}{2l^{}1}(l_{}^{}{}_{}{}^{2}m^2)}\left[\lambda _{(l^{}1)m}\lambda _{lm}\right]_a^b)\times \frac{2\pi }{(l^{}l)(l^{}+l+1)},\hfill \end{array}$$
(C29)
and all other cases are given by

$$\begin{array}{cc}\hfill A_l^m& =A_{l1}^m+\sqrt{\frac{2l1}{2l+3}\frac{(l+1)^2m^2}{l^2m^2}}A_{(l+1)(l1)}^m\sqrt{\frac{2l+1}{2l3}\frac{(l1)^2m^2}{l^2m^2}}A_{l(l2)}^m(l>m)\hfill \\ \hfill A_m^m& =A_{m1}^{m1}+\frac{2\pi }{2m+1}\left[x\lambda _{mm}\lambda _{mm}\right]_a^b\hfill \end{array}$$
(C30)
with $`A_1^1[W]0`$. As only a few operations are needed to evaluate each of these terms, calculating the whole coupling matrix costs $`O(l_{\mathrm{max}}^3)`$ operations and only requires about 10 minutes for $`l_{\mathrm{max}}=2048`$ on a fast work–station.

## D The method of characteristic functions

### 1 Densities of functions of random variates

How do we approach computing the probability density function (pdf) of a (set of) continuous random variate(s) which can itself be written as a function of one or more random variates whose pdfs are known? The usual expression for changing the measure is

$$p(f(x))=p(x)\left|\frac{f}{x}\right|^1$$
(D1)
where $`\left|\frac{f}{x}\right|`$ is the Jacobian of the transformation. However this expression is inconvenient if the transform is not monotonic or the number of transformed variates is not the same as the number of original variates.

In order to attack this more general situation we replace Eq. (D1) by a more powerful integral expression. The idea here is to define an ”atomic” pdf which allows for only one event. For a continuous variate this is the delta function. The sum of these ”atoms” over all possible outcomes, weighted with the probability that the underlying variates produce this outcome gives the resulting pdf. In symbols, if the $`n`$ underlying variates $`X=\{x_i,i=1,\mathrm{},n\}`$ have the joint pdf $`p(X)`$, then the pdf of $`z=f(X)`$ is

$$p(z)=\delta \left(zf(X)\right)p(X)𝑑X,$$
(D2)
where $`dX`$ denotes $`_{i=1}^{i=n}dx_i`$.

### 2 The Method of characteristic functions

In particular we are often interested in the distribution of a linear combination of $`n`$ independent quantities whose distributions are known. Computing these distributions is now straightforward – we can apply Eq. (D2) with the simplification that $`f(X)`$ is linear and p(X) factorises owing to the independence of the variates. The trick is to perform the resulting convolution in Fourier space by application of the Fourier convolution theorem. This is far easier to do than computing $`n`$ convolutions, at least if $`n>2`$. The Fourier transform of a pdf of a variate is called its characteristic function — hence the name of the method. All major pdfs have tabulated Fourier transforms. This makes this method very simple to apply, up to the final step of inverse transforming to obtain the required pdf. This may be tabulated too. If not, extending the integrand to the complex plane and using contour integration perhaps after expanding into partial fractions may lead to success. In any case, at least only one integral needs to be evaluated compared to $`n`$ for the direct convolution.

Even if the inverse transform does not succeed, all is not lost. A key fact about characteristic functions allows us to gain as much information as we like about the statistical properties of $`z`$. By formally computing the Maclaurin series of the Fourier transform which defines the characteristic function $`c(k_z)`$ of a variate $`z`$, we find that it is the generating function for its statistical moments,

$$c(k_z)=\underset{n=0}{\overset{n=\mathrm{}}{}}z^n\frac{(ik_z)^n}{n!}.$$
(D3)
Therefore, the moments can be obtained by straightforward differentiation of $`c(k_z)`$. Equivalently, and sometimes more conveniently, the natural logarithm of the characteristic function generates the cumulants of the distribution. Moments and cumulants are algebraically simply related. Once these moments or cumulants are known there is a host of methods for using them to approximate $`p(z)`$ . In particular, when $`n`$ is large, the large number theorem states that $`p(z)`$ approaches a Gaussian with mean $`z`$ and variance $`z^2`$ (under weak conditions on the $`p_i(x_i)`$).