text:
# Directional Emission from a Microdisk Resonator with a Linear Defect

## I Introduction

The idea to use a microdisk geometry as an alternative to the Fabry-Perot cavity in a resonator design for a semiconductor laser was introduced a decade ago. The advantage of this geometry is that the losses for the whispering-gallery modes of a circular resonator are governed by evanescent leakage and, thus, can be very low. Namely, for a mode with a maximal angular momentum $`M=nk_0R`$, where the $`n`$ is the effective refraction index, $`R`$ is the resonator radius, and $`k_0`$ is the wave number of the radiation, the quality factor, $`Q`$, with exponential accuracy is given by

$$\mathrm{ln}Q=2k_0R\left[n\mathrm{ln}\left(n+\sqrt{n^21}\right)\sqrt{n^21}\right].$$
(1)
The value of the effective refraction index is determined by the disk thickness and the indexes of an active and surrounding passive layers. In the pioneering paper Ref.  the effective index was $`n2`$, while $`k_0R`$ for the smallest microdisk was $`6`$. Then Eq. (1) yields $`Q510^5`$. Experimentally measured values of $`Q`$ are much smaller, $`Q150`$, Ref. . The discrepancy is partially due to a prefactor neglected in Eq. (1), but primarily due to the absorption in the active layer. With such a high $`Q`$-value the lasing threshold for a microdisk resonator is very low. For the same reason the output power is also low, which is not desirable. Another serious drawback of the microdisk geometry is that the angular dependence of the output intensity is $`I(\psi )\mathrm{cos}^2(M\psi )`$, whereas applications require a directed emission. In order to remedy these drawbacks two proposals were put forward

(i) to extract the light out of the resonator by using two parallel disks. The first disk with high $`Q`$ contains a multiple quantum well structure in which the light is generated. The second passive disk coupled to the laser contains an opening serving as a leakage source. The shape of the opening determines the directionality of the output light.

(ii) to couple the light out by introducing either an identation in the form of the “tip of the egg” or corrugation on the circumference of the disk.

A radical solution for increasing the output, and, to a certain extent, directionality, by deforming the shape of the disk seem to devaluate the attempts to extract light from a perfectly circular microdisk. This solution relied on the fact that deformation causes a qualitative change in the light-ray dynamics, so that the whispering-gallery trajectory of a ray becomes unstable. As a result, the ray eventually impinges on the boundary at an angle smaller than the critical angle, $`\mathrm{sin}^1(1/n)`$. This leads to a refractive escape. The improvement of the directionality of the output light from a wave-chaotic resonator was studied theoretically in a great detail. The results of calculations for both “bouncing ball” and “bow-tie” modes and $`nk_0R100`$ can be roughly summarized as follows. In each $`90^{}`$-quadrant the output light is concentrated within total angular interval of about $`60^{}`$ with a strong peak of a width $`30^{}`$ and a large number of narrow satellites.

In the present paper we suggest an alternative approach for improving both the directionality and the output efficiency of a circular microdisk. This improvement can be achieved by introducing a properly oriented linear defect away from the circumference. Proposed geometry is illustrated in Fig. 1. The reason why the linear defect causes directional emission from a microdisk is the following. The field of a whispering-gallery mode “tunnels” towards the defect line, which then assumes a role of the secondary source. Since the source is extended, it emits a secondary light beam which is weakly divergent. The divergence is minimal when this secondary light beam is emitted in the radial direction, i.e. in the direction normal to the disk boundary. It is convenient to characterize the position and orientation of the defect by two parameters, namely $`r_0k_0^1`$ \- radial distance from the edge to the circumference, and $`d`$ \- the minimal distance from the defect line to the disk center. As it will be shown below, the optimal orientation of the defect, for which the direction of the secondary beam is radial, is determined by the condition $`d=(Rr_0)/\sqrt{2}`$. Under this condition the directionality of the output light is maximal. Below we will demonstrate that, with exponential accuracy, the radiative losses caused by the defect are given by

$$\mathrm{ln}Q=\frac{2^{5/2}}{3}\left(\frac{r_0}{R}\right)^{3/2}(nk_0R).$$
(2)
These losses dominate over the evanescent losses Eq. (1) if $`r_0R`$. The angular dependence of the defect-induced emission is a lorentzian, which under the optimal condition $`d=(Rr_0)/\sqrt{2}`$, has the form

$$I(\psi )=\frac{1}{\left(\psi \frac{\pi }{4}\right)^2+2n^2\left(\frac{r_0}{R}\right)},$$
(3)
with the width which is also governed by the ratio $`r_0/R`$. Note, that although Eqs. (2), (3) apply only for $`k_0r_01`$, this ratio can still be quite small as long as $`k_0R`$ is large.

The paper is organized as follows. In Sec. 2 we derive Eqs. (2), (3) within the scalar diffraction theory. In Sec. 3 we discuss the limits of applicability of the theory and provide numerical estimates.

## II Angular Dependence of the Output Light

Neglecting the difference between TE and TM polarizations, the field of a whispering-gallery mode in a microdisk represents a solution of the two-dimensional Helmholtz equation

$$_M(\rho ,\varphi )\mathrm{cos}\left(M\varphi \right)J_M(nk_0\rho ),$$
(4)
where $`\rho `$ and $`\varphi `$ are polar coordinates, $`M`$ is the angular momentum, and $`J_M`$ is the Bessel function. We assume that $`M`$ is close to the maximal value $`nk_0R`$. Then the field Eq. (4) is localized at the boundary $`\rho =R`$ within a narrow ring of a width $`\delta \rho R/(nk_0R)^{2/3}R`$. At smaller $`\rho `$ the field falls off towards the center of the disk as

$$_M\mathrm{cos}\left(M\varphi \right)\mathrm{exp}\left[\frac{2^{3/2}}{3M^{1/2}}\left(nk_0r\right)^{3/2}\right],$$
(5)
where $`r=R\rho `$ is the distance from the boundary.

Within the scalar diffraction theory the emitted field caused by the presence of a defect is determined by the Fresnel-Kirchhoff diffraction integral in which the source is the field $`_M(\rho ,\varphi )`$ taken at $`\rho =\rho (\varphi )`$, where $`\rho (\varphi )`$ describes the defect profile. In the case of a linear defect (Fig. 1) we have $`\rho (\varphi )=d/\mathrm{cos}\varphi `$. It is convenient to introduce instead of $`\rho `$ a variable $`x`$ which is the distance along the defect (Fig. 1). The relation between $`\rho `$ and $`x`$ is the following

$$\rho =\left[d^2+\left(\sqrt{(Rr_0)^2d^2}x\right)^2\right]^{1/2}=Rr_0x\frac{\sqrt{(Rr_0)^2d^2}}{Rr_0}.$$
(6)
In the second equality we have used the fact that $`xR`$. Substituting Eq. (6) into Eq. (5) we obtain

$$_M(x)\mathrm{exp}\left[\frac{2^{3/2}}{3}\left(\frac{r_0}{R}\right)^{3/2}nk_0R\right]\mathrm{exp}(ax)\mathrm{cos}[M\varphi (x)],$$
(7)
where

$$a=2^{1/2}nk_0R\left(\frac{r_0}{R^3}\right)^{1/2}\frac{\sqrt{(Rr_0)^2d^2}}{Rr_0}.$$
(8)
In Eq. (7) we assumed that $`xr_0`$. Indeed, the relevant values of $`x`$ are $`a^1`$. Then the condition $`xr_0`$ can be rewritten as $`r_0a(nk_0R)(r_0/R)^{3/2}1`$. We see that this condition is equivalent to the requirement that the asymptotic Eq. (5) is valid at $`r=r_0`$. The first $`x`$-independent factor in Eq. (7) determines the dependence of the output field on the defect position, $`r_0`$. The expression Eq. (2) immediately follows from this dependence.

The form of the function $`\varphi (x)`$ in Eq. (7) can be easily established from Fig. 1

$$\mathrm{tan}\varphi =\frac{\sqrt{(Rr_0)^2d^2}x}{d}.$$
(9)
Now we are in position to write the expression for the intensity of the outgoing light in the direction $`\psi `$. It is given by the following integral along the defect

$$I(\psi )\left|_0^{\mathrm{}}𝑑xe^{ax}\mathrm{cos}\left[M\varphi (x)\right]_\pi ^\pi 𝑑\phi \mathrm{exp}\left[ink_0b(x,\phi )ik_0R\mathrm{cos}(\psi \phi )\right]\right|^2.$$
(10)
The internal integral over $`\phi `$ is a standard Fresnel-Kirchhoff integral. Parameter $`b`$ in the exponent is the distance from the source on the defect to the exit point (Fig. 1)

$$b^2(x,\phi )=R^2+d^2\mathrm{cos}^2\left[\varphi (x)\right]2Rd\mathrm{cos}\left[\varphi (x)\right]\mathrm{cos}\left[\phi \varphi (x)\right].$$
(11)
It is convenient to express the distance $`b`$ directly through $`x`$ and $`\phi `$, which can be done using Eq. (9)

$$b^2(x,\phi )=R^2+(Rr_0)^22R\left(d\mathrm{cos}\phi +\sqrt{(Rr_0)^2d^2}\mathrm{sin}\phi \right)+2x\left(R\mathrm{sin}\phi \sqrt{(Rr_0)^2d^2}\right)+x^2.$$
(12)
Recall now, that the values of $`x`$ in the integral Eq. (10) are small $`xa^1r_0`$. It can also be seen from Fig. 1 that the outgoing ray is normal to the boundary when $`\mathrm{cos}\phi =d/(Rr_0)`$. This suggests that the difference

$$\delta =\phi \mathrm{cos}^1\left(\frac{d}{Rr_0}\right)$$
(13)
is a small parameter. In other words, the major contribution to the Fresnel-Kirchhoff integral comes from small $`\delta 1`$.

The integrand in Eq. (10) is a rapidly oscillating function. This allows to expand the phase of the oscillations

$$\mathrm{\Phi }(x,\phi )=M\varphi (x)+k_0\left[nb(x,\phi )R\mathrm{cos}(\psi \phi )\right]$$
(14)
in terms of $`x`$ and $`\delta `$

$$\mathrm{\Phi }(x,\phi )=A_xx+A_{xx}x^2+2A_{x\delta }x\delta +A_{\delta \delta }\delta ^2.$$
(15)
As it was already stated in the Introduction, the maximal directionality of the outgoing light is achieved for the position of the defect $`d=(Rr_0)/\sqrt{2}`$. To demonstrate this, we introduce a dimensionless deviation from the optimal defect position

$$\mathrm{\Delta }(d)=\frac{d}{R}\frac{Rr_0}{\sqrt{2}R}.$$
(16)
We will see that the width of the function $`I(\psi )`$ increases dramatically with $`\mathrm{\Delta }`$. Rather involved but straightforward calculations yield the following expressions for the coefficients in the expansion Eq. (15)

