text:
# RX J004717.4-251811: The first eclipsing X-ray binary outside the Local GroupThis work is based on observations obtained with XMM-Newton, an ESA Science Mission with instruments and contributions directly funded by ESA Member States and the USA (NASA).

## 1 Introduction

Luminous X-ray binaries (XRBs) in the Milky Way were among the first X-ray sources detected and optically identified. They are classified as low mass X-ray binary (LMXBs) and high mass X-ray binary (HMXBs) systems according to the mass of the donor star. The accreting collapsed object can be a neutron star or a black hole in either class.

XRBs may show orbital variability (eclipses, absorption dips, regular outbursts), pulsations due to rotation of the compact object (neutron star), long term variability (which may be explained due to precession of an accretion disk or free precession of a neutron star) and flares. Examples for eclipsing XRBs are Her X-1, LMC X-4, SMC X-1 with pulsation periods of 1.24 s, 0.71 s, 13.5 s, orbital periods of 1.7 d, 1.4 d, 3.9 d, full X-ray eclipse durations of at least 0.137, 0.144, 0.147 in phase, long term periods of 35 d, 30 d, 45 d – 60 d, and showing no flares, flares, no flares, respectively (Liu et al. 2000, and references therein). Eclipsing XRBs up to now were only detected within the Local Group and the total number remains small. The first examples of such systems (Cen X-3 and Her X-1: Schreier et al. 1972b; Tananbaum et al. 1972, respectively), were detected in the Milky Way with the UHURU satellite, then in the SMC (SMC X-1: Schreier et al. 1972a) and LMC (LMC X-4: Li et al. 1978; White 1978), and finally one in M33 (Peres et al. 1989). Combining X-ray pulse arrival time information and orbital parameters of pulsating eclipsing XRBs with radial-velocities of the optical companion were essential to determine neutron star masses (see e.g. van Paradijs & McClintock 1995).

In this paper we report on the first detection of an eclipsing XRB outside the Local Group in NGC~253. This prototypical nearby edge-on starburst galaxy was observed by many X-ray satellites. X-ray point sources could be resolved with the high spatial resolution instruments aboard Einstein, ROSAT, Chandra, and XMM-Newton. Fabbiano & Trinchieri (1984, hereafter FT84) reported complex X-ray emission from eight point sources as well as diffuse emission from nuclear area, plume and disk using the Einstein HRI detector. Vogler & Pietsch (1999, hereafter VP99) analyzed X-ray point sources in the bulge, disk and halo of NGC~253 in ROSAT PSPC and HRI observations. The number of point-like X-ray sources increased to 30 sources, 13 of which are time variable. RX~J004717.4-251811 (E2 of FT84, X17 of VP99) is one of the four Einstein sources detected again with ROSAT. During most of the ROSAT observations a second source of similar brightness (X18) was visible just 15″ to the south. RX~J004717.4-251811 shows strong time variability within the ROSAT observations and also compared to the Einstein HRI detection. This time variability and the ROSAT PSPC hardness ratio classification made RX~J004717.4-251811 a good XRB candidate.

NGC~253 was observed by XMM-Newton as a performance verification target in June 2000. First results on bright point sources were discussed in Pietsch et al. (2001, hereafter P01). Of special interest was the spectrum and time variability of the ROSAT source X33 30″ south of the nucleus, most likely a black hole XRB. The pair of ROSAT sources (X17/X18) is clearly visible in Fig. 1 of P01, however deferred for later analysis, as these sources at the time could not be distinguished by the SAS source detection algorithms.

NGC~253 was again observed by us with XMM-Newton in December 2000 within the XMM-Newton survey scientist guaranteed time program. We noticed that during this observation the count rate of RX~J004717.4-251811 abruptly changed from very low to high. This behavior most likely represents an egress from an eclipse by the companion star in this XRB. Triggered by this finding we analyzed archival Chandra, XMM-Newton, ROSAT, and Einstein observations. This led to the detection of a similar behavior during Chandra observation 790. Data from all the satellites were then used to further constrain the binary orbit parameters and long term time variability of the source. In the following we report on these findings.

## 2 Observations

For the detailed analysis of RX~J004717.4-251811 we used data from XMM-Newton, Chandra, ROSAT HRI, and Einstein. Table 1 summarizes the observatory (col. 1), observation identification (2), observation dates (3) and life time (4), extraction radius Re used for XMM-Newton and Chandra count rates and light curves and spectra (5), RX~J004717.4-251811 count rates (6), hardness ratios for XMM-Newton EPIC PN and Chandra ACIS S (7), and luminosities in the 0.5–2.4 keV band (8), an energy band covered by all observatories. Special events are marked under comments (9). As hardness ratio (HR) we use the ratio of the counts in the 1.5–5.0 keV band to the counts in the 0.5–1.5 keV band.

During the XMM-Newton observations (Jansen et al. 2001) the EPIC PN and MOS instruments (Strüder et al. 2001; Turner et al. 2001) were operated in the full frame mode resulting in a time resolution of 73.4 ms and 2.6 s, respectively. The medium filter was in front of the EPIC MOS1 camera in all observations, for MOS2, the thin filter was used during the observations in June 2000 and the medium filter in December. The PN camera operated with the medium filter during observation 0125960101 and with the thin filter for the other observations. These filter changes do not influence the RX~J004717.4-251811 observations significantly as source counts are mainly detected above 0.7 keV (see below). In this energy band the transmission of the thin and medium filter do not differ significantly. We used all EPIC instruments for the imaging and position determination of RX~J004717.4-251811. For the timing and spectral investigations we concentrated on the EPIC PN camera which gives about twice the number of photons and checked the results for consistency and at the beginning of the observing window using the MOS cameras. While the XMM-Newton point spread function (PSF) normally requires extraction radii $`>30`$″ to encircle $`>80\%`$ of the source photons, we had to restrict ourselves to 10″ to avoid contamination from source VP99 X18 (see Fig. 3).

