text:
# Microlensing neutron stars

## 1 Introduction

The determination of neutron star masses so far has been limited to binary systems. Particularly, the masses of radio pulsars in binary systems lie in an impressively narrow range $`M_{\mathrm{PSR}}=1.35\pm 0.04M_{\mathrm{}}`$ (Thorsett et al. T93 (1993); Thorsett & Chakrabarty T99 (1999)). The spread of masses based on the measurement of X-ray binary systems is somewhat larger (Nagase N89 (1989); Reynolds, Bell & Hilditch R92 (1992); Heap & Corcoran H92 (1992); Reynolds et al. R97 (1997); Ash et al. A99 (1999)), but the measurements are consistent with the mass value quoted above. However, it is not clear whether these measurements are representative of all neutron stars. In this paper we investigate the chances of estimating the mass of neutron stars by gravitational microlensing.

The galaxy hosts $`10^9`$ neutron stars (Timmes, Woosley & Weaver TWW96 (1996)), of which $`10^5`$ are pulsars (Narayan Narayan (1987)). Thus the fraction of neutron stars which are also pulsars is only $`10^310^4`$. Moreover, only a small fraction of pulsars appears to be in binary systems (of $`10^3`$ known radio pulsars, only $`10^2`$ are known to be in binary systems). The range of observed masses is surprisingly narrow in view of the theoretical stability considerations, which allow neutron star masses in the range $`0.1M_{\mathrm{}}`$ to $`3M_{\mathrm{}}`$ (see, e.g., Weber Weber (1999)). The predicted initial mass function from supernova simulations provides a more narrow mass range within $`1M_{\mathrm{}}`$ to $`2M_{\mathrm{}}`$, depending on the supernova type (Timmes, Woosley & Weaver TWW96 (1996)).

The neutron star mass function is of great importance not only for astrophysical purposes, but also for fundamental physics. On the one hand, neutron stars are the only known place in nature where theoretical ideas on the phase diagram of quantum chromodynamics at high density may be tested \[for recent reviews see Schwarz (S98 (1998)), Weber (Weber (1999)) and Alford, Bowers & Rajagopal (ABR00 (2000))\]. On the other hand, the equation of state of neutron stars is important to predict the wave pattern of gravitational radiation emitted in mergers of neutron star binaries, one of the most important sources expected for gravitational wave interferometers like LIGO111 Laser Interferometer Gravitational Wave Observatory http://www.ligo.caltech.edu/ and VIRGO222http://www.virgo.infn.it/.