$`A_{x\delta }={\displaystyle \frac{nk_0R}{2r_0}}\left(1{\displaystyle \frac{4\mathrm{\Delta }^2}{14\mathrm{\Delta }^2}}\right){\displaystyle \frac{nk_0R}{2r_0}}\left(14\mathrm{\Delta }^2\right),`$ (17)
$`A_{\delta \delta }=k_0R\left[1{\displaystyle \frac{nR}{r_0}}\left(1{\displaystyle \frac{4\mathrm{\Delta }^2}{14\mathrm{\Delta }^2}}\right)\right]k_0R\left[1{\displaystyle \frac{nR}{r_0}}\left(14\mathrm{\Delta }^2\right)\right],`$ (18)
$`A_{xx}={\displaystyle \frac{nk_0}{4r_0}},A_x=a{\displaystyle \frac{\psi \pi /4}{\delta _\psi }},`$ (19)
where the parameter $`\delta _\psi `$ in the expression for $`A_x`$ is defined as

$$\delta _\psi =n\left(\frac{2r_0}{R}\right)^{1/2}\left[\frac{14\mathrm{\Delta }^2}{14n^2\mathrm{\Delta }^2}\right]^{1/2}n\left(\frac{2r_0}{R}\right)^{1/2}\left[1+2(n^21)\mathrm{\Delta }^2\right].$$
(20)
With the use of the expansion Eq. (15), the Fresnel-Kirchhoff integral can be easily evaluated yielding

$$I(\psi )\left|_0^{\mathrm{}}𝑑x\mathrm{exp}\left[x\left(aiA_x\right)+ix^2\left(A_{xx}\frac{A_{x\delta }^2}{A_{\delta \delta }}\right)\right]\right|^2.$$
(21)
The remaining integral over $`x`$ is of the Fresnel-type. However, it cannot be reduced to the special functions $`Ci(u)`$ and $`Si(u)`$, which describe the diffraction from a semi-infinite plane. This is because the linear term in the exponent contains a contribution $`ax`$ which is real. For this reason, it is convenient to introduce a new variable $`z=ax`$ in the integral (21). Upon substituting the coefficients (17)-(19) into Eq. (21) we arrive at the final result

$$I(\psi )\left|_0^{\mathrm{}}𝑑z\mathrm{exp}\left[z\left(1i\frac{\psi \pi /4}{\delta _\psi }\right)+iz^2\frac{n+1}{4n^2k_0r_0}F(\mathrm{\Delta })\right]\right|^2,$$
(22)
where the function $`F(\mathrm{\Delta })`$ is defined as

$$F(\mathrm{\Delta })=1+\frac{8nR\mathrm{\Delta }^2}{(n+1)r_0}.$$
(23)
As in Eqs.(17)-(19), we kept only the leading $`\mathrm{\Delta }^2`$ term in the definition of $`F`$. Now we can substantiate the statement that the optimal directionality of the emission is achieved at $`\mathrm{\Delta }=0`$. Indeed, $`z^2`$-term in the exponent of Eq. (22) leads to the broadening and oscillations of the angular dependence, $`I(\psi )`$. At small $`\mathrm{\Delta }`$ we have $`F1`$; then the $`z^2`$-term is multiplied by a small factor $`(k_0r_0)^11`$ and, thus, can be neglected. Then we immediately recover the lorentzian Eq. (3). On the other hand, for a general position of the defect we have $`\mathrm{\Delta }1`$, and $`FR/r_0`$. Then the $`z^2`$-term acquires a much larger coefficient $`2R/(nk_0r_0^2)`$, resulting in the loss of the directionality of the output light. This is illustrated in Fig. 2. It is seen that significant broadening and sideback oscillations set in already at small values of $`\mathrm{\Delta }`$. In particular, for $`\mathrm{\Delta }=0.3`$ the broadening is 60 percent.

## III Conclusion

Let us first discuss the validity of the assumptions used in the above calculation

(a) $`I(\psi )`$ was calculated within the scalar diffraction theory using Fresnel-Kirhhof approach. Note, that for a circular geometry, $`I(\psi )`$ can be calculated exactly by solving the scalar wave equation and treating defect as a perturbation. Then the expression for $`I(\psi )`$ is given by a sum over angular momenta of the leaking modes. Fresnel diffraction corresponds to replacing this sum by an integral. The accuracy of such a replacement is determined by the next term in the Poisson expansion, which contains an exponential factor $`\mathrm{exp}\left[2^{3/2}\pi nk_0(r_0R)^{1/2}\right]`$. Thus, the condition of validity of the Fresnel- Kirhgof approach is $`r_01/(k_0^2R)`$, which is not restrictive at all.

(b) According to Eq. (3), the full width at half maximum (FWHM) is equal to $`2\delta _\psi =2n(2r_0/R)^{1/2}`$. This equation was derived under the assumption that the defect is located far enough from the circumference of the disk, i.e. $`r_0\delta \rho R/(nk_0R)^{2/3}`$. It is possible to derive a more general expression for $`I(\psi )`$, that is valid for $`r_0\delta \rho `$, when the asymptotics Eq. (5) is not yet applicable. Derivation is based on the integral representation of the Bessel function and yields

$$I(\psi )\frac{2}{(\pi \gamma )^{1/2}}_0^{\mathrm{}}𝑑s\frac{e^{\gamma s^2}}{(1+s)^2+\left(\frac{\psi \pi /4\text{ }}{\text{ }\delta _\psi }\right)^2},$$
(24)
where the parameter $`\gamma `$ is defined as

$$\gamma =2^{1/2}nk_0R\left(\frac{r_0}{R}\right)^{3/2}.$$
(25)
It is seen that the condition $`r_0\delta \rho `$ corresponds to $`\gamma 1`$. Then we immediately recover the Lorentzian Eq. (3). At moderate $`\gamma `$, the FWHM is given by $`2C(\gamma )\delta _\psi `$, where the function $`C(\gamma )`$ is plotted in Fig. 2, inset. It is seen that within the wide interval $`1\gamma 10`$ the broadening factor $`C(\gamma )`$ changes very slowly. Then the FWHM can be expressed in terms of $`\gamma `$ as $`2^{4/3}nC(\gamma )\left(\gamma /nk_0R\right)^{1/3}`$, which is also a slow function of $`\gamma `$. Choosing for concreteness $`\gamma =1`$, we find for FWHM the expression $`3.35n/(nk_0R)^{1/3}`$.

We now turn to the numerical estimates. Three types of microdisk semiconductor lasers have been described in the literature so far. The lasers for wavelengths $`\lambda 1.5\mu m`$ have $`M`$-values reported are rather low ($`10M70`$) and $`n2.5`$. For this $`n`$ and maximal $`M=70`$ the FWHM is $`116^{}`$. Lasers for $`\lambda 0.8\mu m`$ have $`n3.1`$ and also rather small $`M`$ ($`30M300`$). With maximal $`M=300`$ we get $`89^{}`$ for FWHM. Nitride-based lasers operating at $`\lambda 0.4\mu m`$ have much higher $`M`$-values ($`200M600`$) and $`n2.8`$. This yields FWHM of $`64^{}`$. Microdisk lasers based on non-crystalline materials (polymer and dye solution) have also been reported. For this materials $`n1.8`$ is smaller and the values of $`M`$ ($`930`$ in  and 3000 in ) are high. Both factors tend to narrow $`I(\psi )`$. Namely, for $`M=1000`$ the FWHM of $`34^{}`$ can be achieved.

Let us discuss the physical meaning of the optimal condition, $`d=(Rr_0)/\sqrt{2}`$. As it is seen from Eq. (7), the phase of the whispering-gallery mode changes along the defect. As the defect plays a role of a source of the outgoing light, this change, $`\varphi (x)`$, is equivalent to the rotation of the line of the constant phase by an angle $`\mathrm{sin}^1[d/(Rr_0)]`$. Then, under the optimal condition, the line of the constant phase is perpendicular to the radial line drawn through the edge of the defect (Fig. 1). In other words, under the optimal condition, the defect can be replaced by a constant phase line at distance $`r_0`$ from the circumference that is parallel to the circumference. Clearly, the angular width of the far field emitted by this line is minimal for this parallel orientation.

Note in conclusion, that in  the improvement of the output characteristics of microdisk laser, achieved by introducing the deformation, is due to the fact that when the disk is deformed, the light rays are unable to stay within a whispering-gallery trajectory, and experience refractive escape in course of the chaotic motion. In the present paper we considered a perfectly circular microdisk with a defect. A point-like defect at some distance away from the boundary would be unable to couple out all the whispering-gallery modes, since it will not be able to affect the modes having a node at the defect position. Our main message here is that no whispering-gallery mode can evade an extended defect and will be directed out of the resonator as a result of scattering by this defect.


format: text

text:
# Filtering techniques for the detection of Sunyaev-Zel’dovich clusters in multifrequency CMB maps

## 1 INTRODUCTION

22footnotetext: Currently at Instituto de Física de Cantabria, Fac. de Ciencias, Av. de los Castros s/n, 39005-Santander, Spain

The detection and characterisation of the Sunyaev-Zel’dovich (SZ) effect is currently an area of considerable interest in millimetre and sub-millimetre astronomy. The SZ effect distorts primordial cosmic microwave background (CMB) radiation in the direction of galaxy clusters due to inverse Compton scattering of CMB photons by the hot intracluster plasma. In the context of the ESA Planck mission (and other CMB experiements with high sensitivity and high angular resolution), one must correct for this distortion (and for other foreground contaminants) in order to study the primordial anisotropies in the CMB. Nevertheless, the SZ effect itself is of considerable cosmological interest. The effect in individual clusters can be used to study the intracluster medium. Moreover, since the amplitude of the effect does not depend on redshift, it can also be used as a means for detecting clusters that would otherwise be unobservable. Perhaps most importantly for cosmology, a catalogue of the SZ effects in a large number of clusters, could be used to place constraints on cosmological parameters independently of those resulting from analysis of primordial CMB anisotropies. Indeed, a goal of the forthcoming Planck mission will be to produce a full-sky catalogue of SZ effects containing several tens of thousands of galaxy clusters. It is therefore crucial to have a robust and efficient method for detecting and extracting the SZ effect from multifrequency maps of the CMB.

The issue of component separation using multifrequency CMB maps has been thoroughly studied in the literature. Proposed methods include Wiener filtering (WF, Tegmark & Efstathiou 1996, Bouchet et al. 1997), the maximum-entropy method (MEM, Hobson et al. 1998, 1999), fast independent component analysis (FastICA, Maino et al. 2001), Mexican Hat Wavelet analysis (MHW, Cayón et al. 2000, Vielva el al. 2001a), matched filter analysis (MF, Tegmark & Oliveira-Costa 1998) and adaptive filter techniques (AF, Sanz et al. 2001, Herranz et al. 2001a, Herranz et al. 2001b). A non-parametric Bayesian approach to detecting SZ clusters has also recently been proposed by Diego et al. (2001b). A comparison between these methods is difficult because they have different specific purposes, use different sets of assumptions and the quality of the separation varies under different circumstances. For example, the MHW is designed to detect compact sources with a Gaussian profile (such as point sources convolved with a Gaussian beam) whereas other methods are better suited to deal with diffuse components such as dust and synchrotron emission. This suggests the possibility of combining several of these methods in order to improve the component separation, for example MEM $`+`$ MHW (Vielva et al. 2001b).

In general, a component separation method that uses all the available information will be more powerful than one that does not assume any prior knowledge about the data. The MEM, for example, produces excellent results when the power spectra and the frequency dependences of the components are well known. Unfortunately, if these assumptions about the data are incorrect, errors may arise that would affect the separation of several (or all) components. This is particularly dangerous in methods that perform the separation of all the components simultaneously (WF, MEM, FastICA). The opposite approach is to use a robust method that makes as few assumptions as possible about the data and aims to separate out just a single physical component. An example of the latter approach is the non-parametric SZ detection method given by Diego et al. (2001b), in which only the well-known frequency dependences of the SZ effect and the CMB are assumed. In general, however, such methods are less powerful than ones which assume a greater degree of (correct) prior knowledge about the data and the physical component of interest. Clearly, some compromise between robustness and effectiveness must be made when proposing a component separation method.

