text:
# Chandra Observation of the Merging Cluster Abell 2034

## 1. Introduction

Clusters of galaxies are bright X-ray sources, with typical bolometric luminosities in the range from $`10^{43}`$ to $`10^{45}`$ erg s-1. In the standard hierarchical model of cluster formation, larger clusters are formed by hierarchical accretion of smaller clusters. Occasionally, two clusters of roughly equal mass will merge. These major mergers are expected to have significant effects on the intracluster medium (ICM). Mergers will drive shocks with Mach numbers $`2`$, which can significantly heat the ICM. While strong, centrally condensed cooling flows are anticorrelated with major mergers, merger shocks are too weak to completely disrupt cooling flows at the centers of clusters. Consequently, cold cores of gas have been observed to be moving at high velocities through the ICM of several clusters e.g. Abell 2142, Abell 3667; Markevitch et al., 2000; Vikhlinin et al., 2001a. In principle, these thermal effects should all be visible in X-rays: the shocks as abrupt increases in the entropy of the X-ray emitting gas, and cold cores as regions of high density, cold gas.

The survival of these cold cores during the merger process has been one of the most significant discoveries in the area of cluster research made by Chandra. Prior observations did not have the spatial resolution to resolve these cold cores of gas from the surrounding hotter gas. The leading edges of these cold cores, dubbed “cold fronts,” are contact discontinuities, where the pressure gradient from the hot gas to the cold gas is continuous. The abruptness of the temperature and density change over the small distance of these fronts means that transport processes (diffusion) across the fronts must be suppressed (Vikhlinin et al., 2001b; Ettori & Fabian, 2000).

Abell 2034 is a moderate redshift, high X-ray luminosity cluster at $`z=0.113`$ (Struble & Rood, 1999). Previous observations of the cluster with ASCA and ROSAT have found, respectively, a cooling flow-corrected temperature of 9.6 keV (White, 2000) and a bolometric luminosity of $`2.2\times 10^{45}`$ erg s-1 (David, Forman, & Jones, 1999). White (2000) found a best-fit single temperature of 7.6 keV. Optically, it is quite rich (ACO richness class 2; Abell, Corwin, & Olowin, 1989). It has a cD galaxy which is offset $`1\mathrm{}`$ from the X-ray centroid and $`2.5\mathrm{}`$ from the optical centroid of the cluster. The relative displacements of the cD galaxy, average galaxy position, and center of the X-ray gas distribution suggest that the cluster is out of equilibrium. Indeed, non-equilibrium features were seen in the cluster in a pointed ROSAT observation; the unpublished ROSAT PSPC image is presented below (Figure 2). It showed an excess of emission to the south of the cluster center which was suspected to be a merger shock. Also visible in the ROSAT image was a surface brightness discontinuity to the north of the cluster center, opposite the suspected merger shock. The detection of a radio relic near the position of the northern discontinuity (Kempner & Sarazin, 2001) strengthens the case for the cluster being out of equilibrium due to a recent or ongoing merger.

Despite these pieces of evidence for a merger, Abell 2034 has also been reported to have a large cooling flow with a cooling rate of $``$90–580 $`M_{\mathrm{}}`$ yr-1 (White, 2000). Usually, large cooling flows are strongly anti-correlated with major cluster mergers (Buote & Tsai, 1996). Thus, ASCA and ROSAT observations did not lead to a consistent picture of the dynamical state of this cluster. Higher spatial resolution observations with spatially resolved spectra and a broad X-ray band (given the high temperature) were needed.

This makes Abell 2034 an interesting target for a Chandra observation, particularly since it is located at a sufficiently large redshift that most of the center of the cluster can be imaged with the ACIS-I camera. In this paper, we present the results of such an observation. All errors are quoted at 90% confidence unless otherwise stated. We assume $`H_0=50`$ km s-1 Mpc-1 and $`q_0=0.5`$ throughout this paper. At the redshift of Abell 2034, $`1\mathrm{}160`$ kpc.

## 2. Observation and Data Reduction

Abell 2034 was observed on 2001 May 5 in a single 54 ksec exposure. The data were taken in Very Faint (VF) mode using the four ACIS-I chips and the S3 chip. The focal plane temperature was $``$120 C, the frame time was 3.2 s, and only events with ASCA grades of 0,2,3,4, and 6 were analyzed. The S3 chip was used to check the background during the observation; no background flares were observed. A binned image, uncorrected for background or exposure, is shown in Figure 1a.

Chandra data taken in VF mode retain the full 5$`\times `$5 pixel event islands, which allow for more thorough rejection of particle background than do data with only the 3$`\times `$3 pixel event islands111http://cxc.harvard.edu/cal/Links/Acis/acis/Cal\_prods/vfbkgrnd/index.html. We applied this additional filtering to our data and to the blank-sky backgrounds. We determined the background using the blank sky background files included in the Chandra Calibration Database222http://asc.harvard.edu/caldb/. We used the period D background files, which include data taken in VF mode.

The ancillary response files (arfs) used in the spectral analysis were corrected for the reduction in transmission of the CCD window at low energies using the “apply\_acisabs” tool in CIAO. At the time of the observation, the reduction in sensitivity due to this effect was about 20% at 0.8 keV. Due to additional uncertainty in the calibration at low energies, we restricted our analysis to the range 0.8–9.0 keV.

## 3. X-ray Image

The raw Chandra image, binned by a factor of 4 and uncorrected for background or exposure, is shown in Figure 1a. At its moderate redshift, Abell 2034 fills much of the ACIS-I detector on Chandra. The cluster was positioned on the detector so that the aimpoint would be about halfway between the cD galaxy and the centroid of the X-ray emission as determined from the ROSAT image. The image shows several interesting features in the cluster itself, and a number of point sources.

A smoothed image of the diffuse emission is shown in Figure 1b. The image was adaptively smoothed using the csmooth tool from the CIAO software package333http://asc.harvard.edu/ciao/. Blank sky background and exposure images were smoothed using the same kernel, and were then used to correct the smoothed raw image.

