text:
# New insights on the thermal dust from the far-infrared to the centimeter

## 1 Introduction

To accurately characterize dust emissivity properties represents a major challenge of nowadays astronomy. In the submillimeter domain, large grains at thermal equilibrium (e.g. desert90 ) dominate the dust emission. This thermal dust is characterized by a temperature and a spectral dependence of the emissivity which is usually simply modelled by a spectral index. The temperature, density and opacity of a molecular cloud are key parameters which control the structure and evolution of the clumps, and therefore, star formation. The spectral index ($`\beta `$) of a given dust grain population is directly linked to the internal physical mechanisms and the chemical nature of the grains.

It is generally admitted from Kramers-König relations that 1 is a lower limit for the spectral index. $`\beta `$ = 2 is particularly invoked for isotropic crystalline grains, amorphous silicates or graphitic grains. However, it is not the case for amorphous carbon, which is thought to have a spectral index equal to 1. Spectral indices above 2 may exist, according to several laboratory measurements on grain analogs. Observations of the diffuse interstellar medium at large scales favour $`\beta `$ around 2 (e.g. boulanger96 , dunne01 ). In the case of molecular clouds, spectral indices are usually found to be between 1.5 and 2. However, low values (0.2-1.4) of the spectral index have been observed in circumstellar environments, as well as in molecular cloud cores.

## 2 The dust spectral index as seen by PRONAOS

PRONAOS (PROgramme NAtional d’Observations Submillimétriques) is a French balloon-borne submillimeter experiment ristorcelli98 . Its effective wavelengths are 200, 260, 360 and 580 $`\mu `$m, and the angular resolutions are 2 in bands 1 and 2, 2.5 in band 3 and 3.5 in band 4. The data analyzed here were obtained during the second flight of PRONAOS in September 1996, at Fort Sumner, New Mexico. This experiment has observed various phases of the interstellar medium, from diffuse clouds in Polaris bernard99  and Taurus stepnik03  to massive star-forming regions in Orion (ristorcelli98 , dupac01 ), Messier 17 dupac02 , Cygnus B, and the dusty envelope surrounding the young massive star GH2O 092.67+03.07 in NCS. The $`\rho `$ Ophiuchi low-mass star-forming region has also been observed, as well as the edge-on spiral galaxy NGC 891 dupac03b .

We fit a modified black body law to the spectra: $`I_\nu =ϵ_0B_\nu (\lambda ,T)(\lambda /\lambda _0)^\beta `$, where $`I_\nu `$ is the spectral intensity (MJy/sr), $`ϵ_0`$ is the emissivity at $`\lambda _0`$ of the observed dust column density, $`B_\nu `$ is the Planck function, $`T`$ is the temperature and $`\beta `$ is the spectral index. In most of the areas, we use either only PRONAOS data or PRONAOS + IRAS 100 $`\mu `$m data. We restrain the analysis to all fully independent (3.5 side) pixels for which both relative errors on the temperature and the spectral index are less than 20%. This procedure allows to reduce the degeneracy effect between the temperature and the spectral index. Dupac et al. (2001, 2002) have shown by fitting simulated data that this artificial anticorrelation effect was small compared to the effect observed in the data.

We present in Fig. 1 the spectral index - temperature relation observed. The temperature in this data set ranges from 11 to 80 K, and the spectral index also exhibits large variations from 0.8 to 2.4. One can observe an anticorrelation on these plots between the temperature and the spectral index, in the sense that the cold regions have high spectral indices around 2, and warmer regions have spectral indices below 1.5. In particular, no data points with $`T>`$ 35 K and $`\beta >`$ 1.6 can be found, nor points with $`T<`$ 20 K and $`\beta <`$ 1.5. This anticorrelation effect is present for all objects in which we observe a large range of temperatures, namely Orion, M17, Cygnus and $`\rho `$ Ophiuchi. It is also remarkable that the few points from other regions are well compatible with this general anticorrelation trend. The temperature dependence of the emissivity spectral index is well fitted by a hyperbolic approximating function.

Several interpretations are possible for this effect: one is that the grain sizes change in dense environments, another is that the chemical composition of the grains is not the same in different environments and that this correlates to the temperature, a third one is that there is an intrinsic dependence of the spectral index on the temperature, due to quantum processes such as two-level tunneling effects. Additional modeling, as well as additional laboratory measurements and astrophysical observations, are required in order to discriminate between these different interpretations. More details about this analysis and the possible interpretations can be found in dupac03a .

## 3 The Galactic dust emission as seen by COBE and WMAP

We combine data sets from COBE (DIRBE data points and FIRAS spectrum) and WMAP in order to investigate the Galactic dust spectral energy distribution on the whole sky. We convolve all data to the FIRAS 7o angular resolution, in order to obtain a consistent data set. This is done by co-adding higher resolution data from DIRBE and WMAP into 7 degree-sided pixels. The resulting spectrum is shown in Fig. 2 for the Galactic center area. It is well fitted by a 21 K, $`\beta `$=1.86 one-dust component model. These values are well compatible with the PRONAOS temperature - spectral index relation (see Section 2). More results and conclusions from this COBE - WMAP joint analysis will soon be available.

## Index


format: text

text:
# SAURON dynamical modeling of NGC 2974

## 1. Introduction

NGC 2974 is a nearby field elliptical galaxy classified as E4 with detected HI and H$`\alpha `$ gas apparently distributed in a disk (Kim et al. 1988; Buson et al. 1993). An earlier study by Cinzano & van der Marel (1994) suggested the galaxy harbors an embedded stellar disk not visible in the photometry. The galaxy was observed with the integral-field spectrograph SAURON (Bacon et al. 2001) mounted on the 4.2m WHT on La Palma, as a part of the survey of a representative sample of nearby E, S0 and Sa galaxies (de Zeeuw et al. 2002). With this individual study we explore the dynamical structure of the object, give a quantitative analysis of the embedded disk structure and investigate the distribution and dynamics of the gas in NGC 2974. A similar analysis is being performed on the galaxies of the sample, providing constraints on galaxy formation scenarios.

## 2. Stellar dynamical model

The SAURON observations of emission-lines are consistent with the gas being distributed in a circular disk, except in the inner 3<sup>′′</sup> where a bar-like perturbation was detected (Emsellem, Ferruit, & Goudfrooij, 2003). By assuming the gas is in a thin disk and fitting to the two-dimensional kinematics we are able to constrain the inclination of the disk to $`i60^{}`$. Adopting that inclination we construct an axisymmetric model for the galaxy stellar body using a three-integral orbit superposition Schwarzschild modeling technique adapted for integral field data (e.g. Cretton et al. 1999, Cappellari et al. 2002, Verolme et al 2002).

The best fitting model (Fig. 1) is able to reproduce all the features of the observed stellar kinematics. Analyzing the orbital distribution of the model we find two clearly distinct kinematic components. One component consists of high angular momentum orbits, producing a flattened density distribution, and the other component is rounder and rotates slowly. We associate the former component with a stellar disk and the latter with the spheroid. It is possible to perform a disk-bulge decomposition in phase space and estimate their mass fraction which turns out to be about 1:3.

We construct a dynamical model for the gas with the same potential and inclination used for the stellar model. Solving the Jeans equations in the limit of the asymmetric drift approximation we are able to reproduce the large scale gas kinematics providing an independent consistency test of the potential used for the stellar modeling.

## References

Bacon, R. et al. 2001, MNRAS, 326, 23

Buson, L.M. et al. 1993, A&A, 280, 409

Cappellari, M. et al. 2002, ApJ, 578, 787

Cinzano, P., & van den Marel, R.P. 1994, MNRAS, 270, 325

Cretton, N., de Zeeuw, P.T, van der Marel, R.P, & Rix, H. 1999, ApJS, 124, 383

Emsellem, E., Ferruit, P., & Goudfrooij, P. 2003, A&A, in press.

Emsellem, E. et al. 2003, MNRAS, submitted.

Kim, D.W. et al. 1988, ApJ, 330, 684

Verolme, E.K. et al. 2002, MNRAS, 335, 517

de Zeeuw, P.T. et al. 2002, MNRAS, 329, 513


text:
# Propagation of a mutually incoherent optical vortex pair in anisotropic nonlinear media

## 1 Introduction