Filtering techniques (such as MHW, MF and AF) are single component separation methods that use some of the available information (i.e. the spatial structure of the component to be detected and the power spectrum of the combination of the other ‘background’ components), and generally reduce inaccuracies in the separation that may appear due to error propagation in methods that separate all components simultaneously. The MHW assumes a specific shape (a Gaussian) for the component to be separated (‘sources’) and, given the power spectrum of the background (that can be directly estimated from the data), finds the optimal scale of the filter. The MF generalises by allowing more general spatial profiles for the component (usually discrete objects) to be separated. Adaptive filters put additional emphasis on the characteristic scale of the sources in order to reduce further the number of spurious detections. Although one usually assumes spherical symmetry of the sources to be detected, it is not a general requirement of filtering method and filters can be easily generalised to detect asymmetric features.

In this paper, we discuss the generalisation of filtering techniques, in particular MF and AF, to the case of multiple data maps corresponding to different frequency channels. Multifrequency information can be used both to increase the signal of the sources and to reduce the contribution of the background (‘noise’). This information can be used prior to the filtering of the images (by using the correlations between the different channels to find an optimal combination of channels that maximises the signal-to-noise ratio of SZ clusters) or can be included directly in the construction of the filters. The structure of this paper is as follows. In section 2 we describe the formal aspects of the methods that make use of multifrequency information proposed in this paper. In section 3 we describe the simulations we made to test the multifilters. Section 4 is devoted to summarise and discuss our results.

## 2 Multifrequency filtering

There are two different approaches we can follow to include the frequency dependence of a signal in the filtering of multi-channel data. On the one hand, we can filter each channel separately, but carefully taking into account the cross-correlation between the different channels and the frequency dependence of the signal in order to obtain an output set of filtered maps that, when added to each other, is optimal for the detection. This philosophy gives birth to the *multifilter method*. On the other hand, one can use the information about correlations and frequency dependence before filtering in order to find the optimal combination of channels that maximizes the signal-to-noise ratio of the sources, and then use a filter on the optimally combined map. This leads to the design of a *combination method* plus a single filtering. These two methods are discussed in detail below. We first, however, define our model of the multifrequency observations.

### 2.1 Model of the observations

Let us consider a set of 2-dimensional images (maps) with data values defined by

$$y_\nu (\stackrel{}{x})=f_\nu s_\nu (x)+n_\nu (\stackrel{}{x}),x|\stackrel{}{x}|,\nu =1,\mathrm{},N,$$
(1)
where $`N`$ is the number of maps (or number of frequencies). The first term on the RHS represents the signal in each frequency map due to the thermal SZ effect in clusters, and the generalised noise $`n_\nu (\stackrel{}{x})`$ corresponds to the sum of the other emission components in the map.

For a map of any reasonable size, we would expect the presence of several SZ clusters. To illustrate more clearly the construction of multifilters, however, we assume that the signal is due to a single SZ cluster at the origin $`\stackrel{}{x}=0`$ of the map (the generalisation to several clusters distributed across the map is straightforward). In particular, we assume a spherically-symmetric $`\beta `$-profile for the cluster electron number density

$$n_e(r)[1+(x/r_c)^2]^{\frac{3}{2}\beta },$$

where $`r_c`$ is the core radius and we adopt the standard value $`\beta =2/3`$. One trivially obtains that the two-dimensional microwave decrement from such a cluster has the form $`A\tau (x)`$, where $`A`$ is the amplitude of the effect and $`\tau (x)`$ is the spatial template

$$\tau (x)=\frac{1}{[1+(\frac{x}{r_c})^2]^{1/2}}.$$
(2)
At each observing frequency this template is convolved with the corresponding antenna beam, which we assume to be a 2-D Gaussian of dispersion $`\theta _\nu `$, to produce the convolved profile $`\tau _\nu (x)`$. Thus, in (1), $`s_\nu (x)A\tau _\nu (x)`$. The quantity $`f_\nu `$ in (1) is the frequency dependence of the SZ effect, normalised such that $`f_{\nu _f}=1`$ at the fiducial frequency $`\nu _f`$. Hence $`A`$ is the true amplitude of the SZ effect at the frequency $`\nu _f`$.

The background $`n_\nu (\stackrel{}{x})`$ is modelled as a homogeneous and isotropic random field with average value $`n_\nu (\stackrel{}{x})=0`$ and cross-power spectrum $`P_{\nu _1\nu _2}(q)`$ ($`q|\stackrel{}{q}|`$) defined by

$$n_{\nu _1}(\stackrel{}{q})n_{\nu _2}^{}(\stackrel{}{q^{}})=P_{\nu _1\nu _2}(q)\delta _D^2(\stackrel{}{q}\stackrel{}{q^{}}),$$

where $`n_\nu (\stackrel{}{q})`$ is the Fourier transform of $`n_\nu (\stackrel{}{x})`$ and $`\delta _D^2`$ is the 2-D Dirac distribution.

### 2.2 The mutlifilter method

We now consider the first of two approaches to the filtering of mutlifrequency data, namely the multifilter. In fact, multifilters themselves can be constructed in several different ways. We consider here the two main possibilities, which are the scale-adaptive multifilter and the matched multifilter.

#### 2.2.1 Scale-adaptive multifilter (SAMF)

The idea of an optimal scale-adaptive filter was recently proposed by Sanz et al. (2001), which we now generalise to multifrequency data. One begins by introducing, for each observing frequency, a circularly-symmetric function $`\psi _\nu (x)`$. From this function, one can define a family of filter functions

$$\mathrm{\Psi }_\nu (\stackrel{}{x};R_\nu ,\stackrel{}{b})=\frac{1}{R^2}\psi _\nu \left(\frac{|\stackrel{}{x}\stackrel{}{b}|}{R_\nu }\right),$$
(3)
where the 3 parameters $`R_\nu `$ and $`\stackrel{}{b}`$ define a scaling and a translation respectively. For any particular values of these parameters, we define the filtered map at the $`\nu `$th observing frequency by

$$w_\nu (R_\nu ,\stackrel{}{b})=𝑑\stackrel{}{x}y_\nu (\stackrel{}{x})\mathrm{\Psi }_\nu (\stackrel{}{x};R_\nu ,\stackrel{}{b}).$$
(4)
and the ‘total’ filtered map as

$$w(R_1,\mathrm{},R_n,\stackrel{}{b})\underset{\nu }{}w_\nu (R_\nu ,\stackrel{}{b}).$$
(5)
The convolution (4) can be written as a product in Fourier space, in the form

$$w_\nu (R_\nu ,\stackrel{}{b})=𝑑\stackrel{}{q}e^{i\stackrel{}{q}\stackrel{}{b}}y_\nu (\stackrel{}{q})\psi _\nu (R_\nu q),$$
(6)
where $`y_\nu (\stackrel{}{q})`$ and $`\psi _\nu (q)`$ are the Fourier transforms of $`y_\nu (\stackrel{}{x})`$ and $`\psi _\nu (\stackrel{}{x})`$, respectively. A simple calculation then shows that the expectation value of the filtered field at the origin $`\stackrel{}{b}=0`$, is given by

$$w_\nu (R_\nu ,\stackrel{}{0})=2\pi Af_\nu 𝑑qq\tau _\nu (q)\psi _\nu (R_\nu q),$$
(7)
where the ensemble-average is over realisations of the background emission $`n_\nu (\stackrel{}{x})`$ and the limits in the integrals go from $`0`$ to $`\mathrm{}`$. Similarly, one finds that the variance of the total filter field (5) is given by

$$\sigma _w^2(R_1,\mathrm{},R_n)w^2(R_1,\mathrm{},R_n,\stackrel{}{b})w(R_1,\mathrm{},R_n,\stackrel{}{b})^2$$

$$=2\pi \underset{\nu _1,\nu _2}{}𝑑qqP_{\nu _1\nu _2}(q)\psi _{\nu _1}(R_{\nu _1}q)\psi _{\nu _2}(R_{\nu _2}q).$$
(8)
We choose the filter functions $`\psi _\nu (q)`$ to optimise the detection of the cluster at the origin, taking into account that the source has a bell shape with a single characteristic scale $`R_{o\nu }`$ in each map. We define the optimal scale-adaptive multifilter (SAMF) to be that which ensures the following conditions are satisfied:

1. $`w(R_{o1},\mathrm{},R_{on},\stackrel{}{0})`$ is an unbiased estimator of the amplitude of the source, so $`w(R_{o1},\mathrm{},R_{on},\stackrel{}{0})=A`$;
2. the variance of $`w(R_1,\mathrm{},R_n,\stackrel{}{b})`$ has a minimum at the scales $`R_{o1},\mathrm{},R_{on}`$, i.e. it is an $`\mathrm{𝑒𝑓𝑓𝑖𝑐𝑖𝑒𝑛𝑡}`$ estimator;
3. $`w_\nu (R_\nu ,\stackrel{}{b})`$ has a maximum at $`(R_{o\nu },\stackrel{}{0})`$.

The multifilter satisfying these conditions is given by the matrix equation

$$\stackrel{~}{\psi }(q)=P^1(\alpha F+G),$$
(9)
where we have introduced the column vectors $`\stackrel{~}{\psi }(q)=[\psi _\nu (R_\nu q)]`$, $`F=[f_\nu \tau _\nu ]`$, and $`G=[\mu _\nu \beta _\nu ]`$, where

$$\mu _\nu f_\nu \tau _\nu \left(2+\frac{d\mathrm{ln}\tau _\nu }{d\mathrm{ln}q}\right).$$

Also, $`P^1`$ is the inverse of the matrix $`P[P_{\nu _1\nu _2}(q)]`$, and the quantities $`\alpha `$ and $`\beta _\nu `$ are given by the components

$$\alpha =(A^1)_{00},\beta _\nu =(A^1)_{\nu 0},$$
(10)
where $`A`$ is the $`(1+n)\times (1+n)`$ matrix with elements

$$A_{00}𝑑\stackrel{}{q}F^tP^1F,A_{0\nu }𝑑\stackrel{}{q}\mu _\nu (F^tP^1)_\nu ,$$
(11)
$$A_{\nu 0}𝑑\stackrel{}{q}\mu _\nu (P^1F)_\nu ,A_{\nu \nu ^{}}𝑑\stackrel{}{q}\mu _\nu \mu _\nu ^{}(P^1)_{\nu \nu ^{}}.$$
(12)
The quantity $`\stackrel{~}{\psi }(q)`$ in (9) is called scale-adaptive multifilter, extending the concept considered in our previous paper (Sanz et al. 2001) for a single map. Also, from (8), we see that the variance can be written as

