text:
# Can the Two-Higgs-Doublet Model Survive the Constraint from the Muon Anomalous Magnetic Moment as Suggested ?

## I Introduction

Physicists have been delighted by the extraordinary success of the standard model (SM), while many got frustrated by the lack of experimental clues for the construction of the theory beyond. The most recent measurement on the muon anomalous magnetic moment $`a_\mu `$ BNL  revealed a plausible deviation as large as $`3\sigma `$ from the SM prediction, as suggested by the papers in Ref.amudata . While such a scenario is only a favorable, rather than unquestionably established, conclusion from the analyses, it has been taken by some physicists as a strong suggestion for physics beyond the SM. 111 This $`3\sigma `$ deviation was derived using the $`e^+e^{}`$ data. If the $`\tau `$ decay data are used, the deviation would be reduced to about $`0.9\sigma `$ amudata , which, however, has model dependence and thus less reliable.  Since many extensions of the SM are capable of giving rise to such a deviation, theorists, typically, would like to check the constraints imposed on the parameter space of a specific model. Such constraints are particularly interesting, because they are likely to give information not only on excluded regions but also on the predictions for where else should other evidences on the model be expected. For example, one expects to find a lower bound on the mass of a new particle playing the role on generating the extra contribution to $`a_\mu `$. The information has especially strong implications on models that have a small number of parameters which are already stringently constrained by various precision electroweak data. We present here such a case study, illustrating how far such a $`3\sigma `$ deviation can take us.

One of the simplest extensions of the SM is the two-Higgs-doublet model (2HDM) hunter , which adds one Higgs doublet in addition to the one required in the SM. A generic 2HDM allows flavor-changing neutral currents (FCNC), which can be avoided by restricting the couplings of the doublets, say, by imposing an ad hoc discrete symmetryglas . The most popular version, known as model II, has one Higgs doublet coupled to the down-type quarks (and charged leptons) and the second doublet to the up-type quarks. The physical content of the Higgs sector (assuming no CP violation) includes a pair of CP-even neutral Higgs bosons $`H`$ and $`h`$, a CP-odd neutral boson $`A`$, and a pair of charged-Higgs bosons $`H^\pm `$. The model fits in well with the criteria mentioned above for the $`a_\mu `$ result to have a very strong impact. We focus on the 2HDM II in this paper, answering the question of to what extent the model can survive a requirement of generating the $`3\sigma `$ deviation in $`a_\mu `$, as suggested.

The 2HDM II has been extensively studied in literature and tested experimentally. One of the most stringent tests is the radiative decay of $`B`$ mesons, specifically, the inclusive decay rate of $`bs\gamma `$, which has the least hadronic uncertainties. In the 2HDM, the rate of $`bs\gamma `$ can be enhanced substantially for large regions in the parameter space of the mass $`m_{H^\pm }`$ of the charged-Higgs boson and $`\mathrm{tan}\beta (=v_2/v_1`$, where $`v_1`$ and $`v_2`$ are the vacuum expectation values of the down- and up-sector Higgs doublets, respectively). An earlier analysis has already put a constraint on the charged-Higgs boson mass at $`m_{H^\pm }>380`$ GeVcdgg  (see also Ref.nlo ; misiak ). Updating the constraint while asking for the model to give rise to the $`a_\mu `$ deviation as suggested already imposes a strong and specific mass hierarchy between the pseudoscalar and the charged scalar. Electroweak precision data also have strong implications on the model.

The 2HDM can explain the muon anomalous magnetic moment deviation with a light pseudoscalar boson $`A`$ contributing via a 2-loop Barr-Zee-type diagram darwin ; ours . Our interest here is in updating a previous analysisours  and extending it to a comprehensive treatment of all the relevant constraints. The OPAL Collaboration opal  has recently published an update on their search for the Higgs bosons within the 2HDM framework. Since their result is more stringent than before, by combining the updated constraints from $`a_\mu `$ and other precision measurements, such as the $`\rho `$ parameter, $`R_b`$, and the $`bs\gamma `$ rate, together with this new OPAL result and a study on Yukawa processes from DELPHI delphi , we are able to limit the 2HDM to a tiny window of parameter space with an, perhaps, uncomfortably large value of $`\mathrm{tan}\beta `$. All in all, we will piece together a story of the very stringently constrained 2HDM, almost to the extent of killing it altogether. We note that there have been previous analyses on the 2HDM using the electroweak precision data maria-ew ; chan ; maria , to which we are partially in debt. Our study here can be considered an update, with a presentation along a different line.

We are here talking about the intricate interplay of a few stringent constraints on the overall parameter space of the model. It is not a simple matter to illustrate results on the space of a large number (six, as discussed in the next section) of parameters. In our opinion, the best way to do it may depend on how each of the constraints really works. Our presentation of the constraints follows what we consider the most efficient way to appreciate the overall results. It may be not very conventional, but is considered particularly illustrative. Wherever explicit plots are shown, we are typically plotting two- parameter fits of one or more constraints, based on the usual $`\chi ^2`$-analysis. We will show $`\chi ^2<4`$ regions, corresponding to $`2\sigma `$ deviation limits, which we consider as ”solution” regions — where the 2HDM survives the particular constraints. We will also show regions where the model fits better than the SM, wherever appropriate.

The organization is as follows. In the next section we briefly describe the 2HDM (II) and the relevant parameters used in our analysis. In Sec. III, we look at the $`bs\gamma `$ rate and the $`B^0\overline{B^0}`$ mixing, which require a heavy charged Higgs boson. In Sec. IV, we discuss the Higgs-sector contributions to $`a_\mu `$ and $`R_b`$, and show the strongly complementary character of the two data and how that plays off in the 2HDM. We give our fits to the data, with the exclusions from DELPHI and OPAL further imposed. In Sec. V, we add the consideration of the $`\rho `$ parameter. We conclude in Sec. VI.

## II Parameter Space of The Two Higgs Doublet model II

We take the parameter space of the model as given by a set of six Higgs-sector parameters

$$m_h,m_H,m_A,m_{H^\pm },\mathrm{tan}\beta ,\text{and}\alpha .$$

Here, the first four are masses of the physical Higgs states. The last one, $`\alpha `$, is the real scalar Higgs mixing angle as defined in Ref.hunter . We would like to emphasize that we are taking this set of six parameters as mutually independent experimental parameters. From the theoretical point of view, one has parameters in the scalar potential from which the above can be derived. However, without supersymmetry or anything else to avoid the hierarchy problem, the masses for the Higgs states suffer from quadratic divergences. The tree-level relations among parameters are modified substantially by loop corrections which depend on the renormalization approach and the cut-off imposed. To stay away from such uncertainties, we do not discuss the scalar potential here, except noting that there are enough degrees of freedom in the model in general, and especially with the loop corrections taken into account, to allow us to take the above six parameters as mutually independent. Here we only consider the 2HDM-II without CP-violation, and are not interested in couplings among the Higgs states. It is then obvious that we do not have to consider more than the six parameters. It should be noted that any substantial CP violation in the Higgs sector that may largely invalidate our analysis here is ruled out by the electron electric dipole moment constraint.

The four physical masses are direct experimentally measurable quantities. The other two parameters, $`\mathrm{tan}\beta `$ and the angle $`\alpha `$, come into the game as effective couplings. The Yukawa couplings of $`h,H`$, and $`A`$ to up- and down-type quarks are given by, with a common factor of $`igm_f/2M_W`$,

|  | $`t\overline{t}`$ | $`b\overline{b}`$ | $`\tau ^{}\tau ^+`$ |
| --- | --- | --- | --- |
| $`h`$: | $`\mathrm{cos}\alpha /\mathrm{sin}\beta `$ | $`\mathrm{sin}\alpha /\mathrm{cos}\beta `$ | $`\mathrm{sin}\alpha /\mathrm{cos}\beta `$ |
| $`H`$: | $`\mathrm{sin}\alpha /\mathrm{sin}\beta `$ | $`\mathrm{cos}\alpha /\mathrm{cos}\beta `$ | $`\mathrm{cos}\alpha /\mathrm{cos}\beta `$ |
| $`A`$: | $`i\mathrm{cot}\beta \gamma _5`$ | $`i\mathrm{tan}\beta \gamma _5`$ | $`i\mathrm{tan}\beta \gamma _5`$ |

while the charged Higgs $`H^{}`$ couples to $`t`$ and $`\overline{b}`$ via

$$\overline{b}tH^{}:\frac{ig}{2\sqrt{2}M_W}[m_t\mathrm{cot}\beta (1+\gamma _5)+m_b\mathrm{tan}\beta (1\gamma _5)].$$

From the perspective of our study here, the couplings given above may be considered as defining implicitly the two parameters $`\mathrm{tan}\beta `$ and $`\alpha `$. Other relevant couplings in our study are those to gauge bosons, as given by,

$`hZZ`$ $`:`$ $`igM_Z{\displaystyle \frac{\mathrm{sin}(\beta \alpha )}{\mathrm{cos}\theta _W}}g^{\mu \nu }`$
$`HZZ`$ $`:`$ $`igM_Z{\displaystyle \frac{\mathrm{cos}(\beta \alpha )}{\mathrm{cos}\theta _W}}g^{\mu \nu }`$
$`hAZ`$ $`:`$ $`g{\displaystyle \frac{\mathrm{cos}(\beta \alpha )}{2\mathrm{cos}\theta _W}}(pp^{})^\mu `$
$`HAZ`$ $`:`$ $`g{\displaystyle \frac{\mathrm{sin}(\beta \alpha )}{2\mathrm{cos}\theta _W}}(pp^{})^\mu `$
$`H^+H^{}Z`$ $`:`$ $`ig{\displaystyle \frac{\mathrm{cos}2\theta _W}{2\mathrm{cos}\theta _W}}(pp^{})^\mu ,`$
where $`p(h,H,H^+)`$ and $`p^{}(A,H^{})`$ are the 4-momenta going into the vertex.

## III Requirement of a heavy charged Higgs Boson

It has been well appreciated that $`B`$ physics bars the 2HDM from admitting a relatively light charged-Higgs state. We first review the constraints here. The first important constraint comes from the inclusive $`BX_s\gamma `$ result, and the second one comes from the $`B^0\overline{B^0}`$ mixing. The essential point here is that without a direct source of FCNC, the charged Higgs mediates the only significant contributions to flavor-changing processes, in addition to the $`W^\pm `$-mediated SM process. The experimental data then allows us to bound the charged Higgs mass independent of the other Higgs states.

The detail description of the effective Hamiltonian approach can be found in Refs. buras ; wise . Here we present the highlights that are relevant to our discussions. The effective Hamiltonian for $`BX_s\gamma `$ at a factorization scale of order $`O(m_b)`$ is given by

$$_{\mathrm{eff}}=\frac{G_F}{\sqrt{2}}V_{ts}^{}V_{tb}[\underset{i=1}{\overset{6}{}}C_i(\mu )Q_i(\mu )+C_{7\gamma }(\mu )Q_{7\gamma }(\mu )+C_{8G}(\mu )Q_{8G}(\mu )].$$
(1)
The operators $`Q_i`$ can be found in Ref.buras , of which the $`Q_1`$ and $`Q_2`$ are the current-current operators and $`Q_3Q_6`$ are QCD penguin operators. $`Q_{7\gamma }`$ and $`Q_{8G}`$ are, respectively, the magnetic penguin operators specific for $`bs\gamma `$ and $`bsg`$. Here we also neglect the mass of the external strange quark compared to the external bottom-quark mass. There have been more recent analyses misiak ; gambino ; urban  on $`bs\gamma `$ involving the NLO and other corrections, but the LO treatment here is sufficient for our purpose to put a lower bound on the charged Higgs mass $`m_{H^\pm }`$, which is then used in the central part of our analysis.

The decay rate of $`BX_s\gamma `$ normalized to the experimental semileptonic decay rate is given by