Isotropic media with a repulsive nonlinearity support localized vortex solitons characterized by a dislocation of the wavefront at the point where the field amplitude vanishes . Vortex solitons have been observed in several types of defocusing nonlinear media, including superfluids, superconductors, optical Kerr media and Bose Einstein condensates. Stability of the solitons depends on the structure of the nonlinearity and the dimensionality of the medium. Planar one-transverse dimensional \[(1+1)D\] solutions are modulationally unstable in bulk, two-transverse dimensional \[(2+1)D\] media . Kerr type optical media with a cubic, isotropic, and local nonlinearity support stable (2+1)D vortex solitons . These vortex solitons have attracted considerable attention theoretically and have recently been studied in detail in atomic condensates. A Kerr type nonlinearity is, however, a simplified idealized model of a nonlinear response. The bulk photorefractive nonlinear medium used in the experiments reported here exhibits a nonlinearity that is both anisotropic and nonlocal. This leads to qualitative differences in the spatial dynamics of vortex beams.

One of the most basic phenomena is the propagation of an isolated vortex of unit topological charge. Such a vortex can form a stable soliton in isotropic defocusing media. Theory and experiments have shown that unit charged vortices become delocalized when they propagate in photorefractive media. The nonlinearity in a photorefractive crystal with an externally applied electric field, or photogalvanic response, is anisotropic and nonlocal. The nonlinearity induced by a localized, circular beam is roughly 3 times stronger along the direction of the applied field than in the perpendicular direction (we will refer to the direction of strongest nonlinearity as the $`\widehat{z}`$ axis). Due to the anisotropy a vortex with an initially azimuthally symmetric core profile focuses perpendicular to $`\widehat{z},`$ and stretches along $`\widehat{z}`$, so that the major axis of the core coincides with the direction of greatest material nonlinearity. Simultaneously the elongated vortex starts to rotate due to its phase structure. The direction of rotation is uniquely determined by the vortex charge; changing its sign changes the sign of the rotation. Eventually the rotation is stopped by the anisotropy so that the major axis of the vortex is aligned at some angle with respect to $`\widehat{z}.`$ The stretching, however, proceeds unchecked so that the vortex becomes more and more delocalized. Both theory and experiment show that delocalization of the vortex core is generic in anisotropic media, and not dependent on a specific choice of parameters. The implication is that localized optical vortex solutions and, in particular, vortex solitons of unit topological charge do not exist in these media.

An additional hallmark of anisotropy is the nonlinear decay of a charge $`n`$ vortex into $`n`$ unit charge vortices. The decay, although expected on energetic grounds, does not occur in isotropic media where a high charge vortex is metastable, but was observed using a photorefractive crystal as the nonlinear medium. Upon breakup of the input high-charge vortex the resulting charge-one vortices repel each other and form an array aligned perpendicular to $`\widehat{z}`$. The decay is driven by the intrinsic anisotropy of the medium and takes place for any core profile of the input vortex field. The theory developed in Ref.  suggests that the decay of high-charge vortices is possible in local isotropic media provided some anisotropy is introduced in the problem via, e.g., initial boundary conditions.

Given the instability of unit and high charge vortices in anisotropic media it is natural to ask if there exist self-bound field configurations that remain localized under propagation in anisotropic media. Previous experiments have demonstrated that the effects of anisotropy are considerably weakened for a counterrotating pair of vortices with zero net topological charge. When such a pair is aligned perpendicular to $`\widehat{z},`$ it forms a bound state that translates parallel to $`\widehat{z},`$ in a manner that is similar to the translational motion of a counterrotating point vortex pair in fluid dynamics. The translational motion may, however, be inconvenient in the context of optical processing applications. Furthermore, the initial orientation of the pair is crucial for its subsequent evolution. A pair aligned along $`\widehat{z}`$ annihilates due to diffraction.

In this paper we propose, and demonstrate experimentally a novel technique of creating bound vortex structures in anisotropic media that are free of the above mentioned limitations. The main idea is based on the fact that in the prevailing majority of nonlinear media the material response is slow as compared to the frequency of the optical field, and therefore is a function of the time-averaged intensity of this field. If the light field consists of several features separated by frequencies that are fast on the time scale of the material response, the material response will be a function of the sum of these features intensities, whereas the cross-terms oscillating at high frequencies can be neglected. This fact has been known in photorefractive nonlinear optics and used successfully for implementing several optical processing devices , as well as (1+1)D solitary structures .

We propose to create a bound vortex pair consisting of two copropagating counterrotating vortices sitting on top of each other that are mutually incoherent at the time scales of the nonlinearity. This incoherence can be achieved either by separating the carrier frequencies of the vortices by an amount that is larger than the inverse relaxation time of the medium, or profiling time histories of the input vortex fields such that their time overlap integrated over the relaxation time of the medium is zero. The mutual incoherence of the vortex fields removes the translational motion due to the linear coherence of the vortex pair. At the same time each of the vortices individually tries to rotate in opposite directions which effectively controls the anisotropy induced stretching of a single vortex.

## 2 Theory

In the theoretical analysis we use the set of equations describing propagation of an electromagnetic field $`B(\stackrel{}{r})`$ in a photorefractive self-focusing or self-defocusing medium as developed in . In our case the field consists of two temporal features. For definiteness assume that these features are separated by a frequency shift $`\mathrm{\Omega }`$ such that $`\mathrm{\Omega }\tau 1`$, where $`\tau `$ is the characteristic response time of the nonlinearity:

$$B(\stackrel{}{r},t)=\left[B_+(\stackrel{}{r},t)+B_{}(\stackrel{}{r},t)\mathrm{exp}(i\mathrm{\Omega }t)\right]\mathrm{exp}(ikxi\omega t).$$
(1)
In the steady state the equations governing spatial evolution of the amplitudes $`B_\pm `$ take the form

$`\left[{\displaystyle \frac{}{x}}{\displaystyle \frac{i}{2}}^2\right]B_+(\stackrel{}{r})=i{\displaystyle \frac{\phi }{z}}B_+(\stackrel{}{r}),`$ (2a)
$`\left[{\displaystyle \frac{}{x}}{\displaystyle \frac{i}{2}}^2\right]B_{}(\stackrel{}{r})=i{\displaystyle \frac{\phi }{z}}B_{}(\stackrel{}{r}),`$ (2b)
$`^2\phi +\phi \mathrm{ln}(1+|B_+|^2+|B_{}|^2)={\displaystyle \frac{}{z}}\mathrm{ln}(1+|B_+|^2+|B_{}|^2).`$ (2c)
Here $`=\widehat{y}(/y)+\widehat{z}(/z)`$ and $`\phi `$ is the dimensionless electrostatic potential induced by the light with the boundary conditions $`\phi (\stackrel{}{r}\mathrm{})0`$. The dimensionless coordinates $`(x,y,z)`$ are related to the physical coordinates $`(x^{},y^{},z^{})`$ by the expressions $`x=\alpha x^{}`$ and $`(y,z)=\sqrt{k\alpha }(y^{},z^{})`$, where $`\alpha =(1/2)kn^2r_{\mathrm{eff}}E_{\mathrm{ext}}`$. Here $`k`$ is the wave number of light in the medium, $`n`$ is the index of refraction, $`r_{\mathrm{eff}}`$ is the effective element of the electro-optic tensor, and $`E_{\mathrm{ext}}`$ is the amplitude of the external field directed along the $`z`$ axis far from the beam. The normalized intensity $`I=|B_+|^2+|B_{}|^2`$ is measured in units of saturation intensity $`I_d,`$ so that the physical beam intensity is given by $`I\times I_d.`$ The minus sign on the right hand side of Eqs. (2a,b) corresponds to a self-defocusing nonlinearity.

Numerical solutions of Eqs. (2) describing nonlinear evolution of a single vortex ($`B_{}=0`$) for different values of the applied voltage (nonlinearity) are shown in Fig. 1. The numbers on the frames are the values of this voltage in volts. Superimposed on the images are the equal intensity contour lines visualizing distortions of the vortex core for different values of the nonlinearity. The size of the frames is about $`200\mu \mathrm{m}`$ .

The input field was taken to be

$$B_+(x=0,\stackrel{}{r})=\sqrt{I_{\mathrm{in}}}r\mathrm{exp}(r^2/d_G^2+i\theta )$$
(3)
where $`\theta `$ is the azimuthal angle, $`r=\sqrt{y^2+z^2}`$ and $`d_G`$ is the diameter of the Gaussian beam. The parameters of the calculation were the following: the wavelength $`\lambda =0.63\mu \mathrm{m}`$, the refractive index $`n=2.3`$, the length of the nonlinear medium $`l=2\mathrm{cm}`$, the effective electrooptic coefficient $`r_{\mathrm{eff}}=130\mathrm{pm}/\mathrm{V}`$, $`d_G=115\mu \mathrm{m}`$ ($`135\mu \mathrm{m}`$ Full Width at Half Maximum) and $`I_{\mathrm{in}}=1`$.

