text:
# Beam deconvolution in noisy CMB maps

## 1 Introduction

Many experiments have been designed to measure Cosmic Microwave Background (CMB) anisotropies at small angular scales. Recent and new generation of experiments make use of multi-frequency and multi-beam instruments at a focus of a meter class telescope. Since not all the feeds can be placed along the optical axis of the telescope, the majority of them are necessary off-axis and the corresponding beams show more or less relevant optical aberrations (see e.g. Page et al. 2003 for a recent discussion of the main beam shape in the context of the WMAP project), according to the experiment optical design. For example, in the contex of the ESA Planck project 111http://astro.estec.esa.nl/Planck/ (see e.g. Bersanelli et al. 1996, Tauber 2000), significant efforts have been carried out to significantly reduce the main beam distortions produced by optical aberrations (see e.g. Villa et al. 1998, Mandolesi et al. 2000a). On the other hand, even optimizing the optical design, a certain level of beam asymmetry can not be completely eliminated (see e.g. Sandri et al. 2002, 2003).

If not taken into account in the data analysis, the main beam distortion introduces a systematic effect in the data (Burigana et al. 1998, Mandolesi et al. 2000b) that affects the reconstructed map quality and, in particular, the recovery of the angular power spectrum of the CMB anisotropy (Burigana et al. 2001, Arnau et al. 2002, Fosalba et al. 2002). The main topic of this paper is beam deconvolution in this type of experiments with the aim of remove the main beam distortion effect in the recovery of the angular power spectrum of CMB temperature anisotropy. We wish to reverse the weighted average performed by a non-circular rotating beam in the presence of a significant level of uncorrelated instrumental noise. The true angular power spectrum ($`C_{\mathrm{}}`$ quantities before beam smoothing) should be recovered from the deconvolved maps, at least, in a large enough interval ($`\mathrm{}_{min}`$, $`\mathrm{}_{max}`$).

A preliminary work about beam deconvolution was presented in Arnau & Sáez (2000). In that paper, two methods for beam deconvolution were considered. One of them (hereafter method I) is based on the Fourier transform, and the other one (method II) uses the Jacobi algorithm for solving algebraic systems of linear equations. Applications of these methods in very simple cases were presented. The first method was applied in the case of elliptical non-rotating beams in the presence of a very low level of Gaussian instrumental noise. The other method was only used for a spherical beam in the total absence of noise. More realistic situations must be considered. This is the main goal of this paper.

The formalism of our approach to deconvolution is presented in Sect. 2.

Beam deconvolution can be only studied using simulations. Map making requires a pixelisation, and the accuracy of the angular power spectrum obtained from pixelised maps strongly depends on experiment sensitivity and resolution and on the sky coverage. In the case of small angular scales ($`\mathrm{}\genfrac{}{}{0pt}{}{>}{}\mathrm{\hspace{0.17em}100}`$), the angular power spectrum can be well estimated using small squared maps with edges lesser than $`20^{}`$ (Sáez et al. 1996). In this case, a good map making algorithm and an appropriate power spectrum estimator were described in (Arnau et al. 2002). In a first step (first part of Sect. 3), we work with this type of simulations by assuming a simple elliptical main beam shape. An observational strategy involving repeated measures in the same pixel but without a detailed reference to a specific experiment is adopted at this stage.

Afterwards, we apply the method to more complicated simulations carried out by using the HEALPix 222http://www.eso.org/science/healpix/ (Hierarchical Equal Area and IsoLatitude Pixelization of the Sphere) package by Gòrski et al. (1999) to pixelise the maps and compute the angular power spectrum from them. The beam size, its asymmetry, the variance of the instrumental noise, and the beam rotation associated to an observational strategy simulating the Planck observations at 100 GHz are considered in the second part of Sect. 3.

Finally, we discuss our results and draw the main conclusions in Sect. 4.

We work in the framework of an inflationary flat universe (adiabatic fluctuations) with dark energy (cosmological constant) and dark matter (cold). In this $`\mathrm{\Lambda }`$CDM model, the density parameters corresponding to baryons, dark matter, and vacuum, are $`\mathrm{\Omega }_b=0.05`$, $`\mathrm{\Omega }_d=0.25`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, respectively, and the reduced Hubble constant is $`h=0.65`$. No reionization is considered at all. All the simulations are based on the CMB angular power spectrum corresponding to this model, which has been computed with the CMBFAST code by Seljak & Zaldarriaga (1996).

## 2 Beam smoothing and deconvolution

Let us begin with an asymmetric non-rotating beam which smoothes a map $`T`$ to give another map $`T_s`$. In the continuous formalism, we can write:

$$T_s(\theta ,\varphi )=B(\theta \theta ^{},\varphi \varphi ^{})T(\theta ^{},\varphi ^{})𝑑\mathrm{\Omega }^{},$$
(1)
where the beam is described by function $`B`$. If the beam centre points towards a point with spherical coordinates ($`\theta `$, $`\varphi `$), the weight associated to another direction ($`\theta ^{}`$, $`\varphi ^{}`$) is a function of the form $`B(\theta \theta ^{},\varphi \varphi ^{})`$.

In the absence of rotation, function $`B`$ only depends on the differences $`\theta \theta ^{}`$ and $`\varphi \varphi ^{}`$ and, consequently, Eq. (1) is a mathematical convolution. In such a case, the Fourier transform can be used (as it was explained by Arnau and Sáez 2000) to perform beam deconvolution. Nevertheless, if the asymmetrical beam rotates (as a result of the observational strategy), the function describing the beam is of the form $`B=B(\theta \theta ^{},\varphi \varphi ^{},\theta ,\varphi )`$. Since the beam is different (distinct orientations) when its centre points towards different points in the sky, a new dependence on the angles $`\theta `$ and $`\varphi `$ has appeared. With a $`B`$ function involving this dependence, Eq. (1) is not a mathematical convolution and the method I, which is based on the Fourier transform, does not work.

In practice, non-circular beams rotate due to the observational strategy and, moreover, the effect of this rotation is not negligible in many cases \[see Arnau et al. 2002 for an estimation and Burigana et al. 2001 for an application to Planck LFI (Low Frequency Instrument, Mandolesi et al. 1998)\]. In this situation, method I cannot be used; however, as we are going to show along the paper, method II works.

Let us assume a certain pixelisation and an asymmetric beam which smoothes maps of the CMB sky. We first consider that only one observation per pixel is performed. The beam could have either the same orientation for all the pixels or different orientations for distinct pixels; in both cases, the beam smoothes the sky temperature $`T`$ to give $`T_s`$ according to the relation:

$$T_s^i=\underset{i=1}{\overset{M}{}}B_{ij}T^j,$$
(2)
where $`M`$ is the total number of pixels in the map. Quantity $`B_{ij}`$ is the beam weight corresponding to pixel $`j`$ when the beam centre points inside pixel $`i`$. Eq. (2) can be seen as a system of $`M`$ linear algebraic equations, in which, the independent terms $`T_s^i`$ are the observed temperatures, and the M unknowns are the true sky temperatures $`T^j`$. The Jacobi method can be tried in order to solve this system. The solution would be the deconvolved map with temperatures $`T^j`$. No problems with pixel dependent beam orientations (rotation) are expected a priori; nevertheless, rotating asymmetric beams lead to a $`B_{ij}`$ matrix which is more complicated than that corresponding to a circular beam (studied in Arnau & Sáez 2000); by this reason, we are going to verify that Jacobi method works for any beam, first in the absence of noise and, then, when there is an admissible noise level. In matrix form, Eq. (2) can be written as follows: $`T_s=BT`$.

If we now consider that each pixel is observed N times either with an unique beam and different orientations per pixel or with various non-circular rotating beams (as it occurs in projects as Planck where there are various beams for each frequencies), then, we can write N matrix equations (one for each measure) of the form $`T_s^{(\alpha )}=B^{(\alpha )}T`$, where index $`\alpha `$ ranges from $`1`$ to $`N`$. The average temperature assigned to pixel $`i`$ is $`T_a^i=(1/N)_{\alpha =1}^NT_s^{(\alpha )i}`$ and the above system of N matrix equations leads to

$$T_a=B_aT,$$
(3)
where matrix $`B_a`$ describes the average beam, which can be calculated as follows:

$$B_a^{ij}=\frac{1}{N}\underset{\alpha =1}{\overset{N}{}}B^{(\alpha )ij};$$
(4)
hence, for a given experiment involving various measurements per pixel, the average beam (4) might be calculated at each pixel and, then, we can try to use the Jacobi method to solve the system of linear equations (3).

We see that, in the absence of noise, method II could work in the most general case, in which various beams move according to the most appropriate observational strategy. For a multi-beam experiment one should in principle simulate the effective scanning strategy and the convolution with the sky signal for each beam and then apply the formalism described above by taking into account the various resolutions and shapes of the beams. On the other hand, the power of this method is that it works independently of the small differences between the resolutions and shapes of the various beams at the same frequency in a given experiment. Therefore, considering the data from a single average beam, instead of the data from the whole set of beams, but with the sensitivity per pixel obtained by considering the whole set of receivers at the given frequency in the case in which the noise is taken into account, allows to reduce the amount of data storage and simplify the analysis without introducing a significant loss of information about the accuracy of the method.

Instrumental uncorrelated Gaussian noise makes beam deconvolution more problematic; nevertheless, as we demonstrate in this work (Sects. 3.2 and 3.4), method II works for experiments with a level of noise similar to that of Planck, or better, since the effect of deconvolution on the noise can be quite accurately understood with Monte Carlo simulations.

## 3 Applying method II

In this section, the deconvolution method II is applied in various cases using appropriate simulations. First, in cases A and B, method II is applied to deconvolve a set of small sky patches with regular pixelisation. To make this part of the work almost independent of the detailed scanning strategy of the considered CMB anisotropy experiment, we adopted an observational strategy involving multiple observations of a given pixels and only roughly mimicking that of Planck. The beam shape is assumed to be elliptical. Case A does not involve any noise. Case B is identical to case A except for the presence of noise. Method II is then applied to deconvolve larger sky patches but using the HEALPix package for the sky pixelization and the computation of the angular power spectrum from coadded and deconvolved maps, simulating the Planck observational strategy, and assuming one of the beam shapes simulated in the past year for LFI, both in the absence of noise (case C) and in the presence of noise (case D).

### 3.1 Case A: noiseless, small patches

An elliptical beam of the form

$$B(\theta \theta ^{},\varphi \varphi ^{})=W__Ne^{\left[\frac{(\theta \theta ^{})^2}{2\sigma _\theta ^2}\frac{(\varphi \varphi ^{})^2}{2\sigma _\varphi ^2}\right]}$$
(5)
is assumed, as in Burigana et al. (1998) and Arnau et al. (2002).

It is also assumed that quantities $`\sigma _\theta `$ and $`\sigma _\varphi `$ obey the relations $`\sigma _\theta \sigma _\varphi =\sigma ^2`$ and $`\sigma _\theta /\sigma _\varphi =1.3`$, where $`\sigma =4.54^{}`$ ($`\theta _{_{FWHM}}10.68^{}`$). With this choice, the elliptical beam (5) mimics the 100 GHz Planck beams for some locations of the detectors over the focal plane.