$$\frac{\mathrm{\Gamma }(BX_s\gamma )}{\mathrm{\Gamma }(BX_ce\overline{\nu }_e)}=\frac{|V_{ts}^{}V_{tb}|^2}{|V_{cb}|^2}\frac{6\alpha _{\mathrm{em}}}{\pi f(m_c/m_b)}|C_{7\gamma }(m_b)|^2,$$
(2)
where $`f(z)=18z^2+8z^6z^824z^4\mathrm{ln}z`$. The Wilson coefficient $`C_{7\gamma }(m_b)`$ is given by

$$C_{7\gamma }(\mu )=\eta ^{\frac{16}{23}}C_{7\gamma }(M_W)+\frac{8}{3}\left(\eta ^{\frac{14}{23}}\eta ^{\frac{16}{23}}\right)C_{8G}(M_W)+C_2(M_W)\underset{i=1}{\overset{8}{}}h_i\eta ^{a_i},$$
(3)
where $`\eta =\alpha _s(M_W)/\alpha _s(\mu )`$. The $`a_i`$’s and $`h_i`$’s can be found in Ref. buras . The coefficients $`C_i(M_W)`$ at the leading order in 2HDM II are given by

$`C_j(M_W)`$ $`=`$ $`0(j=1,3,4,5,6),`$ (4)
$`C_2(M_W)`$ $`=`$ $`1,`$ (5)
$`C_{7\gamma }(M_W)`$ $`=`$ $`{\displaystyle \frac{A(x_t)}{2}}{\displaystyle \frac{A(y_t)}{6}}\mathrm{cot}^2\beta B(y_t),`$ (6)
$`C_{8G}(M_W)`$ $`=`$ $`{\displaystyle \frac{D(x_t)}{2}}{\displaystyle \frac{D(y_t)}{6}}\mathrm{cot}^2\beta E(y_t),`$ (7)
where $`x_t=m_t^2/M_W^2`$, and $`y_t=m_t^2/m_{H^\pm }^2`$. The Inami-Lim functionsjapan  are given by

$`A(x)`$ $`=`$ $`x[{\displaystyle \frac{8x^2+5x7}{12(x1)^3}}{\displaystyle \frac{(3x^22x)\mathrm{ln}x}{2(x1)^4}}],`$ (8)
$`B(y)`$ $`=`$ $`y[{\displaystyle \frac{5y3}{12(y1)^2}}{\displaystyle \frac{(3y2)\mathrm{ln}y}{6(y1)^3}}],`$ (9)
$`D(x)`$ $`=`$ $`x[{\displaystyle \frac{x^25x2}{4(x1)^3}}+{\displaystyle \frac{3x\mathrm{ln}x}{2(x1)^4}}],`$ (10)
$`E(y)`$ $`=`$ $`y[{\displaystyle \frac{y3}{4(y1)^2}}+{\displaystyle \frac{\mathrm{ln}y}{2(y1)^3}}].`$ (11)
The most recent experimental data on $`bs\gamma `$ rate has been reported ichep-b , giving

$$B(bs\gamma )|_{\mathrm{exp}}=3.88\pm 0.36(\mathrm{stat})\pm 0.37(\mathrm{sys})_{0.28}^{+0.43}(\mathrm{theory}).$$

The most updated SM prediction is neubert

$$B(bs\gamma )|_{\mathrm{SM}}=(3.64\pm 0.31)\times 10^4,$$

which agrees very well the data. Both the experimental data and the SM prediction have been extrapolated to the total branching ratio. Therefore, there is only a little room for new physics contributions. The constraint on new physics contribution is, explicitly,

$$\mathrm{\Delta }B(bs\gamma )B(bs\gamma )|_{\mathrm{exp}}B(bs\gamma )|_{\mathrm{SM}}=(0.24_{0.59}^{+0.67})\times 10^4,$$
(12)
where we have added the various errors of the experimental data in quadrature. (Note that the theory error quoted in the experimental data is larger than the one quoted by the SM prediction. We take the more conservative value.)

The quantity that parameterizes the $`B^0\overline{B^0}`$ mixing is

$$x_d\frac{\mathrm{\Delta }m_B}{\mathrm{\Gamma }_B}=\frac{G_F^2}{6\pi ^2}|V_{td}^{}|^2|V_{tb}|^2f_B^2B_Bm_B\eta _B\tau _BM_W^2\left(I_{WW}+I_{WH}+I_{HH}\right),$$
(13)
whereabbott

$`I_{WW}`$ $`=`$ $`{\displaystyle \frac{x}{4}}[1+{\displaystyle \frac{39x}{(x1)^2}}+{\displaystyle \frac{6x^2\mathrm{log}x}{(x1)^3}}],`$
$`I_{WH}`$ $`=`$ $`xy\mathrm{cot}^2\beta [{\displaystyle \frac{(4z1)\mathrm{log}y}{2(1y)^2(1z)}}{\displaystyle \frac{3\mathrm{log}x}{2(1x)^2(1z)}}+{\displaystyle \frac{x4}{2(1x)(1y)}}],`$
$`I_{HH}`$ $`=`$ $`{\displaystyle \frac{xy\mathrm{cot}^4\beta }{4}}[{\displaystyle \frac{1+y}{(1y)^2}}+{\displaystyle \frac{2y\mathrm{log}y}{(1y)^3}}],`$
with $`x=m_t^2/M_W^2`$, $`y=m_t^2/m_{H^\pm }^2`$, $`z=M_W^2/m_{H^\pm }^2`$, and the running top mass $`m_t=m_t(m_t)=166\pm 5`$ GeV. The experimental value is pdg2002

$$x_d=0.755\pm 0.015.$$
(14)
We use the following input parameters pdg2002 : $`|V_{tb}V_{td}^{}|=0.0079\pm 0.0015`$, $`f_B^2B_B=(198\pm 30\mathrm{GeV})^2(1.30\pm 0.12)`$, $`m_B=5279.3\pm 0.7`$ MeV, $`\eta _B=0.55`$, and $`\tau _B=1.542\pm 0.016`$ ps. Note that the value of $`|V_{tb}V_{td}^{}|`$ is in fact determined by the measurement of $`x_d`$. Now we can use the data to constrain the new contribution from the charged-Higgs boson.

The two constraints discussed above are quite stringent, giving a lower bound on $`m_{H^\pm }`$ close to $`500\text{GeV}`$ at 95% C.L. for intermediate and large values of $`\mathrm{tan}\beta `$. The result is illustrated in Fig. 1. The heavy charged Higgs mass means that the state pretty much decouples, playing a little role in the contributions of the 2HDM to quantities like $`a_\mu `$ and $`R_b`$, which we turn to in the next section.

## IV $`𝒂_𝝁`$ Vs $`𝑹_𝒃`$

