text:
# X-ray Scattering Halos from the Galactic Center: Implications for Diffuse Emission Around Sgr A^∗

## 1 Introduction

Our Galaxy’s central supermassive black hole provides an unprecedented opportunity for studying the accretion physics of these objects, which are thought to power active galactic nuclei. Since we must view the Galactic Center (GC) through a large column of gas and dust, observations are restricted to wavelengths in and longer than the near infrared (NIR) and to X-rays harder than about 2 keV.

Near infrared observations reveal a dense stellar cluster (Genzel et al. 2003a, and references therein), whose inner members have high proper motions, consistent with the presence of a black hole of mass $`3\times 10^6M_{}`$ (Schödel et al. 2002; Ghez et al. 2003a). Its location, as derived from stellar orbits, coincides with the compact, nonthermal, weakly variable, linearly and circularly polarized radio source Sgr A$`^{}`$(Balick & Brown 1974; Bower et al. 2003, and references therein) and a flaring X-ray point source discovered by the Chandra X-ray Observatory (Baganoff et al. 2001; 2003, hereafter B01 and B03, respectively). The X-ray emission from this source is extended on scales comparable to the Bondi radius, $`0.04`$ pc. Recently a flaring NIR source has been detected at this position (Genzel et al. 2003b; Ghez et al. 2003c).

The spatial and spectral structure of the X-ray emission may help to constrain the diverse set of theoretical models that have been proposed for the accretion onto the black hole (e.g., Falcke, & Markoff 2000; Melia, Liu, & Coker 2000; Liu & Melia 2002; Quataert 2002; Yuan, Markoff, & Falcke 2002; Yuan, Quataert, & Narayan 2003; Quataert 2003). It has been argued that the low luminosity of Sgr A$`^{}`$with respect to its expected classical Bondi accretion luminosity is evidence for very inefficient accretion mechanisms (e.g. Narayan 2002, and references therein). Models that attempt to explain the arcsecond scale diffuse emission via accretion onto a clustered population of unresolved neutron stars have also been proposed (Pessah & Melia 2003). Nayakshin & Sunyaev (2003) proposed that the X-ray flares may result from the interaction of stars with a cold gas disk orbiting close to Sgr A$`^{}`$.

In this paper we consider the scattering and absorption of X-rays from the Galactic Center. In particular we emphasize the importance of scattering in producing extended emission — “halos” — on arcsecond scales. Since much of the intervening dust is located close to the GC, these halos are much more concentrated than those around sources observed along more typical Galactic lines of sight. We assess the contribution made by the concentrated scattering halos of Galactic center sources to the diffuse emission seen within $`10\mathrm{}`$ of Sgr A$`^{}`$. The halos are so concentrated that they can dominate the intensity of the Chandra point spread function (PSF) inside 1″, so they affect determinations of the point source flux and thus intrinsic luminosity. They also affect the estimates of the intensity of emission due to diffuse gas close to the Bondi radius, and thus estimates of the gas density and the black hole’s Bondi accretion rate.

We consider the distribution of gas and dust towards the Galactic center in §2. We then use X-ray observations and a model for how dust and gas absorb and scatter X-rays to calculate the intrinsic X-ray luminosity and spectrum of Sgr A$`^{}`$3). This involves a correction for dust scattering that requires a calculation of the angular intensity profiles of scattering halos as a function of energy. In §4 we calculate the total intensity of the scattering halo, including examples of the delayed halos from flares, and compare to observations. We discuss the implications of our results and conclude in §5.

## 2 The Gas and Dust Distribution to the Galactic Center

To calculate X-ray scattering we require a model for the spatial distribution of dust along the midplane of the Galaxy from the Sun to the Galactic center. The total column of gas and dust also affects the amount of absorption. We take the distance from the Sun to the Galactic center to be 8 kpc (Reid 1993). We assume a uniform metallicity and dust-to-gas ratio, taken to be 1.5 times the value for the local ISM, consistent with Galactic abundance gradients measured in H II regions (Giveon, Morisset, & Sternberg 2002).

### 2.1 Spatial Distribution

Based on studies of the radial distribution of the atomic and molecular gas (e.g. Binney & Merrifield 1998, Fig. 9.19) we approximate the distribution of gas from the main Galactic disk via a Gaussian located at $`R=4`$ kpc and with a HWHM of 2 kpc. We denote the fraction of the total column that is in this outer component by $`f_o`$. This approximation underestimates the local ISM density close to the Sun, but this material makes little contribution to inner X-ray scattering halos, on scales of arcseconds.

A second component of the total column is the material in the inner few hundred parsecs of the Galaxy, and we denote the fraction of material in this region by $`f_i`$. The properties of this gas have been constrained by data from IRAS (Infrared Astronomical Satellite) and COBE (Cosmic Background Explorer) observations (Launhardt, Zylka, & Mezger 2002). A gas disk with typical densities in its diffuse component of $`n_\mathrm{H}300\mathrm{cm}^3`$ is inferred to extend out to $`R200300\mathrm{pc}`$. An inner, denser, inclined, circumnuclear disk is seen as a separate dynamical component (Güsten et al. 1987; Marr, Wright, & Backer 1993) extending from $`R210\mathrm{pc}`$. However, most of this material does not intersect our line of sight to Sgr A$`^{}`$. The distribution of gas in these regions is very inhomogeneous: within the $`200`$ parsec disk are giant molecular clouds, such as the “50 km/s cloud” and the “20 km/s cloud”, which lie close to the line of sight to Sgr A$`^{}`$. Vollmer, Zylka, & Duschl (2003) estimate that the former is $`05\mathrm{pc}`$ closer to us than Sgr A$`^{}`$, while the latter is spread out over a larger distance with its densest region $`2550\mathrm{pc}`$ closer than Sgr A$`^{}`$. The centroids of these clouds are within $`510\mathrm{pc}`$ of our Sgr A$`^{}`$sight line. There is also evidence for clumps of hot molecular gas in the circumnuclear disk on scales $`10\mathrm{pc}`$ (McGary, Coil, & Ho 2001), and even inside the central parsec (Herrnstein & Ho 2002), so that the total gas column to a particular location in the central region may show significant variation. In the study of Herrnstein & Ho (2002) structure was seen on scales all the way down to the beam size of 10″$`0.4\mathrm{pc}`$.

As a first attempt to describe this complicated density structure, we shall approximate the radial distribution of the inner component with a Gaussian centered at various distances, $`R_i`$, towards us from Sgr A$`^{}`$, and with a HWHM equal to $`R_i`$. As a fiducial case we take $`R_i=100`$ pc, but we also consider models with $`R_i=50`$ and 25 pc. In all cases we assume a small central cavity that is free of dust. We take the radius of this cavity to be 2 pc, and then adjust the normalizations of the Gaussians so that the total column from 2 pc to 8 kpc is equal to the adopted overall totals. Individual gas clumps and clouds may have a column density comparable with the total through the more diffuse regions (below). To allow for the possibility of such an enhancement we consider models with a range of column densities in the smoothly distributed component.111Note that the scattering halos due to material inside the inner few parsecs are typically smaller than the size of the point spread function (PSF) of Chandra, and thus do not significantly affect the observable halo profile.

### 2.2 Estimates of the Total Column Density

The total columns of gas and dust to the GC are uncertain. One of the most common methods used to measure these quantities is by estimation of the K band extinction by simultaneously fitting near-IR photometry of GC stars with spectra of a particular stellar type and with a certain amount of reddening (Becklin et al. 1978; Lebofsky, Rieke, & Tokunaga 1982; Sellgren et al. 1987; Rieke, Rieke, & Paul 1989, hereafter RRP). This method requires knowledge of the extinction law in the near-IR (e.g. Rieke & Lebofsky 1985) and accurate models for the intrinsic stellar spectra. The main discriminant for the stellar spectra are the strengths of the CO absorption bands: they become stronger for decreasing temperature and increasing luminosity. $`\mathrm{H}_2\mathrm{O}`$ absorption at the edges of the K band can help distinguish between M giants and K supergiants, because the strength of this absorption decreases with increasing luminosity. Typically, solar metallicity stars have been used as standards for late-type spectra (e.g. Kleinmann & Hall 1986). One possible source of systematic error is the effect of elevated metallicities, which would tend to increase the strength of the CO and $`\mathrm{H}_2\mathrm{O}`$ features for a given luminosity and surface temperature.

An alternative to fitting detailed stellar spectra is to assume that certain, particularly blue, sources have Rayleigh-Jeans spectra: this method was used by RRP for the sources making up the IRS16 complex, which have subsequently been identified as hot, emission line stars (e.g. Eckart et al. 1995; Krabbe et al. 1995; these observations also show that IRS13, discussed below, is a complex dominated by hot He $`\mathrm{I}`$ stars), however the derived dust columns are still somewhat sensitive to the choice of extinction law (see below).

A third method to measure extinction is from the ratio of Br$`\alpha `$ to Br$`\gamma `$ recombination lines. The difficulties here are (1) the need to assume a particular theoretical ratio for the intrinsic line ratio, e.g. based on case B recombination theory at a particular density and temperature, and (2) the measurement of the line ratios in what is a very complicated and crowded region of the sky.

The relatively bright infra-red sources closest to Sgr A$`^{}`$are the IRS16 complex, IRS13, and IRS7, which have projected separations from Sgr A$`^{}`$of $`0.1`$, 0.13, and 0.20 pc. The values of $`A_K`$ adopted by RRP for these sources were 3.4, 3.5, 4.1, respectively, where the value for IRS16 is an average of the 7 components reported by RRP. The extinction law used by RRP has $`A_\lambda \lambda ^\beta `$ between $`1\mu m`$ and $`3\mu m`$, with $`\beta =1.61`$. The Milky Way dust models of Weingartner & Draine (2001, hereafter WD01) with $`R_VA_V/E(BV)=3.1,4.0,5.5`$ have $`\beta 1.73,1.25,1.25`$. If the extinction at K is derived from the degree of reddening of a Rayleigh-Jeans spectrum, then the values of $`A_K`$ derived from WD01 models are $`0.93,1.29,1.29`$ times that derived by RRP. Similar corrections may be expected for determinations of $`A_K`$ from fitting more realistic stellar spectra. Draine (2003a) has discussed the extinction per H nucleon in regions of different $`R_V`$. Using $`A_I/N_\mathrm{H}`$ from Draine (2003a) and the grain size distributions of WD01, we find $`A_K/N_\mathrm{H}=6.04,7.18,8.19\times 10^{23}\mathrm{mag}\mathrm{cm}^2`$. In the fiducial, “local Milky Way metallicity” case of the WD01 model with $`R_V=3.1`$ the effective H column densities for the three GC sources described above are thus $`5.2,5.4,6.3\times 10^{22}\mathrm{cm}^2`$. For the models with $`R_V=4.0,5.1`$ the column densities are different by factors of 1.09 and 0.92, respectively. If the effective metallicity is greater than solar by a certain factor, then these estimates would be reduced by the same factor, assuming a constant dust to metals ratio and for fixed grain properties.