Gravitational lensing \[see, e.g., Schneider, Ehlers & Falco (Schneider (1992))\] by now has been proven to be a very powerful technique to study mass distributions. For dark or low-luminosity stellar mass objects, gravitational microlensing \[for reviews see Paczyński (Paczynski (1996)) and Zakharov & Sazhin (ZS98 (1998))\] provides a tool to measure their masses. Gravitational microlensing towards the galactic bulge and towards the galactic spiral arms has been studied by various groups (OGLE333The Optical Gravitational Lens Experiment http:// astro.princeton.edu/~ogle/, Udalski et al. Udalski1 (1994), Udalski2 (2000), Wozniak et al. Wozniak (2001); MACHO444http://wwwmacho.mcmaster.ca/, Alcock et al. Alcock1 (1997), Alcock2a (1999), Alcock2 (2000), Popowski et al. Popowski (2001); EROS555http://eros.in2p3.fr/ Derue et al. Derue1 (1999), Derue2 (2001)). The optical depth of the bulge is $`\tau =(1.4\pm 0.3)\times 10^6`$ (Popowski et al. Popowski (2001)) and drops to $`\tau =0.45_{0.11}^{+0.24}\times 10^6`$ (Derue et al. Derue2 (2001)) in the spiral arms. About $`800`$ candidates for gravitational microlensing have been observed.

In this work we show that gravitational microlensing of neutron stars in the galaxy can in principle be observed with present-day techniques (photometric lensing). We find that the rate of observing a neutron star microlensing event is about $`10^7`$ per year per observed light source in the direction of the galactic bulge; the optical depth is $`10^8`$. Correspondingly, if we assume that the pulsar number density is proportional to the neutron star number density, the rate of observing a pulsar as a microlens is lower by a factor $`10^310^4`$. As the number of monitored stars is $`10^7`$ in each experiment, we expect that no such observation has been made yet. A data comparison between the ATNF Pulsar Catalogue666http://wwwatnf.atnf.csiro.au/research/pulsar/ catalogue and known microlensing events towards the galactic bulge and towards the galactic spiral arms (OGLE, Wozniak et al. Wozniak (2001); MACHO, Alcock et al. Alcock1 (1997), Alcock2a (1999), Alcock2 (2000); EROS, Derue et al. Derue1 (1999), Derue2 (2001)) shows that indeed none of the known pulsars has been acting as an observed microlens so far.

Once a microlensing event of appropriate duration (a few days to a few months) has been found the challenge is to figure out whether a neutron star or a different compact object has been acting as lens. This task may be accomplished by means of astrometric microlensing (Høg, Novikov & Polnarev HNG (1995), Miyamoto & Yoshi MY (1995), Walker W (1995), Miralda-Escudé M (1996) and Boden, Shao & Van Buren BSB (1997)). Astrometric lensing measures the shift of the centre of light of the source star. While photometric lensing is most efficient if the source passes through the Einstein ring, astrometric lensing is most efficient if the source passes close to, but outside, the Einstein ring. The most important background will be main sequence stars that are not properly resolved (Gould G00 (2000)). White or brown dwarfs as well as black holes are candidates for similar microlensing events due to their similar masses. As has been shown by Gould (G00 (2000)), upcoming astrometric space missions (like GAIA777http://astro.estec.esa.nl/GAIA/ and SIM888Space Interferometry Mission http://sim.jpl.nasa.gov/) combined with ground-based observations should be able to pick out a handful of neutron stars. The observation of astrometric lensing will also be possible with the help of long-baseline optical interferometry, and has been studied in some detail for the ESO Very Large Telescope Interferometer (VLTI999http://www.eso.org/projects/vlti/) by Delplancke et al. (DGR (2001)). The ESO VLTI offers yet another new possibility; it might be possible for the first time to resolve the double images of microlensed stars. For events that are long enough to observe the parallax effect of the light-curve, a measurement of the lens mass will be possible (Delplancke et al. DGR (2001)).

In order to improve the lensing probability for neutron stars we propose a new search strategy that looks for the microlensing event of a known nearby high proper motion pulsar (sweeping a large angle on the sky). In that case the distance and velocity are constrained from the dispersion measure or via parallax measurement and a measurement of the pulsar mass should be possible. Thus gravitational microlensing might enable us to determine the mass of stand-alone pulsars. The cross section of astrometric lensing is larger than that of photometric lensing, which additionally enhances the probability of detection. Thus astrometric measurements around high velocity proper motion pulsars in dense fields might allow a mass determination of pulsars. The potential of forthcoming astrometric space missions for astrometric microlensing of ordinary stars has been investigated in detail by Paczyński (P98 (1998)), Gould (G00 (2000)) and Salim & Gould (SG (2000)) for SIM and by Belokurov & Evans (BE (2001)) for GAIA.

In section 2 we use two pulsar population models of our galaxy to estimate the lensing rate and the optical depth towards different directions. In section 3 we describe the catalogue comparison of pulsars and gravitational microlensing events. We suggest a search strategy to increase the probability to find a neutron star that acts as a gravitational microlense in section 4 and end with a short conclusion.

## 2 Lensing rate and optical depth

Let us first investigate the possibility that a neutron star is detected by photometric microlensing.

### 2.1 Lensing rate

A simple estimate for the rate $`p`$ of gravitational microlensing events per observed background star is obtained easily. We consider a neutron star at distance $`r`$ that passes by the line of sight to a background star at distance $`D_{}`$ within the Einstein radius

$$R_\mathrm{E}=\sqrt{\frac{2R_\mathrm{S}}{D_{}}r(D_{}r)},$$
(1)
where $`R_\mathrm{S}`$ denotes the Schwarzschild radius. The lensing rate is proportional to the solid angle that is swept by the lens, which depends on the (projected) velocity $`v`$ and the Einstein radius $`R_\mathrm{E}`$ of the lens. Moreover, the lensing rate is proportional to the distance $`D_{}`$ of the source star and to the number density $`n`$ of lenses (neutron stars), which we take as constant for this simple estimate. Let us for the moment assume that all neutron stars are distributed homogeneously in a cylinder with radius $`10`$ kpc and height $`1`$ kpc and that the number of neutron stars in our Galaxy is $`N_{\mathrm{NS}}=10^9`$. By putting $`D_{}=10`$ kpc, $`r=D_{}/2`$, $`M_{\mathrm{NS}}=1M_{\mathrm{}}`$ and $`v=200\mathrm{km}/\mathrm{s}`$ we find

$$pD_{}R_\mathrm{E}nv10^7\mathrm{yr}^1.$$
(2)
Presently, a few tens of millions of stars have been monitored over some years by various collaborations. This suggests that about a handful of by now observed lenses should be neutron stars. This rough estimate is consistent with the estimate given in Gould (G00 (2000)). The difficulty is to tell which of the observed microlensing events was due to a neutron star. There are only about $`10^5`$ galactic pulsars (Narayan Narayan (1987)), which gives a pulsar lensing rate of the order $`10^{11}\mathrm{yr}^1`$. Typically, a present-day experiment will therefore not find microlensing events that are due to pulsars.

In order to calculate the lensing rate and the optical depth more precisely we have to take into account that the number density of neutron stars is not constant within our galaxy. There are several models for the pulsar number density $`n_{\mathrm{PSR}}`$ in our galaxy and we assume that the neutron star number density $`n`$ is proportional to it

$$n(R,z)=\frac{N_{\mathrm{NS}}}{N_{\mathrm{PSR}}}n_{\mathrm{PSR}}(R,z),$$
(3)
where $`R`$ is the radial distance of the pulsar/neutron star to the galactic centre in the galactic plane and $`z`$ is the height.

Following Hartman et al. (Hartman (1997)) we use two different models for the radial dependence of the number density and a model for the $`z`$-dependence, which was given by Lyne et al. (Lyne (1998)). The two models from Hartman resemble the models from Narayan (Narayan (1987)) and Johnston (Johnston (1994)). Figure 1 shows the number densities of the four models in the galactic plane plotted against the distance from the galactic centre.

We use the assumption (Narayan Narayan (1987)), that the number density can be factorised

$$n(R,z)2\pi R\mathrm{d}R\mathrm{d}z=N_{\mathrm{NS}}\left[n_R(R)\mathrm{\hspace{0.17em}2}\pi R\mathrm{d}R\right]\left[n_z(z)\mathrm{d}z\right].$$
(4)
The first model from Hartman is exponential in $`R`$,

$$n_{\mathrm{H1}}(R)=\frac{1}{2\pi R_\mathrm{W}^2}\mathrm{exp}\left(\frac{R}{R_\mathrm{W}}\right),$$
(5)
where $`R_\mathrm{W}=5\mathrm{kpc}`$. There is some evidence for a deficit of pulsars in the galactic centre which was shown by Johnston (Johnston (1994)). The second model from Hartman takes this into account

$$n_{\mathrm{H2}}(R)=\frac{c_{\mathrm{H2}}}{2\pi R_\mathrm{w}^2}\mathrm{exp}\left(\frac{\left(RR_{\mathrm{max}}\right)^2}{2R_\mathrm{w}^2}\right),$$
(6)
where $`R_\mathrm{w}=1.8`$ kpc and $`R_{\mathrm{max}}=3.5`$ kpc. The normalisation constant is $`c_{\mathrm{H2}}0.204`$ for the given choice of $`R_{\mathrm{max}}`$. For Gaussian distribution around the centre $`c_{\mathrm{H2}}=1`$. For the $`z`$-dependence we use the Gaussian

$$n_z(z):=\frac{1}{\sqrt{2\pi }\sigma }\mathrm{exp}\left(\frac{1}{2}\frac{z^2}{\sigma ^2}\right)$$
(7)
with $`\sigma =0.45\mathrm{kpc}`$ (Lyne et al. Lyne (1998)).

We transform the equations (5), (6) and (7) into spherical coordinates $`(r,\theta ,\varphi )`$ with the sun at the origin via

$`z`$ $`=`$ $`r\mathrm{sin}\theta ,`$ (8)
$`R^2`$ $`=`$ $`r^2\mathrm{cos}^2\theta +R_{\mathrm{SC}}^22rR_{\mathrm{SC}}\mathrm{cos}\theta \mathrm{cos}\varphi .`$ (9)
The distance from the sun to the centre of the galaxy is given by $`R_{\mathrm{SC}}=8.5\mathrm{kpc}`$.

Hansen and Phinney (Hansen (1997)) suggested a Maxwell distribution for the kick velocities of neutron stars:

$$f(v)=\sqrt{\frac{2}{\pi }}\frac{v^2}{\sigma _v^3}\mathrm{exp}\left(\frac{v^2}{2\sigma _v^2}\right)$$
(10)
where $`\sigma _v=190\mathrm{kms}^1`$, which corresponds to a mean velocity at birth of $`300\mathrm{kms}^1`$. For pulsars that are older than $`10^7`$ years they find a mean velocity of $`198\pm 53\mathrm{kms}^1`$. We can now calculate the lensing rate, which is given by

$$p(D_{},\theta ,\varphi )=N_{\mathrm{NS}}\pi v_{}_0^D_{}R_\mathrm{E}(r,D_{})n(r,\theta ,\varphi )dr.$$
(11)
It only depends on the part of the velocity vector that lies in the lens plane. Therefore we have to calculate the mean velocity projection on the lens plane $`v_{}`$ which is given by

$$v_{}=\frac{1}{\pi }_{\frac{\pi }{2}}^{\frac{\pi }{2}}v\mathrm{cos}\alpha \mathrm{d}\alpha =\frac{2v}{\pi }.$$
(12)
Figures 2 and 3 show the lensing rate as a function of the galactic longitude \[$`\varphi (l/180\mathrm{°})\pi `$\] and as a function of the distance to the light source, respectively. It is most promising to search for neutron stars that act as microlense in the direction of the galactic bulge, but it is also reasonable to look into the disk at galactic longitudes below $`50\mathrm{°}`$.

### 2.2 Optical depth

Analogous to (2) we estimate the order of the optical depth,

$$\tau nD_{}R_\mathrm{E}^2\frac{pR_\mathrm{E}}{v}p\mathrm{\hspace{0.17em}0.1}\mathrm{yr}10^8.$$
(13)
A more thorough estimate for the optical depth is obtained from

$$\tau (D_{},\theta ,\varphi )=_0^D_{}\pi R_\mathrm{E}^2(r,D_{})n(r,\theta ,\varphi )dr.$$
(14)
Figures 4 and 5 show the optical depth as a function of the galactic longitude and as a function of the distance to the light source, respectively.

### 2.3 Typical timescale

Consider a lens of mass $`M`$, passing the line of sight of a distant source at a perpendicular velocity $`v_{}`$. The typical time scale $`t_{\mathrm{var}}`$ of the microlensing event is given by

$`t_{\mathrm{var}}`$ $``$ $`{\displaystyle \frac{R_\mathrm{E}}{v_{}}}`$ (15)
$``$ $`25\mathrm{d}\left({\displaystyle \frac{M}{M_{\mathrm{}}}}\right)^{1/2}\left({\displaystyle \frac{D}{1\mathrm{kpc}}}\right)^{1/2}\left({\displaystyle \frac{200\mathrm{km}/\mathrm{s}}{v_{}}}\right),`$
where $`D=r(D_{}r)/D_{}`$. For $`M=1.4M_{\mathrm{}}`$ we find that the duration of a neutron star microlensing event might range from a few days for close and fast neutron stars up to a few months for neutron stars close to the galactic bulge.

The distribution of timescales is estimated in figure 6. We plot

$`{\displaystyle \frac{\mathrm{d}p(t_{\mathrm{var}})}{\mathrm{d}\mathrm{ln}t_{\mathrm{var}}}}`$ $`=`$ $`{\displaystyle \frac{2N_{\mathrm{NS}}}{t_{\mathrm{var}}^2}}`$ (16)
$`\times `$ $`{\displaystyle _0^D_{}}R_\mathrm{E}^3(r,D_{})f\left[{\displaystyle \frac{R_\mathrm{E}(r,D_{})}{t_{\mathrm{var}}}}\right]n(r,\theta ,\varphi )dr`$
at $`\theta =\varphi =0`$ for $`D_{}=3.25`$ kpc and $`D_{}=7.5`$ kpc, respectively. We use the velocity distribution (10), since we expect young neutron stars at small scale height.

As can be seen from the comparison of the timescale distributions for sources close to the galactic centre and sources at moderate distances, the width of the timescale distribution is dominated by the sources close to the galactic centre. In both models for the radial distribution the timescale distribution peaks at durations of $`20`$ days.

### 2.4 Astrometric microlensing

Astrometric microlensing makes use of the shift of the centroid of the combined (unresolved) images of the light source. While the light amplification used in the photometric effect drops outside the Einstein ring as

$$A1+2\left(\frac{\theta _\mathrm{E}}{\theta _{\mathrm{sl}}}\right)^2,$$

the centroid shift is given by

$$\mathrm{\Delta }\theta \frac{\theta _\mathrm{E}^2}{\theta _{\mathrm{sl}}}.$$

Here $`\theta _\mathrm{E}R_\mathrm{E}/r`$ and $`\theta _{\mathrm{sl}}`$ is the angle between the source and the lense. For a nearby neutron star we may assume $`D_{}r`$, thus $`\theta _\mathrm{E}9(M/M_{\mathrm{}})^{1/2}(100\text{ pc}/r)^{1/2}`$ mas, giving rise to a typical centroid shift of

$$\mathrm{\Delta }\theta 81\left(\frac{M}{M_{\mathrm{}}}\right)\left(\frac{100\text{ pc}}{r}\right)\left(\frac{1\mathrm{}}{\theta _{\mathrm{sl}}}\right)\mu \text{as}.$$
(17)
GAIA will achieve an astrometric accuracy of $`11\mu `$as at $`V=15`$ mag, which drops to $`160\mu `$as at $`V=20`$ mag (Perryman et al. Pea (2001)), while SIM will resolve $`4\mu `$as for sources as faint as that magnitude (Boden, Unwin & Shao BUS (1997)). With the help of VLTI an astrometric accuracy of the order of $`10\mu `$as can be expected (Delplancke, Górski & Richichi DGR (2001)) with a limiting magnitude $`K=1920`$ (available from the middle of 2004). Thus it will be possible, at least in principle, to measure the centroid shift of background stars separated less than a few $`\mathrm{}`$ from a nearby neutron star.

For a survey mission like GAIA, the optical depth and timescale of lensing may be estimated as above, simply correcting for the larger cross section due to the astrometric effect. If one adopts an accuracy of $`100\mu `$as, and aims at measuring the centroid shift with an error of $`10\%`$, $`\theta _{\mathrm{sl}}/\theta _\mathrm{E}`$ might be as big as $`10`$, thus increasing the cross section from photometric lensing by a factor of $`100`$ and the time scale by a factor of $`10`$. However, the sampling of GAIA will be sparse compared with ground-based ongoing microlensing experiments and only a small fraction of all astrometric microlensing events will be detected (Belokurov & Evans BE (2001)).

On the other hand SIM will focus on astrometry of preselected targets. For photometrically-detected microlensing events it will be possible to detect the centroid shift and to obtain high accuracy mass determinations.

## 3 Catalogue comparison

As a consistency check of the above estimates of the photometric lensing rates we investigate whether any pulsar has acted as a microlens thus far.

The positions of 690 candidates for gravitational microlensing events towards the galactic bulge and the spiral arms (OGLE101010 Catalogue of 520 Microlensing Events in the Galactic Bulge:
http://astro.princeton.edu/~wozniak/dia/lens/, Wozniak et al. Wozniak (2001); MACHO111111http://wwwmacho.mcmaster.ca/Data/bulgefts.html, Alcock et al. Alcock1 (1997), Alcock2a (1999), Alcock2 (2000); EROS121212alerts 1998 - 2000, http://www-dapnia.cea.fr/Phys/Spp/ Experiences/EROS/alertes.html, Derue et al. Derue1 (1999), Derue2 (2001)) have been compared with the positions of 1308 radio pulsars from the ATNF Pulsar Catalogue131313http://wwwatnf.atnf.csiro.au/research/pulsar/ catalogue, which is based on the Princeton Pulsar Catalogue141414 http://pulsar.princeton.edu/pulsar/catalog.shtml (Taylor, Manchester & Lyne Taylor (1993)), the Parkes Multibeam Pulsar Survey151515 http://wwwatnf.atnf.csiro.au/research/pulsar/ pmsurv/, (Manchester et al. Manchester (2001)) and the Swinburne Intermediate Latitude Pulsar Survey161616http://astronomy.swin.edu.au/pulsar/ (Edwards et al. Edwards (2001)). Table 1 lists the number of records of the individual groups. Some of the microlensing events have been detected by more than one group. We did not account for this in the table.

All positions are given in the J2000 coordinate system. However, not all of the measurements were made at the same time. Therefore we have to consider a change of the pulsar position that will lead to a deviation angle $`\alpha vt/r`$, where $`t`$ is the time between the pulsar detection and the microlensing event, $`v_{}`$ is the (projected) pulsar velocity and $`r`$ the pulsar distance. For a conservative estimate that includes pulsars that are close, fast, and have been detected long ago, we take $`r=100`$ pc, $`v=500`$ km/s, $`t=10`$ yrs, which gives $`\alpha 10\mathrm{}`$. We found no candidate pair with an angular separation below this value. In a second step we looked for the closest pair, which is PSR J1807-2715 and the OGLE microlensing event sc18 2757. The angular separation is $`0.08\mathrm{deg}`$. According to the Priceton Pulsar Catalogue PSR J1807-2715 is at an estimated distance of $`9.67`$ kpc and no proper motion is reported. Thus it is highly unlikely that it can be responsible for the OGLE event sc18 2757. These findings are in agreement with our estimate, which predicts an extremely small lensing rate for pulsars, i.e. $`10^4`$/yr for $`10^7`$ monitored stars.

## 4 A new search strategy

We have shown that probably some of the gravitational microlensing events towards the galactic bulge are due to neutron stars. For a specific microlensing event with a duration of a few days to a few months, one is interested in finding out whether the lens is a neutron star or a different object. There is not very much hope that with present day experiments a microlensing event will be detected whose lens is a known pulsar. Our catalogue comparison has confirmed this result. To raise the low probability it is clear that a better search strategy is needed. With the prospects of astrometric lensing, we can expect that high precision measurements of masses will become possible, but still the problem of identifying the nature of the lense will remain, especially for masses of $`13M_{\mathrm{}}`$, since many possible objects are stable in that mass range.

A unique identifier of a neutron star is the possibility that it is observed as a radio pulsar. Instead of choosing a field of background stars we suggest to choose a pulsar and to monitor the background stars around it. The solid angle $`\mathrm{\Omega }`$ that is swept during $`\mathrm{\Delta }t`$ by the Einstein-ring of a moving Schwarzschild-lens is given by

$$\mathrm{\Omega }=\frac{1}{4\pi r^2}\left(\pi R_\mathrm{E}^2+2R_\mathrm{E}v_{}\mathrm{\Delta }t\right)\frac{R_\mathrm{S}}{2r}\left(1+\frac{2\mathrm{\Delta }t}{\pi t_{\mathrm{var}}}\right),$$
(18)
for $`rD_{}`$. Our strategy therefore must be to observe distant background stars in the direction of a close, fast-moving neutron star. Table 2 shows some high proper motion pulsars from the Princeton Pulsar Catalogue. We include all pulsars whose proper motion have been measured to be $`|\mu |>100`$ mas/yr and keep only those with proper motion estimates better than $`10\%`$. We expect that pulsars at low galactic latitude have the highest probablility of acting as a gravitational lens. This expectation is confirmed by the result of a search for appropriate background stars, based on the USNO-A2.0 (Monet et al. Monet (1998)) and GSC2.2.1171717http://www-gsss.stsci.edu/gsc/gsc2/GSC2home.htm catalogues. These catalogues contain $`526,280,881`$ resp. $`455,851,237`$ sources and both cover the entire sky. The USNO-A2.0 catalogue above declination $`30\mathrm{°}`$ (only those plates are relevant in our context, see below) goes as deep as $`21`$ in $`O`$ (blue plates) and $`20`$ in $`E`$ (red plates), whereas the GSC2.2.1 catalogue is limited by $`F=18.5`$ in the red and $`J=19.5`$ in the blue. We search in a radius of $`15\mathrm{}`$. The number of hits are shown in the last column of table 2 for USNO-A2.0/GSC2.2.1. This procedure selects three pulsars as interesting candidates for lensing.

The only candidate at high galactic latitude, PSR J0437-4715, is a binary pulsar and the USNO-A2.0 object 0375-01566415 (epoch 1981) is separated by only $`0.5\mathrm{}`$ after the measured pulsar proper motion is used to extrapolate to epoch 1981. Thus the USNO object seems to be identical with the binary companion identified by Bell, Bailes & Bessell BBB (1993), probably a white dwarf with $`R=20.1`$, since they find no other star brighter than $`R=23.8`$ within $`6\mathrm{}`$ of the radio position.

The closest object (USNO-A2.0 object 1500-09063417, GSC2.2.1 object N013123248640) to PSR J2225+6535 is separated by $`14.7\mathrm{}`$ and thus is too far away from the pulsar for significant lensing.

The most promising candidate is the pulsar PSR J1932+1059 (B1929+10). This pulsar is at a (dispersion measure) distance of about $`0.17\mathrm{kpc}`$ (Princeton Pulsar Catalogue) and moves $`105`$ mas/yr. It lies towards the galactic disk ($`l=47.4\mathrm{°},b=3.9\mathrm{°}`$) and is found to have $`7`$ objects brighter than red magnitude $`20`$ within $`15\mathrm{}`$. From the $`9`$ hits shown in table 3, $`2`$ objects are counted twice (see caption of table). In figure 7 (top panel) we map $`9\mathrm{}`$ around the pulsars radio position in 1991. The crosses indicate its estimated position in the years 2001, 2011 and 2021. The pulsar moves away from the objects 1,5,6 and towards 2,7,8 and remains at an almost constant distance to 3. At epoch 2001(2011,2021) objects 5 and 7 are separated from the pulsar by $`4.8\mathrm{}(5.7\mathrm{},6.8\mathrm{})`$ and $`5.9\mathrm{}(4.9\mathrm{},4.0\mathrm{})`$, respectively.

In 1994, Pavlov et al. (Pavlov (1996)) identified the optical counterpart of PSR J1932+1059 using the Hubble Space Telescope. The offset from the estimated radio position at epoch 1994.5 was $`0.39\mathrm{}`$. As a by-product they found 6 additional stars within $`4\mathrm{}`$ of the radio position. We show their positions (from table 5 of Pavlov et al. Pavlov (1996)) in the bottom panel of figure 7, which is a magnification of the top panel by a factor of 2. The pulsar is almost in the middle between HST objects 1 and 2 and approaches HST object 5. At 2010 (2015) the pulsar will be separated from HST objects 1, 2 and 5 by $`2.66\mathrm{}(3.16\mathrm{})`$, $`1.74\mathrm{}(1.95\mathrm{})`$ and $`1.81\mathrm{}(1.35\mathrm{})`$, respectively. According to Pavlov et al. (Pavlov (1996)), objects 1–3, 5 and 6 are perhaps very distant white dwarfs strongly reddened by interstellar extinction.

The angular scale of the Einstein radius of PSR J1932+1059 is $`8.2`$ mas, thus there is no photometric lensing effect. However, there is a chance to detect the astrometric effect. Let us first consider the possibility of detecting lensing with the help of the objects from table 3. These objects are considerably brighter than the six HST objects and therefore astrometric measurements should be simpler. According to equation (17) the images of 5 and 7 move between 2010 and 2020 by $`2\mu `$as and $`3\mu `$as. The pulsar and the objects 5 and 7 are to a good approximation on a straight line; the centroid shift is directed along the same line. The accuracy of SIM will be about $`4\mu `$as, thus the effect seems to be just below what can be detected.

Turning to the HST objects the numbers are much more encouraging. Now the pulsar moves away from 1 towards 5 and stays at an almost constant distance to 2. Between 2010 and 2015 the image of 1 will move away from the pulsar by $`4\mu `$as, the image of 2 will move by $`11\mu `$as almost parallel to the proper motion of the pulsar and the image of 5 will move by $`12\mu `$as towards the pulsar. Thus the angular separation between 1 and 5 should decrease by $`8\mu `$as. If SIM is able to resolve the rather faint stars 1,2 and 5 a measurement of the astrometric lensing of a pulsar should be possible. The best strategy would be to observe the field a few times per year during the lifetime of SIM.

In the above estimates the uncertainties from the distance estimate are most crucial. We used a distance of $`170`$ pc based on the dispersion measure. This is not very reliable due to inhomogeneities within the interstellar medium. Several attempts to measure the parallax of PSR J1932+1059 gave rise to a confusing situation. Salter et al. (Sea (1979)) obtained a parallax distance of $`3080`$ pc, while Backer & Sramek (BS (1982)) found a lower limit of $`250`$pc. Allowing for the range $`30300`$ pc implies an uncertainty of an order of magnitude for the predicted image shifts.

## 5 Conclusions

We have argued that measuring the mass of stand-alone pulsars will be possible in the near future by means of gravitational microlensing.

The same strategy could also be used to measure the mass of isolated neutron stars, which are not seen as radio pulsars. Isolated neutron stars \[for recent reviews see Treves et al. (Treves (2000)) and Popov et al. (Popov (2000))\] may be detected as soft X-ray sources due to accretion in the case of old neutron stars or due to energy released in the cooling of young neutron stars. Several candidates have been found by the ROSAT181818http://wave.xray.mpe.mpg.de/rosat satellite. For two of those candidates, possible optical counterparts have been identified \[RX J185635-3754 by Walter & Mathews (Walter (1997)) and RX J0720-3125 by Kulkarni & van Kerkwijk (Kulkarni (1998))\], which makes them especially interesting. After the submission of the first version of this work, Paczyński (P01 (2001)) pointed out that astrometric lensing of RX J185635-3754 might be possible with the HST in 2003.

In order to increase the list of possible candidates it is important to obtain as many proper motion and parallax measurements of pulsars and isolated neutron stars as possible. Candidates towards the galactic bulge should be especially considered with priority, since the chance of finding a large number of sources for lensing is maximised in that way.

We conclude that measuring neutron star (pulsar) masses by means of gravitational microlensing is possible. We suggest to observe PSR J1932+1059 and three stars around it with SIM over several years to obtain such a mass measurement.

Note added in proof. We learned after acceptance of this work that the possible mass determination of isolated pulsars by gravitational lensing has also been investigated in J. E. Horvath, 1996, MNRAS 278, L46. Recently, a new parallax distance to PSR J1932+1059, $`r=331\pm 10`$ kpc, has been published by W. F. Brisken et al. astro-ph/0204105; thus the astrometric effect is a factor of 2 smaller than estimated in this work. Still, it should be possible to detect the astrometric lensing of PSR J1932+1059 with SIM.

###### Acknowledgements.

We would like to thank W. Kegel, D. Lorimer, B. Paczyński and F. Weber for very helpful conversations and references to the literature. This research has made use of the VizieR catalogue access tool and the Aladin Interactive Sky Atlas, CDS, Strasbourg, France. D. J. S. would like to thank the Austrian Academy of Sciences for financial support.


format: text

text:
# Detection of the 2175Å Extinction Feature at 𝑧=0.83

## 1 INTRODUCTION

Just as the properties of galaxies and stars evolve with cosmic time, so must the properties of the interstellar medium. In particular, with changes in metallicity and ionizing backgrounds we might expect the properties of dust to change with redshift. Locally we can make precise measurements of extinction curves as a probe of the dust grains by comparing the spectral energy distributions (SED) of reddened and unreddened stars of the same spectral type (the pair method, Massa, Savage & Fitzpatrick (1983)), but the need to measure the spectra of individual stars limits this approach to the Galaxy, its satellites and its nearest neighbors.

The average Galactic extinction curve is well modeled over a broad wavelength range as a function of only the ratio of total to selective extinction, $`R_V`$ (e.g. Cardelli, Clayton & Mathis 1989, hereafter CCM; Fitzpatrick 1999). These extinction curves include a spectral feature, the 2175Å bump, and appear to apply to a wide range of interstellar dust environments in the Galaxy. It is not, however, a universally valid law. While some lines of sight in the Magellanic clouds show similar extinction laws, others, like the 30 Doradus region of the Large Magellanic Cloud (LMC) show weaker 2175Å features and steeper far UV extinction curves (e.g., Fitzpatrick, 1986). The most extreme deviations from the Galactic extinction curves are found in the bar of the Small Magellanic Cloud (SMC, Gordon & Clayton 1998). Some newer studies (e.g., Clayton, Gordon & Wolff, 2000) have found lines of sight with weak 2175Å features and steep far UV extinctions in the Galaxy as well. M31 is the most distant galaxy to which the pair method can be applied; the M31 extinction curve seems similar to that of the Galaxy but may have a weaker 2175Å bump (Bianchi et al., 1996).

The studies of extinction in other galaxies are in many cases based on the examination of the attenuation induced by dust in the SED of the galaxies. But the physical properties of dust and its geometrical distribution affect in a complex way the intrinsic SED of the galaxy stellar population, which must be also modeled. The extinction curves inferred from these indirect methods are strongly model-dependent and, in many cases, not unique due to model degeneracies (Bianchi et al., 1996; Gordon, Calzetti & Witt, 1997). In M33, Gordon et al. (1999) noticed that a depression found in the SED should be related to the 2175Å feature, according to a radiative transfer model including attenuation by dust. A very weak 2175Å feature is also found by Rosa & Benvenuti (1994) by comparing spectra of ionizing star clusters in HII regions of M101 with population synthesis models. However, the modeling of the SED of 30 starburst galaxies (Gordon, Calzetti & Witt, 1997) implies that the dust has extinction properties similar to that of the SMC, lacking a 2175Å feature.

Some efforts have also been made in AGNs to study the effects of extinction by comparing the observed SEDs with conveniently selected templates. Crenshaw et al. (2001) have compared the spectra of the Seyfert 1 galaxies NGC 3227 and NGC 4151 (assumed unreddened), and found that the extinction in NGC 3227 increases very steeply toward the far UV (steeper than in the SMC) and that there is no evidence of the 2175Å bump. In QSOs, the comparison of the observed spectra with composite spectra has not revealed traces of the 2175Å feature (Pitman & Clayton, 2000; Vanden Berk, 2001). All the above-mentioned studies suggest that the presence of a significant 2175Å bump in the average extinction curve of the Milky Way is rare (Pitman & Clayton, 2000). In any case, the validity of this or any other conclusion about the behavior of the extinction beyond the Local Group is limited by the lack of a direct method to obtain the extinction curve. This limitation has serious consequences for the study of galaxies (see Gordon et al., 2000, and references therein) and in cosmology (see Falco, Kochanek & Muñoz, 1998, and references therein).

There are now two techniques which can extend the measurement of precise extinction curves into the realm of cosmology. One is based on the use of the relationship between the variations in color and absolute magnitude in SNIa light curves (Riess, Press & Kirshner, 1996). This technique has been used to estimate several reddening ratios from nearby supernovae which are in agreement with the MW average extinction curve. The other is an extension of the standard pair method to obtain extinction curves in gravitational lens galaxies. The method compares two different images of the same lensed QSO which follow different paths through the lens galaxy (Nadeau et al., 1991). As in the pair method each image has the same intrinsic spectrum but suffers from a different amount of extinction along each path through the lens galaxy. Using broad-band optical and NIR photometry, Falco et al. (1999) measured differential extinctions in 23 gravitational lens galaxies. In spite of the coarse spectral sampling, Falco et al. (1999) could fit the average Galactic extinction law to several lens galaxies, finding values for $`R_V`$ in the range 0.64 to 7.2. This range includes some extreme values for $`R_V`$ whose physical significance is not obvious as the CCM law was inferred from values of $`R_V`$ in the 2.5 to 5.5 range. Falco et al. (1999) also derived relatively accurate dust redshifts following the technique proposed by Jean & Surdej (1998). These studies relied on broad band photometry, which can lead to problems distinguishing between extinction and gravitational microlensing and cannot provide complete extinction curves over broad wavelength ranges. These problems can be overcome using spectra of the individual images, particularly using integral field spectroscopy (IFS) to obtain a 2D array of spectra of the lens and its images (Mediavilla et al., 1998).

In this article, we use IFS to search for the 2175Å feature in a cosmologically distant galaxy by selecting a gravitational lens where the lens galaxy shows significant extinction and the redshifted wavelength of the 2175Å feature, $`2175(1+z_l)`$Å is below the Lyman break in the spectrum of the source quasar at $`912(1+z_s)`$Å. We performed the experiment on the two-image, $`1\stackrel{}{\mathrm{.}}11`$ separation quasar lens SBS 0909+532 (Kochanek et al., 1997), where the $`z_l=0.83`$ lens and $`z_s=1.38`$ source redshifts (Oscoz et al., 1997; Lubin et al., 2000) satisfy the wavelength criterion, and Falco et al. (1999) found there was significant differential extinction $`\mathrm{\Delta }E(BV)=0.20\pm 0.03`$ based on the image flux ratios measured in broad band imaging.

## 2 DATA ANALYSIS

### 2.1 Observations and Data Reduction

The spectra of SBS 0909+532 were taken on January 18, 2001 with the 2D-fiber system INTEGRAL (Arribas et al., 1998), at the 4.2-m William Herschel Telescope (WHT). INTEGRAL links the Nasmyth focus of the WHT with the WYFFOS spectrograph. We used the SB1 fiber bundle that consists of a rectangular array of 175 fibers, each $`0\stackrel{}{\mathrm{.}}45`$ in diameter, covering an area of $`7\stackrel{}{\mathrm{.}}9\times 6\stackrel{}{\mathrm{.}}7`$ on the sky, and a ring with diameter $`90\mathrm{}`$ formed by 30 fibers of the same size as those in the inner array. The distance between adjacent fibers is $`0\stackrel{}{\mathrm{.}}5`$ in the central rectangle. At the entrance of the spectrograph, the fibers are aligned to form a pseudo-slit.

We took three exposures of 1800s each. We used the 300lines/mm grating, R300B, that yields a linear dispersion of 6.0Å pixel-1 covering the wavelength range from $`3400`$ to $`9200`$Å. The seeing, measured from a guide star was $`1\mathrm{}`$.

For our data reduction, we used the IRAF111IRAF is distributed by the National Optical Astronomy Observatories, which is operated by the Association of Universities for Research in Astronomy, Inc. (AURA) under cooperative agreement with the National Science Foundation software package. The procedure included bias subtraction, cosmic-ray rejection, spectra extraction, throughput correction, wavelength calibration, and sky subtraction. All these steps were applied to each one of the three individual frames and to an average frame obtained by combining them. For additional details on reduction and analysis of data taken with optical fibers, see Arribas, Mediavilla & Rasilla (1991).

In Figure 1 we show the 2D distribution of spectra from the average frame over a spectral range from 6000 to 7200 Å including the MgII$`\lambda `$2798 emission line. The positions of the two compact components A and B are also shown.

### 2.2 Continuum Maps

From the 2D collection of spectra corresponding to the averaged frame (a subset of them around components A and B is shown in Figure 1) we interpolated 20 continuum maps covering all the observed spectral range. The central wavelengths and widths of the continua are included in Table 1. We fitted and removed the strong emission lines present in the QSO spectra.

In Figure 2 we present six of these continuum maps corresponding to the averaged frame centered at 8727Å, 7233Å, 5739Å, 4447Å, 4006Å, and 3487Å. The two compact components of SBS 0909+532 appear well separated in the reddest map (8727Å, seeing$`1\stackrel{}{\mathrm{.}}05`$). Notice the continuous fading of image B from the red to the blue maps. In fact, in the bluest map (3487Å) the B image is no longer visible. This decrease in the B/A flux ratio towards the blue is explained by dust absorption in the lensing galaxy (Falco et al., 1999).

### 2.3 A and B Spectra

Because of differential atmospheric refraction (DAR), the location of any of the components in the telescope focal plane changes with wavelength. In fact, we measured a displacement of about $`0\stackrel{}{\mathrm{.}}6`$ (greater than the diameter of a fiber face) between the bluest and reddest continuum maps. Although the displacement is modest, it implies that none of the spectra of Figure 1 can be exactly identified with one of the components throughout the full spectral range. This is a general problem of spectroscopy in a large spectral range that admits an a posteriori correction when 2D spectroscopy is available (Arribas et al., 1999). Following the procedure described in this reference, we calculated a DAR model from the location of component A in the 20 continuum maps distributed along the full spectral range. Next, from the original spectra we interpolated a new 2D set of spectra corrected for DAR (see Arribas et al. 1999 for the details). From these spectra we obtained a map in the 3500-9000Å spectral range. We located the centroid of the A component in this map and used the coordinates offset of B from A given in Lehár et al. (2000) to infer the location of the B component. Finally, we interpolated at these locations the A and B spectra (see Figure 3) from the ensemble of spectra corrected for DAR. The changes from the uncorrected spectra are modest, as expected.

## 3 RESULTS

### 3.1 Differential Extinction Curve from Continuum Images

To derive the B/A flux ratios from the ensemble of continuum maps covering the full observed spectral range, we fit the maps using PSFs derived from the continuum images of a star observed on the same night. We used as fitting parameters the separation between the two QSO images, their relative intensity, the sky background, and the blurring of the PSF due to the seeing. We followed this procedure with each one the three independent exposures. After subtracting the quasar model, the pattern of residuals show no emission from the galaxy at any wavelength, so we conclude that our photometric model with only two point-like sources reproduces successfully the observed continua. Acording to the HST data (Lehár et al., 2000) the lens galaxy is not detected in the $`V`$ band and its integrated magnitude is only 5% of component A flux in the $`I`$ band. The location of the two quasar components determined by our photometric model is in excellent agreement with the results obtained with the HST by Lehár et al. (2000) which strengthens the validity of our photometric decomposition.

In Table 1 and Figure 4 we present the mean magnitude differences, $`m_Bm_A`$, derived from the PSF fitting. To estimate the errors we have used the dispersion between the results obtained from each independent exposure. We include in this Figure the $`m_Bm_A`$ differences from Lehár et al. (2000) and Kochanek et al. (1997) corresponding to broad band data; the agreement is excellent.

### 3.2 Differential Extinction Curve from Emission Line Images

Proceeding in the same way as for the continua, we have also computed the $`m_Bm_A`$ differences for the emission in the MgII$`\lambda `$2798 and CIII\]$`\lambda `$1909 lines. To estimate and subtract the continuum below each line, we fitted a straight line to adjacent points. The magnitude differences in the emission lines are approximately 0.3 mag above that of the continuum ($`m_Bm_A`$, see Figure 4). These differences are probably created by microlensing from the stars in the lens galaxy which can differentially magnify the continuum and emission line regions of the quasar because of their differing physical sizes (Schneider, Ehlers & Falco, 1988). While there are too few emission lines to estimate an extinction curve, the emission lines are consistent with the extinction properties estimated from the continuum. We have also tried to compute the $`m_Bm_A`$ difference from the CIV\]$`\lambda `$1548, but the maps were too noisy and the results of the PSF fitting were not consistent.

### 3.3 Differential Extinction Curve from A and B Uncontaminated Spectra

The method followed in Section 3.1 could be repeated with narrower continua to obtain the differential extinction curve with improved spectral sampling that could even match the spectral resolution of the data. However, the spectral sampling that can be achieved with this method is, in practice, limited by the noise in the data and the PSF.

To improve the spectral resolution we can try, alternatively, and only because the lens galaxy is undetected, to obtain the extinction curve directly from the $`A`$ and $`B`$ spectra after correcting them for the cross-contamination induced by the seeing. We can estimate the contribution of $`A`$ to the spectrum measured at the position of $`B`$ by measuring the spectrum at position $`C`$ located at the same distance from $`A`$ as in the image $`B`$ but in the opposite direction (see Figure 1). Similarly we obtain spectrum $`D`$ as an estimate of the contribution to $`A`$ from $`B`$. Our estimate of the seeing-corrected spectra are the differences $`A_u=AD`$ and $`B_u=BC`$.

We can check the results by comparing the spectral ratio to the estimates made by directly fitting the continuum maps. They agree reasonably well, with a $`0.15`$ mag drift between the blue and the red wavelengths. The drift appears to be due to small asymmetries in the PSF. We can improve the agreement by defining $`A_u=A(1ϵ)D`$ and $`B_u=B(1+ϵ)C`$ with $`ϵ=0.1`$ to correct for the asymmetries. Figure 5 shows our final estimate of the magnitude differences $`[m_{B_u}m_{A_u}](1/\lambda _0)`$ after smoothing by a $`\sigma =30`$Å Gaussian. This curve is in good agreement with the difference magnitudes obtained from the continuum images at all wavelengths except in the vicinity of the emission lines. As stated in the previous section, this could arise from microlensing that is not affecting the lines to the same degree as the continuum. This implies that, in some spectral regions (like the one including the CIV$`\lambda `$1548 line in SBS 0909+532) observations made with narrow-band filters could be problematic for obtaining the extinction curves since the contribution from the continuum and the strong emission lines must be separated. (Notice also the presence of the Mg II$`\lambda \lambda `$2796,2803 absorption doublet arising from the lens galaxy).

## 4 DISCUSSION

Assuming that any microlensing effects are achromatic, the magnitude difference (Figure 4) is a simple function of the magnification ratio $`\mathrm{\Delta }M`$ in magnitudes, the differential extinction $`\mathrm{\Delta }E(BV)`$ and the extinction law $`R(\lambda )`$ (Falco et al., 1999),

$$m_B(\lambda )m_A(\lambda )=\mathrm{\Delta }M+\mathrm{\Delta }E(BV)R\left(\frac{\lambda }{1+z_l}\right).$$
(1)
We fitted the $`m_Am_B`$ data from the continua using equation (1) and the CCM parameterization of the extinction curve $`R(\lambda )`$. In a first step, the redshift is restricted to the value determined directly from the absorption lines of the lens galaxy ($`z=0.83`$). The best fit is obtained for $`R_V=2.1\pm 0.9`$, $`\mathrm{\Delta }M=0.2\pm 0.2`$ and differential extinction $`\mathrm{\Delta }E(BV)=0.21\pm 0.02`$. A visual inspection of the fit (see Figure 4) shows reasonable agreement with the CCM extinction curve for the MW ($`\chi ^2/dof2`$). In spite of the uncertainties in $`R_V`$, the fit of the 2175Å bump is good in both shape and central wavelength. When, in a second step, changes not only in $`R_V`$ but also in $`z`$ were allowed, we found that $`R_V`$ was very poorly determined, but we obtained a useful estimate for the redshift, $`z=0.88\pm 0.02`$. The close agreement with the spectroscopic redshift of the lens galaxy ($`z=0.83`$, Oscoz et al., 1997; Lubin et al., 2000) demonstrates the possibilities of the dust-redshift technique (Jean & Surdej, 1998; Falco et al., 1999) and reinforces the reliability of the extinction law.

The detection of a significant 2175Å feature in the extinction curve of a galaxy at $`z=0.83`$ is very noteworthy. Many extragalactic extinction studies have suggested that the 2175Å feature is rare (see Pitman & Clayton, 2000, and references therein). However, these studies were based mainly on starburst galaxies and AGNs, where the dust is subject to processing by radiation and shocks, similar to the environments of star formation regions of the LMC or the SMC, where the 2175Å bump is weak. Thus, the detection of this feature in the extinction curve of the SBS 0909+532 lens (a normal early-type galaxy) and its absence in active galaxies follows naturally under the hypothesis of a connection between activity and suppression of the 2175Å feature.

It has been proposed that the presence of the 2175Å feature is anti-correlated with the FUV extinction, with the feature weakening as the FUV extinction steepens (Clayton, Gordon & Wolff, 2000). Accordingly, the rather low value $`R_V=2.1`$ (although not rare in the extragalactic domain, see Falco et al. 1999 and references therein) that implies a steep FUV rise would be inconsistent with the strong bump (e.g. HD 210121, Clayton, Gordon & Wolff, 2000). Within the uncertainty of the determination of $`R_V`$ the average value for the MW ($`R_V=3.1`$) would also be acceptable and in better accordance with the presence of a significant bump. Additional UV and IR data would be needed to better answer this question.

## 5 CONCLUSIONS

The application of the standard pair method to spectroscopic observations of the gravitationally lensed QSO SBS 0909+532 has allowed us to obtain the first determination of an UV extinction law beyond the Local Group. Our conclusions are:

1. We have detected a significant 2175Å bump in a $`z=0.83`$ galaxy. So far, most studies of extinction in the extragalactic domain have been based on starburst and AGN galaxies which seem to exhibit SMC-like dust properties lacking a 2175Å feature. However, the lens galaxy of SBS 0909+532 seems to contain dust like that of the Galaxy.
2. The Cardelli, Clayton & Mathis (1989) law provides a reasonable fit to the SBS 0909+532 extinction curve with $`R_V=2.1\pm 0.9`$. This value includes, at $`1\sigma `$, the average MW law ($`R_V=3.1`$). However, values of $`R_V2.1`$ will be difficult to reconcile with the proposed experimental correlation between weak bumps and steep FUV extinction.
3. The lens redshift estimated from the dust, $`z_{dust}=0.88\pm 0.02`$, agrees well with the redshift inferred from the absorption lines of the lens galaxy.

We know of about 70 gravitational lenses, many of which show signs of differential extinction (Falco et al., 1999). By applying our techniques to a wider sample of these lenses we should be able to provide a broad outline of the properties and potentially the evolution of extinction curves with cosmic time.

The WHT is operated on La Palma by the Isaac Newton Group (ING) of Telescopes in the Spanish Observatory Roque de los Muchachos of the Instituto de Astrofísica de Canarias (IAC). We would like to thank the support astronomer A. García Pérez and the night assistant N. Mahoney. This work was supported by the P6/88 project of the IAC. V. Motta acknowledges the support of a MUTIS Fellowship from the Agencia Española de Cooperación Internacional (AECI). E.E. Falco & C.S. Kochanek were partially supported by NASA through the grant GO-8252 from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc.


text:
# Journal reference: J. Phys. A: Math. Gen. 35 (2002) 6995-7002 Addenda after Eqs. (34) and(44) Percolation transition in the Bose gas II.

## 1 Introduction

It was perhaps Feynman who first emphasised the importance of long permutation cycles in the description of the $`\lambda `$-transition in liquid helium and, more generally, of Bose-Einstein condensation (Feynman 1953). Some years ago the present author gave a mathematically precise formulation of this idea by introducing the notion of cycle percolation in Bose systems (Sütő 1993, hereafter referred to as I). In I it was suggested that this concept could serve as an alternative characterization of Bose-Einstein condensation (BEC). The examples of the ground state of a system of bosons and the three dimensional ideal Bose gas were treated. For the latter it was shown that BEC indeed implies cycle percolation, but the proof of the absence of percolation in the absence of BEC was missing. Since then, this kind of description has gained some new interest (e.g. Bund and Schakel 1999, Schakel 2000, Martin 2001, Ueltschi 2002). Therefore, it seems to be worthwhile to complete the earlier proof and to make the claim stronger by showing the occurrence of an infinity of macroscopic cycles in the condensate. As it will become clear, our conclusions also apply to the imperfect Bose gas.

Let us recall the principal definitions and relevant results of I. Let $`\mathrm{\Lambda }`$ be a cube of side $`L`$ of the $`d`$ dimensional Euclidean space, and denote $`H_{\mathrm{\Lambda },N}`$ the Hamiltonian of $`N`$ (interacting) bosons in $`\mathrm{\Lambda }`$, taken, e.g., with periodic boundary conditions. By making explicit the summation over the permutations of $`N`$ particles, necessary for symmetrization, the corresponding canonical partition function reads

$`Q_{\mathrm{\Lambda },N}=\mathrm{Tr}_{_{\mathrm{sym}}^N}e^{\beta H_{\mathrm{\Lambda },N}}={\displaystyle \frac{1}{N!}}{\displaystyle \underset{gS_N}{}}\mathrm{Tr}_^NU(g)e^{\beta H_{\mathrm{\Lambda },N}}`$
$`={\displaystyle \underset{\{n_j\}:{\scriptscriptstyle jn_j}=N}{}}\left({\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \frac{1}{j^{n_j}n_j!}}\right)\mathrm{Tr}_^NU(g_{\{n_j\}})e^{\beta H_{\mathrm{\Lambda },N}}`$ (1)
Here $`S_N`$ is the group of permutations of $`\{1,2,\mathrm{},N\}`$ and $`U(g)`$ is the unitary representation of $`g`$ in the $`N`$-particle Hilbert space $`^N`$ where $`=L_2(\mathrm{\Lambda })`$. In the second line of (1) the summation goes over the conjugation classes of $`S_N`$. Each $`g`$ in the same class giving the same contribution, see (I.2.5), $`g_{\{n_j\}}`$ can be any permutation with $`n_j`$ cycles of length $`j`$. Because each term contributing to the sum is positive, cf equation (I.2.18),

$$P_{\mathrm{\Lambda },N}(g)=\frac{1}{N!Q_{\mathrm{\Lambda },N}}\mathrm{Tr}_^NU(g)e^{\beta H_{\mathrm{\Lambda },N}}$$
(2)
can be interpreted as the probability of the occurrence of the permutation $`g`$ in the canonical ensemble. Then, we can introduce random variables, for example $`\xi _1`$, the length of the cycle containing 1 (i.e. particle no 1) and $`n_j`$, the number of cycles of length $`j`$. Constancy of the probability (2) within a conjugation class means that the probability is sensitive only to the number and length of the cycles and not to their contents. Thus, 1 is assigned to any cycle with a probability proportional to the length of the cycle, which yields the conditional probability

$$P_{\mathrm{\Lambda },N}(\xi _1=j|n_j=m)=mj/N$$
(3)
and, as a consequence,

$$P_{\mathrm{\Lambda },N}(\xi _1=j)=\frac{j}{N}\underset{m}{}mP_{\mathrm{\Lambda },N}(n_j=m)\frac{jn_j_{\mathrm{\Lambda },N}}{N}.$$
(4)
Clearly, $`_{j=1}^NP_{\mathrm{\Lambda },N}(\xi _1=j)=1`$. However, if we first perform the thermodynamic limit $`L\mathrm{}`$, $`N\mathrm{}`$ in such a way that $`N/L^d\rho >0`$, and then sum over the limiting probabilities $`P_\rho (\xi _1=j)`$, we can state only that

$$\underset{j=1}{\overset{\mathrm{}}{}}P_\rho (\xi _1=j)1.$$
(5)
We speak about cycle percolation, in an obvious analogy with site or bond percolation, if in this relation the strict inequality holds. In I we conjectured that this happens if and only if there is BEC. Half of the conjecture, namely that BEC implies cycle percolation was shown for the free Bose gas. Now we give a complete proof of a somewhat stronger assertion.

In the theorem below by condensate we mean the ensemble of particles in the one-particle ground state, when BEC occurs. Throughout the paper we work with periodic boundary conditions, so these are particles with zero wavevector. Furthermore, a macroscopic cycle is a cycle of positive density, i.e. a cycle containing a non-vanishing fraction of the total number of particles.

Theorem In the perfect and mean-field Bose gases Bose-Einstein condensation occurs if and only if there is cycle percolation. In the condensate the number of macroscopic cycles is (countable) infinite. Moreover, there are no finite cycles in the condensate and no macroscopic cycles outside it.

In the following section we prove the theorem in three steps. First, we recall the proof presented in I for the existence of cycle percolation when there is condensation. The second and third parts contain the new results: the absence of percolation in the absence of condensation and the analysis of the cycles in the condensate.

## 2 Proof of the Theorem

### 2.1 Proof of cycle percolation in case of condensation

For the ideal Bose gas

$$H_{\mathrm{\Lambda },N}=T_{\mathrm{\Lambda },N}=\frac{\mathrm{}^2}{2m}\underset{i=1}{\overset{N}{}}\mathrm{\Delta }_i,$$
(6)
the N-particle kinetic energy operator with, as we suppose in this paper, periodic boundary conditions in $`\mathrm{\Lambda }`$. In rewriting (1), we make use of the fact that $`U`$ is a representation of $`S_N`$ and, therefore, $`U(g_{\{n_j\}})`$ can be decomposed into a product of $`n_j`$ commuting factors, each corresponding to a cycle. Because now $`\mathrm{exp}(\beta T_{\mathrm{\Lambda },N})`$ also factorizes, the trace will factorize according to the cycles, and the partition function reads

$$Q_{\mathrm{\Lambda },N}=\underset{\{n_j\}:{\scriptscriptstyle jn_j}=N}{}\underset{j=1}{\overset{N}{}}\frac{1}{n_j!}\left(\frac{1}{j}\mathrm{Tr}e^{j\beta T_{\mathrm{\Lambda },1}}\right)^{n_j}$$
(7)
with the trace over $``$. The key to the proof is provided by the bounds

$$\left(\frac{L}{\lambda _B\sqrt{j}}1\right)^d<\mathrm{Tr}e^{j\beta T_{\mathrm{\Lambda },1}}<\left(\frac{L}{\lambda _B\sqrt{j}}+1\right)^d.$$
(8)
These are direct consequences of the inequalities

$$\frac{1}{2}\sqrt{\frac{\pi }{\alpha }}1<_1^{\mathrm{}}e^{\alpha x^2}dx<\underset{n=1}{\overset{\mathrm{}}{}}e^{\alpha n^2}<_0^{\mathrm{}}e^{\alpha x^2}dx=\frac{1}{2}\sqrt{\frac{\pi }{\alpha }}$$
(9)
applied with $`\alpha =\pi j\lambda _B^2/L^2`$ where $`\lambda _B=\mathrm{}\sqrt{2\pi \beta /m}`$, the thermal de Broglie wavelength. Now equation (4) with (7) yields

$$P_{\mathrm{\Lambda },N}(\xi _1=j)=\frac{1}{N}\frac{Q_{\mathrm{\Lambda },Nj}}{Q_{\mathrm{\Lambda },N}}\mathrm{Tr}e^{j\beta T_{\mathrm{\Lambda },1}}<\frac{1}{N}\mathrm{Tr}e^{j\beta T_{\mathrm{\Lambda },1}}.$$
(10)
The inequality holds because $`Q_{\mathrm{\Lambda },Nj}<Q_{\mathrm{\Lambda },N}`$ for periodic boundary conditions. (See (I.6.16) for Dirichlet boundary condition.) Taking the thermodynamic limit and using the upper bound in (8), we find

$$P_\rho (\xi _1=j)\frac{1}{j^{d/2}\rho \lambda _B^d}.$$
(11)
This leads to

$$\underset{j=1}{\overset{\mathrm{}}{}}P_\rho (\xi _1=j)\frac{1}{\rho \lambda _B^d}\underset{j=1}{\overset{\mathrm{}}{}}\frac{1}{j^{d/2}}=\frac{g_{d/2}(1)}{\rho \lambda _B^d}$$
(12)
where

$$g_\alpha (z)=\underset{j=1}{\overset{\mathrm{}}{}}\frac{z^j}{j^\alpha }.$$
(13)
According to the well-known condition (e.g. Huang 1987), there is BEC in the free Bose gas if and only if the right-hand side of (12) is less than 1, which implies cycle percolation.

### 2.2 Proof of the absense of percolation if there is no condensation

Here we refer to the strong equivalence of ensembles in the absence of BEC. As we show below, for any fugacity $`z<1`$ and any fixed $`j`$ the grand-canonical probability distribution of $`n_j/L^d`$ becomes degenerate in the thermodynamic limit, i.e. will be concentrated onto a single value. As a consequence, the canonical probability distribution of $`n_j/L^d`$ will also be degenerate and concentrated on the same value, provided that the particle density $`\rho `$ corresponds to $`z`$. So it will be possible to compute $`P_\rho `$ from the grand-canonical probabilities $`P_z`$. First we show the absence of cycle percolation in the absence of BEC in the grand-canonical ensemble.

The grand-canonical partition function

$`Z_{\mathrm{\Lambda },z}`$ $`=`$ $`{\displaystyle \underset{N=0}{\overset{\mathrm{}}{}}}z^NQ_{\mathrm{\Lambda },N}={\displaystyle \underset{N=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{z^N}{N!}}{\displaystyle \underset{gS_N}{}}\mathrm{Tr}_^NU_N(g)e^{\beta H_{\mathrm{\Lambda },N}}`$
$`=`$ $`{\displaystyle \underset{\{n_j\}_{j1}:{\scriptscriptstyle n_j}<\mathrm{}}{}}\left({\displaystyle \underset{j}{}}\left({\displaystyle \frac{z^j}{j}}\right)^{n_j}{\displaystyle \frac{1}{n_j!}}\right)\mathrm{Tr}_{^{{\scriptscriptstyle jn_j}}}U_{{\scriptscriptstyle jn_j}}(g_{\{n_j\}})e^{\beta H_{\mathrm{\Lambda },{\scriptscriptstyle jn_j}}}`$
can also be considered as the generator of a probability distribution over $`S_N`$. Since any $`gS_N`$ is naturally embedded into $`S_N^{}`$ for any $`N^{}>N`$ (by setting $`g(i)=i`$ for $`i>N`$), it is more precise to say that now probabilities are assigned to the representations $`U_N(g)`$,

$$P_{\mathrm{\Lambda },z}(U_N(g))=\frac{z^N}{N!Z_{\mathrm{\Lambda },z}}\mathrm{Tr}_^NU_N(g)e^{\beta H_{\mathrm{\Lambda },N}}$$
(15)
or, as the second line of (2.2) shows it, to terminating sequences of non-negative integers, each $`\{n_j\}`$ with $`n_j<\mathrm{}`$ corresponding to a conjugation classe of $`S_N`$ for $`N=jn_j`$. As in the case of the canonical ensemble, we will consider $`n_j`$ to be random variables. For a general Hamiltonian any finite number of them has a coupled joint distribution. However, if for the ideal Bose gas we replace (7) into equation (2.2), we arrive at

$$Z_{\mathrm{\Lambda },z}=\underset{j=1}{\overset{\mathrm{}}{}}\mathrm{exp}\left\{\frac{z^j}{j}\mathrm{Tr}e^{j\beta T_{\mathrm{\Lambda },1}}\right\},$$
(16)
showing that the $`n_j`$ are independent random variables with probabilities

$$P_{\mathrm{\Lambda },z}(n_j=n)=\mathrm{exp}\left\{\frac{z^j}{j}\mathrm{Tr}e^{j\beta T_{\mathrm{\Lambda },1}}\right\}\frac{1}{n!}\left(\frac{z^j}{j}\mathrm{Tr}e^{j\beta T_{\mathrm{\Lambda },1}}\right)^n.$$
(17)
The $`L`$-dependence of (17) leads to the asymptotic degeneracy of the grand-canonical distribution of $`n_j/L^d`$ for $`z<1`$:

$$\underset{L\mathrm{}}{lim}P_{\mathrm{\Lambda },z}(\frac{n_j}{L^d}x)=\mathrm{\Theta }(x\rho _j)$$
(18)
where $`\mathrm{\Theta }`$ is the Heaviside function and

$$\rho _j=\frac{z^j}{j^{1+d/2}\lambda _B^d}.$$
(19)
Equation (18) is the manifestation of the following large deviation principle. Let $`0<a<b`$, then

$$\underset{L\mathrm{}}{lim}\frac{1}{L^d}\mathrm{ln}P_{\mathrm{\Lambda },z}(a\frac{n_j}{L^d}b)=\underset{axb}{inf}I_j(x).$$
(20)
where the rate function $`I_j`$ is given by

$$I_j(x)=\rho _jx+x\mathrm{ln}\frac{x}{\rho _j}.$$
(21)
To prove (20) and (18), we extend the right-hand-side of (17) to real positive values of $`n`$ and define

$$f_L(x)=\mathrm{exp}\left\{\frac{z^j}{j}\mathrm{Tr}e^{j\beta T_{\mathrm{\Lambda },1}}\right\}\frac{1}{(xL^d)!}\left(\frac{z^j}{j}\mathrm{Tr}e^{j\beta T_{\mathrm{\Lambda },1}}\right)^{xL^d}$$
(22)
for $`x>0`$. Then, by using (8) and Stirling’s formula, asymptotically, as $`L\mathrm{}`$, we find

$$f_L(x)\frac{1}{\sqrt{2\pi xL^d}}e^{L^dI_j(x)}.$$
(23)
Equations (18) and (20) follow from (23) by noting that $`I_j(x)`$ is strictly convex for $`x>0`$ and $`I_j(\rho _j)=0`$ is its minimum. For the distribution of $`n_j`$ this implies that asymptotically $`n_j`$ is concentrated to an $`O(L^{d/2})`$-neighbourhood of $`\rho _jL^d`$.

As in the canonical case, $`\xi _1`$ is the length of the cycle containing 1. In finite volumes

$`P_{\mathrm{\Lambda },z}(\xi _1=j)`$ $`=`$ $`{\displaystyle \underset{\{m_i\}:{\scriptscriptstyle m_i}<\mathrm{}}{}}P_{\mathrm{\Lambda },z}(\xi _1=j|n_i=m_i,i1)P_{\mathrm{\Lambda },z}(n_i=m_i,i1)`$ (24)
$`=`$ $`{\displaystyle \underset{\{m_i\}:{\scriptscriptstyle m_i}<\mathrm{}}{}}{\displaystyle \frac{jm_j/L^d}{im_i/L^d}}{\displaystyle \underset{i}{}}P_{\mathrm{\Lambda },z}(n_i/L^d=m_i/L^d).`$
With our large-deviation result then

$$P_z(\xi _1=j)=\underset{L\mathrm{}}{lim}P_{\mathrm{\Lambda },z}(\xi _1=j)=\frac{j\rho _j}{_ii\rho _i}=\frac{z^j}{j^{d/2}g_{d/2}(z)}.$$
(25)
Since the denominator is finite for $`z<1`$, we find $`_{j=1}^{\mathrm{}}P_z(\xi _1=j)=1`$ indeed. Now for $`z<1`$ $`\rho `$ and $`z`$ are related through

$$g_{d/2}(z)=\rho \lambda _B^d.$$
(26)
Substituting this expression into (25) we obtain $`P_\rho `$ and the absence of cycle percolation in the canonical ensemble.

Two remarks are in order here:

1. As we have mentioned at the beginning of this section, the asymptotic degeneracy (18) implies that the canonical distribution of $`n_j/L^d`$ is also asymptotically degenerate. Knowing this, we can derive $`P_\rho (\xi _1=j)`$ without passing through $`P_z(\xi _1=j)`$. Indeed, in equation (4) for $`j=1,\mathrm{},N`$ the averages $`n_j_{\mathrm{\Lambda },N}`$ can be approximated by the set $`\{n_j\}`$ belonging to the largest term of (7), which can be found by conditional maximization. The asymptotically valid result is

$$n_j\frac{Nz_N^j}{j^{1+d/2}\rho \lambda _B^d}$$
(27)
Here $`z_N`$ appears via a Lagrange multiplier and satisfies

$$\underset{j=1}{\overset{N}{}}\frac{z_N^j}{j^{d/2}}=\rho \lambda _B^d.$$
(28)
If $`\rho <\rho _c(\beta )=g_{d/2}(1)/\lambda _B^d`$ then $`z_N<1`$. Writing equation (28) for $`N`$ and $`N+1`$ and taking the difference one can see that for any fixed $`\rho >0`$, $`z_N>z_{N+1}`$. So for $`\rho <\rho _c`$, $`z_Nz<1`$ which, hence, is just the solution for $`z`$ of equation (26). In sum,

$$P_{\mathrm{\Lambda },N}(\xi _1=j)\frac{z_N^j}{j^{d/2}\rho \lambda _B^d}P_\rho (\xi _1=j)=\frac{z^j}{j^{d/2}\rho \lambda _B^d}.$$
(29)
2. Comparing the expression (25) of $`P_{z=1}(\xi _1=j)`$ with the right member of the inequality (11), we see that the inequality saturates at $`\rho =\rho _c`$. In fact, for $`\rho \rho _c`$ there is equality in (11): Although for $`\rho >\rho _c`$ there is no asymptotic equivalence between the canonical and grand canonical distributions of $`n_j/L^d`$ (reflecting the fact that the canonical and grand canonical Gibbs states are related through the nondegenerate Kac density, see Cannon 1973 and Lewis and Pulè 1974), the former still becomes asymptotically degenerate. So the strong equivalence of the canonical ensemble with a microcanonical one (with $`n_j`$ fixed) persists. The largest term of (7) is given by (27). From (28) it follows that $`z_N>1`$ and converges to 1, therefore

$$P_\rho (\xi _1=j)=\frac{1}{j^{d/2}\rho \lambda _B^d}\mathrm{if}\rho >\rho _c.$$
(30)
### 2.3 Existence of an infinity of macroscopic cycles

While $`_{j=1}^{\mathrm{}}P_\rho (\xi _1=j)<1`$ means that infinite cycles occur with positive probability, this inequality does not inform us about the asymptotic size and number of large cycles. Below we show that in the domain of condensation macroscopic cycles are in abundance, their number is infinite, and each of them is associated with the condensate.

Let, therefore, $`d3`$ and fix a $`\rho >\rho _c(\beta )=g_{d/2}(1)/\lambda _B^d`$. Thus, $`\rho _0\rho \rho _c>0`$. Let, moreover, $`N=N(L,\rho )=[\rho L^d]`$. Writing the trace in the basis of the eigenstates of $`T_{\mathrm{\Lambda },N}`$, the partition function reads

$$Q_{\mathrm{\Lambda },N}=\frac{1}{N!}\underset{K}{}e^{\beta _kN_k(K)\epsilon _k}\underset{g}{}\delta _{gK,K}.$$
(31)
Here $`K=(k_1,\mathrm{},k_N)`$, $`k_i\mathrm{\Lambda }^{}=(2\pi /L)^d`$, $`\epsilon _k=(\mathrm{}^2k^2/2m)`$, $`N_k(K)`$ is the number of occurrence of $`k`$ in the sequence $`K`$ and $`gK=(k_{g^1(i)})_{i=1}^N`$. From (31) we can infer the probability of a pair $`(g,K)`$,

$$P_{\mathrm{\Lambda },N}(g,K)=\frac{e^{\beta _kN_k(K)\epsilon _k}\delta _{gK,K}}{N!Q_{\mathrm{\Lambda },N}}$$
(32)
and that of $`K`$, $`P_{\mathrm{\Lambda },N}(K)=_gP_{\mathrm{\Lambda },N}(g,K)`$.

The asymptotic distribution of the random variables $`N_k`$ was studied in detail by Buffet and Pulè (1983). They showed that the canonical distribution of $`N_k/L^d`$ becomes asymptotically degenerate, namely it converges to the Dirac delta concentrated at $`\rho _0`$ if $`k=0`$ and at 0 if $`k0`$. From their analysis it follows that $`\xi _1=O(N)`$ is possible only if $`k_1=0`$. Moreover, let

$$A_{L,ϵ}=\left\{K(\mathrm{\Lambda }^{})^N\right||N_0(K)/L^d\rho _0|<ϵ,N_k(K)/L^d<ϵ\text{ for }k0\},$$
(33)
then

$$\underset{L\mathrm{}}{lim}P_{\mathrm{\Lambda },N}(A_{L,ϵ})=1\text{ for any }ϵ>0.$$
(34)
We prove equation (34) by showing that the probability of the complement of $`A_{L,ϵ}`$ goes to zero.

$$P_{\mathrm{\Lambda },N}(A_{L,ϵ}^c)P_{\mathrm{\Lambda },N}(A_{L,ϵ}(k=0)^c)+\underset{k0}{}P_{\mathrm{\Lambda },N}(A_{L,ϵ}(k)^c)$$

where

$$A_{L,ϵ}(k=0)=\{K:|N_0(K)/L^d\rho _0|<ϵ\}$$

and for $`k0`$

$$A_{L,ϵ}(k)=\{K:N_k(K)/L^d<ϵ\}.$$

Now $`P_{\mathrm{\Lambda },N}(A_{L,ϵ}(k=0)^c)`$ goes to zero as $`L`$ tends to infinity. With $`x_k=e^{\beta \epsilon (k)}`$

$`P_{\mathrm{\Lambda },N}(A_{L,ϵ}(k)^c)={\displaystyle \frac{_{m=ϵL^d}^Nx_k^m_{\{n_q\}_{qk}:{\scriptscriptstyle n_q}=Nm}_{qk}x_q^{n_q}}{_{\{n_q\}:{\scriptscriptstyle n_q}=N}_qx_q^{n_q}}}`$
$`x_k^{ϵL^d}{\displaystyle \frac{_{m=ϵL^d}^N_{\{n_q\}_{qk}:{\scriptscriptstyle n_q}=Nm}_{qk}x_q^{n_q}}{_{\{n_q\}:{\scriptscriptstyle n_q}=N}_qx_q^{n_q}}}x_k^{ϵL^d}`$
because in the numerator we have a part of the sum of the denominator, over those sets of occupation numbers for which $`n_q=N`$, $`n_k=0`$ and $`n_0ϵL^d`$ (observe that $`x_0=1`$). But

$$\beta \epsilon (k)=\frac{\beta \mathrm{}^2k^2}{2m}=\pi \lambda _B^2n^2/L^2$$

where $`n=(n_1,\mathrm{},n_d)`$, $`n^2=_{i=1}^dn_i^2`$ with $`n_i`$ integers. Therefore

$`{\displaystyle \underset{k0}{}}P_{\mathrm{\Lambda },N}(A_{L,ϵ}(k)^c)`$ $``$ $`{\displaystyle \underset{nZ^d\{0\}}{}}\mathrm{exp}\{\pi ϵ\lambda _B^2L^{d2}n^2\}`$
$`=`$ $`\left({\displaystyle \underset{i=\mathrm{}}{\overset{\mathrm{}}{}}}\mathrm{exp}\{\pi ϵ\lambda _B^2L^{d2}i^2\}\right)^d1`$
$`=`$ $`\left[1+2{\displaystyle \underset{i=1}{\overset{\mathrm{}}{}}}\mathrm{exp}\{\pi ϵ\lambda _B^2L^{d2}i^2\}\right]^d1`$
$``$ $`\left[1+(\pi ϵ\lambda _B^2L^{d2})^{1/2}\right]^d10`$
with $`L`$ going to infinity if $`d>2`$.

The crucial property of the probabilities (32) is that they are the same for all permutations leaving $`K`$ invariant. Due to this uniform distribution and to (34), the discussion of the ground state cycle percolation in I becomes relevant and we can derive the analogues of equations (I.5.16) and (I.5.18) (notice the different meaning of $`d`$ – the number of spin states – in I). Related classical results can also be found in Feller 1968.

We fix $`ϵ<\rho _0/2`$ and choose a $`K`$ in $`A_{L,ϵ}`$. Then $`N_0(K)>N_k(K)`$ for every $`k0`$. We consider a $`j`$ such that $`N_0(K)j>\mathrm{max}_{k0}N_k(K)`$. Thus, cycles of length $`j`$ can occur only among particles with wavevector 0. The conditional expectation value of $`n_j`$, provided that $`K`$ is given, is

$$n_j_K\underset{m}{}mP_{\mathrm{\Lambda },N}(n_j=m|K)=\frac{\underset{m}{}m|\{gS_N|gK=K,n_j=m\}|}{|\{gS_N|gK=K\}|}.$$
(35)
Here $`|\{\mathrm{}\}|`$ means the number of elements of $`\{\mathrm{}\}`$. Now the number of permutations leaving $`K`$ invariant is $`_kN_k(K)!`$, so we obtain

$$n_j_K=n_j_{N_0(K)}\frac{1}{N_0(K)!}\underset{m}{}m|\{gS_{N_0(K)}|n_j=m\}|=\frac{1}{j}.$$
(36)
Indeed, for any $`Mj`$ the uniform distribution $`P_M(g)=1/M!`$ over $`S_M`$ gives rise to

$$n_j_M=\underset{m}{}mP_M(n_j=m)=\frac{1}{M!}\underset{m}{}m|\{gS_M|n_j=m\}|=\frac{1}{j};$$
(37)
in order to see this it suffices to recall that the probability of a conjugation class $`\{n_i\}`$ of $`S_M`$ ($`_iin_i=M`$) is

$$P_M(\{n_i\})=\underset{i=1}{\overset{M}{}}\frac{1}{n_i!}\left(\frac{1}{i}\right)^{n_i}.$$
(38)
From (36) we obtain

$$P_{\mathrm{\Lambda },N}(\xi _1=j|K)=\frac{jn_j_{N_0(K)}}{N_0(K)}\delta _{k_1,0}=\frac{\delta _{k_1,0}}{N_0(K)}.$$
(39)
Only particles with the same wave vector can be in the same cycle. Therefore, $`k_1=0`$ implies $`\xi _1N_0(K)`$ and

$$P_{\mathrm{\Lambda },N}(\xi _1>j|K)=\left(1\frac{j}{N_0(K)}\right)\delta _{k_1,0}.$$
(40)
Fix now $`0<x\rho _0/\rho `$. Equations (34) and (40) imply

$$\underset{L\mathrm{}}{lim}P_{\mathrm{\Lambda },N}\left(\frac{\xi _1}{N}>x\right)=\frac{\rho _0}{\rho }x.$$
(41)
Indeed, by choosing $`ϵ<x\rho `$,

$`\underset{L\mathrm{}}{lim}P_{\mathrm{\Lambda },N}\left({\displaystyle \frac{\xi _1}{N}}>x\right)`$ $`=`$ $`\underset{L\mathrm{}}{lim}{\displaystyle \underset{KA_{L,ϵ}}{}}P_{\mathrm{\Lambda },N}\left({\displaystyle \frac{\xi _1}{N}}>x|K\right)P_{\mathrm{\Lambda },N}(K)`$ (42)
$`=`$ $`\underset{L\mathrm{}}{lim}{\displaystyle \underset{KA_{L,ϵ}}{}}\left(1{\displaystyle \frac{xN}{N_0(K)}}\right)\delta _{k_1,0}P_{\mathrm{\Lambda },N}(K)`$
$``$ $`\left(1{\displaystyle \frac{x\rho }{\rho _0\pm ϵ}}\right)\underset{L\mathrm{}}{lim}{\displaystyle \underset{KA_{L,ϵ}}{}}\delta _{k_1,0}P_{\mathrm{\Lambda },N}(K)`$
$`=`$ $`\left(1{\displaystyle \frac{x\rho }{\rho _0\pm ϵ}}\right)P_\rho (k_1=0)`$
$`=`$ $`\left(1{\displaystyle \frac{x\rho }{\rho _0\pm ϵ}}\right){\displaystyle \frac{\rho _0}{\rho }}.`$
The upper and lower bounds holding for all $`ϵ<x\rho `$, we can send $`ϵ`$ to zero and find (41).

Complementary to this result is

$$P_\rho (\xi _1=j|k_1=0)\underset{L\mathrm{}}{lim}P_{\mathrm{\Lambda },N}(\xi _1=j|k_1=0)=0\text{for any fixed }j$$
(43)
whose analogue also existed in the ground-state cycle percolation, see (I.5.13). Equation (43) can be obtained with the following argument. If we restrict the probability distribution $`P_{\mathrm{\Lambda },N}(g|K)`$ to the permutations of the $`N_0(K)`$ particles with 0 wavevector, the restricted distribution is still uniform. Therefore, by (37), equation (36) extends to cycles of an arbitrary length among the zero wavevectors. As a consequence, if $`KA_{L,ϵ}`$ is chosen so that $`k_1=0`$, equation (39) also extends to all $`jN_0(K)`$ and together with (34) implies (43). Thus, the full probability (30) comes from $`k_10`$.

The number of infinite cycles of positive density is (countable) infinite. The number of those with density larger than $`x\frac{\rho _0}{\rho }`$ can be obtained from equations (34) and (36):

$$\underset{L\mathrm{}}{lim}\underset{j=xN}{\overset{N_0}{}}n_j_{\mathrm{\Lambda },N}=\mathrm{ln}\frac{\rho _0}{\rho x}.$$
(44)
This number can be arbitrarily large if $`x`$ is sufficiently small. For $`m>\mathrm{ln}\rho \mathrm{ln}\rho _0`$ the expected number of cycles of density between $`e^{(m+1)}`$ and $`e^m`$ is

$$\mathrm{ln}\frac{e^{m+1}\rho _0}{\rho }\mathrm{ln}\frac{e^m\rho _0}{\rho }=\mathrm{ln}e=1.$$

The intervals $`[e^{(m+1)},e^m)`$ are disjoint, their number is infinite and on average there belongs one infinite cycle to each interval.

Finally, we note that the above results apply without any further ado to the imperfect Bose gas (cf Huang 1987). The Hamiltonian of this latter reads

$$H_{\mathrm{\Lambda },N}=T_{\mathrm{\Lambda },N}+\frac{a}{2L^d}N^2$$
(45)
with some $`a>0`$, so probabilities and averages in the canonical ensemble coincide with those of the perfect Bose gas.

Acknowledgment This work was supported by the Hungarian Scientific Research Fund under grant no T 30543.

References
Buffet E and Pulè J V 1983 J. Math. Phys. 24 1608
Bund S and Schakel A M J 1999 Mod. Phys. Lett. B 13 349
Cannon J T 1973 Commun. Math. Phys. 29 89
Feller W 1968 An introduction to probability theory and its applications I.
——— (New York: John Wiley) Ch. X. 6
Feynman R P 1953 Phys. Rev. 91 1291
Huang K 1987 Statistical Mechanics (New York: Wiley)
Lewis J T and Pulè J V Commun. Math. Phys. 36 1
Martin Ph 2001 private communication
Schakel A M J 2001 Phys. Rev. E 63 026115
Sütő A 1993 J. Phys. A: Math. Gen. 26 4689
Ueltschi D 2002 In and out of equilibrium: Physics with a probability flavor
———– ed V Sidoravicius (Basel: Birkhäuser) p 363


text:
# KiloHertz quasi–periodic oscillations difference frequency exceeds inferred spin frequency in 4U 1636–53

## 1 Introduction

In recent years kiloHertz quasi–periodic oscillations (kHz QPOs) have been found in power spectra of more than 20 low–mass X–ray binaries (LMXBs) using observations made with the Rossi X–ray Timing Explorer (RXTE) satellite (?). In most sources the kHz QPOs are found in pairs. The high frequencies of these QPOs suggest that at least one of the two peaks reflects the orbital motion of matter close to the neutron star. Another high frequency phenomenon was detected with RXTE: during some type I X–ray bursts nearly coherent oscillations which slightly drift in frequency, the so–called burst oscillations, have been observed in the power spectra of ten sources (for a review see ?). The X–ray flux–modulation is consistent with being due to the changing aspect of a drifting hot spot on the surface of the neutron star. Therefore, these burst oscillations are thought to reflect the spin frequency of the neutron star. In high luminosity ‘Z’–source LMXBs QPOs at low frequencies were found with the EXOSAT and GINGA satellites (for a review see ?). Observations with the RXTE satellite have also revealed low–frequency QPOs in a number of atoll sources (e.g. ?; ?)

The atoll source 4U 1636–53 is one of the few sources that displays all of the high–frequency QPO phenomena that have been observed in LMXBs; two kHz QPOs have been observed (?; ?) as well as a sideband to the lower of the two kHz QPOs (?) and burst oscillations at $``$581 Hz (?; ?). ? presented evidence that these oscillations are in fact the second harmonic of the neutron star spin of $``$290.5 Hz. However, using another dataset these findings were not confirmed (?).

Since the frequency difference between the two kHz QPO peaks (the peak separation, $`\mathrm{\Delta }\nu `$) is nearly equal to (half) the burst oscillation frequency a beat frequency mechanism was proposed for the kHz QPOs (?). In such a model the higher frequency QPO is attributed to the orbital frequency of clumps of plasma at a special radius near the neutron star while the lower frequency QPO is due to the beat between the orbital frequency and the stellar spin frequency. A specific model incorporating the beat frequency mechanism, the sonic–point beat frequency model, is now one of the leading models explaining kHz QPOs in LMXBs (?). However, subsequent findings complicated the picture. As was first shown in Sco X–1, the peak separation decreases as the kHz QPO frequencies increase (?). Furthermore, in 4U 1636–53 the peak separation was found to be less than half the burst oscillation frequency (?). These findings could be explained within the sonic point beat frequency model (?) by taking into account the effect of the inward velocity of the plasma on the emerging kHz QPO frequencies. A firm prediction of this revised sonic–point beat frequency model is that for prograde spinning accretion disks the kHz QPO peak separation is smaller than the neutron star spin frequency. This follows from the detailed description of the plasma flow patterns in this model which produce a number of corrections to the observed frequencies relative to the true orbital and beat frequency, all of which are expected to make $`\mathrm{\Delta }\nu `$ smaller (?). Indeed, until now, in various sources the peak separation was always found to be consistent with, or less than (half) the burst oscillation frequency.

In this Letter we show that in recent observations of 4U 1636–53 the kHz QPO peak separation is significantly larger than half the burst oscillation frequency. Using other observations ? had shown earlier that the peak separation was significantly smaller than half the burst oscillation frequency. This means that in order for the sonic–point beat frequency to explain both observations further changes will have to be made in the description of the flow pattern. The only change in the flow pattern previously described that could produce the observed change in the separation frequency is that in which the accretion disk changes from prograde to retrograde when the lower kHz QPO moves from $``$750 Hz to $``$800 Hz.

## 2 Observations and analysis

4U 1636–53 was observed in September 2001 and January 2002 with the proportional counter array (PCA; ?) onboard the RXTE satellite (?). Due to the reduced duty cycle of the PCA detectors only a subset of the 5 detectors is operational for most of the time. In this analysis we only used data in which 4 of the 5 detectors were operational. A log of the observations is given in Table 1. In total $``$28.4 ksec of data were used. Data were obtained using an Event–mode in which information on detected X–ray photons is stored onboard and telemetered to the ground on an event–to–event basis. The time resolution of the data is $``$125$`\mu `$s and RXTE’s effective energy range (2–60 keV) is covered by 64 channels.

Using this data we calculated power density spectra of segments of 64 s up to a Nyquist frequency of 4096 Hz in one energy band covering the total PCA energy range. We used a dynamical power spectrum displaying consecutive power spectra rebinned to a frequency resolution of 2 Hz and a time resolution of 128 s to visualize the time evolution of the lower kHz QPO for all observations except for the power spectra of observation 60032-01-21-00 which were rebinned to a time resolution of 162s (segment 7) and 256s (segment 8), and the power spectra of observation 60032-01-22-00, which were treated separately. We traced the frequency of the lower kHz QPO and fitted the power spectrum in a window of 200 Hz centered on the traced frequency with a Lorentzian and a constant. All peaks in each power spectrum of 128 s, 162 s, and 256 s had a significance of $`3\sigma `$ or more. The frequency of the bin associated with the peak of the fitted Lorentzian was taken to be the frequency of the QPO in that power density spectrum. Finally, we used the shift–and–add method described by ? to shift the lower kHz QPO to the same (arbitrary) frequency in each power spectrum before averaging. This assures that at least for timescales longer than 256 s the determined peak separation is not influenced by the fact that the two kHz QPOs have a different frequency–amplitude relation (?).

We averaged the aligned power spectra and fitted the average power spectrum from 256–1500 Hz with a constant, to represent the power introduced by the Poisson counting noise, and two Lorentzians to represent the two kHz QPOs. Errors on the fit parameters were determined using $`\mathrm{\Delta }\chi ^2=1.0`$ ($`1\sigma `$ single parameter). The observation 60032-01-22-00 was treated separately since we were not able to find a more than 3$`\sigma `$ lower kHz QPO using data stretches shorter than or equal to 256 s. Instead we averaged the power spectra in this observation without shifting the data. We note that by doing so we may have artificially introduced deviations in the kHz QPO peak separation if the two kHz QPOs moved in frequency during the observation for reasons outlined above.

The fits of the average power spectra were good, with a reduced $`\chi ^2`$ slightly larger than 1 for 146 degrees of freedom. We found that, combining all observations, except 60032-01-22-00, the frequency separation between the two kHz QPO peaks was 323.3$`\pm `$4.3 Hz, significantly larger than 290.5 Hz the putative spin frequency according to an analysis of ?, but still significantly smaller than 581 Hz, the frequency of the strong burst oscillations found by ? and ? (see Figure 1). We then divided the data set into three parts according to the frequency of the lower kHz QPO and measured the peak separation in each of the parts. For each part we also averaged the power spectra without shifting in order to measure the average frequency of the lower kHz QPO and that of a low–frequency QPO changing in frequency from $``$24–42 Hz. In Table 2 we present our results on the kHz QPOs. The fractional rms amplitude of the low–frequency QPO changed from 7.8$`\pm `$0.6 per cent to 4.3$`\pm `$0.7 per cent (2–60 keV) while its frequency increased from 23.6$`\pm `$2.4 Hz to 42$`\pm `$4 Hz and its FWHM was consistent with being constant at 26$`\pm `$9 Hz. In Figure 2 we plot the measured peak separation ($`\mathrm{\Delta }\nu `$) as a function of the frequency of the lower and upper kHz QPO (left and right panel, respectively). The squares are our new measurements using the shift–and–add method, the diamond is our measurement using observation 60032-01-22-00, and the dots are the measurements of ? based on RXTE data obtained from 1996–2000 (see also ?; ?). Clearly, $`\mathrm{\Delta }\nu `$ changes from more than half the burst oscillation frequency to less than that. The average $`\mathrm{\Delta }\nu `$ of the three leftmost points is 327.2$`\pm `$4.2 Hz and that of the dots is 247.1$`\pm `$1.9 Hz making the jump in $`\mathrm{\Delta }\nu \mathrm{\hspace{0.17em}70.1}\pm 4.6`$ Hz.

## 3 Discussion

We found that in 4U 1636–53 the kHz QPO frequency separation ($`\mathrm{\Delta }\nu `$) changes as a function of the frequency of the lower kHz QPO from 323.3$`\pm `$4.3 Hz to 242.4$`\pm `$3.6 Hz (Fig. 2 left); the maximum peak separation is significantly larger than 290.5 Hz (half the burst oscillation frequency of $``$581 Hz; see Section 1), while the minimum peak separation is significantly smaller than that. The change in $`\mathrm{\Delta }\nu `$ as a function of upper kHz QPO frequency is abrupt; the data are consistent with a sudden jump at $`\nu =1050`$ Hz (Fig. 2 right). We note that similar jumps in $`\mathrm{\Delta }\nu `$ as a function of the upper kHz QPO have also been seen in 4U 1728–34 and 4U 1608–52, although in those cases $`\mathrm{\Delta }\nu `$ always remained smaller or equal to (half) the burst oscillation frequency (?; ?; ?; ?). This is the first time that the kHz peak separation has been shown to be significantly larger than the inferred neutron star spin frequency and also the first time that $`\mathrm{\Delta }\nu `$ has been seen to vary between less and more than (half) the burst oscillation frequency.

In previous work burst oscillations detected in the dipping burster 4U 1916–05 at $``$272 Hz also seemed to occur at too low a frequency for $`\mathrm{\Delta }\nu `$ (?). However, in that source the kHz QPOs are weak and sometimes very broad, therefore whole segments of several ksec of data had to be averaged to obtain a detection (?). This means that there is no certainty that the kHz QPOs were detected simultaneously or that their apparent separation was not affected by frequency shifts combined with changes in their power ratio. Therefore, the measured $`\mathrm{\Delta }\nu `$ may not reflect the true peak separation. In the current work, using the shift–and–add method we excluded the possibility that effects changing $`\mathrm{\Delta }\nu `$ on timescales longer than 256 s have influenced the observed $`\mathrm{\Delta }\nu `$. On timescales shorter than 256 s, changes in the rms amplitudes could bias the peak separation for reasons outlined above. However, since our measurements show that the lower kHz QPO rms amplitude increased while the rms amplitude of the upper kHz QPO decreased this effect will have made the observed $`\mathrm{\Delta }\nu `$ smaller than the true peak separation assuming the rms amplitude–frequency relation is the same for timescales shorter than 256 s. Fast changes in the upper kHz QPO frequency, associated with the 5–6 Hz QPO in Sco X–1 ($``$20 Hz; ?) could also bias $`\mathrm{\Delta }\nu `$ measurements. However, those $``$5–6 Hz QPOs are only found when atoll and Z sources are in a high–luminosity state (e.g. ?; ?). Changes in the kHz QPO frequencies associated with milliHz QPOs present in 4U 1636–53 and 4U 1608–52 in a small luminosity interval (?) are so small (less than a Hz; ?) that it is unlikely that they alter our findings significantly, although we note that such a milliHz QPO was present during some of the observations.

Systematic variations in the profile of the kHz QPOs could in principle cause the kHz QPO separation frequency to become large at low kHz QPO frequencies and small at high frequencies since the centroid frequencies have been measured assuming a Lorentzian peak profile. Since our fits were good (reduced $`\chi ^2`$1) and no systematic residuals were apparent, we estimate that shifts in the centroid frequency due to asymmetries in the peak profiles could at most be as large as the FWMH of the peak divided by its significance, i.e. negligible in case of the lower kHz QPO but up to $``$10 Hz in case of the upper kHz QPO. For power spectra with a lower kHz QPO in the frequency range 650–750 Hz, the average data$``$model residuals within one FWHM are less then 0.15 per cent of the power in the peak. Furthermore, there is no evidence for a sudden change of the QPO profile which could produce the observed abrupt decrease in $`\mathrm{\Delta }\nu `$. When using Gaussians to model the kHz QPO peaks the fits were also good (reduced $`\chi ^2`$1) and we obtained results consistent with those using Lorentzians. We conclude that it is unlikely that our results can be explained by changing asymmetries in the peak profile of the kHz QPOs.

The principal motivation for a beat frequency model is the relative closeness of the inferred spin frequency to $`\mathrm{\Delta }\nu `$. We found that $`\mathrm{\Delta }\nu `$ is consistent with being distributed symmetrically around half the burst oscillation frequency, emphasizing the apparent connection between the inferred spin and $`\mathrm{\Delta }\nu `$. The present formulation of the sonic–point beat–frequency model can explain the observed decrease in $`\mathrm{\Delta }\nu `$ as a function of kHz QPO frequency (?). However, all versions of the sonic–point beat–frequency model described so far involving a prograde spinning accretion disk (?; ?) predict that $`\mathrm{\Delta }\nu `$ is always less than or equal to the spin frequency of the neutron star for Keplerian orbital frequencies larger than the neutron star spin frequencies; since our results for 4U 1636–53 show that $`\mathrm{\Delta }\nu `$ can be larger than the neutron star spin frequency this prediction is now falsified.

? showed that it is still unclear whether the burst oscillations at $``$581 Hz reflect the spin frequency of the neutron star or its second harmonic as suggested by ?. Therefore, one might argue that the spin frequency is really at 581 Hz (not 290.5 Hz) and hence that $`\mathrm{\Delta }\nu `$ is less than the spin frequency. This would, however, take away the principal motivation of beat–frequency models, since the neutron star spin frequency and the kHz QPO peak separation are in that case completely different.

That in the sonic–point beat–frequency model the neutron star spin and $`\mathrm{\Delta }\nu `$ are not exactly equal is explained by the way in which the beat–frequency interaction produces the observed frequencies. The key element is the radial motion of the plasma clumps near the sonic radius (?). This makes the observed lower kHz QPO frequency somewhat larger than the beat–frequency since it squeezes the spatial separation between peaks of enhanced mass flow. It also makes the upper kHz QPO frequency smaller than the orbital frequency since it makes the clump’s footpoint move backward in a corotating frame. It might be possible to find other effects in the plasma flow patterns besides those discussed by ? (all of which make $`\mathrm{\Delta }\nu `$ smaller) which may lead to larger values of $`\mathrm{\Delta }\nu `$ at lower kHz QPO frequencies. Fully relativistic simulations of the gas flow and radiation transport reported by ? have not, so far, shown evidence for this. If, however, the sense of rotation of (part of) the accretion disk would change from prograde to retrograde when $`\mathrm{\Delta }\nu `$ changes from less to more than the spin frequency, then our results could still be in accordance with current formulations of the sonic–point beat–frequency model. This also requires that in our observations for $`\mathrm{\Delta }\nu `$ larger than 290.5 Hz the lower kHz QPO reflects the (near) Keplerian orbital frequency, whereas for those points in Figure 2 where $`\mathrm{\Delta }\nu `$ is smaller than half the burst oscillation frequency, the Keplerian orbital motion is represented by the upper kHz QPO.

We detected a low–frequency QPO between $``$25–40 Hz in 4U 1636–53 simultaneously with the kHz QPO pair. The relation between this low–frequency QPO and the lower kHz QPO is similar to the one found by ? for 4U 1728–34. It has been suggested that such low–frequency QPOs in atoll sources are related to the Z source Horizontal Branch Oscillations (HBO; e.g. ?; ?). According to the magnetospheric beat–frequency model these HBO originate near the magnetospheric radius (?; ?). The frequency of the low–frequency QPO in 4U 1636–53 increased as the frequency of the lower kHz QPO increased. Therefore, if the disk near the sonic radius is indeed retrograde within our observations the retrograde annulus must extend at least up to the magnetospheric radius if the sonic–point and magnetospheric beat–frequency model both apply. Partially retrograde spinning disks have been proposed as an explanation for observed torque reversals (?) in slowly rotating neutron stars with a dipole magnetic field strength of 1011-12 Gauss (?). Modelling by ? shows that the accretion disk can get inclinations of more than 90 degrees. It is not clear if their work is also applicable to the case of low–magnetic field LMXBs. Clearly, allowing disks with alternating rotation senses complicates the sonic–point beat frequency interpretation of the kHz QPOs considerably.

## Acknowledgments

PGJ is supported by EC Marie Curie Fellowship HPMF–CT–2001–01308. MK is supported in part by a Netherlands Organization for Scientific Research (NWO) grant. We would like to thank Cole Miller for useful discussions and the referee for comments which helped improve the manuscript.