The most recent data on the $`a_\mu `$ indicates amudata  (see also the footnote # 1)

$$\mathrm{\Delta }a_\mu a_\mu ^{\mathrm{exp}}a_\mu ^{\mathrm{SM}}=(33.9\pm 11.2)\times 10^{10}.$$
(15)
The result shows a $`3\sigma `$ deviation to be explained by new physics. Adopting the view that the $`a_\mu `$ problem is real and demands new physics contributions, we will see that it has a strong and definite implication on the Higgs spectrum of the 2HDM. In fact, we had performed an analysis ours  along the line for the earlier data. The major point is that a light pseudoscalar, together with a large $`\mathrm{tan}\beta `$ value, is required to explain the positive $`\mathrm{\Delta }a_\mu `$ contribution via a two-loop Barr-Zee diagram BZ . A real scalar contributes in the negative direction. To avoid a cancellation, the real scalar mass has to be heavy. We will show here that such a mass splitting is strongly disfavored by the allowed contribution to $`R_b`$, for which the experimental data agrees well with the SM prediction. The two constraints are hence strongly complementary.

It has been emphasized in Ref.darwin ; ours  that for the Higgs boson mass larger than about $`3\text{GeV}`$, the dominant Higgs contributions to $`a_\mu `$ actually come from the two-loop Barr-Zee diagram with a heavy fermion ($`f`$) running in the upper loop. A $`m_f^2/m_\mu ^2`$ factor could easily overcome the $`\alpha /4\pi `$ loop factor. In our calculation here, we include all one-loop contributions and all two-loop Barr-Zee-type contributions with an internal photon and a third-family fermion running in the loop. The latter diagrams with the bottom and tau loops are strongly enhanced by $`\mathrm{tan}\beta `$. If the internal photon was replaced by a $`W^\pm `$ or a $`Z^0`$, the contributions will be much suppressed (see edm , for examples). The $`W^\pm `$ case is in particular strongly suppressed, partly as a result of the fact that the Higgs boson has to be $`H^\pm `$, the mass of which we have shown above to be heavier than 500 GeV. The only other important diagrams of the Barr-Zee type are the SM diagrams, such as the one with the $`W^\pm `$ replacing the heavy fermion. We neglect the small “extra” contributions from such diagramswy  because of the small difference between the Higgs boson mass used here and that used in the Ref. M .

Explicitly, we first write the fermion couplings of a neutral Higgs mass eigenstate $`\varphi ^0`$ as

$$^{\overline{f}\varphi ^0f}=\lambda _f\frac{m_f}{v}\overline{f}\varphi ^0f+i\gamma _5A_f\frac{m_f}{v}\overline{f}\varphi ^0f,$$
(16)
where $`\lambda _f\frac{m_f}{v}`$ and $`A_f\frac{m_f}{v}`$ are the effective scalar and pseudoscalar couplings explicitly given in Sec. II, and $`v=246\text{GeV}`$. The two-loop photon Barr-Zee diagram contribution from $`\varphi ^0`$, with a heavy fermion $`f`$ running in the second loop, is given by

$$\mathrm{\Delta }a_\mu ^\varphi =\frac{N_c^f\alpha _{\text{em}}}{4\pi ^3v^2}m_\mu ^2𝒬_f^2\left[A_\mu A_fg\left(\frac{m_f^2}{m_\varphi ^2}\right)\lambda _\mu \lambda _ff\left(\frac{m_f^2}{m_\varphi ^2}\right)\right],$$
(17)
where

$`f(z)={\displaystyle \frac{1}{2}}z{\displaystyle _0^1}𝑑x{\displaystyle \frac{12x(1x)}{x(1x)z}}\mathrm{ln}{\displaystyle \frac{x(1x)}{z}},`$
$`g(z)={\displaystyle \frac{1}{2}}z{\displaystyle _0^1}𝑑x{\displaystyle \frac{1}{x(1x)z}}\mathrm{ln}{\displaystyle \frac{x(1x)}{z}};`$ (18)
$`N_c^f`$ represents the number of color degrees of freedom in $`f`$, and $`𝒬_f`$ its electric charge. Here, we have three scalars and no CP violating mixing is assumed. The real scalars $`h`$ and $`H`$ give negative contributions only (from the second part) while a pseudoscalar $`A`$ gives a positive contribution only (from the first part). The diagrams with the $`b`$ and $`\tau `$ loops are $`\mathrm{tan}\beta `$ enhanced.

We plot in Fig. 2 the $`2\sigma `$ range of solution to $`a_\mu `$ on the plane of the pseudoscalar mass $`m_A`$ versus $`\mathrm{tan}\beta `$, considering only the pseudoscalar contribution. Note that while the region below the solution band is excluded, the solution in the region above the band may be admissible when the $`a_\mu `$ is compensated by some negative contributions from the real scalar(s). We also superimpose on the plot the excluded region from the DELPHI study on Higgs Yukawa processes in the $`4b`$ and $`2b2\tau `$ final states delphi . We can see that one obtains lower bounds on $`m_A`$ and $`\mathrm{tan}\beta `$ as $`26\text{GeV}`$ and $`30`$ respectively. 222The plot is similar to the one given in Ref. maria , in which some more constraints are superimposed. We include here only the important ones. Note that the Tevatron exclusion region claimed in the paper is not used here. The exclusion result was from Ref.roco , which is an analysis based on the minimal supersymmetric standard model. The result should not be directly applicable to the present case of 2HDM. We do not find a similar study on the Tevatron data based on the 2HDM. We had given in Ref. ours  plots of the solution regions on the $`m_h`$-$`m_A`$ plane for the old $`a_\mu `$ data with specific values of $`\mathrm{tan}\beta `$ and the scalar mixing angle $`\alpha `$. We will present similar results here with, however, the complementary $`R_b`$ constraint included. We will see that not much area of the parameter space can survive the combination of both $`a_\mu `$ and $`R_b`$ data.

The current $`R_b`$ measurement is given by lep-ew

$$R_b^{\mathrm{exp}}=0.21646\pm 0.00065.$$

With $`R_b^{\mathrm{SM}}=0.215768`$, we have

$$\mathrm{\Delta }R_bR_b^{\mathrm{exp}}R_b^{\mathrm{SM}}=0.000692\pm 0.00065.$$
(19)
The $`\mathrm{\Delta }R_b`$ contributions in the 2HDM are given, for example, by formulas in Ref.denner . The charged Higgs contribution is always negative. On the other hand, there is a window of parameter space for the neutral Higgs contributions to be positive. From our discussions above, we need a light pseudoscalar and have to live with a heavy charged Higgs. We are therefore more interested in the neutral Higgs contributions.

Let us first focus on the contributions from the lighter Higgs bosons, the pseudoscalar $`A`$ and the real scalar $`h`$, pushing the other scalar $`H`$ to the heavy-mass limit together with $`H^\pm `$. The scenario will actually be well justified by our discussion on the constraint of the $`\rho `$ parameter in the next section. The smallness of $`\mathrm{\Delta }R_b`$ contributions admitted here generally disfavors a large mass splitting between $`m_h`$ and $`m_A`$. This is in contrast to the requirement of a positive contribution to $`a_\mu `$. Since the SM result (with $`M_{H_{SM}}`$ at $`115\text{GeV}`$) now represents a $`3\sigma `$ deviation in $`a_\mu `$, better fits to the combined $`a_\mu `$-$`R_b`$ data are possible from the 2HDM. We illustrate some such fits in Fig. 3. In the figure, we take the case of $`\mathrm{tan}\beta =58`$ and check various values of the Higgs mixing angle $`\alpha `$. Here, and in the discussion below unless specifically stated otherwise, we stick to $`m_{H^\pm }=500\text{GeV}`$ and $`m_H=1\text{TeV}`$. The exact value of $`m_H`$ does not matter at all here, one, however, should note that the charged Higgs boson still gives a contribution of $`2.32\times 10^4`$ to $`R_b`$, which is about $`\frac{1}{3}\sigma `$ in strength. Bearing this in mind, it is easy to estimate from our plots the slight shift in each of the admissible region as $`m_{H^\pm }`$ is being pushed towards the decoupling limit.

Each of the plots in Fig. 3 gives $`\pm 2\sigma `$ limits for $`a_\mu `$ and $`R_b`$ fits, with darker shaded regions indicating the solution of interest defined by a total $`\chi ^2`$ of $`4`$ or less. Also marked in the plots are regions with a total $`\chi ^2`$ less than the SM value of $`10.3`$ (the sum of $`a_\mu `$ and $`R_b`$). The purpose of showing the area with a total $`\chi ^2<\chi ^2(\mathrm{SM})`$ is to indicate the region of parameter space that can fit the $`a_\mu `$ and $`R_b`$ better than the SM, other than the decoupling limits. From now on, we concentrate on the dark area of a total $`\chi ^2<4`$ as a valid solution to the $`a_\mu `$ and $`R_b`$ data. We can see that there are no solutions for $`\frac{\pi }{8}<\alpha <\frac{\pi }{4}`$. While a larger magnitude of the $`\alpha `$ looks more favorable, the best solution, with higher Higgs masses ($`m_h`$ especially), stay close to $`|\mathrm{sin}\alpha |=1`$, inclining more towards the negative sign. The most favorable range is around $`\frac{\pi }{2}<\alpha <\frac{3\pi }{8}`$ \[cf. plots (c) and (d)\]. The plots in Fig. 3 illustrate well the trend of the changes in the $`a_\mu `$-$`R_b`$ solution regions with variations in $`\alpha `$, which is quite generic for $`\mathrm{tan}\beta `$ around and larger than 50. A higher $`m_h`$ solution is preferred as the light mass solutions are easily killed by searches at LEP, including the DELPHI exclusion used in Fig. 2 and a particularly focused 2HDM analysis from OPALopal , to be discussed below. In fact, as we will illustrate below, if any part of the solutions to the $`a_\mu `$-$`R_b`$ fits survive the exclusions from LEP, it is more or less the part with a high enough $`m_h`$ value. 333In these plots, we used the standard $`R_b`$ formulas for 2HDM as available in the reference denner  quoted. The formulas have not taken into consideration the possible tree-level decays of the $`Z^0`$ boson into a pair of Higgs states. Such very light Higgs regions are typically not of interest to theorists, as they would be easily excluded by direct searches. By the same token, we use the formulas as they are and show the $`a_\mu `$-$`R_b`$ fits without bothering about the tree-level modifications of $`R_b`$ from such Higgs channels, leaving such fake solutions to be taken care of by the exclusions from LEP searches to be discussed below. We thank A. Sanda for alerting us to explicitly address the issue here.

For instance, the DELPHI exclusion we used in Fig. 2 obviously kills quite a part of the above solutions. At $`\mathrm{tan}\beta =58`$, the lower bound on $`m_A`$ is actually about $`40\text{GeV}`$. It is particularly interesting to check the case of a relatively small $`\mathrm{tan}\beta `$. We show in Fig. 4 the case of $`\mathrm{tan}\beta =40`$, in which the DELPHI almost kills all solutions. In the figure, only a tiny window survives at $`\mathrm{sin}\alpha =\frac{3\pi }{8}`$. This is, however, at a $`m_h`$ (as well as $`m_A`$) value too small to survive the OPAL exclusion discussed below. From the same study by DELPHI, the Yukawa processes were also used to impose bounds on the $`m_h`$delphi . Note that the plots in the Ref. delphi  give bounds on Higgs masses versus the $`b\tau `$ coupling enhancement factor, which is simply $`\mathrm{tan}\beta \sqrt{B(Ab\overline{b},\tau ^+\tau ^{})}`$ for the pseudoscalar case, but has an added $`|\mathrm{sin}\alpha |`$ dependence for the case of the scalar $`h`$. In the range of $`m_A`$ values of interest, however, the $`m_h`$ bound is typically superseded by the OPAL exclusion. Note that the $`a_\mu `$ solution with $`m_h`$ at a few GeV is also inadmissible, from the consideration of $`\mathrm{{\rm Y}}`$-decayours  and otherwise. Figure 4 also illustrates that a positive value of $`\alpha `$, around $`\frac{3\pi }{8}`$ tends to give $`a_\mu `$-$`R_b`$ solutions with the largest $`m_A`$. This is mainly a result of the $`a_\mu `$ constraint. At a fixed $`m_h`$, the contribution of the scalar through the 2-loop Barr-Zee diagram to $`a_\mu `$ is suppressed by $`|\mathrm{sin}\alpha |`$, and thus allows a larger $`m_A`$. A slight asymmetry in the cases of positive and negative $`\alpha `$ values comes in as a result of the different way the 1-loop contributions go.

We should also point out that we find no solutions to the $`a_\mu `$-$`R_b`$ fit at all for much lower $`\mathrm{tan}\beta `$. The $`a_\mu `$ solutions shown in Fig. 2 simply produce a $`m_h`$-$`m_A`$ splitting too large to accommodate the $`R_b`$ constraint. Therefore, the solutions to the $`a_\mu `$-$`R_b`$ fits start to emerge as $`\mathrm{tan}\beta `$ approaches a large enough value, not much below 40. To get solutions with high enough masses for the Higgs bosons so as to survive the exclusion limits from the LEP searches, one will have to get to a higher and higher $`\mathrm{tan}\beta `$ value. To check the details, we first turn to the powerful exclusions from OPAL.

Based also on the LEP data, OPAL has been publishing Higgs-search analyses specifically focused on the 2HDM. Here, we use the most recent results available opal ; pamela . The results exclude a region at the lower left corner of the $`m_h`$-$`m_A`$ plane for each specific value of the mixing angle $`\alpha `$ (four explicitly shown). As presented in Ref.opal , however, the results are not $`\mathrm{tan}\beta `$ specific. We show in Fig. 5 the solution regions from $`a_\mu `$-$`R_b`$ with the OPAL and the above mentioned DELPHI excluded regions superimposed. Here, we show the cases of $`\mathrm{tan}\beta =50`$ and $`58`$ for two $`\alpha `$ values, $`\frac{\pi }{4}`$ and $`\pm \frac{\pi }{2}`$. The latter are chosen as they are among the $`\alpha `$ values for which the OPAL paper gives the explicit exclusion. The exclusion for non-specific $`\alpha `$ values only presents a substantially weakened result to be of interest here. The two $`\alpha `$ values are also close to or within the range of value for an optimal $`a_\mu `$-$`R_b`$ fit, as illustrated in Fig. 3 above.

In fact, the OPAL exclusions are based on the $`\mathrm{tan}\beta `$ value in the range $`158`$, given only at four values of $`\alpha `$, which we adopted here for the plots. An excluded point is one excluded at all value of $`\mathrm{tan}\beta `$ within the range. The excluded ranges may hence be extended at each specific value of $`\mathrm{tan}\beta `$, with more detailed analysis of the datapamela2 . In fact, one would expect the enhanced couplings at a larger $`\mathrm{tan}\beta `$ generally push the exclusion regions towards higher masses. One can see in the plots that the $`a_\mu `$-$`R_b`$ solution regions are largely cut off by the presented OPAL exclusions in general. Actually, nothing survives in the illustrated plots for $`\mathrm{tan}\beta =50`$ \[cf. plots (a) and (b) of Fig. 5\]; and it is also quite obvious that the same is true for the case of $`\mathrm{tan}\beta =40`$. Figures 4 and 5 together clearly illustrate the general trend of how the increase in $`\mathrm{tan}\beta `$ gives better results. While it does not do much for the case of $`\alpha =\frac{\pi }{4}`$, at $`\alpha =\pm \frac{\pi }{2}`$, the admissible $`m_h`$ and $`m_A`$ values are pushed to be high enough to escape the OPAL exclusion at $`\mathrm{tan}\beta >50`$ \[cf. plot (c) of Fig. 5\]. However, we are not bold enough to say for sure if no solution survives at $`\mathrm{tan}\beta =50`$ though. The optimal case giving the largest $`m_h`$ value for the $`a_\mu `$-$`R_b`$ is likely to be around $`\alpha =\frac{3\pi }{8}`$, at which the present result of OPAL did not show its greatest strength.

We pick $`\mathrm{tan}\beta =58`$ for the above detailed illustration of the $`a_\mu `$-$`R_b`$ fits because it gives the most favorable case within the limit of a direct application of the OPAL exclusion. There is apparently a surviving window in the parameter space. Combining our present study with a refined version of the OPAL analysis to focus on our $`a_\mu `$-$`R_b`$ solution regions, especially the upper blots as shown in plots (a) and (c) of Fig. 5, will certainly be very interesting. Exclusions have to be checked for each specific value of $`\mathrm{tan}\beta `$ and that of $`\alpha `$. Such a study will further narrow down the survival parameter space regions of the 2HDM, especially with further improved exclusion limits. We are told that the LEP data actually allows such an improvementpamela2 . A very tantalizing question is if the current constraints are actually strong enough to kill the model altogether!

The OPAL analysis is limited to $`\mathrm{tan}\beta `$ value at or below $`58`$, while the DELPHI analysis stops at $`100`$. The very large $`\mathrm{tan}\beta `$ region is theoretically unfavorable and may provide practical problems due to the much enhanced $`b`$-quark Yukawa coupling, which signals a breakdown of the perturbative treatment. Nevertheless, we will include some results from such uncomfortably large $`\mathrm{tan}\beta `$ values, and urge the OPAL group to push on a bit further in their analysis.

For $`\mathrm{tan}\beta >58`$, we again illustrate some results for $`\alpha =\frac{\pi }{4}`$ and $`\pm \frac{\pi }{2}`$ in Fig. 6, in which we still put in the (no longer exactly valid) OPAL exclusion from the $`\mathrm{tan}\beta =158`$ scan for reference. For the $`\alpha =\frac{\pi }{4}`$ case, the $`a_\mu `$-$`R_b`$ solution regions never rise above $`m_h=60\text{GeV}`$ (as also for the smaller $`\mathrm{tan}\beta `$ cases), and partially excluded by DELPHI. There seem to be good enough reasons to believe that such low Higgs mass regions should be excluded by the available LEP data if an analysis along the OPAL line is performed. The $`\alpha =\pm \frac{\pi }{2}`$ case is better. The surviving region seen at $`\mathrm{tan}\beta =58`$ moved further to the right, towards larger $`m_A`$, and further up a bit for larger $`m_h`$. The rise in $`m_h`$ actually more or less saturates at $`\mathrm{tan}\beta =100`$, and falls for even larger values of $`\mathrm{tan}\beta `$. Such a region still lives at the boundary of the DELPHI exclusion, but would still have a surviving part even if the OPAL exclusion can still be imposed. The region is, however, shrinking with increases in $`\mathrm{tan}\beta `$, due to a more fine-tuned $`a_\mu `$-$`R_b`$ solution. One should also bear in mind that the (upward) shifting, and hence slight enlargement, of this solution region with a further increase in the $`m_{H^\pm }`$ value. For this purpose, we give the charged-Higgs contribution to $`R_b`$ in Fig. 7.

Let us summarize our results so far. We have seen that combining the suggested requirement of producing a definite positive contribution to $`a_\mu `$ while keeping a limited deviation from the SM $`R_b`$ result is an extremely stringent constraint on the parameter space of the 2HDM. When the available direct experimental search results from the LEP experiments are further implemented, there is at most a tiny window of parameter space that can survive. The apparently surviving region from the above discussions is restricted to very large values of $`\mathrm{tan}\beta `$. In fact, it may be already uncomfortably large, inviting the problem of the perturbativity of the Yukawa coupling of the $`b`$ quark. As for the mixing angle $`\alpha `$, it is being pushed close to the $`\frac{\pi }{2}<\alpha <\frac{3\pi }{8}`$ region. All these are based on a strong mass splitting between the charged Higgs and the light scalars, with the pseudoscalar $`A`$ lying typically below $`80\text{GeV}`$ and the scalar $`h`$ below $`140\text{GeV}`$.

Recall that we have essentially decoupled the heavy Higgs boson $`H`$ in the above analysis. We promised to justify this as a physics requirement. We first noted that a not too heavy $`H`$ would more or less add to the effect of the other real scalar $`h`$ in the contributions to $`a_\mu `$ and $`R_b`$. So, we expect it to ask for a larger mass splitting when fitting $`a_\mu `$ is concerned, but a smaller mass splitting to fit $`R_b`$. In another word, further tightening of the apparent solution window. In the section below, we will show that fitting another precision EW parameter, the $`\rho `$ parameter, actually does require sending $`m_H`$ to a very large value, indeed well beyond $`m_{H^\pm }`$.

## V The $`\rho `$ parameter constraint

The parameter $`\rho `$ was introduced to measure the relation between the masses of $`W^\pm `$ and $`Z^0`$ bosons. In the SM $`\rho M_W^2/M_Z^2\mathrm{cos}^2\theta _W=1`$ at tree-level. However, the $`\rho `$ parameter receives contributions from the SM corrections and from new physics. The deviation from the SM prediction is usually described by the parameter $`\rho _0`$ defined bylangacker

$$\rho _0\frac{M_W^2}{\rho M_Z^2\mathrm{cos}^2\theta _W},$$
(20)
where the $`\rho `$ in the denominator absorbs all the SM corrections, including the corrections from the top quark and the SM Higgs boson. By definition, $`\rho _0=1`$ in the SM. Sources of new physics that contribute to $`\rho _0`$ can be written as

$$\rho _0=1+\mathrm{\Delta }\rho _0^{\mathrm{new}},$$
(21)
where $`\mathrm{\Delta }\rho _0^{\mathrm{new}}=\mathrm{\Delta }\rho ^{\text{2HDM}}\mathrm{\Delta }\rho ^{\text{SM-Higgs}}`$ in our case. Note that since the two-doublet Higgs sector (in the 2HDM) is employed here to replace the SM Higgs, the latter contribution to $`\mathrm{\Delta }\rho `$ has to be subtracted out.

The most recent reported value of $`\rho _0`$ is pdg2002

$$\rho _0=1.0004\pm 0.0006,\text{(with }M_{H_{\mathrm{SM}}}\text{ fixed at 115 GeV)}.$$
(22)
In terms of new physics the constraint becomes:

$$\mathrm{\Delta }\rho _0^{\mathrm{new}}=0.0004\pm 0.0006.$$
(23)
In 2HDM $`\mathrm{\Delta }\rho `$ receives contributions from all Higgs bosons given by hunter ; maria-ew

$`\mathrm{\Delta }\rho ^{\text{2HDM}}`$ $`=`$ $`{\displaystyle \frac{\alpha _{\mathrm{em}}}{4\pi \mathrm{sin}^2\theta _WM_W^2}}[F(m_A,m_{H^+})+\mathrm{cos}^2(\beta \alpha )[F(m_{H^+},m_h)F(m_A,m_h)]`$ (24)
$`+\mathrm{sin}^2(\beta \alpha )[F(m_{H^+},m_H)F(m_A,m_H)]]`$
$`+\mathrm{cos}^2(\beta \alpha )\mathrm{\Delta }\rho ^{\text{SM}}(m_H)+\mathrm{sin}^2(\beta \alpha )\mathrm{\Delta }\rho ^{\text{SM}}(m_h),`$
where

$`F(x,y)`$ $`=`$ $`{\displaystyle \frac{1}{8}}x^2+{\displaystyle \frac{1}{8}}y^2{\displaystyle \frac{1}{4}}{\displaystyle \frac{x^2y^2}{x^2y^2}}\mathrm{log}\left({\displaystyle \frac{x^2}{y^2}}\right)=F(y,x),`$
$`\mathrm{\Delta }\rho ^{\text{SM}}(M)`$ $`=`$ $`{\displaystyle \frac{\alpha _{\mathrm{em}}}{4\pi \mathrm{sin}^2\theta _WM_W^2}}\left[3F(M,M_W)3F(M,M_Z)+{\displaystyle \frac{1}{2}}(M_Z^2M_W^2)\right].`$ (25)
Let us take a closer look into the implication of the formulas above. First of all, we note that $`\mathrm{\Delta }\rho ^{\text{SM}}(M)`$ has a negative value with magnitude increasing with $`M`$. As to be expected from above, the value is about $`0.0004`$ at $`M=115\text{GeV}`$. It has a relatively mild variation, and does not go beyond $`0.005`$ even as $`M`$ gets to $`10\text{TeV}`$. In direct contrast, the other contributions to $`\mathrm{\Delta }\rho ^{\text{2HDM}}`$ in the above formula are very sensitive to the masses involved. The $`F(x,y)`$ function is always positive, vanishes only at $`x=y`$, and increases with a faster and faster rate with the splitting between $`x`$ and $`y`$. In a typical scenario that is of interest here, we expect the pseudoscalar to be the lightest Higgs state with a quite heavy charged Higgs. That makes the contribution from the first term \[involving $`F(m_A,m_{H^+})`$\] large; indeed of order $`0.01`$ for $`m_{H^+}`$ satisfying the lower bound from $`bs\gamma `$ and $`B\overline{B}`$ mixing. To get the required almost zero value of $`\mathrm{\Delta }\rho ^{\text{2HDM}}`$, we need some negative contributions from the terms involving $`[F(m_{H^+},m_h)F(m_A,m_h)]`$ and $`[F(m_{H^+},m_H)F(m_A,m_H)]`$. And the solution is obviously a fine-tuned one. Consider a case of degenerated Higgs real scalars, $`m_H=m_h=M`$. The $`(\beta \alpha )`$ dependence in the formula is removed as the sine-square and cosine-square parts are combined. We need $`[F(m_{H^+},M)F(m_A,M)]`$ to have a negative value close in magnitude to that of $`F(m_A,m_{H^+})`$. The former is obviously positive roughly when $`M`$ is closer to $`m_A`$ than $`m_{H^+}`$, and negative when it is the other way round; and it can be larger than $`F(m_A,m_{H^+})`$ only for $`M>m_{H^+}`$. The $`(\beta \alpha )`$ dependence comes back in the generic situation, with a mass splitting between the two real scalars. Obviously, at least one of them, by definition $`H`$, has to be heavy. However, they cannot be both heavy, because the $`R_b`$ constraint does not allow only a light pseudoscalar giving a substantial contribution. Hence, this additional requirement of a limited splitting between $`m_A`$ and $`m_h`$ clearly suggests a large $`m_H`$, typically larger than $`m_{H^+}`$. The larger the $`m_H`$ value, the smaller the $`\mathrm{sin}^2(\beta \alpha )`$ is required. Admissible solutions, though fine-tuned, can be obtained so long as the required $`\mathrm{sin}^2(\beta \alpha )`$ falls into the legitimate interval. A large splitting between $`m_H`$ and $`m_h`$ also resulted in a very narrow range of admissible $`(\beta \alpha )`$ values, and hence the $`\alpha `$ values at a fixed $`\mathrm{tan}\beta `$ of interest.

Our numerical results corroborate well with the above analytical discussions. We illustrative our discussion with a plot in Fig. 8. Here, we take a “surviving” solution point to the $`a_\mu `$-$`R_b`$ fits and perform a further fitting together with the $`\rho `$ parameter by varying $`m_H`$ and $`m_{H^+}`$. The extremely fine-tuned nature of the solution is well-illustrated by the very narrow $`\chi ^2`$ bands. Moreover, $`m_H`$ is always more than twice of $`m_{H^+}`$ for an $`|\mathrm{sin}\alpha |`$ larger than $`0.8`$. Note that $`\mathrm{sin}\alpha `$ of $`0.8`$ and $`0.92`$ roughly correspond to an $`\alpha `$ value of $`\frac{3\pi }{10}`$ and $`\frac{3\pi }{8}`$, respectively. The basic features remain if other solution points are taken, hence, we refrain from showing more plots.

## VI Conclusions

The requirement for the 2HDM II to give rise to the suggested $`+3\sigma `$ deviation in the muon anomalous magnetic moment demands a light pseudoscalar, preferably around $`40\text{GeV}`$ for $`\mathrm{tan}\beta 58`$. The FCNC constraints, in particular the $`bs\gamma `$, require a heavy charged Higgs beyond $`400\text{GeV}`$. This large mass splitting in the Higgs mass spectrum is difficult to be accommodated by the precision EW measurements. In particular, the $`\rho `$ parameter constraint then admits only very fine-tuned solutions, with cancellation from opposite contributions good to one in a hundred, favoring heavy real scalars. The fine-tuned nature of the solutions, while making many physicists uncomfortable with the model, is not in itself a good enough reason to pronounce the death of the model. If one has reasons to be confident about the correctness of the model, one would say that the available constraints are just strong enough to pin down for us the values of the unknown model parameter. A good example of such a situation is given by the pinning down of the top mass value from precision EW data prior to the experimental discovery of the top quark.

Nevertheless, we have not had much of a reason to believe in the correctness of the 2HDM. In fact, in our analysis here, we focus more on the simultaneous fits to $`\mathrm{\Delta }a_\mu `$ and $`R_b`$. Having only the pseudoscalar contribution dominated the corrections to $`R_b`$ is fatal. Hence, we require a relatively light $`m_h`$, while pushing $`m_H`$ to way beyond $`m_{H^+}`$ in order to satisfy the $`\rho `$ parameter. Even then, not much of a solution to the $`a_\mu `$-$`R_b`$ fits survives the direct search exclusion. Here, we require a total $`\chi ^2`$ of 4 or smaller for the $`a_\mu `$-$`R_b`$ fits to claim a good solution. For roughly $`\mathrm{tan}\beta <40`$, no solution to $`a_\mu `$-$`R_b`$ fits survives the DELPHI exclusion. In fact, there is no solution to the $`a_\mu `$-$`R_b`$ fit for a $`\mathrm{tan}\beta `$ value quite a bit smaller than 40. For $`\mathrm{tan}\beta `$ value from around 40 to a bit beyond 50, solutions surviving the DELPHI exclusions exist, but only to be killed by the OPAL exclusions. For even larger $`\mathrm{tan}\beta `$, the $`a_\mu `$-$`R_b`$ solutions surviving both the DELPHI and OPAL exclusions started to emerge. This is very much restricted to an $`\alpha `$ value in the range $`\frac{\pi }{2}<\alpha <\frac{3\pi }{8}`$. Solution regions shrink fast outside the range as the $`m_h`$ values given by the $`a_\mu `$-$`R_b`$ fits drop towards OPAL exclusion bound.

In summary, under the strong restriction of the available constraints, we show that only a very tiny window of apparent solutions exist close to the limit of $`\mathrm{tan}\beta 58`$. If the OPAL group could tailor their analysis of the LEP data to focus on the apparent solution window as shown here, we would be able to have a more definite conclusion. It looks to us that the solution window will be shut down quite substantially. So, we have “almost enough” constraints to kill the model altogether. For $`\mathrm{tan}\beta `$ beyond the $`58`$ limit, our hands are tied at the moment by the unavailability of the strong LEP exclusion results as presented by OPAL. The very large $`\mathrm{tan}\beta `$ values certainly make many of us uncomfortable though. It is theoretical undesirable as limited by the blowing up of the bottom Yukawa couplings.

## Acknowledgments

We thank D. Chang, C. Kao, W.-Y. Keung, and A. I. Sanda for valuable discussions. We are also in debt to P. Ferrari from OPAL and M. Boonekamp from DELPHI for clarification about their results. This research was supported in part by the National Science Council grant number NSC91-2112-M-008-042 (O.K.) and by the NCTS.

Otto Kong thanks the hospitality of National Center for Theoretical Sciences, Taiwan during the early phase of the study, and that of KEK Theory Group, Japan where the final draft is finished.


format: text

text:
# A high resolution foreground cleaned CMB map from WMAP

## I INTRODUCTION

The release of the first results BennettWMAP ; BennettMission ; BennettForegs ; barnes1-03 ; barnes2-03 ; hinshawpower03 ; hinshaw03 ; jarosik1-03 ; jarosik2-03 ; kogut03 ; komatsu03 ; limon03 ; page1-03 ; page2-03 ; Page3-03 ; peiris03 ; Spergel03 ; verde03  from the Wilkinson Microwave Anisotropy Probe (WMAP) constituted a major milestone in cosmology, laying a solid foundation upon which to found the cosmological quest in coming years. Although much of the attention in the wake of the WMAP release has focused on “sexy” issues like the power spectrum and its cosmological implications, the primary stated science goal of WMAP is to produce maps. Indeed, one of the qualitatively new types of research made possible by WMAP involves taking advantage of this spatial information by cross-correlating the maps with other cosmological templates such as galaxy Peiris00 , x-ray Boughn02 ; DiegoSilk03 , infrared and lensing maps Song02 , which can reveal interesting signals ranging from the Late Integrated Sachs-Wolfe effect to the SZ effect and lensing Cooray01 ; Hirata02 .

Many such future studies will be looking for signals of modest statistical significance, so it is important to quantify and minimize unwanted signals in the map due to foreground contamination and detector noise. Accurately understanding the foreground signal is also important for the interpretation of the WMAP early reionization detection BennettWMAP ; Page3-03 ; Spergel03 ; Haiman03 ; Holder03 ; Hui03  and for the interpretation of the low WMAP quadrupole BennettWMAP ; hinshawpower03 ; Spergel03 , since Galactic foregrounds are most important on large angular scales wiener ; polarforgs . The WMAP team has already performed a careful foreground analysis BennettForegs  combining the five frequency bands into a single foreground-cleaned map, shown in Figure 1 (top). Given the huge effort that has gone into creating the spectacular multifrequency maps, it is clearly worthwhile to subject them to an independent foreground analysis. This is the purpose of the present paper, which we will argue not only corroborates the findings of the WMAP team, but also makes some further improvements that we believe are useful.

The main goal of this paper is to remove foregrounds, not to understand or model them. For reviews of foreground modeling issues, see, e.g., BennettForegs ; wiener ; foregpars ; polarforgs ; forgbook ; forgx ; Bouchet1 ; Bouchet2 ; Bouchet3 ; tucci00 ; Baccigalupi00 ; Burigana02 ; bruscoli02 ; tucci02 ; giardino02  and references therein. There is a rich literature of techniques for foreground removal. The work most closely related to the present paper is that done in preparation for the Planck mission Bouchet1 ; Bouchet2 ; Bouchet3 , developing multipole-based cleaning techniques and testing them on simulations.

## II ANALYSIS AND RESULTS

### II.1 The problem

A key goal of the CMB community is to measure the function $`x(\widehat{𝐫})`$, the true CMB sky temperature in the sky direction given by the unit vector $`\widehat{𝐫}`$. The WMAP team have observed the sky in five frequency bands denoted K, Ka, Q, V and W, centered on the frequencies of 22.8, 33.0, 40.7, 60.8 and 93.5 GHz, respectively, producing five corresponding maps (Figure 2, left) that we will refer to as $`𝐲_i`$, $`i=1,\mathrm{},5`$. In practice, each of these maps are discretized into $`n=12\times 512^2=3`$,$`145`$,$`728`$ HEALPix111The HEALPix package is available from http://www.eso.org/science/healpix/.  pixels healpix1 ; healpix2 , so $`𝐲_i`$ is an $`n`$-dimensional vector. However, since the maps are more than adequately oversampled relative to the beam resolution, it is equivalent and often simpler to think of them simply as five smooth functions $`y_i(\widehat{𝐫})`$. Conversely, we will occasionally write the true sky as a pixelized vector $`𝐱`$.

These maps are related to the true sky $`𝐱`$ by the affine relation

$$𝐲_i=𝐀_i𝐱+𝐧_i,$$
(1)
where the matrix $`𝐀_i`$ encodes the effect of beam smoothing and the additive term $`𝐧_i`$ is the contribution from detector noise and foreground contamination, collectively referred to as “junk” below since it complicates the recovery of $`𝐱`$. Defining the larger matrices and vectors

𝐀(
A1A5),𝐲(
y1y5),𝐧(
n1n5),formulae-sequence𝐀
A1A5formulae-sequence𝐲
y1y5𝐧
n1n5{\bf A}\equiv\left(\begin{tabular}[]{c}${\bf A}_{1}$\\
$\vdots$\\
${\bf A}_{5}$\end{tabular}\right),\quad{\bf y}\equiv\left(\begin{tabular}[]{c}${\bf y}_{1}$\\
$\vdots$\\
${\bf y}_{5}$\end{tabular}\right),\quad{\bf n}\equiv\left(\begin{tabular}[]{c}${\bf n}_{1}$\\
$\vdots$\\
${\bf n}_{5}$\end{tabular}\right),\quad (2)
we can rewrite the system of equations given by (1) as

$$𝐲=\mathrm{𝐀𝐱}+𝐧,$$
(3)
a set of linear equations that would be highly over-determined by a factor five if it were not for the presence of unknown junk $`𝐧`$.

Foreground removal involves inverting this overdetermined system of noisy linear equations. The most general linear222In addition to simplicity and transparency, linear methods have the advantage that the noise properties of $`\stackrel{~}{𝐱}`$ can be readily calculated from those of the input maps. estimate of $`\stackrel{~}{𝐱}`$ of the true sky $`𝐱`$ can be written

$$\stackrel{~}{𝐱}=\mathrm{𝐖𝐲}$$
(4)
for some $`n\times (5n)`$ matrix $`𝐖`$. We will require that the inversion leaves the true sky unaffected, i.e., that the expected measurement error $`\stackrel{~}{𝐱}𝐱`$ is independent of $`𝐱`$. Bennett et al. BennettForegs  refer to this property as the inversion having unit response to the CMB. Methods involving maximum-entropy reconstruction or some form of smoothing typically lack this property. Substituting equation (3) into equation (4) shows that this requirement corresponds to the constraint

$$\mathrm{𝐖𝐀}=𝐈.$$
(5)
### II.2 The mathematically optimal solution

Which choice of $`𝐖`$ gives the smallest rms errors from foregrounds and detector noise combined? Physically different but mathematically identical problems were solved in a CMB context by wright96 ; tegmark97 , showing that the minimum-variance choice is

$$\stackrel{~}{𝐱}=[𝐀^t𝐍^1𝐀]^1𝐀^t𝐍^1𝐲,$$
(6)
where $`𝐍\mathrm{𝐧𝐧}^t`$. For an extensive discussion of different methods proposed in the literature and the relations between them, see foregpars . Although optimal, this method is unfortunately unfeasible for the WMAP case, for two reasons. First, it requires the inversion of the $`(5n)\times (5n)`$ matrix $`𝐍`$. Although the detector noise contribution to this matrix is close to diagonal for WMAP, the foreground contribution is certainly not. Second, it requires knowing the matrix $`𝐍`$, which has many more components than there are pixels in the map.

### II.3 The WMAP team solution

In producing their internal linear combination (ILC) map, the WMAP team adopt a simpler approach BennettForegs , first smoothing all five maps to a common angular resolution of $`1^{}`$ and then performing the cleaning separately for each pixel. The smoothing implies that $`𝐀_i=𝐈`$ (if we redefine $`𝐱`$ to be the true sky smoothed to $`1^{}`$), and equation (4) reduces to the simple form

$$\stackrel{~}{𝐱}(\widehat{𝐫})=w_i𝐲_i(\widehat{𝐫})$$
(7)
for five weights $`w_i`$ that according to equation (5) must sum to unity. The WMAP team chose the weights that minimize the rms fluctuations in the cleaned map $`\stackrel{~}{𝐱}`$, using 12 separate weight vectors for 12 disjoint sky regions.

Although this method works well, it can be improved by allowing the weights to depend on angular scale (i.e. on harmonic number $`ł`$) as well as on Galactic latitude. This has two advantages. First, the angular resolution is limited not by that of the lowest resolution channel (the K-band FWHM is $`49.2^{}`$), but rather by that of the highest resolution channel (the W-band FWHM is $`12.6^{}`$). Second, as shown by wiener , letting the weights depend on angular scale produces a cleaner map by taking into account that the frequency dependence of the unwanted signals varies with scale: at large scales, galactic foregrounds dominate, whereas at small scales, detector noise dominates.

### II.4 Our solution

In this paper, we will take an approach intermediate between the two described above, aiming to approximate the optimal method while staying within the realm of numerical feasibility. Our method is essentially that of Tegmark & Efstathiou wiener , implemented to make it blind and free of assumptions both about the CMB power spectrum and about foreground and noise properties. The only assumption, which is crucial, is that the CMB has the Blackbody spectrum determined by the COBE/FIRAS experiment Mather94 . The method strictly respects the requirement of equation (5), which is most easily seen by verifying that each of the steps described below do so individually. The gist of our method is to combine the five WMAP bands with weights that depend both on angular scale and on distance to the Galactic plane (we subdivide the sky into 9 regions of decreasing overall cleanliness). We begin by describing the angular scale separation, then turn to the spatial subdivision in Section II.6.

We perform our cleaning in multipole space as in wiener , and go back and forth between spherical harmonic coefficients

$$a_\mathrm{}m^iY_\mathrm{}m(\widehat{𝐫})^{}𝐱_i(\widehat{𝐫})𝑑\mathrm{\Omega }$$
(8)
and real-space maps

$$x_i(\widehat{𝐫})=a_\mathrm{}m^iY_\mathrm{}m(\widehat{𝐫})$$
(9)
using the HEALPix package healpix1 ; healpix2  with $`\mathrm{}_{\mathrm{max}}=1024`$. Since this employs a spherical harmonic transform algorithm using fast Fourier transforms in the azimuthal direction Muciaccia97 , each transformation takes only about a minute on a Linux workstation. Each cleaned map is defined by

$$a_\mathrm{}m\underset{i=1}{\overset{5}{}}w_ł^i\frac{a_\mathrm{}m^i}{B_ł^i}$$
(10)
for some set of five-dimensional weight vectors $`𝐰_ł`$, where $`B_ł^i`$ is the beam function for the $`i^{th}`$ channel from page2-03 . (There are 4 W-band maps, 2 V-band maps and 2 Q-band maps; we combine these into single maps at each frequency by straight averaging and therefore average the corresponding beam functions as well.) When computing our final cleaned map in real space, we multiply $`a_\mathrm{}m^i`$ by $`B_ł^5`$ in equation (9) so that it has the beam corresponding to the highest-resolution WMAP band.

To gain intuition for the weight vectors $`𝐰_ł`$ that specify a cleaned map, we have plotted them for four interesting cases in figures 3, 4, 5 and 6. Figure 3 corresponds to the weighting used by the WMAP team for the region away from the galactic plane, and is simply independent of $`ł`$. To recover their published internal linear combination map shown in Figure 1 (top), one simply applies these weights after first multiplying each $`a_\mathrm{}m^i`$ by a Gaussian beam with FWHM=1 in equation (9).

The main drawback of this weighting is that it neglects that there is a tradeoff between foregrounds and detector noise which depends strongly on angular scale. Diffuse foregrounds are most important on large scales where detector noise is negligible, warranting large negative and positive weights to aggressively subtract foregrounds. Detector noise, on the other hand, is most important on small scales, both because of its Poissonian nature ($`𝐂_l`$ roughly constant in the observed map) and because the beam correction in equation (10) causes it to blow up exponentially on scales smaller than the angular resolution KnoxNoise ; wiener . In the limit $`ł\mathrm{}`$, the best weighting is therefore $`𝐰_ł(0,0,0,0,1)`$ regardless of what the foregrounds are doing, as illustrated in Figure 4, since the W-band has the best resolution.

We choose to minimize the total unwanted power from foregrounds and noise combined, separately for each harmonic $`ł`$ as in wiener , as opposed to only for the combination corresponding to the $`1^{}`$ pixel variance as in BennettForegs . As seen in Figure 5, one expects such a weighting to combine features from the two previous figures: rather aggressive subtraction using all five channels at low $`ł`$, more cautious subtraction using only the higher-resolution channels on intermediate scales, and all the weight on the W-band at extremely high $`ł`$. In particular, it is crucial to downweight the K-band when cleaning on the scales of the acoustic peaks, otherwise the acoustic peaks in the resulting map will be dominated by K-band noise.

The constraint equation (5) that we leave the CMB untouched corresponds to the requirement that the weights sum to unity ($`_i𝐰_ł^i=1`$) for each $`ł`$, i.e., that

$$𝐞𝐰_ł=1,$$
(11)
where $`𝐞=(1,1,1,1,1)`$ is a column vector of all ones. Minimizing the power $`|a_\mathrm{}m|^2`$ in the cleaned map of equation (10) subject to this constraint gives wiener ; foregrounds ; foregpars

$$𝐰_ł=\frac{𝐂_ł^1𝐞}{𝐞^t𝐂_ł^1𝐞},$$
(12)
where $`𝐂_ł`$ is the $`5\times 5`$ matrix-valued cross-power spectrum

$$𝐂_ł^{ij}a_{\mathrm{}m}^{i}{}_{}{}^{}a_\mathrm{}m^j.$$
(13)
As an example, Figure 6 shows the weights we obtain for the 3rd cleanest of the 9 sky regions shown in Figure 8 below. We see that just as forecast in Figure 5, and as in the WMAP team weighting of Figure 3, the 61 GHz V channel is “the breadwinner”, getting a large positive weight on large scales since it has the lowest overall foreground level. The 94 GHz W channel gets a negative weight to subtract out dust, and the three lower frequency channels are used to subtract out synchrotron, free-free and any other emission dominating at low frequencies. In cleaner sky regions, weights get less aggressive in the sense of acquiring smaller absolute values. In particular, we recover weights similar to those of the WMAP team (Figure 3) on large scales for the Kp2 sky cut defined and used by BennettForegs .

### II.5 Blind analysis and the power spectrum matrix

We compute the power spectrum matrix $`𝐂_l`$ in practice using the method of Hivon ; a similar approach was used by the Boomerang Ruhl02  and WMAP hinshawpower03  teams. This simply consists of expanding the masked sky patch in question in spherical harmonics and then correcting for window function effects. Our only variation is that we do not invert the window matrix to obtain anticorrelated band power estimates with delta function window functions. Rather, we simply divide by the quadratic estimators by the area of the window $`\delta 𝐓_ł^2`$ window functions, which asymptotes to the unbiased minimum variance estimators of cl  on scales much smaller than the sky patch analyzed. An example of our window functions is shown in Figure 7.

As an example, Figure 9 shows the measured power spectrum for the V-band, the cleanest of WMAP’s frequency bands.

One fact worth emphasizing is that our weighting scheme of equation (12) is totally blind, assuming nothing whatsoever about the CMB power spectrum, the foregrounds, the WMAP detector noise or external templates. The only assumption is that the CMB spectrum is the blackbody that the WMAP team have modeled it as, so that the CMB contributes equally to all five channels — otherwise the vector $`𝐞`$ would be replaced by some other constant vector. We see that there is no need to model the CMB, the foregrounds or the noise, since all we need for computing the optimal weights is the total power spectrum matrix $`𝐂_l`$, containing the total contribution from CMB, foregrounds and noise combined — and this can be measured directly from the data.

We can decompose $`𝐂_l`$ as a sum of two terms,

$$𝐂_ł=𝐂_ł^{\mathrm{junk}}+𝐂_ł^{\mathrm{cmb}}=𝐂_ł^{\mathrm{junk}}+C_ł^{\mathrm{cmb}}\mathrm{𝐞𝐞}^t,$$
(14)
where the second term is the CMB contribution and the first term is the contribution from noise and foregrounds. Note that if we keep $`𝐂_ł^{\mathrm{junk}}`$ fixed and change $`C_ł^{\mathrm{cmb}}`$, the weights given by equation (12) stay the same. The easiest way to see this is to note that the quantity we are minimizing is $`|a_\mathrm{}m|^2=𝐰^tC_ł𝐰=𝐰^t𝐂_ł^{\mathrm{junk}}𝐰+C_ł^{\mathrm{cmb}}(𝐞𝐰)^2=𝐰^t𝐂_ł^{\mathrm{junk}}𝐰+C_ł^{\mathrm{cmb}}`$, so the CMB power is just an additive constant that does not affect the optimal weighting. This means that our method is blind to assumptions about the underlying (ensemble-averaged) CMB power spectrum. Although we will return below in Section II.7 to the issue of how to determine what fraction of the power $`𝐂_ł`$ comes from each of the two terms in equation (14), it is important to remember that this affects only the physical interpretation, not our cleaning method and the maps we produce.

### II.6 Subdividing the sky

To minimize the variance in our cleaned map, we should take advantage of all ways in which the unwanted signals (noise and foregrounds) differ from the CMB in their contribution to the covariance matrix $`𝐍`$ in equation (6). Above we exploited their different dependence on angular scale $`ł`$. Unlike the CMB, foregrounds are not an isotropic Gaussian random field. Rather, their variance differs dramatically between clean and dirty regions of the sky. It is therefore desirable to subdivide the sky into a set of regions of increasing cleanliness and perform the cleaning separately for each one wiener . One then expects our method described above to settle on more aggressive weights for the dirtier regions, where foregrounds are much more of a concern than noise. A second advantage of such a subdivision is that the frequency dependence of the foregrounds is likely to differ between very dirty and very clean regions, again resulting in different optimal weights. The WMAP team used the latter argument to motivate their subdivision of the sky into 12 regions, and convincingly demonstrated that foreground spectra indeed did vary across the sky, notably for synchrotron radiation where the spectrum was found to steepen towards increasing galactic latitudes.

The WMAP team have created a set of sky-masks of increasing cleanliness based solely the K-band map. Although these masks are undoubtedly fine in practice, the procedure used in creating them is not blind, since it rests on the assumption that all dirty regions are dirty in the K-band. In particular, if a foreground manifests itself only at higher frequencies (imagine, say, a blob with localized dust emission without detectable synchrotron or free-free emission), it would go unnoticed. A second minor drawback of this K-band approach is that random CMB fluctuations affect the masks at a low level. In other words, the mask was based on the upper left map in Figure 2 which, as opposed to the upper right map, contains CMB fluctuations.

To preserve the blind nature of our method, we therefore create sky masks with a different procedure. We first form four difference maps W-V, V-Q, Q-K and K-Ka, thereby obtaining maps guaranteed to be free of CMB signal that pick up any signals with a non-CMB spectrum. We then form a combined “junk map” by taking the largest absolute value of these four maps at each point in the sky (Figure 8, top). Finally, we create disjoint sky regions based on contour plots of this map. We use cuts that are roughly equispaced on a logarithmic scale, corresponding to thresholds of 30000, 10000, 3000, 1000, 300 and 100$`\mu `$K (Figure 8, bottom). We emphasize that we have found no evidence whatsoever for any actual problems with the WMAP team masks, and opt to use our own simply to preserve the blind and CMB-independent nature of our analysis. As a cross-check, we also repeated our entire analysis using the WMAP masks Kp0, Lp2 and Kp12, obtaining similar results.

We followed the WMAP team procedure in the details of converting the contours into the masks shown in Figure 8 (bottom): we downsampled the junk map to HEALPpix resolution 64, imposed the cuts, went back up to HEALPix resolution 512, performed Gaussian smoothing with FWHM=$`2^{}`$ on the $`\{0,1\}`$-valued mask and imposed a cutoff of 0.5. The WMAP team reported strong spatial variations of foreground spectra in the innermost parts of the galactic plane, and therefore subdivided this into 11 disjoint regions. We were unable to reproduce this procedure since they did not specify which these regions were. Instead, we merely lopped off three spatially disconnected islands in the two dirtiest regions as their own separate masks, as illustrated in Figure 8 (bottom), leaving 9 separate masks in total.

Our multipole-based cleaning is nonlocal in the spirit of equation (4), and although it guarantees that the CMB signal is preserved separately for each pixel, this is of course not the case for foregrounds. To avoid mirages of foreground emission from the Galactic plane leaking up to high latitudes, we clean the galaxy “from inside out”, i.e., clean the dirtiest region first, the second dirtiest region second, etc. To be specific, we start by defining five temporary maps, initially set equal to the five WMAP channels (Figure 2, left). We then repeat the following cleaning procedure nine times, once corresponding to each region $`i=1,\mathrm{},9`$:

1. Compute the power spectrum matrix $`𝐂_ł`$ and the optimal weights $`𝐰_ł`$ for the $`i^{th}`$ region only.
2. Expand the five temporary maps in spherical harmonics and compute a cleaned all-sky map using the weights from step 1.
3. Replace the $`i^{th}`$ region of the temporary maps by the corresponding region in the cleaned map from step 2 (smoothed to the resolution appropriate for that channel, of course, by multiplying by $`B_ł^i`$ in $`ł`$-space).

Visually, one thus sees the foreground contamination gradually being cleaned off from Figure 2 as this iteration proceeds, starting in the inner Galactic plane and proceeding outward. At the end of the nine iterations, the five temporary maps equal the desired cleaned map at the K-, Ka-, Q-, V-, and W-band resolutions, respectively.

### II.7 Interpretation of the cleaned maps

Figure 1 (bottom) shows our final cleaned map and Figure 10 shows its power spectrum in sky regions of varying cleanliness. We plot all maps after 21 Gaussian smoothing (giving a net FWHM of 24) to prevent them from being undersampled by the pixels in the image. The reader interested in using this map can download the corresponding 13 Megabyte HEALPix fits file from the web333Our cleaned maps are available for download at http://www.hep.upen.edu/~max/wmap.html, together with a high-resolution version of this paper.. To use this map, it is important to be clear on how to interpret it.

First of all, it is a sum of CMB, foreground and detector noise fluctuations. Although it was constructed by minimizing its power spectrum $`C_ł^{\mathrm{clean}}=𝐰_ł^t𝐂_ł𝐰_ł`$, its power nonetheless gives a lower bound on the CMB power spectrum, since the minimization was performed subject to the constraint that the CMB be preserved; $`𝐞𝐰=1`$444 The only way in which its power could be biased low would be if random fluctuations in our estimate of the $`𝐂_ł`$ matrix conspired to remove power. Although we found no indication of this actually happening, we computed our weights $`𝐰_ł`$ using a heavily smoothed version of $`𝐂_ł`$-matrix as a precaution. Specifically, we smooth over at least $`\mathrm{\Delta }ł=10`$ or 100 $`(ł,m)`$-modes, whichever is larger, obtaining an $`ł`$-dependence of $`𝐂_ł`$ with no visible trace of random fluctuations..

So what fraction of the power seen in Figure 10 is due to CMB, foregrounds and detector noise, respectively? Let us first get some rough estimates from the figures, then present more quantitative limits. Since this paper is focused on minimizing foregrounds, not on physically modeling them, the interested reader is referred to the detailed foreground study of the WMAP team BennettForegs  for more information.

The WMAP team cleverly eliminated the average noise contribution by using only cross-correlations between different differencing assemblies (DAs) to measure the power spectrum, and used a combination of sky cutting and foreground subtraction (with external templates) to minimize the foreground contribution. The best fit cosmological model Spergel03 ; verde03  to their measured CMB power spectrum hinshawpower03  is shown for comparison in figures 9 and 10. For $`ł\mathrm{}<100`$, it is seen to agree well with the lower envelope of our curves in Figure 9, suggesting that foregrounds are subdominant in the cleanest parts of the 61 GHz sky. Figure 10 shows no noticeable excess power at $`łł\mathrm{}<100`$ due to foregrounds in the cleaned map in any of the four cleanest sky regions from Figure 8, which together cover all but the very innermost Galactic plane.

The slight power deficit on the very largest scales has two causes: one is the low quadrupole, to which we return in Section III below, which pulls down neighboring band power estimates because the band-power window, illustrated in Figure 7, contains small off-diagonal contributions. The second cause is that we have not corrected for the effects of monopole and dipole removal, which pulls down the power estimates on the scale of the sky patch in question (our method produces an unbiased CMB map regardless of what $`𝐂_ł`$ we use in equation (12), so this merely raises the variance in the cleaned map above the optimal level). On smaller scales, detector noise starts to dominate, and is seen to push our curves way above the CMB curve. A worthwhile future project for further quantifying foregrounds would be to repeat our analysis with $`𝐂_ł`$ estimated in a way that removes the detector noise contribution. For Q, V and W bands, which each have more than one DA, this can be done using cross-correlations. For the $`(1,1)`$ and $`(2,2)`$ components of $`𝐂_ł`$ (the K and Ka power spectra), it would involve subtracting the noise power using the WMAP team’s noise model.

We can, however, give some quantitative limits even based on our measured $`𝐂_ł`$ alone. Grouping the five coefficients $`a_\mathrm{}m^i`$ into a five-dimensional vector $`𝐚_{łm}`$, equation (13) becomes simply $`𝐂_ł=𝐚_{łm}^{}𝐚_{łm}^t`$, and the cleaned procedure can be written $`a_\mathrm{}m=𝐰𝐚_{łm}`$. Using equation (12) shows that the power in the cleaned map is

$$C_ł^{\mathrm{clean}}=\frac{1}{𝐞^t𝐂_ł^1𝐞}.$$
(15)
By subtracting our cleaned map from the input map, we obtain maps showing in Figure 2 (right). These are guaranteed to be free from CMB power, since the cleaned map gave a lower limit on the CMB. Most likely, we have subtracted some foreground power too, so the five maps should be interpreted as placing lower limits on the foreground power. (As a toy example, imagine synchrotron, free-free and dust emission tracing each other perfectly; the sum of their three spectra can then be written as a constant, which our method will interpret as CMB, plus a non-negative residual, which our method will interpret as foregrounds). Using equation (14), the covariance matrix of these five CMB-free “junk maps” is

$$𝐂_ł^{\mathrm{junk}}𝐂_ł^{\mathrm{nocmb}}𝐂_łC_ł^{\mathrm{clean}}\mathrm{𝐞𝐞}^t=𝚷_ł𝐂_ł𝚷_ł^t,$$
(16)
where the projection matrix

$$𝚷_ł𝐈\mathrm{𝐞𝐰}_ł^t$$
(17)
satisfies $`𝚷_ł^2=𝚷_ł`$, $`𝚷_ł𝐞=0`$ and $`𝚷_ł^t𝐰=0`$ and can be interpreted as projecting out the CMB component. In this notation, the maps in Figure 2 (right) are defined by simply $`𝚷_ł𝐚_{łm}`$. The inequality $`𝐂_ł^{\mathrm{junk}}𝐂_ł^{\mathrm{nocmb}}`$ refers only to the diagonal elements of these matrices. The five diagonal elements of $`𝐂_ł^{\mathrm{junk}}`$ are plotted in Figure 11 for our cleanest sky region from Figure 8.

How much of $`𝐂_ł`$ can possibly be due to CMB? In other words, forgetting for a moment about the weighting that gave equation (15), how large can we make $`C_ł^{\mathrm{cmb}}`$ in equation (16) before the covariance matrix $`𝐂_ł^{\mathrm{junk}}`$ gets unphysical properties? First of all, no covariance matrix can have negative eigenvalues, so we we must stop increasing $`C_ł^{\mathrm{cmb}}`$ once the smallest eigenvalue of $`𝐂_ł^{\mathrm{junk}}`$ drops to zero. In fact, $`𝚷^t𝐰=0`$ implies that when using our optimal weighting and hence equation (15), $`𝐂_ł^{\mathrm{junk}}`$ has a vanishing eigenvalue corresponding to the vector $`𝐰`$, and it is easy to show that this alternative method for estimating $`C_ł^{\mathrm{cmb}}`$ is equivalent to our original method, being simply an alternative derivation of equation (15).

Let us now make a second assumption: that $`𝐂_ł^{\mathrm{junk}}`$ cannot have any negative elements. The noise covariance matrix is guaranteed to have this property, since the absence of correlations between bands implies that it is diagonal. The foreground covariance matrix is also guaranteed to have this property if foreground emission is indeed emission, i.e., if foregrounds can make only positive contributions to the sky maps. Pure absorption at all frequencies likewise give only positive correlations. This assumption, made also in the Maximum-Entropy analysis of BennettForegs  is likely to be valid WMAP frequencies, since the only known exception is the thermal SZ effect, and it changes sign only outside of the WMAP frequency range (around 217 GHz). In summary, we therefore have two separate limits on the CMB power spectrum:

$`C_ł^{\mathrm{cmb}}`$ $``$ $`C_ł^{\mathrm{clean}},`$ (18)
$`C_ł^{\mathrm{cmb}}`$ $``$ $`\underset{i,j}{\mathrm{min}}(𝐂_ł)_{ij}.`$ (19)
$`C_ł^{\mathrm{clean}}`$ is also the actual variance of the cleaned map, so if the second limit is lower than first, then we know that the difference cannot be due to CMB.

To gain intuition, consider the following two examples of the junk + CMB decomposition of equation (14) for the simple case of only two frequency bands:

$`𝐂_ł=\left(\begin{array}{cc}\text{5}& \text{3}\\ \text{3}& \text{2}\end{array}\right)`$ $`=`$ $`\left(\begin{array}{cc}\text{4}& \text{2}\\ \text{2}& \text{1}\end{array}\right)+\left(\begin{array}{cc}\text{1}& \text{1}\\ \text{1}& \text{1}\end{array}\right)`$ (26)
$`𝐂_ł=\left(\begin{array}{cc}\text{2}& \text{1}\\ \text{1}& \text{2}\end{array}\right)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(\begin{array}{cc}\text{1}& \text{-1}\\ \text{-1}& \text{1}\end{array}\right)+{\displaystyle \frac{3}{2}}\left(\begin{array}{cc}\text{1}& \text{1}\\ \text{1}& \text{1}\end{array}\right)`$ (33)
$`=`$ $`\left(\begin{array}{cc}\text{1}& \text{0}\\ \text{0}& \text{1}\end{array}\right)+\left(\begin{array}{cc}\text{1}& \text{1}\\ \text{1}& \text{1}\end{array}\right)`$ (38)
In the first case, equation (15) gives $`C_ł^{\mathrm{cmb}}=1`$, i.e., the largest contribution we can possible attribute to the CMB (the second term) is $`\mathrm{𝐞𝐞}^t`$ — if we scaled up the last term, then $`𝐂_ł^{\mathrm{junk}}`$ (the first term) would acquire an unphysical negative eigenvalue. Moreover, we see that $`𝐂_ł^{\mathrm{junk}}`$ looks like the covariance matrix of a perfectly respectable foreground, with a spectrum such that it contributes twice as high rms fluctuations to the first channel as to the second, and with perfect correlation between the two (the dimensionless correlation coefficient is $`r=1`$). Such a foreground would be completely removed by our method, and so we cannot place any lower limit on the foreground contribution to the cleaned map in this case. In the second case, equation (15) gives $`C_ł^{\mathrm{cmb}}=3/2`$, so this will be the variance in the cleaned map. However, $`𝐂_ł^{\mathrm{junk}}`$ (the first term) in equation (33) exhibits an unphysical anticorrelation between the two bands which neither noise nor foregrounds could produce, and equation (19) tells us that the CMB power $`C_ł^{\mathrm{cmb}}1`$, corresponding to the alternative decomposition on the last line. This means that the junk map must have a power contribution of at least $`1/2`$ which is not due to CMB.

This difference, which places a lower limit on the amount of non-CMB signal in our cleaned map, is shown as the black line in Figure 11 for our actual five-band case. We see that although we get an interesting lower limit for $`ł>100`$, presumably dominated by detector noise, the residual foreground level is consistent with zero at low $`ł`$.

### II.8 Wiener filtering

Figure 1 (bottom) shows our cleaned map of the CMB sky, and this is the map that should be used for cross-correlation analysis with other data sets and other scientific applications. For visualization purposes, however, we can do better. The “best-guess” map of what the CMB looks like, in the sense of minimizing the rms errors, is the Wiener filtered map defined by Wiener49

$$a_\mathrm{}m^{\mathrm{Wiener}}=\frac{C_ł^{\mathrm{cmb}}}{C_ł^{\mathrm{clean}}}a_\mathrm{}m.$$
(39)
This staple signal processing technique, multiplying by “signal over signal plus noise”, has additional attractive properties; for instance, it constitutes the best-guess (maximum posterior probability) map in the approximation of Gaussian fluctuations. Examples of recent applications of Wiener filtering to CMB and galaxy mapping include Rybicki92 ; Lahav94 ; Fisher95 ; Zaroubi95 ; Bond95 ; Bunn96 ; saskmap . Our resulting Wiener filtered map is shown in Figure 1 (middle), using the best fit model from the WMAP team Spergel03 ; verde03  as our estimate of $`C_ł^{\mathrm{cmb}}`$ in the numerator of equation (39). For the denominator, we take the larger of our measured $`C_ł^{\mathrm{clean}}`$ and $`C_ł^{\mathrm{cmb}}`$, so that the ratio is guaranteed to be $`1`$. The result is not an unbiased CMB map. Rather, equation (39) shows that each multipole gets multiplied by a number between 0 and 1, so features with high signal-to-noise are left unaffected whereas features that are not statistically significant become suppressed. This means that the features that you see in Figure 1 (middle) are likely to be real CMB fluctuations, having signal-to-noise exceeding unity.

The cleaned map at the bottom of Figure 1 reveals some residual galactic fluctuations on very small angular scales, caused mainly by the fact that no other channels have fine enough angular resolution to help clean the W-band map for very large $`ł`$. We Wiener filter each of the nine sky regions from Figure 8 separately, so $`C_ł^{\mathrm{clean}}`$ is much higher near the Galactic plane. This is why the Galactic contamination is imperceptible in the Wiener filtered map: equation (39) automatically suppresses fluctuations in regions with large residual foregrounds.

## III Discussion

We have performed an independent foreground analysis of the WMAP maps to produce a cleaned CMB map. The only assumption underlying our method is that the CMB contributes equally to all five channels. This assumption rests on very solid ground BennettWMAP . The basic reason for this is that the COBE/FIRAS determination of the CMB spectrum Mather94  is based on the absolute CMB signal, which is about $`10^5`$ times larger than the fluctuations that we have considered in this paper.

Figure 1 shows that our map agrees very well with the internal linear combination (ILC) map from the WMAP team on the scales $`\mathrm{}>1^{}`$ where they can be compared. This is yet another testimony to the high signal-to-noise in the WMAP data and to the fact that unpolarized CMB foregrounds are manageable: the basic spatial features of the CMB are insensitive to the details of the foreground removal method used.

The basic advantage of our map is illustrated in Figure 12, which compares its all-sky power spectrum with that of the WMAP team ILC map. Both power spectra have had beam effects removed here; in Figure 1, our map is shown at the W-band resolution and the ILC map is shown at $`1^{}`$ resolution as released.

The first thing to note from Figure 12 is that foregrounds appear highly subdominant in both maps, since their power spectra essentially coincide with the WMAP CMB power spectrum on large scales where noise becomes unimportant. Second, as expected, the main improvement in our map is seen to be on the smaller scales where noise is important, gaining a factor of 30% at the first acoustic peak and about a factor of two at the second peak where noise from the low frequency channels is beginning to exponentially dominate the ILC map.

We hope that our map will prove useful for a variety of scientific applications. For cross-correlation with external maps, its lowered noise power should be particularly advantageous for pulling out small-scale signals, for instance those associated with lensing and SZ clusters. Rather that attempting detailed modeling of the residual noise and foreground fluctuations in our map, a simple way to place error bars on such correlations will be repeating the analysis with a suite of rotated and flipped versions of our map as in, e.g., 19ghz .

Table 1 – Measurements of the CMB quadrupole and octopole.

| Measurement | $`\delta T_2^2`$ $`[\mu `$K$`{}_{}{}^{2}]`$ | p-value | $`\delta T_3^2`$ $`[\mu `$K$`{}_{}{}^{2}]`$ |
| --- | --- | --- | --- |
| Spergel et al. model | 869.7 |  | 855.6 |
| Hinshaw et al. cut sky | 123.4 | 1.8% | 611.8 |
| ILC map all sky | 195.1 | 4.8% | 1053.4 |
| Cleaned map all sky | 201.6 | 5.1% | 866.1 |
| Cosmic quadrupole | 194.2 | 4.7% |  |
| Dynamic quadrupole | 3.6 |  |  |

Let us close by returning from small to large angular scales. The surprisingly small CMB quadrupole has intrigued the cosmology community ever since if was first observed by COBE/DMR Smoot92 , and simulations by Spergel03 ; verde03  show that the low value observed by WMAP is sufficiently unlikely to warrant serious concern. The WMAP team measured the quadrupole using only the part of the sky outside of their Galactic cut, and stress that the dominant uncertainty in its value is foreground modeling. While we agree with this assessment, it should be borne in mind that noise variance and beam issues are completely negligible on these huge angular scales, so this does not imply that the foreground uncertainties are large compared to the signal itself. Indeed, Figure 12 suggests that foregrounds are subdominant to the intrinsic CMB signal even without any Galaxy cut if a foreground-cleaned map is used. While the reader may feel disturbed by the clearly visible Galactic residuals in Figure 1 (in both the top and bottom maps), it is important to bear in mind that these signals are present in only a tiny fraction of the total sky area and therefore contribute little to the total power spectrum.

To quantify the impact of the Galactic plane, we computed a sequence of all-sky power spectra for both our map and the ILC map. We first used the unadulterated maps, then repeated the calculation after repacing the dirtiest of the nine regions from Figure 8 with zeroes, then zeroed our the 2nd dirtiest region as well, etc. This zeroing procedure obviously biases the measurements by removing power, by an amount related to the zeroed area, but provides a powerful test of how sensitive the results are to Galactic plane details. We also band-pass filtered the resulting maps to produce spatial plots of the quadrupole, octopole and hexadacapole. We found that zeroing out the dirtiest parts of the Galactic plane had a negligible effect on both the power spectrum and on the spatial structure of these lowest multipoles. Specifically, we could zero out all but the three cleanest regions Figure 8 (everything with $`T>1`$mK) without the quadrupole or octopole changing substantially. The spatial morphology of the quadrupole, octopole and hexadecapole for the all-sky analysis agrees well between our map and the ILC map, again showing insensitivity to galaxy modeling details (in particular, the ILC map is likely to have less contamination in the Galactic plane due to more subdivisions there).

Although more detailed foreground modeling would be needed to rigorously quantify the foreground contribution to low multipoles, let us, encouraged by the above-mentioned tests, tentatively assume that this contribution is unimportant and perform an all-sky analysis of our cleaned map. The resulting power spectrum is shown in Figure 13, which is simply a blow-up of the leftmost part of the previous figure, and the corresponding angular correlation function is shown in Figure 14.

Table 1 summarizes the quadrupole and octopole results. We see that although the quadrupole is still low, it is not quite as low as that from the cut-sky WMAP team analysis of hinshawpower03 . Moreover, our map has a quadrupole virtually identical to the WMAP team ILC map despite the differences in foreground modeling, further supporting our hypothesis that the quadrupole is not strongly affected by foregrounds. The second column in Table 1 shows the probability of the quadrupole in our Hubble volume being as low as observed if the best-fit WMAP team model from Spergel03 ; verde03  is correct. This is computed for the all-sky case where $`\delta T_2^2`$ has a $`\chi ^2`$-distribution with 5 degrees of freedom, so the probability tabulated is simply $`1\gamma [5/2,(5/2)T_2^2/855.6\mu \mathrm{K}^2]/\mathrm{\Gamma }(5/2)`$, where $`\gamma `$ and $`\mathrm{\Gamma }`$ are the incomplete and complete Gamma functions, respectively555 The actual probability is slightly larger for the cut-sky case where there are fewer effective degrees of freedom, although Spergel03  show with Monte Carlo simulations that there it is still a disturbingly small for the measured WMAP quadrupole — indeed, they obtain a consistency probability of $`0.7\%`$ using Markov chains.. We see that the statistical significance of the low quadrupole problem drops substantially with our all-sky analysis, below the 95% significance level, requiring merely a one-in-twenty fluke.

To understand why the all-sky quadrupole is larger than the cut-sky one measured by the WMAP team, we plot the lowest three multipoles of our cleaned map in Figure 15, all on the same temperature scale. Several features are noteworthy. First of all, the quadrupole is low with an interesting alignment. A generic quadrupole has three orthogonal pairs of extrema (two maxima, two minima and two saddle points). We see that the actual CMB quadrupole has its strongest pair of lobes, apparently coincidentally, fall near the Galactic plane. Applying a Galaxy cut therefore removes a substantial fraction of the quadrupole power. The saddle point is seen to be close to zero. In other words, there is a preferred axis in space along which the observed quadrupole has almost no power.

Second, the observed quadrupole is the sum of the cosmic quadrupole and the dynamic quadrupole due to our motion relative to the CMB rest frame - see KamionKnox02  for a detailed discussion. Since this motion is accurately known from the CMB dipole measurement Smoot92 ; BennettWMAP , the dynamic quadrupole can and should be subtracted when studying the cosmic contribution. Figure 15 shows the dynamic quadrupole approximated by $`(v/c)^2(\mathrm{cos}\theta ^21/3)`$, where $`v369km/s`$ towards $`(l,b)(264,48)`$ is the velocity of the Solar System relative to the CMB BennettWMAP  and $`\theta `$ is the angle relative to this velocity vector. This is a small correction with peak-to-peak amplitude $`(v/c)^24\mu \mathrm{K}`$, and Table 1 shows that it reduces the cosmic quadrupole slightly. This approximation is crude since the dynamic quadrupole in fact depends on frequency KamionKnox02 , and therefore may have been either over- or under-estimated in our foreground cleaned map \- we have shown it here merely to illustrate its spatial orientation and give an order-of-magnitude estimate of its importance.

Third, although the overall octopole power is large, not suppressed like the quadrupole, it too displays the unusual property of a preferred axis along which power is suppressed. Moreover, this axis is seen to be approximately aligned with that for the quadrupole. The reason that our measured octopole in Figure 13 is larger than that reported by the WMAP team is therefore, once again, that much of the power falls within the Galaxy cut. In contrast, the hexadecapole is seen to exhibit the more generic behavior we expect of an isotropic random field, with no obvious preferred axis.

How significant is this quadrupole-octopole alignment? As a simple definition of preferred axis for an arbitrary sky map, smalluniverse  computes the unit vector $`\widehat{𝐫}`$ around which the angular momentum dispersion $`(\widehat{𝐫}Ł)^2=_mm^2|a_{łm}(\widehat{𝐫})|^2`$ is maximized. Here $`a_{łm}(\widehat{𝐫})`$ denotes the spherical harmonic coefficients of the map in a rotated coordinate system with its $`z`$-axis in the the $`\widehat{𝐫}`$-direction. The preferred axes $`\widehat{𝐫}_2`$ and $`\widehat{𝐫}_3`$ for the quadrupole and octopole, respectively, are

$`\widehat{𝐫}_2`$ $`=`$ $`(0.1145,0.5265,0.8424),`$
$`\widehat{𝐫}_3`$ $`=`$ $`(0.2578,0.4207,0.8698),`$ (40)
i.e., both roughly in the direction of $`(l,b)(110^{},60^{})`$ in Virgo. The angle between these two axes is merely $`10.3^{}`$, and the probability that a random axis falls inside a circle of radius $`10.3^{}`$ around the quadrupole axis is simply the area of this circle over the area of the half-sphere, about $`1/62`$. In other words, if the CMB is an isotropic Gaussian random field, then a chance alignment this good requires about a 1-in-62 fluke. This issue is discussed in greater detail in smalluniverse .

What does this all mean? Although we have presented these low multipole results merely in an exploratory spirit, and more thorough modeling of the foreground contribution to $`ł=2`$ and $`ł=3`$ is certainly warranted, it is difficult not to be intrigued by the similarities of Figure 13 with what is expected in some non-standard models, for instance ones involving a flat “small Universe” with a compact topology as in Starobinsky ; Stevens ; bb1 ; bb2 ; Levin01 ; Rocha  and one of the three dimensions being relatively small (of order the Horizon size or smaller). This could have the effect of suppressing the large-scale power in this particular spatial direction in the same sense as is seen in Figure 13. As so often in science when measurements are improved, WMAP has answered old questions and raised new ones.

We thank the the WMAP team for producing such a superb data set and for promptly making it public via the Legacy Archive for Microwave Background Data Analysis (LAMBDA)666The WMAP data are available from http://lambda.gsfc.nasa.gov. . Support for LAMBDA is provided by the NASA Office of Space Science. We thank Krzysztof Górski and collaborators for creating the HEALPix package healpix1 ; healpix2 , which we used both for harmonic transforms and map plotting. Thanks to Ed Bertschinger, Gary Hinshaw, Jim Peebles, Ned Wright and Matias Zaldarriaga for helpful comments. This work was supported by NSF grants AST-0071213, AST-0134999, & AST-0205981, and by NASA grants NAG5-10763 & NAG5-11099. MT was supported by a David and Lucile Packard Foundation fellowship and a Cottrell Scholarship from Research Corporation.