Stars such as S0-2, which are very close to Sgr A$`^{}`$in both the plane of the sky and in line of sight dimension, are much fainter than the brighter IR sources discussed above, so the continuum slopes of their published spectra (Ghez et al. 2003b) are not yet reliably determined (A. Ghez 2003, private communication). With this caveat in mind, the K band extinction can be estimated assuming an intrinsic Rayleigh-Jeans spectrum (which is a reasonable assumption) and assuming a K band continuum spectrum $`F_\lambda \lambda ^\gamma `$ with $`\gamma =5.06.5`$. These values yield $`A_K4.65.4`$ for the RRP normalization, about 50% greater than the value estimated using the same technique towards IRS16.

The above range of values for the extinction appear to be fairly typical of the GC region on larger scales. Using the same reddening law as RRP, Figer et al. (1998) find $`A_K=3.2\pm 0.5`$ towards the Quintuplet star cluster, located 30 pc in projection from Sgr A$`^{}`$. Lebofsky & Rieke (1987) found the extinction over a 24 square arcminute region around the GC was spatially variable, but never less than $`A_K3.4`$. Sellgren et al. (1987) (using $`\beta =1.9`$ so that values of $`A_K`$ are about a factor of 0.85 of those derived using $`\beta =1.61`$ of RRP) reported that their extinction results were consistent with a uniform component with $`A_K2`$ and a variable component due to material near the GC with $`A_K01`$. These results are consistent with the recent model of Launhardt et al. (2002), who find $`A_K1.7`$ ($`A_V15`$) due to dust in the outer Galactic disk and a similar value for extinction due to dust in the GC region itself.

We note that the clumps reported by Herrnstein & Ho (2002) are estimated (but with large uncertainties) to have $`N_\mathrm{H}2\times 10^{24}\mathrm{cm}^2`$, which would correspond to extinctions a factor of $``$30 larger than typically reported, if they contain typical dust grain populations.

As we shall describe later, estimates of the total column density have also been made by fitting models for the X-ray emission from Sgr A$`^{}`$to observations. Porquet et al. (2003) found $`N_\mathrm{H}1.9\times 10^{23}\mathrm{cm}^2`$, corresponding to $`A_K119.9`$ for $`R_V=3.35.5`$. B01 found $`N_\mathrm{H}=5.3_{1.1}^{+0.9}\times 10^{22}\mathrm{cm}^2`$, i.e. about a factor of 3 or 4 smaller than Porquet et al. (2003).

Given the uncertainties in the estimates of the total gas and dust columns along the line of sight to the GC, and since it is clear that there is some variation towards different objects in this region, we shall consider a range of different density distributions and total column densities. As an average for the gas along the line of sight to the GC, we adopt a metallicity of 1.5 times solar and a dust-to-gas mass ratio of 1.5 times the local interstellar value, which is taken to be 0.008. This metallicity enhancement is consistent with ISO (Infrared Space Observatory) observations of Galactic H II regions (Giveon, Morisset, & Sternberg 2002). We assume depletions into dust grains of 67% for C, 20% for O, 90% for Mg, Si, and Fe, and 0% for N, Ne, and S (WD01). For the fiducial model we set $`A_K=3.5`$, which for the WD01 grain model with $`R_V=3.1`$ and the above metallicities yields $`N_{\mathrm{H},\mathrm{tot}}=3.90\times 10^{22}\mathrm{cm}^2`$. For this model we also set $`f_i=f_o=0.5`$. We shall also consider a model with $`A_K=7.0`$, so that $`N_{\mathrm{H},\mathrm{tot}}`$ is twice the above value, and which has $`f_i=0.75`$ and $`f_o=0.25`$: i.e. we have added material only to the inner component.

The $`R_V=3.1`$ dust model with $`\overline{Z}=1.5Z_0`$ has C/H = 90 ppm in polycyclic aromatic hydrocarbons (PAHs) and consists of a mixture of carbonaceous grains and amorphous silicate grains. As discussed by Li & Draine (2001), the carbonaceous grains have the properties of PAH molecules when they contain $`10^4`$ C atoms, and the properties of graphite particles when they contain $`10^5`$ C atoms. By altering the size distributions of the carbonaceous and silicate particles, the dust model appears to be able to reproduce observed extinction curves in various Galactic regions, and in the Large and Small Magellanic Clouds. This dust model, when illuminated by starlight, produces infrared emission consistent with the observed emission spectrum of the interstellar medium (Li & Draine 2001, 2002).

For the amorphous silicate grains and the carbonaceous grains we employed the dielectric functions recently estimated by Draine (2003b) and calculated the scattering using Mie theory. The X-ray absorption and scattering properties of this model are consistent with observations towards Nova Cygni 1992 that probe typical dust in the local diffuse ISM (Draine & Tan 2003). In Figure 1 we show the energy dependence of the optical depth (due to photoelectric absorption by gas and extinction by dust) of our models of gas and dust toward the GC. The contribution of the dust extinction due to scattering is also shown.

Lutz et al. (1996) reported extinction towards the GC in the 4-8 $`\mu \mathrm{m}`$ region of the spectrum exceeding that expected for standard graphite-silicate mixes (e.g. Draine 1989), which may indicate that additional species are present, such as iron/metal sulphides or metal oxides, or carriers related to the production of the 3 $`\mu \mathrm{m}`$ ‘ice’ absorption features that are seen in the spectra. From an analysis of such absorption features, Chiar et al. (2000) found evidence for enhanced molecular material along the line of sight to Sgr A$`^{}`$as compared to the Quintuplet star cluster.

A shift in the grain size distribution toward larger grains, with increased values of $`R_V`$, is expected in the dense, cold conditions of the cores of giant molecular clouds. Since these conditions may apply to a significant amount of the gas and dust towards Sgr A$`^{}`$, we shall also consider the effect of using models with $`R_V=5.5`$ for the entire dust column. Dust models with larger grains yield lower values of $`N_{\mathrm{H},\mathrm{tot}}`$ for a given $`A_K`$ and metallicity. For example, the $`R_V=5.5`$ model of WD01 gives columns that are a factor of 0.70 smaller than the $`R_V=3.1`$ case.

Given an observed X-ray spectrum we can use the models shown in Figure 1 to estimate the form of the intrinsic spectrum (§3). However, we shall also see that when the flux from a particular source is defined as all photons arriving from within a certain angular radius (e.g. on arcsecond scales for Chandra), then a substantial fraction of the photons that are scattered can remain in this source region and must be accounted for.

## 3 The Intrinsic X-Ray Spectrum of Sgr A$`^{}`$

Utilizing the high spatial resolution ($`0.5\mathrm{}`$) of Chandra, B01 and B03 reported the first detection of an X-ray source that can be convincingly associated with Sgr A$`^{}`$. B03 fit the observed spectrum (derived from counts within 1.5$`\mathrm{}`$ of the source from the first 40.3 ks observation) with an absorbed power-law model, with $`(d\dot{N}/dE)_0=A_0(E/\mathrm{keV})^\mathrm{\Gamma }\mathrm{ph}\mathrm{cm}^2\mathrm{s}^1\mathrm{keV}^1`$ with $`\mathrm{\Gamma }=2.7_{0.9}^{+1.3}`$ and $`N_\mathrm{H}=9.8_{3.0}^{+4.4}\times 10^{22}\mathrm{cm}^2`$ (90% confidence interval). However, these calculations did not account for dust scattering (B01), leading to an overestimation of the column density necessary to produce the observed soft X-ray cutoff. B01 reported results from a second observation with effective exposure time of 35.4 ks, which included a $`10`$ ks flare. From a joint analysis of the quiescent emission from the two epochs and the flare emission from the second epoch, B01 derived $`\mathrm{\Gamma }_q=2.2_{0.7}^{+0.5}`$, $`\mathrm{\Gamma }_f=1.3_{0.6}^{+0.5}`$, and $`N_\mathrm{H}=5.3_{1.1}^{+0.9}\times 10^{22}\mathrm{cm}^2`$, where the subscripts $`q`$ and $`f`$ refer to the quiescent and flaring states, respectively. Models of the spectra of an absorbed optically thin thermal plasma (Raymond & Smith 1977) yielded similar column densities with $`kT=1.9_{0.5}^{+0.9}\mathrm{keV}`$ (B03). Systematic errors due to the energy dependence of the enclosed energy fraction inside 1.5$`\mathrm{}`$ and of charge transfer inefficiency effects could lead to an overestimation of $`\mathrm{\Gamma }`$ by $`0.20.3`$, an underestimation of the intrinsic X-ray luminosity by $`20\%`$, and an overestimation of $`N_\mathrm{H}`$ (B03).

XMM-Newton observed a bright 3 ks flare from Sgr A$`^{}`$(Porquet et al. 2003) and from these observations were derived $`\mathrm{\Gamma }_f=2.5\pm 0.3`$, and $`N_\mathrm{H}=20\pm 3\times 10^{22}\mathrm{cm}^2`$. The spectrum could also be fit with models based on bremsstrahlung continuum, thermal blackbody, and optically thin plasma emission.

The method of estimating the intrinsic spectrum of a source typically involves counting all the photons arriving inside a given circular area around a source (a radius of $`1.5\mathrm{}`$ was used for the Chandra observations of the quiescent state), and then subtracting off a background flux, estimated from a surrounding annulus. Given an assumed dust and gas column, a correction is then made to estimate the intrinsic flux that would be received if there was no absorption. However, for gas and dust distributed close to the source, much of the scattering halo is located very close to the PSF, and needs to be accounted for.

Consider a source observed at 4 keV. For the model with $`A_K=3.5`$ that has $`f_o=f_i=0.5`$, we have $`\tau _{\mathrm{gas}}=0.19`$, $`\tau _{\mathrm{dust},\mathrm{ext}}=0.45`$, and $`\tau _{\mathrm{dust},\mathrm{sca}}=0.32`$. Figure 2 shows the radial intensity profiles of dust-scattered halos of a monochromatic source at 4 keV for the adopted models for the density distribution to the GC, i.e. for $`A_K=3.5`$ and $`7.0`$. The halo profile was convolved with the on-axis energy-dependent Chandra PSF, derived from the CIAO (Chandra Interactive Analysis of Observations) software. We see that the scattering halos can dominate over the PSF inside 1″. In this figure we also show the dependence on the scattering halo on the grain size distribution, as parameterized via $`R_V`$: a model with a larger value of $`R_V`$, corresponding to larger average grain sizes, produces a halo that is more concentrated inside a few arcseconds. The halos at several different energies are compared in Figure 3: at lower energies the halo intensities are stronger relative to the PSF.