The output intensity distribution in the absence of the nonlinearity is given by the frame labeled $`0\mathrm{V}`$ and corresponds to an annular ring having approximately $`79\mu \mathrm{m}`$ and $`267\mu \mathrm{m}`$ internal and external diameters, respectively. The frames corresponding to nonzero nonlinearity demonstrate focusing of the vortex along $`\widehat{y}`$ and its stretching along $`\widehat{z}.`$ Also clearly seen is the rotation and alignment of the vortex. This rotation is charge-dependent and changes sign if the charge of the vortex is changed from plus to minus one. The magnitude of all these effects is directly proportional to the value of the applied voltage.

Numerical solutions of Eqs. (2) describing nonlinear evolution of the mutually incoherent vortex pair for different values of the applied voltage (nonlinearity) are shown in Fig. 2. The input field was taken to be

$$B_\pm (x=0,\stackrel{}{r})=\sqrt{I_{\mathrm{in}}}r\mathrm{exp}(r^2/d_G^2\pm i\theta )$$
(4)
All parameters were the same as in the case of a single vortex (Fig. 1). The parameter $`I_{\mathrm{in}}`$ was again chosen such that the total maximum input intensity was equal to one (the maximum intensity of each of the constituents of the pair was 0.5). Figure 2 demonstrates some considerable reduction in the magnitude of anisotropy effects as compared to the case of a single vortex (Fig.1). In particular, the degree of ellipticity of the vortex core at high voltages is several times smaller in Fig. 2 than in Fig. 1. The horizontal intensity dip appearing at high intensities on the Gaussian beam and passing from left to right is several times smaller in Fig.2 than in Fig. 1. Figure 2 also confirms that the vortex pair remains stationary and does not translate with respect to the Gaussian beam.

The above stability of the vortex pair is due to its phase structure. To prove this point we have carried out numerical calculations for the input field consisting of two temporal features analogous to Eqs. (4) and (1) but without azimuthal phase dependence

$$B_\pm (x=0,\stackrel{}{r})=\sqrt{I_{in}}r\mathrm{exp}(r^2/d_G^2)$$
(5)
The field (5) has an input intensity distribution that is identical to the above discussed cases but lacks the topological phase structure present in the case of Eqs. (3) and (4). Figure 3 shows results of the calculations. All parameters are the same as in Figs. 1 and 2. In sharp contrast to Figs. 1 and 2 the output distribution of the field does not contain any intensity zeros because of the diffraction. The output intensity for zero applied voltage is a bright ring with a smaller local maximum in the center that is about two times weaker than the ring. Increasing nonlinearity results in the transfer of energy from the ring to the center of the beam and the appearance of two intensity minima on the right and left. The intensity in these minima at the highest value of the applied voltage ($`700\mathrm{V}`$) is about 0.2 of the maximum intensity in the center.

## 3 Experiment

The experimental setup is shown in Fig. 4. A 30 mW He-Ne laser beam ($`\lambda =0.63\mu \mathrm{m}`$) was passed through a variable beam splitter, and a system of lenses controlling the size of the beam waist. Diffraction of the beam by a diffractive element with a fringe dislocation created unit charge vortices of opposite signs in the first order diffracted beams to either side of the transmitted beam. One of the beams was phase modulated by reflection from a mirror mounted on a piezo-electric transducer driven by a sawtooth voltage, such that $`\mathrm{\Omega }\tau 1.`$ The beams were then recombined and directed into a photorefractive crystal of SBN:60 doped with 0.002% by weight Ce. The beams propagated perpendicular to the crystal $`\widehat{c}`$-axis ($`=z`$ axis), and were polarized along it. The crystal measured $`19\mathrm{mm}`$ along the direction of propagation, and was $`5\mathrm{mm}`$ wide along the $`\widehat{c}`$-axis. The experimentally measured value of the relevant component of the electrooptic tensor was found to be equal to $`r_{33}=130\mathrm{pm}/\mathrm{V}`$. A variable dc voltage was applied along the $`\widehat{c}`$-axis to control the value of the nonlinearity coefficient. The crystal was illuminated by a source of incoherent white light to increase the level of the effective saturation intensity. Images of the beams at the output face of the crystal were recorded with a CCD camera.

Figure 5 presents experimental output intensity distributions of the light beam with two embedded overlapping mutually incoherent vortices for different values of the applied voltage. The numbers on the frames give the applied voltage in volts. The size of the frames is about $`200\mu `$m. The output intensity distribution in the frame with no applied voltage corresponds to an annular ring with internal and external diameters of about $`80`$ and $`260`$ $`\mu `$m, respectively. Increasing the voltage results in the focusing of the vortex core and its stretching along the direction of the anisotropy which is clearly seen on the last frame corresponding to the $`700`$ V applied voltage. This stretching, however, is much smaller than that for a single vortex for similar parameters. Experimental data on the distortions of a single vortex are given in Ref. . Comparison with the results of the theoretical description shows good quantitative agreement between the experiment and the theory regarding the degree of ellipticity of the vortex core. Positions of intensity lobes (maxima) on the Gaussian beam above and below the vortex core in the theoretical Fig. 3 also are in agreement with those on the experimental Fig. 5.

## 4 Conclusions

In summary we presented an experimental and theoretical study of the propagation of a mutually incoherent pair of vortices. Incoherent coupling between the oppositely rotating vortices creates a symmetric attractive potential that arrests the rapid spreading and decay of an isolated vortex.

M. S. was supported by the A. P. Sloan Foundation. We are grateful to E. Eilertsen and E. Rasmussen at Risø National Laboratory for preparation of the holographic mask.

## References


text:
# Monitoring Chandra observations of the quasi-persistent neutron-star X-ray transient MXB 1659–29 in quiescence: the cooling curve of the heated neutron-star crust

## 1 Introduction

Neutron stars in low-mass X-ray binaries accrete matter from solar mass companions. Among those systems, the sub-group of neutron-star transients spend most of their time in quiescence during which hardly any or no accretion occurs. However, these transients sporadically become very bright ($`>`$$`10^{3638}`$ erg s-1) owing to a huge increase in the accretion rate onto their neutron stars. During those outbursts, these sources can be readily studied with the available X-ray instruments, but obtaining high quality quiescent data remains a challenge. In spite of this, several systems have now been studied in detail: they typically exhibit 0.5–10 keV luminosities of $`10^{3233}`$ erg s-1 and their spectra are usually dominated by a soft component which can be described by a thermal model. This emission is thought to be due to the cooling of the neutron star which has been heated during the outbursts (Brown, Bildsten, & Rutledge 1998; Campana et al. 1998a).

Most neutron-star transients are active for only weeks to months, but several systems have remained active for years and even decades (the ’quasi-persistent’ neutron-star transients; Wijnands et al. 2003). Wijnands et al. (2001) realized that those systems are excellent targets to study the effects of accretion on the behavior of neutron stars by observing them in quiescence. The accreting material is expected to have a larger effect on the neutron stars in such systems than on the neutron stars in short-duration transients (Wijnands et al. 2001; Rutledge et al. 2002). In the latter systems, the crust is only marginally heated during the outbursts and will quickly return to thermal equilibrium with the core after the end of the outbursts. In the quasi-persistent transients, however, the crust is heated to high temperatures and becomes significantly out of thermal equilibrium with the core (Rutledge et al. 2002). After the end of the prolonged outbursts, it will cool until it returns to equilibrium with the core. The exact cooling time depends on the thermal conductivity of the crust, the core cooling processes, and the accretion history of the source.

KS 1731–260 was the first quasi-persistent transient to be studied in detail in quiescence. It was observed using Chandra shortly after the end of its $``$12.5 year outburst (Wijnands et al. 2001) and it was found to have a luminosity of $`10^{33}`$ erg s-1 (for a distance $`d=7`$ kpc; 0.5–10 keV). Half a year later it was observed with XMM-Newton and it was found that its luminosity had decreased by a factor of 2–3 (Wijnands et al. 2002b). Using the cooling curves calculated by Rutledge et al. (2002), this drop in brightness can be explained if the neutron star has a large crustal conductivity and enhanced core cooling processes. In September 2001, a second quasi-persistent neutron-star transient (MXB 1659–29) turned off after having accreted for $``$2.5 years. Wijnands et al. (2003) obtained a Chandra observation of this source within a month after the end of its outburst and detected it at a luminosity of $`34\times 10^{33}`$ erg s-1 (0.5–10 keV; $`d=10`$ kpc). Several years before this outburst, the source was observed with ROSAT, but could not be detected (Verbunt 2001). The flux upper limit was $``$10 times lower than the Chandra flux (Oosterbroek et al. 2001; Wijnands 2002). Wijnands et al. (2003) concluded that during the Chandra observation the observed radiation was due to a hot crust and not associated with the core.