We simulate squared $`14.6^{}\times 14.6^{}`$ patches, with 256 (128) nodes per edge; thus, our pixel size is $`\mathrm{\Delta }=3.43^{}`$ ($`\mathrm{\Delta }=6.86^{}`$). These sizes are allowed by HEALPIX and, consequently, this choice will facilitate some comparisons. We use seventy five of these regions covering about the forty per cent of the sky. With this coverage and $`\mathrm{\Delta }=3.43^{}`$ ($`\mathrm{\Delta }=6.86^{}`$), the angular power spectrum can be estimated (from simulated maps) with good accuracy from $`\mathrm{}_{min}=100`$ to $`\mathrm{}_{max}=10800/\mathrm{\Delta }3100`$ ($`\mathrm{}_{max}1550`$). See Sáez & Arnau (1997) for details about partial coverage. In the theoretical model under consideration (see Sect. 1), the CMB temperature is a Gaussian homogeneous and isotropic statistical two dimensional field. In such a case, a certain method proposed by Bond & Efstathiou (1987) can be used to make the $`14.6^{}\times 14.6^{}`$ maps used in this paper. This method is based on the following formula:

$$\frac{\delta T}{T}=\underset{s_1,s_2=N}{\overset{N}{}}D(\mathrm{}_1,\mathrm{}_2)e^{i(\theta \mathrm{}_1+\varphi \mathrm{}_2)},$$
(6)
where $`\mathrm{}_1=2\pi s_1/\mathrm{\Lambda }`$, $`\mathrm{}_2=2\pi s_2/\mathrm{\Lambda }`$, and $`\mathrm{\Lambda }`$ stands for the angular size of the square to be mapped. This equation defines a Fourier transform from the position space ($`\theta `$, $`\varphi `$) to the momentum space ($`\mathrm{}_1`$, $`\mathrm{}_2`$). The Gaussian quantities $`D(\mathrm{}_1,\mathrm{}_2)`$ have zero mean, and their variance is proportional to $`C_{\mathrm{}}`$, where $`\mathrm{}=(\mathrm{}_1^2+\mathrm{}_2^2)^{1/2}`$. Since $`\delta T/T`$ is real, the relation $`D(\mathrm{}_1,\mathrm{}_2)=D^{}(\mathrm{}_1,\mathrm{}_2)`$ must be satisfied. From given $`C_{\mathrm{}}`$ coefficients, the above $`D(\mathrm{}_1,\mathrm{}_2)`$ quantities can be easily calculated and, then, according to Eq. (6), a Fourier transform leads to the map. Sáez et al. (1996) used this map making algorithm to get very good simulations of $`20^{}\times 20^{}`$ squared regions.

In the case of small squared maps, the above map making method suggests the power spectrum estimator used in Sects. 3.1 and 3.2 and also in Arnau et al. (2002). Given one of these maps $`\delta T/T(\theta ,\varphi )`$, an inverse Fourier transform leads to quantities $`D(\mathrm{}_1,\mathrm{}_2)`$ and, then, the average $`|D(\mathrm{}_1,\mathrm{}_2)|^2`$ can be calculated on the circumference $`\mathrm{}^2=\mathrm{}_1^2+\mathrm{}_2^2`$. Some interpolations are necessary to get the $`D(\mathrm{}_1,\mathrm{}_2)`$ values at the points located on the circumference. The resulting average is proportional to $`C_{\mathrm{}}`$, where $`\mathrm{}`$ is the radius of the circumference.

For $`\mathrm{\Delta }=3.43^{}`$ ($`\mathrm{\Delta }=6.86^{}`$), the beam average can be restricted to a square with seventeen (nine) nodes per edge. Outside this square, beam weights given by Eq. (5) appear to be negligible in this context 333The signal entering at angles larger than $`1^{}`$ from the beam centre direction produces the so-called straylight contamination, dominated by the Galactic emission (see e.g. Burigana et al. 2001, 2003), a systematic effect different from the main beam distortion effect considered in this work..

Our elliptical beam is rotating while it covers a given $`14.6^{}\times 14.6^{}`$ patch. In order to simulate beam rotation (see Fig. 1), the squared patch is located with random orientation (angle $`\alpha `$) in the plane ($`\theta `$ , $`\varphi `$ ), and a different beam orientation is assigned to each pixel of the patch. If Q is the centre of a certain pixel, we find a point P on the $`\theta `$-axis which is the centre of an auxiliary circumference with radius $`85^{}`$ which passes by Q and, then, the beam orientation – in the pixel under consideration – is fixed by assuming that the major axis of the elliptical beam is tangent – at point $`Q`$ – to the auxiliary circumference. The distance from the patch centre, C, to the $`\theta `$-axis is random, but it is constrained to be smaller than $`74^{}`$ in order to ensure the existence of an auxiliary circunference passing by the centre of every pixel.

For each patch, the sky (T field) is simulated using either 256 or 128 nodes per edge and, afterwards, the beam described above is used to get the smoothed map ($`T_s`$). The pixel temperatures of the $`T_s`$ map are the independent terms of Eq. (2) and, moreover, the terms of the $`B`$ matrix can be built up when necessary using Eq. (5) and beam orientation. Taking into account that all the terms of the B diagonal are identical to the central weight of the beam $`b`$, the $`n+1`$ iteration of the Jacobi method can be written as follows:

$$T^{(n+1)}=T^{(n)}+b^1T_sb^1BT^{(n)},$$
(7)
where $`T^{(n)}`$ is the previous one. At zero order, we take $`T^{(0)}=T_s`$.

Since the map is a $`256\times 256`$ ($`128\times 128`$) square and the beam is another $`17\times 17`$ ($`9\times 9`$) square (see Fig. 2), when the beam centre points towards a pixel located outside the ninth (fifth) row or column (counting from the nearest boundary), there are no map temperatures to be weighted. In practice, for CMB maps, we have verified that the following procedure works very well: Write an equation for every internal node where the beam average is well defined and, then, solve the resulting system, which has $`240\times 240`$ ($`120\times 120`$) equations and the same unknowns. The remaining temperatures (external points) are used when required by beam smoothing, but they are not altered along the iterative process.

Fig. 3 shows the main results obtained in case A. Top (bottom) left panel shows quantities $`\mathrm{}(\mathrm{}+1)C_{\mathrm{}}/2\pi `$ in units of $`\mu K^2`$ before smoothing (continuous line) and after deconvolution (pointed line) for $`\mathrm{\Delta }=3.43`$ ($`\mathrm{\Delta }=6.86`$). Both curves are indistinguishable except at the largest $`\mathrm{}`$ values included in the Figure. The relative deviations between the dotted and dashed lines of each panel are given in the corresponding right panels. The relative error introduced by deconvolution – in the absence of noise – is smaller than $`5\%`$ (0.5 %) for $`\mathrm{}1900`$ ($`\mathrm{}1480`$) in the case $`\mathrm{\Delta }=3.43`$ ($`\mathrm{\Delta }=6.86`$); the deviations grow beyond the sixth (fifth) acoustic peak.

### 3.2 Case B: noisy, small patches

In this case, beam, rotation, coverage and pixelizations are identical to those of case A; however, there is instrumental uncorrelated Gaussian noise with $`\sigma __N=9\mu K`$ ($`\sigma __N=4.5\mu K`$), in antenna temperature, for $`\mathrm{\Delta }=3.43`$ ($`\mathrm{\Delta }=6.86`$), just the noise expected by combining all the beams of Planck at 100 GHz. A joint treatment of the impact of main beam distortions and of correlated $`1/f^\alpha `$ type noise (see e.g. Seiffert et al. 2002) and other kind of instrumental systematics (see e.g. Mennella et al. 2002) is out of the scope of this paper. On the other hand, this does not represents a crucial limitation, since blind destriping algorithms can strongly reduce the impact of these effects (see e.g. Delabrouille 1998, Maino et al. 1999, Mennella et al. 2002) also in the presence of optical distortions (Burigana et al. 2001) and, possibly, of non negligible foreground fluctuations (Maino et al. 2002). The system to be solved has the form:

$$T_s^i=\underset{i=1}{\overset{M}{}}B_{ij}T^j+N^i,$$
(8)
where $`N^i`$ is the noise at pixel $`i`$. Using matrices, this equation can be written as follows:

$$T_s=B(T+B^1N).$$
(9)
This last equation is formally identical to the matrix form of Eq. (2) and it can be solved in the same way – using Jacobi method – to find the map $`T+B^1N`$. After applying this method, some numerical error $`E`$ is expected and, consequently, the numerical solution of system (9) is of the form

$$T^{}=T+B^1N+E.$$
(10)
In general, $`T^{}`$ is different from $`T`$ (sky temperature before smoothing); hence, the angular power spectrum extracted from the map $`T^{}`$ is different from that of the unconvolved sky, which is to be extracted from map $`T`$. Results are shown in Fig. 4, which has the same structure as Fig. 3. We see that the spectra before smoothing and after deconvolution (which are obtained from maps $`T`$ and $`T^{}`$, respectively) separate at middle $`\mathrm{}`$ values. In the right panels, we can verify that the deviation produced by deconvolution – in the presence of the assumed level of noise – is smaller than five per cent for $`\mathrm{}1100`$ ($`\mathrm{}950`$) in the case $`\mathrm{\Delta }=3.43`$ ($`\mathrm{\Delta }=6.86`$). The relative deviations rapidly increases for greater $`\mathrm{}`$ values. From the comparison of Figs. 3 and 4 it follows that the presence of noise with $`\sigma __N=9\mu K`$ ($`\sigma __N=4.5\mu K`$) has important consequences for beam deconvolution.

Fortunately, the angular power spectrum of the map $`N^{}=B^1N`$ can be estimated and subtracted from that of the $`T^{}`$ map. In order to do that, we first solve the matrix equation:

$$N^{}=BN^{}.$$
(11)
This equation can be also solved using the Jacobi method. The independent terms are the temperatures $`N^{}`$ corresponding to a new noise realization different from $`N`$. We can now extract the $`C_{\mathrm{}}`$ quantities from the map $`N^{}`$. New numerical errors should appear when we use the Jacobi method in the presence of noise. On the other hand, we can try a Monte Carlo approach to evaluate the deconvolution effect on the noise. We can take various $`N^{}`$ noise realizations to get an average spectrum of the corresponding $`N^{}`$ maps; in this way, the effect of the noise variance in the estimate of the angular power spectrum of $`N^{}`$ is strongly reduced (forty noise realizations suffice). When we subtract this spectrum from that of $`T^{}`$, namely, when we correct the $`T^{}`$ spectrum taking into account noise effects, results are much better than those showed in Fig. 4 444We observe that an analogous approach can be pursued also in the presence of correlated noise, provided that the noise properties can be known from laboratory measures and/or directly reconstructed from the data (Natoli et al. 2002). Of course, in this context, destriping (or, possibly, map-making, see e.g. Natoli et al. 2001) should be previously applied both to the data and to the simulated pure noise data.. These new results are presented in Fig. 5. The structure of this figure is identical to that of Figs. 3 and 4. As in Fig. 3, the range of $`\mathrm{}`$ values – in the left panels – has been appropriately chosen to include the region where the displayed curves separate significantly. The deconvolution is better than that of Fig. 4, where no correction for the noise has been considered. According to the right panels of Fig. 5, the relative deviation produced by deconvolution plus correction is smaller than five per cent for $`\mathrm{}1500`$ ($`\mathrm{}1300`$) in the case $`\mathrm{\Delta }=3.43`$ ($`\mathrm{\Delta }=6.86`$). For $`\mathrm{\Delta }=3.43`$ ($`\mathrm{\Delta }=6.86`$), deconvolution works very well up to the end of the fifth (fourth) acoustic peak. For equivalent levels of noise in different pixels, we see that – as expected – deconvolution has recovered more $`C_{\mathrm{}}`$ coefficients for $`\mathrm{\Delta }=3.43`$; hence, we can say that results corresponding to $`\mathrm{\Delta }=3.43`$ are sensibly better than those of $`\mathrm{\Delta }=6.86`$. On account of this fact and also for the sake of briefness, pixels with $`\mathrm{\Delta }=6.86`$ are not considered in cases C and D below.

### 3.3 Case C: application to Planck in the absence of noise

The selected orbit for Planck is a Lissajous orbit around the Lagrangian point L2 of the Sun-Earth system (see e.g. Bersanelli et al. 1996). The spacecraft spins at 1 r.p.m. and the field of view of the two instruments – LFI and HFI (High Frequency Instrument, Puget et al 1998) – is about $`10^{}\times 10^{}`$ centered at the telescope optical axis (the so-called telescope line of sight, LOS) at a given angle $`\alpha `$ from the spin-axis direction, given by a unit vector, $`s`$, chosen to be pointed in the opposite direction with respect to the Sun. In this work we consider values of $`\alpha 85^{}`$, as adopted for the baseline scanning strategy. The spin axis will be kept parallel to the Sun–spacecraft direction and repointed by $`2.5^{}`$ every $`1`$ hour (baseline scanning strategy). Hence Planck will trace large circles in the sky and we assume here, for simplicity, 60 exact repetitions of the set of the pointing directions of each scan circle. A precession of the spin-axis with a period, $`P`$, of $`6`$ months at a given angle $`\beta 10^{}`$ about an axis, $`f`$, parallel to the Sun–spacecraft direction (and outward the Sun) and shifted of $`2.5^{}`$ every $`1`$ hour, may be included in the scanning strategy, possibly with a modulation of the speed of the precession in order to optimize data transmission (Bernard et al. 2002). The quality of our deconvolution code is of course almost independent of the details of these proposed scanning strategies, and we assume here the baseline scanning strategy for sake of simplicity.

The code implemented for simulating Planck observations for a wide set of scanning strategies is described in detail in Burigana et al. (1997, 1998) and in Maino et al. (1999). In this study we do not include the effects introduced by the Planck orbit, to be currently optimized, by simply assuming Planck located in L2, because they are fully negligible in this context.

We compute the convolutions between the antenna pattern response and the sky signal as described in Burigana et al. (2001) by working at $`3.43^{}`$ resolution and by considering spin-axis shifts of $`2.5^{}`$ every hour and 7200 samplings per scan circle. We simulate 11000 hours of observations (about 15 months) necessary to complete two sky surveys with the all Planck beams.

With respect to the reference frames described in Burigana et al. (2001), following the recent developments in optimizing the polarization properties of LFI main beams (see e.g. Sandri et al. 2003), the conversion between the standard Cartesian telescope frame $`x,y,z`$ and the beam frame $`x_{bf},y_{bf},z_{bf}`$ actually requires a further angle $`\psi _B`$ other than the standard polar coordinates $`\theta _B`$ and $`\varphi _B`$ defining the colatitude and the longitude of the main beam centre direction in the telescope frame. Appendix A provides the transformation rules between the telescope frame and the beam frame, as well as the definition of the reference frames adopted in this work.

The orientation of these frames as the satellite moves is implemented in the code. For each integration time, we determine the orientations in the sky of the telescope frame and of the beam frame, thus performing a direct convolution with the sky signal by exploiting the detailed main beam response in each considered sky direction. The detailed main beam shape and position on the telescope field of view adopted in this application is that computed in the past year for the feed LFI9 (Sandri et al. 2002) which shows an effective FWHM resolution of $`10.68^{}`$ and deviations from the symmetry producing a typical ellipticity ratio of 1.25. Such values of resolution and asymmetry parameter are in the range of them that it is possible to reach with a 1.5 m telescope like that of Planck by optimizing the optical design (see e.g. Sandri et al. 2003). Although our deconvolution method is largely independent of the details of the considered beam shape, it is interesting to exploit its reliability under quite realistic conditions.

The CMB anisotropy map has been projected into the HEALPix scheme (Gòrski et al. 1999) starting from the angular power spectrum of the assumed $`\mathrm{\Lambda }`$CDM model (see Sect. 1).

To make the application of the deconvolution code easier and the system solution possible without large RAM requirement and in a reasonable computational time 555In the current implementation, about 19 hours of computation are required to deconvolve a single patch with $`1024^2`$ pixels on an 64 bit alpha digital unix machine with single cpu at 533 MHz and 1 Gb RAM. we implemented a code that identify in the simulated time ordered data (TOD) all and only the beam centre pointing directions in an equatorial patch (in ecliptic coordinates) of $`1024\times 1024`$ pixels with a $`3.43^{}`$ side ($`n_{side}=1024`$). We keep the exact information on the beam centre pointing direction and the beam orientation (defined for instance by an angle between the axis $`x_{bf}`$ and the parallel in the beam centre pointing direction) as computed by our flight simulator. All the samples of the TOD within the same pixel are identified and restored in contiguous positions. At this aim, we take advantage from the nested, hierarchical ordering of the HEALPix. This is quite simple in the current simplified simulation, but it will require the development of efficient and versatile tools to manage the more general case in which the all samples from the experiment multi-beam array are considered, particularly for the ecliptic polar patches, which pixels are observed many and many times because of the Planck scanning strategy. In the context of the Planck project, this effort will be pursued by taking advantage from the development of Planck Data Model (see e.g. Lama et al. 2003).

¿From the simulated TOD, possibly restored as described above, we extract a map of a patch of simply coadded data and a map of a patch deconvolved by applying method II. The latter map can be then symmetrically smoothed with a beam FWHM of $`10.68^{}`$ by using the HEALPix tools for comparison with the former one, obtained from the convolution with the simulated asymmetric beam and taking into account the scanning strategy. Of course, from the input map we can extract the same sky patch.

We consider four different patches covering an equatorial region of $`28.3`$ % of the sky (analogously to the case of small patches, see Sect. 3.1, avoiding the boundary regions of the four patches slightly reduces the originally considered, $`33.3`$ %, sky coverage).

All the above maps are inverted with the anafast code of HEALPix to extract the correponding angular power spectra. The result is shown in Fig. 6. Of course, all the angular power spectra are in strict agreement at multipoles $`\genfrac{}{}{0pt}{}{<}{}\mathrm{\hspace{0.17em}200}`$, where the main beam distortion effect is negligible for a beam with a FWHM of about $`10^{}`$ and a reasonable ellipticity. Note the very good agreement between the power spectrum of the input map and that derived with method II: the difference becomes significant only at the seventh acoustic peak \[compare the solid (black) line with the dashed (green) line\]. Note the power excess at high $`\mathrm{}`$ introduced by the beam distortions, compared with the power spectrum derived from the deconvolved map subsequently symmetrically smoothed \[compare the dash-three dots (fuxia) with the dash-dots (blue)\]. Also, the power spectrum derived from coadded map when divided by the window function corresponding to the symmetric equivalent beam, $`\mathrm{exp}[(\sigma \mathrm{})^2]`$, significantly exceeds that of the input map at $`\mathrm{}`$ larger than the fourth acoustic peak \[compare the dots (red) with the solid line (black)\]. We find a similar disagreement even by varying the assumed value of the symmetric beam width $`\sigma `$: an improvement on a limited range of multipoles results in a worsening on a different range of multipoles.

This demonstrates that a kind of deconvolution is necessary to remove the main beam distortion effect at very high multipoles and that method II works well in the absence of noise. The impact of instrumental noise is discussed in the next subsection.

### 3.4 Case D: application to Planck in the presence of noise

Analogously to the Case B (see Sect. 3.2) we have simulated four patches (of $`1024^2`$ pixels, as in the previous section) of instrumental uncorrelated Gaussian noise with $`\sigma _N=9\mu `$K (in terms of antenna temperature) for a pixel of $`3.43^{}`$, as appropriate to the global Planck sensitivity at 100 GHz; for simplicity we have assumed a uniform noise.

This realization of noise map has been added to the map of signal and then method II has been applied to deconvolve the map of signal plus noise, as described in Sect. 3.2. We produced also four realizations of pure noise maps to be deconvolved in the same way. Finally, we generated another map of noise to be superimposed to the coadded map obtained from the convolution with the simulated beam including scanning strategy, for comparison.

We computed the angular power spectrum of the four maps of pure noise and of the four maps of pure noise deconvolved with method II. Fig. 7 (left panel) compares the averages of the four realizations of these power spectra and their relative variance (right panel). As evident, deconvolution increases the noise: a rough approximation of the ratio between of the noise angular power spectrum after deconvolution and before deconvolution is given by $`2(\mathrm{FWHM}/\mathrm{\Delta })^2\mathrm{exp}[(\sigma \mathrm{}/2)^2+(\sigma \mathrm{}/2)^6]`$ where, as usual, $`\mathrm{FWHM}=\sqrt{8\mathrm{l}\mathrm{n}2}\sigma `$ ($`=10.68^{}`$) and $`\mathrm{\Delta }`$ is the pixel side ($`=3.43^{}`$). On the other hand, the relative variance of these power spectra is almost similar.

In spite of the relatively large increase of the noise power, we find that method II results to work quite well in removing the effect of main beam distortions up to the end of the fifth acoustic peak, when the average power spectrum of the deconvolved pure noise maps is subtracted to the power spectrum of the deconvolved noisy map (see the dashed (green) line in Fig. 8).