The fraction of enclosed flux in the scattering halo, relative to that in the PSF is shown in Figure 4. These results allow one to correct for the contribution that a scattering halo makes to the direct flux from a source, given the angular scale inside which the direct flux of a particular observation is estimated, i.e. the observed flux inside an angle $`\theta `$ is the sum of the unscattered flux and the part of the scattered flux that lies inside $`\theta `$.

For absorption and scattering with our fiducial density models, we calculate the intrinsic spectra of Sgr A$`^{}`$in its quiescent and flaring states, based on the observed fluxes and appropriate extraction radii (Figure 5 and Table 1). From these results we see that the uncertainty in the intrinsic X-ray luminosity of Sgr A$`^{}`$is set by the uncertainty in the extinction, as parameterized by $`A_K`$. This also sets the uncertainty in the slope of the spectrum, e.g. when fit by power-law models. The quiescent state 2-10 keV luminosity is a few$`\times 10^{33}\mathrm{ergs}\mathrm{s}^1`$. The two flares reported by B01 and Porquet et al. (2003) are about a factor of 10 to 30 times more luminous in this energy range, when averaged over the flare durations. We estimate that the peak luminosity of the flares is a factor of 1.3 and 2 times greater than the mean values quoted in the table. We tend to favor a harder intrinsic spectrum for the flare than Porquet et al., who find their best fitting power-law model has $`\mathrm{\Gamma }=2.5\pm 0.3`$ (and $`N_\mathrm{H}=2\times 10^{23}\mathrm{cm}^2`$), however at least part of this discrepancy is due to the fact that our model uses a much smaller column of gas: $`3.9,7.8\times 10^{22}\mathrm{cm}^2`$ for the $`A_K=3.5,7.0`$ cases respectively. Part of this difference is due to our use of abundances that are a factor of 1.5 times greater than the local interstellar values, and part may be due to differences in the grain models, including our preference to normalize to the observed extinction in the K band, rather than the canonical visual extinction of 30 mag. Our results are broadly consistent with those of B01.

## 4 The X-ray Scattering Halo of Sgr A$`^{}`$

X-ray scattering by dust creates diffuse “halos” of emission around astronomical sources (Overbeck 1965; Martin 1970; Hayakawa 1973), and these can be used to constrain models of the dust (e.g. Catura 1983; Mathis, & Lee 1991; Predehl & Klose 1996; Smith, Edgar, & Schafer 2002; Draine & Tan 2003, and references therein). X-ray scattering is most efficient at small angles, $`1`$ degree, and if the dust and gas have a fairly uniform distribution along the line of sight to a source, then this angle is also representative of the size of the observed scattering halo. The diffuse emission directly associated with Sgr A$`^{}`$extends only over a few arcseconds (B03). It has been argued that any diffuse emission created by dust-scattering would be fairly uniform on these scales, and therefore removed during the process of background subtraction: B03 estimated the background from the mean intensity between scales 3.5″to 10″.

However, dust that is close to the source tends to produce a more intense and concentrated scattering halo. This is relevant to the case of Sgr A$`^{}`$because much of the line-of-sight dust is thought to be within a 100 pc or so of the Galactic center. The median scattering angle, $`\theta _{s,50}360\mathrm{}(\mathrm{keV}/E)`$, and dust at a distance $`R`$ from a point source produces a scattered halo with 50% of the scattered power within $`\theta _{h,50}(R/D)\theta _{s,50}2.25\mathrm{}(R/50\mathrm{p}\mathrm{c})(8\mathrm{k}\mathrm{p}\mathrm{c}/D)(\mathrm{keV}/E)`$, where $`D`$ is the distance from the observer to the point source, and multiple scattering is neglected (Draine 2003).

For the adopted models of the source spectrum in the quiescent and flaring states, Figure 6 shows the calculated azimuthally-averaged intensity profile of the scattering halo. As in the monochromatic case, the halo profile was convolved with the on-axis Chandra PSF, derived from the CIAO data analysis software.222Smith et al. (2002) have shown that the PSF at angles greater than $`10\mathrm{}`$ is best derived empirically from study of astronomical sources. In this study we restrict our attention to scales inside $`10\mathrm{}`$, where the PSF is well-determined from theoretical modeling of the telescope’s optics. Note that the energy dependence of the PSF can be seen in Figure 3.

The strength of the scattering halo is relatively more important for the quiescent state, because of its softer spectrum. However, in all cases the halo intensity can become significant relative to the PSF on scales of about an arcsecond or so.

In Figure 7 we show the effect of changing the location of the inner component of dust that is near Sgr A$`^{}`$. As the dust is brought closer to Sgr A$`^{}`$, the halo becomes more concentrated, becoming more and more point-like, and thus indistinguishable from the intrinsic source after convolution with the (Chandra) PSF.

In Figure 8 we compare the theoretical radial intensity profiles, which are the sum of the PSF and the dust-scattered halo, with the observed emission from and around Sgr A$`^{}`$in its quiescent state (B03). The observations are presented as the mean values of annuli 0.492″wide, which is the pixel size. Note that the dithering of the telescope during the observation means that these profiles are not broadened significantly by the finite pixel size. B03 also corrected for pointing errors, estimated to be $`0.15\mathrm{}`$ by $`0.37\mathrm{}`$. We also show the radial intensity profiles of two other nearby point sources, J174540.9-290014 and J174538.0-290022, observed by B03, normalized to the intensity of the innermost annulus. These sources have comparable spectra to Sgr A$`^{}`$, are separated on the sky from Sgr A$`^{}`$by 18.7″(0.73 pc) and 27.0″(1.05 pc), respectively. Being heavily absorbed they are likely to be in the GC region, and given the angular distribution of sources around Sgr A$`^{}`$, their proximity makes it likely that they are physically close to Sgr A$`^{}`$in all three spatial dimensions. It is likely, although not certain (§2), that their lines of sight are subject to similar extinctions, particularly for contributions from gas and dust that is $`R25\mathrm{pc}`$ towards us from the Galactic center.

We find that dust-scattering of point source emission is important for broadening the PSF of Chandra for typical sources in the vicinity of the GC. The level of enhancement is sufficient to explain the observed profiles of the sources near Sgr A$`^{}`$. The largest enhancement at the arcsecond scale occurs with the inner dust component about 50 pc from Sgr A$`^{}`$, and with dust that has $`R_V=5.5`$, i.e. larger grains. However, even with the model with the largest dust column, corresponding to $`A_K=7.0`$ magnitudes, the scattering halo appears to be insufficient to explain the intensity of the diffuse emission within a few arcseconds of Sgr A$`^{}`$by factors of several. The total dust column would have to be increased by similar factors above the value of the $`A_K=7.0`$ model in order to explain the entire diffuse emission as being due to dust-scattered emission from an unresolved, steady point source.

We note that our model intensity profiles are greater than the data in the central annulus by about a factor of 1.8. This is because we have normalized the models to produce the observed flux, which was estimated inside an extraction radius of $`1.5\mathrm{}`$. The models underpredict the intensity beyond about 0.6″, and so to compensate have larger central intensities. If there are instrumental effects that lead to broadening, such as pointing errors, that have not been corrected for in the data analysis of B03, then these would cause the intensity of the central annulus to appear artificially smaller. If there are no such effects (as suggested by the profiles of the two nearby point sources) so that the observed intensity of the central annulus is accurate, then our point source models would have to be reduced in flux by a factor of about 1.8. This would lead to a similar reduction in the halo intensities and the contribution they make to the extended emission.

We now make a crude estimate for the contribution of the scattering halos from sources surrounding Sgr A$`^{}`$to the diffuse background inside 10″. From Figure 4 we estimate that at 4 keV for the $`A_K=7.0`$ model the scattering halo contributes about 20% of the PSF flux at angles 1.5-10″. B03 report that the sources inside 10″have a combined flux of $`1.4\pm 0.1\times 10^2\mathrm{counts}\mathrm{s}^1`$. Thus over the 40.9 ks of the observation, the number of expected counts per pixel (0.492″) is about 0.1. The $`A_K=3.5`$ model would predict about half of this value. These estimates are much smaller than the 1.2 counts per pixel measured by B03 for the strength of the diffuse background, so we conclude that the scattering halos from the innermost GC sources make a relatively minor contribution to the local diffuse emission.

### 4.1 Time Delay of Scattered Photons

An additional effect that can boost the observed intensity is the delayed scattering of enhanced flux from an earlier flare. B03 observed that the flux at the start of their observation was rapidly declining and three times greater than the mean value. They concluded that a flare occurred around the time of the start of their observation ($`99.34\%`$ confidence). Since flares from Sgr A$`^{}`$have been observed to occur very frequently (about once per day, Baganoff 2003), here we estimate the effects of such flares on the scattered halos. The time dependence of the scattered halo intensity at various angles can help to constrain the dust and gas distribution near the Galactic center. If a typical strong flare occurred just before or at the start of the reported observation of quiescent emission of B03, it could have contributed significantly to the intensity of the diffuse emission.

Photons that reach us via one or more scatterings will have to travel a longer path length and so will be delayed compared to those arriving directly. If the source is variable, such as a flare, then this time delay can yield information on the distance and dust distribution to a source (Trümper & Schönfelder 1973; Draine & Tan 2003). For single scattering from dust that is a distance $`R=yD`$ along the line of sight to a source at distance $`D`$, the delay is

$`\delta t_1(y,\theta )`$ $`=`$ $`{\displaystyle \frac{D}{c}}\left\{{\displaystyle \frac{y}{\mathrm{cos}\theta }}+(1y)\left[1+\left({\displaystyle \frac{y}{1y}}\right)^2\mathrm{tan}^2\theta \right]^{1/2}1\right\}`$ (1)
$``$ $`{\displaystyle \frac{D}{c}}{\displaystyle \frac{y}{1y}}{\displaystyle \frac{\theta ^2}{2}}=9.7\left({\displaystyle \frac{D}{8.0\mathrm{kpc}}}\right){\displaystyle \frac{y}{1y}}\left({\displaystyle \frac{\theta }{1\mathrm{}}}\right)^2\mathrm{s}.`$ (2)
We expect much of the dust towards Sgr A$`^{}`$to be close to the source: for $`R_i=100,50,25`$ pc, $`1y=0.0125,6.25\times 10^3,3.125\times 10^3`$, so the delay $`\delta t_1`$ for $`\theta =1\mathrm{}`$ becomes 760, 1500, 3100 s.