$`\sigma _w^2`$ $`=`$ $`{\displaystyle 𝑑\stackrel{}{q}\stackrel{~}{\psi }^tP\stackrel{~}{\psi }}`$
$`=`$ $`(A_{00})^1\left[1+2{\displaystyle 𝑑\stackrel{}{q}F^tP^1G}\right]+{\displaystyle 𝑑\stackrel{}{q}G^tP^1G}.`$
An interesting special case is when $`P`$ is a diagonal matrix, i.e. there is no cross-correlation between the backgrounds in the different frequency maps, so $`P_{\nu \nu ^{}}=\delta _{\nu \nu ^{}}P_\nu (q)`$. In this case, the multifilter is given by

$$\stackrel{~}{\psi }_\nu (q)=\frac{f_\nu \tau _\nu }{P_\nu }\frac{1}{a_\nu f_\nu ^2\frac{b_\nu ^2}{H_\nu }}\left[1\left(2+\frac{d\mathrm{ln}\tau _\nu }{d\mathrm{ln}q}\right)\frac{b_\nu }{H_\nu }\right],$$
(13)
where

$`a`$ $`=`$ $`{\displaystyle \underset{\nu }{}}f_\nu ^2{\displaystyle 𝑑\stackrel{}{q}\frac{\tau _\nu ^2}{P_\nu }},`$ (14)
$`b_\nu `$ $`=`$ $`f_\nu ^2{\displaystyle 𝑑\stackrel{}{q}\frac{\tau _\nu ^2}{P_\nu }\left(2+\frac{dln\tau _\nu }{dlnq}\right)},`$ (15)
$`H_\nu `$ $`=`$ $`f_\nu ^2{\displaystyle 𝑑\stackrel{}{q}\frac{\tau _\nu ^2}{P_\nu }\left(2+\frac{d\mathrm{ln}\tau _\nu }{d\mathrm{ln}q}\right)^2}.`$ (16)
If, additionally, one assumes that the backgrounds are white noise, i. e. the $`P_\nu `$ are constants and the resolution is the same in all maps, producing the convolved profile $`\tau (r_c,\theta )`$, then

$$f_\nu \stackrel{~}{\psi }_\nu (q)=\frac{\tau }{N}\frac{f_\nu ^2P_\nu ^1}{_\nu f_\nu ^2P_\nu ^1}\left[1\frac{b}{H}\left(2+\frac{d\mathrm{ln}\tau }{d\mathrm{ln}q}\right)\right],$$
(17)
where

$`N`$ $``$ $`{\displaystyle 𝑑\stackrel{}{q}\tau ^2}{\displaystyle \frac{b^2}{H}},`$ (18)
$`b`$ $`=`$ $`{\displaystyle 𝑑\stackrel{}{q}\tau ^2\left(2+\frac{d\mathrm{ln}\tau }{d\mathrm{ln}q}\right)},`$ (19)
$`H`$ $`=`$ $`{\displaystyle 𝑑\stackrel{}{q}\tau ^2\left(2+\frac{d\mathrm{ln}\tau }{d\mathrm{ln}q}\right)^2}.`$ (20)
#### 2.2.2 Matched multifilter (MMF)

If we do not assume condition (iii) in the definition of the SAMF in the previous subsection, one obtains another type of multifilter after minimization of the variance condition (ii) subject to the constraint (i). The multifilter satisfying these less restrictive conditions is given by the matrix equation

$$\stackrel{~}{\psi }(q)=\alpha P^1F,\alpha ^1=𝑑\stackrel{}{q}F^tP^1F,$$
(21)
where $`F`$ is the column vector $`F=[f_\nu \tau _\nu ]`$ and $`P^1`$ is the inverse matrix of the cross-spectrum $`P`$. This will be called matched multifilter extending the concept usually considered for a single map. In this case, the variance of the filtered field is given by

$$\sigma _w^2=𝑑\stackrel{}{q}\stackrel{~}{\psi }^tP\stackrel{~}{\psi }=\alpha .$$
(22)
In the special case where there is no cross-correlation between the backgrounds in the different maps, so $`P_{\nu \nu ^{}}=\delta _{\nu \nu ^{}}P_\nu (q)`$, the multifilter is given by

$$\stackrel{~}{\psi }_\nu (q)=\alpha f_\nu \tau _\nu P_\nu ^1,\alpha ^1=\underset{\nu }{}f_\nu ^2𝑑\stackrel{}{q}\frac{\tau _\nu ^2}{P_\nu }.$$
(23)
If additionally, one assumes that the backgrounds are white noise, i.e. the $`P_\nu `$ are constants, and the resolution is the same in all maps, producing the convolved profile $`\tau (r_c,\theta )`$, then

$$f_\nu \stackrel{~}{\psi }_\nu (q)=\frac{\tau }{N}\frac{f_\nu ^2P_\nu ^1}{_\nu f_\nu ^2P_\nu ^1},N𝑑\stackrel{}{q}\tau ^2.$$
(24)
A similar result has been independently developed by Naselsky et al. (2001).

### 2.3 The combination method

Our second approach to the filtering of mutlifrequency data is to perform a two stage process in which one first constructs some optimal linear combination of the individual frequency maps, and then applies a single filter to this combined map.

The linear combination of the original frequency maps is constructed so that it maximally ‘boosts’ objects with a given spatial template and frequency dependence above the background. In general, the combination map is given by

$$y(\stackrel{}{x})=\underset{\nu }{}c_\nu y_\nu (\stackrel{}{x}),$$
(25)
where $`c_\nu `$ is some set of weights (derived below). From (1), we see that this combined map can be written as

$$y(\stackrel{}{x})=At(x)+ϵ(\stackrel{}{x}),$$
(26)
where the ‘combined’ template $`t(x)`$ and background $`ϵ(\stackrel{}{x})`$ are given simply by

$$t(x)=\underset{\nu }{}c_\nu f_\nu \tau _\nu (x),ϵ(\stackrel{}{x})=\underset{\nu }{}c_\nu n_\nu (\stackrel{}{x}).$$
(27)
The optimal values of the weights $`c_\nu `$ are found by maximising the quotient $`Q=At(\stackrel{}{0})/\sigma _ϵ`$, where $`\sigma _ϵ`$ is the dispersion of the combined noise field $`ϵ(\stackrel{}{x})`$. Thus, one is maximising the ‘detection level’ of the object of interest above the background.

In practice, it is easier (but equivalent) to maximise $`Q^2`$. It is straightforward to show that, in the numerator of $`Q^2`$,

$$[t(\stackrel{}{0})]^2=c^tGc,$$
(28)
where

$$c=(c_\nu ),G=(G_{\nu \nu ^{}}),G_{\nu \nu ^{}}=f_\nu \tau _\nu (0)f_\nu ^{}\tau _\nu ^{}(0).$$
(29)
Similarly, in the denominator of $`Q^2`$,

$$\sigma _ϵ^2=c^tMc,$$
(30)
where

$$c=(c_\nu ),M=(M_{\nu \nu ^{}}),M_{\nu \nu ^{}}=n_\nu (\stackrel{}{x})n_\nu ^{}(\stackrel{}{x}).$$
(31)
Thus one must maximise, with repect to the weights $`c_\nu `$, the quantity

$$Q^2=A^2\frac{c^tGc}{c^tMc}.$$
(32)
The solution is easily found to be the eigenvector $`c`$ associated to the largest eigenvalue $`\lambda `$ of the generalized eigenvalue problem

$$(G\lambda M)c=0.$$
(33)
Once we have constructed the combined template $`t(x)`$, given in (25), with the optimal weights $`c`$ given by (33), we can apply a either a scale-adaptive filter or matched filter to the single combined image given by (26). The appropriate form of the scale-adaptive filter for a single image is given by Sanz et al. (2001). The appropriate form of the matched filter for a single image is given by

$$\stackrel{~}{\psi }(q)=\frac{1}{2\pi a}\frac{t(q)}{P(q)},$$
(34)
where

$$a=𝑑qq\frac{t^2(q)}{P(q)},P(q)=c^tNc,N=[P_{\nu \nu ^{}}(q)],$$
(35)
where $`t(q)`$ and $`P(q)`$ are the Fourier transform of $`t(x)`$ and the power spectrum of the noise $`ϵ`$, respectively.

### 2.4 Comparison of methods

We have introduced two different methods for filtering mutlifrequency data, namely multifilters and the combination method. Moreover, each method can be implemented for two different kinds of filter (scale-adaptive filters and matched filters). This makes four different possible ways to filter the data, and we would expect each method to have advantages and shortcomings with respect to the others.

From the methodological point of view, the main difference between the combination method plus a single filter and the multifilter method is that the first one uses the multi-frequency information to give the optimal starting point for the filter, while multifilters produce the optimal ending point after filtering. In that sense, multifilters are more powerful than a single filter applied to an optimally combined map. The cost of this higher efficiency, however, is an increase in the complexity of the filters and therefore of the computational time required to perform the data analysis. For example, in the application to simulated Planck observations in 10 frequency channels described below, the multifilter method requires $`10`$ times more CPU time than the combination method (because the last one filters only once).

Regarding the type of filter, matched filters give higher gains while adaptive filters reduce the number of errors in the detections (spurious sources). This is due to the fact that the optimal scale-adaptive filter includes a constraint (numbered (iii) above) about the scale of the sources that is not present in matched filter design. This constraint characterises with more precision the sources but also restricts the minimisation of the variance in the filtered maps (that is, the final gain).

## 3 SIMULATIONS

In order to test the multifrequency filtering methods outlined above, we now apply them to simulated Planck observations. These simulations mimic the main features of the future Planck mission, and are realistic in the sense that they include the latest available information about different physical components of emission (CMB, Galactic components, SZ effect and extra-galactic point sources) and that they reproduce the technical specifications of the different Planck channels (pixel sizes, antenna beams and noise levels). Table 1 shows the assumed observational characteristics of the simulated maps. The simulations were performed in patches of the sky of $`12.8^{}\times 12.8^{}`$. However, the method can be easily extended to the full sphere.

The CMB simulation was generated using the $`C_l`$’s provided by the CMBFAST code (Seljak & Zaldarriaga, 1996) for a spatially-flat $`\mathrm{\Lambda }`$CDM Universe with $`\mathrm{\Omega }_m=0.3`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ (Gaussian realization).

The Galactic emission is assumed to originate from four different components: thermal dust, spinning dust, free-free and synchrotron. Thermal dust was simulated using the template given by Finkbeiner et al. (1999). This model assumes that dust emission is due to two grey-bodies: a hot one with a dust temperature $`T_D^{hot}16.2K`$ and an emissivity $`\alpha ^{hot}2.70`$, and a cold one with $`T_D^{cold}9.4K`$ and $`\alpha ^{cold}1.67`$. These quantities are mean values.

For the free-free template we used one correlated with the dust emission in the manner proposed by Bouchet et al. (1996). The frequency dependence of the free-free emission is assumed to vary as $`I_\nu \nu ^{0.16}`$, and is normalised to give an rms temperature fluctuation of $`6.2\mu K`$ at 53 GHZ.

The synchrotron spatial template has been produced using the all sky map of Fosalba & Girardino 333ftp://astro.estec.esa.nl/pub/synchrotron, which is an extrapolation of the 408 MHz radio map of Haslam et al. (1982), from the original 1 deg resolution to a resolution of 5 arcmin. We further extrapolated the small-scale structure to 1.5 arcmin following a power-law power spectrum with an exponent of $`3`$. The frequency dependence is assumed to be $`I_\nu \nu ^{0.9}`$ and is normalized to the Haslam 408 MHz map. We include in our simulations information on the changes of the spectral index as a function of electron density in the Galaxy. This template has been made by combining the Haslam map with the Reich & Reich (1986) map at 1420 MHz and the Jonas et al. (1998) map at 2326 MHz, and can be found in the FTP site refered to above.

