text:
# Small-Scale structure in the Galactic ISM: Implications for Galaxy Cluster Studies

## 1 Introduction

In the course of analyzing X-ray observations of extragalactic objects, corrections are applied for the effects of Galactic absorption. The amount of Galactic absorption can be fit directly from the observations for adequately strong sources with a known underlying spectrum, and this value can be compared to the amount of HI measured from 21 cm line emission. The similarity or difference between these two measures of the absorption column can be used for a variety of purposes, such as to examine whether there is excess absorption within a cluster of galaxies (White et al., 1991; Arabadjis & Bregman, 2000) or whether absorption by molecular material is commonplace within the Galaxy (Arabadjis & Bregman, 1999a). In other cases, one can fix the X-ray absorption column at the 21 cm HI value and examine whether there is an additional emission component in a system, such as the very soft X-ray excess (the EUVE excess) that has been claimed to exist in clusters of galaxies (e.g., Lieu, Bonamente & Mittaz (1999); but see Bowyer, Berghöfer & Korpela (1999) and Arabadjis & Bregman (1999b)).

In these studies, one usually ignores the small-scale structure of the neutral Galactic ISM, assuming that the mean 21 cm column density within a radio telescope beam is uniform across the face of the beam, which is typically 30 for the single-dish radio telescopes that most of the HI values are based upon. Naturally, there is structure to the total neutral column within the radio telescope beam and this has the potential of influencing the analysis of extragalactic X-ray emission as well as absorption studies in the optical and ultraviolet region. Here we examine the magnitude of the variation of the small-scale structure in the neutral Galactic layer and we discuss the implications for the analysis of data. First, we present an idealized model to illustrate the magnitude of the effect (§2), and then we present the HI and X-ray observations that show the magnitude of the variation in the absorption column (§3, 4). Among other effects, we show that by ignoring such structure in the ISM, absorption corrected X-ray spectra of clusters will systematically show an apparent excess soft emission component, of the sort that is seen in some EUVE observations of galaxy clusters (§5).

## 2 An Example of the Importance of Small-Scale Structure

The importance of small-scale structure on absorption corrections can be illustrated by adopting simple models for the fluctuations in $`N_H`$ that are sensible and are consistent with the observations discussed below. In one model, we assume that the column density distribution $`f(N)`$ has a minimum level, $`N_{low}`$, a maximum level, $`N_{high}`$, and behaves as a power-law between these levels so that

$`f(N)=`$ $`N^m`$ $`\mathrm{for}N_{low}NN_{high}`$
$`f(N)=`$ $`0`$ $`\mathrm{for}N<N_{low},N_{high}<N`$ (1)
The mean column density is the 21 cm HI column, which is

$`N_{avg}=N_{21cm}`$ $`=`$ $`{\displaystyle \frac{_{N_{low}}^{N_{high}}fN𝑑N}{_{N_{low}}^{N_{high}}f𝑑N}}`$
$`=`$ $`{\displaystyle \frac{m1}{m2}}\left({\displaystyle \frac{N_{low}^{(m2)}N_{high}^{(m2)}}{N_{low}^{(m1)}N_{high}^{(m1)}}}\right),`$ (2)
and in the case of $`N_{high}=\mathrm{}`$, we have the simple form $`N_{21cm}=\frac{m1}{m2}N_{low}`$. For individual lines of sight out of the Galaxy, one is more likely to measure a column that is typically closer to the median column density rather than the average column density $`N_{avg}`$. The median column density, or the column density most often be sampled, is

$$N_{median}=2^{1/(m1)}(N_{low}^{(m1)}+N_{high}^{(m1)})^{1/(m1)}$$
(3)
and in the limit of $`N_{high}=\mathrm{}`$, $`N_{median}=2^{1/(m1)}N_{low}`$. Then, the ratio of the median to the mean is, in the limit of $`N_{high}=\mathrm{}`$,

$$N_{\mathrm{median}}/N_{21cm}=\frac{(m2)\mathrm{\hspace{0.17em}2}^{1/(m1)}}{m1}.$$
(4)
If the distribution function were infinitely sharply peaked, $`N_{median}=N_{21cm}`$, but if the distribution has some significant width, these two quantities can be measurably different. For example, with a value of $`m=3.5`$, consistent with the data of Wakker, Oosterloo & Putman (2002), $`N_{median}/N_{21cm}=0.792`$, or about a 20% difference between the median and the mean. This implies that individual lines of sight out of the Galaxy will generally return a value for the column that is about 20% lower than the column induced from the 21 cm measurements. As an example of the importance of this in an astrophysical setting, this phenomenon will affect the interpretation of absorption line measurements (single lines of sight) when they are compared to the amount of 21 cm emission (as was first pointed out to us by Wakker). In a comparison of a Ly$`\alpha `$ absorption line to a 21 cm emission line measurement, one would find the systematic effect that the absorption column was less than the emission column. In the determination of metallicities that compare metal absorption lines to 21 cm emission, one would determine a lower than true metallicity more often than a higher metallicity.

These small scale variations in the HI column will affect the absorption of light in the X-ray waveband in significant ways, especially at high optical depths. The transmitted fraction of light at some energy is given by

$$Tr(E)=I_o(E)\frac{e^{\sigma (E)N(x,y)}𝑑x𝑑y}{𝑑x𝑑y}$$
(5)
At low optical depth ($`\tau (x,y,E)=\sigma (E)N(x,y)<<1`$, where $`\sigma (E)`$ is the energy dependent cross section), this just reduces to $`Tr(E)=I_o(E)(1\sigma (E)N_{avg})`$, so it is appropriate to use $`N_{avg}=N_{21cm}`$, as is usually done. However, if there are variations in $`N(x,y)`$ when $`\tau >1`$, the effective X-ray absorption column changes with energy, with the deviation from $`N_{21cm}`$ becoming greater toward lower energies, which is the direction of increasing optical depth. This can be illustrated in Figure 1, where we show the transmitted fractions as a function of energy for the case where $`N_{avg}=N_{21cm}=2\times 10^{20}`$ cm-2 , and where the variation occurs from $`N_{low}=0.82N_{avg}`$ to $`N_{high}=1.31N_{avg}`$, with $`m=3.5`$. At high optical depths, the transmission is underestimated when using $`N_{21cm}`$. The reason is due to the non-linearity in the behavior of the transmission with column, so that the fluctuations toward low density (lower optical depth) permit far more photons to flow through than an equivalent high density fluctuation removes. This leads to a residual flux when fitting a model with constant $`N_{21cm}`$ to the data, and this positive residual is a systematic effect that is prominent in a moderately narrow energy range where $`4\tau 1`$. The absolute value of the residual is small at low energies because the optical depth through the Galaxy is large, while the residual vanishes at larger energies because the optical depth becomes low enough that using a constant $`N_{21cm}`$ is a good approximation. This effect can be rather important in the discussion of the reality of the soft excess X-ray emission, as we discuss below. However, first it is important to understand the magnitude of the fluctuations in $`N`$ so that the importance of the effects can be assessed.

## 3 The Magnitude of Small-Scale Structure in the ISM from 21 cm Studies

Most extragalactic objects that have been studied lie at Galactic latitudes $`b>20`$, or lines of sight that are out of the plane of the Milky Way. For such lines of sight, there are few studies of the small-scale structure in the total HI column, and almost no information on the structure of H2. The type of information needed is the distribution function of column densities within the larger beam of single-dish surveys (i.e., a histogram of NHI). In principle, this can be obtained with synthesis arrays, but such instruments normally resolve away the larger scale structure, showing only small-scale variations. However, synthesis array observations, used in conjunction with single-dish observations (the zero spacing data) can recover all the emission and show the true distribution of NHI. Unfortunately, few such observations have been made out of the plane, although there are many such observations in the plane, which show a rich amount of structure. Out of the plane, two published observations are by Wakker, Oosterloo & Putman (2002), which is toward a high velocity cloud (WW 187) in the direction of NGC 3783, $`l=287`$, $`b=23`$, and by Joncas, Boulanger & Dewdney (1992), which is toward a filament of far infrared cirrus emission at $`l=141`$$`b=30`$. Although these may not be typical regions, we examine them since they provide a starting point for understanding small-scale structure; we will be obtaining 21 cm data toward more typical sight lines in the near future.