Figure 9 shows the time evolution of the scattered halo intensity following two different hypothetical flares, whose properties we have take to be similar to the flare observed by B01 (F1) and that observed by Porquet et al. (2003) (F2). The evolution of intensities depends on the distribution of gas near Sgr A$`^{}`$. Given the actual lightcurve of an observed flare, and with careful subtraction of the quiescent component, one may use the time-varying component of the X-ray halo to constrain the line of sight gas distribution close to the Galactic Center.

To compare to the observations of B03 we find the mean intensity averaged over the $``$40 ks observation period of B03. The results are shown in Figure 10. The delayed arrival of X-ray photons from when Sgr A$`^{}`$was in a more luminous state can boost the mean halo intensity in the following $``$40 ks by a factor of two or so. We stress that for these particular observations, while there is evidence for a flare near the beginning of the observation period, the actual lightcurve is unconstrained. Thus our calculation makes use of representative flare properties, as seen during other observations.

## 5 Conclusions

The large column densities of gas and dust towards the Galactic center strongly absorb and scatter X-rays, so that a careful treatment of these processes is necessary to infer the intrinsic properties of sources in this region, such as their luminosity, spectra, and spatial size. The fact that much of the dust lies close to the GC means that the scattering halos are much more concentrated than the halos seen around typical sources in the outer Galaxy. The same phenomenon also occurs for halos observed around X-ray sources in star-forming regions, since these also tend to contain a large column density of dust near the source.

In this paper we have analyzed the effects of the intervening gas and dust on Galactic center X-ray sources, using a realistic range of models for the total column and spatial distribution. We used reasonable estimates for chemical abundances and depletions of the gas and dust. We calculated the radiative transfer using realistic scattering properties of dust grains, and considered the effects of multiple scattering. We also evaluated the importance of different grain size distributions.

The principal conclusions of our work are the following:

1. The intrinsic luminosities and spectra of emission from Sgr A$`^{}`$in its quiescent state and in two flaring states are shown in Figure 5 and Table 1. These estimates depend on the amount of intervening dust and gas, which we have related to the most relevant observable, the extinction in the K band. The quiescent 2-10 keV luminosity ranges from 3 to $`6\times 10^{33}\mathrm{ergs}\mathrm{s}^1`$, as $`A_K`$ ranges from 3.5 to 7.0 magnitudes. At the same time, the spectral index, $`\mathrm{\Gamma }`$, ranges from 1.8 to 3.3. The flaring states are about 10 to 30 times more luminous, with the flare observed by Chandra having a harder spectral index by about 0.7-1.0 and the flare observed by XMM-Newton having a harder spectral index by about 2.0. The uncertainties can be reduced once a better estimate for $`A_K`$ is available.

2. We show that scattering by dust within $`100`$ pc of the Galactic center can account for the arcsecond-scale broadening of Galactic center X-ray sources seen in Chandra observations by B03.

3. We have calculated the azimuthally averaged intensity profile of the dust-scattered halo around Sgr A$`^{}`$, and shown that it can dominate the Chandra PSF beyond $``$1 arcsec (see Fig. 7). This effect contributes to the diffuse emission seen around Sgr A$`^{}`$.

4. The dust-scattered halos of multiple sources in the GC region overlap, but only account for at most $`10\%`$ of the diffuse emission inside 10″of Sgr A$`^{}`$.

5. Motivated by the detection of enhanced emission at the start of the observation of B03 and by other observations that show that strong flares occur about once per day (Baganoff 2003), we have modeled the effects of flared emission from Sgr A$`^{}`$. For reasonable distributions of the dust, the time delay for scattered emission is of order a thousand to several thousand seconds at 1″, increasing quadratically with the angle. We show how measurements of the evolution of the scattered halo following a flare can help constrain the dust’s spatial distribution. It is possible that delayed scattering from a flare accounts for some of the diffuse emission close to Sgr A$`^{}`$in the observation of B03. To disentangle the quiescent component, data in a period at least $`50`$ ks after a major flare needs to be analyzed.

6. Our fiducial models for the contribution of a dust-scattered halo to the extended emission of Sgr A$`^{}`$can explain up to $`1/3`$ of the observed intensity at $`1\mathrm{}`$, corresponding to the extent of the Bondi radius. Thus earlier estimates of the gas luminosity and density may be overestimated by factors of $`3/2`$ and $`\sqrt{3/2}`$, respectively. Previous determinations of the Bondi accretion rate and luminosity could be overestimated by similar factors. We conclude, like B03, that the emission around Sgr A$`^{}`$is extended. To explain the entire emission as being due to dust scattering would require at least a doubling of the observed column in material at least $``$25 pc from Sgr A$`^{}`$(but not substantially changing the columns to the other nearby GC sources), and/or the presence of an exceptionally strong flare just prior to the observation period. The latter possibility can be tested by looking at the observed halo from Sgr A$`^{}`$in more quiescent periods in existing (but unpublished) data.

We thank Fred Baganoff for providing the radial intensity profiles of several GC sources in electronic form and for helpful explanations about the Chandra observations and data reduction. We thank Diab Jerius, Eliot Quataert and Randall Smith for helpful discussions. This work was supported in part by NSF grant AST-9988126, and in part by NASA grant NAG5-10811. JCT has received support via a Spitzer-Cotsen Fellowship from the Department of Astrophysical Sciences and the Society of Fellows in the Liberal Arts of Princeton University.


format: text

text:
# Resonance Effects in the Interaction of NLS Solitons with Potential Wells

## Abstract

The interaction of nonlinear Schrödinger (NLS) solitons with potential wells with variable shapes is investigated numerically. For fixed initial velocities below the threshold for transmission, the outcome pattern as a function of the width of the potential yields periodically repeating regions of trapping, transmission and reflection. The observed effects are explained by an excitation and a following resonant deexcitation of amplitude (shape) oscillations of the solitons at the boundaries of the well, associated with radiation modes.

Self-localized nonlinear waves (solitons) have been studied in many areas of physics including optics, solid state, molecular, plasma, elementary particles etc. Within integrable models, solitons exhibit remarkable stability - they propagate with constant velocities and shapes and emerge from collisions unchanged except for phase and space shifts. Real physical systems are often described by nonintegrable equations or such containing nonintegrable perturbations. This leads to inelastic soliton interactions with a variety of outcomes. As solitons provide an important mechanism for energy and information transport in nonlinear systems, such interactions have attracted considerable attention (see i.e. Ref. 1 for a review of earlier works on soliton dynamics in nearly integrable systems). Investigations have been focussed on collisions between solitons in nonintegrable models and interactions of solitons with defects and inhomogeneities. In both cases, due to the inelasticity of the interactions, solitons can change their velocities, break into a number of localized and dispersive waves and/or be trapped into bound states. In addition, interesting resonance phenomena have been observed.

Resonace effects in kink-antikink collisions have been studied numerically in some nonintegrable equations including $`\varphi ^4`$, double and modified sine-Gordon and others abl ; camp1 ; peyr . For initial velocities below the threshold for trapping, a sequence of narrow regions of reflection have been obtained. These reflection windows have been explained by a ”two-bounce” resonance mechanism involving excitation of an internal shape mode during the first collision, temporal trapping of the solitons due to loss of kinetic energy, deexcitation of the shape mode during the second (backward) collision and escape of the kinks to infinity (reflection). The resonance condition requires that the time between the two collisions is commensurate with the period of the shape mode. Fine three- and four-bounce resonance structures have also been obtained camp2 . Collisions of vector NLS solitons have been investigated in Yang1  where fractal resonant patterns have been obtained.

Similar effects have been observed in the interactions of solitons with impurities. The latter break the translational symmetry of the unperturbed system and create an effective potential for scattering or capture of the solitons kiv1 . Resonance effects in the kink-impurity interaction have been investigated in kiv2 ; fei1 ; fei2 . It has been shown in particular that kinks can be reflected by an attractive impurity via a ”two-bounce” resonance mechanism, analogous to that of kink-antikink interaction involving the excitation and deexcitation of a localized impurity mode fei1 , or an impurity and a shape mode fei2 . Scattering of NLS solitons from point defects has been studied in kiv3 ; kiv4 ; cao  involving a variety of nonresonant outcomes.

A problem of considerable theoretical and practical importance is the interaction of solitons with extended defects shar ; ting ; frau ; kalb . Trapping of solitons in potential wells and nonclassical behavior for kinetic energies close to the height of potential step have been obtained in kalb , but no resonance phenomena have been observed. The investigations however have been restricted to potential widths comparable to the soliton width. In the present work we investigate in detail the dynamics of bright NLS solitons impinging on potential wells with variable shapes. For fixed initial velocities slightly below the threshold for transmission, the increase of the width of the well yields alternating regions of capture and transmission, and occasionally - narrow reflection windows. The regions of transmission, capture and reflection follow a remarkable periodicity. The observed effects are explained by excitation and a following resonant deexcitation of amplitude (shape) oscilations of the soliton at the boundaries of the well. These oscillations are not associated with true internal shape modes of the soliton, but with dispersive radiation modes.

As a model in our numerical simulations we used the discrete nonlinear Schrödinger equation which describes the dynamics of nonlinear Bose-type excitations in atomic or molecular chains in the presence of defects which change locally the energy:

$`i{\displaystyle \frac{\alpha _n}{t}}=(\alpha _{n+1}+\alpha _{n1}2\alpha _n)2\alpha _n^2\alpha _n+d_n\alpha _n.`$ (1)
In the continuum limit corresponding to wide solitons (compared to the lattice constant) it turns into a perturbed NLS equation:

$$i\frac{\alpha }{t}+\frac{^2\alpha }{x^2}+2\alpha ^2\alpha =d(x)\alpha .$$
(2)
For $`d(x)0`$ (2) possesses a fundamental bright soliton solution:

$`\alpha (x,t)={\displaystyle \frac{1}{L}}\mathrm{sech}{\displaystyle \frac{xvtx_0}{L}}e^{i\left(vx/2\omega t\right)},\omega ={\displaystyle \frac{v^2}{4}}{\displaystyle \frac{1}{L^2}},`$ (3)
where $`L`$ and $`v`$ are the width and the velocity of the soliton.

It is well known, that for $`d_n=0`$ (1) is a nonintegrable discretization of the completely integrable continuum NLS equation. For sufficiently wide solitons however, the discreteness-induced effects are negligible, and the solution (3) is stable on ideal discrete lattices. We checked numerically the stability of (3) with $`L4`$ for very long time intervals. Thus (3) was input as initial condition in the simulations, placed 50 sites away from the defect region to avoid radiation losses due to overlapping. A predictor-corrector method sham  was used, periodic boundary conditions and chains much longer than the defect region in order to eliminate boundary effects. The accuracy of the calculations was controlled through the conservation of the norm (number of particles), which is better than $`10^6`$.