We have also taken into account the possible Galactic emission due to spinning grains of dust, proposed by Draine & Lazarian (1998). This component could be important at the lowest frequencies of Planck (30 and 44 GHz) in the outskirts of the Galactic plane.

The extragalactic point source simulations follow the same cosmological model as that assumed for the CMB simulation and corresponds to the model of Toffolati et al. (1998). The thermal SZ effect was made for the same cosmological model. The cluster population was modelled following the Press-Schechter formalism (Press & Schechter 1974) with a Poissonian distribution in $`\theta `$ and $`\varphi `$. The simulated cluster population fits well all the available X-ray and optical cluster data sets (see Diego et al. 2001a for a discussion). The different components used for the simulation are shown at 300 GHz in figure 1. Figure 2 shows the simulated channels taking into account all the components, as well as the antenna beam effect and the instrumental noise.

The simulation described above is the same used by Diego et al. (2001b). We have chosen this particular simulation in order to compare results with that work under the same conditions.

## 4 RESULTS AND DISCUSSION

### 4.1 Testing the methods

Before applying the methods presented in section 2 to the simulations described in section 3, let us illustrate how they work in a simplified case. For the sake of simplicity, the simulations include the same foregrounds and technical features as the simulations described in section 3, except for the fact that the spinning dust component was not included. Spinning dust is the weakest component and its contribution is almost negligible in most Planck channels. The other main difference with the realistic simulations is the thermal SZ effect itself. Instead of the realistic SZ clusters, we simulated 200 clusters of the same size ($`r_c=0.5`$ pixels), and amplitudes uniformly distributed between 0 and the maximum amplitude of the clusters belonging to the realistic simulation. The simulated clusters have the frequency dependence of the thermal SZ effect and the radial profile

$$\tau (x)=N\left(\frac{1}{\sqrt{r_c^2+x^2}}\frac{1}{\sqrt{r_v^2+x^2}}\right)$$
(36)
where $`N=r_vr_c/(r_v+r_c)`$. Here $`r_c`$ denotes the core radius and $`r_v`$ is a ‘cut scale’ that can be interpreted as the virial radius of the cluster. The profile given by (36) is a modified multiquadric profile that behaves as a $`\beta `$ model for $`x<<r_v`$ and decays quickly for $`x>>r_v`$. In order to simplify the simulations still further, the kinematic SZ effect is not included. In the real case, the contribution of the kinematic SZ effect is expected to be $`30`$ times lower than the contribution of the thermal SZ effect, thus justifying our approximation.

The simulated maps were filtered using the four possible combinations of filters and methods presented in this work. The combination method took about 1 minute of processing in a 700 MHz PC. The multifilter method took about 10 minutes of processing on the same computer. Processing times were slightly higher for the case of adaptive filters than in the case of matched filters. To compute the filters and the combinations of the maps the low resolution channels were re-binned to the resolution of the highest resolution channel (1.5 arcmin). Figure 3 shows the filters (in Fourier space) employed in the multifilter method. The adaptive filters are represented by solid lines, whereas matched filters are represented by dotted lines. Filters, specially those used in low frequency channels, are quite complicated and show several peaks at different wave numbers $`q`$ that correspond to the different scales the filters are trying to identify on the images. Conversely, the filters used after the combination method are rather simple and are shown in figure 4. The differences between matched and adaptive filters in the combination method are small. This indicates that there will be few differences in the results of both filters. In the case of the multifilters the differences are more pronounced, in particular for the 217 GHz channel.

The results from the test can be summarised as follows. The number of detections is higher in the multifilter case. This is not surprising since by definition the multifilter method is more powerful than the combination method. In particular, if one raises the detection threshold, the difference between the two methods increases. As an example, table 2 shows the results of one of the tests. After filtering, the sources were detected by looking for peaks above a $`4\sigma `$ threshold and then compared with a catalogue of the original simulated clusters. By comparison, in the same simulations, if the detection threshold is raised to $`5\sigma `$ (this case is not shown in table 2) the matched multifilters produce 90 detections whereas the matched filter in the combination method gives only 78.

A detected peak is considered a spurious detection if the distance to the closest object in the original catalogue is greater than 1.5 pixels. The number of spurious detections is strikingly low in all cases. Nevertheless, the adaptive filter seems to work better than the matched filter (for example, in table 2 the adaptive filter gives 0 spurious detections instead of 1 with the combination method). Of all the components present in the simulations, point sources are the most likely to produce spurious detections due to the similar scale of sources and SZ clusters. However, the frequency dependence of the SZ effect greatly reduces the probability of contamination.

The position of the sources was determined with errors below the pixel size (1.5 arcmin). However, the error in the determination of the amplitude is quite large ($`35\%`$). Table 2 shows that this error is due to a systematic bias. To explain this bias let us consider the case of matched filter. The filter normalization is given by integral $`a=_0^{\mathrm{}}𝑑qq\tau ^2/P`$, being $`\tau `$ the source profile in Fourier space and $`P`$ the power spectrum of the background. However, when we analyse a patch of the sky we do not have information about the power spectrum at all the wave numbers $`q`$. An image divided into pixels is limited by a minimum value $`q_{min}`$ and a maximum $`q_{max}`$. Therefore, the normalization we calculate is $`a^{}=_{q_{min}}^{q_{max}}𝑑qq\tau ^2/P`$. For the case of images with the same size and pixel scale as our simulations, and a multiquadric profile $`\tau `$ with $`r_c=0.5`$ pixels, the normalisations can differ by $`10\%`$ for a background spectral index $`\gamma =1`$ to $`60\%`$ for a background spectral index $`\gamma =3`$. The case for adaptive filters is more complicated but qualitatively similar. This bias is independent of the source flux and can be calibrated using simulations.

Figure 5 shows the contribution of each channel to the source amplitude estimation for each method. The channels with larger contribution are the 143 GHz and the 353 GHz. This is not surprising since they are the channels with more SZ contribution. The 100 GHz HFI channel has more contribution than the 100 GHz LFI because of its better signal-to-noise ratio. The combination method puts more emphasis in the 100 GHz and less in the 143 GHz channels than the multifilter method. Contributions from the 30 GHz, 217 GHz and 857 GHz channels are negligible. However, that does not mean that these channels do not contribute to the filter construction. For example, if we repeat the analysis with only the 5 channels with more contribution (70, 100 LFI, 100 HFI, 143 and 353 GHz), the number of detections reduces by $`10\%`$ and the number of spurious detections increases to 6 (at $`4\sigma `$ detection threshold) for the combination method case (for the two filters) and to 9 (adaptive filter) and 5 (matched filter) for the multifilter method. Even the 217 GHz channel is important; repeating the analysis with all channels except for the 217 GHz increases the error in the determination of the amplitudes and also increases the number of spurious detections at low $`\sigma `$ thresholds.

We can summarise the conclusions of the test as follows:

* Multifilter method is more powerful in the detection and estimation of cluster parameters.
* Combination method is faster than the multifilter method.
* Multifrequency information reduces the number of spurious detections. Therefore, it is not critical to use an adaptive filter to that end. The matched filter allows one to detect more sources.
* Some channels contribute more than others to the analysis, but all of them carry valuable information; the analysis should include all the available data.

### 4.2 Results for realistic simulations

Taking into account the insights provided by the test presented in the last sub-section, we are now prepared to confront the analysis of realistic simulations. These realistic simulations have been described in section 3. The main difference with respect to the previously performed test is that in the realistic case clusters have different sizes that are not known a priori. Herranz et al. (2001b) proposed a method to accommodate with this problem:

* Choose a trial core radius $`r_c`$ and construct the corresponding filters.
* Convolve the data with the filters, *varying the scales* $`R_\nu `$ *of eq.* (3). This is performed by substituting $`\stackrel{~}{\psi }(q)`$ by $`\stackrel{~}{\psi }(qx)`$, where $`x=R_\nu ^{}/R_\nu `$ is simply a dilation factor.
* In the case of scale-adaptive filters, condition (iii) (see section 2.2.1) implies that, if the assumed value of $`r_c`$ corresponds with the true core radius of the cluster, the coefficients will be maximum when $`x=1`$. If this is not true, the trial $`r_c`$ is discarded.
* In the case of matched filters, condition (iii) no longer holds. However, a similar criterion can be used: the best performance of the filter will occur when the scale of the filter and the scale of the cluster is the same. Therefore, if the signal-to-noise ratio of the coefficients is not maximum for $`x=1`$, the assumed radius $`r_c`$ is not the optimal choice. If this condition is verified for more than one value of $`r_c`$, the one with the highest signal-to-noise ratio is chosen.
* We repeat the process with as many different values of $`r_c`$ as desired.

In Herranz et al. (2001b), the above method has been successfully applied to single frequency maps containing simulated multiquadric profiles and different kind of backgrounds.

The clusters in the realistic simulations have a profile

$$\tau (r)=\frac{2}{\sqrt{1+(r/r_c)^2}}\mathrm{tan}^1\left[\sqrt{\frac{p^2(r/r_c)^2}{1+(r/r_c)^2}}\right],$$
(37)
where $`p=r_v/r_c`$ is the ratio between the virial radius and the core radius of the cluster. This profile is realistic from the point of view of the physics of the Sunyaev-Zel’dovich effect. The profiles (36) and (37) are almost identical for $`r<r_v`$. In order to calculate the filters in Fourier space, it is easier to work with the modified multiquadric (36). Therefore, as a first approximation we will assume that the clusters can be described by the profile (36). As we will see, this is a good approximation.

Following the results of the test in the last subsection, we choose a matched multifilter to perform the analysis. After applying our method to the simulations we detect the clusters by looking for sets of connected pixels above a certain threshold. The maximum of these sets give the position and the amplitude of the sources. At the $`3\sigma `$ level (regions with 5 or more connected pixels) we are able to detect 63 cluster candidates, of which 62 correspond to real clusters. The spurious detection appears in one of the borders of the image and therefore can be considered as an edge effect. The mean error in the determination of the position of the clusters was $`1`$ pixel.

Using the multi-scale analysis we are able to determine the core radii of the clusters with a mean absolute error of 0.30 pixels. The mean bias in the determination of the core radii was $`0.15`$ pixels. Since pixelisation effects are expected to corrupt structures whose typical scale is much smaller than the pixel size, all clusters with correspondingly small core radii will not follow the multi-quadric profile. Most of the clusters have very small core radii and therefore can be considered as point-like sources. To detect these clusters, a Gaussian profile (corresponding to the beam at each channel) was assumed instead of the multiquadric given by (36). The separation between clusters that are considered as point sources and extended sources was set in $`r_c=0.4`$ pixel, since if $`r_c`$ is below this limit the FWHM of the multiquadric profile is of the order of a pixel. All the clusters that were detected as point sources were considered to have the (arbitrary) value $`r_c=0.1`$ pixel. The bias in the determination of the core radius mentioned above is due to the asymmetry introduced when assigning $`r_c=0.1`$ pixel to all clusters whose core radius is less than 0.4 pixel. Figure 6 shows the resulting maps after filtering with different scales (point source in the left, $`r_c=1.0`$ pixel in the middle and $`r_c=2.5`$ pixel in the right side). Small clusters stand out mainly in the left panel, while large clusters are emphasized in the right panel.