## 2 Observations, analysis, and results

Chandra observed MXB 1659–29 twice for $``$27 ksec: on October 15, 2002 (the 2002 observation), and on May 9, 2003 (the 2003 observation). We also used the $``$19 ksec observation performed on October 15–16, 2001 (the 2001 observation; Wijnands et al. 2003). During all observations the ACIS-S3 chip was used. The data were reduced and analyzed using CIAO 3.0. To make use of the latest calibration products, we reprocessed the 2001 observation. A minor background flare occurred during the 2003 observation (factor of $``$2; lasting $``$2 ksec). Its effect on the quality of the source data was negligible and we did not to remove this flare from the data. No flares occurred during the other observations.

For each observation, we extracted the number of source photons, the light curve, and the spectrum, using a circle with a radius of 3′′ as source extraction region and an annulus with an inner radius of $`7^{\prime \prime }`$ and an outer radius of 22′′ as background region. We detected 948$`\pm `$31, 263$`\pm `$16, and 107$`\pm `$10 counts (0.3–7 keV; background corrected) for the 2001, 2002, and 2003 observations, respectively, resulting in corresponding count rates of 0.050$`\pm `$0.002, 0.0097$`\pm `$0.0006, and 0.0039$`\pm `$0.0004 counts s-1. Wijnands et al. (2003) observed an eclipse and dipping behavior during the 2001 observation (similar to the outburst behavior of the source; Lewin 1979; Cominsky et al. 1983; Cominsky & Wood 1984, 1989). To search for eclipsing behavior during the 2002 and 2003 observations, we determined the orbital phase range covered by those observations using the time of the eclipse in the 2001 observation as the reference time. Given the orbital phase range traced during each observation, we expect to see a single eclipse per observation and, as anticipated, we did not detect any photons during the expected eclipse intervals. However, we also found that no photons were detected during several time intervals (of equal duration as the lengths of the eclipses) at different phases of the orbital period. Therefore, without prior knowledge of the eclipsing nature of MXB 1659–29, we could not have concluded that we saw eclipses during the 2002 and 2003 observations. Owing to the limited statistics of the 2002 and 2003 observations, no conclusions can be drawn about possible dipping behavior during these observations.

When extracting the spectra, we used all data, including those taken during the intervals of eclipses and possible dipping behavior. The eclipses could not be removed from the data before extracting the spectra because the uncertainties in the ephemeris presented by Oosterbroek et al. (2001) are sufficiently large so that the exact start and end times of the expected eclipses could not be determined. Instead we decreased the exposure time in the resulting spectral files by 900 s since the eclipse duration during outburst was found to be $``$900 s (Wachter et al. 2000) and Wijnands et al. (2003) reported an eclipse duration of 842$`\pm `$90 seconds for the 2001 observation. Small differences in the eclipse duration might be present between the observations but the expected effects on the resulting fluxes will be marginal. We also did not remove the data obtained during the dipping interval observed in the 2001 observation. Such dipping intervals are likely present during the other two observations but they cannot be identified in the light curves due to limited statistics. For those two observations all data had to be used and to obtain a homogeneous data selection across observations, we included the dipping interval observed during the 2001 observation. Wijnands et al. (2003) found evidence that this dipping behavior is likely due to a change in internal absorption in the system and not due to actual changes in the neutron-star properties. Therefore, the inclusion of the (possible) dipping intervals will likely result in a somewhat higher column density ($`N_\mathrm{H}`$) in the spectral fits than the true interstellar $`N_\mathrm{H}`$ toward the source, but should not significantly impact other source properties.

We grouped the spectra in bins of 15 counts to validate the use of the $`\chi ^2`$ fitting method and simultaneously fitted the three spectra using Xspec (Arnaud 1996). A variety of one-component models111E.g., a power-law model could fit the spectra but with an index of 4.7–5.8 suggesting soft thermal spectra. We also fitted a NSA plus power-law model to determine the upper-limits on the contribution of such a power-law tail to the 0.5–10 keV flux. Those limits are $`<`$20%–25%, $`<`$35%–45%, and $`<`$50%–100%, for the 2001, 2002, and 2003 observation, respectively. The range of upper limits is due to the range assumed in photon indices (between 1 and 2). could fit the individual spectra satisfactorily, but since we expect that the X-rays from MXB 1659–29 are due to the cooling of the neutron-star surface, for this paper we only fit the data using a neutron-star hydrogen atmosphere model (NSA; for weakly magnetized neutron stars; Zavlin et al. 1996). In such models the normalization is given by $`1/d^2`$, with $`d`$ in pc. The distance should be constant between observations and therefore we left the normalization tied among the different spectra (when leaving the normalizations free between observations, we find that they are consistent with each other). We expect the $`N_\mathrm{H}`$ toward the source to be very similar between observations (only minor variations are expected due to variable internal absorption) and this parameter was also tied. We assume a ’canonical’ neutron star with a radius of 10 km and a mass of 1.4 M.

From the fits, we found that the normalization was $`1.4_{0.8}^{+2.2}\times 10^8`$ which yields a source distance of 5–13 kpc. This is consistent with the distance range given in the literature (10–13 kpc; Oosterbroek et al. 2001; Muno et al. 2001). However, we found that the errors on the fit parameters were dominated by the large uncertainties in the normalization and did not allow us to realize the full potential of the data. If the source distance were established through an independent method, we could fix the normalization in the NSA models, resulting in considerably smaller errors on the remaining fit parameters. Therefore, instead of leaving the normalization as a free parameter, we fixed it so that it corresponded to a distance of 5, 10, and 13 kpc, covering the full range of allowed distances obtained when the normalization was a free parameter. To estimate the bolometric fluxes ($`F_{\mathrm{bol}}`$) we extrapolated the model to the energy range 0.01–100 keV which gives approximate bolometric fluxes222We verified that the 0.01–100 keV fluxes approximate $`F_{\mathrm{bol}}`$ by calculating the bolometric luminosity $`L_{\mathrm{bol}}=4\pi \sigma R_{\mathrm{}}^2T_{\mathrm{eff}}^\mathrm{}4`$, with $`\sigma `$ Stefan-Boltzmann constant, $`T_{\mathrm{eff}}^{\mathrm{}}`$ the effective temperature (at infinity), and $`R_{\mathrm{}}`$ the neutron-star radius (at infinity). The 0.01–100 keV fluxes were indeed consistent with the calculated $`F_{\mathrm{bol}}`$. We use the measured fluxes because their errors takes into account the uncertainties in $`N_\mathrm{H}`$ and the $`T_{\mathrm{eff}}^{\mathrm{}}`$ obtained for all observations. The $`L_{\mathrm{bol}}`$ errors are only calculated using the $`T_{\mathrm{eff}}^{\mathrm{}}`$ errors during one specific observation. . To calculate the flux errors, we fixed each free fit parameter (only one at a time) either to its minimum or maximum allowed value. After that we refitted the data and recalculated the fluxes. This process was repeated for each free parameter and the final flux range determined the flux errors. The fit parameters obtained are listed in Table 1.