As can be seen in Figures 1a and 1b, the cluster has significant emission out to a distance of $`5\mathrm{}`$$`6\mathrm{}`$ in all directions. There is a sharp discontinuity in the surface brightness $``$$`3\mathrm{}`$ to the north of the cluster center. We will refer to this feature as the “northern cold front.” It is discussed in detail below (§6). This feature was seen in the earlier ROSAT image (Figure 2; the image was corrected for background and exposure using the ROSAT extended object analysis tools by S. Snowden444ftp://legacy.gsfc.nasa.gov/rosat/software/fortran/sxrb/). This feature was one of the motivations for observing this cluster with Chandra. Based on the ROSAT observations, which provided no useful information on the X-ray spectra near this feature, we had interpreted this feature as a possible merger shock.

Also visible in the raw and smoothed X-ray images is an excess of emission to the south of the cluster center. We will refer to this region as the “south excess;” it is discussed in more detail below (§ 5). This feature was also seen in the earlier ROSAT image (Figure 2), where it appeared as a rather distinctive arc of emission running from the southeast to the northwest. Initially, we suspected that this feature might be a merger shock cutting across the cluster. The extensions of this arc in the ROSAT image are due, at least partially, to unresolved point sources seen in the Chandra image. Because Chandra is able to better resolve these distinct regions of emission and separate out point sources, the south excess appears as a fairly regular region of diffuse emission rather than as the arc which we believed to see in the ROSAT image (Figure 2).

The cluster does not have a particularly sharply peaked core, suggesting that it does not contain a strong cooling flow. The core of the cluster has a somewhat irregular structure. The central arcminute or so is shaped like a curved teardrop, with the brightness peak centered on the broad end of the teardrop, to the north of center. Outside of 1′, the cluster becomes more regular and circular. This structure is clearly visible in the contours in Figure 3.

An optical image of the Digital Sky Survey (DSS) is shown in Figure 3. Note that the cD galaxy which is the brightest cluster member is not located at the peak in the X-ray surface brightness, but is $``$77″ south of the peak. There is a smaller peak in the X-ray surface brightness which does coincide with the central cD galaxy.

A very bright cluster galaxy, also classified as a cD, is located NNW of the X-ray peak near the location of the northern cold front. There are several other bright elliptical galaxies in this neighborhood, and numerous fainter galaxies. In fact, the number of moderately bright to faint galaxies in this region significantly outnumbers that in the region of the more central cD galaxy. The majority of these galaxies are north, or ahead, of the cold front, perhaps suggesting that the collisionless galaxies have moved ahead of the collisional gas, which has been slowed during the merger. Just inside the cold front, there is a significant excess of X-ray emission relative to that expected from a smooth $`\beta `$-model, as we discuss later in §4.1.

There are no bright galaxies and no obvious excess in the galaxy density evident in the DSS image of the region covered by the south excess. Two faint galaxies ($`m_\mathrm{R}17`$) are visible on the western edge of the south excess, and are associated with the wide angle tail (WAT) radio galaxy FIRST J150957.2+332716 and the narrow angle tail (NAT) radio galaxy FIRST J150959.5+332746. The latter is also visible as an X-ray source in the Chandra image.

## 4. Global Cluster Properties

### 4.1. Radial Gas Distribution

The X-ray image of Abell 2034 shows a number of strong features which suggest that the cluster may not be relaxed. These include the northern cold front, the south excess, the lumpy structure near the center, and the displacement of the central cD galaxy from the X-ray peak. Despite these local irregularities, the radial surface brightness profile of the cluster is fairly consistent over most azimuthal angles.

We determined the radial surface brightness profile in 5 sectors where the cluster is approximately circularly symmetric within each sector. The sectors and annuli used to extract the surface brightness profile are shown in Figure 4a. The annuli were selected to have between 1000 and 2000 source counts per annulus in order to produce reasonable spectral fits. This meant that for the purposes of the surface brightness analysis presented here they had a much higher flux than was necessary. In Figure 4b, we show the radial surface brightness profiles in these five sectors. The sectors’ common origin is the approximate centroid of the cluster emission. A $`\beta `$-model fit to the sum of these profiles yields the best-fit core radius of $`r_c=290\pm 10`$ kpc and slope of $`\beta =0.69\pm 0.02`$ (1-$`\sigma `$ errors). We could not find any previously published values for these parameters for comparison. This model provides a reasonably good fit to all of the inner regions of the sectors except for the one encompassing the northern edge, plotted in red in Figure 4b. This particular sector shows an excess of emission immediately inside the edge and a smaller deficit outside, but at smaller and larger radii is reasonably consistent with the profiles in the other sectors.

The southwest sector has a noticeably flatter profile with brighter emission at large radii due to the south excess region. The south excess may in fact be a separate, possibly more distant cluster seen in projection against Abell 2034, as we discuss in detail in §5. In Figure 5, we show the residuals after subtracting a $`\beta `$-model from the Gaussian smoothed image shown in Figure 4a. Since the southeast sector is largely free of emission from the south excess and from other substructure within the cluster, we fit a $`\beta `$-model separately to this sector for use in creating Figure 5. The resultant $`\beta `$-model fit has $`r_c=190`$ kpc and $`\beta =0.54`$, which is a significantly smaller core radius and flatter slope than the best fit to all the sectors including the south excess and internal substructure. As is evident in Figure 5, however, this still oversubtracts slightly the large-scale diffuse emission in the cluster, particularly at radii less than $``$800 kpc. This oversubtraction is probably due to the $`\beta `$-model still being affected by a small contribution from the south excess and from the excess at the very center of the cluster.

### 4.2. Spectral Properties

To test the previous report using data from ASCA of a large cooling flow in Abell 2034 (White, 2000), we extracted a spectrum from a circular region with radius equal to the core radius determined in §4.1. We fit it with a single temperature MEKAL model (Kaastra, 1992; Liedahl, Osterheld, & Goldstein, 1995), and then with a single temperature plus a cooling flow. In both cases we fixed the absorption to the Galactic value of $`1.58\times 10^{20}`$ cm-2 (Dickey & Lockman, 1990). The single temperature fit determined a temperature of $`7.9\pm 0.4`$ keV, and an abundance of $`0.29\pm 0.07`$ times solar, with $`\chi ^2/d.o.f.=1.4`$. This abundance is consistent to within the errors with that determined by White (2000) using data from ASCA. The addition of a cooling flow component, with the minimum temperature set to essentially zero and the maximum temperature and abundance tied to the values of the single-temperature component, did not improve the fit. The resultant cooling rate to low temperature was determined to be $`23_{20}^{+21}M_{\mathrm{}}`$ yr-1. With the low temperature as a free parameter, the fit was too poorly contrained to provide any useful limits. In short, we find evidence for only a small amount of cooling in the center of the cluster. It is likely that the large cooling rate measured by White (2000) was due to the presence of cooler gas in the cold front (see §6) and in the south excess (see §5) rather than in a classical cooling flow in the cluster core.

We extracted spectra in the same sectors as were used in §4.1, and fit them with a single temperature MEKAL model. For the initial spectral fits, we fixed the absorption at the Galactic value. We first allowed only the temperature and normalization to vary, with the abundance fixed at the mean value of 0.29 times solar determined above. We then allowed the abundance to vary and re-fit the spectra. Allowing the abundances to vary did not change the temperatures we measured by a significant amount, and the abundances we determined were not well constrained—in most cases they provided only upper limits. There is no evidence from these spectra, however, for any large abundance gradients across the cluster, and all the abundances we measured were consistent with the value of 0.29 times solar determined above. We therefore use the fits with fixed abundance in our subsequent analysis. The results of these fits are shown in Figure 6.

With the exception of the southwest quadrant, where contamination from the south excess is significant at large radii, we find no evidence for a large scale temperature gradient like that seen in some clusters (Markevitch et al., 1998). The temperature gradients reported by Markevitch et al. (1998), however, were generally at radii $``$1 Mpc, whereas our temperature profiles are only significant out to $``$700 kpc.

### 4.3. Temperature Structure

In Figure 7, we present an adaptively binned temperature map of the cluster. The temperature map was created using the adaptive binning algorithm described in Houck & Denicola (2000) and Houck, Wise, & Davis (in preparation). The algorithm extracts a spectrum at each pixel in the output map from a region which is as small as possible while containing at least 800 counts. The maximum extraction region was 63″$`\times `$63″. The smallest extraction region used was 26$`\stackrel{}{\mathrm{.}}`$6$`\times `$26$`\stackrel{}{\mathrm{.}}`$6. For pixels around which it was not possible to extract 800 counts from a region of this maximum size, no spectrum was extracted and no temperature was calculated. Each spectrum was corrected for background using a spectrum from a matching region in the blank sky background file discussed in §2. Regions containing point sources were eliminated from the event lists used for the source and background spectra. The spectra were grouped to have a minimum of 20 photons per bin, then fit by a single MEKAL model with absorption fixed at the Galactic value, and abundance fixed at $`0.29Z_{\mathrm{}}`$. The negative errors in the temperature map are less than 15% everywhere. The positive errors are larger since the average temperature is close to or above the upper limit of the Chandra spectral band, i.e. $`9`$ keV. Where the temperature is less than $``$8 keV, the positive errors are between 15% and 30%.

The temperature map shows inhomogeneities on scales roughly the size of the binning scale and larger. Since the binning washes out any inhomogeneities on smaller scales, and since the inhomogeneities span a range of scales up to several times the binning scale, it is likely that better photon statistics would reveal variations on even smaller scales. The temperature variations show no correlation with surface brightness. Because of the overall high temperature of the gas relative to the useful upper limit of the Chandra energy range, the excursions in the map to lower temperatures from the ambient temperature of $``$9 keV are statistically significant, while the excursions to higher temperature are not. Thus, the hot spot to the southwest is not significantly hotter than the surrounding gas at 90% confidence. In contrast, the cold spot corresponding to the central brightness peak is cooler than the surrounding gas at 90% confidence. On the other hand, excursions to higher temperature in regions of lower ambient temperature, such as the extended hot region just inside the northern cold front, are significant. This hot region is coincident with the middle part of the cold front, which has a less well defined edge than the rest of the front. As we will discuss in §6, this feature may be related to the development of a hydrodynamic instability in the front. A detailed spectral analysis of this hot region shows that it is significantly hotter than the gas to either side of it: $`9.9_{1.9}^{+3.4}`$ keV compared to $`6.1_{0.9}^{+1.2}`$ keV elsewhere just inside the cold front for a single-temperature MEKAL model with the absorption and abundance fixed as above. We discuss this hot spot in greater detail in §6

### 4.4. Mass Distribution

Combining the surface brightness and temperature profiles from above, we produce two mass profiles of the cluster in each of the five sectors: a profile of the gas mass, and a profile of the total gravitational mass. The gas mass is determined primarily by the surface brightness profile, which is deprojected under the assumption that the cluster is spherically symmetric. The temperatures from the variable-abundance fits were used to derive the emissivity. The gravitational mass is calculated by deriving the pressure in each concentric annulus and applying the condition of hydrostatic equilibrium to derive the necessary mass interior to the given shell. The gravitational and gas masses are plotted in Figure 8.

As can be seen from the figure, the gas masses determined from the different sectors are quite consistent. This should not be surprising given the similarities between the surface brightness profiles. The gravitational mass profiles are not as well determined as those of the gas mass, for several reasons. The first is the scatter in the temperature measurements, which contribute directly to the pressure calculations. The second is that the calculation of the gravitational mass assumes that the cluster is in hydrostatic equilibrium. Since it is undergoing a merger and shows evidence of relative motion between the gas in different parts of the cluster, such as the cold front in sector 2 (the red points in Figure 8), it is almost certainly not in hydrostatic equilibrium. Not only is the gas unlikely to be at rest relative to the gravitational potential, it has also probably been heated by the merger beyond the temperature one would derive simply by applying the hydrostatic condition to a cluster with the combined mass of the original subclusters (Randall, Sarazin, & Ricker, 2002).

Also shown in Figure 8 is a mass profile derived from the $`\beta `$-model fit produced in §4.1, under the assumption that the cluster is isothermal. Here we used the mean temperature from the spectral models in §4.2 with fixed abundance. This model fits the data remarkably well, both for the gas mass and total mass profiles, particularly at large radii. At radii less than the core radius of the $`\beta `$-model ($`r<290`$ kpc), the fit to the gravitational mass is worse. In this region, however, the pressure gradient should be small, and since our method for determining the gravitational mass from the data involves computing pressure differences between points, the accuracy of mass at these radii is suspect.

Using the outermost bin for which we have estimates of both the gas mass and the gravitational mass, the gas mass fraction is $``$13%. The $`\beta `$-model mass profile gives a similar fraction. This fraction is consistent with the typical value of 10–20% (David et al., 1995; Allen et al., 2002, all for $`H_0=50`$ km s-1 Mpc-1).

## 5. South Excess

As mentioned above, there is a faint, large scale excess of emission to the south and southwest on the cluster periphery. In the ROSAT image of the cluster (Figure 2, this excess appeared to be linear, leading us to believe that it was likely to be a merger shock. As mentioned in §3, the Chandra image of the cluster shows this region to be much more diffuse, and not a sharp feature at all. It appears much more round, and probably extends off the south edge of the detector. The ROSAT image, with its large field-of-view, shows that the south excess does indeed extend about 4′ past the edge of the ACIS-I detector. The excess is slightly more than 8′ across, so approximately half of the south excess is outside the field-of-view of the ACIS-I detector.

We fit the south excess with a two-temperature MEKAL model, with one component’s parameters fixed at the values determined at the same distance from the cluster center and to the northeast of the excess, and the other component’s temperature, abundance and normalization allowed to vary. This second component was found to have a temperature of $`4.9_{0.8}^{+1.2}`$ keV and an abundance of $`<0.73Z_{\mathrm{}}`$. Thus, the south excess is cooler than the gas in the main body of the Abell 2034 cluster. The excess in this region therefore cannot be due to shock compression of this hotter gas.

One possible interpretation is that the south excess is a smaller cluster which is merging with Abell 2034, but which has not yet been shock heated by the interaction. This would explain the low temperature of the south excess. However, the flux from the south excess is much too low for a $``$4.9 keV cluster at $`z=0.113`$, if it follows the standard $`L_X`$$`T`$ relation for clusters (e.g. Markevitch, 1998).

Another possible explanation of the south excess is that it was originally a relatively small, cool, and X-ray faint subcluster which has been shock heated to $`kT=4.9`$ keV as a result of a merger with Abell 2034. There are two possible concerns with this explanation. First, there is no evidence for a strong interaction between the south excess and the main cluster in Abell 2034. Second, the observed temperature of the south excess is lower than might be expected for gas shocked while merging with the main cluster, given the very high temperature of Abell 2034. This is true even for the measured upper limit on the temperature of the south excess. However, it is possible that the merger shock was a weak oblique shock resulting from an offset merger (a merger with a large impact parameter), or that the south excess region originated as substructure within Abell 2034 and did not fall through the full potential well of the cluster. It is also possible that while the subcluster was initially shock heated, it has since had time to adiabatically expand as the center of its dark matter potential has moved ahead of the gas. This would require the subcluster either to have a very small mass or to be a quite advanced merger, probably having already fallen through the cluster. Otherwise, the cumulative effect of ram pressure over the lifetime of the merger would be insufficient to decouple the gas from the potential.

Alternatively, we consider the possibility that the south excess is a cluster at a considerably larger redshift than that of Abell 2034. By requiring that this cluster fit the $`L_X`$$`T`$ relation of Markevitch (1998), it must be at a redshift of $`0.30<z<1.25`$. Unfortunately, there are no strong line features in the X-ray spectrum of the south excess region, so it is not possible to determine the redshift associated with this emission directly from the spectrum. There are only a handful of faint galaxies visible in the Digital Sky Survey 2 (DSS-2) red image in the region of the south excess. This would be consistent with the south excess being associated with a moderately distant cluster. Assuming the brightest cluster galaxy had an absolute magnitude of $`M_V=22`$, it would have an apparent R magnitude $`m_R=19.1`$ at a redshift of $`z=0.3`$, and would be only slightly above the sensitivity limit of the DSS-2. The bulk of the galaxies in the cluster would therefore be below the limit of the DSS-2. We assume a rest frame color $`(VR)=0.9`$ and take $`K`$-corrections from Coleman, Wu, & Weedman 1980. We have obtained moderately deep optical and near infrared images of the field, which will be described in a subsequent paper, along with a more detailed analysis of the X-ray emission from this region. Even deeper optical imaging and follow-up spectroscopy should also be done to verify the existence and proposed redshift of this background cluster.

## 6. Northern Cold Front

To the north, the cluster contains a sharp surface brightness discontinuity which was visible in an earlier ROSAT image (Fig. 2). As with other such features seen prior to Chandra, it was initially assumed to be a merger shock. Also like some other features of its kind, it in fact has turned out to be a “cold front” of the type first discovered in Abell 2142 (Markevitch et al., 2000), although the temperature contrast is small relative to Abell 2142 or other similar cases, and the front itself does not appear to be stable. In addition, the density contrast across the front is lower than in other clusters: it is at most a factor of 2, compared to factors of 2–4 observed in other cold fronts e.g. A2142, Markevitch et al., 2000; A3667, Vikhlinin et al., 2001a. The cold front in Abell 2034 may also be less regular than those observed in other clusters. The temperature contrast across the front is also low: the gas inside the surface brightness edge is not significantly cooler than that outside the edge. We conclude therefore that it is a cold front, since the high surface brightness side of the edge is hot hotter than the low surface brightness side, and therefore it cannot be a merger shock. As can be seen in Figure 9a, there is a region of excess emission ahead of the cold front to the north and northwest of the front. Moreover, the gas inside the center of the cold front is hotter than the rest of the gas inside the front (see Figure 9b), a further indication that the front may be unstable. The leading edge of the cold front is roughly circular within $``$$`50^{}`$ of its axis of symmetry. The axis of symmetry is at a position angle of approximately $`2^{}`$ clockwise from North. The cD galaxy associated with the cold front is ahead of the front, $`12^{}`$ west of the axis of symmetry. The projected direction of motion, then, is probably somewhere between these two values.

We extracted surface brightness and temperature measurements from a set of concentric circular annular wedges which were defined to match the curvature of the cold front as well as possible (see Figure 9a). These annuli were divided into three sectors so that the boundaries of the sectors bracketed the hot region in the center of the front. The projected temperature and surface bright profiles are shown in Figure 9b,c. We then deprojected the surface brightness profile separately for each sector under the assumption of spherical symmetry. The low surface brightness of the gas in these regions serves to minimize projection effects in the temperature profile, so the differences between the projected temperature profile and a deprojected one would be minimal. We therefore used the projected profile for simplicity. The deprojected surface brightnesses were then converted to electron densities, which, when combined with the measured projected temperatures, produced pressure profiles in each of the three sectors (Figure 9d). As mentioned above, the density contrast across the front is less than a factor of two. The pressures are consistent with being continuous across the cold front in each sector, and remain continuous out to large radii. Both within the cold front and at large radii past the cD galaxy, the pressures are the same in each sector to within the errors. Just outside the cold front, however, the pressure in the two sectors containing the excess emission is higher than the pressure in the sector without the excess, although not by more than 2-$`\sigma `$. In all three sectors, the pressure drops with increasing distance from the cluster center.

The overpressure in the region of the excess emission, combined with the surprisingly hot gas in the central sector lead us to believe that the front is unstable. While the overpressure is measured at a fairly low significance, the difference in temperatures between the sectors inside the front is significant at greater than 90% confidence. Therefore, we consider a mechanism for creating an instability in the front. As we discussed in §3, the cD galaxy associated with the cold front has moved ahead of the front, as it has not been slowed by the ram pressure forces that affect the ICM. Since the dark matter is also presumably collisionless, it too has moved ahead of the cold front. This creates a situation in which the gradient of the potential is inverted relative to the density and pressure gradients in the gas that used to be at the center of the potential. Thus, the dense gas which was formerly at the center of the cluster is now further out from the center of the potential than is the less dense gas. In the absence of a gravitational potential inside the front, a Rayleigh-Taylor instability should form along the front (e.g. Jones, Ryu, & Tregillis, 1996). Instabilities should be visible if the timescale for their formation is short compared to the timescale of the merger. The timescale for the formation of this instability is dependent on the acceleration due to the potential. Unfortunately, since the dark matter and the X-ray emitting gas are decoupled, it is not possible to reliably measure the potential, for example using the technique described in Vikhlinin & Markevitch (2002)

A further test of the stability of the cold front comes from a measurement of its specific entropy. In a stable cold front centered on the gravitational potential, the gas in the center of the cold core has a lower specific entropy than the surrounding gas. Therefore, a higher specific entropy in the center of the cold front would indicate that cold core has been disrupted. The specific entropy is defined as

$$\mathrm{\Delta }ss_1s_2=\frac{3}{2}k\mathrm{ln}\left[\left(\frac{T_1}{T_2}\right)\left(\frac{\rho _1}{\rho _2}\right)^{2/3}\right],$$
(1)
where the subscripts denote the two regions of interest. Since this estimate of the entropy is only useful in a qualitative sense, we approximate $`\rho _1/\rho _2`$ as $`(I_{X1}/I_{X2})^{1/2}`$, where $`I_X`$ is the X-ray surface brightness in a given region. We find that the hot region has a higher specific entropy than the surrounding gas, with $`\mathrm{\Delta }s/(\frac{3}{2}k)=0.6\pm 0.2`$. Disruption of the core by Rayleigh-Taylor instability could induce turbulence in the gas, which could explain the heating and consequent increase in entropy.

We should note that hydrodynamic instability is by no means the only possible explanation for the increase in entropy at the center of the front. A hydrodynamic shock, such as a merger shock, increases the specific entropy of the shocked gas. It is possible that a small merger shock has heated the gas at the center of the front and left the rest of it untouched, thereby increasing the entropy only in this central region.

Another possible explanation for the heating is that a now-extinct AGN associated with the cD or one of the other nearby galaxies had heated this gas at some point in the past. We find this explanation to be quite unlikely, however, as there is no evidence for radio activity in any of these galaxies, nor is there any evidence for diffuse emission on the scale of the hot region, even at low frequencies. The radio relic, discussed in the next section, covers a much larger area, including the cold gas to either side of the hot region.

## 7. Radio Relic

Kempner & Sarazin (2001) reported a tentative detection of a radio relic in Abell 2034 near the position of the cold front. The radio emitting electrons that produce relics are also expected to produce X-rays by Inverse Compton (IC) scattering off of Cosmic Microwave Background (CMB) photons. The ratio of radio flux to IC flux is equal to the ratio of the energy density in the magnetic field to that in the CMB. Measuring both the radio synchrotron and X-ray IC emission would allow both the magnetic field strength and the energy density in relativistic electrons to be determined.

To measure the strength of Inverse Compton emission, we have extracted the X-ray spectrum from the region of the radio relic, and fit it with the sum of a MEKAL model and a power-law model, with the photon index of the power-law set to $`1+|\alpha |`$, where $`\alpha `$ is the relic’s spectral index as determined by Kempner & Sarazin (2001). This analysis is complicated by the temperature structure discussed in §6. To account for this, we extracted spectra separately from the regions having different temperatures and fit their temperatures separately while requiring a single power law model for the different regions. The spectra are extremely well fit by a single-temperature MEKAL model, with $`\chi ^2/\mathrm{d}.\mathrm{o}.\mathrm{f}.=0.99`$. The addition of a power-law component improves the fit marginally. The power-law component has a 1–10 keV flux of $`7.5_{6.4}^{+4.9}\times 10^{14}`$ ergs cm-2 s-1, which implies a mean amplitude of the magnetic field of $`0.3<B<0.9`$ $`\mu `$G (90% confidence), including the 1-$`\sigma `$ measurement errors in the radio flux (Kempner & Sarazin, 2001). This is consistent with the magnetic field strengths derived from analyses of radio halos in other clusters (e.g. Feretti & Giovannini, 1996; Govoni et al., 2001). The errors on this measurement are admittedly quite large, but improved calibration of the Chandra detectors at energies below 1 keV could reduce the errors on this measurement in a future re-analysis of the data.

## 8. Summary

Our analysis of the Chandra observation of the massive, moderate redshift cluster Abell 2034 has revealed a number of interesting features. We have determined that the temperature of the cluster is fairly constant out to a radius of at least about 800 kpc. This includes the central region of the cluster, which shows no evidence for the large cooling flow that had previously been claimed. In contrast, the cluster shows some strong signatures of an ongoing merger. We have shown that the surface brightness discontinuity on the NE edge of the cluster consistent with a cold front which is in the process of being disrupted by a combination of ram pressure and the decoupling of the gas from the dark matter potential. This decoupling of the gas from the dark matter has also been demonstrated to have occurred in 1E0657-56 (Markevitch et al., 2002). A large concentration of galaxies, including a cD galaxy, is visible just ahead of the cold front. We suggest that these galaxies, including a cD galaxy, are centered on the potential well of a subcluster, and that the gas in the cold front was the cooling core of this subcluster. The collisionless galaxies and dark matter have now moved ahead of the collisional gas, which has been slowed by ram pressure during the merger. This has caused the front to begin to develop Rayleigh-Taylor instabilities. In addition to the cold front and the lack of a centrally condensed cooling flow, other major evidence for the cluster being out of equilibrium includes the significant offset of the central cD from the center of the X-ray emission.

We have made a weakly constrained detection of IC emission from the region of the proposed radio relic. The magnetic field derived from this limit is consistent with other measurements of cluster magnetic fields using radio halos.

There is a region of excess emission at a cooler temperature (4.9 keV at the redshift of Abell 2034) to the south of the main cluster. There is no clear evidence for any interaction between this gas and the main Abell 2034 cluster. The south excess is too faint to be an undisturbed cluster at the redshift of Abell 2034 or closer, if it follows the normal $`L_X`$$`T`$ relation. It could be the remnant of a shocked subcluster, but this too would be too faint, as shock heating would boost its luminosity by an amount that would keep it more or less consistent with the normal $`L_X`$$`T`$ relation. Instead, we suggest that the south excess is a background cluster at an estimated redshift of $`0.3<z<1.25`$. Much more work is needed to confirm both the existence of and distance to this new cluster, including deep optical imaging and spectroscopy. We have begun this follow-up work, which we will discuss in a future paper.

Support for this work was provided by the National Aeronautics and Space Administration through Chandra Award Numbers GO1-2123X, GO2-3159X, and GO2-3164X issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-39073; and by NASA contract NAS8-39073. We thank John Houck for the S-Lang code which served as the basis for the program to make the temperature map. We also thank Roger Chevalier, William Forman, and D. Mark Whittle for helpful comments on the manuscript.


format: text

text:
# A maximum likelihood approach to the destriping technique

## 1 Introduction

One of the major goals of cosmology is to determine the cosmological parameters which describe the structure and evolution of the Universe. In this respect CMB observations are a powerful tool directly probing the early phases of the Universe. Recent results from the WMAP satellite (Bennett et al. bennett03 (2003)) show that high accuracy in such measurements can be achieved with an optimal choice of observing site (the second Lagrange point of the Sun-Earth system, $`L_2`$, for good thermal and environmental stability), careful instrument design and control of systematic effects. The last point is related to the observing strategy adopted, which should be as redundant as possible, with different measurements of the same sky region with different detectors and on different time scales in order to properly control systematics.

Future space missions like Planck111http://astro.estec.esa.nl/SA-general/Projects/Planck/ which are designed to have a signal-to-noise ratio of the order of few tens (far larger than WMAP), require control of systematic effects at the $`\mu `$K level. In this respect several techniques have been developed to treat systematics, to detect and remove them in the best possible way. Burigana et al. (Burigana99 (1997)), Delabrouille (Delabrouille98 (1998)) and Maino et al. (Maino99 (1999), Maino02 (2002)) have considered, in the context of the Planck mission, a simple destriping algorithm to remove systematics like the $`1/f^\alpha `$ noise. Mennella et al. (Mennela02 (2002)) have instead considered destriping when dealing with periodic fluctuations such as those induced by thermal instabilities.

Destriping methods work on time-ordered data (TOD) and produce TOD cleaned of systematics. When TOD is cleaned it is possible to co-add observations on the same region (pixel) of the sky to obtain a sky map which gives a visual impression of the data. Although it is non-optimal, in the sense that it would not necessarily produce the map with the minimum possible variance as instead provided by the Generalized Least Square solution of the map-making problem (see e.g. Natoli et al. natoli01 (2001)), it provides a fast and accurate map-making algorithm. In addition, the analysis of TOD cleaned of systematics is useful for several applications relevant for the Planck data analysis (e.g. in-flight main beam reconstruction (Burigana et al. Bur02 (2002)) and calibration (Bersanelli et al. Ber97 (1997), Piat et al. Piat03 (2003), Cappellini et al. Cap03 (2003)), time series analysis) and for the scientific exploitation of Planck data (e.g. source variability studies (Terenzi et al. Ter02 (2002))).

In this paper we consider the destriping technique in the light of maximum-likelihood analysis and present a general formulation of the destriping technique. We restrict our analysis to $`1/f`$ noise fluctuations. They produce noise which is strongly correlated in time and, when coupled with the observing strategy, will lead to stripes in the final maps that would alter the signal statistics. This is of extreme importance for the CMB which is expected to be a Gaussian random field.

The basic idea in destriping is to model the noise in the TOD by a linear combination of simple arithmetic functions, such as polynomials or Fourier components. The amplitudes of these base functions are determined taking advantage of the redundancy of the scanning strategy for which the same points on the sky are monitored several times during the mission. In its simplest form destriping involves fitting uniform baselines, i.e. one baseline for each elementary scanning period. In order to improve the accuracy of the method, we present here the possibility of fitting several components (base functions).

The destriping method of Burigana et al. (Burigana99 (1997)) and Maino et al. (Maino99 (1999, 2002)) differs from the destriping method of Delabrouille (Delabrouille98 (1998)) in the weights they assign to different map pixels based on the number of measurements falling on that pixel. Our maximum-likelihood analysis presented in this paper leads to a weighting scheme that differs from both of these. Therefore we compare results obtained from all these three methods.

The paper is organized as follows. In Sect. 2 we present the maximum-likelihood approach to the destriping technique, in Sect. 3 we apply it in the case of uniform baselines, and in Sect. 4 we generalize the discussion to arbitrary base functions. We present our conclusions in Sect. 5.

## 2 Destriping – maximum likelihood approach

### 2.1 Maximum likelihood analysis

In the following we present a maximum-likelihood based approach to the destriping problem. We assume that data produced by a generic detector at a given time $`t`$ could be written as:

$$y_t=P_{tp}m_p+n_{t,\mathrm{corr}}+n_t$$
(1)
where $`m_p`$ is the sky signal, assumed to be pixelized, $`P_{tp}`$ is the pointing matrix, $`p`$ is the pixel index, $`n_{t,\mathrm{corr}}`$ is the correlated noise component while $`n_t`$ is the white noise component. The variance of the white noise component is represented by a diagonal matrix $`𝐂_n`$ in the time domain. Eq.(1) could be written in vector form as:

$$𝐲=\mathrm{𝐏𝐦}+𝐧_{\mathrm{corr}}+𝐧.$$
(2)
We model the correlated noise component of the TOD as follows. The TOD is divided in to elementary scanning periods, which we shall here call “rings” (as appropriate for the Planck scanning strategy). For each ring $`j`$ we define a constant offset $`a_j`$, so that

$$y_t=P_{tp}m_p+F_{tj}a_j+n_t.$$
(3)
Here $`F_{tj}`$ equals unity if point $`t`$ lies on ring $`j`$. We write this in matrix form as

$$𝐲=\mathrm{𝐏𝐦}+\mathrm{𝐅𝐚}+𝐧.$$
(4)
We treat both the map and the correlated noise component as deterministic. With these assumptions, we obtain the likelihood function

$$\chi ^2=(𝐲\mathrm{𝐅𝐚}\mathrm{𝐏𝐦})^T𝐂_n^1(𝐲\mathrm{𝐅𝐚}\mathrm{𝐏𝐦}).$$
(5)
If $`N_t`$ is the length (the number of samples) of the TOD stream, $`N_{\mathrm{pix}}`$ is the number of pixels in the map, and $`N_a`$ is the number of unknown amplitudes, then the sizes of the matrices are: $`[F]=(N_t,N_a)`$, $`[P]=(N_t,N_{\mathrm{pix}})`$, $`[𝐂_n]=(N_t,N_t)`$.

We now want to find the maximum likelihood solution for $`𝐚`$. We need to minimize the function in Eq. (5) with respect to both of the unknown variables $`𝐦`$ and $`𝐚`$. First we find the minimum with respect to $`𝐦`$,

$$_m\chi ^2=2𝐏^T𝐂_n^1(𝐲\mathrm{𝐅𝐚}\mathrm{𝐏𝐦})=0.$$
(6)
From this we can solve the map $`𝐦`$,

$$𝐦=(𝐏^T𝐂_n^1𝐏)^1𝐏^T𝐂_n^1(𝐲\mathrm{𝐅𝐚}).$$
(7)
Substituting Eq. (7) back into Eq. (5) we obtain

$$\chi ^2=(𝐲\mathrm{𝐅𝐚})^T𝐙^T𝐂_n^1𝐙(𝐲\mathrm{𝐅𝐚}),$$
(8)
where

$$𝐙=𝐈𝐏(𝐏^T𝐂_n^1𝐏)^1𝐏^T𝐂_n^1.$$
(9)
Here $`𝐈`$ denotes the unit matrix.

The next step would be to minimize $`\chi ^2`$ with respect to $`𝐚`$,

$$_a\chi ^2=2𝐅^T𝐂_n^1𝐙(𝐲\mathrm{𝐅𝐚})=0.$$
(10)
The minimum is given by

$$𝐅^T𝐂_n^1\mathrm{𝐙𝐅𝐚}=𝐅^T𝐂_n^1\mathrm{𝐙𝐲}.$$
(11)
Here we have used the property $`𝐙^T𝐂_n^1𝐙=𝐂_n^1𝐙`$.

We assume from now on $`𝐂_n=diag(\sigma ^2)`$. With this simplification, the minimum of (8) is given by

$$𝐅^T\mathrm{𝐙𝐅𝐚}=𝐅^T\mathrm{𝐙𝐲}.$$
(12)
where

$$𝐙=𝐈𝐏(𝐏^T𝐏)^1𝐏^T.$$
(13)
The effect of $`𝐙`$ acting on a TOD is to subtract from each sample the average of all samples hitting the same pixel. The solution to Eq. (12) is not unique. We may add an arbitrary constant to $`𝐚`$ without changing the value of $`\chi ^2`$. This is equivalent to varying the monopole component of the CMB map, and is irrelevant for anisotropy measurements. To remove this ambiguity we require that the sum of baselines is zero, $`𝐚^T\mathrm{𝟏}=0`$. Here $`\mathrm{𝟏}`$ is a vector with all elements equal to 1. This is equivalent to adding term $`\mathrm{𝟏𝟏}^T𝐚`$ to the left-hand side of Eq.(12).

The solution is now given by

$$𝐚=[𝐅^T\mathrm{𝐙𝐅}+\mathrm{𝟏𝟏}^T]^1𝐅^T\mathrm{𝐙𝐲}.$$
(14)
Matrix $`𝐅^T\mathrm{𝐙𝐅}+\mathrm{𝟏𝟏}^T`$, unlike $`𝐅^T\mathrm{𝐙𝐅}`$, is non-singular, provided that there are enough intersection points between the rings. $`\mathrm{𝟏𝟏}^T`$ denotes a matrix with all elements equal to one. When the amplitude vector $`𝐚`$ has been determined, the CMB map can be computed according to Eq. (7). These equations are the main theoretical result of this paper.

### 2.2 Pointing and beam shape

The matrix $`P`$ spreads the map $`𝐦`$ into TOD. In principle, the beam shape and profile can be incorporated in $`𝐏`$. The elements of $`𝐏`$ then determine the weights that different pixels contribute to a given measurement. In this way beams with arbitrary shapes and profiles can be treated. In this case the matrix $`𝐏^T𝐏`$ is non-diagonal. Due to its large size, in Planck type of missions its inversion would present a major computational burden.

A simpler approach is to consider each sample in the TOD to represent the temperature of the pixel at the center of the beam. Then $`𝐏`$ takes a particularly simple form, consisting of ones and zeros for a full-power measurement like Planck. Each row contains one non-zero element identifying the pixel on which the corresponding measurement falls. Matrix $`𝐏^T𝐏`$ becomes diagonal, the diagonal elements giving the number of hits on each pixel. We follow this approach from here on.

### 2.3 Comparison with earlier work

Equation(8) can be put into form

$$\chi _n^2=\underset{p}{}\frac{_{ik,jl}(a_ia_jy_{ik}+y_{jl})^2d_{ik}^pd_{jl}^p}{2\sigma ^2_{ik}d_{ik}^p}$$
(15)
where index $`p`$ labels pixels, $`i,j`$ scanning rings, and $`k,l`$ samples on a given ring. A combined index $`ik`$ or $`jl`$ identifies a measurement. We define the symbol $`d_{ik}^p`$ so that $`d_{ik}^p=1`$ if measurement $`ik`$ falls into pixel $`p`$, otherwise $`d_{ik}^p=0`$. Due to the factors $`d_{ik}^p`$ and $`d_{jl}^p`$ in the numerator of Eq. (15) only those pixels $`p`$ contribute to the pixel sum which lie on two or more scanning rings, and the sum $`_{ik,jl}`$ is equal to 2 times the sum over all pairs of measurements falling on pixel $`p`$.

Eq. (15) can be compared to Eq. (2) of Maino et al. (Maino99 (1999)) or to Eq. (10) of Burigana et al. (Burigana99 (1997)). The formulae differ in that in Eq. (15) $`\chi _n^2`$ has in the denominator the term $`n_p=_{ik}d_{ik}^p`$, which gives the total number of hits in pixel $`p`$.

Delabrouille (Delabrouille98 (1998)) gives the general form

$$S=\underset{p\mathrm{sky}}{}\underset{\mathrm{pairs}}{}w(p,ik,jl)(y_{ik}y_{jl}a_i+a_j)^2$$
(16)
for the function $`S`$ to be minimized. Here the second sum refers to all pairs $`(ik,jl)`$ that can be formed of the measurements falling onto pixel $`p`$, and $`w`$ is a weight function to be chosen. Based on the fact that pixel $`p`$ contributes $`n_p(n_p1)/2`$ pairs, Delabrouille suggests choosing $`w1/(n_p1)`$. The result of Maino et al. (Maino99 (1999)) corresponds to $`w=\text{const.}`$ and our new result, Eq. (15), to $`w=1/n_p`$.

### 2.4 Circular scanning

In the nominal scanning strategy of Planck the spin axis follows the ecliptic plane. The spin axis is kept anti-solar by repointing it by $`2.5^{}`$ every hour. The spacecraft rotates around the spin axis at a rate of 1 rpm. During one hour Planck scans the same circle on the sky 60 times. As the sky signal is almost time-independent, the data can be coadded to reduce the length of the data stream by a factor of 60. We call this set of 60 circles, and the corresponding segment of the coadded TOD a “ring”. The crossing points of the rings are important calibration points, which allow for the removal of the correlated noise component from the TOD.

The opening angle of the scanning circle varies between 80-90 degrees, depending on the location of the detector on the focal plane. The sampling frequency for the 100 GHz LFI receiver is 108.3 Hz. The instrument then collects 6498 temperature values, or “samples”, at each rotation of the spacecraft, corresponding roughly to a $`3^{}`$ shift between successive samples. A total of 8766 rings builds up one year of observations.

In reality, the angular velocity of the rotation does not stay exactly constant, especially immediately after repointing, and the samples from different circles of the same ring are shifted in position. This will probably require discarding the first few circles of each ring, and resampling or phase binning the rest before performing the coaddition. During the ring the spin axis of the satellite will follow a nutation ellipse with maximum amplitude of 1.5 at the end of the nutation damping phase (van Leeuwen et al. vanL02 (2002)). This may degrade slightly the performance of destriping. We ignore these complications in this paper.

The destriping technique applies particularly well to a Planck-like measurement pattern resulting from the coadding of scanning circles into rings which breaks the stationarity of the data. The $`1/f^\alpha `$ noise component in the coadded TOD is well presented by a piecewise defined function, where each piece consists of a linear combination of a few base functions.

## 3 Uniform baselines

### 3.1 Simulation results

We have carried out simulations of the Planck LFI 100 GHz detector. The underlying CMB map was created by the Synfast code of the HEALPix package222http://www.eso.org/science/healpix, starting from the CMB anisotropy angular power spectrum computed with the CMBFAST code333http://physics.nyu.edu/matiasz/CMBFAST/cmbfast.html (see Seljak & Zaldarriaga cmbfast96 (1996), and references therein) using the cosmological parameters $`\mathrm{\Omega }_{\mathrm{tot}}=1.00`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`\mathrm{\Omega }_\mathrm{b}h^2=0.02`$, $`h=0.7`$, $`n=1.00`$, and $`\tau =0.0`$. We created the input map with HEALPix resolution $`N_{\mathrm{side}}=1024`$ and with a symmetric Gaussian beam with a full width at half maximum (FWHM) of $`10^{}`$. We then formed the signal TOD by picking temperatures from this map. All our output maps have the resolution parameter $`N_{\mathrm{side}}`$=512, corresponding to an angular resolution $`7^{}`$.

The scanning pattern corresponds to the nominal Planck scanning strategy of the 100 GHz LFI detector number 10.444Simulation Software is part of the Level S of the Planck DPCs and is available for Planck collaboration at http://planck.mpa-garching.mpg.de. The angle between the satellite spin axis and the optical axis of the telescope is $`85\mathrm{°}`$. The beam center is pointing towards $`(\theta ,\varphi )=(3\stackrel{}{.}737,126\stackrel{}{.}228)`$. Here $`\theta `$ is the angle from the optical axis and $`\varphi `$ is an angle counted clockwise from the axis pointing from the center of the focal plane towards the satellite spin axis. We assumed no spin axis precession. Our simulated data set consists of 5040 scanning rings, corresponding to 7 months of measurement time. The TOD stream contains 6498 samples on each ring, corresponding to a sampling frequency of $`f_\mathrm{s}=108.3`$ Hz. The sky coverage is 98.5%.

In our simulations we have assumed a symmetric Gaussian beam, and convolved the input map with the beam.

The rings cross at points which are mostly concentrated near the ecliptic poles. We count as crossing points all points where two measurements on different rings fall on the same pixel. Since our pixel resolution ($`7^{}`$) exceeds the repointing angle of the spin axis ($`2.5^{}`$), the crossing points include cases where two nearby rings pass parallel through the same pixel, without actually crossing each other.

We used the Stochastic Differential Equation (SDE) technique to create the instrument noise stream, which we added to the signal TOD.555SDE is one of the two methods in the Planck Level S pipeline for producing simulated instrument noise. The method was implemented by B. Wandelt and K. Górski and modified by E. Keihänen. We generated noise with the power spectrum

$$P(f)=\left(1+\frac{f_\mathrm{k}}{f}\right)\frac{\sigma ^2}{f_\mathrm{s}},(f>f_{\mathrm{min}})$$
(17)
with parameters $`\sigma =4800\mu \text{K}`$ (CMB temperature scale), $`f_\mathrm{k}=0.1\mathrm{Hz}`$, and $`f_{\mathrm{min}}=10^6\text{Hz}`$. The parameter $`f_\mathrm{k}`$ is called the knee frequency. The noise level $`4800\mu \text{K}`$ corresponds to the estimated white noise level of one $`100\text{GHz}`$ LFI detector.

Fig. 1 shows a five-hour section of a coadded noise TOD and the baselines fitted to it. Fig. 2 shows the baseline distribution from a set of 10 simulated 7-month noise TODs.

There was no foreground included in the simulations presented in this paper, but we have also verified our destriping codes on simulated data sets with foreground. We found that the foreground has an insignificant impact on the baselines determined by the destriping, in agreement with the discussion in Maino et al. (Maino02 (2002)). The quality of destriping is also almost independent of the impact of another class of instrumental systematic effects, main beam distortions and straylight, as the temperature differences at crossing pixels are dominated by the noise and only minimally affected by the spurious signals ($`\mu `$K) introduced by optical distortions (see e.g. Burigana et al. 2001).

Fig. 3 shows the input $`C_{\mathrm{}}`$ spectrum and the spectrum derived from the simulated TOD after destriping. We used the Anafast code of the HEALPix package to compute the $`C_{\mathrm{}}`$ spectrum of the destriped map. We subtracted from the derived spectrum an estimate of the noise level $`C_{\mathrm{noise}}=0.197\mu \mathrm{K}^2`$ (estimated as the average of $`C_{\mathrm{}}`$ over $`\mathrm{}=980\mathrm{}1000`$) and corrected the spectrum for the beam shape convolution and pixel convolution. The angular spectrum shown is thus $`\stackrel{~}{C}_{\mathrm{}}=(C_{\mathrm{}}C_{\mathrm{noise}})/(B_{\mathrm{}}^2h_{\mathrm{}}^2)`$, where $`B_{\mathrm{}}=\mathrm{exp}(\sigma _b^2\mathrm{}(\mathrm{}+1)/2)`$, with $`\sigma _b=10^{}/\sqrt{8\mathrm{ln}(2)}`$, is the beam convolution function corresponding to the assumed beam width of $`10^{}`$ (FWHM), and $`h_{\mathrm{}}`$ is the pixel convolution function (provided by the HEALPix package).

Fig. 4 shows the same for a noise level reduced by the factor $`\sqrt{24}`$, corresponding to the combination of 24 detectors. Here the subtracted noise level was $`C_{\mathrm{noise}}=0.197\mu \mathrm{K}^2/24=0.0082\mu \mathrm{K}^2`$.

Note that Figs. 3 and 4 are for illustration only, as this paper does not address the full CMB $`C_{\mathrm{}}`$ estimation problem, and thus we have just used the above crude estimate for $`C_{\mathrm{noise}}`$.

### 3.2 Comparison of different weighting schemes: maps and angular power spectra

We have written a destriping code which allows us to compare the different weight functions discussed in Sect. 2.3. We use the Cholesky decomposition technique to solve the set of linear equations.

We have chosen the root-mean-square (rms) value of the residual noise map (see below) pixels as a figure of merit that we use to compare different destriping methods. This map rms is related to the $`C_{\mathrm{}}`$ spectrum through the relation

$$\mathrm{rms}^2=\underset{\mathrm{}}{}\frac{2\mathrm{}+1}{4\pi f_{\mathrm{sky}}}C_{\mathrm{}}.$$
(18)
The map rms squared is thus a weighted sum of the angular spectrum, with high weight on high multipoles. The sky coverage $`f_{\mathrm{sky}}=0.985`$ enters here because we have computed our rms values over the visited pixels only. When computing $`C_{\mathrm{}}`$ spectra, we have set $`T=0`$ in the remaining pixels.

We compute the residual noise map by taking the difference between the destriped map and the noise-free reference map and subtracting the monopole component. Note that because of the incomplete sky coverage, removing the monopole affects the $`C_{\mathrm{}}`$ spectrum at all $`\mathrm{}`$ (not only $`\mathrm{}=0`$). The reference map is computed by coadding the pure signal TOD into a map of resolution $`N_{\mathrm{side}}`$=512. The expected contribution from white noise to the residual map rms is 220.95 $`\mu `$K.

While the rms of the residual noise map is a natural measure of the CMB map quality, the main scientific interest is perhaps not in the CMB map itself, but rather in its angular power spectrum $`C_{\mathrm{}}`$. It is therefore of interest to see the impact of the destriping methods on the different parts of the $`C_{\mathrm{}}`$ spectrum. The map rms is dominated by the high $`\mathrm{}`$ part, and does not reveal the difference in performance of the various methods in the low $`\mathrm{}`$ part. Thus we have computed the angular power spectra $`C_{\mathrm{}}`$ of the residual noise maps.

Because of the random nature of the noise, the result of a comparison between methods may vary from one noise realization to another. We therefore performed destriping 10 times, with different realizations of instrument noise. The underlying CMB map was kept the same. We then used the average of the 10 residual map rms values and the residual noise map $`C_{\mathrm{}}`$ spectra to compare the methods. We also calculated the standard deviation (std) of the 10 map rms and $`C_{\mathrm{}}`$ values, to see whether the differences between the methods were statistically significant. Thus this average rms approximates the expectation value for the rms with an accuracy of $`\text{std}/\sqrt{10}`$. However, the std itself tells us how much we can expect the residual map rms for a single realization to deviate from this average value.

Table 1 presents a comparison between different weight functions discussed in Sect. 2.3. The corresponding $`C_{\mathrm{}}`$ spectra are shown in Fig. 5. Since from our simulations we have the noise streams available separately, we also computed reference baselines as the average of the noise stream over each ring. This reference baseline can be thought of as the “true” baseline of the noise. For comparison, we then subtracted the reference baselines from the TOD and computed the rms of the resulting map. This represents an ideal situation, where we could determine the baselines exactly. The reference rms is given on the last line of Table 1. Actual residual noise rms values are always larger.

Weight functions $`w=1/n_p`$ and $`w=1/(n_p1)`$ in Eq. (16) give similar results, due to the fact that for most pixels $`n_p1`$. Weight $`w=1/n_p`$ suggested by our maximum likelihood analysis is clearly superior to $`w=1`$. However, the weight $`w=1/(n_p1)`$ gives an even smaller rms, although the difference is very small. The difference of the rms between $`w=1`$ and $`w=1/n_p`$ is significant because the difference is about 10 times larger than the respective std of the rms. Although the difference between $`w=1/n_p`$ and $`w=1/(n_p1)`$ is much less than the std between different noise realizations, it was in the same direction in each realization. Note that we were able to measure this small difference only because we used the same set of random seeds for all weighting schemes.

Since the difference between the weight functions $`w=1/n_p`$ and $`w=1/(n_p1)`$ is so small that it does not show up in Fig. 5, we show just the differences in Fig. 6.

It is well known that maximum-likelihood analysis should provide the minimum-variance solution. Therefore it may at first seem surprising that the maximum-likelihood based weight function did not give the best results. However, the maximum-likelihood solution is the optimal one only if the model used corresponds to reality. Here we have modelled the $`1/f`$ noise component in the TOD by a uniform baseline, which is a simplified model. Further, we have assumed that the baselines are independent from ring to ring. The reason for the maximum-likelihood solution not giving the best result is that the noise model used in the analysis does not exactly correspond to the actual noise properties.

To verify this, we re-generated our input noise in a way that better corresponds to the model assumed in the analysis. We generated the $`1/f`$ noise in the usual way, but then, for each ring, we took the average over the ring, and replaced the original $`1/f`$ contribution to the ring with this average value, on top of which we added white noise. This way we obtained noise which still has a realistic correlation between scanning rings, but consists of baselines + white noise only. The results are shown in Table 2. Now the maximum-likelihood based weight function gives the smallest variance map.

Thus it seems that the slightly better performance of the Delabrouille weighting scheme ($`w=1/(n_p1)`$) is related to the effect of that part of the correlated noise which deviates from uniform baselines.

The scatter in the individual residual noise map $`C_l`$ values from one noise realization to another was larger than the difference between the methods. (See Fig. 7 for the $`C_{\mathrm{}}`$ spectra from the first three noise realizations, using our weighting scheme, $`w=1/n_p`$.) The difference between $`w=1`$ and $`w=1/n_p`$ becomes however statistically significant when the $`C_{\mathrm{}}`$ are binned into larger bins.

## 4 Increasing the number of base functions

As shown by Delabrouille (Delabrouille98 (1998)) and Maino et al. (Maino99 (1999, 2002)), a simple model with uniform baselines and white noise models quite well the coadded noise stream. One can try to further improve the performance of destriping by introducing more base functions. Delabrouille (Delabrouille98 (1998)) found that the addition of a small number of base functions improved the performance of destriping, whereas Maino et al. (Maino02 (2002)) found no significant improvement. The difference between these results could be due to the different noise spectra considered: $`1/f`$ noise for Maino et al. (Maino02 (2002)) and $`1/f^2`$ for Delabrouille (Delabrouille98 (1998)).

We generalize the discussion of Sect. 2 to include arbitrary base functions. We model the correlated noise component of the TOD by a linear combination of base functions as $`𝐧_{\mathrm{corr}}=\mathrm{𝐅𝐚}`$. Here $`𝐅`$ is a matrix, whose columns contain the base functions, and $`𝐚`$ is a vector containing their (unknown) amplitudes. It is convenient to select an orthogonal set of the base functions, so that $`𝐅^T𝐅`$ is diagonal. Eqs. (5)–(13) hold as such for general baselines.

We have studied two sets of base functions: Fourier components and Legendre polynomials. Both form an orthogonal set.

### 4.1 Simulation results

The solution of the general destriping problem involves the solution of a large linear system of equations

$$𝐅^T\mathrm{𝐙𝐅𝐚}=𝐅^T\mathrm{𝐙𝐲}.$$
(19)
Matrix $`𝐀𝐅^T\mathrm{𝐙𝐅}`$ becomes very large if several base functions are fitted. We use the iterative conjugate gradient method (see, e.g., Press et al. NumRes (1992)) to solve the system. The conjugate gradient method only requires multiplication by matrix $`𝐀`$. That can easily be done algorithmically, without actually storing the full matrix $`𝐀`$ at any one time. The conjugate gradient method allows us to perform the destriping in a relatively small memory space.

Note that we do not need to add the term $`\mathrm{𝟏𝟏}^T`$, since the conjugate gradient method has the property that, when solving system $`\mathrm{𝐀𝐱}=𝐛`$, it automatically finds the solution for which $`𝐱_0^T𝐱=0`$, if iteration is started with $`𝐱=0`$, and $`\mathrm{𝐀𝐱}_0=0`$ and $`𝐛^T𝐱_0=0`$. (The amplitude vector which gives unit amplitudes to the uniform baselines and zero amplitudes to other baselines is an example of such a vector $`𝐱_0`$, so the average of the uniform baselines is set to zero. The case of other possible vectors in the null space of $`𝐀`$ is discussed further below.)

We compare four sets of functions:
1. (“Un.”) Uniform baselines only.
2. (“F1”) Uniform baseline + first Fourier frequency, which gives three components for each ring: constant, $`\mathrm{sin}(2\pi f_{\mathrm{sc}}t)`$, and $`\mathrm{cos}(2\pi f_{\mathrm{sc}}t)`$ ($`q=0,1,2`$). Here $`f_{\mathrm{sc}}=1/`$(60 s) is the scanning frequency.
3. (“L1”) Uniform baseline + 1st (linear) Legendre polynomial.
4. (“L2”) Uniform baseline + 1st and 2nd Legendre polynomials.

The first Legendre polynomials are $`L(x,0)=1`$, $`L(x,1)=x`$, and $`L(x,2)=\frac{1}{2}(3x^21)`$, for $`x[1,1]`$.

Again we averaged the residual noise $`C_{\mathrm{}}`$ spectra and the residual noise map rms over 10 noise realizations.

Table 3 and Fig. 8 present our results of fitting several base functions. We find no improvement in the map rms.

The last column of Table 3 gives a reference rms, which was computed as follows. We defined the reference amplitude vector as $`𝐚_0=(𝐅^T𝐅)^1𝐅^T𝐧`$, where $`𝐧`$ is the pure noise stream (i.e. $`𝐧_{\mathrm{corr}}+𝐧`$ of Eq. (4). We coadded a map from the TOD stream, from which we had removed the reference base functions, $`𝐲F𝐚_0`$, and computed the rms of the residual noise map. Fourier components give the lowest reference rms, showing that Fourier components model the noise best of our base function sets. However, when we fit the data, Fourier components give the poorest results, as the first column of Table 3 shows. The worst of the 10 runs gave a rms of 452 $`\mu `$K. The code took a long time to converge, and the final maps contained a very strong and obviously unphysical dipole-like structure. (The large std in the case of Fourier components does not mean that for some realizations the Fourier components would have given a lower rms than the other methods. Rather, the distribution was just very skew; a small number of realizations gave a very large rms, but all 10 realizations gave an rms larger than the average for any of the other methods.)

The problem with Fourier components is related to small eigenvalues of matrix $`𝐀`$. If there exists a map $`𝐦^{}`$ such that

$$\mathrm{𝐏𝐦}^{}=\mathrm{𝐅𝐚}^{}$$
(20)
for some $`𝐚^{}`$, then $`𝐚^{}`$ is an eigenvector of $`𝐀`$ with zero eigenvalue. If $`𝐚`$ is a solution of Eq. (19), then $`𝐚+𝐚^{}`$ is also a solution. In other words, if we can produce the same TOD stream both as a combination of the base functions, and by picking it from a map with our scanning strategy, then it is not possible, without further information, to tell if this TOD component comes from noise or signal. In practice zero eigenvalues are unlikely to happen, but already small but non-vanishing eigenvalues cause similar problems. Eq. (20) then holds approximatively.

The difficulty of fitting Fourier components is related to the symmetry of the nominal scanning strategy of Planck. Suppose the scanning rings are regular circles with centers on the ecliptic plane. If we give all cosine and sine components equal amplitudes $`a`$ and $`b`$ such that the resulting function $`a\mathrm{cos}\varphi +b\mathrm{sin}\varphi `$ ($`\varphi `$ runs from $`0`$ to $`2\pi `$ around the ring) is at maximum at the northernmost point and at minimum at the southernmost point (or vice versa), and coadd a map from this TOD stream, we obtain a meridian symmetric map for which relation (20) holds. This noise component cannot be resolved without further information on noise properties. The same problem also affects fitting 2nd order Legendre polynomials, albeit less severely. We would expect that other, less symmetric, scanning strategies would not be as prone to such problems.

We discuss this problem quantitatively in the following.

### 4.2 Eigenvalue analysis of the destriping problem

Consider the solution of Eq. (19) in light of eigenvalue analysis.

Assume the TOD stream is of the form

$$𝐲=\mathrm{𝐏𝐦}_\mathrm{𝟎}+\mathrm{𝐅𝐚}_0+𝐧$$
(21)
where $`𝐦_0`$ is the actual map, $`𝐚_0`$ is the “true” baseline vector, and $`𝐧`$ represents the remaining noise component. With these assumptions, Eq. (19) becomes

$$𝐅^T\mathrm{𝐙𝐅𝐚}=𝐅^T\mathrm{𝐙𝐅𝐚}_0+𝐅^T\mathrm{𝐙𝐧}.$$
(22)
Let $`\lambda _i`$ and $`𝐮_𝐢`$ be the eigenvalues and corresponding orthonormal eigenvectors of matrix $`𝐀=𝐅^T\mathrm{𝐙𝐅}`$. Because $`𝐀`$ is symmetric and non-negative definite, all eigenvalues are positive or zero, $`\lambda _i0`$, and the eigenvectors form a complete orthogonal basis. The matrix can be presented by its eigenvectors as

$$𝐅^T\mathrm{𝐙𝐅}=\underset{i}{}\lambda _i𝐮_i𝐮_i^T.$$
(23)
We can expand $`𝐚=_ia_i𝐮_i`$ and $`𝐚_0=_ia_i^0𝐮_i`$ and

$$𝐅^T\mathrm{𝐙𝐧}=\underset{i}{}c_i𝐮_i.$$
(24)
We substitute these expansions into (22) to find

$$a_i=a_i^0+c_i/\lambda _i.$$
(25)
Here it must be understood that this relation holds for components for which $`\lambda _i>0`$. For $`\lambda _i=0`$ one can easily show that $`c_i=0`$ and $`a_i`$ remains undefined.

Consider then the statistical properties of coefficients $`c_i`$. Since $`c_i=𝐮_i^T𝐅^T\mathrm{𝐙𝐧}`$, we find, assuming that $`𝐧`$ is white and $`\mathrm{𝐧𝐧}^T=diag(\sigma ^2)`$, that $`c_i=0`$ and

$$c_ic_k=\sigma ^2𝐮_i^T𝐅^T\mathrm{𝐙𝐅𝐮}_k=\sigma ^2\lambda _i\delta _{ik}.$$
(26)
Looking back at Eq. (25) we observe that

$$(a_ka_k^0)(a_ia_i^0)=\frac{\sigma ^2}{\lambda _i}\delta _{ik}.$$
(27)
We see that if one of the eigenvalues is very small, then the inaccuracy in the corresponding component $`a_i`$ is very large. Actually, the vanishing eigenvalues do not pose a problem, since the conjugate gradient algorithm always sets the corresponding amplitude to zero. Problems are caused by moderately small, but non-vanishing, eigenvalues. How small a value must be regarded as zero depends on the floating point accuracy of the computer and on the convergence criterion one has chosen for the conjugate gradient algorithm.

We have seen that fitting uniform baselines only already gives good results.There is no advantage in trying to fit additional components which have a large inaccuracy. Fitting poorly determined components causes more error than leaving them out entirely. We therefore aim to fit only components that correspond to a large eigenvalue, and eliminate small eigenvalue components. In the following we present a practical method to do this. The method presented does not require full determination of eigenvalues or eigenvectors of matrix $`𝐀`$, which would be a computationally expensive task.

### 4.3 A practical method

Consider the following equation:

$$(𝐅^T\mathrm{𝐙𝐅}+ϵ𝐅^T𝐅)𝐚=𝐅^T\mathrm{𝐙𝐲}$$
(28)
where $`ϵ`$ is a small positive constant. An eigenvalue analysis similar to that presented above shows that the solution of Eq. (28) is related to the solution of (19) through

$$a_i^{}=a_i\frac{\lambda _i}{\lambda _i+ϵ\lambda _{max}},$$
(29)
where $`\lambda _{max}`$ is the largest possible eigenvalue of $`𝐀`$. With our chosen normalization $`𝐅^T𝐅=diag(n_b)`$ it is equal to the number of samples on a ring, $`\lambda _{max}=n_b`$.

The effect of the $`ϵ`$ term in Eq. (28) is to wash out components with eigenvalues smaller than $`ϵ\lambda _{max}`$, while the large eigenvalue component remains unaffected, as long as $`ϵ`$ is small. At the limit $`ϵ0`$ the solution of Eq. (28) approaches that of Eq. (19).

We have repeated our computations with this method. Table 4 presents our results for different values of $`ϵ`$. The results were again obtained using 10 different realizations of the instrument noise, over which the average and the standard deviation were calculated. We see that, with values $`10^6ϵ10^3`$, the accuracy of fitting Fourier components is strongly improved with respect to the $`ϵ=0`$ case. Also the required computation time is reduced. For 2nd order Legendre polynomials we also find a clear improvement.

However, the results for multiple base functions are still worse than for uniform baselines only.

Uniform baselines and 1st order Legendre polynomials, which exhibited no problems with $`ϵ=0`$, are unaffected with $`ϵ=10^4`$ or less, but with $`ϵ=10^3`$ or larger the accuracy of fitting them begins to deteriorate.

Fig. 9 shows the residual noise $`C_{\mathrm{}}`$ spectra for the $`ϵ=10^4`$ case.

Depending on the number of base functions, the code took 3-5 seconds per iteration step on one processor of an IBM eServer Cluster 1600 computer. The total computation time varied between 2 and 30 minutes.

To illustrate the use of several base functions we show in Fig. 10 the same 5 hour coadded noise TOD as in Fig. 1, but now with different sets of base functions. Note that the deviations from uniform baselines are exaggerated in this figure. The actual amplitudes of the other base functions are much smaller (by about a factor of 20) than those of the uniform components.

In order to check how the results depend on the knee frequency, we repeated our computations with rescaled noise. We took the $`1/f`$ noise stream, which was originally generated with $`f_\mathrm{k}=0.1`$ Hz, scaled it by a factor 0.5 or 2, and added white noise with the same variance in all cases. This is equivalent to changing the knee frequency by a factor of 0.25 or 4. We thus have results for three knee frequencies: $`f_k=0.025`$Hz, $`f_\mathrm{k}=0.1`$Hz, and $`f_\mathrm{k}=0.4`$Hz.

The obtained residual map rms for $`f_k=0.4`$Hz and $`f_\mathrm{k}=0.025`$Hz are shown in Tables 5 and 6, for different values of $`ϵ`$. The optimal value of $`ϵ`$ seems to depend somewhat on knee frequency, being smaller at higher knee frequencies.

The $`C_{\mathrm{}}`$ spectra for $`ϵ=10^4`$ for knee frequencies $`f_\mathrm{k}=0.1\text{Hz}`$, $`f_\mathrm{k}=0.4\text{Hz}`$, and $`f_\mathrm{k}=0.025\text{Hz}`$ are shown in Figs. 9, 11, 12, respectively. The std of the $`C_{\mathrm{}}`$ for the $`f_\mathrm{k}=0.1\text{Hz}`$ case are shown in Fig. 14.

Looking at the average $`C_{\mathrm{}}`$ spectra of residual noise and their std we see that fitting additional base functions decreases the accuracy of destriping at low $`\mathrm{}`$. However the situation for the map rms values in Table 5 seems more complicated. For $`f_\mathrm{k}=0.1\text{Hz}`$ or smaller, fitting additional base functions does not improve the performance of destriping, but with $`f_\mathrm{k}=0.4\text{Hz}`$ fitting one or two Legendre polynomials besides the constant baselines decreases the map rms, while Fourier components still give inferior results. The improvement in the map rms for $`f_\mathrm{k}=0.4\text{Hz}`$, when fitting Legendre polynomials, comes from the high multipoles. This can be seen from Fig. 13, where we plot the difference between the residual noise $`C_{\mathrm{}}`$ obtained when fitting Legendre polynomials or Fourier components, and when fitting uniform baselines only.

We see that increasing the number of base functions improves the high $`\mathrm{}`$ but worsens the low $`\mathrm{}`$ part of the $`C_{\mathrm{}}`$ spectra. This is true both for Fourier components and for Legendre polynomials. This trend persists for lower $`f_\mathrm{k}`$, but the value of $`\mathrm{}`$ above which we get an improvement goes up and the improvement for those $`\mathrm{}`$ becomes smaller.

Delabrouille (Delabrouille98 (1998)) obtained improved results by fitting several base functions already with $`f_\mathrm{k}=0.1\text{Hz}`$. The difference between our results and his is probably due to differences in the noise model. While we assume $`Pf^1`$, as appropriate for LFI radiometers (Seiffert et al. Seiffert02 (2002)), Delabrouille assumes a noise spectrum of the form $`Pf^2`$ to account also for possible thermal fluctuations and atmospheric noise in ground based and balloon borne bolometer experiments. This leads to more low-frequency noise for a given knee frequency.

The std of the residual noise $`C_{\mathrm{}}`$ influences the accuracy at which the $`C_{\mathrm{}}`$ spectrum of the CMB can be estimated from the noisy data. From Fig. 14 we can see that at low $`\mathrm{}`$ uniform baselines give the best performance in the sense that the $`C_{\mathrm{}}`$ of the residual noise varies the least from one realization to another. At high $`\mathrm{}`$ there is no clear difference between the performances of different sets of base functions.

## 5 Conclusions

We have presented a maximum-likelihood formulation of the destriping approach to the CMB map-making problem, and a rigorous derivation of the destriping algorithm, and we have applied it to the case of the Planck mission.

We have formulated the method in matrix form, which allows us to apply the conjugate gradient technique in such a way that we can handle very large data sets.

We have compared the three different destriping methods, the one derived here and the other two already presented in the literature, using simulated Planck data (one 100 GHz LFI detector). The differences between these methods can be expressed in terms of a weight function, which varies between methods. This function assigns weights to pixels based on the number of observations falling on that pixel.

We found that our new method provides some improvement to the method used in Burigana et al. (Burigana99 (1997)) and Maino et al. (Maino99 (1999, 2002)). However, our new method was not better than the method given by Delabrouille (Delabrouille98 (1998)), although he gives only a heuristic justification for his weight function. The difference between the latter two methods was insignificantly small, but was systematic. That the maximum-likelihood derivation did not lead to the optimal method in practice is due to actual noise properties differing from the idealization used in the derivation. We recommend using either the weight function derived here ($`w=1/n_p`$) or the one given by Delabrouille ($`w=1/(n_p1)`$).

We have tested the possibility of improving the accuracy of destriping by fitting more base functions besides the uniform baseline, but we have found no systematic improvement in the case of instrumental $`1/f`$ noise. (Fitting several base functions may be more beneficial when removing other types of systematics, i.e. periodic fluctuations induced by thermal instabilities.) The optimal selection of base functions seems to depend on the actual spectrum of the noise, and on which multipoles one is mainly interested in. However, the great virtue of the destriping method is its simplicity: it does not require prior information on the noise spectrum. We lose this advantage if we incorporate information on the noise spectrum into the method.

###### Acknowledgements.

This work was supported by the Academy of Finland Antares Space Research Programme grant no. 51433. TP wishes to thank the Väisälä Foundation for financial support. We thank CSC (Finland) and NERSC (U.S.A.) for computational resources. DM and CB warmly thank all the coauthors of our “classical destriping code” and the colleagues that contributed to the related applications. We acknowledge use of the cmbfast code for the computation of the theoretical CMB angular power spectrum. We gratefully acknowledge K. Górski and B. Wandelt for their implementation of the SDE noise generation method. Some of the results in this paper have been derived using the HEALPix package (Górski et al. gorski99 (1999)). We warmly thank the referee for constructive comments on the first version of this paper.