Figure 7 shows the detected amplitudes of the 60 brightest detected clusters vs. the true amplitudes. The brightest sources in the graph were recovered with relative errors of $`30\%`$. As might be expected, the error increases as the true amplitude decreases. For the weakest clusters, the detected amplitude tends to a horizontal line. This bias is due to the incorrect determination of the normalization described in section 4.1, as well as two additional factors: first, only the faint clusters that are overamplified are able to reach the detection threshold; second, the detection is performed by looking for peaks in the data and, hence, those clusters which are by chance enhanced by residual fluctuations of the background (noise) are more likely to be detected. The bias due to the normalization should be the same for all sources with a given scale and therefore could be calibrated from simulations, whereas the second kind of bias affects mainly the weak clusters (that dominate the number counts in this case) and can not be easily calibrated. However, the amplitude estimation could be improved by refining the detection method.

The catalogue of detected sources was complete above 0.17 Jy (at 300 GHz). Below this limit the completeness gradually decreases. There are a few very faint clusters that are detected with this method (fluxes of few tens of mJy), but the practical detection limit is $`4050`$ mJy (at 300 GHz). The 62 detections in our small sky patch of $`(12.8^{})^2`$ corresponds to $`10000`$ detections in $`2/3`$ of the sky (the area not dominated by the Galaxy).

The performance of the other filtering techniques proposed in this work was compared to the performance of the matched multifilter for the realistic Planck simulations. For the comparison, only the point-like clusters (that dominate the number counts) were considered. The scale-adaptive multifilter, as seen in section 4.1, performs worse than the matched multifilter because their main advantage, that is the removal of spurious sources, has already been accomplished by the multi-frequency information included in the filters. Since the clusters one is trying to detect are very faint, it is the overall gain which is most important. At $`3\sigma `$, the scale-adaptive multifilter detects 40 per cent fewer clusters than the matched multifilter. Both matched and scale-adaptive filtering of an optimally combined map are slightly less sensitive than multifilters, producing, in both cases, around 20 per cent fewer detections at $`3\sigma `$ than with the matched multifilter. In other words, in the case of single filters, matched and scale-adaptive perform very similarly (as was expected from figure 4, where both filters look alike).

The results obtained with the matched multifilter are comparable in number of detections and reliability at the $`3\sigma `$-threshold to those obtained by Diego et al. (2001b). To perform their method it is necessary first to clean the maps carefully. This cleaning includes point source removal with a MHW technique (Vielva el al. 2001a), dust substraction using the 857 GHz map and CMB substraction in Fourier space up to a certain wavenumber limit using the 217 GHz map. In order to compare results, we performed the same cleaning method and then applied the matched multifilter. At $`3\sigma `$ detection threshold, we obtained a $`7`$ per cent more detections than before. At $`2\sigma `$, we also detected $`7`$ per cent more clusters and the number of spurious detections decreased by about $`12`$ per cent. This improvement is small given all the cleaning work that is necessary to obtain it. Moreover, it indicates that the multifilter scheme is powerful and robust enough to deal successfully with contaminants such as point sources, dust and the CMB signal.

## 5 Conclusions

The filtering techniques described in this work have proven to be useful tools for the detection of discrete objects in multifrequency data, which are known to have a prescribed frequency dependence. The only other assumption about the objects is their spatial template; any other quantities of interest are directly calculated from the data. We have assumed in this work that the sources have circular symmetry, but the method can be extended to deal with asymmetric structures. Future work will take this point into account. We have applied these techniques to the problem posed by the presence of Sunyaev-Zel’dovich effect in multifrequency CMB maps. The best results are obtained with matched multifilters. The detection level of the method is comparable to other single-component separation methods, but the number of assumptions and previous data processing is reduced to a minimum. Owing to the weakness of the signals that we aim to detect, the determination of the amplitude of the clusters is biased; this is the strongest limitation of the method. Therefore, the method should be complemented with a posteriori analysis (for example data fitting) in order to improve the amplitude estimation of weak clusters. Finally, the filters here presented can be used in other fields of Astronomy, such as X-ray observations, and large-scale structure detection.

## Acknowledgments

We thank Patricio Vielva for useful suggestions and comments. DH acknowledges support from a Spanish MEC FPU fellowship and the hospitality of the Astrophysics Group of the Cavendish Laboratory during a research stay in July 2001. RBB acknowledges financial support from the PPARC in the form of a research grant. RBB thanks the Instituto de Física de Cantabria for its hospitality during a stay in 2001. JMD acknowledges support from a Marie Curie Fellowship of the European Community programme Improving the Human Research Potential and Socio-Economic knowledge under contract number HPMF-CT-2000-00967. We thank FEDER Project 1FD97-1769-c04-01, Spanish DGESIC Project PB98-0531-c02-01 and INTAS Project INTAS-OPEN-97-1192 for partial financial support.


text:
# Detection of 𝑩_𝒔→𝝁⁺𝝁⁻ at the Tevatron Run II and Constraints on the SUSY Parameter Space

## Abstract

A measurement of the branching ratio for the rare decay mode $`B_s\mu ^+\mu ^{}`$ at the Tevatron is an opportunity to test various supersymmetric scenarios. We investigate the prospects for studying this mode in Run II and estimate that CDF would be sensitive to this decay for a branching ratio $`>1.2\times 10^8`$ with 15 $`\mathrm{fb}^1`$(or, if a similar analysis holds for D0, $`>6.5\times 10^9`$ for the combined data). We calculate the branching ratio in minimal supergravity (mSUGRA) parameter space, and find that $`\mathrm{tan}\beta >30`$ can be probed. (This mSUGRA parameter space cannot be probed by direct production of SUSY particles at Run II.) Including other experimental constraints on the mSUGRA parameter space, one finds that CDF $`B_s\mu ^+\mu ^{}`$ measurements would be able to cover the full mSUGRA parameter space for $`\mathrm{tan}\beta =50`$ if the muon $`g_\mu 2`$ anomaly exceeds $`11\times 10^{10}`$, and about half the allowed parameter space for $`\mathrm{tan}\beta `$ = 40. A large branching ratio $`>7(14)\times 10^8`$ (feasible with only 2 $`\mathrm{fb}^1`$) would be sufficient to exclude the mSUGRA model for $`\mathrm{tan}\beta 50(55)`$. Dark matter neutralino-proton detection cross sections are examined in the allowed region, and should be large enough to be accessible to future planned experiments. Combined measurements of $`B_s\mu ^+\mu ^{}`$, the Higgs mass $`m_h`$ and the muon $`g_\mu 2`$ anomaly would be sufficient to determine the $`\mu >0`$ mSUGRA parameters (or show the model is inconsistent with the data). We also briefly discuss the $`B_s\mu ^+\mu ^{}`$ decay in R parity violating models. There, for some models, the branching ratio can be large enough to be detected even for small $`\mathrm{tan}\beta `$ and large $`m_{1/2}`$.

preprint: CTP-TAMU-02-02

1. INTRODUCTION. Supersymmetric (SUSY) extensions of the Standard Model (SM) represent a natural candidate for the new physics expected to occur in the TeV energy domain. One of the difficulties in determining predictions of such models lies in the large number of new parameters the theory implies. Thus the most general low energy model, the Minimal Supersymmetric Standard Model (MSSM), has over 100 free parameters, and even the minimal model based on supergravity grand unification at the scale $`M_G=2\times 10^{16}\mathrm{GeV}`$ (mSUGRA) possesses four new parameters and one sign in addition to the SM parameters. Fortunately SUSY models apply to a large number of different accelerator and cosmological phenomena, and a great deal of effort has been involved in recent years to to use this data to limit the parameter space. Part of the difficulty in doing this resides in the success of the model in not disturbing the excellent agreement of the precision tests of the SM  due to the SUSY decoupling theorems which suppress SUSY contributions at low energies. Historically, the absence of flavor changing neutral currents at the tree level played an important role in the construction of the SM. They represent therefore an important class of phenomena that might show the presence of new physics, since the SM and the SUSY contributions contribute first at the loop level with comparable size. Thus the decay $`bs\gamma `$ has been a powerful tool in limiting the SUSY parameter space. In this paper, we consider the decay $`B_s\mu ^+\mu ^{}`$ within the framework of mSUGRA models and R parity violating models. This process is particularly interesting for several reasons: The SM branching ratio is quite small, $`i.e.`$, $`Br[B_s\mu ^+\mu ^{}]_{\mathrm{SM}}=3.5\times 10^9`$ . The SUSY contribution  has terms that grow as $`\mathrm{tan}^6\beta `$ (where $`\mathrm{tan}\beta =<H_2>/<H_1>`$ and $`H_{1,2}`$ are the two SUSY Higgs bosons) and thus can become quite large for large $`\mathrm{tan}\beta `$. Finally, as we shall show below, the two collider detectors at the Tevatron will be sensitive to this decay for $`\mathrm{tan}\beta \begin{array}{c}>\hfill \\ \hfill \end{array}30`$ in Run II.

2. MSUGRA MODELS. We consider first mSUGRA models with R parity invariance, and combine the existing accelerator constraints on the parameter space ($`e.g.`$, the light Higgs mass $`m_h`$ bound, $`bs\gamma `$$`etc.`$), the cosmological dark matter constraints, the possible muon magnetic moment anomaly, with what might be expected from the Tevatron measurement of the $`B_s\mu ^+\mu ^{}`$ decay. The combined constraints can significantly limit the SUSY parameter space, and thus allow better predictions as to what the models predict at the LHC.

The mSUGRA model  depends on four parameters which we take to be the following: $`m_0`$ (the universal scalar mass at $`M_G`$), $`m_{1/2}`$ (the universal gaugino mass at $`M_G`$), $`A_0`$ (the universal cubic soft breaking mass at $`M_G`$), and $`\mathrm{tan}\beta `$ (at the electroweak scale). In addition, the sign of $`\mu `$ (the Higgs mixing parameter appearing in the superpotential as $`\mu H_1H_2`$) is arbitrary. We consider the parameter range of $`m_0,m_{1/2}<1\mathrm{TeV}`$, $`3<\mathrm{tan}\beta <55`$, and $`|A_0|<4m_{1/2}`$. We take a $`2\sigma `$ bound on the $`bs\gamma `$ decay  of $`1.8\times 10^4<Br[bs\gamma ]<4.5\times 10^4`$, and the LEP bound on the light Higgs of $`m_h>113.5\mathrm{GeV}`$. Since there is still a (2-3) GeV uncertainty in the theoretical calculation of $`m_h`$, we will (conservatively) interpret this to mean $`m_h(\mathrm{theory})>111\mathrm{GeV}`$. For $`\mathrm{tan}\beta \begin{array}{c}>\hfill \\ \hfill \end{array}45`$, results are sensitive to the precise values of $`m_b`$ and $`m_t`$. We assume here that $`m_b(m_b)`$ = 4.25 GeV and $`m_t(pole)`$ = 175 GeV, and include the large $`\mathrm{tan}\beta `$ corrections to the $`b`$-quark Yukawa coupling constant . We assume that any muon $`g2`$ deviation from the SM is due to SUSY . The experimental deviation has now been reduced to a $`1.6\sigma `$ effect , and we take here this deviation to be greater than $`1\sigma `$ below the current central value $`i.e.`$, $`a_\mu ^{\mathrm{SUGRA}}>11\times 10^{10}`$. We will show however, what would happen if this deviation were to change with the new BNL E821 data currently being analyzed. For models with R parity invariance, the lightest neutralino, $`\stackrel{~}{\chi }_1^0`$ is the dark matter particle, and we require that $`0.07<\mathrm{\Omega }_{\stackrel{~}{\chi }_1^0}h^2<0.21`$, in accord with current CMB and other astronomical data . All stau-neutralino co-annihilation channels  are included in the relic density calculations.