In Fig. 9 we report the relative (per cent) errors introduced by beam distortions in the absence of deconvolution, in the presence of deconvolution without applying the subtraction of the average deconvolved noise spectrum and by applying the deconvolution and the subtraction of average deconvolved noise spectrum. As evident, in the last case the power spectrum can be recovered with a good accuracy up to high multipoles (relative errors $`\genfrac{}{}{0pt}{}{<}{}\mathrm{\hspace{0.17em}5},10,15,20`$ % respectively for $`\mathrm{}\genfrac{}{}{0pt}{}{<}{}\mathrm{\hspace{0.17em}1250},1470,1500,1650`$ – see the dashed (green) line in the middle panel – to be compared with errors $`10,20,30`$ % at $`\mathrm{}1200,1250,1300`$ – see the dotted (red) line in the middle panel – and then dramatically increasing with $`\mathrm{}`$ in the absence of deconvolution). The right panel shows what already found for the noiseless case (Sect. 3.3): even by varying the assumed value of the symmetric beam width $`\sigma `$, an improvement in the $`C_{\mathrm{}}`$ recovery can not be reached simultaneously on the whole relevant range of multipoles.

We have considered here for simplicity only equatorial patches. On the other hand, we have verified that the main beam distortion affects the reconstructed power spectrum in a similar way also on polar patches. In fact, even at high ecliptic latitudes a given pixel is preferentially observed with a limited number of beam orientations. In polar patches the Planck sensitivity is significantly better than the average. Therefore, in spite of the more complex data storage, deconvolution will be there less affected by the noise 666For example, the average sensitivity on a polar region of about 25 squared degrees is about 5 times better (i.e. $`2\mu `$K on pixel of $`3.43^{}`$!) than the average full sky sensitivity: then, we expect there a deconvolution quality intermediate between that found here and that found in the previous section..

## 4 Discussion and conclusions

We have presented the basic formalism for a robust and feasible method for deconvolution in noisy CMB maps.

We have implemented and tested this method for two completely different situations: small patches observed with an elliptical beam and with a scanning strategy involving repeated measures of the same pixel but not related to a specific experiment and quite large sky areas observed with a realistic beam and a scanning strategy like that of the Planck satellite. A sensitivity level and a beam resolution typical of the Planck experiment at 100 GHz have been exploited.

After having verified the good accuracy of our deconvolution code to remove the main beam distortion effect in the angular power spectrum estimation in the absence of noise, we applied it to noisy maps. We demonstrate that it is possible to accurately evaluate the effect of deconvolution on pure noise simulated maps so deriving, with Monte Carlo simulations, a good estimation of the average deconvolved noise angular power spectrum to be subtracted from the deconvolved noisy maps.

In practice, for the considered sensitivity, $`9\mu `$K for a pixel of $`3.43^{}`$, and beam resolution, FWHM $`10^{}11^{}`$, our deconvolution code allows to efficiently remove the main beam distortion effect and accurately reconstruct the CMB power spectrum up to the end of the fifth acoustic peak.

Clearly, in the context of the Planck project, the measure of the very high multipole region of the CMB angular power spectrum will be take advantage from the cosmological frequency channels at highest resolution, namely the 217 GHz channel (having a FWHM $`5^{}`$), where Poisson fluctuations from extragalactic sources (see e.g. Toffolatti et al. 1998) are expected to be at a very low level and anisotropies from thermal Sunyaev-Zeldovich effects are, if not exactly null because of possible unbalanced contributions within the bandwidth, certainly very small. On the other hand, the frequency range about 100 GHz is where the global (Galactic plus extragalactic) foreground contamination is expected to be minimum. Therefore, it is extremely relevant to extract at these frequency an accurate estimation of the sky angular power spectrum, cleaned, as better as possible, from the all systematic effects. In addition, the removal of the main beam distortion effect, relevant at large multipoles, greatly helps the comparison between the results obtained at different frequency channels.

Of course, we plan to apply this method in the future also to lowest and highest beam resolutions and in the presence of other kinds of systematic effects. We believe that the results found here are very encouraging, suggesting that the main beam distortion effect, previously reduced by optimizing the optical design, can be further reduced in the data analysis.

###### Acknowledgements.

This work has been partially supported by the Spanish MCyT (project AYA2000-2045). A part of the calculations was carried out on a SGI Origin 2000s at the Centro de Informática de la Universidad de Valencia. It is a pleasure to thank the LFI DPC, were some computations were carried out, and the staff working there for many useful discussions. We warmly thank M. Sandri and F. Villa for having provided us with one of the LFI main beam, simulated at high resolution, adopted in some parts of this work. Some of the results in this paper have been derived using the HEALPix (Gòrski et al. 1999). We thank U. Seljak and M. Zaldarriaga for the use of the CMBFAST code.

## Appendix A Transformation rules between Planck telescope frame and beam frame

Let $`s`$ be the unit vector, choosen outward the Sun direction, of the spin axis direction and $`\widehat{k}`$ that of the direction, $`z`$, of the telescope line of sight (LOS), pointing at an angle $`\alpha 85^{}`$ from the direction of $`s`$.

On the plane tangent to the celestial sphere in the direction of the LOS we choose two coordinates $`x`$ and $`y`$, respectively defined by the unit vector $`\widehat{i}`$ and $`\widehat{j}`$ according to the convention that the unit vector $`\widehat{i}`$ points always toward $`s`$ and that $`x,y,z`$ is a standard Cartesian frame, referred here as telescope frame.

Let $`\widehat{i}_{bf},\widehat{j}_{bf},\widehat{k}_{bf}`$ be the unit vectors corresponding to the Cartesian axes $`x_{bf},y_{bf},z_{bf}`$ of the beam frame; $`\widehat{k}_{bf}`$ defines the direction of the beam centre axis in the telescope frame. The beam frame is defined with respect to the telescope frame by three angles: $`\theta _B`$$`\varphi _B`$$`\psi _B`$ ($`\theta _B`$ and $`\varphi _B`$, two standard polar coordinates defining the direction of the beam centre axis, range respectively from $`0^{}`$, for an on-axis beam, to some degrees, for LFI off-axis beams, and from $`0^{}`$ to $`360^{}`$).

Let $`\widehat{i^{}}_{bf},\widehat{j^{}}_{bf},\widehat{k^{}}_{bf^{}}`$ ($`\widehat{k^{}}_{bf}=\widehat{k}_{bf}`$) be the unit vectors corresponding to the Cartesian axes $`x^{},y^{},z^{}`$ of an intermediate frame, defined by the two angles $`\theta _B`$ and $`\varphi _B`$, obtained by the telescope frame $`x,y,z`$ when the unit vector of the axis $`z`$ is rotated by an angle $`\theta _B`$ on the plane defined by the unit vector of the axis $`z`$ and the unit vector $`\widehat{k}_{bf}`$ up to reach $`\widehat{k}_{bf}`$:

$$\widehat{k^{}}_{bf}=\widehat{k}_{bf}=\mathrm{cos}(\varphi _B)\mathrm{sin}(\theta _B)\widehat{i}+\mathrm{sin}(\varphi _B)\mathrm{sin}(\theta _B)\widehat{j}+\mathrm{cos}(\theta _B)\widehat{k}$$
(12)
$$\widehat{i^{}}_{bf}=[\mathrm{cos}(\varphi _B)^2\mathrm{cos}(\theta _B)+\mathrm{sin}(\varphi _B)^2]\widehat{i}+[\mathrm{sin}(\varphi _B)\mathrm{cos}(\varphi _B)(\mathrm{cos}(\theta _B)1)]\widehat{j}\mathrm{sin}(\theta _B)\mathrm{cos}(\varphi _B)\widehat{k}$$
(13)
$$\widehat{j^{}}_{bf}=[\mathrm{sin}(\varphi _B)\mathrm{cos}(\varphi _B)(\mathrm{cos}(\theta _B)1)]\widehat{i}+[\mathrm{cos}(\theta _B)\mathrm{sin}(\varphi _B)^2+\mathrm{cos}(\varphi _B)^2]\widehat{j}\mathrm{sin}(\theta _B)\mathrm{sin}(\varphi _B)\widehat{k}.$$
(14)
The beam frame is obtained from the intermediate frame through a further (anti-clockwise) rotation of an angle $`\psi _B`$ (ranging from $`0^{}`$ to $`360^{}`$ 777We note that, in other conventions, angles $`\varphi _B^{}`$ and $`\psi _B^{}`$ ranging from $`180^{}`$ to $`180^{}`$ are given, instead of $`\varphi _B`$ and $`\psi _B`$. The angles $`\varphi _B`$ and $`\psi _B`$ here defined are equal to $`\varphi _B^{}`$ and $`\psi _B^{}`$ when they are positive and are given respectively by $`360^{}+\varphi _B^{}`$ and $`360^{}+\psi _B^{}`$ for negative $`\varphi _B^{}`$ and $`\psi _B^{}`$.) around $`\widehat{k}_{bf}`$ and is therefore explicitely given by:

$`\widehat{i}_{bf}`$ $`=`$ $`[\mathrm{cos}(\psi _B)\widehat{i^{}}_{bf,x}+\mathrm{sin}(\psi _B)\widehat{j^{}}_{bf,x}]\widehat{i}+[\mathrm{cos}(\psi _B)\widehat{i^{}}_{bf,y}+\mathrm{sin}(\psi _B)\widehat{j^{}}_{bf,y}]\widehat{j}`$ (15)
$`+[\mathrm{cos}(\psi _B)\widehat{i^{}}_{bf,z}+\mathrm{sin}(\psi _B)\widehat{j^{}}_{bf,z}]\widehat{z}`$
$`\widehat{j}_{bf}`$ $`=`$ $`[\mathrm{sin}(\psi _B)\widehat{i^{}}_{bf,x}+\mathrm{cos}(\psi _B)\widehat{j^{}}_{bf,x}]\widehat{i}+[\mathrm{sin}(\psi _B)\widehat{i^{}}_{bf,y}+\mathrm{cos}(\psi _B)\widehat{j^{}}_{bf,y}]\widehat{j}`$ (16)
$`+[\mathrm{sin}(\psi _B)\widehat{i^{}}_{bf,z}+\mathrm{cos}(\psi _B)\widehat{j^{}}_{bf,z}]\widehat{z},`$
where the bottom index $`x`$ $`(y,z)`$ indicates the component of intermediate frame unit vector along the axis $`x`$ $`(y,z)`$ of the telescope frame, as defined by eqs. (A1–A3).


format: text

text:
# An Unusual Discontinuity in the X-ray Surface Brightness Profile of NGC 507: Evidence of an Abundance Gradient?

## 1 Introduction

Chandra observations of early galaxies, groups of galaxies, and clusters of galaxies have demonstrated that the hot, X-ray emitting coronal gas of these systems is often not resting idly in their gravitating dark matter halos and that many of these systems are far from relaxed. Several different structures have been identified and studied in detail. Chandra observations of shocks and ”cold fronts” associated with sub-cluster mergers have contributed to our knowledge of the formation of structure. Evidence for large scale, coherent, non-hydrostatic motions of the central regions of clusters has also been reported (Vikhlinin et al., 2001; Markevitch et al., 2001; Mazzotta et al., 2001). One of the most interesting, and to some degree unexpected, results from Chandra has been the study of diverse interactions between jets and radio lobes of radio galaxies and the ambient ISM/IGM/ICM (Fabian et al., 2000; Finoguenov and Jones, 2000; McNamara et al., 2000; Jones et al., 2002; Kraft et al., 2003). The study of these various motions and hydrodynamical interactions is of fundamental importance to several astrophysical and cosmological problems. These observations give us insight into the merger process and subsequent relaxation. By studying merger shocks, thermal histories, etc., we learn about the process of the formation of clusters. In addition, understanding the radio lobe/ICM interactions may be critical in understanding the life cycle of clusters as there may be a cyclical relationship between cooling flows, nuclear activity, and nuclear outflow.