The total energy associated with the solution (3) on an ideal lattice is:

$$E_s=_{\mathrm{}}^{\mathrm{}}(\frac{\alpha }{x}^2\alpha ^4)𝑑x=\frac{v^2}{2L}\frac{2}{3L^3},$$
(4)
where the first term describes the kinetic energy of the free quasiparticles, and the second term - the nonlinear energy responsible for the formation of the soliton. The evolutionary pattern depends in general on the interplay between these two energies and the energy of interaction with the defects $`E_d`$

$$E_d=_{\mathrm{}}^{\mathrm{}}d(x)\alpha ^2𝑑x$$
(5)
The effects studied below correspond to the most interesting case of slow solitons with kinetic energy of the order of the energy of interaction with the defect and large nonlinear energy:

$$E_{kin}|E_d||E_{int}|$$

or equivalently:

$$v^2L^22|d|L1.$$
(6)
The large nonlinear energy is required to preserve the integrity of the soliton during the interaction.

Scattering of bright NLS solitons from single point defects has been studied in detail in kiv3 ; kiv4 ; cao . The corresponding interaction energy when the soliton is on top of the defect is:

$$E_d=d_{\mathrm{}}^{\mathrm{}}\delta (x)\alpha ^2𝑑x=\frac{d}{L^2}$$
(7)
When $`E_{kin}|E_d|`$, the solitons are not influenced considerably by the defect, and for $`E_{kin}|E_d|`$ the solitons are reflected even by an attractive defect. The possible outcomes in the case of slow solitons ($`E_{kin}|E_d|`$) and moderate defect strengths are transmission or capture. No resonance reflection windows have been obtained. A natural question arises as to what happens when the defect spreads over several lattice sites. The energy of interaction with $`N`$ consecutive defects when the soliton is in the middle of the defect region is:

$$E_d=d_{N/2}^{N/2}\alpha ^2𝑑x=\frac{2d}{L}\mathrm{tanh}\frac{N}{2L}$$
(8)
One can expect, that for a small number of defects, ($`NL`$), the evolution should be similar to this of a soliton interacting with a single defect with $`N`$-times greater strength. This turns out to be true only for kinetic energies much smaller than the interaction energy. In the case of comparable energies, the delocalization of the defect (which decreases the interaction energy) can change the evolution from capture to transmission (Fig. 1). It is worth noting that the sharper the defect - the stronger the radiation accompanying the interaction.

The focus of the present study lies in the interaction of NLS solitons with potential wells considerably wider than the soliton width. We modelled rectangular potential wells by $`N`$ equal consecutive defects with $`d=0.007`$ and used solitons in the form (3) with widths $`L=5.75`$, which are stable on an ideal lattice. For initial velocities $`v<0.04`$ the solitons are trapped inside the well and for $`v>0.06`$ they pass through it and escape to infinity. For initial velocities in the intermediate region, the outcome pattern as a function of the width of the well exhibits periodically alternating regions of transmission and capture, and occasionally, at the boundaries between them - narrow reflection windows. This is shown on Fig. 2 where we have plotted the final velocity of the soliton as a function of the width of the well for different values of the initial velocity. When trapped, the soliton oscillates back and forth inside the well with zero average velocity which we have plotted as final. The sharp minima with negative final velocity ($`v_f<0`$) correspond to reflection. These reflection ”windows” are extremely sensitive to the initial velocity. The relative widths of the regions of transmission and capture depend on the initial velocity and can be quite different, but the period of repeat remains constant and depends weakly on the velocity.

Fig. 3 illustrates the evolutionary patterns corresponding to capture, transmission and reflection. It is clearly seen, that amplitude (shape) oscillations are excited when the soliton enters the potential well and persist while the soliton is inside the well. Whenever the soliton leaves the defect region, the shape oscillations are almost totally extinguished. This brings forward the explanation of the observed resonance phenomena: when the soliton reaches the defect region it interacts inelastically with the sharp boundary and looses part of its kinetic energy exciting additional modes. These can be internal shape modes of the soliton and/or dispersive modes (radiation). We can distinguish between them by examining the frequency of the shape oscillations. The frequencies of the true internal shape modes can be obtained by adding small perturbation $`\alpha _1(x,t)=\phi _1(x)e^{i\mathrm{\Omega }t}`$ to the unperturb solution (3) and solving the linearized Schrödinger equation for it:

$$\frac{^2\phi _1}{x^2}+(\mathrm{\Omega }+\frac{4}{L^2\mathrm{cosh}^2(x/L)})\phi _1=0.$$
(9)
The frequencies of the shape mode determined from (9) are $`\mathrm{\Omega }_1=0.0735`$ and $`\mathrm{\Omega }_2=0.00094`$. The period of the amplitude oscillations obtained from the numerical data on Fig. 3 is T=208 which yields a frequency of $`\mathrm{\Omega }=2\pi /T=0.030`$. It practically coincides with the internal frequency of the unperturbed soliton (3) $`\omega =0.0296`$. Shape oscillations with the internal soliton frequency have been reported in Chbat  and are explained by interference of the soliton with radiation modes with twice the soliton frequency. This shows that similarly to the case of collision of vector solitons Yang1 , the amplitude oscillations which we observe are due to radiation, emitted during the inelastic interaction of the soliton with the boundary, and not to true internal shape modes. The soliton crosses the defect region accompanied by the radiation and this configuration is stable and weakly damped. When the soliton meets the second boundary, different outcomes are possible depending on the timing. In the general case, as the initial velocities are very small, the reduced kinetic energy of the soliton is not sufficient to overcome the potential barrier of the second boundary, the soliton is reflected from it and eventually gets trapped. However, the interaction of the soliton with the boundary is phase-sensitive and if the time for which the soliton crosses the defect region is commensurate with the period of the shape oscillations, the inelastic interaction with the second boundary can extinguish the shape oscillations, adding their energy back to the kinetic energy of the translational motion and allowing the soliton to overcome the barrier and escape to infinity which leads to transmission. The higher the initial velocity of the soliton - the wider the transmission regions (Fig. 2). In some cases the resonant condition for escape is achieved after the soliton has crossed the defect region twice - in the forward and backward directions. This yields the observed narrow reflection windows which are analogous to the three-bounce resonances observed in camp2 . Due to the radiation losses accompanying the propagation of the soliton inside the well, these higher-order resonances are very sharp, extremely sensitive to the initial velocity and difficult to observe.

The interaction energy when the soliton is at the boundary of the defect region is:

$$E_d=d_0^N\alpha ^2𝑑x=\frac{d}{L}\mathrm{tanh}\frac{N}{L}$$
(10)
For narrow potential wells, comparable to the soliton width, the interaction energy is smaller, leading to weaker shape oscillations and a narrower region of trapping. For sufficiently wide potential wells ($`NL`$, $`\mathrm{tanh}(N/L)1`$), when the soliton does not ”feel” the second boundary, the interaction energy and the excited shape oscillations are constant, which results in the observed periodic outcome pattern. For very wide potential wells, due to the radiation losses accompanying the oscillating soliton, the regions of transmission get narrower and eventually close down.

An increase of the depth of the well leads to wider regions of trapping and narrower regions of transmission (Fig. 4), while the total period remains unchanged. The inelastic interaction of the soliton with the boundary is stronger in this case and a larger portion of the kinetic energy of the soliton is transformed into radiation. A more exact resonance condition is required at the second boundary for the escape of the soliton, which yields narrower regions of transmission.

Contrarily, the change of the shape of the potential well from rectangular to trapezoid leads to wider transmission regions and narrower regions of capture (Fig. 5). The potential in this case is smoother and the perturbation it induces is weaker. Hence a smaller portion of the kinetic energy is transformed into radiation and the resonance condition at the second boundary is more relaxed.

We also checked the dependence of the evolutionary pattern on the initial position of the soliton with respect to the boundary of the defect region. For an initial soliton in the form (3) and a fixed velocity $`v=0.05`$ we obtain a threshold initial distance of 15 lattice sites, above which the soliton passes through the defect region and below it gets trapped. This can be explained by the radiation losses when the initial unperturbed soliton overlaps with the defect region. The chaotic behavior of the outcome with the initial soliton position obtained in kalb  can be attributed to the different type of coupling between the soliton and the shape mode.

###### Acknowledgements.

This work is supported in part by the National Science Foundation of Bulgaria under Grant No. F911.


text:
# Cross-correlation of the CMB and foregrounds phases derived from the WMAP data

## 1 Introduction

The recently-published Wilkinson Microwave Anisotropy Probe (WMAP) data sets (see Bennett et al. 2003 a-c, Hinshaw et al. 2003 a-b) open a new epoch for the CMB investigation. These data serve for the development of more refined technique for future high resolution measurements and for a choice of the realistic cosmological model.

The WMAP has observed the full sky in five frequency bands: K (centered frequency 22.8 GHz), Ka (33.0 GHz),Q (40.7 GHz), V (60.8 GHz), and W (93.5 GHz) and produced five maps represented in $`n=12\times 512^2`$ HEALPix (Górski et al. 1999) pixels. Then, using smoothing of the maps to $`1^o`$ and performing cleaning of the combined multifrequency map by minimization of the rms variance for each pixel, the WMAP team produced the Internal Linear Combination (ILC) map taking 12 optimization coefficients $`w_i,i=1,2,\mathrm{},5`$ for 12 disjoint sky regions. Thus, the ILC map has a minimal variance corresponding to the minimization of the Galaxy and foreground contamination at the range of multipoles $`\mathrm{}100`$.

Tegmark, de Oliveira-Costa and Hamilton (2003) (hereafter TOH) have suggested to use ideologically similar but technically different optimization scheme based on Tegmark and Efstathiou (1996) (hereafter TE) method. The main goal of the TOH method is to combine the K–W band maps into the Foreground Cleaned Map (FCM) using variable weighting coefficients $`w_i(\mathrm{})`$ for each frequency band $`i`$=K–W. In order to minimize the total unweighted power from the Galaxy, foregrounds and noise separately for each harmonic $`\mathrm{}`$, TOH subdivided the sky into 9 disjoint regions and performed cleaning of the 5 frequency maps in each region. Unlike the ILC map, the TOH FCM has smaller power at the multipole range $`\mathrm{}12`$. Both the ILC and TOH FCM are available from the web and can be used for analysis of the foregrounds and the statistical properties of the cleaned maps. For example, Komatsu et al. (2003) tested the non–Gaussianity of the WMAP CMB signal using Minkowski functionals as a statistic. Chiang et al. (2003) have also tested the statistics of the phases for the ILC and the TOH FCM maps and discussed non-Gaussianity of the TOH FCM map. Dineen and Coles (2003) have presented a diagnostic of the Galactic synchrotron contamination and have discovered the cross–correlation between Faraday rotation measure and the TOH FCM map (see also Coles et al. 2003).