This table shows that $`T_{\mathrm{eff}}^{\mathrm{}}`$ and $`F_{\mathrm{bol}}`$ decreased in time (Fig. 3). We fitted the $`T_{\mathrm{eff}}^{\mathrm{}}`$ and $`F_{\mathrm{bol}}`$ curves with an exponential decay function $`y(t)=c_0e^{\frac{tt_0}{\tau }}`$, with $`c_0`$ a normalization constant, $`t_0`$ the start time, and $`\tau `$ the $`e`$-folding time. We found that the other fit parameters were not very sensitive to the value of $`t_0`$, but when $`t_0`$ was left free it had adverse effects on the errors on those parameters. Therefore, we fixed $`t_0`$ to MJD 52159.5 which corresponds to midday September 7, 2001 (the last day MXB 1659–29 was found to be active; Wijnands et al. 2002a) and which can be regarded as an approximation of the time when $`T_{\mathrm{eff}}^{\mathrm{}}`$ and $`F_{\mathrm{bol}}`$ began to decrease. The assumed exponential functions could adequately describe the decrease in $`T_{\mathrm{eff}}^{\mathrm{}}`$ and $`F_{\mathrm{bol}}`$ (Fig. 3; alternative functions did not provide adequate fits). We found that $`\tau `$ and $`c_0`$ for the $`F_{\mathrm{bol}}`$ curve were 289$`\pm `$37, 262$`\pm `$33, and 254$`\pm `$29 days, and 70$`\pm `$9, 48$`\pm `$6, and 43$`\pm `$6 $`\times 10^{14}`$ erg cm-2 s-1, when assuming a distance of 5, 10, or 13 kpc, respectively, in the spectral fits. The corresponding $`\tau `$ and $`c_0`$ for the $`T_{\mathrm{eff}}^{\mathrm{}}`$ curve were 1153$`\pm `$160, 1060$`\pm `$126, and 1055$`\pm `$112 days, and 0.099$`\pm `$0.004, 0.126$`\pm `$0.004, and 0.139$`\pm `$0.004 keV. We saw no evidence that the curves approached a rock-bottom value: we found an upper limit on such a value of 3.5–7.5 $`\times 10^{14}`$ erg cm-2 s-1 for the $`F_{\mathrm{bol}}`$ curve (resulting in bolometric luminosity limits of 2.2–7.0 $`\times 10^{32}`$ erg s-1), and 0.06–0.07 keV for the $`T_{\mathrm{eff}}^{\mathrm{}}`$ curve.

## 3 Discussion

We have presented monitoring Chandra observations of MXB 1659–29 in quiescence. The first observation was taken only a month after the end of its last outburst which lasted 2.5 years; the second and third observations were taken $``$1 and $``$1.5 years after this initial one. Because it is expected that the emission should be dominated by thermal emission from the hot neutron-star crust (see Wijnands et al. 2003), we fitted the data with a NSA model for weakly ($`B<10^{89}`$ G) magnetized neutron stars. We found that $`F_{\mathrm{bol}}`$ decreased by a factor of $``$8 in $``$1.5 years and the rate of decrease followed an exponential decay function. Furthermore, $`T_{\mathrm{eff}}^{\mathrm{}}`$ also decreased and the rate of decrease again followed an exponential decay function. We found that the $`e`$-folding time of the $`T_{\mathrm{eff}}^{\mathrm{}}`$ curve was consistent with four times that of the $`F_{\mathrm{bol}}`$ curve, as expected if the emission is caused by a cooling black body for which the bolometric luminosity is given by $`L_{\mathrm{bol}}=4\pi \sigma R_{\mathrm{}}^2T_{\mathrm{eff}}^\mathrm{}4`$ (see footnote 2): if $`T_{\mathrm{eff}}^{\mathrm{}}`$ decays exponentially, $`L_{\mathrm{bol}}`$ (and thus $`F_{\mathrm{bol}}`$) will also decay exponentially but with an $`e`$-folding time four times smaller than that of $`T_{\mathrm{eff}}^{\mathrm{}}`$, exactly what we observe.

Our results support the suggestion that the crust was heated to high temperatures during the prolonged accretion event, which ended a month before our first observation, and that it is now cooling until it reaches thermal equilibrium with the core. Rutledge et al. (2002) calculated cooling curves for the neutron star in KS 1731–260, assuming different behaviors of the crustal micro-physics and the core cooling processes. Those curves can be used as a starting point to investigate how our results of MXB 1659–29 could be explained. Of those curves, only the one which assumes a large crustal conductivity and the presence of enhanced core cooling processes exhibits a large luminosity decrease in the first two years after the end of the last outburst, suggesting that the neutron star in MXB 1659–29 has similar properties. This conclusion was already tentatively reached by Wijnands et al. (2003) based on a comparison of the luminosity seen during the October 2001 Chandra observation with the significantly lower luminosity upper-limit found with ROSAT. But detailed cooling curves for the neutron star in MXB 1659–29 need to be calculated to fully explore (and exploit) the impact of our observations on our understanding of the structure of neutron stars. The cooling curves calculated by Rutledge et al. (2002) for KS 1731–260 only give us a hint of the behavior of MXB 1659–29 because they depend on the long-term ($`>10^4`$ years) accretion history of the source. For KS 1731–260, this long-term accretion behavior was quite unconstrained due to large uncertainties in the averaged duration of the outbursts, the time-averaged accretion rate during the outbursts, and the time the source spent in quiescence. However, the accretion history of MXB 1659–29 over the last three decades is much better constrained (Wijnands et al. 2003), which will help to reduce the uncertainties in its long-term averaged accretion history allowing for more detailed cooling curves to be calculated for MXB 1659–29. This might help to constrain the physics of the crust better for MXB 1659–29 than for KS 1731–260. The only significant uncertainty left is that of the source distance; however, we found that this only affects the exact values of the bolometric fluxes and the effective temperatures, but not their rate of decay.

Our 0.5–10 keV flux during the May 2003 Chandra observation is still higher than the upper limit found with ROSAT, suggesting that the crust will cool even further in quiescence and that we have not yet reached thermal equilibrium between the crust and core. Further monitoring observations are needed to follow the cooling curve of the crust to determine the moment when the crust is thermally relaxed again. When this occurs, no significant further decrease of the quiescent luminosity is expected and from this bottom level the state of the core can be inferred. As of yet, we have found no evidence that the flux and temperature are reaching a leveling-off value, associated with the temperature of the core, although the limits we obtained are not very stringent.

Jonker, Wijnands, & van der Klis (2004) suggested that the difference in luminosity of MXB 1659–29 between the ROSAT non-detection and the 2001 Chandra observation might be due to differences in residual accretion rate onto the surface. Residual accretion could indeed produce soft spectra (e.g., Zampieri et al. 1995), but to explain the exponential decay we observe for $`F_{\mathrm{bol}}`$ and $`T_{\mathrm{eff}}^{\mathrm{}}`$, the residual accretion rate must also decrease exponentially with a timescale of a year. Although this cannot be completely ruled out, we believe this is unlikely since other neutron-star transients have been observed to reach their quiescent states on timescales of only tens to several tens of days at the end of their outbursts (e.g., Campana et al. 1998b; Jonker et al. 2003) and the variations in accretion rate tend to be more stochastic. Moreover, if the neutron star has a significant magnetic field strength, this might inhibit material from reaching the surface when accreting at the inferred low rates.


text:
# Gravitational lensing as a contaminant of the gravity wave signal in CMB

## I Introduction

CMB polarization is generated by Thomson scattering of photons off free electrons. To generate polarization one needs an anisotropic distribution of photons in the electron rest frame (more specifically, a non-vanishing quadrupole moment) and scattering that couples this angular anisotropy to the polarization (Thomson scattering in this case). These two conditions are satisfied during recombination, when most of the polarization signal is generated. After recombination the free streaming of photons leads to a large quadrupole moment, so if some fraction of the photons is rescattered when the universe is reionized then a new polarization contribution will be generated at the angular scale of horizon at reionization. One often refers to the two contributions as the recombination and reionization components, respectively.

Thomson scattering generates linear polarization only. This is usually expressed in terms of Stokes parameters $`Q`$ and $`U`$, which are coordinate dependent. They can be decomposed into coordinate independent $`E`$ and $`B`$ type polarizations Seljak (1997); Zaldarriaga and Seljak (1997); Kamionkowski et al. (1997a) with opposite parities. To linear order in perturbation theory, primordial scalar (density) perturbations can only generate $`E`$ polarization, while gravitational waves (GWs) can generate both scalar $`E`$ and pseudoscalar $`B`$. If the amplitude of the gravity waves is very small relative to scalars it cannot be isolated from the temperature anisotropies or $`E`$ polarization due to cosmic variance. The $`B`$ polarization is however insensitive to cosmic variance from scalar modes and is limited only by instrument noise, foregrounds and sky coverage. This fact generated attention as a potentially promising tool to detect gravity waves and test inflation Seljak and Zaldarriaga (1997); Kamionkowski et al. (1997b). The amplitude of gravity waves produced during inflation depends on its energy scale: higher energy scales give larger amplitude of gravity waves. In terms of the tensor to scalar power spectrum ratio one has $`r=P_h/P_RV_{}^4`$, where $`P_h`$ is the tensor power spectrum, $`P_R`$ is the scalar curvature power spectrum and $`V_{}^4`$ is the energy density during inflation when the present Hubble scale exited the horizon. $`V_{}`$ has units of energy and has been termed the “energy scale” of inflation. We do not know this energy scale, but one of the possibilities is the scale of grand unification theories (GUTs) at $`V_{}10^{16}`$ GeV. At this energy scale the gravity wave contribution is sufficiently large to be detectable. In this paper we define the tensor to scalar ratio in terms of their primordial power spectra, rather than the quadrupole moments as often defined. This has the advantage of relating the tensor to scalar ratio directly to the inflationary predictions independent of the cosmological parameters, which affect the CMB anisotropy spectra. For typical parameters we find $`C_2^T/C_2^Sr/2`$.