NGC 507 is the dominant, massive elliptical galaxy of a nearby ($`z=0.01646`$) group/poor cluster (the so-called Pisces cluster). Previous ROSAT observations have shown that NGC 507 is one of the most X-ray luminous early type galaxies ($`L_x`$ 1043 ergs s-1) in the local universe with evidence of a large cooling flow (cooling rate of 30-40 $`M_{}`$ yr-1) at the center (Kim and Fabbiano, 1995; Paolillo et al., 2003). This group has an unusually large mass-to-light ratio ($``$125), similar to that observed in the X-ray overluminous elliptical galaxies (OLEGs) (Vikhlinin et al., 1999). NGC 507 is also known to be an FR I radio galaxy (Fanti et al., 1986).

In this letter, we present results from a Chandra/ACIS-I observation of the hot ISM of this group. We detect a sharp discontinuity in the surface brightness profile 55 kpc E and SE of the center of the galaxy spanning a 125 arc. This discontinuity has been seen in XMM-Newton observations of NGC 507 as well (Fabbiano et al., 2002). A more detailed discussion of the larger scale emission and the central cooling flow peak based on ROSAT/HRI and Chandra/ACIS-S observations is presented elsewhere (Paolillo et al., 2003). Initially it appears that this discontinuity was another example of a ”cold-front” as seen in clusters of galaxies. Such discontinuities have recently been detected in clusters of galaxies and are commonly attributed to either infalling sub-clumps or non-hydrostatic motions of the central core (e.g. A2142 (Markevitch et al., 2000), A3667 (Vikhlinin et al., 2001), A1795 (Markevitch et al., 2001)). However, there are significant problems with this interpretation of the discontinuity observed in NGC 507. In particular, the temperature jump across the discontinuity is small (0.1-0.2 keV), so that the discontinuity should not appear sharp. In addition, if the sharp boundary is maintained by motion of the galaxy, a large, nearly supersonic velocity of the galaxy with respect to the larger scale halo is required. Spectral analysis indicates that there is a significant jump in elemental abundance across the discontinuity. The X-ray discontinuity is coincident with a radio lobe detected in the NVSS survey. We suggest that the discontinuity in the surface brightness is due to a combination of the near-supersonic expansion of the radio lobe and to turbulent entrainment of high abundance material from the center of the galaxy to more distant regions caused by this expansion.

Throughout this paper, we use $`H_0`$=75 km s-1 Mpc-1, $`\mathrm{\Omega }`$=0.3 and $`\mathrm{\Lambda }`$=0.7. The observed redshift of NGC 507 then corresponds to a distance of 65.5 Mpc, and $`1^{}`$=18.8 kpc.

## 2 Observations

The central region of the nearby (z=0.01646, d=65.5 Mpc) NGC 507 group was observed twice by Chandra. The first observation, taken Oct 11, 2000, was made with the ACIS-S instrument, and the second, taken Jan 8, 2002, with the ACIS-I in very faint mode. The particle backgrounds in both observations were examined for periods of flares, and the average levels compared with the nominal ACIS background rates. We found that the first observation was contaminated by a steady but anomalously high background, so we have used only data from the second observation in all analyses presented in this letter. No flares or periods of high background were found in the second dataset. Very faint mode filtering was applied to this second event file to reduce the background to the lowest possible level. A background event file was created from the standard ACIS background data appropriate for the date of the observation, and very faint mode event filtering applied to these data as well. We also have obtained archival data from the XMM-Newton observation, and have used these data to constrain spectral parameters as described below.

A raw X-ray image in the 0.5 to 2.0 keV band is shown in Figure 1. An adaptively smoothed, background subtracted, exposure corrected image of NGC 507 is shown in the Figure 2. The complex azimuthal structure in the X-ray emission can be seen in both figures. The large scale coronal emission can be seen in Figure 2 with the X-ray peak centered on NGC 507. X-ray emission from the less massive NGC 508 is detected to the N. Some structure in the emission is seen to the west of NGC 507, and is most likely related to interaction between the radio lobe and the ISM.

## 3 Analysis

One of the most striking features of Figure 1 is the sharp discontinuity in the surface brightness approximately $`120^{\prime \prime }`$ to the east of the central peak spanning an angle of approximately $``$125 from the north through the east to the south. The radial surface brightness profile in the 0.5 to 2.0 keV band of the 60 wedge shown in Figure 1 is plotted in Figure 3. The discontinuity appears to be sharpest to the SE, but in creating the surface brightness profile in Figure 3, we chose the wedge to the E and ENE to avoid uncertainties associated with the CCD chip gap toward the SE and NGC 508 to the N.

The temperature of the gas around NGC 507 was determined by two independent methods. First, we divided a 120 wedge (twice the angle of the sector shown in Figure 1 to improve statistics) into eight regions, three to the west of the discontinuity, one straddling the discontinuity, and four to the east. A single temperature absorbed APEC model was fit to both the Chandra and the XMM-Newton data (independently) with the abundance as a free parameter. The absorption was held fixed at the galactic value ($`N_H`$=5.5$`\times `$1020 cm-2). The lower limit of the spectral fits was restricted to 0.8 keV to avoid uncertainties in the ACIS-I response due the build up of contaminants, although a correction also was applied to the response to account for this effect.

The best fit Chandra and XMM-Newton/MOS temperatures and abundances with 90% confidence intervals as a function of distance from the center of NGC 507 are plotted in Figure 4. Figure 4 shows only a small (0.1-0.2 keV) difference in the temperature between the region interior and exterior of the discontinuity. The abundance increase from 0.3$`\pm `$0.1 to 0.6$`\pm `$0.1 towards the center is statistically significant at the 90% confidence level. If the abundance is held constant, the temperature difference between the interior regions and the exterior regions increases by a only small ($``$0.15 keV) amount.

We also created a temperature maps from both data sets using the algorithm described in Churazov et al. (1996). The map created from the XMM-Newton/MOS data is shown in Figure 5 and supports the conclusion of the spectral fitting that there is little change in temperature just across the discontinuity, but a small temperature decrease toward the center. Interestingly the temperature is lowest in the region to the S/SE where the discontinuity appears sharpest.

## 4 Interpretation

If the abundance is roughly constant across the discontinuity, there must be a large pressure jump due to the large density discontinuity because the observed temperature jump across the boundary is small. We consider three possible interpretations to explain the observed surface brightness feature. In the first case, we consider the possibility that the edge is caused by motion of NGC 507 relative to the larger scale group dark matter potential. Second, we consider a scenario where the pressure across the boundary is balanced by an unseen relic radio lobe left over from an earlier (and more powerful) epoch of nuclear activity. Finally, we investigate the possibility that the edge is due to the subsonic expansion of the relatively weak lobe currently observed lying interior to the discontinuity. This expansion has pushed higher elemental abundance gas out from the central region. The appearance of the discontinuity is then a consequence both of the compression of the gas due to expansion of the lobe and to the increase in the abundance across the discontinuity.

One interpretation commonly invoked to explain the sharp discontinuities observed in galaxy clusters is that the central object is moving with a significant velocity relative to the larger halo (Vikhlinin et al., 2001; Markevitch et al., 2001). These are the so-called cluster ”cold fronts”. There are several problems with this scenario for NGC 507. First, the implied velocity must be very large. If the pressure difference between the interior and exterior regions is balanced by ram pressure, the implied velocity is approximately 450 km s-1 (Mach 0.7) relative to the ambient medium. Such a large velocity relative to the dark matter halo is not reasonable since NGC 507 is the largest galaxy in the group and lies at or near the center defined by the outer X-ray isophotes. The large scale gravitating mass of the NGC 507 group can be estimated from the parameters of the $`\beta `$-profile assuming the larger scale gas is in hydrostatic equilibrium. If we consider NGC 507 to be an undamped harmonic oscillator in this dark matter potential moving with a velocity of 450 km s-1, it should reach a distance of $``$75 kpc from the center before falling back. Since there are no other known instances of the largest galaxy in a group or cluster being so far from the center of the gravitational potential, we consider this possibility to be unlikely.

A second significant complication with this interpretation arises from hydrodynamic considerations. The pressure gradient is due largely to a density difference, not a temperature difference. There is no evidence for supersonic motions (i.e. there is no hot, shock heated shell in front of the gas which would be moving supersonically). This implies that the motion is subsonic (consistent with the above analysis), but since the pressure difference is due only to differences in density, the boundary should not appear sharp, as there must be a continuous and gradual pressure and density transition between the two regions. If the pressure difference were due to a temperature discontinuity (such as in A3667 (Vikhlinin et al., 2001)), a sharp boundary could be observable because the surface brightness is proportional to $`n^2`$. It is conceivable that a much smaller region of the galaxy is actually moving, and the discontinuity is material that is piled up in front of the stagnation point. A large velocity would still be required to contain this material.

For the sharp boundary to exist in the absence of motion, additional pressure support could exist beyond the discontinuity. An unseen relic radio lobe could provide the needed pressure support beyond the discontinuity. NGC 507 is a known radio source and has been classified as an FR I radio galaxy (Fanti et al. 1986). Data from two radio observations are shown in Figure 6 overlaid onto an adaptively smoothed Chandra image. These radio observations detected a double lobe structure lying along the east-west axis of the galaxy. The eastern (weaker) radio lobe is contained inside the discontinuity. The western lobe is contained within some of the complex emission to the west and may be responsible for the complex morphology (Forman et al., 2001). No evidence for radio emission beyond the X-ray discontinuity has been reported.

Is it feasible that such an unobserved, relic lobe is providing pressure support beyond the observed discontinuity? Relativistic electrons with $`\gamma `$5$`\times `$102 in a magnetic field of 50 $`\mu `$G (a typical value for the lobes of FR II radio galaxies) would radiate at $``$10 MHz, which is unobservable. The density of such electrons can be estimated by requiring the pressure of the lobe to balance the pressure difference across the discontinuity. Assuming that the lobe has a $``$ 50 kpc radius to cover the observed discontinuity, the total energy in the lobe would be about 5$`\times `$1060 ergs. This number is a factor of a few larger than the total energy in one of the lobes of the FR II radio galaxy Cyg A (Carilli et al., 1991), and would imply an energy deposition rate of approximately 1046 ergs s-1 for 107 years. The luminosity of IC scattered CMB photons would be considerable ($``$109 L), but peaking in the EUV ($``$30 eV), again difficult to observe. The existence of such an energetic but otherwise unobservable radio relic is not ruled out by this simple energy argument, but is highly unlikely.