Recently Naselsky et al. (2003) proposed the Phase Cleaning Method (PCM) for decomposition of the CMB signal and foregrounds using the WMAP data. The main idea of the PCM is to minimize the variance of the combined K–W maps and cross–correlations of the CMB and foreground phases simultaneously. The result of PCM in its application to the WMAP data reproduces well the power spectrum of the best-fit WMAP $`\mathrm{\Lambda }`$CDM model at the range $`\mathrm{}<50`$. The main target of our paper is to extend the resolution of the derived PCM CMB map up to $`\mathrm{}100`$, to perform the the cross-correlation analysis (circular and linear statistics and “friend–of–friend” analysis) of the CMB and foreground phases and to estimate their interconnection. We have included the ILC , TOH FCM and the PCM maps in our analysis and have discovered significant deviation of the phase statistics from the expected for uniformly random phases. We have shown that such a kind of non-Gaussianity arises owing to the contamination of the foregrounds signal in the ILC, TOH FCM and the PCM maps.

## 2 Phase cleaning method and the CMB extraction

The fluctuations of measured CMB plus foregrounds radiation on a sky sphere can be expressed as a sum over spherical harmonics:

$$\mathrm{\Delta }T(\theta ,\phi )=T(\theta ,\phi )T=\underset{\mathrm{}=0}{\overset{\mathrm{}}{}}\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}a_\mathrm{}mY_\mathrm{}m(\theta ,\phi ),$$
(1)
where $`a_\mathrm{}m`$ are the coefficients of expansion, $`T`$=2.73K, and $`\mathrm{\Delta }T`$=0. Homogeneous and isotropic CMB Gaussian random fields (GRFs), as a result of the simplest inflation paradigm, possess $`a_\mathrm{}m`$ modes whose real and imaginary parts are independently distributed. The statistical properties of this field are completely specified by its angular power spectrum $`C_{\mathrm{}}`$,

$$a_\mathrm{}m^{cmb}(a^{cmb})_\mathrm{}^{^{}}m^{^{}}^{}=C_{\mathrm{}}^{cmb}\delta _{\mathrm{}\mathrm{}^{^{}}}\delta _{mm^{^{}}},$$
(2)
and random phases

$$\mathrm{\Psi }_\mathrm{}m^{cmb}=\mathrm{tan}^1\frac{Im(a_\mathrm{}m^{cmb})}{Re(a_\mathrm{}m^{cmb})},$$
(3)
which are uniformly distributed at the range $`0,2\pi `$. For the foregrounds, the signal is obviously non-Gaussian. So, for the combined CMB + foregrounds signal we define

$$a_\mathrm{}m=|a_\mathrm{}m|\mathrm{exp}(i\mathrm{\Psi }_\mathrm{}m),$$
(4)
where $`|a_\mathrm{}m|`$ is the modulus and $`\mathrm{\Psi }_\mathrm{}m`$ is the phase of each $`\mathrm{},m`$ harmonic. In practice, each of the WMAP K–W maps is decomposed into set of $`𝐚^j`$ where $`𝐚a_\mathrm{}m`$ and index $`j=1,\mathrm{2..5}`$ corresponds to K, Ka, …, W bands, using the HEALPix code (Górski et al. 1999).111In addition we used the GLESP pixelization scheme (Doroshkevich et al. 2003a) in order to control the accuracy of the $`𝐚^{(j)}`$ harmonics estimation.

The basic idea of the PCM is to generalize the TOH and TE minimization scheme including also the minimization of the cross-correlations between derived CMB signal and foregrounds (Naselsky et al. 2003). The method does not require any galactic cut–offs or disjoint regions. Following Naselsky et al. (2003), we consider the combinations of the WMAP maps: Ka–Q, Ka–V and Q–V, for which the higher correlations between foreground phases takes place. 222We did not include W and K bands to the separation procedure because of the peculiar phases (Naselsky et al. 2003).. For each pairs of the maps, we introduce weighting coefficients $`𝐰^{(j)}(\mathrm{})`$ similar to TOH and TE methods and minimize the variance of the derived map per each mode $`\mathrm{}`$. Neglecting the instrumental for $`\mathrm{}100`$ noise and taking into account the beam shape deconvolution, we consider the signal as a superposition of the CMB and the foregrounds.

From all $`𝐚^{(j)}`$, we have found the phases $`\mathrm{\Psi }_\mathrm{}m^{(j)}𝚿^{(j)}`$ (see Eq (2)–(4)). Each set of the $`𝐚^{(j)}`$ coefficients is defined by a combination of the different foreground coefficients $`_k𝐚_k^{(j)}=𝐆^{(j)}`$ and the CMB signal $`𝐚^{cmb}`$, where index $`k`$ marks the synchrotron, free-free and dust emission. For each combination of the maps, the derived CMB map is defined as

$$a_\mathrm{}m^M=\underset{j=1}{\overset{2}{}}𝐰^{(j)}(\mathrm{})a_\mathrm{}m^{(j)}=a^{cmb}+\underset{j=1}{\overset{2}{}}𝐰^{(j)}G^{(j)},$$
(5)
where

$`𝐰^{(1)}(\mathrm{})={\displaystyle \frac{_m\left\{|𝐚^{(2)}|\left(|𝐚^{(2)}||𝐚^{(1)}|\mathrm{cos}(𝚵^{(2)}𝚵^{(1)})\right)\right\}}{_m|𝐚^{(1)}𝐚^{(2)}|^2}},`$
$`𝐰^{(2)}(\mathrm{})={\displaystyle \frac{_m\left\{|𝐚^{(1)}|\left(|𝐚^{(1)}||𝐚^{(2)}|\mathrm{cos}(𝚵^{(2)}𝚵^{(1)})\right)\right\}}{_m|𝐚^{(1)}𝐚^{(2)}|^2}},`$ (6)
and $`𝚵^{(j)}`$ is the phase of $`j`$-channel.

The phases of the reconstructed CMB are related to the foregrounds amplitudes and phases as follows (Naselsky et al. 2003)

$$𝚿^M=\xi +\mathrm{arcsin}\frac{_j𝐰^{(j)}|𝐆^{(j)}|\mathrm{sin}(𝚿^{(j)}\xi )\mathrm{cos}(𝚿^{(M)})}{_j𝐰^{(j)}|𝐆^{(j)}|\mathrm{cos}𝚿^{(j)}+|𝐚^{cmb}|\mathrm{cos}\xi },$$
(7)
where $`\xi `$ is the true CMB phase and $`𝚿^{(j)}`$ are the foreground phases. As is mentioned in Naselsky et al. (2003), the phases of the reconstructed CMB signal $`𝚿^M`$ would have correlations with those of foregrounds. The main goal of the PCM is to minimize such correlations by minimizing of the weighting variance

$$V=\frac{1}{2\pi (2\mathrm{}+1)}\underset{m}{}\frac{|𝐚^{cmb}|^2}{C_{\mathrm{}}}_0^{2\pi }𝑑\xi \left(𝚿^M\xi \right)^2min$$
(8)
using Eq.(7) and weighting coefficients $`𝐰^j`$ in a form

$`𝐰^{(1)}={\displaystyle \frac{_m\left\{|𝐆^{(2)}|\left[|𝐆^{(2)}||𝐆^{(1)}|\mathrm{cos}(𝚿^{(2)}𝚿^{(1)})\right]\right\}}{_m|𝐆^{(1)}𝐆^{(2)}|^2}}`$
$`𝐰^{(2)}={\displaystyle \frac{_m\left\{|𝐆^{(1)}|\left[|𝐆^{(1)}||𝐆^{(2)}|\mathrm{cos}(𝚿^{(2)}𝚿^{(1)})\right]\right\}}{_m|𝐆^{(1)}𝐆^{(2)}|^2}}.`$ (9)
Note that minimization of the weighting variance Eq.(8) is equivalent to minimization of the error bars of the CMB reconstruction $`_m|a_\mathrm{}m^Ma_\mathrm{}m^{cmb}|^2\mathrm{min}`$, where $`a_\mathrm{}m^{cmb}`$ is the true CMB (Naselsky et al. 2003).

Practical implementation of the PCM includes the iteration scheme in which the first step corresponds to minimization of the variance of the derived map in each mode $`\mathrm{}`$ with the choice of coefficients Eq.(6). This step of optimization reconstructs the pre-CMB and pre-foreground map for each pair of channels where the foreground maps are a simple subtraction $`a_\mathrm{}m`$ coefficients of the K–W signals and derived pre-CMB $`a_\mathrm{}m`$. Then, using the moduli and phases of the pre-foregrounds we perform the next iteration, using Eqs.(9), (5) and so on. This iteration scheme is stable and reproduce quite well the CMB signal from each pair of K–W maps after two steps of iteration.

However, the CMB maps reconstructed from each pair of K–W maps are slightly different because of residues of the foregrounds. To minimize such residues, as the last step of the PCM we use the so called MIN-MAX filter for each pixels of the map. The MIN-MAX filter compares signals, $`\mathrm{\Delta }T_p`$, in each pixel $`p`$ in a set of maps. In this method, we chose the minimal amplitude of $`\mathrm{\Delta }T_p`$ in the pixel as the CMB signal, i.e. $`L_{\mathrm{min}}(\{\mathrm{\Delta }T_p^{(i)}\})\mathrm{\Delta }T_p^{\mathrm{min}}`$, where $`|\mathrm{\Delta }T_p^{\mathrm{min}}|=\mathrm{min}\{|\mathrm{\Delta }T_p^{(i)}|\}`$. To estimate residues from point sources, we find the maximal amplitude of the signal in each pixel $`L_{\mathrm{max}}(\{\mathrm{\Delta }T_p^{(i)}\})\mathrm{\Delta }T_p^{\mathrm{max}}`$, where $`|\mathrm{\Delta }T_p^{\mathrm{max}}|=\mathrm{max}\{|\mathrm{\Delta }T_p^{(i)}|\}`$, for all the pre-CMB maps produced by the method Eq.(9). For true CMB, the difference $`\mathrm{\Delta }T_p^{\mathrm{max}}\mathrm{\Delta }T_p^{\mathrm{min}}`$ is equal to zero and we consider $`\mathrm{\Delta }T_p^{\mathrm{max}}\mathrm{\Delta }T_p^{\mathrm{min}}`$ as a measure of deviation from the CMB map.