A future satellite mission dedicated to $`B`$ type polarization has been identified as one of the NASA Einstein probes to be built over the next decade. One of the outstanding questions regarding such a mission is what are the required angular resolution and sensitivity to maximize the science output and at what level do systematics swamp the improvements in these. It has been pointed out Zaldarriaga and Seljak (1998) that gravitational lensing leads to a generation of $`B`$ polarization even if none was present in the early universe. This could limit the extraction of gravity wave signal if unaccounted for Hu (2002); Lewis et al. (2002). One can try to reconstruct the gravitational lensing potential using the non-Gaussian information present in the CMB data to improve the limits (the large-scale lensing $`B`$-modes exhibit higher-order correlations with small-scale polarization whereas inflationary GW $`B`$-modes do not). Recent work applied quadratic estimators Hu and Okamoto (2002) to argue that using these estimators leads to an order of magnitude improvement for the no noise experiment Kesden et al. (2002); Knox and Song (2002). This work has been interpreted as providing a fundamental limit to the gravity wave extraction due to the lensing. However, it is important to note that these papers do not rule out the possibility that better reconstruction methods may be constructed. Indeed, in a recent work we have shown that better estimators are indeed possible Hirata and Seljak (2003a). We have demostrated explicitly that one can improve upon the no noise limit of the quadratic estimator. Indeed, in the idealized case of no noise, perfect resolution and lensing by a single scalar deflection potential, the lensing reconstruction can be achieved exactly and the lensing contamination can be removed completely. It is easy to see why this is so: lensing displaces photons, so one can write the final $`Q`$ and $`U`$ polarization in terms of the initial $`Q`$ and $`U`$ and lensing deflection angle $`𝐝`$. In the absence of gravity waves (null hypothesis) initial $`B`$ vanishes and ignoring lensing rotation (i.e. taking $`\times 𝐝=\mathrm{𝟎}`$) the deflection angle can be written in terms of a scalar field with one degree of freedom at each point, $`𝐝=𝚽`$. In this case for $`N`$ pixels there are 2$`N`$ observables (two at each point, $`Q`$ and $`U`$ or $`E`$ and $`B`$), and 2$`N`$ unknowns, initial $`E`$ and $`\mathrm{\Phi }`$. The number of unknowns thus equals the number of equations, so one can solve it exactly in the absence of noise. It is of course possible that there are degenerate modes that cannot be reconstructed; however, it has been shown that the fraction of modes that are degenerate is small and may even be zero. (See Appendix B of Ref. Hirata and Seljak (2003a) for details.)

The idealized case discussed above is unrealistic, since in the real world noise always limits the achievable sensitivity (in Ref. Hirata and Seljak (2003a) we found that the lensing rotation is not a limitation at the sensitivity levels that can be achieved in foreseeable future). At the same time, if the detector noise is large gravitational lensing is not limiting the GW detection anyways. As the detector noise is lowered our ability to clean the lensing contamination improves as well and if the scaling between the detector and residual lensing noise is linear then lensing may never be the dominating source of noise. The relevant question regarding the lensing contamination is thus not whether it provides a fundamental limit (which remains an open question), but rather how much does it degrade the gravity wave sensitivity for a given instrument noise and angular resolution. We will phrase this question in terms of a $`B`$-mode noise power spectrum: the minimum detectable tensor-to-scalar power spectrum ratio $`r`$ that can be observed then scales linearly with it.

The lensing degradation issue is particularly interesting in the context of the required noise and angular resolution of a future CMB polarization satellite. There are important instrument and cost tradeoffs that need to be included in the design of such a mission. For example, since the bulk of the gravity wave signal is at large scales one could devise a high sensitivity low angular resolution instrument costing significantly less than the equivalent high angular resolution mission. However, in this case one would not be able to reconstruct the lensing potential, making the lensing contamination more significant. The goal of this paper is to provide some guidance to these considerations, obtaining the lensing degradation factors as a function of detector noise and resolution. We will use both quadratic method Hu and Okamoto (2002) and the maximum-likelihood method Hirata and Seljak (2003a), which improves upon quadratic in the low noise regime. We will ignore other sources of contamination such as foregrounds and instrument-related issues, which should be included in the full consideration of pros and cons of a mission design. Unless otherwise specified, we will use the fiducial cosmology of Ref. Hirata and Seljak (2003a).

## II Lensing and $`B`$-modes

In this section, we briefly review techniques for lensing reconstruction from CMB polarization data. We then arrive at the main objective of this paper: to calculate, via simulations, the remaining “noise” level (including lensing residuals) after the lensing reconstruction has been performed and used to remove lensing $`B`$-modes. We compute $`w_{P,eff}^{1/2}`$, the white noise power spectrum of the combined instrumental $`B`$-mode noise and lensing residuals that remains after lens cleaning.

### II.1 Lensing reconstruction

While unique features of weak lensing effect on CMB power spectrum, such as smoothing of the peaks and transfer of power to small scales Seljak (1996), can be used to deduce its presence, it is the nongaussian signatures that allow for the lensing reconstruction Bernardeau (1997). There are several lensing reconstruction methods proposed in the literature. Simple quadratic estimators Zaldarriaga and Seljak (1999) have been shown to be near optimal for the Wilkinson Microwave Anisotropy Probe (WMAP) and possibly Planck Zaldarriaga (2000), but more efficient quadratic estimators have been shown to improve upon these for low noise, high angular resolution experiments Hu (2001); Hu and Okamoto (2002). Recently we used likelihood techniques to develop an iterative estimator Hirata and Seljak (2003b, a), which further improves upon these if detector noise is sufficiently low, specially for the polarization data. The corresponding reconstruction errors have been computed as a function of noise and angular resolution in the absence of gravity wave contribution Hirata and Seljak (2003a).

We use only polarization information (not temperature) for our lensing reconstruction in order to reduce the computation time required for the simulations. This is obviously a conservative assumption. But note that – although temperature anisotropies are potentially useful for measuring the convergence power spectrum, and for mapping the convergence on large angular scales ($`l200`$) – most of the lensing-induced $`B`$-modes come from smaller-scale convergence modes that are only accessible using polarization data Hu and Okamoto (2002). Additionally, the small-scale temperature anisotropies are contaminated by secondary processes such as the kinetic Sunyaev-Zeldovich/Ostriker-Vishniac effect and scattering during inhomogeneous reionization, which have the same spectral dependence as primary CMB fluctuations and hence can significantly degrade lensing reconstruction from temperature Santos et al. (2003). By comparison, CMB polarization is expected to be essentially free of secondary scattering contamination Hu (2000); Santos et al. (2003).

Gravitational lensing re-maps the primary CMB polarization field according to:

$$\stackrel{~}{P}(\widehat{𝐧})=P\left(\widehat{𝐧}2^2\kappa (\widehat{𝐧})\right),$$
(1)
where $`P`$ is the primary polarization, $`\stackrel{~}{P}(\widehat{𝐧})`$ is the lensed polarization in direction $`\widehat{𝐧}`$, and $`\kappa `$ is the lensing convergence field. The lensing has two effects on polarization that are of interest here. One is the generation of $`B`$-modes in the lensed polarization field; the other is the generation of an anisotropic two-point function:

$$E_{\stackrel{}{l}_1}B_{\stackrel{}{l}_2}=\frac{1}{\pi ^{1/2}l^2}C_{l_1}^{EE}\stackrel{}{l}\stackrel{}{l}_1\mathrm{sin}(2\alpha )\kappa _\stackrel{}{l},$$
(2)
where $`\stackrel{}{l}=\stackrel{}{l}_1+\stackrel{}{l}_2`$, $`\alpha `$ is the angle between $`\stackrel{}{l}_1`$ and $`\stackrel{}{l}_2`$, and higher-order terms in $`\kappa `$ have been neglected. (The average is taken over CMB realizations only, not including noise.) The $`EB`$ quadratic estimator developed by Ref. Hu and Okamoto (2002) is obtained by taking a minimum-variance linear combination of $`EB`$ products subject to the constraint that the expectation value is $`\kappa _\stackrel{}{l}`$:

$$\widehat{\kappa }_\stackrel{}{l}^{EB,(\mathrm{quad})}=\frac{A_{EB}(l)}{2}\underset{\stackrel{}{l}_1}{}\frac{C_{l_1}^{EE}E_{\stackrel{}{l}_1}B_{\stackrel{}{l}_2}}{(\stackrel{~}{C}_{l_1}^{EE}+N_{l_1}^{EE})(\stackrel{~}{C}_{l_2}^{BB}+N_{l_2}^{BB})},$$
(3)
where $`A_{EB}`$ is a normalization factor, $`N_{l_1}^{EE}`$ is the $`E`$-mode polarization noise power spectrum, and $`\stackrel{}{l}_2=\stackrel{}{l}\stackrel{}{l}_1`$. Ref. Hu and Okamoto (2002) gives similar quadratic estimators using the $`EE`$ products; in this paper we use an $`EE+EB`$ quadratic estimator that is a minimum-variance weighting of the $`EE`$ and $`EB`$ products (although almost all of the information for low-noise experiments comes from $`EB`$). For the purpose of cleaning out the lensing $`B`$-mode, it is best to Wiener-filter the quadratic estimator Knox and Song (2002) since this minimizes the mean squared difference between the estimated and “true” convergence maps:

$$\widehat{\kappa }_\stackrel{}{l}^{WF}=\frac{C_l^{\kappa \kappa }}{C_l^{\kappa \kappa }+\sigma _l^{\kappa \kappa }}\widehat{\kappa }_\stackrel{}{l}^{(\mathrm{quad})},$$
(4)
where $`\sigma _l^{\kappa \kappa }`$ is the power spectrum of the noise in the estimator $`\widehat{\kappa }_\stackrel{}{l}^{(\mathrm{quad})}`$. Lens cleaning consists of applying a de-lensing operation that inverts the mapping of Eq. (1).

The quadratic estimator using polarization has excellent performance, in particular it can ultimately recover the convergence map with signal-to-noise ratio exceeding unity out to $`l1000`$ if the instrument noise is sufficiently low. Nevertheless, as is shown in Ref. Hu and Okamoto (2002), the accuracy of quadratic lens reconstruction is fundamentally limited by the $`B`$-mode cosmic variance. Put another way, the estimation of one convergence mode $`\kappa _\stackrel{}{l}`$ is contaminated by $`B`$-mode power generated by the other covergence modes $`\kappa _\stackrel{}{l}^{}`$ Kesden et al. (2003).

The iterative reconstruction algorithm Hirata and Seljak (2003a) avoids this problem. While Ref. Hirata and Seljak (2003a) introduces the iterative lensing reconstruction using the likelihood function, the algorithm used can also be thought of as follows: for low-noise experiments, the dominant source of uncertainty in the lens reconstruction is the above-mentioned cross-talk among different convergence modes. Once an estimate of the lensing field is available using Eq. (4), the CMB can be de-lensed using the estimated convergence field. Since the de-lensed CMB map has most of the lensing-induced $`B`$-mode removed, the cross-talk among the convergence modes is reduced; therefore, re-application of the quadratic estimator (with different weighting) results in an even lower-noise convergence map. The reader is referred to Ref. Hirata and Seljak (2003a) for implementation details. The iterative procedure is only useful if the instrument noise is below the lensing-induced $`B`$ signal (5.3 $`\mu \mathrm{K}\mathrm{arcmin}`$ in the fiducial model at low $`l`$). However, as the instrument noise is reduced, the uncertainty in the iterative lensing reconstruction is also reduced; indeed, in the absence of foregrounds and the rotational component of the deflection field, it is not known whether there is any fundamental limitation to the reconstruction accuracy that can be obtained via the iterative estimator.

Since both the quadratic and iterative estimators assume absence of gravity waves one must be careful when applying them to the case with gravity waves. As discussed in the introduction, in the absence of noise one is simply solving for lensing assuming $`B=0`$, which is of course violated if gravity waves are present. The most general approach is to solve for the tensor $`B`$ modes, combined scalar and tensor $`E`$ modes and lensing potential $`\mathrm{\Phi }`$ simultaneously. However, since most of the lensing signal is at high $`l`$, while most of the intrinsic gravity wave signal is at low $`l`$, we can simplify and ignore the information from low $`l`$ in the lensing reconstruction, but only use low $`l`$ in $`B`$ mode cleaning. We have tried cutoffs of $`l=50`$ and $`l=150`$ as well as no cutoff and found no significant difference between them. This is not surprising, since number of modes not used in the reconstruction is in both cases small compared to the total number of modes. Note that $`l=150`$ corresponds roughly to the scale below which the gravity waves become negligible, while the reionization peak only contributes for $`l50`$.

While analytical (Fisher matrix) expressions for noise spectra exist for both quadratic and iterative methods, they may not be fully reliable in either case. The Fisher matrix tends to significantly underestimate the noise in the iterative method because the lensing reconstruction error is non-Gaussian and realization-dependent Hirata and Seljak (2003a). The analytical estimate for the quadratic estimator error estimate based on the approximation of Gaussian $`E`$ and $`B`$ modes Hu and Okamoto (2002) is more reliable, but still does not include all terms in the covariance matrix. Some of these are analytically tractable and have been shown to increase the covariance by up to 20% Cooray and Kesden (2003); Kesden et al. (2003), but other terms (specifically higher order terms in the deflection angle $`𝐝`$) have not yet been computed. The safest approach is thus to use numerical Monte Carlo simulations, which by construction include all of the terms.

### II.2 Simulations

We compute the post-cleaning $`B`$-mode noise power spectrum via simulation as follows. First a simulated pure $`E`$ Gaussian CMB polarization field is generated using the unlensed $`C_l^{EE}`$ power spectrum from cmbfast Seljak and Zaldarriaga (1996) for the fiducial $`\mathrm{\Lambda }`$CDM cosmology of Ref. Hirata and Seljak (2003a). Then a Gaussian convergence field $`\kappa `$ is generated, and the polarization $`Q`$ and $`U`$ fields are re-mapped according to Eq. (1). We next add noise to the polarization field with power spectrum:

$$C_l^{EE}(\mathrm{noise})=C_l^{BB}(\mathrm{noise})=w_P^1\mathrm{exp}\frac{l(l+1)\theta _{FWHM}^2}{8\mathrm{ln}2},$$
(5)
where $`w_P^{1/2}`$ is the instrument noise (typically measured in $`\mu \mathrm{K}\mathrm{arcmin}`$) and $`\theta _{FWHM}`$ is the full width at half maximum of the instrument’s beam. The quadratic and iterative lens reconstruction algorithms are then applied as described in Ref. Hirata and Seljak (2003a), with the modification that we remove $`l<150`$ modes (technically we have set all rows and columns corresponding to $`l<150`$ modes to zero in the $`\sigma `$-matrices of Ref. Hirata and Seljak (2003a)). Once the estimated convergence field $`\widehat{\kappa }`$ has been determined, it is used to “de-lens” the CMB polarization map, thereby yielding a map of the primary polarization field. We then compute the $`B`$-mode power by averaging $`(B_𝐥^{(\mathrm{res})})^2`$ over the $`l<150`$ modes in the simulation. Each simulation is run on a $`2048\times 2048`$ square grid with periodic boundary conditions and grid spacing of 1 arcmin, corresponding to a total area of 1165 square degrees. The residual $`B`$-mode power spectra quoted here are determined by averaging over 4 such simulations. We define the “effective noise” by:

$$w_{P,eff}^1=C_l^{BB}(\mathrm{residual});$$
(6)
this has units of $`(\mu \mathrm{K}\mathrm{arcmin})^2`$ and represents the white noise from combined instrument noise and lensing residuals that limits detectability of the GW signal. Fig. 1 shows an input $`B`$ polarization map assuming $`C_2^T/C_2^S=0.012`$ (upper left), no cleaning map (upper right), quadratic cleaning map (lower right) and iterative cleaning map (lower right).

### II.3 Results