As a final possibility, we consider a scenario where the inflation/expansion of the known radio lobe is creating the discontinuity. Such a scenario can explain all the observed features, and has recently been investigated theoretically by Reynolds, Heinz, and Begelman (2000) to understand Chandra observations of cavities and cool shells seen around the lobes of several other FR I radio galaxies. We suggest that there is a low surface brightness radio lobe interior to the discontinuity, expanding to the S/SE. Such a lobe is hinted at in the NVSS radio map (see Figure 6). The relationship between the discontinuity and the smaller scale radio lobe presented in Fanti et al. (1986) is not clear. In this scenario, the subsonic expansion of the radio lobe, whether due to ongoing nuclear activity or by buoyancy if the central source has shut off, gently pushes the lower temperature, higher abundance material from the cooling center of the galaxy out to more distant regions (Churazov et al., 2001; Nulsen et al., 2002).

We hypothesize that the sharp appearance of the discontinuity is due to the combination of the motion of the radio lobe and to a significant abundance discontinuity across the interface. At the observed $``$ 1.2-1.3 keV temperature of NGC 507, the emissivity of the gas is a strong function of the heavy element abundance, and there is an abundance gradient from the halo to the core (see Figure 4). Material entrained from the central regions of NGC 507 would have higher abundance than that in the halo and would naturally give a sharp boundary if the emission were due to a cap surrounding a radio lobe. That is, the sharp abundance discontinuity (shown in Figure 4) would naturally produce a discontinuity in the surface brightness, and at the same time maintain pressure equilibrium without the need for a large temperature jump.

To demonstrate the effects of the enhanced elemental abundance interior to the discontinuity, we consider a model for the surface brightness where the emissivity of the gas is determined by the extreme values of temperature and abundance in Figure 4. That is, the gas interior to the discontinuity has an abundance more than twice that of the material beyond the discontinuity, but with a 20% lower temperature. We assume that the gas density is described by a single $`\beta `$-profile, and that there is no pressure discontinuity across the surface brightness discontinuity. The emissivity of the interior material is enhanced by a factor of 2.16 due to its higher abundance and lower temperature. We have overplotted the surface brightness profile from this model onto Figure 3.

Such an idealized model closely approximates the observed surface brightness profile and can account for most, but not all, of the emission interior to the discontinuity. The remaining difference in the surface brightness (between the data and the model) can be attributed to compression of the gas by the inflation of the radio lobe. We have increased the density, and therefore pressure, of the interior region to match the surface brightness profile, and estimate that an expansion velocity of the lobe on the order of 100 km s-1 is required to match the profile in detail and balance the increased pressure. We therefore conclude that a model combining the effect of the enhanced abundances and compression from subsonic motion can easily match the observed surface brightness profile. Further detailed radio observations are required to constrain the relationship between the lobe and this discontinuity.

Chandra has now detected more than a dozen examples of active galaxies ”blowing bubbles” into the ICM (Blanton et al. (2001); McNamara et al. (2000) among many others). In most cases, there is no evidence of shocks being driven into the IGM/ICM. This is clearly the case with NGC 507 as there is no evidence of hotter gas associated with either the discontinuity or any of the radio components. The inflation of the lobe has however disrupted the gas in the central region and may provide an efficient mechanism to mix the higher abundance material at the center with the lower abundance material (presumably much older) in the halo.

## 5 Conclusions

We have observed a sharp discontinuity in the X-ray surface brightness profile in a Chandra observation of NGC 507. We conclude that this discontinuity is most likely not a ”cold front” as has been observed in many clusters, and this is not due to non-hydrostatic motion of the group core. It is also highly improbable that this discontinuity is supported by an unseen, relic radio lobe from an earlier epoch of galaxy activity. The appearance of this discontinuity is probably due to two factors. First, the subsonic expansion/inflation the radio lobe is compressing the ambient ISM. Second, there is an abundance gradient across the boundary that creates a discontinuity in the surface brightness. We suggest that the inflation of the radio lobe has pushed up and/or entrained higher abundance material from the central regions within NGC 507 into the group halo. It is unclear whether this is a common phenomenon or not because it is most readily observable only in relatively cool ($`k_BT<`$ 2 keV) galaxies and clusters where the emissivity of the gas is a strong function of the abundance. This type of disruption/mixing of the central, high abundance region of galaxy groups and clusters by jets/lobes may be common, but would be difficult to detect in high temperature systems where abundance variations would not produce changes in surface brightness but instead would require detailed abundance maps for their detection. There is at least one other object, A576, that shows evidence for a surface brightness discontinuity that is at least partially attributable to a gradient in the elemental abundance (Kempner and David, in preparation). If this is indeed a common phenomenon, the infall velocity of cool subclusters or the expansion velocity of radio lobe/ICM interactions would be overestimated, perhaps by as much as a factor of a few.

This work was supported by NASA contracts NAS8-38248, NAS8-39073, the Chandra X-ray Center, and the Smithsonian Institution.


text:
# Observations of H i absorbing gas in compact radio sources at cosmological redshifts

## 1 Introduction

The energy output from AGN is commonly agreed to be powered by accretion onto a super-massive black hole, resulting from the infall of gas transported from the host galaxy down to the central regions (e.g., Rees rees84 (1984)). Recent detailed and statistical studies of Seyfert galaxies (e.g., Cid Fernandes et al. cid01 (2001)) are beginning to chart the evolutionary sequence that leads to the production of an active galactic nuclei in late-type (Sd through S0) galaxies and the effect of that evolution on the host galaxy (e.g., Storchi-Bergmann et al. storchi01 (2001)). The picture that is emerging is one of a close tie between the galaxy environment and galaxy-galaxy interactions leading to a temporal sequence of circumnuclear starbursts and accretion onto the nuclear black hole (e.g., Sanders et al. sanders88 (1988)).

However, while this picture is becoming well defined for Seyferts, very little is known regarding the situation in powerful radio galaxies. The mostly radio-quiet AGN in a large nearby galaxy sample studied by Ho et al. (ho97 (1997)) often reside in early-type disk galaxy hosts, many of which have spiral structures, and the amount of gas needed to sustain the low-luminosity activity for the expected lifetime is readily available. But looking at radio-loud AGN, Martel et al. (martel99 (1999)) have found that at least 89% of the hosts of 46 low-redshift ($`z<0.1`$) 3CR radio galaxies are ellipticals. At larger redshifts, Dunlop et al. (dunlop03 (2003)) have shown that the host galaxies of radio-quiet and radio-loud quasars, as well as radio galaxies, are massive ellipticals. These may of course themselves be merger products of late type galaxies. True spiral hosts of powerful radio galaxies seem to be very rare (see also Ledlow et al. ledlow98 (1998)).

The predominance of elliptical hosts is often thought to be puzzling, since these are usually assumed to be gas poor. It should be noted first that major galaxy merger products are likely places where gas has lost angular momentum and can feed into the central engine. Apart from that, Walsh et al. (walsh89 (1989)) have shown that in ellipticals the amount of long-wavelength IR emission, indicative of the presence of dust and gas, is correlated with the occurrence and strength of a central radio source. More recently, it has been shown that early-type galaxies can contain gas in both molecular form (see e.g. Rupen rupen97 (1997); Knapp & Rupen knapp96 (1996)), and in atomic form (see e.g. Bregman et al. bregman92 (1992); Oosterloo et al. oosterloo99 (1999)), and this gas could act as a fuel reservoir for the central engine. Indeed, central structures of gas are frequently observed in nearby radio sources; for instance, 17 of 19 nearby FRI radio galaxies observed with the HST by Verdoes Kleijn et al. (verdoes99 (1999)) have dust structures in their centres. Such dust disks are also reported in molecular (e.g in 4C12.50, Evans et al. evans99 (1999)) and atomic lines, e.g. in NGC 4261 (Van Langevelde et al. vanlangevelde00 (2000)) and NGC 1052 (Vermeulen et al. 2003a ). Furthermore, observations of nearby radio-loud AGN by van Gorkom et al. (vangorkom89 (1989)) have directly shown, in about 30% of the galaxies studied, the presence of cold gas visible through the H i $`\lambda =21`$ cm line of neutral hydrogen in absorption towards the central radio source. In a recent study, Morganti et al. (morganti01 (2001)) have found evidence that in their radio galaxy sample the H i detections are related to circumnuclear tori.

At cosmological redshifts, the galactic medium is largely unexplored. But the recent upgrade of the Westerbork Synthesis Radio Telescopes (WSRT) has opened up a new avenue of research. The array of 14 telescopes, with its new multi-channel correlator (DZB), offers good sensitivity as well as a superior capability to reject external interference compared to single dish telescopes. The wide band UHF-high receivers, spanning 725–1200 MHz, enable studies of the occurrence and kinematics of atomic gas seen by means of the redshifted H i line in absorption against radio sources almost out to $`z=1`$. We here present the combined results of a number of H i line searches performed at the WSRT. These projects had a variety of goals, all related to the properties of cold gas in compact radio sources with a total extent less than a few arcseconds, within the appropriate redshift range, in the Northern sky, and with a flux density above a few hundred mJy (in order to have a reasonable opacity detection threshold). A majority of the sources in this paper are classified as Gigahertz Peaked Spectrum (GPS) or Compact Steep Spectrum (CSS) sources, believed to be intrinsically compact (e.g. Fanti et al. fanti95 (1995), O’Dea odea98 (1998)). Further, a subset of the sources in this paper are end-on classical doubles (including both quasars and galaxies). As a result of coordination between the present authors, we here present the complete observational results for the sources in the original programmes for which usable spectra were obtained, plus some statistics on the absorber depth and kinematics. In a follow-up paper, using a more tightly defined sample of GPS and CSS sources including data from the literature, Pihlström et al. (pihlstrom03 (2003)) discuss H i absorption in the context of other properties, and show, for example, the existence of an anti-correlation between absorption depth and linear size of the radio source.

## 2 WSRT observations and data analysis

This paper presents the combined results of a number of different WSRT projects, performed between 1997 and 2001, using the UHF-high receivers. While there were slight variations in the observational and data reduction methods used, the description below adequately describes the essentials. The data reduction was carried out using the NRAO AIPS and Caltech DIFMAP packages, typically switching back-and-forth to make use of their respective strong points.

Typically, each target was first observed for a few hours, in each of two orthogonal linear polarisations, with 128 spectral channels covering a 10 MHz wide observing band (the maximum available when the surveys started), centred at the frequencies predicted for H i based on the optical redshifts (these have some uncertainty, as discussed in Section 3). Doppler tracking was used; it ensures that spectra taken separately at arbitrary times of the day and year can be directly combined pixel-by-pixel. Most spectra were Hanning smoothed online, chiefly to combat the spectral rippling effects of strong narrow-band interference somewhere in the observing band. As a result, the effective spectral resolution in the initial observations is about 150 kHz.