The observation by Joncas, Boulanger & Dewdney (1992) is toward a gas and dust complex in Ursa Major that has conspicuous IR-cirrus emission and is also known to have some molecular gas. However, the mean HI column in this direction is 3.4$`\times 10^{20}`$ cm-2 , which is fairly typical for sight lines at this Galactic latitude. The observation was made with the Synthesis Radio Telescope of Dominion Radio Astrophysical Observatory, including zero spacing data and complementary observations in the continuum; the resulting resolution was 1 and the field of view was 2.6. Since the sensitivity is dropping toward the edges, and the absolute calibration suffers, we have chosen to look at the variations in the central 1 region, smoothed to 2.3 and formed a histogram of the column densities (Figure 2). There is a fairly sharp cutoff in N(HI) for values below 1.7$`\times 10^{20}`$ cm-2 , but there is a distribution that extends to $`5.9\times 10^{20}`$ cm-2 , a range of almost a factor of four (0.54 in the log).

The observation by Wakker is of a HVC projected upon a background AGN (NGC 3783) in the Southern Hemisphere, for the purpose of understanding the optical and ultraviolet absorption line results toward NGC 3783. The synthesis observations were made with the Australia Telescope Compact Array, supplemented with data from the Parkes Multibeam survey, so that HI structures were not resolved out of the image. The map that they produced was not of all of the Galactic HI in that region, but only the HI in the velocity range surrounding cloud, or 190-280 km s-1. Although it would be better to have the entire Galactic column, such a map was not produced, so we have analyzed the existing map, which should give us some insight into HI variations. We have used a square box 1 on a side located in the central region in order to determine the column density histogram (Figure 3), which shows that there is a sharp drop in the column below $`0.65\times 10^{20}`$ cm-2 . The HI distribution is moderately flat to $`1.2\times 10^{20}`$ cm-2 , declining toward higher values ($`16\times 10^{20}`$ cm-2 in a way that may be described as a broken power-law of slopes -3 to -4 (or with an exponential model). A detailed analysis of this distribution is given in Wakker, Oosterloo & Putman (2002).

Although there are few HI maps such as the above, the entire sky has been mapped in the far infrared with DIRBE and IRAS data by Schlegel, Finkbeiner & Davis (1998). These authors have produced maps of E(B-V), which are, on average, proportional to the HI column density, and with 2.67 resolution. We sampled five randomly chosen 1 square regions located at $`l,b`$ = $``$141.3, 37.3; 14.0, 74.8; 112.8, 47.7; $``$44.345.4; and 52.4, 55.9, and analyzed the column density distribution, using the relationship of $`N(HI)=8.18\times 10^{21}E(BV)+2.54\times 10^{19}`$ cm-2 , which was determined by us using other lines of sight out of the Galaxy. It is evident that the full range is at least 100% relative to the minimum value (Figure 4). However, we can better characterize the range by using the quartile points, referenced to the median column density, which in percentage, is $`100(N(75\%)N(25\%))/N_{median}`$, where $`N(75\%)`$ and $`N(25\%)`$ are the columns at the 75% and 25% quartile locations. Using this measure, we find that for these selected regions, the range is fairly narrow, 14.6-19.2% with an average value of 16.5%.

Unfortunately, it is nearly impossible to determine whether the degree of fluctuations seen in the HI is the same as in the IRAS-DIRBE data at a resolution of 1. The HI data from Wakker show a greater degree of HI fluctuations than does the IRAS-DIRBE data, but this may be because the HI map does not include HI at all velocities. In the region studied by Joncas, the HI and IRAS-DIRBE variations are similar in magnitude, but this region was chosen because of the prominent IRAS structure in the field, which guarantees substantial variation. One might expect the IRAS-DIRBE variations to be smaller than the HI variations just because the ”beam” of the IRAS-DIRBE data is an order of magnitude larger than the HI radio synthesis beam of about 1 (so small scale variations are washed out). This is an issue where more data would be extremely valuable.

There is additional evidence for small-scale structure in the interstellar medium, such as from the studies that are based upon optical and ultraviolet absorption measurements. The studies of NaI absorption systems toward globular cluster stars show that there are variations in the equivalent width by a factor of a few on scales of 1-4 (Langer, Prosser & Sneden, 1990) as well as on scales of $`10^{\prime \prime }`$ (Andrews, Meyer & Lauroesch, 2001). Absorption line studies toward binary stars show that there is significant structure in the ISM on a scale of 0.1′′ (Laroesch, Meyer, and Blades 1999). In addition, HI absorption studies of the 21 cm line show that 10-15% of cold HI occurs on scales below 0.1′′. Frail et al. (1994) study the changes in 21 cm absorption toward pulsars with known velocities, detecting changes in the absorption line over a timescale of a few years, which corresponds to a change of 10-100 mas. On comparable angular scales, Faison et al. (1998) used VLBI studies toward radio AGNs and find variations in the HI optical depth as large as a factor of two. Unfortunately, neither the optical absorption lines (NaI, KI) or the 21 cm absorption line is a true linear indicator of the gas column because the optical metal lines can be sensitive to changes in the ionization state while the 21 cm optical depth is inversely proportional to the temperature. Nevertheless, these studies are consistent with the presence of structure in the interstellar column on relatively small scales.

## 4 Small-Scale Variation From X-Ray Observations

The absorption of soft X-rays depends upon the column of intervening gas as well as the energy of the observation, since the approximate relationship for the optical depth is $`\tau NE^3`$, where $`N`$ is the column of X-ray absorbing gas (usually expressed as an equivalent HI column) and $`E`$ is the energy of the photons. For a column density of $`N=2\times 10^{20}`$ cm-2 , $`\tau =1`$ at 0.24 keV, while for a column of $`4\times 10^{20}`$ cm-2 , $`\tau =1`$ occurs at 0.32 keV, so if there are significant fluctuations in the X-ray absorption column, we should see a “drop-out” of the soft emission (the ratio of soft to hard emission decreases). At these column densities, the absorption of X-rays is due primarily to HeI, with HI being the second most important absorber. Other components of the ISM are less important as absorbers (warm ionized gas or hot X-ray emitting gas), so the X-ray column is a fairly good measure of the neutral column.

This effect is, in principle, possible to detect against a relatively smooth extended extragalactic X-ray emitting object such as a cluster of galaxies. In practice, it is possible to search for this effect only in clusters of high surface brightness and for the correct range of Galactic column density, at least when using ROSAT data. The ROSAT PSPC data has rather limited energy resolution so that the data can only be separated into a two or three independent spectral regions. One of these divisions is due to the presence of carbon in the optical path, which leads to a sharp absorption edge at 0.284 keV that becomes progressively less important at higher energies. Operationally, this leads to a local minimum in the detection efficiency near 0.5 keV, so we have divided the PSPC data into two bands, 0.2-0.5 keV and 0.5-2.0 keV (PI channels 20-51 and 52-201). Snowden et al. (1994) have developed procedures to extract diffuse emission in various energy bands, corrected for point sources, instrumental background, and local background sources (solar scattered X-rays, etc.). Here we use those procedures to form the ratio of band R2+R3 (the 0.2-0.5 keV band) to R4+R5+R6+R7 (0.5-2.0 keV).

In order to detect variations in the ratio of these two bands, there must be sufficient photons in the soft band for statistical reasons and that implies a mean Galactic absorption column of less than about $`6\times 10^{20}`$ cm-2 (otherwise there are no soft photons) but the absorption in the soft band should be significant ($`\tau =1`$) so that modest variations in $`\tau `$ lead to detectable variations in the X-ray ratio. For example, if the mean Galactic column is $`3\times 10^{20}`$ cm-2 and there is a region where the column rises to $`4\times 10^{20}`$ cm-2 , our hardness ratio will decrease by 30% while if there is a region where the column rises to $`5\times 10^{20}`$ cm-2 , the hardness ratio will drop to half its original value. From the above discussion, variations of this sort are not likely, although the probabilities of finding such variations can depend strongly on region. For the HI data of Joncas, Boulanger & Dewdney (1992) and Wakker, Oosterloo & Putman (2002) about 20% of the data are above the mean by 4/3 and 10% of the data are above the mean by a factor of 5/3 or greater. However, if one uses the variations in the DIRBE/IRAS data for these two regions and the five other regions (above), typically only 1% of the area has an intensity value that is above the mean by 4/3 or greater. Since the variation in the DIRBE/IRAS data is less than that in the HI for the region observed by Wakker, Oosterloo & Putman (2002), this latter value of 1% may be an underestimate. This suggests that finding significant variation in our hardness ratio is likely if there are several hundred independent positions.

We analyzed five clusters as given in Table 1, all of which are bright nearby systems that are fairly symmetrical in their inner 6. The X-ray temperatures and the bolometric luminosities are taken from Wu, Xue & Fang (1999), for H0 = 75 km s-1 Mpc-1, and the HI columns are those given by Dickey & Lockman (1990).

\[EDITOR: INSERT TABLE 1 HERE.\]

Because galaxy clusters have radial temperature gradients, which would lead to radial softness variation, we examined softness variations azimuthally, either in circular annuli, or in the case of clusters with elliptical isophotes, in elliptical annuli. We avoided the central region and used four annular regions with radii (or semi-major axis) of 1-2 (annulus 1), 2-3 (annulus 2), 3-4.25 (annulus 3) and 4.25-5.5 (annulus 4). We also examined larger annuli, but the surface brightness was generally too low to be of use. In addition, each annular region was divided into 24 regions of 15 each, with region 1 being from position angle 0 to 15, with the position angle increasing east from north as is the usual convention. This leads to 96 independent positions per cluster.

The signature that we are looking for in the X-ray emission is a local hardening of the emission due to the soft photons being absorbed by an excess in the Galactic HI. This will lead to a decrease in the hardness ratio with the change due to decrease in the number of soft photons rather than an increase in the number of hard photons. We discuss each of the systems in turn.

Abell 119: The hardness ratio is shown as a function of annular ring and region (position angle) in Figure 5. There are three positions for which the hardness ratio is more than 3$`\sigma `$ from the mean of the sample (annulus 1, region 2; annulus 2, region 19; annulus 3, region 8). In each case, the deviation is in the sense that the region is harder than the mean and it is the soft band that has decreased relative to it neighbors rather than an increase in the hard band. These are precisely the expectations for regions with excesses in their Galactic HI columns. Also, the next four deviant points in the sample (about 2.5$`\sigma `$; regions 6, 11, and 22 in annulus 3 and region 4 in annulus 4) are also harder than the average and due to a decrease in the number of soft photons.

Abell 2029: There is one region that deviates by more than 3$`\sigma `$ and is in annulus 1 at region 10, but this is due to a positive deviation in the hard flux rather than a decrease in the soft flux, so this does not meet our criteria.

Abell 2142: In this system, we used elliptical isophotes with an eccentricity of 0.82 and an orientation of the major axis of $``$58. There is a strong point source that could not be removed and it affects the first three regions in annuli 2, 3, and 4 (at PA 0-45), so they were removed from the analysis. The photon statistics were good for this long observation, so the annuli were divided into 30 regions, each 12 wide. There are two regions with deviations more than 3$`\sigma `$ (annulus 2, region 26; annulus 3, region 8), both regions being harder than the average and the hardness is due to the lack of soft photons rather than an excess of hard photons (Figure 6).

For the remaining two systems, AWM4 and Abell 3376, there was no evidence of deviations in the hardness ratio beyond 3$`\sigma `$, although the numbers of photons were less than the other clusters, so the statistical significance is poorer.

To summarize, there were about 600 locations in which we could have seen fluctuations in the hardness ratio, although about half those were of lower statistical quality than the others. Six locations deviated from the mean hardness ratio in their annuli by more than 3$`\sigma `$, all toward the region of fewer soft photons. About one such deviation is expected from statistical fluctuations, so this represents an excess. Five of these six regions were harder because of a paucity of soft photons, as would be expected if there were a positive excess in the Galactic neutral hydrogen column. We cannot rule out that these statistically meaningful fluctuations in the hardness ratio is intrinsic to the cluster. To determine whether the fluctuations are due to variations in the Galactic gas will require an accurate HI map in this region. In addition, we will be able to more accurately define the hardness variations with future X-ray observations, especially with XMM-Newton.

Although point sources are largely removed in the data processing, we comment on the effect by point sources that are unresolved and therefore would not be removed. In the ROSAT band, the typical point source is of comparable hardness to the thermal Bremsstrahlung X-ray emission from a galaxy cluster. Also, the surface brightness of a galaxy cluster is much greater than that of the unresolved X-ray background sources. Taken together, we expect that unresolved point sources will not influence our analysis.

## 5 Applications

The fluctuations in the foreground absorption column can have a number of effects, each of which should be modeled. These fluctuations can lead to variable absorption and reddening, which will affect the properties of an HR diagram of a star cluster. Another example occurs in absorption line studies (as discussed above), since the mean neutral gas column inferred from 21 cm emission studies is likely to be higher than the median line of sight, which is more relevant for individual lines of sight. A third example, which we explore in more detail, is the modeling of galaxy clusters, since this will be a systematic effect. This systematic effect either leads one to obtain an X-ray absorption column lower than the mean column (which would help explain the results of Arabadjis & Bregman (1999b), or alternatively, it will lead to a flux residual that might be attributed to a new ”soft” component of the X-ray emission.

The presence of a new soft X-ray component was first presented by Lieu et al. (1996) from data obtained with the EUVE satellite, combined with ROSAT X-ray data. The instrument used by the EUVE for this detection was the DS telescope (Boyer & Malina, 1991), which is sensitive to photons up to about 0.19 keV, which is the soft X-ray region. Soft X-rays can be absorbed by Galactic gas, and for a moderately low column density of $`2\times 10^{20}`$ cm-2 , the optical depth to continuous absorption is unity at 0.24 keV and $`\tau =3`$ at 0.16 keV, so the X-rays that are detected in the EUVE band lie in the heavily absorbed region of the X-ray spectrum. Lieu et al. (1996) took the X-ray spectrum of some bright galaxy clusters and fit an absorption spectrum with a single column density and then compared the flux in the EUVE band to the measured value. They found that the EUVE emission was brighter than anticipated, leading them to conclude that there was a very soft X-ray excess in the emission that needed to be explained.

The validity of this emission has been questioned from two other works on the subject. Arabadjis & Bregman (2000) pointed out that the same soft excess should be visible in the ROSAT PSPC as well, and possibly easier to detect. The effective areas of the PSPC and the EUVE DS telescope are shown in Figure 7, which includes not only the raw collecting area of both instruments, but the effective collecting are due to absorption by Galactic gas with a column density of $`2\times 10^{20}`$ cm-2 , which is fairly typical of the sight lines involved. One finds that the effective EUVE bandpass is 0.144-0.186 keV (FWHM), which can be a great aid because although the PSPC has more effective collecting area in this energy range, its energy resolution is significantly poorer than the EUVE bandpass of 0.04 keV. In this EUVE bandpass, the optical depth is about 3, so if the Galactic column has been overestimated slightly, it makes a large difference in the predicted flux. Galactic columns out of the plane, such as these, are difficult to determine accurately, and Arabadjis & Bregman (1999b) showed that in nearly all cases, the soft X-ray excess would vanish if one adopted a column that was less than 3$`\sigma `$ from the published values. The only exception was the Coma cluster, for which the Galactic column is particularly small. For low HI columns, the corrections to the 21 cm observations are the greatest, so it is possible that the HI column is particularly low and that there is no soft X-ray excess.

The validity of the soft X-ray excess was also questioned by Bowyer, Berghöfer & Korpela (1999), who showed that certain data processing issues could complicate the ability to analyze the data. In particular, they showed that flat-fielding was a key issue, so they constructed flat fields from many long observations, and when this was applied, soft excesses vanished in Abell 2199 and Abell 1795, two clusters that Lieu, Bonamente & Mittaz (1999) claimed to have excesses. However, they find that the Coma cluster still has an excess in its soft X-ray emission, the same result as Arabadjis & Bregman (1999b), working from the ROSAT data.

Here we show that even if one has measured the mean Galactic 21 cm column correctly, fluctuations in that HI column lead to a soft X-ray excess. As an example, when we use the model described in §2, the residual flux has a peak at about 0.19 keV and a substantial amount of flux falls in the EUVE band (Figure 8). In the EUVE band, the fractional excess is claimed to be about 20-50% above the expected values, and that is approximately the amount that we predict from our model of column density fluctuations (Figure 9).

Ideally, one should use the detailed column density map associated with the particular clusters in question, but 21 cm data of this sort does not exist. Consequently, we cannot demonstrate that this effect is responsible for the extraction of a soft X-ray excess; we merely suggest it as a possibility.

## 6 Discussion and Final Comments

We have shown the potential implications of variations in the Galactic absorption column in trying to interpret a variety of observations. The greatest shortcoming in determining the true importance of this effect is the lack of accurate data on the HI variations on the sky at small angular scale. Such observations need a combination of synthesis array observations in combination with single dish data in order to include the ”zero spacing” data. Aside from the two observations discussed above, data of this kind have been largely confined to the Galactic plane, where the Dominion Radio Astronomy Observatory has undertaken a survey from $`74^{}<l<147^{}`$, but in the narrow latitude band $`3.5^{}<b<5.5^{}`$. This is far from the region where most of the galaxy clusters lie, typically $`20^{}<|b|`$.

Another issue that needs to be examined is the degree to which the HI fluctuations in these fields follows the DIRBE/IRAS maps of Schlegel, Finkbeiner & Davis (1998). If it could be shown that the HI generally follows the DIRBE/IRAS data to good accuracy, then one could use that data as a proxy for the HI data. Currently, the situation is unclear since in the WW 187 field, there is a poor correlation between the HI column and the DIRBE/IRAS data, while in the IRAS filament mapped by Joncas, the correlation is much better.

The are a variety of effects that will be clarified by future HI observations, so we strongly encourage observers to obtain such data.

We are especially indebted to Bart Wakker and Gilles Joncas for sending us copies of their data for us to analyze. Also, we would like to thank Thomas Berghöfer for his patience in answering our many questions about EUVE data. We would like to thank Jay Lockman, Jimmy Irwin, Renato Dupke, Morton Roberts, and Stu Bowyer for their advice and comments. Partial support for this work has been provided by NASA through LTSA grant NAG5-10765.


format: text

text:
# Electroweak Corrections to 𝐞⁺𝐞⁻→𝜈𝜈̄𝐇 and 𝐞⁺𝐞⁻→𝐭𝐭̄𝐇This work was supported in part by the Swiss Bundesamt für Bildung und Wissenschaft and by the European Union under contract HPRN-CT-2000-00149.

## 1 Introduction

One of the most important future challenges in particle physics is the understanding of the mechanism of electroweak symmetry breaking and the discovery of the up to now only missing particle of the Electroweak Standard Model (SM), the Higgs boson. The mass of the Higgs boson is expected to be in the range between the current lower experimental bound of $`114.4\mathrm{GeV}`$ and $`1\mathrm{TeV}`$, where a light Higgs boson (with $`M_\mathrm{H}100`$$`200\mathrm{GeV}`$) is favoured by the global fit of the SM to electroweak precision data. While the LHC will find a SM Higgs boson in the full mass range up to $`1\mathrm{TeV}`$ if it exists and has no exotic properties, a complete determination of the Higgs interactions is only possible in the clean environment of an $`\mathrm{e}^+\mathrm{e}^{}`$ collider.

In this work we focus on the reactions $`\mathrm{e}^+\mathrm{e}^{}f\overline{f}\mathrm{H}`$ where $`f`$ is a neutrino or a top quark, which belong to the most interesting Higgs production processes at future $`\mathrm{e}^+\mathrm{e}^{}`$ colliders.

## 2 Radiative corrections to $`𝐞^+𝐞^{}\nu \overline{\nu }𝐇`$

At an $`\mathrm{e}^+\mathrm{e}^{}`$ collider the two main Higgs production mechanisms in the SM are Higgs radiation off Z bosons, so-called Higgs strahlung, and Higgs production via WW fusion. Both mechanisms are present in the reaction $`\mathrm{e}^+\mathrm{e}^{}\nu _l\overline{\nu }_l\mathrm{H}`$ where $`l`$ can be an $`\mathrm{e}`$, $`\mu `$, or $`\tau `$. The l.h.s. of Figure 1 shows the different contributions to the lowest-order cross section as a function of the centre-of-mass (CM) energy $`\sqrt{s}`$ for $`M_\mathrm{H}=150\mathrm{GeV}`$ (see Ref. Denner:2003yg  for details). While the Higgs-strahlung contribution to the total cross section rises sharply at threshold to a maximum of a few tens of GeV above $`\sqrt{s}=M_\mathrm{Z}+M_\mathrm{H}`$ and falls off as $`1/s`$, the WW-fusion channel, which is only present in the reaction $`\mathrm{e}^+\mathrm{e}^{}\nu _\mathrm{e}\overline{\nu }_\mathrm{e}\mathrm{H}`$, dominates the cross section well above the ZH threshold and grows as $`\mathrm{ln}s`$ in the high-energy limit. The difference between the total cross section and the sum ZH+WW is the interference between both channels which is relatively small.

The $`𝒪(\alpha )`$ electroweak corrections to the process $`\mathrm{e}^+\mathrm{e}^{}\mathrm{ZH}`$ have been calculated in Refs. Fleischer:1982af ; Kniehl:1991hk ; Denner:1992bc , and a Monte Carlo algorithm for the calculation of the real photonic corrections to this process was described in Ref. Berends:dw . The electroweak corrections to the full process $`\mathrm{e}^+\mathrm{e}^{}\nu \overline{\nu }\mathrm{H}`$ have attracted a lot of interest recently. Analytical results for the one-loop corrections to this process have been studied in Ref. Jegerlehner:2002es ; however, no numerical results are given there. The contributions of fermion and sfermion loops in the Minimal Supersymmetric Standard Model (MSSM) have been evaluated in Refs. Eberl:2002xd ; Hahn:2002gm . Complete calculation of the full $`𝒪(\alpha )`$ electroweak corrections to $`\mathrm{e}^+\mathrm{e}^{}\nu \overline{\nu }\mathrm{H}`$ in the SM have been performed in Refs. Denner:2003yg ; Belanger:2002ik . These calculations agree within $`0.3\%`$, which is of the same order as the integration error of Ref. Belanger:2002ik .

In the following we briefly summarize some results of our calculation Denner:2003yg  of the $`𝒪(\alpha )`$ electroweak corrections, where we included also corrections from initial-state radiation (ISR) beyond $`𝒪(\alpha )`$ in the structure-function approach. The calculation is done in the $`G_\mu `$ scheme which absorbs the corrections proportional to $`m_\mathrm{t}^2/M_\mathrm{W}^2`$ in the fermion–W-boson couplings and the running of $`\alpha (Q^2)`$ from $`Q^2=0`$ to the electroweak scale. The numerical evaluation of the virtual corrections is particularly complicated due to the appearance of pentagon diagrams. Therefore, we have applied the approach of Ref. Denner:2002ii  which avoids the appearance of inverse Gram determinants in the reduction of tensor 5-point functions. The soft and collinear singularities are treated both in the dipole subtraction method following Refs. Dittmaier:1999mb ; Roth:1999kk  and in the phase-space slicing method following closely Ref. Bohm:1993qx . Two completely independent Monte Carlo programs have been constructed; one applies the multi-channel approach similar to Refs. Roth:1999kk ; Denner:1999gp ; Dittmaier:2002ap , the second uses VEGAS Lepage:1977sw .

The relative corrections are shown on the r.h.s. of Figure 1. The ISR corrections vary strongly in the region of the $`\mathrm{ZH}`$ threshold but are nearly flat for energies above $`400\mathrm{GeV}`$. They are always negative since the lowest-order cross section is continuously rising. The fermionic corrections reach a maximum of about $`6\%`$ in the region where the $`\mathrm{ZH}`$-production channel dominates and are small above $`500\mathrm{GeV}`$ where the WW-fusion process is most important. The non-ISR bosonic corrections exhibit a minimum of about $`7\%`$ at $`300\mathrm{GeV}`$ and are between $`0\%`$ and $`4\%`$ elsewhere. Near the $`\mathrm{ZH}`$ threshold the electroweak corrections are dominated by the $`𝒪(\alpha )`$ ISR corrections from ZH production, while for higher energies, where the main contributions come from WW fusion, the corrections become flat reaching relative corrections of about $`10\%`$ above $`500\mathrm{GeV}`$.

## 3 Radiative corrections to $`𝐞^+𝐞^{}𝐭\overline{𝐭}𝐇`$

Another interesting Higgs production process is the reaction $`\mathrm{e}^+\mathrm{e}^{}\mathrm{t}\overline{\mathrm{t}}\mathrm{H}`$. If the Higgs-boson mass is not too large, i.e. $`M_\mathrm{H}100`$$`200\mathrm{GeV}`$, the Higgs boson is produced mainly through Higgs radiation off top quarks, while emission from intermediate Z bosons plays only a minor role. Therefore, this process can be used for the determination of the top-quark Yukawa coupling, $`g_{\mathrm{t}\overline{\mathrm{t}}\mathrm{H}}`$, which is by far the largest one in the SM ($`g_{\mathrm{t}\overline{\mathrm{t}}\mathrm{H}}0.5`$). If the Higgs boson is light, i.e. $`M_\mathrm{H}120\mathrm{GeV}`$, a precision of around $`5\%`$ can be reached at an $`\mathrm{e}^+\mathrm{e}^{}`$ collider operating at $`\sqrt{s}=800\mathrm{GeV}`$ with a luminosity of $`Ldt1000\mathrm{fb}^1`$ Baer:1999ge . Combining the $`\mathrm{t}\overline{\mathrm{t}}\mathrm{H}`$ channel with information from other Higgs-production and decay processes an even better accuracy can be obtained in a combined fit Battaglia:2000jb .

The $`𝒪(\alpha _s)`$ corrections to the total cross section within the SM have been calculated for the dominant photon-exchange channel in Ref. Dawson:1998ej , while the full set of diagrams has been evaluated in Ref. Dittmaier:1998dz . Recently, considerable progress has been achieved in the calculation of the electroweak corrections to $`\mathrm{e}^+\mathrm{e}^{}\mathrm{t}\overline{\mathrm{t}}\mathrm{H}`$. Results for the electroweak $`𝒪(\alpha )`$ corrections in the SM have been presented in Refs. You:2003zq ; Belanger:2003nm ; Denner:2003ri . While the results of Refs. Belanger:2003nm ; Denner:2003ri  agree well, those of Ref. You:2003zq  differ at large CM energies and close to threshold.

In the last part of this article we show some results of our calculation of the $`𝒪(\alpha )`$ electroweak corrections and of the $`𝒪(\alpha _s)`$ corrections. Although the virtual corrections are much more complex than for $`\mathrm{e}^+\mathrm{e}^{}\nu \overline{\nu }\mathrm{H}`$, we were able to perform the calculation using the same computational techniques as in the former case. Results for total cross sections and various distributions have been presented in Ref. Denner:2003ri .

On the l.h.s. of Figure 2 we show the lowest-order cross section and the cross section including both electroweak and QCD corrections as a function of the CM energy for $`M_\mathrm{H}=150\mathrm{GeV}`$. Away from the kinematic threshold at $`\sqrt{s}=2m_\mathrm{t}+M_\mathrm{H}`$ the total cross section is typically of the order of a few $`\mathrm{fb}`$ and becomes maximal at an energy of about $`800\mathrm{GeV}`$. The relative corrections are presented on the r.h.s. of Figure 2. While the weak bosonic corrections are around $`+10\%`$ close to threshold and fall off rapidly with increasing CM energy, the fermionic corrections are about $`+10\%`$ and depend only weakly on the CM energy. For energies above $`600\mathrm{GeV}`$, the fermionic and weak bosonic contributions partially cancel. The QED corrections, which include both the complete photonic and higher-order ISR corrections, are about $`40\%`$ at threshold and rise to a few per cent at $`1.5\mathrm{TeV}`$. In the electroweak corrections, both QED and weak contributions partially compensate each other. The QCD corrections are positive and rather large in the threshold region, where soft-gluon exchange between in the $`\mathrm{t}\overline{\mathrm{t}}`$ system leads to a Coulomb-like singularity. In the region above threshold, the QCD corrections decrease and even turn negative for energies $`\stackrel{>}{}800\mathrm{GeV}`$.

In summary, for both processes $`\mathrm{e}^+\mathrm{e}^{}\nu \overline{\nu }\mathrm{H}`$ and $`\mathrm{e}^+\mathrm{e}^{}\mathrm{t}\overline{\mathrm{t}}\mathrm{H}`$ we find corrections that are typically of the order of $`\pm 10\%`$ and thus important ingredients in the theoretical predictions for future $`\mathrm{e}^+\mathrm{e}^{}`$ colliders.


text:
# Asymfast, a method for convolving maps with asymmetric main beams

## I Introduction

With the increasing in accuracy and angular scale coverage of the recent Cosmic Microwave Background (CMB) experiments, a major objective is to include beam uncertainties when estimating cosmological parameters bridle , and in particular the asymmetry of the beam arnau . As seen in Cosmosomas cosmosomas , Boomerang boom , Maxima maxima , Archeops archeops , WMAP wmap  and anticipated for Planck, the systematics errors are dominated at the smaller scales (higher multipoles, $`\mathrm{}`$) by the uncertainties on the reconstruction of and deconvolution from the beam pattern. So the beam must be considered more realistically, in particular by rejecting the assumptions of Gaussianity and/or symmetry. It leads to a better modelization of the beam to convolve with and a better estimation of its representation in the harmonic space, its transfer function $`B_{\mathrm{}}^{eff}`$.

The use of simulated sky maps takes an important part in the analysis of astrophysical data sets to study possible systematic effects and noise contributions, and also, to compare the data to theoretical predictions or to observations from other instrumental setups. This is one of the most challenging problems in data analysis : obtaining a simple and accurate model for the instrumental response. In many cases, Monte-Carlo approaches are favoured as they are in general simpler than the analytic ones. However, they need a great number of simulations which often requires too much time of execution for the available computing facilities. A large amount of this time is spent in the convolution of the simulated data by the instrumental response or beam pattern. For an asymmetric beam, the convolved map at a given pointing direction on the sky would depend both on the relative orientation of the beam on the sky and on the shape of the beam pattern. This makes brute-force convolution particularly painful and slow (e.g. buriganaLFI ).

Therefore, either we work in the spherical harmonic space using a general and accurate convolution algorithm (see wandelt  for a fast implementation), either we model the beam pattern by a series of easy-to-deal-with functions and compute for each a fast convolution in the harmonic space.

The work of wandelt  presents how to convolve exactly two band limited but otherwise arbitrary functions on the sphere - which can be the 4-$`\pi `$ beam pattern and the sky descriptions. At each point on the sphere, one computes a ring of different convolution results corresponding to all relative orientations about this direction. To allow subsequent interpolation at arbitrary locations it is sufficient to discretize each Euler angle describing the position on the sphere into $`𝒪\left(L\right)`$ points where L measures the larger of the inverse of the smallest length scale of the sky or beam. The method rapidity depends on the scanning strategy : in $`𝒪\left(L^3\right)`$ for constant-latitude scans and in $`𝒪\left(L^4\right)`$ for other strategy - even if factorization may be found in some cases. This method is used for Planck simulations in case of elliptical Gaussian beams for the polarized channels whansen  and is one of the most promising techniques to decrease side-lobe effects on the WMAP results barnes . While it is efficient for a certain class of observational strategies whansen , it may be difficult to implement in the more general case of non-constant latitude scanning strategies challinor .

For the second method based on the modelling the beam pattern, several solutions have been proposed either circularizing the beam (page  or wu ) or assuming an elliptical Gaussian beam (souradeep  or pablo ). Here we propose an alternative method, Asymfast, which easily can account for any main beam shape. In general the main beam maps can be obtained from point-sources like the planets Jupiter and Saturn (see chiang  for a general view, beamplanck , beamlfi  for HFI, LFI beams and page  for WMAP beams). The beam maps are then fitted with the appropriate model from which the beam transfer function can be computed. Moreover the beam modelization can be used as input for deconvolution methods like burigana  or optimal map-making iterative methods as MapCUMBA dore .

By contrast to wandelt , Asymfast deals only with main beams which are decomposed onto a sum of 2D symmetric Gaussians and does not take into account far-side lobes. In the following, we will simply call Gaussian a 2D symmetric Gaussian function. The sky is then easily convolved by each Gaussian and the resulting convolved sub-maps are combined into a single map which is equivalent to the sky convolved by the asymmetric original beam. This method can be used for any observational strategy and is easy to implement. As the goal of Asymfast is only to deal with main beams, we shall not directly compare with the wandelt  method in this paper.

The use of 4-$`\pi `$ beam is beyond the scope of Asymfast and so Asymfast and the wandelt  method are not directly compared in this paper.

We describe our approach in Section II and the simulations we used to check the accuracy and performance of the method in Section III. Section IV describes the symmetric expansion of the beam. In Section V, we compare the accuracy of our method with respect to the elliptical Gaussian and the brute-force approach. Section VI discusses the time-computing efficiency of the different convolution method considered. Finally, in Section VII, we describe a method to estimate the effective beam transfer function. This technique has been successfully used for the determination of the Archeops main beam in section VIII.

## II Method

Asymfast approaches any asymmetric beam by a linear combination of Gaussian functions centered at different locations within the original beam pattern.

The convolution is performed separetly for each Gaussian and then the convolved maps are combined into a single map. The map convolution with the Gaussian beams is computed in the spherical harmonic space. This allows us to perform the convolution in a particularly low time consuming way. By contrast, the brute-force convolution by an asymmetric beam needs to be performed in the real space for each of the time samples so that the relative orientation of the beam on the sky is properly taken into account for each pointing direction.

The Asymfast method can be described in five main steps :

1. The beam is decomposed into a weighted sum of $`N`$ Gaussians. The number of Gaussians is chosen by minimazing residuals according to user-defined precision (see Section IV).
2. The initial map is oversampled111we consider oversampled sub-maps to reduce the influence of the pixelisation. by a factor of 2 and convolved with each of the Gaussian functions. The sky map is decomposed into $`a_\mathrm{}m`$ coefficients in the harmonic space. Then, $`N`$ sub-maps are computed, each of them smoothed with the corresponding Gaussian sub-beam, by multiplying the $`a_\mathrm{}m`$ coefficients of the original sky map by the current approximation of the transfer functions of the Gaussian beams $`B_{\mathrm{}}^N=\mathrm{exp}\left[\mathrm{}\left(\mathrm{}+1\right)\sigma _N^2\right]`$.
3. The sub-maps are deprojected into timelines using the scanning strategy of the corresponding sub-beam position.
4. The $`N`$ timelines are stacked into a single one weighted by their sub-beam amplitude.
5. The timeline obtained this way is projected onto the sky with the scanning strategy corresponding to the center of the beam pattern and we obtain the convolution of the original sky map with the fitted beam along the scan.

The HEALPix healpix  package is used to store maps (ring description), to compute the decomposition of the sky map in spherical harmonics (anafast) and to reconstruct maps from the $`a_\mathrm{}m`$ coefficients (synfast).

## III Beams and scanning strategy simulations

The method was checked using realistic simulations of a sky observation performed by an instrument with asymmetric beam pattern and complex scanning strategy. To keep the time consumption reasonnable within our computing capabilities without degrading the quality of the results, we have chosen beams with FWHM between 40 and 60 arcmin sampled by steps of 4 arcmin (square beam map of 60$`\times `$60 pixels corresponding to 4$`\times `$4 deg).

We consider two sets of simulations : a first one corresponding to a quasi-circular beam to which we have added 1 % (simulation 1a) and 10 % (simulation 1b) of noise, and a second one corresponding to an irregular beam to which we have also added noise in the same way (simulations 2a and 2b). These two sets are obtained from a sum of a random number of elliptical Gaussians with random positions within 0.4 degree around the center and random FWHMs and amplitudes. Then, the beams are smoothed with a 4 arcmin width Gaussian and white noise is added.

The scanning strategy is assumed to be a set of consecutive meridians with sampling of 3 arcmin over each meridian and a lag of 3 arcmin between two meridians. This strategy is roughly similar to the Planck observing plan. It allows an efficient check of the method as, close to the equator, all beams are parallel and the effect of the orientation is maximum. However, close to the pole, the high number of beams in different direction per pixel makes the effective beam more circular. In order to save time, we keep only half an hemisphere of the sky, so we use about 6.5 million data samples. The number of hits per pixel is highly variable from equator to pole, ranging from 1 to about 8000 near the pole.

For high-resolution instruments, like Planck, we should consider maps with resolution of the order of 1 arcmin. However, this would require too much computing time for brute-force convolution. Instead, in this paper, we consider maps of the sky with pixel size of $``$24 arcmin which are stored in the HEALPix format ($`N_{side}=256`$). We have checked that the results obtained with Asymfast considering lower resolution maps with larger beams can be generalized to higher resolution maps with narrower beams.

## IV Symmetric Gaussian expansion of the beam

We model the original beam pattern using $`N`$ Gaussians, $`i\{1,N\}`$, of the form

$$g(x_i,y_i,\sigma _i;x,y)=\frac{1}{2\pi \sigma _i^2}\mathrm{exp}\left[\frac{(xx_i)^2+(yy_i)^2}{2\sigma _i^2}\right]$$
(1)
The beam is fitted with a weighted sum of the $`N`$ Gaussians :

$$b(x,y)=\underset{i=1}{\overset{N}{}}A_ig(x_i,y_i,\sigma _i;x,y)$$
(2)
so we have $`4N`$ free parameters corresponding to width ($`\sigma _i`$, FWHM = $`\sqrt{8ln2}\sigma _i`$), amplitude ($`A_i`$) and center position of each Gaussian ($`x_i,y_i`$).

The optimal value of $`N`$ depends on the required precision. We perform the fit using 1 to 10 Gaussians (which is typically enough to attain residuals of the order of the noise level) and compute the quadratic deviation to the original beam pattern

$$S(N)=\frac{1}{n_p}\underset{p=1}{\overset{n_p}{}}\left[b_N^{fit}(p)b(p)\right]^2.$$
(3)
where $`n_p`$ represents the total number of pixels, $`b_N^{fit}`$ corresponds to the fitted beam with $`N`$ Gaussians and $`b`$ to the original beam pattern.

We consider now the two sets of simulations (simulations 1a, 1b, 2a and 2b) described in Section III. The fit is performed using a least square fit in the non-linear case following the algorithm described in brandt . Figure 1 presents the main results for these two sets of simulations : the first two rows for simulations 1a and 1b, and the last two rows for simulations 2a and 2b. In each case, we have plotted, from left to right :

* the initial simulated beam pattern
* the initial simulated beam pattern + noise
* the quadratic deviation $`S`$ as a function of the number of Gaussians $`N`$ (the dotted line represents the number of Gaussians chosen for the final beam pattern model)
* the best-fit model for the beam pattern
* the residual map after substraction of the input noise in percents
* the histogram of the residual map including input noise, fitted to a Gaussian (in red) with residuals to the Gaussian (in blue)

The distribution of $`S(N)`$ is modelled by a decreasing exponential, $`exp(\tau N)`$, plus a constant $`k`$. The number of Gaussians $`N`$ is chosen such that $`N`$ is the smaller value verifying

$$S(N)k<\frac{t}{\tau }$$
(4)
where $`t`$ is a threshold defined by the user for the required precision.

Figure 1 shows that a small number of Gaussians is enough for a good fit (typically less than 10 Gaussians give less than 2 % residuals). In some cases, the algorithm does not converge and therefore no data point is plotted on the figure. Lack of convergence may appear for some combination of number of Gaussians. This is generally solved by using a smaller or larger number of Gaussians in the fit.

The efficiency of the fit is illustrated by the distribution of the residuals which are very close to the white noise distribution centered at zero with the dispersion corresponding to the level of the input noise. Note also that the reconstructed beams are marginally sensitive to the noise level and that the apparent decomposition in three different peaks (case 1b) leads to residuals with amplitude only of the level of the noise.

Once the Gausssians parameters have been estimated, map simulations using any pointing strategy can be done quickly and precisely.

## V Accuracy of the Asymfast method

In this section, we compare the accuracy of Asymfast to that of other common approaches including the modelling of the beam pattern either by a single circular Gaussian or by an elliptical Gaussian. We use as a reference the brute-force convolution with the true beam pattern at each pointing direction.

For this purpose, we have convolved, using the different methods described before, a map containing 26 point sources with the same amplitudes uniformly distributed over half an hemisphere. The level of accuracy for each of the former methods of convolution is estimated by computing the quadratic deviation of the convolved map obtained for that particular method with respect to the convolved map obtained by brute-force convolution.

Table 1 represents the quadratic deviation for each of the methods as a percentage of the maximum of the original point-sources map. This quantity can be interpreted as the percentage of spurious noise introduced by the convolution method. In the case of quasi-circular beams (simulations 1a and 1b), Asymfast is 1.5 (resp. 6) more accurate than the elliptical (resp. circular) approximation. For more realistic irregular beams (simulations 2a and 2b), Asymfast is about 4 (resp. 10) times more accurate than the elliptical (resp. circular) approach. Moreover Asymfast depends neither on the shape of the beam pattern nor on the noise level.

Figure 2 shows, for the two sets of simulations considered (simulations 1a, 1b on the top panel and simulations 2a, 2b and the bottom panel), the histogram of the percentage of the quadratic deviation residual maps for the Asymfast (solid line), elliptical (dashed line) and circular (dotted line) approximations. The figure confirms the results shown in Table 1, indicating that Asymfast is a very good approximation to both quasi-circular and irregular beams (residuals smaller than 5%). By contrast, the elliptical approximation can only be safely used in the case of quasi-circular beams (residuals about 20% for simulations 2). The circular approximation, as expected, is very poor in any case (residuals of the order 40%).

In conclusion, the Asymfast results are very close to the standard solution while the circular and elliptical Gaussian approximations are very poor for realistic beam patterns. In addition, the standard approach, which works in real-space, necessarily uses a beam pattern with fixed resolution. By contrast, Asymfast works in the spherical harmonic space where the resolution is only limited by the pixelisation of the sky map. Moreover, Asymfast can obtain full resolution in frequency with no mixing. Obtaining the same frequency resolution with brute-force convolution requires a full-sky beam pattern which increases considerably the computing time.

## VI Time-computing efficiency of Asymfast vs brute-force convolution

As discussed in the following section, the determination in the spherical harmonic space of the effective circular transfer function of an asymmetric beam pattern, $`B_{\mathrm{}}^{eff}`$, requires a large number (of the order of 1000) of accurate Monte-Carlo simulations of fake data appropriately convolved by the instrument beam. Therefore, any algorithm used for this purpose need to be much faster than the brute-force convolution procedure to keep the computing time reasonable.

To quantify the time-computing efficiency of Asymfast, we consider $`n_{simu}=1000`$ Monte-Carlo simulations of a convolved full-sky map for a single Planck-like detector. We use a squared map of $`n_{beam}=3600`$ pixels to describe the beam pattern and simulated timelines for 12 months of Planck observations which correspond to $`n_{scan}=\mathrm{5\; 10}^9`$ samples. With Asymfast, we model the beam pattern using $`N=10`$ Gaussians. For each simulation, the total number of pixels222We have considered the HEALPix convention healpix  on the full-sky map is $`n_{pix}=12N_{side}^2`$ for three values of $`N_{side}`$ : 256, 512 and 1024 considered.

For these simulations, Table 2 shows the computing-time, in terms of arbitrary CPU-units, for circular beam convolution in the harmonic space using the HEALPix package, and for asymmetric beam convolution using Asymfast or brute-force convolution algorithms.

We observe that Asymfast is more than 50 times faster than the brute-force convolution for the high resolutions ($`N_{side}=512`$ ou $`1024`$). This is because the convolution in the spherical harmonic space used by Asymfast is much faster than in real space. In addition, Asymfast is particularly efficient for Monte-Carlo purposes because the beam modelization and the computing of the pointing directions corresponding to each of the $`N`$ sub-beam are performed only once.

Furthermore, for high resolutions, the Asymfast computing-time is just a factor of 2 larger than the computing-time needed for the convolution of $`N`$ circular Gaussians in spherical harmonic space. The extra computing time with respect to the convolution of $`N`$ circular Gaussians comes from the projection and the deprojection operations performed by Asymfast.

Finally, the difference in computing-time between Asymfast and the brute-force convolution also increases with the beam map resolution as the second method depends linearly on beam map number of pixels while Asymfast is very marginally sensitive to it, only at the beam fit step.

## VII Application to CMB analysis : estimation of the beam transfer function $`B_{\mathrm{}}`$

In this section, we describe how Asymfast can be used to estimate this effective circular transfer function for realistic asymmetric beams. A beam is described in the spherical harmonic space by a set of coefficients $`b_\mathrm{}m`$ for each scanning orientation at each pointing position. The Gaussianity of the primordial fluctuations of the CMB is a key assumption of modern cosmology, motivated by simple models of inflation. So the angular power spectrum of the CMB should be a circular quantity and it is common to just consider an effective circular beam transfer function, $`B_{\mathrm{}}^{eff}`$. The “circularization”of the beam may be obtained by several ways : assuming a Gaussian beam (the easiest way, leading to a simple analytic description), computing the optimal circularly symmetric equivalent beam (wu , applied to Maxima) or fitting the radial beam profile by a sum of Hermite polynomials (page , applied to WMAP). The analytic expression of the beam transfer function of elliptical Gaussian beams has also been studied in details by pablo  and by souradeep  who applied it to Python V.

With Asymfast, the “circularization”of the beam depends on the scanning strategy as it is computed on the map after convolution by the asymmetric main beam. So it is comparable to or more precise than other methods, depending on how the scans intersect each other.

For a circular Gaussian beam (angles $`\theta `$ and $`\varphi `$ corresponding to the beam direction),

$$b(\theta ,\varphi )=\frac{1}{2\pi \sigma ^2}\mathrm{exp}\left[\frac{\theta ^2}{2\sigma ^2}\right],$$
(5)
the beam transfer function in case of a pixel scale much smaller than the beam scale is approximated by white

$`B_{\mathrm{}}`$ $`=`$ $`\mathrm{exp}\left[{\displaystyle \frac{1}{2}}\mathrm{}(\mathrm{}+1)\sigma ^2\right].`$ (6)
Note that in this case, the beam transfer function does not depend on the orientation of the beam on the sky and is given by a simple analytical expression.

For a moderate elliptical Gaussian beam pattern, we can compute an analytic approach pablo  to the beam transfer function by introducing a small perturbation to Eq. 5 so that

$$b^{ell}(\theta ,\varphi )=\frac{1}{2\pi \sigma ^2}\mathrm{exp}\left[\frac{\theta ^2}{2\sigma ^2}f(\varphi )\right]$$
(7)
where $`f(\varphi )`$ describes the deviation from circularity.

For general irregular asymmetric beams, the previous approximation is not well adapted and the orientation of the beam on the sky has to be taken into account. For this purpose, we estimate an effective beam transfer function, $`B_{\mathrm{}}^{eff}`$, from Monte-Carlo simulations which use the Asymfast method for convolution as follows :

1. We convolve CMB simulated maps by the beam pattern using Asymfast.
2. For each simulation, we estimate a $`B_{\mathrm{}}^{eff}`$ inverting the equation master

$$\stackrel{~}{C}_{\mathrm{}}=M_{\mathrm{}\mathrm{}^{}}B_{\mathrm{}^{}}^{eff}{}_{}{}^{2}C_{\mathrm{}^{}}$$

where $`M_{\mathrm{}\mathrm{}^{}}`$ is the coupling kernel matrix that takes into account the non-uniform coverage of the sky map, $`\stackrel{~}{C}_{\mathrm{}}`$ is the pseudo-power spectrum computed on the convolved map and $`C_{\mathrm{}}`$ is the input theorical model.
3. We compute the effective transfer function of the beam, by averaging the $`B_{\mathrm{}}^{eff}`$ obtained for each of the simulations.

We have tested this method on simulations of pure circular Gaussian beams for which we have obtained a beam transfer function fully compatible with Eq. 6 as expected.

Figure 3 shows the estimate of the effective beam transfer function $`B_{\mathrm{}}^{eff}`$ for the two sets of beam simulations discussed in the previous sections (simulations 1a, 1b, 2a, 2b) using a Monte-Carlo of 25 simulations. The results of this Monte-Carlo (solid line) is compared to the analytic transfer function of a circular approximation of the beam (semi-dashed line) and to the elliptical approximation computed from pablo  (dashed line). The top panel represents the beam transfer functions whereas the bottom panel shows the differences to the Asymfast estimation for both analytical approximations (circular, semi-dashed line and elliptical, dashed line). The Asymfast effective beam transfer function shows some irregularities at very low $`\mathrm{}`$ due to the cosmic variance.

As expected, for the quasi-circular beam (simulations 1a and 1b), the estimate transfer functions for the three methods are similar within 10 %. The differences observed between our approach and the elliptical approximation are mainly due to the fact that the latter does not take into account the scanning strategy. Due to the complex but realistic scanning strategy used in our simulations, we expect that, in general, a given position on the sky will be observed with different relative beam orientations and therefore the effective beam will appear more circular.

By contrast, for a more irregular beam pattern (simulations 2a and 2b), the differences between the three methods are much larger (from 10 % at low $`\mathrm{}`$ up to 70 % at high $`\mathrm{}`$). The beam transfer function obtained using the elliptical approximation follows better the one obtained using Asymfast. We expect that the complex scanning strategy would make the effective beam more circular. However as the beam pattern is very asymmetric, the effective beam will contain complex highly irregular structures which cannot be mimicked by an oriented elliptical beam. Therefore, the largest differences are found at high resolution.

Note that the results presented here are also valid for higher resolution beams. For this paper, we have considered low resolution beams to be able to directly compare the Asymfast convolution to standard brute-force convolution.

## VIII Application to Archeops

Asymfast has been applied to the main beam of the Archeops bolometers which are identical to the Planck-HFI ones. The beam shapes were measured on Jupiter beamarcheops  and are for most of them moderately elliptical, the multimode ones being rather irregular. Figure 4 shows the study of the Archeops 217 GHz photometer used for the first Archeops CMB analysis archeops . On the first row we show the $`1\times 1`$ deg with 1 arcminute pixel maps of the main beam, at left the initial map, at right the reconstructed one with 10 Gaussians. The second row presents the residual map and its histogram. The latter is fitted with a Gaussian (in red), the difference to it is shown in blue. The Monte-Carlo estimation using Asymfast described in Section VII (solid line) is compared to the analytic transfer function of a circular approximation of the beam (semi-dashed line), to the elliptical approximation to order 4 from pablo  (dashed line) and to the estimation obtained by simulation used in archeops  (dotted line). They are shown on the third row figure. The bottom panel shows the differences to the Asymfast beam transfer function estimation for the other three approaches.

The beam is well reproduced as the residual map do not exhibit any structure. Moreover the dispersion of the distribution of the residuals is compatible to the noise level of the initial map. The deviation to the Gaussian and to the elliptical models arises from $`\mathrm{}200`$, i.e. the top of the first acoustic peak and culminates around $`\mathrm{}500600`$, i.e. at the second peak location. Therefore the fine structure of main beams with equivalent FWHM of 13 arcminutes is already determinant for the measurement of the second acoustic peak with Archeops. Because of the irregularity of the beam, we are more sensitive to higher resolution structures.

## IX Conclusions

Asymfast is a fast and accurate convolution procedure particularly well-adapted to asymmetric beam patterns and complex scanning strategies which are often used in CMB observations. Asymfast can both produce convolved maps from input timelines and compute, from Monte-Carlo simulations, an accurate circular approximation to the transfer function, $`B_{\mathrm{}}^{eff}`$, of any asymmetric beam pattern. The computing time needed to obtain a convolved map is dominated by the HEALPix software computing time. So it scales as $`𝒪\left(n_{pix}^{3/2}\right)`$ where $`n_{pix}`$ is the number of pixels of the map, with a multiplying factor depending on the number of Gaussians.

Asymfast models any general beam pattern by a linear combination of circular 2D Gaussians, permitting an accurate reconstruction of the instrumental beam, with residuals smaller than 1 % (compared to 4 % for an elliptical Gaussian model). In addition, Asymfast convolution is at least a factor of 50 faster than the brute-force convolution algorithm for full-sky maps of 12.5 million pixels and even faster at higher resolution. This allows us to perform a large number of Monte-Carlo simulations in a reasonable computing time to estimate accurately the effective circular transfer function of the beam pattern.

By constrast to other modelling techniques like page  and pablo , Asymfast can be used with non-circular and non-elliptical beam patterns. Furthermore Asymfast approximates the main beam pattern while wandelt  uses an exact 4-$`\pi `$ beam description, nevertheless it can be equally easily applied to any general scanning strategy while the feasibility of wandelt  strongly depends on the former challinor .

Note that Asymfast is a general convolution algorithm which can be also used successfully in many other astrophysical areas to reproduce the effects of asymmetric beam patterns on sky-maps and to compare observations from independant instruments which requires the cross-convolution of the datasets. Any circular functions with analytic description in the harmonic space may be used instead of 2D symmetric Gaussians.

###### Acknowledgements.

The authors would like to thank F.-X. Désert for fruitfull discussions. The Healpix package was used throughout the data analysis healpix . We also want to thank the anonymous referee for pointing out some issues which helped us to improve the manuscript.