The results for the error power spectrum of such a reconstruction are shown in Table 1 for a variety of noise and beam levels, for both quadratic and iterative methods. Since the spectra are close to white noise for $`l<150`$, we only show the amplitudes and not the full spectrum. Note that the noise amplitude shown is the total noise and contains both lensing and detector noise contributions. Also note that for the 20’ beam (and for larger beams), cleaning can actually make the $`B`$-mode noise worse because the de-lensing operation transfers power from high-$`l`$, unresolved CMB modes down to low $`l`$. In principle this problem can be circumvented by Wiener-filtering the CMB prior to the de-lensing operation; we have not implemented this because lens cleaning is not useful for such wide beams anyway.

The results can be divided into high, intermediate and low noise regimes. For high detector noise, $`w_P^{1/2}>5\mu \mathrm{K}\mathrm{arcmin}`$, lensing is a minor contributor to the total noise. Note that this noise level is still a factor of 100 (in power) lower than expected Planck polarization noise, so clearly any discussion of lensing induced noise in $`B`$ is relevant only for a post-Planck CMB mission dedicated to polarization. For $`w_P^{1/2}=6\mu \mathrm{K}\mathrm{arcmin}`$ a $`10^{}`$ beam results in a total rms noise of $`7.5\mu \mathrm{K}\mathrm{arcmin}`$ and $`7^{}`$ beam in $`7.3\mu \mathrm{K}\mathrm{arcmin}`$ for either method. Without cleaning the combined lensing and detector noise would be $`8\mu \mathrm{K}\mathrm{arcmin}`$, so cleaning hardly improves anything at all. Improvements appear when the detector noise drops below the lensing noise, which for our model is $`5.3\mu \mathrm{K}\mathrm{arcmin}`$. In the intermediate range ($`2`$$`5\mu \mathrm{K}\mathrm{arcmin}`$) the quadratic estimator is very similar to the iterative method in terms of the residual noise. For example, for $`w_P^{1/2}=3\mu \mathrm{K}\mathrm{arcmin}`$ the residual noise is $`4.5\mu \mathrm{K}\mathrm{arcmin}`$ for $`7^{}`$ beam, a factor of 2 in power greater than the detector noise alone and a factor of 2 lower than no lens cleaning/large beam case. Going to a $`4^{}`$ beam marginally improves upon this.

Finally, in the low noise regime the iterative method clearly outperforms quadratic method. The quadratic method bottoms out roughly at $`w_{P,eff}^{1/2}=2\mu \mathrm{K}\mathrm{arcmin}`$, which is a factor of 7 improvement over the no-cleaning lensing noise of $`w_{P,eff}^{1/2}=5.3\mu \mathrm{K}\mathrm{arcmin}`$. This bottoming out of the quadratic method has lead to the suggestion that lensing noise limit may be fundamental and cannot be improved upon Knox and Song (2002); Kesden et al. (2002). However, this conclusion is only valid for the quadratic method, which is not the optimal method in the low detector noise regime. We find that the iterative method always reduces the overall noise as the detector noise is decreased, at least over the range tested with our simulations (which should cover the range of interest for the next generation CMB satellite dedicated to polarization). At the lowest noise and smallest beam ($`0.25\mu \mathrm{K}\mathrm{arcmin}`$, $`2^{}`$) tested in our simulation the lensing noise is reduced by more than a factor of 40. Further improvements are likely if the detector noise is reduced below $`0.25\mu \mathrm{K}\mathrm{arcmin}`$, but the iterative method becomes very computationally expensive and we have not explored these very low detector noise cases here.

For a given noise and lensing induced $`B`$-mode power spectrum the resulting uncertainty on $`r`$ is:

$$\sigma _r^2=f_{sky}w_{P,eff}^2\underset{l}{}\frac{2l+1}{2}\left(\frac{C_{Tl}^{BB}}{r}\right)^2,$$
(7)
where $`C_{Tl}^{BB}`$ is the GW power spectrum of $`B`$ modes (Fig. 2), and $`w_P`$ is the inverse noise variance per solid angle per polarization. Since $`r`$ is merely supplying the normalization of the tensor power spectrum, $`C_{Tl}^{BB}/r`$ is fixed by the background cosmology. We are interested in large scales only, so the noise spectrum has been approximated as a constant and taken out of the sum. It is easy to see from this expression that the limit on $`r=T/S`$ is proportional to the noise weight per solid angle $`w_{P,eff}^1`$. Therefore, the noise degradation factors are the same as $`r`$ degradation factors. For partial-sky coverage, Eq. (7) must be modified to take into account sky cuts; while $`\sigma _rf_{sky}^{1/2}`$ for the recombination peak on degree scales, the reionization peak present at $`l<20`$ exhibits a much more complicated dependence on the survey geometry due to cross-leakage of $`E`$ and $`B`$ modes induced by, e.g. the Galactic Plane cut Lewis et al. (2002). Nevertheless, it remains true even in the presence of sky cuts that $`\sigma _rw_{P,eff}^1`$ if only the pure $`B`$-modes are used for GW searches. (But note that for the reionization peak, the uncertainty on $`r`$ can be non-Gaussian due to the small number of modes, which complicates hypothesis tests for a GW contribution Lewis et al. (2002).)

## III Discussion and Conclusions

While the degradation factors can be computed independently of the theoretical spectrum, the actual achievable values of $`r`$ depend on it. It is worth considering the reionization and recombination peaks separately. While the reionization peak gives typically higher signal to noise if $`f_{sky}=1`$, it depends sensitively on the Thomson scattering optical depth $`\tau `$ due to reionization, which is still rather uncertain (although this may improve with future observations of the reionization-induced $`TE`$ correlation, which was recently detected by WMAP Bennett et al. (2003); Kogut et al. (2003)). In addition, incomplete sky coverage and foregrounds are particularly worrisome on large scales, so the reionization peak may be more difficult to observe than the recombination peak at $`l100`$.

As an example, for $`l<20`$ and optical depth $`\tau =0.17`$ one finds $`\sigma _r=5\times 10^7`$ for $`f_{\mathrm{sky}}=1`$ and $`w_{P,eff}^{1/2}=0.8\mu \mathrm{K}\mathrm{arcmin}`$, which is the lowest noise level found in our simulations. The energy scale of inflation scales as $`r^{1/4}`$, so in this case the minimum energy scale one can detect at $`3\sigma `$ is $`V_{}10^{15}`$ GeV. For full-sky coverage, the signal-to-noise scales somewhat less rapidly than $`\tau ^2`$ and reducing the optical depth to $`\tau =0.07`$ increases $`\sigma _r`$ by a factor of 3. Additional reduction will be caused by incomplete sky coverage, which will cause leakage of $`E`$ into $`B`$ Lewis et al. (2002); Bunn et al. (2003); this degradation will be worst for models with late reionization because this pushes the $`B`$ reionization peak to the low multipoles where sky-cut effects are most severe. For the recombination peak at $`l>20`$ one has $`\sigma _r=10^5`$ for $`w_{P,eff}^{1/2}=0.8\mu \mathrm{K}\mathrm{arcmin}`$, which is a factor of 20 worse than reionization peak if $`\tau =0.17`$ and a factor of 7 worse if $`\tau =0.07`$. The corresponding energy scale that can be detected at $`3\sigma `$ is $`V_{}2.3\times 10^{15}`$ GeV.

To summarize, in this paper we present a detailed numerical study of how well can one clean $`B`$-type polarization of the contamination caused by gravitational lensing. We find that quadratic methods are able to reduce the noise power by up to a factor of 7, while our iterative method is able to reduce the noise significantly beyond that, at least a factor of 40 in our simulations. With this method we do not find any fundamental lower limit caused by gravitational lensing, in the sense that over the range of parameters considered, reducing the instrument noise always leads to a reduction in lensing noise as well. However, the scaling is sublinear, so at low detector noise levels lensing noise dominates over instrument noise. With the noise cleaning one can achieve tensor to scalar ratios as low as $`r10^6`$ and possibly even lower, which should allow us to differentiate between different models of structure formation with high precision. In particular, inflationary models with an energy scale significantly below $`10^{15}`$ GeV, as well as cyclic/ekpyrotic models Khoury et al. (2001); Steinhardt and Turok (2002); Boyle et al. (2003) predict that the gravity waves should be negligible and so their predictions could be falsifiable with the future CMB polarization studies.

###### Acknowledgements.

U.S. is supported by Packard Foundation, Sloan Foundation, NASA NAG5-1993 and NSF CAREER-0132953. C.H. is supported by a fellowship from the NASA Graduate Student Researchers Program (GSRP).