3. $`B_s\mathbf{}\mu ^\mathbf{+}\mu ^{\mathbf{}}`$ DECAY. The branching ratio for $`B_s\mu ^+\mu ^{}`$ is given in  which we write in the form

$`Br[B_s\mu ^+\mu ^{}]`$ $`=`$ $`{\displaystyle \frac{2\tau _BM_B^5}{64\pi }}f_{B_s}^2\sqrt{1{\displaystyle \frac{4m_l^2}{M_B^2}}}`$ (2)
$`\left[\left(1{\displaystyle \frac{4m_l^2}{M_B^2}}\right)\left|{\displaystyle \frac{(C_SC_S^{})}{(m_b+m_s)}}\right|^2+\left|{\displaystyle \frac{(C_PC_P^{})}{(m_b+m_s)}}+2{\displaystyle \frac{m_\mu }{M_{B_s}^2}}(C_AC_A^{})\right|^2\right]`$
where $`f_{B_s}`$ is the $`B_s`$ decay constant, $`M_B`$ is the $`B`$ meson mass, $`\tau _B`$ is the mean life and $`m_l`$ is the mass of lepton. $`C_S`$, $`C_S^{}`$, $`C_P`$, $`C_P^{}`$ include the SUSY loop contributions due to diagrams involving the particles such as stop, chargino, sneutrino, Higgs etc.. For large $`\mathrm{tan}\beta `$, the amplitude has terms that grow like $`\mathrm{tan}^3\beta `$ as can be seen in the example of Figure 1. Thus at large $`\mathrm{tan}\beta `$, the dominant contribution to $`C_S`$ is given approximately by

$`C_S`$ $``$ $`{\displaystyle \frac{G_F\alpha }{\sqrt{2}\pi }}V_{tb}V_{ts}^{}\left({\displaystyle \frac{\mathrm{tan}^3\beta }{4\mathrm{sin}^2\theta _W}}\right)\left({\displaystyle \frac{m_bm_\mu m_t\mu }{M_W^2M_A^2}}\right){\displaystyle \frac{\mathrm{sin}2\theta _{\stackrel{~}{t}}}{2}}\left({\displaystyle \frac{m_{\stackrel{~}{t}_1}^2\mathrm{log}\left[\frac{m_{\stackrel{~}{t}_1}^2}{\mu ^2}\right]}{\mu ^2m_{\stackrel{~}{t}_1}^2}}{\displaystyle \frac{m_{\stackrel{~}{t}_2}^2\mathrm{log}\left[\frac{m_{\stackrel{~}{t}_2}^2}{\mu ^2}\right]}{\mu ^2m_{\stackrel{~}{t}_2}^2}}\right)`$ (3)
where $`m_{\stackrel{~}{t}_{1,2}}`$ are the two stop masses, and $`\theta _{\stackrel{~}{t}}`$ is the rotation angle to diagonalize the stop mass matrix. We need to muliply the above expression by $`1/(1+ϵ_b)^2`$ to include the SUSY QCD corrections. $`ϵ_b`$ is proportional to $`\mu \mathrm{tan}\beta `$ . We have $`C_P=C_S,C_S^{}=(m_s/m_b)C_s`$ and $`C_P^{}=(m_s/m_b)C_P`$. The operators are given in ref.. In our numerical calculation, we use all complete one loop contributions to the branching ratio.

4. DETECTION OF $`B_s\mathbf{}\mu ^\mathbf{+}\mu ^{\mathbf{}}`$ at TEVATRON. We consider now the possibility of detecting the decay $`B_s\mu ^+\mu ^{}`$ by the CDF and D0 detectors at the Tevatron in Run II. Both detectors have been upgraded with excellent tracking and muon detector systems . The dimuon trigger is the key to collect the $`B_s\mu ^+\mu ^{}`$ decays.

In order to estimate the limits on $`Br[B_s\mu ^+\mu ^{}]`$ detection, we use the 95$`\%`$ C.L. limit on $`Br[B_s\mu ^+\mu ^{}]`$ published by CDF. Thus our discussion is based on the CDF detector, although both CDF and D0 detectors should have a similar perfromance.

In the Run I analysis, CDF observed one candidate that was consistent with $`B_s\mu ^+\mu ^{}`$ with an estimate of 0.9 backound (BG) events in 98 $`\mathrm{pb}^1`$ . The primary Run-I selection variables and cut values were that at least one muon track with $`c\tau L_{xy}M_B/P_T^{\mu \mu }>100\mu `$m, $`IP_T^{\mu \mu }/[P_T^{\mu \mu }+\mathrm{\Sigma }P_T]>0.75`$ for the muon pair, and $`\mathrm{\Delta }\mathrm{\Phi }<0.1`$ rad. Here, $`L_{xy}`$ is the transverse decay length; $`p_T^{\mu \mu }`$ is the transverse momentum of the dimuon system. $`\mathrm{\Sigma }P_T`$ is the scalar sum of the transverse momenta of all tracks, excluding the muon candidates, within a cone of $`\mathrm{\Delta }R\sqrt{(\mathrm{\Delta }\eta )^2+(\mathrm{\Delta }\varphi )^2}=1`$ around the monetum vector of the muon pair. The $`z`$ coordinate of each track along the beam line  must be within 5 cm of the primary vertex. $`\mathrm{\Delta }\mathrm{\Phi }`$ is an opening azimuthal angle between $`P_T^{\mu \mu }`$ and the vector pointing from the primary vertex to the secondary vertex (the reconstructed $`B`$-meson decay position). As a conservative estimate, CDF took the one event as signal to calculate 95% C.L. limit of signal events ($`N_1^{95\%}`$5.06 events) and had set a limit of $`Br<2.6\times 10^6`$. In the analysis, the selection efficiency ($`ϵ`$) for signal events and the rejection power ($``$) for background events (pass a baseline selection ) are estimated to be $`ϵ_1=0.45`$ and $`_1=440`$ by using a sample of like-sign dimuon events ($`5<M_{\mu \mu }<6\mathrm{GeV}/c^2`$).

A dimuon trigger in Ref. will improve the acceptance for signal events by a factor of 2.8. The trigger will soon be tested using the Run IIa (2 $`\mathrm{fb}^1`$) data. This will allow us to modify the trigger design for the higher luminosity expected in Run IIb (15 $`\mathrm{fb}^1`$). In our analysis, we assume that the dimuon trigger can be designed by maintaing the acceptance for signal events. We expect to improve the acceptance for signal events by a factor of 2.8 . If we assume the factor 2.8 to be the same for BG events, then we would observe 51 (386) events in 2 (15) $`\mathrm{fb}^1`$ with the same cuts as in Run I. Therefore, CDF has to require a set of tighter cuts to obtain the best possible upper limit.

We first need to understand the background contents and expected improvement by the new Run II detector. Two types of backgrounds must be taken into account: (i) non-$`b`$ backgrounds comming from the primary vertex; (ii) $`b`$ background events, such as the gluon-spliting $`b\overline{b}`$ events.

The most important feature for $`B`$ decays is the displaced vertex. One way to reduce prompt background is to require a minimum decay length $`L_{xy}`$. However, two tracks can appear to form a secondary vertex if one of two tracks originates from the primary vertex and the other has an impact parameter ($`d`$). Therefore, the requirement of a minimal impact parameter of individual tracks can further clean up the sample. It has been shown for example, in the Run-I analysis for $`B^0K^0\mu ^+\mu ^{}`$ events , that a tight impact parameter cut on significance for individual track ($`d/\sigma _d>2`$) significantly improve the background rejection even with $`L_{xy}>100\mu `$m. One has then $`ϵϵ_1\times 0.43`$ and $`_1\times 190`$.

We indeed observed the similar numbers of OS and LS events in the Run-I final selection for $`B_s\mu ^+\mu ^{}`$ events, while we should expect more OS events than LS events if the cuts were tight enough for non-$`b`$ background events. This indicates the final dimuon event candidates in Run I analysis are actually dominated by the non-$`b`$ background sources. Thus a higher track impact parameter is neccessary to reduce the non-$`b`$ backgrounds. We would expect larger reduction with good efficiency even after the $`L_{xy}`$ cut. In other words, we would have been able to set better limits if we had included the impact parameter for the cut optimization. The silicon vertex detector (SVX-II) will provide us much better reduction for the non-$`b`$ background than Run I. The non-$`b`$ background will not be a problem in Run II.

In Run II, the most severe background will be the two muons from gluon-spliting $`b\overline{b}`$ events. Since both particles are $`b`$ quarks, the impact parameter does not help. Both $`b`$ and $`\overline{b}`$ also go in the same direction, so that cut on $`L_{xy}`$ does not help either. However, $`\mathrm{\Delta }\mathrm{\Phi }`$ is still usefull to remove the background events. Furthermore, in Run II, we can use $`\mathrm{\Delta }\mathrm{\Theta }`$ in $`r`$-$`z`$ view since we have $`z`$-strips in SVX-II.

There is some room to improve the isolation cut. We can form a new isolation parameter by only using the tracks with large impact parameter. This new isolation cut will work to reject the $`b\overline{b}`$ rather than non-$`b`$ background. Furthermore, we can search for tracks with large impact parameter on the opposite side of the dimuon candidates to make sure that the $`b`$ and $`\overline{b}`$ go to the opposite side.

Therefore, CDF could improve the BG rejection by a factor of 200-400 with further reduction of the signal efficiency by a factor of 2-3. Based on these facts, we now consider two cases to evaluate Run II limits as a function of luminosity.

In the first case (Case A), we naively assume new tighter cuts in Run II, described above, will gain additional BG rejection power of 450 for additonal efficiency of 0.45, or $`_2=450^{0.45/ϵ_2}`$. This gives us

$`{\displaystyle \frac{ϵ_2}{ϵ_1}}`$ $`=`$ $`{\displaystyle \frac{1}{1+\mathrm{log}(_2/_1)/\mathrm{log}(450)}}`$ (4)
If we could optimize the BG rejection in Run IIa (2 $`\mathrm{fb}^1`$) to be $`_251_1`$ with $`ϵ_20.61ϵ_1`$ (from Eq. 4), then we would expect one BG event in 2 $`\mathrm{fb}^1`$. Thus, with an assumption of the same size of the total systematic uncertainty in Run II as in Run I, we can extrapolate the 95% C.L. limit to be $`Br<7.7\times 10^8`$ for 2 fb<sup>-1</sup> using $`N_1^{95\%}`$.

In the second case (Case B), we simply assume the Run-II background rejection could be improved (without loosing the signal efficiency) to keep the expected BG events in 2 $`\mathrm{fb}^1`$ at the level of Run I ($`i.e.`$, 0.9 events). If we would observe one event in 2 $`\mathrm{fb}^1`$, then we could set the limits by scaling the Run-I $`Br`$ limit down by the luminosity (2000 $`\mathrm{pb}^1`$/98 $`\mathrm{pb}^1`$) and the acceptance by (2.8/1.0). Thus we obtain $`Br<4.6\times 10^8`$. This would certainly be the optimistic scenario, but it would be a goal of this analysis in Run IIa. Here, the systematic uncertainty in Run II is assummed to be the same as in Run I.