The reason for such a filter is quite obvious. The signal in each pixel is a superposition of the CMB signal and (small) residues of the foregrounds. If correlations between the CMB signal and foregrounds are minimal then we expect that all the deviations of the $`\mathrm{\Delta }T_p^{(i)}`$ in the pixel are caused by the residues. In our case, two approximations of the CMB map restored by Ka–V and Q–V channels (Naselsky et al. 2003) have been used and $`\mathrm{\Delta }T^{\mathrm{min}}`$ has been taken as the CMB map.

The result of the PCM application to the WMAP data at the range $`\mathrm{}100`$ is shown in Fig. 1 for the CMB signal. In this Figure, the differences between the ILC and PCM maps (top right), TOH FCM and PCM (bottom left) and between TOH Wiener and PCM maps (bottom right) are also shown. As is shown in Fig. 1, the main difference is connected with sources in the Galactic plane and the harmonic $`a_{51}`$ which is probably related to the galactic emission. The difference shows foreground components being absent in the PCM map.

For the simplest estimator of the power spectrum, $`C(\mathrm{})=(2\mathrm{}+1)^1_m|a_\mathrm{}m|^2`$, we show in Fig. 2 the power spectrum for the ILC, FCM, PCM and TOH Wiener filtered maps. As one can see from Fig.2 the PCM map reproduces well the spectrum for the TOH Wiener filtered map.

## 3 Correlations of the cleaned signals and foregrounds for the ILC, FCM and PCM maps

Three cleaned maps (ILC, FCM and PCM) are produced from different methods of the foreground component separation. For the perfect separation of the CMB signal and foregrounds we can expect small random correlations between their phases. However, these methods provide the approximate separation only. Therefore, these remaining cross–correlations between phases of the cleaned signal, $`\varphi _s`$, and the foregrounds, $`\psi _f`$, can be used to characterize the degree of separation achieved.

To do this, here we will use the simplest phase–phase circular (Fisher, 1993) correlation coefficients. We will use also the linear correlation coefficients and more refined “friend–of–friend” statistics discussed, for example, by Roeder (1992) neglecting circularity of the phases. Indeed, equivalence of the phases $`\varphi _s\&\varphi _s\pm 2\pi `$ and $`\psi _f\&\psi _f\pm 2\pi `$, changes their correlation functions and cumulants and the results obtained with the “friend–of–friend” statistics. To suppress this effect we compared our result with 200 random realizations of Poissonian process prepared in the same manner (see Sec. 3.3).

For five frequency channels K–W, the maps are taken from the WMAP web site333http://lambda.gsfc.nasa.gov/product/map/m\_products.cfm, and all the phases are obtained by the spherical harmonics decomposition using the HEALPix (Górski et al. 1999) and GLESP (Doroshkevich et al. 2003a) codes. We consider separately two ranges of multipoles, $`2\mathrm{}50`$ and $`51\mathrm{}100`$. For the ILC, FCM and PCM CMB maps, the derived foregrounds are determined as differences between the signal ($`S`$) and cleaned map ($`C`$) $`F=SC`$ for each ($`\mathrm{},m`$). In addition, for the ILC cleaned map, WMAP own foregrounds given in the same web–site were also used. Examples of difference of $`S`$, $`ILC`$ and WMAP own foregrounds for V and W WMAP bands are shown on Fig. 3.

### 3.1 Circular cross-correlation of the phases.

Following Fisher(1993) we define the statistics

$$M_{sp}=\frac{1}{\mathrm{}_{max}}\underset{m}{\overset{\mathrm{}_{max}}{}}\mathrm{exp}[ip(\varphi _m\varphi _m)],$$

$$M_{fp}=\frac{1}{\mathrm{}_{max}}\underset{m}{\overset{\mathrm{}_{max}}{}}\mathrm{exp}[ip(\psi _m\psi _m)],$$

$$\varphi _m=\mathrm{tan}^1\left(\underset{m}{\overset{\mathrm{}_{max}}{}}\mathrm{sin}\varphi _m/\underset{m}{\overset{\mathrm{}_{max}}{}}\mathrm{cos}\varphi _m\right),$$

$$\psi _m=\mathrm{tan}^1\left(\underset{m}{\overset{\mathrm{}_{max}}{}}\mathrm{sin}\psi _m/\underset{m}{\overset{\mathrm{}_{max}}{}}\mathrm{cos}\psi _m\right),$$

$$R_{sf}(\mathrm{})=\frac{1}{\mathrm{}}\underset{m=1}{\overset{\mathrm{}}{}}\mathrm{cos}(\varphi _m\psi _m),$$
(10)
$$R_{sf}=\frac{1}{\mathrm{}_{max}\mathrm{}_{min}+1}\underset{\mathrm{}=\mathrm{}_{min}}{\overset{\mathrm{}_{max}}{}}R_{sf}(\mathrm{})$$

where $`M_{sp}`$ and $`M_{fp}`$ are the $`p`$-th trigonometric moments of the samples, $`\varphi `$ and $`\psi `$ are corresponding mean directions, $`R_{sf}(\mathrm{})`$ is the circular cross-correlation coefficient in each mode $`\mathrm{}`$ and $`r_{sf}`$ is the mean circular cross-correlation coefficient for all phases. For $`m=0`$ and for all $`\mathrm{}`$ phases $`\varphi (\mathrm{},0)=\psi (\mathrm{},0)=0`$ and here we neglect them.

For the K–W bands the circular coefficients, $`R_{sf}(\mathrm{})`$, are plotted in Fig. 4 for the ILC cleaned signal and its own and derived foregrounds and for PCM signal and derived foregrounds. For the first three channels, these coefficients are quite moderate and do not exceed the random scatter ($`1\sigma `$) obtained from 200 random realizations. Note that, for all the bands, the shape of the functions $`R_{sf}(\mathrm{})`$ are quite similar to each other what reflects strong correlation of phases in all the foregrounds (Naselsky et al. 2003). As one can see from Fig. 4, for the V & W channels the cross–correlations of both the PCM and ILC phases with the derived foregrounds seems to be quite significant.

The same tendency is seen from the estimations of the mean coefficients, $`R_{sf}`$, for $`2\mathrm{}50`$ listed in Table 1 for all three cleaned signal and the own WMAP and three derived foregrounds. For first three channels these coefficients are small but they become significant for channels V and W and derived foregrounds.

Correlations of phases of different foregrounds with the ILC, FCM and PCM cleaned signals are plotted in Figs. 5–7. For the free–free foregrounds, correlations are quite moderate and they exceed the random scatter only for the 35 – 36 harmonics. However, for the synchrotron and dust foregrounds we see the significant (95%) correlations for the 21 – 23 harmonics for all cleaned signal maps. In all the cleaned maps this range corresponds to minima of the power (see Fig.2̇).

### 3.2 Linear cross-correlation of the phases

For $`\mathrm{}_{min}\mathrm{}\mathrm{}_{max}`$, the linear correlation coefficient between phases of the cleaned signal, $`\varphi _s`$, and foregrounds $`\psi _f`$, is defined as follows:

$$r_{sf}(\mathrm{}_{min},\mathrm{}_{max})=\frac{\varphi _s\psi _f\varphi _s\psi _f}{\sigma _\varphi \sigma _\psi },$$
(11)
$$\varphi _s=\frac{1}{N_{\mathrm{}}}\underset{\mathrm{}=\mathrm{}_{min}}{\overset{\mathrm{}_{max}}{}}\underset{m=1}{\overset{\mathrm{}}{}}\varphi _s(\mathrm{},m)\pi ,$$

$$\sigma _\varphi ^2=\frac{1}{N_{\mathrm{}}}\underset{\mathrm{}=\mathrm{}_{min}}{\overset{\mathrm{}_{max}}{}}\underset{m=1}{\overset{\mathrm{}}{}}\varphi _s^2(\mathrm{},m)\pi ^2/3,$$

$$\varphi _s\psi _f=\frac{1}{N_{\mathrm{}}}\underset{\mathrm{}=\mathrm{}_{min}}{\overset{\mathrm{}_{max}}{}}\underset{m=1}{\overset{\mathrm{}}{}}\varphi _s(\mathrm{},m)\psi _f(\mathrm{},m),$$

$$N_{\mathrm{}}=(\mathrm{}_{max}\mathrm{}_{min}+1)(\mathrm{}_{max}+\mathrm{}_{min})/2.$$

For the foreground phases, $`\psi _f`$, the mean value, $`\psi _f`$, and the variance, $`\sigma _\psi `$, are defined with similar relations. Here $`N_{\mathrm{}}`$ is the number of phases in the sample under consideration. For $`m=0`$ and for all $`\mathrm{}`$, phases $`\varphi _s(\mathrm{},0)=\psi _f(\mathrm{},0)=0`$ and here we neglect them.

The linear correlation coefficients for $`2\mathrm{}50`$ and $`51\mathrm{}100`$ are listed in Table 2 for the ILC, FCM and PCM cleaned signals and five frequency channels. First three columns characterize the correlation with the ILC own foregrounds. Six other columns represent these coefficients for the same signals and derived foregrounds. For six samples the correlations between the signal and derived foregrounds, $`r_{sf}(\mathrm{}_{min},\mathrm{}_{max})`$, increase with the channel frequency and are comparable for a given frequency. The random scatter of the coefficient determined by averaging of 200 random realizations is $`r_{sf}(\mathrm{}_{min},\mathrm{}_{max})0.06`$ at $`68\%`$ CL. As is seen from Table 1, for high frequency channels V and W the measured $`r_{sf}(\mathrm{}_{min},\mathrm{}_{max})`$ exceeds this value. This fact indicates the noticeable correlation between the cleaned signal and derived foregrounds in channels V and W and the limited precision of separation methods.

For all the channels and WMAP own foregrounds, the coefficients $`r_{sf}(\mathrm{}_{min},\mathrm{}_{max})`$ are less than the random value that indicates high efficiency of corrections of own foregrounds. However, this result depends upon the foregrounds rather than the cleaned signal. For the FCM and PCM cleaned signal and the same foregrounds, the coefficients, $`r_{sf}(\mathrm{}_{min},\mathrm{}_{max})`$ are also listed in the same Table. For all samples these coefficients $`r_{sf}(\mathrm{}_{min},\mathrm{}_{max})`$ do not exceed the random value. Therefore, the small correlations between the own foregrounds of ILC and the cleaned signals is determined by the properties of foregrounds.

By definition, this coefficient is more sensitive to $`\mathrm{}_{max}`$ because the majority of phases used comes from $`\mathrm{}\mathrm{}_{max}`$. However, it decreases when we correlate phases with different $`\mathrm{}`$.

### 3.3 Linear correlation coefficient of the phases per each multipole $`\mathrm{}`$