The combined initial target list incorporated about 80 sources. Despite the WSRT’s impressive ability to suppress external interference, by virtue of its being an interferometric array, and having a multi-channel correlator, we have encountered a number of frequency intervals in which the level of the interference was always too high to obtain any useful data, even after trying to observe at different times, with different telescope orientations (hour angles), or with slightly shifted frequency bands. In the end, we pursued the observations for 57 sources to the point where the spectra showed either a line detection or gave a reasonable upper limit. This often involved follow-up observations done with narrower, offset bandwidths and/or twice the number of frequency channels (at the price of fewer correlated baselines) in the pursuit of line candidates.

Bright external sources such as 3C48, 3C147, and 3C286 were used as calibrators. Time-limited or baseline specific interference was first removed, before using them for standard complex gain and bandpass-calibration. No attempt was made to obtain polarisation information; the data correlated in the two orthogonal linear polarisations were simply added after separate gain calibration, to increase the signal-to-noise of the final spectra. Furthermore, since our primary interest was on line optical depths, no particular care was expended on setting the overall flux density scale, but we believe that the flux densities shown in Figures 1 and 2 are accurate to better than 5%.

After external calibration, interference was removed from the target spectra to the extent possible. Figures 1 and 2 show that this was not always completely successful: some interfering signals, particularly when their strength did not dominate in the observing band, could not be adequately discriminated by their temporal or spatial characteristics. After editing to the extent possible, those spectral channels without any line features were averaged, and this continuum dataset was then used in iterative cycles of self-calibration alternated with image cleaning and/or model-fitting to the uv-data. The compact targets are all unresolved to the WSRT. Characterisation of the continuum sources in the field was not a goal in itself, and was pursued only to the level needed to find adequate self-calibration gain factors, and to ensure that confusing components would not disturb our spectral results. The self-calibration complex gains found for the continuum were subsequently applied to all spectral line channels. Final spectra were then produced by coherently integrating all phase-calibrated visibilities per spectral channel; this method produces an optimal signal-to-noise ratio given that in all cases the target, at the phase centre, is by far the dominant source in the field.

We find H i absorption to be present in the spectra of 19 of the 57 sources (33%); these are shown in Fig. 1. Table 1 gives the particulars of the detected absorption lines, determined by Gaussian profile fitting of the peak optical depths, $`\tau _{\mathrm{peak}}`$, and FWHM line widths, $`\mathrm{\Delta }V`$; the column densities were computed from these assuming uniform coverage of the entire radio source and a spin temperature $`T_{\mathrm{sp}}=100`$ K, using $`N_{\mathrm{HI}}=1.82\times 10^{18}T_{\mathrm{sp}}\tau _{\mathrm{peak}}\mathrm{\Delta }V`$. Both assumptions mean that the listed column depths are lower limits: the covering factor might be below unity, and if the H i absorption were to arise in the pc-scale vicinity of the AGN, conditions might be such that $`T_{\mathrm{sp}}8000`$ K (Maloney et al. maloney96 (1996)). For the other 38 sources, we show the spectra in Fig. 2. For Table 2, we have used a FWHM of 100 km s-1 as a plausible width for an absorption line in these objects (see Section 4.1 and Fig. 5), in order to be able to derive $`2\sigma `$ upper limits to the line depths and H i column depths, based on the continuum flux densities and r.m.s. noise levels of the spectra.

## 3 Detection uncertainties

We expect that many of the optical redshifts used to centre the initial observing bands are uncertain at the level of a few hundred km s-1; particularly at the higher redshifts, they are typically based on emission lines which can be quite broad, and may have a velocity centroid which is offset from the systemic velocity of the host galaxy. With a 10 MHz wide band, we have spanned a velocity range of approximately $`\pm 1250`$ km s-1 at $`z=0.2`$, increasing to approximately $`\pm 2000`$ km s-1 at $`z=0.85`$. Also note that an uncertainty of $`\pm 0.001`$ in $`z`$, which may be typical for the optical redshifts in the literature, corresponds to a frequency uncertainty for H i of nearly $`\pm 1`$ MHz at $`z=0.2`$ but only about $`\pm 0.4`$ MHz at $`z=0.85`$. Thus, (optical) redshift uncertainties give the greatest (radio) frequency uncertainties at the lowest redshifts. In view also of interpreting the offsets discussed in Section 4.2, it would be valuable to obtain accurate systemic redshifts from high sensitivity, high resolution optical spectroscopy, for all of the sources we have detected in H i. Nevertheless, Fig. 6, discussed in Section 4.2, shows that in the 19 sources with detections, the H i velocity distribution around the optical redshift tails off rapidly towards $`\pm 500`$ km s-1. We therefore believe that our spanned bandwidth was adequate to cover the vast majority of the potential lines, and quite likely all of them.

In view of the large range in flux densities amongst the targets, it was impractical to achieve a uniform sensitivity to line optical depth (opacity). It is difficult to quantify what this means for an individual non-detection, but the following limited statistical analysis is indicative. Figure 3 shows the individual $`2\sigma `$ upper limits plotted along with the line opacities for the detected lines, while Fig. 4 shows their cumulative distributions. Firstly, we point out that from Fig. 3 there appears to be no evidence for a strong correlation of optical depth with FWHM, which means that adopting a somewhat arbitrary width of 100 km s-1 in order to compute upper limits is reasonable (the limits would have been somewhat tighter had we adopted a larger width). It is also clear that the distribution of detected line strengths is different from that of the upper limits: many of the limits are below most of the actual detections, meaning that these limits are tight enough to be very significant. A KS test shows that the cumulative distribution of the detections is different from that of the 2$`\sigma `$ limits at the 98% confidence level. It is also clear from Fig. 3, however, that about half of the 2$`\sigma `$ limits leave something to be desired: while they rule out the presence of prominent absorption (roughly: above $`\tau =1`$%), they still allow the presence of shallower lines analoguous to those seen in some of the weaker detections (roughly: below $`\tau =0.5`$%). A significant investment in telescope time would be needed in most cases in order to constrain this potential incompleteness.

Nevertheless, the detections that we do have are already an important next step in studying the cold neutral medium in cosmologically distant radio sources. Below, we briefly discuss some of the line properties. In a follow-up paper (Pihlström et al. pihlstrom03 (2003)), the line properties are related to other source parameters. VLBI observations in order to image the neutral gas distribution are already being pursued for some of the most prominent absorbers (e.g., Vermeulen rcviau (2003), Vermeulen et al. 2003b ).

## 4 Line Properties

### 4.1 Line widths

Figure 5 shows the distribution of the FWHM of the Gaussians fitted to the absorption lines. This plot confirms that H i lines in these sources often are broad, with a FWHM of 100–200 km s-1 or even more. We find that the primary lines (i.e. the deepest Gaussians, with the highest optical depth in each source) are probably somewhat broader than the secondary lines; the mean FWHMs, of 156 km s-1and 72 km s-1, respectively, differ at a confidence level of 97%, according to a Student’s $`t`$-test.

Often, H i absorption is interpreted as being due to a circumnuclear disk; for example VLA 21 cm H i absorption imaging of a sample of Seyfert and starbursts displayed almost exclusively sub kpc scale rotating disks aligned with the outer galaxy disk (Gallimore et al. gallimore99 (1999)). Data for radio loud sources are also consistent with disks on sub kpc scales, for example in Cygnus A (Conway conway99 (1999)) and 1946+708 (Peck & Taylor peck01 (2001)). Those examples all have line widths of the order of 100 km s-1, comparable to the detections reported in this paper. If the wide lines arise on sub kpc scales close to the nucleus, the narrower lines could be additional absorption occurring on larger scales within the galaxy, provided that the line of sight does not cross a substantial amount of differential rotation. Alternatively, the narrower lines could be from single cohesive clouds, with only moderate internal velocity dispersion, as is observed in 3C236 (Conway & Schilizzi conway00 (2000)).

### 4.2 Inflowing and outflowing gas ?

It was pointed out by van Gorkom et al. (1989) that in all 4 out of 29 nearby radio galaxies ($`z<0.13`$) with detected H i absorption, and also in 4 other nearby systems, the H i velocity is positive, i.e. infalling with respect to the systemic velocity of the galaxy. The radio structure (or, in the case of 3C236, at least the absorbed part) of those nearby galaxies with H i absorption is very compact, and many of these objects are now classified as GPS or CSS sources. Thus, they potentially match well with the types of sources that constitute the majority of our sample.

As Figure 6 shows, we now find H i lines at both positive and negative velocities with respect to the optical redshift. Also, while typical velocities in the lower redshift sources of van Gorkom (1989) were in the range 0 to $`+200`$ km s-1, with only 1 value around $`+400`$ km s-1, 10 of our 19 main (highest opacity) lines have a velocity more than 200 km s-1 different from the optical redshift; the main line velocities range between $`v=1420`$ km s-1 and $`v=+318`$ km s-1. While, as already discussed in Section 3, inaccuracies in the optical redshift may play a role in causing some of the scatter, we believe the larger H i velocities are likely to be significant. We see an intriguing hint that the velocity distribution may be skewed towards outflow rather than infall. This effect can be stated more dramatically than it appears through the binning of Fig. 6: 13 of the main line velocities are negative as compared to 6 positive; of these, 7 have $`v<200`$ km s-1 as compared to 3 having $`v>200`$ km s-1, and of these, 3 have $`v<600`$ km s-1 as compared to none at $`v>600`$ km s-1. The mean centroid velocities are $`v=223`$ km s-1 for the main lines and $`v=58`$ km s-1 for the secondary lines, respectively. However, Student’s $`t`$ tests at a confidence level of 95% show that the mean velocities of the main lines and the secondary lines indeed do not differ significantly from zero, and also not from each other, possibly only because the sample size is too small.

Thus, we find the opposite of the results obtained for nearby radio sources by van Gorkom et al. (vangorkom89 (1989)) in both the sign and the magnitude of the velocities. Speculatively, we think that in these higher luminosity sources, interaction between the radio jets and the surrounding inner galactic medium could lead to significant motions of the gas, such as found in 3C236 by Conway & Schilizzi (conway00 (2000)). Other evidence for such interactions is often seen as the so-called radio-optical alignment effect. The kinematics of the aligned gas seen in CSS sources shows velocity offsets of a few hundred km s-1 (both positive and negative) with respect to the nuclear velocity (e.g., O’Dea et al. odea02 (2002)). Alternatively, in some cases we may be looking at the line-of-sight rotational motion of a neutral component of broad-line gas, seen in front of a parsec-scale radio source, as has been discussed for lower redshift galaxies for example by Morganti et al. (morganti01 (2001)).

## 5 Summary

H i absorption searches have been carried out with the WSRT for a sample of mostly compact radio sources with redshifts out to nearly $`z=1`$, giving observing frequencies which fall in the UHF-high band (725–1200 MHz). Out of 57 sources for which adequate spectra could be obtained, 19 are found to have associated H i in absorption, with opacities ranging from 16% down to 0.2%; below 1% opacity, our statistics may well be incomplete. The typical line width is $`150`$ km s-1, but the lines can also be several times narrower or wider. Column densities ranging from $`1\times 10^{19}`$ to $`3.3\times 10^{21}`$ are obtained by assuming uniform coverage of the radio source and a spin temperature $`T_{\mathrm{sp}}=100`$ K. In contrast to nearby radio sources, we find not only positive but even more negative H i velocities, up to more than 1000 km s-1(as compared to the optical redshifts). Perhaps, inflow and outflow can both occur, or there is interaction with the radio source, or we are seeing rotating broad-line gas.