We repeat the same arguement for different luminosities. Figure 2 shows 95% C.L. limits on $`Br[B_s\mu ^+\mu ^{}]`$ at CDF in Run II as a function of integrated luminosity for Cases A and B. For 15 $`\mathrm{fb}^1`$ in case A, CDF would be sensitive to $`Br>1.2\times 10^8`$ and the combined CDF and D0 data (30 $`\mathrm{fb}^1`$) would be sensitive to $`Br>6.5\times 10^9`$.

5. RESULTS FOR MSUGRA. We examine first the parameter region for the mSUGRA model that would be accessible to CDF at Run II with 15 $`\mathrm{fb}^1`$ of data. Figure 3 shows the $`Br[B_s\mu ^+\mu ^{}]`$ as a function of $`m_{1/2}`$ for $`A_0`$ = 0, $`m_0`$ = 300 GeV. One sees that with a sensitivity of $`Br[B_s\mu ^+\mu ^{}]>1.2\times 10^8`$ for 15 $`\mathrm{fb}^1`$, the Tevatron Run II can probe the $`B_s\mu ^+\mu ^{}`$ decay for $`\mathrm{tan}\beta >30`$. Further, a search for this decay would sample much higher regions of $`m_{1/2}`$ than a direct search at Run II for SUSY particles which is restricted to $`m_{1/2}<250\mathrm{GeV}`$ . As $`m_0`$ increases, the branching ratio goes down. However, this dependence becomes less significant for large $`m_{1/2}`$, where $`m_0`$ as large as 800 GeV can be sampled for large $`m_{1/2}`$.

In Figure 4 the contours of $`Br[B_s\mu ^+\mu ^{}]`$ are plotted in the $`m_0`$-$`m_{1/2}`$ plane for $`\mathrm{tan}\beta `$ = 50, $`A_0`$ = 0. We combine now this result with the other experimental constraints. Thus the shaded region to the left is ruled out by the $`bs\gamma `$ constraint, and the shaded region on the right hand side is disallowed if $`a_\mu ^{\mathrm{SUGRA}}>11\times 10^{10}`$. The narrow shaded band in the middle is allowed by the dark matter constraint. We note that independent of whether the astronomically observed dark matter is SUSY in origin, the dark matter allowed region for mSUGRA cannot significantly deviate from this shaded region, for below the narrow shaded band, the stau would be lighter than the neutralino (leading to charged dark matter), while above the band, mSUGRA would predict more neutralino dark matter than is observed.

Using our estimate that $`Br>1.2\times 10^8`$ can be observed with 15 $`\mathrm{fb}^1`$, we see that almost the entire parameter space allowed by the $`a_\mu `$ constraint can be probed in Run II for $`\mathrm{tan}\beta `$ = 50. Note that an observed $`Br[B_s\mu ^+\mu ^{}]>7\times 10^8`$, possible with only 2 $`\mathrm{fb}^1`$(see Figure 2), would be sufficient to rule out the mSUGRA model for $`\mathrm{tan}\beta 50`$. In Figure 4 we also show the expected dark matter detector cross section for Milky Way dark matter (the short solid lines). They are of a size that can be observed by future planned dark matter detectors such as GENIUS, Cryoarray, ZEPLIN IV and CUORE. In Figure 5 we plot the same information for $`\mathrm{tan}\beta `$ = 40, $`A_0`$ = 0. We see here about half the parameter space can be scanned by the CDF detector (and the whole parameter space if a similar analysis holds for the D0 detector). See Figure 2. We note, further, that if $`A_0`$ = 0, a simultaneous measurement of both $`Br[B_s\mu ^+\mu ^{}]`$ and $`a_\mu `$ would essentially determine the mSUGRA parameters, as the $`m_0`$ allowed region at fixed $`m_{1/2}`$ is very narrow due to the dark matter constraint. The effect of varying $`A_0`$ is shown in Figure 6, where the allowed region for $`A_0=2m_{1/2}`$, $`\mathrm{tan}\beta `$ = 40 in the $`m_0`$-$`m_{1/2}`$ plane is plotted. The effect is to tilt (and narrow) the allowed dark matter band. The entire allowed parameter space can again be probed.

In Figures 4, 5 and 6, we have also drawn lines for various light Higgs masses (vertical dotted lines). A measurement of $`B_s\mu ^+\mu ^{}`$, $`a_\mu `$ and $`m_h`$ would then effectively determine the parameters of mSUGRA for $`\mu >0`$ by requiring that they intersect with the dark matter allowed band at a point. (If no choice of parameters allowed this, mSUGRA would be ruled out.) The Tevatron Run II should be able to either rule out a Higgs mass or give evidence for its existence at the 3 $`\sigma `$ level over the entire allowed mass range of SUSY light Higgs masses. Alternatively, the LHC’s determination of $`m_h`$ or the gluino mass (to determine $`m_{1/2}`$) would fix the parameters of mSUGRA.

The BNL E821 experiment should shortly have a more accurate value for $`a_\mu `$, with errors reduced by a factor of two or more. We note from the above figures the importance this result might have. Thus if $`a_\mu `$ increases, the $`a_\mu `$ bound moves downward, encroaching further on the allowed part of the parameter space, and a value of $`a_\mu \begin{array}{c}>\hfill \\ \hfill \end{array}50\times 10^{10}`$ would eliminate the mSUGRA model . However, if $`a_\mu `$ decreases significantly (but is still positive), the mSUGRA model would predict a heavy SUSY particle spectra closer to the TeV region, having significant effects on accelerator and dark matter detection physics. An accurate determination of $`a_\mu `$ corresponds to a line from upper left to lower right (or more precisely a band when errors are included) running parallel to the $`a_\mu <11\times 10^{10}`$ boundary, and cutting through the allowed dark matter band which runs from lower left to upper right. Thus these two experiments are complementary for determining the mSUGRA parameters.

6. R PARITY VIOLATING MODELS. We consider briefly here the case when R parity is broken. The general renormalizable R parity violating superpotential has the form

$$W_{R/_p}=\kappa _iL_iH_2+\lambda _{ijk}L_iL_jE_k^c+\lambda _{ijk}^{}L_iQ_jD_k^c+\lambda _{ijk}^{\prime \prime }U_i^cD_j^cD_k^c$$
(5)
The first three terms are lepton number violating, and the last term is baryon number violating. For the $`B_s\mu ^+\mu ^{}`$ decay, we need to consider lepton violating terms, and so we set $`\lambda ^{\prime \prime }`$ to zero (to prevent rapid proton decay). Unlike the R parity conserving models, the SUSY contribution to $`B_s\mu ^+\mu ^{}`$ can now occur at the tree level, and so can be considerably larger. For example, if the $`\lambda `$ and $`\lambda ^{}`$ terms are present in the theory, we can have

$$C_S=\frac{\lambda _{i23}^{}\lambda _{i22}^{}}{2m_{\stackrel{~}{\nu }}^2},C_P=\frac{\lambda _{i23}^{}\lambda _{i22}^{}}{2m_{\stackrel{~}{\nu }}^2}.$$
(6)
In Figure 7, we plot the $`Br`$ as a function of $`m_{1/2}`$ in the $`R`$ parity violating scenarios. We have chosen $`\lambda _{i23}^{}=\lambda _{i22}=0.02`$ and $`\mathrm{tan}\beta =10`$. These small couplings are not restricted by any another physical process. We see that in this case smaller values of $`\mathrm{tan}\beta `$ can be probed and the $`Br`$ is large compared to the R parity conserving case. The contributions of $`C_{S,P}`$ to the newly discovered decay mode $`B_sK\mu \mu `$ are small. For example, for $`m_{1/2}=300`$ GeV and $`m_0=300`$ GeV, $`Br[B_sK\mu \mu ]`$ gets a contribution of $`3.30\times 10^7`$. (This calculation, however, involves a large QCD uncertainty.) The observed branching fraction for this mode is $`(0.99_{0.30.14}^{+0.4+0.13})\times 10^6`$.

In contrast to the R conserving case, the R parity violating scenarios are not infested with tau’s and therefore has the potential to be observed at Run II directly. For example, let us consider the production of $`\chi _1^\pm \chi _2^0`$ and let us assume that $`\lambda _{i22},\lambda _{i23}^{}`$ are present. In this case we can have several interesting final states such as $`6l+2`$jets+missing $`E_T`$ (from neutrino) and $`4l+\mathrm{\hspace{0.17em}4}`$jet+missing $`E_T`$ from the production of $`\stackrel{~}{\chi }_1^\pm `$-$`\stackrel{~}{\chi }_2^0`$.

7. CONCLUSIONS. We have investigated the prospects for studying the rare decay mode $`B_s\mu ^+\mu ^{}`$ by the CDF detector in Run II at the Tevatron. We have analyzed the background for this process and find that a $`Br>1.2\times 10^8`$ for 15 fb<sup>-1</sup> can be probed at Run II. (This is nearly a factor of 100 improvement over the Run I bound.) In this background analysis, we have not only extrapolated the Run I limit, but also introduced new selection cuts to reduce the background. The improvement of the Run II detector has also been taken into account. If a similar analysis can be performed for the D0 detector, the combined data would increase the sensitivity for detecting this decay mode to $`Br>6.5\times 10^9`$. The LHC should be able to reach the sensitivity at the level of SM branching ratio.

The $`B_s\mu ^+\mu ^{}`$ is an important process for SUSY searches for new physics, as the Standard Model prediction of the branching ratio is quite small ($`3.5\times 10^9`$), and the SUSY contribution increases for large $`\mathrm{tan}\beta `$ as $`\mathrm{tan}^6\beta `$. For the mSUGRA model, the above sensitivity implies that Run II could probe a region of parameter space for $`\mathrm{tan}\beta >30`$, a region which could not be probed by a direct search at Run II. We have combined the expectations for $`B_s\mu ^+\mu ^{}`$ for mSUGRA with other experimental bounds on the parameter space. Thus a large branching ratio, $`i.e.`$, $`>7\times 10^8(14\times 10^8)`$ would be sufficient to eliminate the mSUGRA model for $`\mathrm{tan}\beta 50(55)`$. A measurement of $`B_s\mu ^+\mu ^{}`$ , the muon $`g2`$ anomaly and the light Higgs mass combined with the astronomical bounds on cold dark matter would essentially determine the mSUGRA model, allowing predictions of all other sparticle masses and cold dark matter neutralino-proton cross sections. If the muon $`g2`$ anomaly remains positive these cross sections are $`\begin{array}{c}>\hfill \\ \hfill \end{array}10^9`$ pb, and should then be accessible to future planned detectors such as GENIUS, Cryoarray, ZEPLIN IV and CUORE.

A simple model of R parity violation was also considered, and here the $`Br[B_s\mu ^+\mu ^{}]`$ can be large for both small and large $`\mathrm{tan}\beta `$.

ACKNOWLEDGEMENTS This work was supported in part by National Science Foundation grant PHY-0101015 and in part by Department of Energy grant DE-FG03-95ER40917. We should like to thank D. Toback for valuable comments.