Three Chandra ACIS S observations of NGC~253 (see Table 1) were obtained from the Chandra Data Archive (http://asc.harvard.edu/cgi-gen/cda). The instrument was operated in the full frame mode (3.2 s time resolution). RX~J004717.4-251811 is positioned in the back-illuminated CCD chip S3 during observation 969 and 383, in the front-illuminated chip S2 during 790. Due to the larger off-axis position of RX~J004717.4-251811 in the latter observation, the Chandra PSF urged us to use a larger extraction radius.

The deep space orbits of the satellites XMM-Newton and Chandra led to long observation times of RX~J004717.4-251811 without interruption. The low earth orbits of the Einstein and ROSAT satellite on the other hand, led to observations split in many short intervals of typically less than 1 500 s.

NGC~253 was observed with the ROSAT PSPC and HRI in several observing runs from Dec 1991 to Dec 1998. We restricted our analysis on 75.5 ks HRI observations as only this ROSAT detector fully resolves RX~J004717.4-251811 from VP99 X18. The HRI observation 601113h (July 1998, 2 ks) could not be used for the investigations due to attitude problems.

RX~J004717.4-251811 can not be resolved from nearby sources with the Einstein IPC. However, the source was clearly detected with the HRI (FT84) during a bright state. The nearby source VP99 X18 was not detected. For additional timing analysis we extracted the NGC~253 HRI data from the Einstein Data Archive at HEASARC.

The data analysis was performed using tools in the SAS v5.3.3, CIAO v2.2, EXSAS/MIDAS 1.2/1.4, and FTOOLS v5.1 software packages, the imaging application DS9 v2.1, the timing analysis package XRONOS v5.19 and spectral analysis software XSPEC v11.2.

For the time variability investigations all RX~J004717.4-251811 event times were corrected to solar system barycenter arrival times.

### 2.1 Time variability

The EPIC PN as well as the MOS light curves for XMM-Newton observation 0110900101 show egresses from low to high state at solar system barycenter corrected MJD $`51892.146\pm 0.006`$ ( for EPIC PN see Fig. 1) with residual emission during low state which partly may be contaminated from emission collected from the PSF wing of VP99 X18. During the high state the average HR was 0.86 and constant within the errors (average HR 0.86, see Table 1), the HR during the low state (0.5) indicates a softer spectrum and/or less absorption. During each of the XMM-Newton observations 0125960101 and 0125960201 the flux was constant within the errors (see Table 1) and the source was a factor of 1.9 brighter than in the 0110900101 low state, a factor of more than five below that in the 0110900101 high state. The HRs during these observations on the other hand, appear to be quite similar to the 0110900101 low state.

The light curve for the Chandra observation 790 shows an egress from low to high state at solar system barycenter corrected MJD $`51539.276\pm 0.006`$ (Fig. 2). During the high state the HR defined in the same energy bands as for XMM-Newton EPIC, was 0.91 and again constant within the errors. The residual emission during the low state clearly originates from RX~J004717.4-251811 as is demonstrated in the right part of Fig. 3 and corresponds to a count rate of 1.4 ct ks-1. The HR again is significantly smaller than during the high state. The low state flux was about a factor of 2 below that in the XMM-Newton low state. During Chandra observation 969, RX~J004717.4-251811 flux and HR were constant within the errors and the source about a factor of 1.4 brighter with similar average HR than that measured in the high state of 790. During Chandra observation 383 the RX~J004717.4-251811 was a factor of 12 below that during 969, but a factor of 2.3 higher than in the 790 low state, and showed a small HR (see Table 1). However, one has to keep in mind that the Chandra ACIS S count rates and HRs of RX~J004717.4-251811 can not directly be compared to the others as they originate from a CCD with a different energy response function. The same is even more true if one wants to compare Chandra and XMM-Newton count rates and HRs.

Figure 3 shows images (1′ to a side) of the RX~J004717.4-251811 area in the 0.5–5.0 keV band. For XMM-Newton all EPIC cameras are combined for the high state of observation 0110900101. RX~J004717.4-251811 is the source to the north, VP99 X18 the one to the south. The Chandra ACIS S images of observation 790 are divided in a low and high state image. RX~J004717.4-251811 is clearly detected also during the low state.

We searched for pulsations within the bright parts of the XMM-Newton and Chandra observations in the frequency range 10-3–3 Hz. The strongest signal was found at 29.5 s in the XMM-Newton observation 0110900101 which, following the Rayleigh Z$`{}_{\mathrm{n}}{}^{}{}_{}{}^{2}`$ method (Buccheri et al. 1983) as described in Haberl & Zavlin (2002) yields only a 71% confidence level ($`1\sigma `$).

During the ROSAT HRI observations covering 6 years, 110 photons were detected within 6″ from RX~J004717.4-251811 in an effective observing time of 75.5 ks. Due to the low ROSAT HRI background less than 10 of these counts may be due to field background. The source photons were not distributed evenly over the observing time but clearly showed strong intensity changes and even short term flares as is demonstrated by a Kolmogorov-Smirnov test against constancy (see Fig. 4).

### 2.2 Improved position

VP99 give the position of RX~J004717.4-251811 with an error of 2$`\stackrel{}{.}`$7 which is mostly determined by the included systematic error of 2$`\stackrel{}{.}`$5. The much bigger number of photons detected with the XMM-Newton EPIC and Chandra ACIS S detectors allows us to determine a significantly improved source position. Systematic errors can be reduced by adjusting the positions using the well determined optical (0$`\stackrel{}{.}`$2 error using USNO-B positions, Monet et al. 2003) and X-ray positions of several QSOs in the field of view.

For Chandra ACIS S observation 790 the VP99 source X22 identified with a QSO (optical blue magnitude of 19.9, redshift 1.25) is detected in chip S3 close to the center of the field of view. During observation 969 it is positioned at a larger off-axis angle in chip S4. For XMM-Newton EPIC besides VP99 X22 the sources X4 and X58 are clearly detected (QSOs with optical blue magnitudes of 20.3 and 18.5 and redshifts of 1.02 and 0.66, respectively). Averaging the resulting positions from Chandra observations 790 and 969 as well as from the individual EPIC instruments in observation 0110900101 we get a highly improved position for RX~J004717.4-251811: $`\alpha `$=00h47m17$`\stackrel{s}{.}`$51, $`\delta `$=$``$25°18′11$`\stackrel{}{.}`$2 (J2000), with an remaining error radius of 0$`\stackrel{}{.}`$3 mainly determined by systematics. See Fig. 5 for an overlay of the improved position on a deep optical blue image.

### 2.3 Energy spectra

As discussed in Sect. 2.1, HR changes of RX~J004717.4-251811 in different brightness regimes indicate spectral variability which may be either caused by reduced absorption during the low intensity states or due to changes of the stiffness of the source spectrum. Unfortunately, the source count rates do not allow a more detailed investigation. For the brighter states of RX~J004717.4-251811, that is during Chandra observation 969 as well as during the high state of Chandra observation 790, and XMM-Newton observation 0110900101 more than 600 counts were collected from the source which allowed us a rough spectral analysis. RX~J004717.4-251811 photons were extracted using the same areas as for the time variability analysis and spectra were grouped in channels with at least 30 source plus background photons. Due to the limited number of photons we just fitted absorbed one component spectral models: power law (POWL), bremsstrahlung (BREMS), and thin thermal plasma with element abundance fixed to solar (Raymond-Smith, THPL). The resulting integration times, raw count rates and degrees of freedom (DOF) are listed in Table 2 together with the results of the fits.

The absorbed POWL model fits show the lowest $`\chi ^2/\nu `$ values. Absorbing column and photon index $`\mathrm{\Gamma }`$ of all observations coincide within the errors. However, only for the Chandra observations the fit is formally acceptable. Several factors may explain the unacceptably high $`\chi ^2/\nu `$ of the POWL fit for the XMM-Newton EPIC PN spectrum of RX~J004717.4-251811. To the one hand, there may be contributions to the spectrum from the nearby source X18. On the other hand, to get sufficient statistics for the fit we accepted all photons in the extraction region, not rejecting events on or close to a CCD column with slightly increased offset which passes through the extraction region. The BREMS spectra exhibit marginally higher $`\chi ^2/\nu `$ values than the POWL fits and find high temperatures (around 7 to 10 keV). The absorbing column needed in addition to the POWL and BREMS models is $`2\times 10^{21}`$ cm-2, significantly above the absorption expected within the Milky Way in this direction ($`1.3\times 10^{20}`$ cm-2, Dickey & Lockman 1990). Neither the Chandra nor the XMM-Newton data can be described by a THPL spectrum.

The best fit POWL model parameters of Chandra observation 969 are used in Table 1 for conversion from count rates to intrinsic source luminosities for observations of all observatories. We are aware that the smaller HR values during low intensity states indicate a higher low energy flux and therefore a different spectrum. This may lead to an overestimation of the intrinsic low state luminosities in the table.

## 3 Orbital period determination

The intensity changes of RX~J004717.4-251811 from low to high in the XMM-Newton observation 0110900101 and in the Chandra observation 790 are not resolved by the 1000 s time resolution of our plots. Shorter time bins do not help as the source is not bright enough. In the following we interpret the intensity changes as egresses from eclipse of an XRB. The times of egresses have to be separated by an integer number of periods and therefore give a first constrain on the possible orbital period p of RX~J004717.4-251811:

p = ($`352.870\pm 0.012`$) d / n ; n = 1,2,3,… (1)

without exactly knowing the number n of periods between the two egresses. The relative period error is 3.3$`\times 10^5`$.

The Chandra observation 790 is the longest continuous observation of RX~J004717.4-251811 covering 0.5092 d, however only part of an orbit. This together with Eq. 1 constrains the number of possible periods to 692 (covering the range 0.5099 d and 352.870 d).

During Chandra observation 969 the flux of RX~J004717.4-251811 was constant within the errors and even brighter than during the 790 bright state. The source therefore undoubtedly was out of eclipse. If we use the conditions that the non-eclipse observations are not allowed to include orbit phases covering the minimum eclipse duration (0.14 d during the Chandra observation) and assuming that the eclipse phase is longer than 0.13. This constrains the number of possible periods to 457 (between 0.5129 d and 70.574 d).

We can try to further reduce the number of acceptable periods with the help of the other XMM-Newton and Chandra observations and in addition use the sparse ROSAT and Einstein HRI photons. This selection however has to be taken with care as it is hampered by the limited knowledge of the luminosity state of RX~J004717.4-251811 during these observations and the limited photon statistics.

XMM-Newton observations 0125960101 (duration 0.451 d) and 0125960201 (0.164 d) were separated by 0.174 d. RX~J004717.4-251811 was again constant within the errors during each observation and also from observation to observation, however only a factor of 1.9 brighter than during eclipse of 0110900101 with comparable hardness ratio values. As the residual flux during eclipse may vary from eclipse to eclipse similar to the source flux out of eclipse we can not decide if the observations are out of eclipse (I) or if one eclipse covers both observations (II). The condition (I) reduces the number of allowed periods to 141 (between 0.6416 d and 70.574 d). If we assume for condition (II) a maximum eclipse duration of 0.25 in phase (similar to that reported for Cen X-3), the period has to be longer than 3.156 d, the two observations at the right phase and Chandra 969 not within the same phase interval. The condition (II) reduces the number of allowed periods to 10 (between 4.1031 d and 32.079 d).

Also for Chandra observations 383 we can not decide if the source was in or out of eclipse as the source is only a factor of two brighter than during eclipse of 790 and shows similar hardness ratio. For 23 of the 141 period candidates of condition (I) and 3 of the 10 for (II) the source would be in eclipse during Chandra observation 383.

The detection of two eclipse egresses within the few XMM-Newton and Chandra observations of less than 0.5 d duration strongly favors a short orbital period of RX~J004717.4-251811. If we constrain the period to less than 10 d there only remain 118 periods for condition (I) and 5 for (II) to be investigated. During 18 (1) of those RX~J004717.4-251811 would be in eclipse during Chandra observation 383.

To further reduce the number of acceptable periods we folded the ROSAT HRI counts detected from the RX~J004717.4-251811 area over the period ranges allowed by the candidate periods. We checked if the source could have remained in low intensity (eclipse) for at least the duration of the measured longest low states (0.14 d) or for 0.13 of the candidate period and if this low intensity time was compatible in phase with the Chandra eclipse egress for period candidates derived with condition (I). For condition (II) the eclipse had to last as long as indicated by the XMM-Newton observations. Many of the candidate periods could be rejected using these criteria. Some periods can not be rejected because the expected eclipse phase is not covered by ROSAT exposure. Some periods however show only background counts at the expected eclipse phases even when the ROSAT exposure is high. All periods selected under condition (II) are not acceptable using these criteria. Under condition (I), we get 14 candidate periods with a lower error range. During three of the periods, RX~J004717.4-251811 is in eclipse during Chandra observation 383.

In a last step, the candidate periods were checked for consistency with the Einstein observations (collected during a RX~J004717.4-251811 high state in July 1979, see Table 1). Only seven periods, two with Chandra observation 383 in eclipse, pass this test (see Table 3). The periods above 5 days are less probable as the chance to detect two eclipse egresses during the XMM-Newton and Chandra observations would have been low. From the remaining periods the ones with allowed eclipse duration longer than 0.13 in phase seem to be more probable.

Taking all the selection criteria into account 1.470243 d seems to be the best period candidate followed by 3.207928 d. Additional observations of RX~J004717.4-251811 in the high state are urgently needed to finally decide on the correct orbital period.

## 4 Discussion

In the last sections we have shown that the variety of observations of RX~J004717.4-251811 can be explained in an eclipsing XRB scenario. RX~J004717.4-251811 is the first eclipsing XRB identified outside the Local Group.

The length of eclipse and the luminosity during its high state of $`2\times 10^{38}`$ erg s-1, i.e. at the Eddington limit of a 1.4 M neutron star, are most easily explained by a HMXB with a neutron star as companion. Similar luminosities were reported for the eclipsing HMXBs SMC X-1 and LMC X-4 during the high state. Detection of pulsations in the tens of second range or shorter are expected in such systems. For a faint source like RX~J004717.4-251811 such short pulsations are difficult to detect as time delays of photon arrival times due to the motion in the binary orbit will smear out the periodic signal.

Accreting XRBs normally show pulsation averaged X-ray energy spectra above 1 keV that can be modeled by power laws with high energy cut-offs. The photon indices of the power law spectra range between 0.8 and 1.5 for most pulsars, but some sources show rather soft spectra (see the review of Nagase 1989). The steep power law slope of 1.9 in the brighter states puts RX~J004717.4-251811 in this special category.

For most eclipsing XRB residual emission during eclipse was measured. The emission can be explained by re-processing of primary X-rays originating from the bright area at the magnetic poles of the neutron star, in an extended accretion disk corona (which is not fully occulted) or by scattering in the companion atmosphere/stellar wind. Residual emission of up to $``$10% was reported (Ebisawa et al. 1996; Haberl 1991) depending on system geometry and wind density. The RX~J004717.4-251811 residual emission of $``$5% is well within these limits.

As discussed in Sect. 3, RX~J004717.4-251811 was most likely out of eclipse during the XMM-Newton observations in December 2000 and the Chandra observation in August 2001. During these observations the source luminosity was up to a factor of ten lower than during the high states (luminosities close to Eddington) indicating clear long term time variability. During the Einstein observations in July 1979 and also during the ROSAT observations in January 1995 the source was in high state, while during the rest of the ROSAT observations the intensity was significantly lower (see Fig. 4). Similar long term variability was detected in several other low and high mass XRBs where it even may manifest itself in a long term superorbital period (e.g. Her X-1, LMC X-4, SMC X-1) typically on the scale of tens of days. Due to the up to now sparse time coverage of RX~J004717.4-251811 observations, such long term periodicities were not detected, but at present also can not be ruled out. The Her X-1 and LMC X-4 long term periodicities correspond to $`21`$ orbital periods. If a RX~J004717.4-251811 long term variability follows a similar pattern, the most probable orbital period of 1.47 d would lead to $`31`$ d as period of the long term variability.

During some ROSAT observation intervals of typically 1000 s duration there are some indications of short term flaring (see Fig. 4). Similar flares up to now were not detected during XMM-Newton and Chandra observations. Short term flaring was detected in several HMXBs (e.g. LMC X-4, Pietsch et al. 1985).

The most promising orbital periods of RX~J004717.4-251811 do not allow us to decide if the optical companion is a high or low mass star. LMC X-4 and SMC X-1 are HMXB systems with orbital periods of 1.4 d and 3.9 d, Her X-1 with an intermediate period (1.7 d) on the other hand is a LMXB. The X-ray luminosity of RX~J004717.4-251811 in the bright state at or above the Eddington limit and possible flaring (see above) may indicate a HMXB as the LMXB Her X-1 only reaches about one third of the Eddington luminosity (Forman et al. 1972; Howarth & Wilson 1983) and does not show flares.

If we therefore assume that RX~J004717.4-251811 is a HMXB we can use the optical magnitudes of the known systems in the Magellanic Clouds to determine the optical brightness in NGC~253. The optical brightnesses of LMC X-4 and SMC X-1 – V of 14.0 mag and 13.3 mag (Liu et al. 2000) and distances of 50 kpc and 59 kpc, respectively (see review by van den Bergh 1999) – correspond to (21.4–22.6) mag at the assumed distance of 2.58 Mpc to NGC~253. If this distance has to be corrected to 3.94 MPc as recently indicated by Karachentsev et al. (2003), not only the estimated optical magnitudes would have to be increased by 0.9 mag but also the derived X-ray luminosities for RX~J004717.4-251811 would rise by a factor of 2.3. The expected optical magnitude may have to be increased by up to one optical magnitude if the absorption measured in the X-ray spectrum in the high state ($`N_\mathrm{H}`$=1.9$`\times 10^{21}`$ cm-2, see above) is due to extinction within NGC~253 and not due to absorption in the near surrounding of the X-ray emitting area within the binary system (Predehl & Schmitt 1995). Additional extinction within NGC~253 may well be present as RX~J004717.4-251811 seems to be positioned at the edge of a dark cloud (see Fig. 5). A preliminary analysis of the image indicates that stars of visual magnitude 23 should just be visible. There is a faint object at the east edge of the error box which could be the optical counterpart. However, if RX~J004717.4-251811 is a LMXB similar to Her X-1 (14 mag at a distance of 6 kpc, see e.g. Howarth & Wilson 1983), it would be fainter than 27 mag at the distance of NGC~253, very difficult to be separated from the many other stars of similar brightness in the NGC~253 surrounding, and not be detectable in current optical imaging. Only deeper time resolved imaging and spectroscopy can lead to a secure identification of RX~J004717.4-251811as a HMXB and solve the high mass – low mass ambiguity.

## 5 Conclusions and outlook

Starting from an XMM-Newton observation in which RX~J004717.4-251811 changed state from low to high, we found a similar transition in a Chandra observation. Interpreted as eclipse egresses we determined candidate orbital periods for the source. With the help of archival XMM-Newton, Chandra, ROSAT and Einstein observations we selected as most probable orbital periods 1.47029 d and 3.20793 d. RX~J004717.4-251811 is one of the few eclipsing X-ray binary systems and the first detected outside the Local Group. No pulsation period was found. The source is clearly variable on time scales of tens of days and also short time flaring (1000 s time scale) is indicated. The energy spectrum of the source during the high state can be described by a hard power law spectrum with absorption significantly higher than the foreground absorption in the Milky Way. Using the improved source position we indicate a possible optical identification. RX~J004717.4-251811 is most likely a HMXB system similar to SMC X-1 and LMC X-4.

Two long observations of NGC~253 accepted for the XMM-Newton and Chandra program (140 ks and 85 ks, respectively) hopefully will detect RX~J004717.4-251811 in the high state and allow us to finally decide on the orbital period of the source and enable sensitive neutron star rotation searches from seconds to tens of seconds. The XMM-Newton observation would cover one full binary orbit if the 1.5 d period is correct and finally clarify the correct orbital solution. A search for sub-second periods will only be possible in dedicated observations of RX~J004717.4-251811 with higher time resolution (that could not fully cover the galaxy).

###### Acknowledgements.

The XMM-Newton project is supported by the Bundesministerium für Bildung und Forschung / Deutsches Zentrum für Luft- und Raumfahrt (BMBF/DLR), the Max-Planck Society and the Heidenhain-Stiftung.


format: text

text:
# Noiseless scattering states in a chaotic cavity

(11 February 2003)

## Abstract

Shot noise in a chaotic cavity (Lyapunov exponent $`\lambda `$, level spacing $`\delta `$, linear dimension $`L`$), coupled by two $`N`$-mode point contacts to electron reservoirs, is studied as a measure of the crossover from stochastic quantum transport to deterministic classical transport. The transition proceeds through the formation of fully transmitted or reflected scattering states, which we construct explicitly. The fully transmitted states contribute to the mean current $`\overline{I}`$, but not to the shot-noise power $`S`$. We find that these noiseless transmission channels do not exist for $`N\sqrt{k_FL}`$, where we expect the random-matrix result $`S/2e\overline{I}=1/4`$. For $`N\sqrt{k_FL}`$ we predict a suppression of the noise $`(k_FL/N^2)^{N\delta /\pi \mathrm{}\lambda }`$. This nonlinear contact dependence of the noise could help to distinguish ballistic chaotic scattering from random impurity scattering in quantum transport.

Shot noise can distinguish deterministic scattering, characteristic for particles, from stochastic scattering, characteristic for waves. Particle dynamics is deterministic: A given initial position and momentum fixes the entire trajectory. In particular, it fixes whether the particle will be transmitted or reflected, so the scattering is noiseless. Wave dynamics is stochastic: The quantum uncertainty in position and momentum introduces a probabilistic element into the dynamics, so it is noisy.

The suppression of shot noise in a conductor with deterministic scattering was predicted many years ago from this qualitative argument Bee91 . A better understanding, and a quantitative description, of how shot noise measures the transition from particle to wave dynamics in a chaotic quantum dot was put forward by Agam, Aleiner, and Larkin Aga00 , and developed further in Ref. Hen02 . The key concept is the Ehrenfest time $`\tau _E`$, which is the characteristic time scale of quantum chaos Zas81 . The noise power $`S\mathrm{exp}(\tau _E/\tau _D)`$ was predicted to vanish exponentially with the ratio of $`\tau _E`$ and the mean dwell time $`\tau _D=\pi \mathrm{}/N\delta `$ in the quantum dot (with $`\delta `$ the level spacing and $`N`$ the number of modes in each of the two point contacts through which the current is passed). A recent measurement of the $`N`$-dependence of $`S`$ is consistent with this prediction for $`\tau _E<\tau _D`$, although an alternative explanation in terms of short-range impurity scattering describes the data equally well Obe02 .

The theory of Ref. Aga00  introduces the stochastic element by means of long-range impurity scattering and adjusts the scattering rate so as to mimic the effect of a finite Ehrenfest time. Here we take the alternative approach of explicitly constructing noiseless channels in a chaotic quantum dot. These are scattering states which are either fully transmitted or fully reflected in the semiclassical limit. They are not described by random matrix theory Jal94 . By determining what fraction of the available channels is noiseless, we can deduce a precise upper bound for the shot-noise power. A random matrix conjecture for the remaining noisy channels gives an explicit form of $`S(N)`$. We find that the onset of the classical suppression of the noise is described not only by the Ehrenfest time, but by the difference of $`\tau _E`$ and the ergodic time $`\tau _0`$, which we introduce and calculate in the paper. The resulting nonlinear dependence of $`\mathrm{ln}S`$ on $`N`$ may help to distinguish between the competing explanations of the experimental data Obe02 .

We illustrate the construction of noiseless scattering states for the two-dimensional billiard with smooth confining potential $`U(x,y)`$ shown in Fig. 1. The outer equipotential defines the area in the $`xy`$ plane which is classically accessible at the Fermi energy $`E_F=p_F^2/2m`$ (with $`p_F=\mathrm{}k_F`$ the Fermi momentum). The motion in the closed billiard is chaotic with Lyapunov exponent $`\lambda `$. We assume the billiard to be connected at $`x=0`$ and $`x=L`$ by two similar point contacts to leads of width $`W`$ extended along the $`\pm x`$-direction.

The beam of electrons injected through a point contact into the billiard has cross-section $`W`$ and transverse momenta in the range $`(p_W,p_W)`$. The number of channels $`Np_WW/\mathrm{}`$ in the lead is much smaller than the number of channels $`Mp_FL/\mathrm{}`$ supported by a typical cross-section of the billiard. While $`W/L1`$ in general, the ratio $`p_W/p_F`$ depends on details of the potential near the point contact. If $`p_W/p_F1`$ one speaks of a collimated beam. This is typical for a smooth potential, while a hard-wall potential typically has $`p_Wp_F`$ (no collimation). We define $`r_{\mathrm{min}}=\mathrm{min}(W/L,p_W/p_F)`$ and $`r_{\mathrm{max}}=\mathrm{max}(W/L,p_W/p_F)`$.

The classical phase space is four-dimensional. By restricting the energy to $`E_F`$ and taking $`x=0`$ we obtain the two-dimensional section of phase space shown in Fig. 2. The accessible values of $`y`$ and $`p_y`$ lie in a disc-shaped region of area $`𝒜=Nh`$ in this surface of section. Up to factors of order unity, the disk has width $`r_{\mathrm{min}}`$ and length $`r_{\mathrm{max}}`$ (if coordinate and momentum are measured in units of $`L`$ and $`p_F`$, respectively). In Fig. 2 one has $`r_{\mathrm{min}}r_{\mathrm{max}}`$. Each point in the disc defines a classical trajectory that enters the billiard (for positive $`p_x`$) and then leaves the billiard either through the same lead (reflection) or through the other lead (transmission). The points lie in narrow bands, which we will refer to as “transmission bands” and “reflection bands”.

It is evident from Fig. 2 that the area $`A_j`$ enclosed by a typical transmission (or reflection) band $`j`$ is much less than $`𝒜`$. For an estimate we consider the time $`t(y,p_y)`$ that elapses before transmission. Let $`t_j`$ be the dwell time averaged over the starting points $`y,p_y`$ in a single band. The fluctuations of $`t`$ around the average are of the order of the time $`t_WmW/p_W`$ to cross the point contact, which is typically $`t_j`$. As we will see below, the area of the band decreases with $`t_j`$ as

$$A_j𝒜_0\mathrm{exp}(\lambda t_j)\mathrm{if}t_j1/\lambda ,t_W.$$
(1)
The prefactor $`𝒜_0=𝒜r_{\mathrm{min}}/r_{\mathrm{max}}`$ depends on the degree of collimation. In Ref. Sil02  the symmetric case $`r_{\mathrm{min}}=r_{\mathrm{max}}`$ was assumed, when $`𝒜_0=𝒜`$.

We now proceed to the construction of fully transmitted scattering states. To this end we consider a closed contour $`𝒞`$ within a transmission band $`j`$. The starting points on the contour define a family of trajectories that form a flux tube inside the billiard (see Fig. 1). The semiclassical wave function

$$\psi (x,y)=\underset{\sigma }{}\sqrt{\rho _\sigma (x,y)}\mathrm{exp}[i𝒮_\sigma (x,y)/\mathrm{}]$$
(2)
is determined as usual from the action $`𝒮_\sigma `$ and density $`\rho _\sigma `$ that solve the Hamilton-Jacobi and continuity equations

$$|𝒮|^2=2m(E_FU),(\rho 𝒮)=0.$$
(3)
The action is multivalued and the index $`\sigma `$ labels the different sheets. Typically, there are two sheets, one originating from the upper half of the contour $`𝒞`$ and one from the lower half.

The requirement that $`\psi `$ is single-valued as one winds around the contour imposes a quantization condition on the enclosed area,

$$_𝒞p_y𝑑y=(n+1/2)h.$$
(4)
The increment $`1/2`$ accounts for the phase shift acquired at the two turning points on the contour. The quantum number $`n=0,1,2,\mathrm{}`$ is the channel index. The largest value of $`n`$ occurs for a contour enclosing an area $`A_j`$. The number of transmission channels $`N_j`$ within band $`j`$ is therefore given by $`A_j/h`$, with an accuracy of order unity. In view of Eq. (1) we have

$`N_j(𝒜_0/h)\mathrm{exp}(\lambda t_j),\mathrm{for}t_j<\tau _E,`$ (5a)
$`N_j=0,\mathrm{for}t_j>\tau _E.`$ (5b)
The time

$$\tau _E=\lambda ^1\mathrm{ln}(𝒜_0/h)=\lambda ^1\mathrm{ln}(Nr_{\mathrm{min}}/r_{\mathrm{max}})$$
(6)
above which there are no fully transmitted channels is the Ehrenfest time of this problem.

By decomposing one of these $`N_j`$ scattering states into a given basis of transverse modes in the lead one constructs an eigenvector of the transmission matrix product $`tt^{}`$. The corresponding eigenvalue $`𝒯_{j,n}`$ equals unity with exponential accuracy in the semiclassical limit $`n1`$. Because of the degeneracy of this eigenvalue any linear combination of eigenvectors is again an eigenvector. This manifests itself in our construction as an arbitrariness in the choice of $`𝒞`$.

We observe in Fig. 1 that the spatial density profile $`\rho (x,y)`$ of a fully transmitted scattering state is highly non-uniform. The flux tube is broad (width of order $`W`$) at the two openings, but is squeezed down to very small width inside the billiard. A similar effect was noted Sil02  in the excited states of an Andreev billiard (a cavity connected to a superconductor). Following the same argument we estimate the minimal width of the flux tube as $`W_{\mathrm{min}}L\sqrt{N_j/k_FL}`$.

The total number

$$N_0=\underset{j}{}N_j=N_0^{\tau _E}P(t)𝑑t$$
(7)
of fully transmitted and reflected channels is determined by the dwell-time distribution $`P(t)`$ Bau90 . Fig. 3 shows this distribution in our billiard. One sees three different time scales. The narrow peaks represent individual transmission (reflection) bands. They consist of an abrupt jump followed by an exponential decay with time constant $`t_W`$. These exponential tails correspond to the borders of the bands, where the trajectory bounces many times between the sides of the point contact. If we smooth $`P(t)`$ over such short time intervals, an exponential decay with time constant $`\tau _D=\pi \mathrm{}/N\delta `$ is obtained (inset). The decay starts at the so called “ergodic time” $`\tau _0`$. There are no trajectories leaving the cavity for $`t<\tau _0`$. So the smoothed dwell-time distribution has the form.

$$P(t)=\tau _D^1\mathrm{exp}[(\tau _0t)/\tau _D]\theta (t\tau _0),$$
(8)
with $`\theta (t)`$ the unit step function.

In order to find $`\tau _0`$ we consider Fig. 4, where the section of phase space along a cut through the middle of the billiard is shown (line $`b`$ in Fig. 1). It is convenient to measure the momentum and coordinate along $`b`$ in units of $`p_F`$ and $`L`$. The injected beam crosses the section for the first time over an area $`O_{\mathrm{initial}}`$ of size $`r_{\mathrm{max}}\times r_{\mathrm{min}}=hN/p_FL`$. (Fig. 4 has $`r_{\mathrm{min}}r_{\mathrm{max}}`$, but the estimates hold for any $`r_{\mathrm{min}}<r_{\mathrm{max}}<1`$.) Further crossings consist of increasingly more elongated areas. The fifth crossing is shown in Fig. 4. The flux tube intersects line $`b`$ in a few disjunct areas $`O_j`$, of width $`r_{\mathrm{min}}e^{\lambda t}`$ and total length $`r_{\mathrm{max}}e^{\lambda t}`$. (Due to conservation of the integral $`𝐩𝑑𝐫`$ enclosing the flux tube, the total area $`_jO_j`$ decreases only when particles leave the billiard.) The typical separation of adjacent areas is $`(r_{\mathrm{max}}e^{\lambda t})^1`$. To leave the billiard (through the right contact) without further crossing of $`b`$ a particle should pass through an area $`O_{\mathrm{final}}r_{\mathrm{max}}\times r_{\mathrm{min}}`$. This is highly improbable endend  until the separation of the areas $`O_j`$ becomes of order $`r_{\mathrm{max}}`$, leading to the ergodic time

$$\tau _0=\lambda ^1\mathrm{ln}r_{\mathrm{max}}^2.$$
(9)
The ergodic time varies from $`\tau _0\lambda ^1`$ for $`r_{\mathrm{max}}1`$ to $`\tau _0=\lambda ^1\mathrm{ln}(k_FL/N)`$ for $`r_{\mathrm{min}}r_{\mathrm{max}}`$. The overlap of the areas $`O_j`$ and $`O_{\mathrm{final}}`$ is the mapping of the transmission band onto the surface of section $`b`$. It has an area $`p_FLr_{\mathrm{min}}^2e^{\lambda t}=𝒜(r_{\mathrm{min}}/r_{\mathrm{max}})e^{\lambda t}`$, leading to Eq. (1).

Substituting Eq. (8) into Eq. (7) we arrive at the number $`N_0`$ of fully transmitted and reflected channels,

$`N_0=N\theta (\tau _E\tau _0)\left[1e^{(\tau _0\tau _E)/\tau _\mathrm{D}}\right],`$ (10)
$`\tau _E\tau _0=\lambda ^1\mathrm{ln}(N^2/k_FL).`$ (11)
There are no fully transmitted or reflected channels if $`\tau _E<\tau _0`$, hence if $`N<\sqrt{k_FL}`$. Notice that the dependence of $`\tau _E`$ and $`\tau _0`$ separately on the degree of collimation drops out of the difference $`\tau _E\tau _0`$. The number of noiseless channels is therefore insensitive to details of the confining potential. An Ehrenfest time $`\mathrm{ln}(N^2/k_FL)`$ has appeared before in connection with the Andreev billiard Vav02 , but the role of collimation (and the associated finite ergodic time) was not considered there.

Eqs. (5) and (8) imply that the majority of noiseless channels group in bands having $`N_j1`$, which justifies the semiclassical approximation. The total number of these noiseless bands is $`(NN_0)/\lambda \tau _D`$, which is much less than both $`NN_0`$ and $`N_0`$. Because of this inequality the relatively short trajectories contributing to the noiseless channels are well separated in phase space from other, longer trajectories (cf. Fig. 2).

The shot-noise power $`S`$ is related to the transmission eigenvalues by But90

$$S=2e\overline{I}g^1\underset{k=1}{\overset{N}{}}𝒯_k(1𝒯_k),$$
(12)
with $`\overline{I}`$ the time-averaged current and $`g=_k𝒯_k`$ the dimensionless conductance. The $`N_0`$ fully transmitted or reflected channels have $`𝒯_k=1`$ or $`0`$, hence they do not contribute to the noise. The remaining $`NN_0`$ channels contribute at most 1/4 per channel to $`Sg/2e\overline{I}`$. Using that $`g=N/2`$ for large $`N`$, we arrive at an upper bound for the noise power $`S<e\overline{I}(1N_0/N)`$.

For a more quantitative description of the noise power we need to know the distribution $`P(𝒯)`$ of the transmission eigenvalues for the $`NN_0`$ noisy channels, which can not be described semiclassically. We expect the distribution to have the same bimodal form $`P(𝒯)=\pi ^1𝒯^{1/2}(1𝒯)^{1/2}`$ as in the case $`N_0=0`$ Jal94 . This expectation is motivated by the earlier observation that the $`N_0`$ noiseless channels are well separated in phase space from the $`NN_0`$ noisy ones. Using this form of $`P(𝒯)`$ we find that the contribution to $`Sg/2e\overline{I}`$ per noisy channel equals $`_0^1𝒯(1𝒯)P(𝒯)𝑑𝒯=1/8`$, half the maximum value. The Fano factor $`F=S/2e\overline{I}`$ is thus estimated as

$`F`$ $`=`$ $`\frac{1}{4},\mathrm{for}N\sqrt{k_FL},`$ (13a)
$`F`$ $`=`$ $`\frac{1}{4}(k_FL/N^2)^{N\delta /\pi \mathrm{}\lambda },\mathrm{for}N\sqrt{k_FL}.`$ (13b)
This result should be compared with that of Ref. Aga00 : $`F^{}=\frac{1}{4}(k_FL)^{N\delta /\pi \mathrm{}\lambda }`$. The ratio $`F^{}/F=\mathrm{exp}[(2N\delta /\pi \mathrm{}\lambda )\mathrm{ln}(N/k_FL)]`$ is always close to unity (because $`N\delta /\pi \mathrm{}\lambda N/k_FL1`$). But $`F\frac{1}{4}`$ and $`F^{}\frac{1}{4}`$ are entirely different for $`N\sqrt{k_FL}`$, which is the relevant regime in the experiment Obe02 . There the $`N`$-dependence of the shot noise was fitted as $`F=\frac{1}{4}(1t_Q/\tau _D)=\frac{1}{4}(1\mathrm{constant}\times N)`$, where $`t_Q`$ is some $`N`$-independent time. Eq. (13) predicts a more complex $`N`$-dependence, a plateau followed by a decrease as $`\mathrm{ln}FN\mathrm{ln}(N^2/k_FL)`$, which could be observable if the experiment extends over a larger range of $`N`$.

We mention two other experimentally observable features of the theory presented here. The reduction of the Fano factor described by Eq. (13) is the cumulative effect of many noiseless bands. The appearance of new bands with increasing $`N`$ introduces a fine structure in $`F(N)`$, consisting of a series of cusps with a square-root singularity near the cusp. The second feature is the highly nonuniform spatial extension of open channels, evident in Fig. 1, which could be observed with the STM technique of Ref. Top01 . From a more general perspective the noiseless channels constructed in this paper show that the random matrix approach may be used in ballistic systems only for sufficiently small openings: $`N\sqrt{k_FL}`$ is required. For larger $`N`$ the scattering becomes deterministic, rather than stochastic, and random matrix theory starts to break down.

This work was supported by the Dutch Science Foundation NWO/FOM.


text:
# Nucleon participants or quark participants?

## I Introduction

Charged particle multiplicity densities near mid rapidity in high energy nuclear collisions depend strongly on collision centrality. In order to better understand this dependence and to disentangle pure nuclear effects, this density is often plotted as per nucleon participant pair phobos . Participants are the nucleons that have encountered at least one inelastic collision. The number of participants at a given centrality is usually calculated in the Glauber model either analytically or with the help of a Monte-Carlo algorithm. The charged particle density per participant increases with centrality. In the 5% most central Au+Au collisions it is about 20% larger than that in semi-peripheral collisions (50–70% centrality region), and it is about 50% larger compared to $`pp`$ collisions at the same energy.

The reason for this increase in number of produced particles per participant at midrapidity is still not well understood. The most common explanation of the phenomena involves particle production in hard processes. Hard processes have much smaller cross-sections than soft collisions and depend differently on collision centrality. They scale with the number of binary collisions (the number of collisions the incident nucleon would experience if it were not altered at all while passing through the nucleus). The number of binary collisions increases with centrality faster than the number of participants; this results in an increase of particle production per participant nucleon as centrality increases. In such approaches the particle density is often presented simply as kharzeev :

$$\frac{dN_{ch}}{dy}\alpha N_{Npart}+(1\alpha )N_{Ncoll},$$
(1)
where the parameter $`\alpha `$ is the relative fraction of particles produced in the soft collisions, and $`(1\alpha )`$ is the relative fraction produced in hard collisions. With proper parameters, this fits the data fairly well; see Ref. phobos . Note, however, that in such models, the relative contribution of hard processes is expected to increase with the collision energy. The data seem inconsistent with such an energy dependence.

In the approach proposed in this paper, both nuclei and single nucleons are considered as a superposition of constituent quarks (also often called as “dressed” quarks or valons); there are three such dressed quarks per nucleon. The concept of the constituent quarks has been known for many years (see Refs. constq-first ; hwa1982  and references therein). The constituent quark approach is able to explain many features of hadron-hadron and hadron-nucleus collisions constq-kniga . QCD calculations support the statement that inside a nucleon there are three objects of the size of 0.1 – 0.3 fm (see Ref. constq-shuryak  and references therein); for some recent works using or discussing constituent quarks, see also Refs. constq-recent ; hwa2002 .

In the constituent quark picture, a $`NN`$ collision looks like a collision of two light nuclei. Most often only one $`qq`$ pair interacts, with other quarks being spectators. Only part of the entire nucleon energy is spent for particle production at midrapidity (as $`\sqrt{s_{qq}}\sqrt{s_{NN}}/3`$). The quark spectators form hadrons in the nucleon fragmentation region. In the case of $`AA`$ collisions, more than one quark per nucleon interacts due to the large nucleus size and the possibility for quarks from the same projectile nucleon to interact with different target nucleons. The goal of this study is to find the number of produced particles per participant quark (pair) and to check for its centrality dependence.

## II Calculations of the number of participants

We calculate mean number of nucleon/quark participants using a Monte-Carlo based implementation misko  of the nuclear overlap model nom . We use Woods-Saxon nuclear density profile

$$n_A(r)=\frac{n_0}{1+\mathrm{exp}[(rR)/d]},$$
(2)
with parameters: $`n_0=0.17\mathrm{fm}^3,R=(1.12A^{1/3}0.86^{1/3})\mathrm{fm},d=0.54\mathrm{fm}.`$

In the nuclear overlap model, the mean number of participants in $`A+B`$ collision at impact parameter $`b`$ is given by

$`N_{part,AB}={\displaystyle d^2sT_A(\stackrel{}{s})\{1[1\frac{\sigma _{NN}T_B(\stackrel{}{s}\stackrel{}{b})}{B}]^B\}}`$
$`+{\displaystyle d^2sT_B(\stackrel{}{s})\{1[1\frac{\sigma _{NN}T_A(\stackrel{}{s}\stackrel{}{b})}{A}]^A\}},`$ (3)
where $`T(b)=_{\mathrm{}}^{\mathrm{}}𝑑zn_A(\sqrt{b^2+z^2})`$ is the thickness function; then $`(1\sigma _{NN}T_A(b)/A)^A`$ is the probability for a nucleon to pass through the nucleus without any collision. We use the inelastic $`NN`$ cross section $`\sigma _{NN}=41`$ mb at $`\sqrt{s_{NN}}=130`$ GeV.

The number of participating nucleons for a given centrality can be determined directly using the web interface misko . In order to calculate the number of participating quarks we downloaded FORTRAN code and modified it by increasing the density three times ($`n_0^q=3n_0=51\mathrm{fm}^3`$) and changing $`\sigma _{NN}`$ to $`\sigma _{qq}`$. For the quark-quark cross-section in our calculation, we use two values $`\sigma _{qq}=\sigma _{NN}/9=4.56`$ mb and a somewhat arbitrary value of $`\sigma _{qq}=6`$ mb; the latter was used mostly to illustrate the sensitivity of the results to the value of $`\sigma _{qq}`$. The choice of $`\sigma _{qq}=6`$ mb is not unreasonable since at RHIC energies, approximately 1.2 – 1.3 quarks per nucleon can participate in a single $`NN`$-collision constq-kniga . In principle, $`\sigma _{qq}`$ could be probably as high as 8 mb based on the early estimates of $`(r_q/R_N)^21/5`$ hwa1982 .

Fig. 1a shows the number of the nucleon and quark participants vs. impact parameter. Fig. 1b presents the centrality dependence of the ratio of $`N_{qpart}/N_{Npart}`$. Smooth curves are the polynomial fits to the Monte-Carlo results to smooth out the statistical fluctuations. The ratio $`N_{qpart}/N_{Npart}`$ is used later for the renormalization of the particle pseudorapidity distributions from per nucleon participant to per quark participant.

Fig. 2 presents the comparison of our calculation of $`N_{Npart}`$ in the nuclear overlap model with PHOBOS calculations phobos  based on the HIJING model. The number of participants is plotted as a function of a given fraction of the total cross section. Open circles represent PHOBOS calculations. The nuclear overlap model results (using the Woods-Saxon density profile, the same as used in the HIJING model) are shown by solid symbols. Note a small deviation of our calculations from that of PHOBOS in very central region; this is not important for the conclusion of the current study.

The PHOBOS Collaboration presents their results on centrality dependence of the charged particle pseudorapidity density by plotting it vs. the number of the nucleon participant pairs. In this paper, we continue to use the same quantity for the centrality characterization, but note that in the constituent quark picture, the number of the nucleon participants no longer has the meaning of the number of the particle production sources.

## III Results

The centrality dependence of the charged particle multiplicity per participant pair is shown in Fig. 3. The results per nucleon participant pair are in the upper part of the figure and the results per quark participant pair are shown in the lower part. The original PHOBOS results phobos  on $`dN_{ch}/d\eta `$ per nucleon participant pair (calculated using HIJING) for $`\sqrt{s_{NN}}=130`$ GeV and 200 GeV are shown in solid symbols. In open symbols we also show the same results renormalized for the number of the nucleon participant pairs from our calculations using the nuclear overlap model. Our main result – the centrality dependence of the charged particle pseudorapidity density per quark participant pair – is presented in the same plot.

We observe no, or even slightly decreasing, dependence of $`(dN_{ch}/dy)/N_{qpart}`$ on centrality, with the ratio being dependent only on the energy of the collision. The slight decrease in particle production at midrapidity with increasing centrality could be due either to low values of the constituent quark inelastic cross section used in our calculations or to parton saturation effects.

Note that in the constituent quark picture, $`(dN_{ch}/dy)/N_{qpart}`$ as a function of centrality depends very weakly on the collision energy, as the change in the inelastic cross section is probably less than 5% between $`\sqrt{s_{NN}}=130`$ GeV and 200 GeV.

Hard processes scale with the number of binary collisions. Although it was not necessary to include the contribution of hard processes into our calculation in order to describe the centrality dependence of the charged particle density at midrapidity, we have calculated the number of binary collisions as well: see Fig. 4. Note that the number of binary collisions per participant has a much weaker centrality dependence in the constituent quark approach than it has in the the nucleon participant model.

## IV Summary

We have shown that the particle multiplicity density near mid rapidity scales linearly with the number of constituent quark participant pairs. The experimentally observed increase of $`dN_{ch}/d\eta `$ per nucleon participant pair with centrality in this picture is explained by the relative increase in the number of interacting constituent quarks in more central collisions.

We thank D. Mis̀kowiec for making his code available to us. This work was supported in part by U.S. DOE Grant No. DE-FG02-92ER40713 and Wayne State University Office of the Provost, Global Education Grant.


text:
# Effects of a Collective Spin Resonance Mode on the STM Spectra of D-Wave Superconductors

## Abstract

A high-energy spin resonance mode is known to exist in many high-temperature superconductors. Motivated by recent scanning tunneling microscopy (STM) experiments in superconducting Bi2Sr2CaCu2O8+δ, we study the effects of this resonance mode on the local density of states (LDOS). The coupling between the electrons in a $`d`$-wave superconductor and the resonance mode produces high-energy peaks in the LDOS, which displays a two-unit-cell periodic modulation around a nonmagnetic impurity. This suggests a new means to not only detect the dynamical spin collective mode but also study its coupling to electronic excitations.

A prominent feature in the excitation spectrum of the high-$`T_c`$ superconductors is the “41 meV” collective spin resonance mode, seen by inelastic neutron scattering experiments in most of the cuprate families Rossat-Mignod91 ; Mook93 ; Fong96 ; He01 . The physics of this resonance mode – including its microscopic origin, its connections with other physical properties, as well as its role on superconductivity itself – has been the subject of considerable debate. Given the recent developments of the atomic resolution scanning tunneling microscopy (STM) Hoffman02a ; Hoffman02b ; Howald02a ; Howald02b ; McElroy02 , it is timely to address the possible manifestation of this resonance mode in the local density of states (LDOS). A number of theoretical works Polkovnikov02 ; Podolsky02 ; Han02 ; Kivelson02  have addressed the effect of related spin fluctuations on the LDOS. These works focused on the pinning of the spin fluctuations by impurity: a dynamical spin mode centered around the wavevector $`𝐐`$ leads to a $`2𝐐`$ spatial modulation in the low-energy LDOS. This result is smoothly connected to what happens in the case of a static spin-density-wave ordering Zachar98 ; Zhu02 . However, for the resonance mode – which is sharply peaked at $`𝐐=(\pi ,\pi )`$ – such effects would not be manifested \[since $`2𝐐=(2\pi ,2\pi )`$ is equivalent to $`(0,0)`$\]. There are also works about quasiparticle scattering from nomagnetic impurities Wang02 ; Zhang02 .

In this paper, we show that the coupling of $`d`$-wave quasiparticles to the resonance mode does produce spatial modulations in the LDOS around an impurity. The feature is located at relatively high energies, $`\pm (\mathrm{\Delta }_0+\mathrm{\Omega }_0)`$, where $`\mathrm{\Delta }_0`$ is the maximum superconducting energy gap and $`\mathrm{\Omega }_0`$ the resonance energy. In addition, the wavevector of the LDOS modulation is close to $`𝐐=(\pi ,\pi )`$. Our predictions could be observable by the STM experiments. Such STM studies represent a new means to characterize the coupling between the electronic excitations and the resonance mode. The STM feature we discuss relates to the “peak-dip-hump” structure of the angle resolved photoemission spectroscopy (ARPES) Dessau91 ; Ding96 ; Abanov02 ; Kee02 ; Eschrig00 ; Abanov00 ; the inference about the electron-spin coupling from the ARPES and related spectroscopies is a topic of recent controversy Kee02 ; Abanov02  and we hope that the STM studies we propose will shed new light on this important issue.

We start with a model Hamiltonian describing two-dimensional electrons coupled to a collective spin mode and in the presence of a single-site impurity: $`=_{BCS}+_{sp}+_{imp}`$. Here the BCS-type Hamiltonian for a uniform $`d`$-wave superconductor is given by $`_{BCS}=_{𝐤,\sigma }(\epsilon _𝐤\mu )c_{𝐤\sigma }^{}c_{𝐤\sigma }+_𝐤(\mathrm{\Delta }_𝐤c_𝐤^{}c_𝐤^{}+\mathrm{\Delta }_𝐤^{}c_𝐤c_𝐤)`$, where $`c_{𝐤\sigma }^{}`$ ($`c_{𝐤\sigma }`$) creates (annihilates) a conduction electron of spin $`\sigma `$ and wavevector $`𝐤`$, $`\epsilon _𝐤`$ is the normal state energy dispersion for the conduction electrons, $`\mu `$ the chemical potential, and $`\mathrm{\Delta }_𝐤=\frac{\mathrm{\Delta }_0}{2}(\mathrm{cos}k_x\mathrm{cos}k_y)`$ the $`d`$-wave superconducting energy gap. The coupling between the electrons and the resonance mode is modeled by an interaction term $`_{sp}=g_i𝐒_i𝐬_i`$, where the quantities $`g`$, $`𝐬_i`$, and $`𝐒_i`$ are the coupling strength, the electron spin operator at site $`i`$, and the operator for the collective spin degrees of freedom, respectively. The dynamics of the collective mode will be specified below. The impurity scattering is given by $`_{imp}=U_0_\sigma c_{0\sigma }^{}c_{0\sigma }`$, where without loss of generality we have taken a single-site impurity of strength $`U_0`$ located at the origin, $`𝐫_i=0`$. By introducing a two-component Nambu spinor operator, $`\mathrm{\Psi }_𝐤=(c_𝐤,c_𝐤^{})^T`$, the matrix Green’s function for the $`d`$-wave BCS Hamiltonian $`_{BCS}`$ is determined by $`G_0^1(𝐤;i\omega _n)=\left(\begin{array}{cc}i\omega _n\xi _𝐤& \mathrm{\Delta }_𝐤\\ \mathrm{\Delta }_𝐤& i\omega _n+\xi _𝐤\end{array}\right),`$ where $`\xi _𝐤=\epsilon _𝐤\mu `$ and $`\omega _n=(2n+1)\pi T`$ is the fermionic Matsubra frequency. We have also assumed that the $`d`$-wave pair potential is real.

For a homogeneous system, where only the inelastic scattering of quasiparticles from the collective mode occurs, we calculate the self-energy to the second order in the coupling constant (see Fig. 1b):

$$\widehat{\mathrm{\Sigma }}(𝐤;i\omega _n)=\frac{3g^2T}{4}\underset{𝐪}{}\underset{\mathrm{\Omega }_l}{}\chi (𝐪;i\mathrm{\Omega }_l)G_0(𝐤𝐪;i\omega _ni\mathrm{\Omega }_l),$$

where $`\chi (𝐪;i\mathrm{\Omega }_l)`$ is the dynamical spin susceptibility $`\chi _{ij}(\tau )=T_\tau (S_i^x(\tau )S_j^x(0))`$ and $`\mathrm{\Omega }_l=2l\pi T`$ the bosonic Matsubra frequency. The dressed Green’s function is:

$$\underset{¯}{G}_0^1(𝐤;i\omega _n)=\left(\begin{array}{cc}i\omega _n\xi _𝐤\mathrm{\Sigma }_{11}& \mathrm{\Delta }_𝐤\mathrm{\Sigma }_{12}\\ \mathrm{\Delta }_𝐤\mathrm{\Sigma }_{21}& i\omega _n+\xi _𝐤\mathrm{\Sigma }_{22}\end{array}\right).$$
(1)
The corresponding real-space dressed Green’s function $`\underset{¯}{G}_0(i,j;i\omega _n)`$ is obtained through a Fourier transform with respect to $`𝐫_i𝐫_j`$. For the $`d`$-wave pairing symmetry, one can show that the local Green’s function, $`\underset{¯}{\overset{^}{g}}_0(i\omega _n)=\underset{¯}{G}_0(i,i;i\omega _n)`$ is diagonal. In the presence of a single-site impurity at $`𝐫_i=0`$ with potential strength $`U_0`$, the site-dependent Green’s function can be written in terms of the $`T`$-matrix:

$`G(i,j;E)`$ $`=`$ $`\underset{¯}{G}_0(i,j;E)`$ (2)
$`+{\displaystyle \underset{lm}{}}\underset{¯}{G}_0(i,l;E)\widehat{T}_{lm}(E)\underset{¯}{G}_0(m,j;E).`$
Due to the vertex corrections induced by the coupling to the collective modes (Fig. 1a), the $`T`$ matrix in general contains site-off-diagonal terms. We will first carry out the calculation without the vertex corrections, in which case $`\widehat{T}_{lm}=\widehat{T}\delta _{l,0}\delta _{m,0}`$, with $`\widehat{T}^1=U_0^1\sigma _3\underset{¯}{\overset{^}{g}}_0`$, where $`\sigma _3`$ is the $`z`$-component of the Pauli matrix. The LDOS at the $`i`$-th site, summed over two spin components, is

$`\rho (𝐫_i,E)={\displaystyle \frac{2}{\pi }}\text{Im}G_{11}(i,i;E+i\gamma ),`$ (3)
where $`\gamma =0^+`$.

Up to now, our discussion and formulation are quite general and can be used to study the effects of any dynamic mode once the susceptibility $`\chi `$ is known. We treat the susceptibility in a phenomenological form (based on the inelastic neutron scattering observations), see also Eschrig00 :

$$\chi (𝐪;i\mathrm{\Omega }_l)=\frac{\delta _{𝐪,𝐐}}{2}\left[\frac{1}{i\mathrm{\Omega }_l\mathrm{\Omega }_0}\frac{1}{i\mathrm{\Omega }_l+\mathrm{\Omega }_0}\right],$$
(4)
where we denote the wavevector $`𝐐=(\pi ,\pi )`$ and the mode energy by $`\mathrm{\Omega }_0`$. This form is especially suitable for the optimally doped YBa2Cu3O6+y (YBCO) compounds in the superconducting phase, where the observed neutron resonance peak is almost resolution-limited in energy and fairly sharp in wavevector. The resonance peak in BSCCO is broadened in both energy and wavevector. The finite width in the wavevector space might be important for the ARPES lineshape in general and in particular the understanding of the ARPES spectra away from the $`M`$ points \[$`𝐤=(\pi ,0)`$ and symmetry-related points\] of the Brillouin zone Eschrig00 , but should not change the qualitative conclusion of our work: the LDOS effects we will discuss arise from the fact that the dominant effects of the resonance mode on the single-electron spectral functions occur near the $`M`$ points which is expected to remain to be the case beyond our simplified form for the susceptibility. In addition, given that the peak in BSCCO is still quite sharp in energy, we expect that the main effect of the broadening in energy of the resonance mode is to extend the bias window for the LDOS feature we will discuss. We have also neglected the incommensurate peaks seen in the inelastic neutron scattering experiments in YBCO (the part that disperses “downward” away from the resonance peak) Arai99 ; Fong00 ; Brinckmann99 ; Kao00 , since their spectral weight is significantly smaller than that of the resonance mode. For the normal-state energy dispersion, we use $`\epsilon _𝐤=2t(\mathrm{cos}k_x+\mathrm{cos}k_y)4t^{}\mathrm{cos}k_x\mathrm{cos}k_y`$, where $`t`$ and $`t^{}`$ are the nearest and next-nearest neighbor hopping integral. Unless specified explicitly, the energy is measured in units of $`t`$. We choose $`t^{}=0.2`$ to model the band structure of the hole-doped cuprates. Since the maximum energy gap for most of the cuprate superconductors at the optimal doping is about $`30\text{meV}`$ while the resonance mode energy is in the range between $`35`$ and $`47`$ meV, we take $`\mathrm{\Delta }_0=0.1`$ and $`\mathrm{\Omega }_0=0.15`$ (i.e., $`1.5\mathrm{\Delta }_0`$). The chemical potential is tuned to give an optimal doping value 0.16. To mimic the intrinsic life time broadening, in our numerical calculation we take $`\gamma `$ of Eq. (3) to be $`0.04\mathrm{\Delta }_0`$. A system size of $`N_x\times N_y=1000\times 1000`$ is taken in the numerical calculation.

In the absence of impurities, the density of states is the summation of the spectral function, $`A_𝐤(E)=\frac{2}{\pi }\text{Im}\underset{¯}{G}_{0,11}(𝐤;E+i\gamma )`$, over all wavevectors $`𝐤`$. Fig. 2(a) shows the spectral function at an $`M`$ point of the Brillouin zone. Without the electron-mode coupling, as is well known, the spectral function is peaked at the maximum gap edges $`\pm \mathrm{\Delta }_0`$. As the electron-mode coupling is switched on, new peaks emerge at the energies $`\pm E_1\pm (\mathrm{\Delta }_0+\mathrm{\Omega }_0)`$. (For simplicity, we have neglected the broad “background” part of the single-electron spectral function in our consideration.) The peaks in $`A_𝐤`$ originate from the poles of the Green’s function $`\underset{¯}{G}_{0,11}`$. Note that the weight of the peak at $`\mathrm{\Delta }_0`$ is larger than at $`\mathrm{\Delta }_0`$ because the van Hove singularity is below the Fermi energy. Since the spectral weight of the spin resonance mode \[i.e., $`\text{Im}\chi (𝐪;\omega )`$\] is peaked at $`𝐐`$, the feature of the quasiparticle self-energy is the strongest around the $`M`$ points of the zone because they are connected by $`𝐐`$. The singularity in the quasiparticle self-energy causes additional poles in the Green’s function. As the coupling constant $`g`$ increases, these peaks are shifted to higher energies and, in addition, their spectral weight is enhanced; simultaneously, the weight of the superconducting coherent peaks is reduced to obey the sum rule. The shift of states due to inelastic scattering is expected in DOS and is also expected for scattering off local mode Balatsky02 . Fig. 2(b) plots the density of states and clearly shows that the high energy peaks still occur around $`\pm E_1`$. In other words, the contributions from near the $`M`$ points dominate the wavevector summation for the density of states, reflecting the flat nature of the normal state band near this point. Furthermore, the highly asymmetrical structure in the DOS at energies $`E_1`$ and $`E_1`$ comes from the singular structure in the quasiparticle self-energy. These results, for the clean case, are consistent with earlier studies of the ARPES Dessau91 ; Ding96 ; Abanov02 ; Kee02 ; Eschrig00  and DOS Abanov00 .

We are now in a position to address the LDOS in the presence of a single nonmagnetic impurity. For concreteness, we take the on-site potential $`U_0=100\mathrm{\Delta }_0`$. Fig. 3 shows the LDOS directly at the impurity site, as well as at its nearest neighbor. The near-zero energy resonant state triggered by the quasiparticle scattering from the impurity Balatsky95  is robust against the electron-mode coupling. Our key new results are two-fold. First, the impurity modifies the shape of the spectral features at $`\pm E_1`$, which can now be either a dip or a peak. Second, and more importantly, these high energy features at $`\pm E_1`$ exhibit a spatial dependence. At the impurity site, the LDOS displays a dip at $`E_1`$ but a peak at $`+E_1`$. The behavior is reversed at the site closest to the impurity.

To explore this spatial variation of the LDOS in more detail, we have calculated the LDOS, with the energy fixed at $`E_1`$, in the vicinity of the impurity with and without the mode coupling. In the absence of mode coupling, we obtain the results (not shown) which are consistent with previous studies Balatsky95 : the LDOS exhibits a Friedel oscillation along the diagonals of CuO2 plane but no other non-trivial features. When the mode coupling is turned on, in addition to the Friedel oscillation along the diagonals, the LDOS displays a new type of modulation with a period of $`2a`$ in the wide regions along the bond direction. This new modulation can be more easily seen Fig. 4(a) when the pre-dominant Friedel oscillation is filtered away. It occurs in four disconnected triangles in the field of view.

Also in order to highlight the new modulation, we find it useful to perform a filtered Fourier transform, $`\rho (𝐪,E)=_i^{}\mathrm{e}^{i𝐪𝐫_i}\rho (𝐫_i,E)`$, where $`_i^{}`$ denotes a summation over all the sites in four triangles in Fig. 4(a). The resulting Fourier-transformed image note3  is given in Fig. 4(b), which unambiguously shows that the new feature induced by the coupling to the spin resonance mode has a spatial modulation wavevector $`(\pi ,\pi )`$.

This new type of LDOS modulation with a wavevector close to $`𝐐=(\pi ,\pi )`$ reflects the dominance of the collective mode effect near the $`M`$ points of the Brillouin zone. It follows from Eqs. (2,3), with the form of the $`T`$-matrix under consideration and when the chosen field of view has an inversion symmetry with respect to the impurity site, that the Fourier transformed LDOS for any finite $`𝐪`$ is,

$`\rho (𝐪,E)={\displaystyle \frac{2}{\pi }}\text{Im}{\displaystyle 𝑑𝐩\left[\underset{¯}{G}_0(𝐩+𝐪;E)\widehat{T}(E)\underset{¯}{G}_0(𝐩,E)\right]_{11}}.`$ (5)
At positive energies, the poles shown in Fig. 2 for $`\underset{¯}{G}_0(𝐩;E)`$ are located Dessau91 ; Ding96 ; Abanov02 ; Kee02 ; Eschrig00  at $`E_𝐩`$ and $`E_{𝐩𝐐}+\mathrm{\Omega }_0`$, respectively. (Here, $`E_𝐩\sqrt{\xi _𝐩^2+\mathrm{\Delta }_𝐩^2}`$.) Likewise, $`\underset{¯}{G}_0(𝐩+𝐪;E)`$ is the sum of two poles, one at $`E_{𝐩+𝐪}`$ and the other at $`E_{𝐩+𝐪𝐐}+\mathrm{\Omega }_0`$. At $`E=E_1`$, Eq. (5) is dominated by the convolution between the poles at $`E_{𝐩𝐐}+\mathrm{\Omega }_0`$ and $`E_{𝐩+𝐪𝐐}+\mathrm{\Omega }_0`$. This term in turn is dominated by the contributions corresponding to when both $`𝐩𝐐`$ and $`𝐩+𝐪𝐐`$ are near to the $`M`$ points, leading to a $`\rho (𝐪,E_1)`$ that is peaked near $`(\pi ,\pi )`$.

Similar phase-phase considerations show that the vertex correction terms lead to a similar momentum dependence. In the presence of vertex corrections, the $`T`$-matrix satisfies an integral equation in the wavevector space. The vertex correction to the $`T`$-matrix to the same order ($`g^2`$) of our calculation for the self-energy is shown in Fig. 1(a). It involves a wavevector convolution of a form similar to that given in Eq. (5). We therefore expect note2  an additive contribution to $`\rho (𝐪,E_1)`$ that is also peaked near $`(\pi ,\pi )`$.

Finally, we remark on issues which go beyond the idealized model we have considered so far: (i) In the presence of a gap inhomogeneity at the nanoscale, the mode signature would appear at $`\mathrm{\Delta }_0+\omega _0`$, with a slight smearing. Here, $`\mathrm{\Delta }_0`$ is the spatially averaged superconducting gap. The fact that a well-defined peak-dip-hump structure appears in the break-junction tunneling spectrum Zasadzinski01  implies that the smearing is not too large; (ii) The mode signature depends on the detailed band structure. In the cuprates, the band is flat near the antinodal points (i.e., $`(\pi ,0)`$ etc.), and a mode with a momentum significantly different from $`(\pi ,\pi )`$ will have a weaker effect compared to that of the $`(\pi ,\pi )`$ mode we addressed; (iii) Tunneling matrix elements need to be taken into account in order to understand the detailed spatial variation of the LDOS as observed around zinc impurities in BSCCO Pan00 ; Martin02 . Such a filtering effect, however, will not affect our conclusion on the momentum of the LDOS modulation.

To summarize, we have studied the effects of the magnetic resonance mode on the tunneling spectrum in the presence of a nonmagnetic impurity. The LDOS around the impurity displays resonant features at relatively high energies \[close to $`\pm (\mathrm{\Omega }_0+\mathrm{\Delta }_0)`$\], which modulates in space with a wavevector close to $`(\pi ,\pi )`$. Our prediction can be tested straightforwardly by operating the existing high resolution STM at a relatively high energy window. Such studies should shed considerable new light on the physics of the spin resonance mode, in particular its coupling to the electronic excitations.

We are especially grateful to J. C. Davis for helpful conversations at the early stage of this work. The authors have benefited considerably from discussions with Ar. Abanov, Y. Bang, A. V. Chubukov, K. Damle, and M. Norman. This work was supported by the US DOE (JXZ and AVB), by TcSAM and the NSF Grant No. DMR-0090071 (JS and QS). JXZ also acknowledges the hospitality of the Rice University, where part of this research was carried out.