Linear correlation coefficient of phases of cleaned signal and foregrounds, $`r_{sf}(\mathrm{})`$, can also be found for each $`\mathrm{}`$ in the same manner as it was been done for the circular cross– correlation coefficient in Sec. 3.1. In the case, the coefficient is defined by the same relations (11) neglecting the summation over $`\mathrm{}`$ and for $`N_{\mathrm{}}=\mathrm{}`$. It allows one to characterize properties of each mode separately and to determine harmonics with maximal $`r_{sf}(\mathrm{})`$. However, the small sample statistics of phases at moderate $`\mathrm{}`$ increases significantly its random scatter. Note also that, as is seen from the definition, this coefficient characterizes the sample used. Thus, in contrast with the circular cross–correlation coefficient, the direct averaging $`r_{sf}(\mathrm{})`$ over $`\mathrm{}`$ does not reproduce results obtained in previous Section.

The coefficients $`r_{sf}(\mathrm{})`$ are plotted in Fig. 8 for channels V and W for two foregrounds of the ILC cleaned signal and for the PCM cleaned signal. For each $`\mathrm{}`$ the error bars are found for 200 random realizations and are well consistent with the expected behavior $`\mathrm{}^{1/2}`$.

For channels K–V and for all three samples, the correlation coefficients do not exceed the random scatter but the disposition of their maxima and minima are very similar. For the channel W we see significant differences between functions $`r_{sf}(\mathrm{})`$ for different samples. More pronounced correlations are seen for $`4\mathrm{}10`$, $`\mathrm{}=1315`$, $`\mathrm{}=1827`$, $`\mathrm{}=3236`$ and $`\mathrm{}=4550`$ As is seen from the comparison with Fig. 2, some of these harmonics correspond to local minima in $`C_{\mathrm{}}`$. This fact indicates that, perhaps, these peculiarities are caused by the same sources as the Galactic dust emission.

### 3.4 “Friend–of–Friend” statistics

To compare the foregrounds and the foreground-cleaned signal with random Poissonian process we use a more refined technique, the “friend–of–friend” approach for the 2D phase diagramme. In these diagrammes each point has coordinates $`x_i=\varphi _s(\mathrm{},m),y_i=\psi _f(\mathrm{},m)`$ where $`\varphi _s(\mathrm{},m)`$ and $`\psi _f(\mathrm{},m)`$ are the phases of the cleaned signal and the foregrounds, respectively.

For each diagramme, we find the fraction of clusters, $`f_M`$, with a given richness, $`M`$, in wide range of linking lengths, $`r_{lnk}`$, which determines the maximal point separation within a cluster. Here we use the normalized linking length,

$$l_{lnk}=\pi r_{lnk}^2n,n=N_{phs}/4\pi ^2,$$
(12)
where $`r_{lnk}`$ is the actual distance between two points, and $`N_{phs}`$ is the number of points used. As was discussed in White (1979) and Borgani (1996), the fractions $`f_M`$ depend on the correlation functions (or cumulants) of all orders that provides high sensitivity to deviations of the sample from the Poissonian–like one. This sensitivity is weak for $`M=1\&\mathrm{\hspace{0.17em}2}`$ and fast increases for larger $`M`$. But the random scatter of $`f_M`$ increases for larger $`M`$ as well, and the analysis becomes ineffective for $`M7`$.

For the truly Poissonian sample of points, the correlations functions of all orders are equal to zero and these fractions can be approximately described as follows:

$$f_M(l_{lnk})[1\mathrm{exp}(l_{lnk})]^M\mathrm{exp}(l_{lnk}),$$

where the first term describes the probability to find $`M`$ separation between points $`l_{lnk}`$ while the last term gives probability to find larger separation for the $`M+1`$ point. The scatter of $`f_M`$ is determined by averaging of $`f_M`$ over 200 random realizations of diagrammes analyzed in the same manner.

To find clusters with a given richness we use the Minimal Spanning Tree technique which, firstly, allows one to connect all points within a unique tree. Further on, rejection of edges of a tree larger than a chosen linking length transforms the tree to the system of clusters of all richness for a given linking length. The Minimal Spanning Tree technique was proposed in Barrow et al. (1985) and van de Weygaert (1991) and was applied for analysis of rich galaxy catalogues (see, e.g., Doroshkevich et al. 2001, 2003b).

Firstly, we are testing the hypothesis that the distribution of phases are Poissonian-like and correlations are negligible for both foregrounds and cleaned signal. Main results of our analysis are presented in Figs. 9,10 where for all five frequency channels variations of fractions $`f_M(l_{lnk}),1M7`$ are plotted versus the linking length, $`l_{lnk}`$. As is expected, the functions $`f_M(l_{lnk})`$ are quite similar to Poissonian ones for $`M=1`$ and $`M=2`$, and even for $`M=3`$ the deviations are moderate. However, for $`M4`$ the deviations from the Poissonian samples become significant, especially for V and W channels. This result illustrates the expected non-Gaussian character of the cleaned signal and its correlations with the foregrounds.

This result is expectable because of non-Gaussian character of the foregrounds signals. However, comparing the phase diagramme for the cleaned signal and foregrounds with 100 diagrammes prepared for the foreground and random phases we can quantify their divergences using $`\chi ^2`$ statistics. These results are shown in Table 3. If for $`M=1`$ and $`M=2`$ we see small $`\chi ^2N`$ where $`N`$ is the number of linking lengths used then for $`M3`$, we have $`\chi ^2N`$. Larger ratio $`\chi ^2/N`$ is found again for higher frequency channels V and W. This fact indicates that all the cleaned maps contain some residues from the dust emission.

## 4 Conclusion

In this paper we have proposed new methods of investigation the statistical properties of the signal derived from the WMAP K–W bands foreground cleaned maps. With this approach we show significant correlations between the phases of the cleaned maps and the foregrounds. We have compared three cleaned maps, ILC, TOH FCM and PCM, and have investigated variations of the linear cross-correlation coefficients versus the multipole index, $`\mathrm{}`$. We suggest also more sophisticated “friend–of–friend” statistics. All the methods described in Sec. 3 show significant cross-correlations between the phases of the cleaned maps and the foregrounds which manifest themselves more clearly for the W band of the WMAP. We have pointed out that some of peculiarities of the ILC and PCM power spectrum, for example, local minima at $`\mathrm{}6`$, and $`\mathrm{}2025`$ are accompanied by significant cross-correlations between phases of the foregrounds and cleaned signal. Such correlations are important indicators for investigation and detection of possible non-Gaussianity of the CMB signal which could be a byproduct of the component separation methods.

## Acknowledgments

This paper was supported by Danmarks Grundforskningsfond through its support for the establishment of the Theoretical Astrophysics Center. We thank Max Tegmark et al. for providing their processed maps and making them public with openness. We thank Igor Novikov and Lung-Yih Chiang for useful discussions. We thank TAC CMB collaboration for used GLESP code. We also acknowledge the use of HEALPix package (Górski et al. 1999)444http://www.eso.org/science/healpix/ to produce $`a_\mathrm{}m`$ from WMAP maps and some figures.

## 5 Appendix

In this section we would like to introduce an analytical approach for investigation of the cross-correlation coefficient $`r_{sf}(\mathrm{})`$ between foregrounds and derived CMB phases per each mode $`\mathrm{}`$. We define $`r_{sf}(\mathrm{})`$ as

$$r_{sf}^k(\mathrm{})\frac{3}{4\pi ^4l}\underset{m=1}{\overset{\mathrm{}}{}}_\pi ^\pi _\pi ^\pi 𝑑𝚿^M𝑑\mathrm{\Psi }^k𝚿^M\mathrm{\Psi }^k$$
$`(A1)`$
where $`k`$ marks the foreground derived from the WMAP frequency band. Here we assume that the distributions of both phases are close to homogeneous ones with $`\mathrm{\Psi }`$=0, $`\mathrm{\Psi }^2\pi ^2/3`$ and we average a correlation coefficient over $`m`$ for a given $`\mathrm{}`$. For the phase of the cleaned CMB signal we will use Eq.(7), which reflects directly the cross-correlation between the derived CMB and foregrounds phases. Then, using the definition $`\mathrm{\Psi }^j=\alpha _{jk}\mathrm{\Psi }^k`$, where $`\alpha _{jk}`$ is the cross-correlation coefficient between the foreground phases in $`j`$-th and $`k`$-th channels, from Eq.(A1) we obtain

$$r_{sf}^k(\mathrm{})=\frac{3}{\sqrt{2}\pi \mathrm{}}\underset{m=1}{\overset{\mathrm{}}{}}\underset{j}{}\frac{𝐰^{(j)}|𝐆^{(j)}|}{|𝐚^{cmb}|}\frac{J_{\frac{3}{2}}(\pi \alpha _{jk})}{(\alpha _{jk})^{\frac{1}{2}}}$$
$`(A2)`$
Let us assume that cross-correlation between different foreground phases is very strong and corresponding coefficients are $`\alpha _{ik}1`$555see Naselsky et al. (2003) for details of the phase correlations.. Then, from (A2) we obtain

$$r_{sf}^k(\mathrm{})=\frac{3}{\pi ^2\mathrm{}}\underset{m=1}{\overset{\mathrm{}}{}}\frac{1}{|𝐚^{cmb}|}\underset{j}{}𝐰^{(j)}|𝐆^{(j)}|$$
$`(A3)`$
For the PCM method, $`𝐰^{(j)}`$ coefficients are the functions of $`\mathrm{}`$, but not $`m`$, while ratio $`|𝐆^{(j)}|/|𝐚^{cmb}|`$ depends on $`m`$ for a given value $`\mathrm{}`$.

From Eq.(5) we get

$$|𝐚^M𝐚^{cmb}|=\left|\underset{j}{}𝐰^{(j)}𝐆^{(j)}\right|$$
$`(A4)`$
and, because the foregrounds seem to be highly correlated,

$$\left|\underset{j}{}𝐰^{(j)}𝐆^{(j)}\right|=sign\left(\underset{j}{}𝐰^{(j)}|𝐆^{(j)}|\right)\underset{j}{}𝐰^{(j)}|𝐆^{(j)}|,$$

$$r_{sf}^k(\mathrm{})=\frac{3}{\pi ^2\mathrm{}}\underset{m=1}{\overset{\mathrm{}}{}}\frac{|𝐚^{(M)}𝐚^{cmb}|}{|𝐚^{cmb}|}sign(|𝐚^{(M)}||𝐚^{cmb}|)(A5)$$

Thus, the cross-correlation coefficient $`r_{sf}^k(\mathrm{})`$ simply describes the accuracy of the reconstruction of the CMB signal.