###### Acknowledgements.

The WSRT is operated by ASTRON (The Netherlands Foundation for Research in Astronomy) with support from the Netherlands Foundation for Scientific Research (NWO). We are very grateful for the dedication of the ASTRON staff who worked on the WSRT upgrade, as well as those who assisted in performing the survey observations. This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration in the United States of America.


text:
# Superfluid response in monolayer high-𝑇_𝑐 cuprates

## Abstract

We have studied the doping dependence of the in-plane and out-of-plane superfluid density, $`\rho ^s(0)`$, of two monolayer high-$`T_c`$ superconductors, HgBa2CuO4+δ and La2-xSrxCuO4, using the low frequency $`ac`$-susceptibility and the muon spin relaxation techniques. For both superconductors, $`\rho ^s(0)`$ increases rapidly with doping in the under- and optimally doped regime and becomes nearly doping independent above a critical doping, $`p_c0.20`$.

Measurements of the magnetic penetration depth have been important in probing the order parameter and in testing theories of high-$`T_c`$ superconductors (HTS) . In hole doped HTS, the low temperature dependence of the in-plane penetration depth $`\lambda _{ab}(T)`$ is linear and doping independent, indicating the presence of nodes in the superconducting energy gap . The $`c`$-axis penetration depth $`\lambda _c`$ is a key parameter for some theories describing the mechanism of high temperature superconductivity . It is sensitive to the electromagnetic anisotropy of the system and has been used to test the pairing symmetry and properties of interlayer dynamics .

In a recent study of the spin and charge response of HTS, it was found that both the superfluid density $`\rho ^s(0)`$ $``$ $`\lambda ^2(0)`$ and the muon spin relaxation, $`\mu `$SR, rate show dramatic changes at a critical doping $`p_c0.20`$, slightly above optimal doping, in pure and Zn-doped La2-xSrxCuO4 (La-214) and Bi2.1Sr1.9Ca1-xYxCu2O8+y (Bi-2212) at zero temperature . The sharp changes in the superfluid density with the disappearance of a spin glass phase transition near $`p_c`$ suggested a change in symmetry of the ground state. The existence of such a special doping has been demonstrated in many other physical quantities  and the $`\rho ^s(0)`$ and $`\mu `$SR data could be linked to the presence of a quantum phase transition at $`p_c`$, that is in turn related to the opening of the normal state pseudogap.

To elucidate further the changes in the ground state across the phase diagram of HTS we have studied the doping dependence of the zero-temperature in-plane and out-of-plane superfluid responses, $`\rho _{ab}^s`$ and $`\rho _c^s`$, for two monolayer HTS materials: HgBa2CuO4+δ (Hg-1201) and La-214. This study allows us to determine both the in- and out-of-plane responses as a function of doping and to perform a direct comparison between two simple HTS with different degrees of disorder. We find that in both systems, the superfluid density is strongly doping dependent below $`p_c`$ and shows abrupt changes around $`p_c`$. For Hg-1201 the effect is sharper and there is actually a peak in the superfluid density at $`p_c`$.

The Hg-1201 samples were prepared in Houston using a method similar to that described in ref. . Their doping level can be continuously varied from very under- to heavily over-doped regime by adding or removing oxygen. Unlike La-214, where the doping is varied by Sr substitution for La, which may cause pronounced disorder effects, the variation of oxygen concentration in Hg-1201 is known to induce little lattice disorder . The Hg-1201 samples were characterised using magnetisation and thermoelectric power measurements. The doping level $`p`$ was determined by both thermopower  and the universal relation $`T_c`$=$`T_{c,max}`$\[1-82.6($`p`$-0.16)2\] . The La-124 samples were synthesised in Cambridge using solid-state reaction followed by oxygenation. Effort was taken to ensure high purity and homogeneity. All La-214 powders were dried, reacted, ground, milled, re-pressed and re-sintered at least four times. The phase purity was verified by powder x-ray diffraction as well as extensive transport and thermodynamic measurements. No signal of impurities or inhomogeneity was captured in micro-analytical spectroscopic studies . The $`T_c`$ values as well as lattice parameters of these samples were in good agreement with published data. In La-214 $`p`$ is taken as equal to the Sr concentration. The heat capacity anomalies and $`ac`$ -susceptibility transitions are sharp.

We have measured the magnetic penetration depth $`\lambda `$ using the low-field $`ac`$-susceptibility technique for grain-aligned powders . The superfluid density is inversely proportional to the square of the in-plane penetration depth. To determine the in-plane and $`c`$-axis penetration depths separately, the grains were magnetically aligned in a static field of 12T at room temperature. X-ray powder diffraction scans  for both La-214 and Hg-1201 samples showed that more than $`90\%`$ of the grains had their CuO2 planes aligned. The $`ac`$-susceptibility measurements were performed down to 1.2K with a home-made susceptometer using miniature coils with $`H_{ac}`$=1-3G rms at $`f`$=333Hz. The absence of weak links among grains was confirmed by the linear response of the signal with $`H_{ac}`$ from 0.3 to 3G rms and $`f`$ from 33 to 333Hz. We also used a commercial susceptometer to confirm some of our findings. Taking the grains to be approximately spherical, as indicated by scanning electron microscopy, the data were analysed using London’s model . The $`ac`$-susceptibility data were also confirmed by standard transverse field $`\mu `$SR experiments performed on unaligned powders at 400G .

Figure 1 shows the data for (a) $`T_c`$, $`\lambda _{ab}^2(0)`$ and (b) $`\lambda _c^2(0)`$ for La-214. The $`T_c`$ and $`\lambda _{ab}^2(0)`$ data were published in an earlier paper and are included here for comparison . $`\rho _{ab}^s(0)`$ is nearly doping independent in the overdoped regime, but drops fast below $`p=0.19`$. This suppression of the superfluid density for $`p<0.19`$ was previously discussed in terms of a competition between quasi-static magnetic correlations and superconductivity . It has also been linked to the strong reduction in entropy as well as condensation energy associated with the opening of the normal state pseudogap .

$`\rho _c^s(0)`$ shows similar behaviour as its in-plane counterpart. However, in contrast to the nearly linear doping dependence of $`\rho _{ab}^s(0)`$ on $`p`$, $`\rho _c^s(0)`$ shows a stronger doping dependence below $`p_c`$ corresponding to $`1/\lambda _c^2`$ $``$ $`p^n`$ with $`n2.7`$. This difference in the doping dependence between $`\rho _{ab}^s(0)`$ and $`\rho _c^s(0)`$ is probably associated with the unconventional interlayer coupling of electrons in high-$`T_c`$ oxides, and is worthy of further theoretical and experimental investigation.

Figure 2(a) shows the doping dependence of $`T_c`$ and $`\rho _{ab}^s(0)`$ for Hg-1201. Similar to La-214, $`\rho _{ab}^s(0)`$ is relatively doping independent in the overdoped regime and shows a sharp drop below 0.19. A similar $`p`$ dependence of $`\rho _{ab}^s`$ has been found for Bi-2212 , and recently also reported for Y0.8Ca0.2Ba2Cu3O7-δ (Ca:Y-123) and $`Tl_{0.5y}Pb_{0.5+y}Sr_2Ca_{1x}Y_xCu_2O_7`$ (Pb:Tl-2212) . The maximum of $`\rho _{ab}^s(0)`$ is located at $`p_c`$ for all high-$`T_c`$ compounds. It suggests that the observed doping dependence of $`\rho _{ab}^s(0)`$ below $`p_c`$ is common to all HTS compounds and is not due to a structural transition or inhomogeneity. The relatively doping independent $`\rho _{ab}^s(0)`$ for $`p>p_c`$ in La-214 and Hg-1201 is in agreement to Bi-2212  but seems to differ from the data for $`Tl_2Ba_2CuO_{6+\delta }`$, Ca:Y-123 and Pb:Tl-2212 . The mechanism causing this difference is unknown and is certainly worth further investigation. Nevertheless, it is clear that the maximum of $`\rho _{ab}^s(0)`$ is located at $`p_c`$ for all high-$`T_c`$ cuprates.

Figure 2(b) shows the doping dependence of the $`c`$-axis superfluid density for Hg-1201. A sharp change from large to low superfluid response is also observed around $`p_c`$. This is the sharpest change in $`\rho _c^s(0)`$ ever being reported and together with the observed peak at $`p_c`$ could be related to its tetragonal crystal structure and the fact that Hg-1201 is more ordered than La-214. It is worth noting that a significantly weaker glassy response has been observed in Hg-1201 . We may speculate that this observation suggests that the sharper changes near $`p_c`$ may be linked to a quantum critical point for which disorder causes smoothing of the doping dependence of various physical properties and associated phase transitions.

The interlayer distance between the CuO2 planes may be a key parameter for optimal $`T_c`$. This has been emphasised by Uemura recently . Indeed, for the same in-plane superfluid density, $`T_c`$ is higher if the interlayer distance is shorter. Therefore, the interlayer coupling seems to be essential for obtaining higher $`T_c`$. The observed variation of $`T_c`$ cannot be explained by the simple Kosterlitz-Thouless transition where $`T_c`$ is solely determined by the 2D superfluid density. The similar doping dependence of $`\rho _c^s(0)`$ to $`\rho _{ab}^s(0)`$ observed here supports this view and indicates the fundamental role of the $`c`$-axis electrodynamics to the overall superconducting condensation. As a matter of fact, $`\lambda _c(0)`$ for both monolayer cuprates studied here is small above $`p_c`$ and agrees, for example, with the interlayer tunneling model of Anderson and Chakravarty . Large superfluid response above $`p_c`$ seems to occur in connection with a crossover from two-dimensional to three-dimensional transport, as suggested by the doping dependence of the anisotropy in $`\lambda `$ (Fig. 3) and the associated behaviour of the anisotropy of the normal state resistivity .

In summary, for the two monolayer high-$`T_c`$ cuprates, La-214 and Hg-1201, both the in-plane and $`c`$-axis superfluid response remain relatively constant above $`p_c`$ but drop rapidly below $`p_c`$. We have found a peak in $`\rho ^s(0)`$ at $`p_c`$ for Hg-1201 indicating the strongest superconductivity at the point where the spin glass phase transition (the glass transition temperature $`T_g`$ vesus $`p`$ curve) vanishes and the normal state gap extrapolates to zero . The rapid change and peak may be due to a change in the superconducting ground state. Furthermore, we have observed that the doping dependence of $`\rho _c^s(0)`$ in La-214 follows a power law of approximately 2.7.

C.P. thanks S. Chakravarty and J. W. Loram for enlightening discussions, D.N. Basov for an earlier collaboration and discussions on the subject and The Royal Society for financial support. T.X. acknowledges support from the National Natural Science Foundation of China.