text:
# A Black Hole in the Superluminal source SAX J1819.3-2525 (V4641 Sgr)1

## 1 Introduction

SAX J1819.3-2525 was discovered as a relatively faint X-ray source independently with the Wide Field Cameras on BeppoSAX on 1999 February 20 (in ’t Zand et al., 1999, 2000) and with the Proportional Counter Array on the Rossi X-ray Timing Explorer on 1999 February 18 (Markwardt, Swank, & Marshall, 1999). The source had two rapid and bright X-ray flares around 1999 September 15. The first flare had a peak intensity of about 4.5 Crab in the 2-12 keV X-ray band (Smith, Levine, & Morgan, 1999) and 5 Crab in 20-100 keV X-ray band (McCollough & Finger, 1999). The second peak, which came about 0.8 days later, reached intensities of 12.2 and 8 Crab in the two bands, respectively. A few days after the giant flares the source was no longer detected in X-rays \[e.g. Wijnands & van der Klis (2000)\].

SAX J1819.3-2525 was detected as a bright radio source (0.4 Jy) on 1999 September 16 (Hjellming et al., 2000), about 16 hours after the second X-ray flare. The radio source was marginally resolved in this first observation, indicating high expansion velocities. Thereafter, the source was detected at radio wavelengths for about three weeks. Hjellming et al. (2000) modelled the radio images and the detailed radio light curves as a combination of the ejection of a relativistic, freely expanding jet and a subsequent ejection of a more slowly decaying, optically thin jet segment. The radio observations did not resolve the moving components of the jets, so as a result one must make an assumption about when the jets were ejected in order infer the proper motion. If the jets were ejected during the brief 4.5 Crab X-ray flare, the apparent proper motion is 0.22 arcseconds per day. Hjellming et al. (2000) concluded the most likely time of the jet ejection was during the rise of the 12 Crab X-ray flare, in which case the apparent proper motion would be 0.36 arcseconds per day. These high proper motions obviously imply large velocities; an apparent proper motion of $`0.22`$ arcseconds per day corresponds to an apparent velocity of $`1.28cd`$, where $`c`$ is the velocity of light and $`d`$ is the distance in kpc.

SAX J1819.3-2525 was already known as an optical variable star (Goranskij, 1978), which was discovered 21 years before the first known X-ray activity. The source had a $`2`$ mag flare in 1978 June, and otherwise showed variability with an amplitude of up to $`1`$ mag about a mean level of $`B14.2`$ (Goranskij, 1990). In spite of differences in the classifications and variability behavior, Goranskij’s variable was confused with GM Sgr, a long-period variable star discovered by Luyten (1927). As a result, many of the papers discussing the optical counterpart of SAX J1819.3-2525 refer to GM Sgr. Hazen et al. (1999) located the plate material used by Luyten and found that Luyten’s variable GM Sgr is actually one arcminute south of Goranskij’s variable star, which is the star at the precise radio coordinates of the X-ray transient. Goranskij’s variable has been given the new variable star designation V4641 Sgr, and we shall use this name for the remainder of the paper. Goranskij (1990) assembled a quiescent light curve using 345 plates taken over a span of about 30 years with the Crimean 40 cm astrograph. He noted a periodicity at a frequency of 2.7151 cycles per day, corresponding to a double-wave period of 0.7365483 days. However, owing to severe aliasing problems, Goranskij (1990) could not rule out 1 day aliases of his preferred period. We demonstrate below that the orbital period is 2.8713 days, which is an alias of the period noted above (i.e. $`1/(0.52.8173)0.71`$ cycles per day).

A network of optical observers had noted an increase in the variability of V4641 Sgr just prior to the giant X-ray flares (Kato et al., 1999). About six days before these flares the source showed a $`1`$ mag modulation with a period of $`2.5`$ days. The optical brightness peaked at $`V=8.8`$ on 1999 September 15.8, and subsequently decayed rapidly to its mean quiescent level within two days. An optical spectrum obtained 1999 September 16.25 showed very strong Balmer emission lines (Djorgovski et al., 1999). The Balmer lines made an abrupt transition from emission to absorption between 1999 September 17 and 19 (Garcia, McClintock, & Callanan, 1999). A spectrum obtained 1999 September 30 (Wagner, 1999) showed only strong Balmer absorption lines and interstellar features. At the time, this was quite confusing since the quiescent counterpart was thought to be an early K-type star based on the spectrum of GM Sgr shown in Downes et al. (1995). However, in hindsight, the Balmer absorption line spectrum makes sense since we now know that Downes et al. (1995) did not observe V4641 Sgr, which in fact contains a late B star. It now appears that the quiescent optical state was reached by 1999 September 19 (i.e. about four days after the X-rays were undetectable).

It was initially thought that V4641 Sgr was relatively nearby based on two arguments. First, if the intrinsic velocity of the jet(s) observed in V4641 Sgr is similar to the velocities observed in GRO J1655-40 and GRS 1915+105 (Mirabel & Rodriguez, 1999), then the large observed proper motion of the jet observed in V4641 Sgr would imply $`d500`$ pc. Second, a K-star in a $`0.7`$ day orbit with an apparent $`V`$ magnitude of $`13.7`$ would have a distance of $`1`$ kpc. We demonstrate below that $`d7.4`$ kpc at the 90% confidence level. Thus the apparent velocity of the 0$`\stackrel{}{\mathrm{.}}`$25 jet(s) seen several hours after the beginning of the X-ray event was $`9.5c`$, assuming the jets were ejected at the start of the giant X-ray flare.

The galactic “microquasars” are excellent laboratories for the study of relativistic jets since they evolve orders of magnitude more quickly than the jets in quasars evolve (Mirabel & Rodriguez, 1999). V4641 Sgr is of special interest because of its very rapid X-ray flaring behavior and its extreme superluminal jets, and because the source is optically bright in quiescence ($`V13.7`$). In this paper we report the results of our spectroscopic observations of V4641 Sgr. The observations and basic data reductions are summarized in Section 2. In Section 3 we establish the basic orbital parameters of the system. The modelling of a high resolution spectrum and the derivation of the secondary star properties are discussed in Section 4. We discuss the astrophysical parameters of the binary system in Section 5 and discuss the implications of our results in Section 6. The paper ends with a short summary in Section 7.

## 2 Observations

We obtained a total of 45 spectra of V4641 Sgr between 1999 September 17 and 1999 October 16 using the FAST spectrograph (Fabricant et al., 1998) on the 1.5 m telescope at the Fred L. Whipple Observatory (FLWO) on Mount Hopkins, Arizona. The 300 g mm-1 grating was used with a Loral $`512\times 2688`$ CCD. About half of the observations were taken with a 1$`\stackrel{}{\mathrm{.}}`$5 wide slit, yielding a spectral resolution of $`4.1`$ Å FWHM. The other observations were taken with a 5$`\stackrel{}{\mathrm{.}}`$0 slit, yielding a spectral resolution of $`6.0`$ Å FWHM. The spectra cover the wavelength range 3650-6800 Å. The exposure times were from three to ten minutes and the spectra were generally obtained under mostly clear skies with seeing between 1 and 2 arcseconds. The spectrograph slit was rotated to maintain an approximate alignment with the parallactic angle, and a wavelength calibration lamp was observed after each exposure.

As we noted above, the source was active until about 1999 September 18. All of the spectra from September 17 and 18 show a broad H$`\alpha `$ emission line, with a full width at zero intensity of up to $`5000`$ km s-1. In contrast, the higher Balmer lines were always in absorption. By September 19, the H$`\alpha `$ line was also in absorption, and apparently the source was near or at its quiescent level. In the present paper we will discuss 14 spectra obtained in quiescence (i.e. between 19 September and 16 October) that were taken through the 1$`\stackrel{}{\mathrm{.}}`$5 wide slit.

A total of forty-seven additional spectra of the source were obtained 2000 June 4-7 using the FORS1 instrument on Antu, which is the first 8.2 m telescope at the European Southern Observatory, Paranal. Thirty-seven of the spectra were taken with the 600B grism and a 0$`\stackrel{}{\mathrm{.}}`$7 wide slit; this combination gives a spectral resolution of 4 Å FWHM and a wavelength coverage of 3366-5751 Å. The other ten spectra were taken with the 600R grism and a 0$`\stackrel{}{\mathrm{.}}`$7 wide slit, yielding a spectral resolution of 3.2 Å FWHM and a wavelength coverage of 5142-7283 Å. The seeing was relatively poor for the VLT ($`1\stackrel{}{\mathrm{.}}5`$), and clouds were present on the night of June 5-6. The exposure times were from 1.5 to 4 minutes depending on the conditions. An atmospheric dispersion corrector was used, so the slit was kept at the default north-south direction. Following the standard procedure at Paranal, the flat-field and wavelength calibration exposures were obtained during the daytime hours with the telescope pointed at the zenith.

Ten additional spectra of the source were obtained 2000 June 16 and 17 with the 3.5 m New Technology Telescope (NTT) at the European Southern Observatory, La Silla, using the blue arm of EMMI. The instrumental configuration consisted of the Tek $`1024\times 1024`$ CCD, a 1$`\stackrel{}{\mathrm{.}}`$0 wide slit, and the grating #11, which is a high efficiency holographic grating with 3000 “grooves” per mm (Wilson et al., 1991). The spectral resolution is 0.43 Å FWHM and the wavelength coverage is 3870.6-4032.8 Å. The slit was aligned with the parallactic angle, and a thorium-argon lamp was observed for wavelength calibration. The weather was clear June 16 and a total of eight 30-minute exposures were obtained. The conditions on June 17 were quite poor, and only two 30-minute exposures were taken.

Finally, seven additional spectra were obtained 2000 August 12 and 13 with the 4.2 m William Herschel Telescope (WHT) at the Observatorio de Roque de los Muchachos, la Palma, Spain. The instrumental configuration was the blue arm of the ISIS double beam spectrograph, the R1200B grating, and the EEV12 $`4096\times 2048`$ CCD. We used a 0$`\stackrel{}{\mathrm{.}}`$8 wide slit, yielding a spectral resolution of 0.58 Å FWHM and a coverage of 3680-4500 Å. The slit was aligned with the parallactic angle, and a copper-argon lamp was observed for wavelength calibration. The weather was photometric, and the seeing was generally 1$`\stackrel{}{\mathrm{.}}`$1 or less.

We reduced all of the spectra using IRAF111IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation.. The standard tasks were used to apply the bias and flat-field corrections, and to extract and wavelength-calibrate the spectra. The bright night sky emission line at $`5578`$ Å was used to make small adjustments to the wavelength scales of the FLWO and VLT spectra. The shifts required to align this feature to a common wavelength (5578.0 Å) were generally less than $`15`$ km s-1, although the VLT spectra from the end of the night of June 7 required shifts on the order of 80 km s-1.

## 3 Orbital Period and Spectroscopic Elements

### 3.1 Spectroscopic Period

We measured the radial velocities of the secondary star using the fxcor task within IRAF, which is an implementation of the technique of Tonry & Davis (1979). Synthetic spectra with $`T_{\mathrm{eff}}=10500`$, $`\mathrm{log}g=3.5`$, and $`V_{\mathrm{rot}}\mathrm{sin}i=123`$ km s-1 (see below) were used as templates for the five different sets of spectra (i.e. FLWO, VLT “blue”, VLT “red”, NTT, and WHT). The cross correlation regions covered the available wavelengths between 3700 and 6847 Å, excluding the interstellar Ca II H and K and Na D lines, the diffuse interstellar bands near 4330 and 5780 Å, and a telluric feature near 6280 Å. The final velocities for the FLWO, the VLT blue, and WHT spectra were insensitive to the exact starting and ending points of the cross correlation regions and to the function used to find the centroid of the cross correlation peak. On the other hand, the velocities of the VLT red spectra and the NTT spectra were sensitive to these parameters since these spectra contain very few lines. We can estimate approximate velocity corrections for the VLT red and NTT spectra by using observations of an A0 template star observed both with the VLT and NTT. However, we find that including the corrected velocities did not improve the resulting spectroscopic parameters since the corrected velocities all have relatively large uncertainties. We therefore excluded the VLT red and NTT velocities from the analysis presented below.

To search for the spectroscopic period we computed a three-parameter sinusoid fit to the velocities for a range of trial periods between 0 and 5 days and recorded the values of the reduced $`\chi ^2`$ for the fits. The free parameters at each trial period are the velocity semiamplitude $`K_2`$, the epoch of maximum velocity $`T_0`$(spect), and the systemic velocity $`\gamma `$. The error bars on the individual velocities were scaled by a factor of 2.65 to yield $`\chi _\nu ^2=1.0`$ at the minimum. Fig. 1a shows the resulting $`\chi _\nu ^2`$ vs. $`P`$ curve. The best fit is at a period of $`P=2.81678\pm 0.00056`$ days ($`1\sigma `$ error). The next best fit occurs at $`P=2.699`$ days, where $`\chi _\nu ^2=2.27`$. This alias period and the others are clearly ruled out by inspection of the folded velocity curves and we adopt the following spectroscopic elements: $`P_{\mathrm{spect}}=2.81678\pm 0.00056`$ days, $`K_2=211.0\pm 3.1`$ km s-1, $`\gamma =107.4\pm 2.9`$ km s-1, and $`T_0(\mathrm{spect})=\mathrm{HJD}2,451,442.523\pm 0.052`$. The resulting optical mass function is then $`f(M)=2.74\pm 0.12M_{}`$ ($`1\sigma `$ errors). Fig. 2 shows the velocities and the best fitting sinusoid, and the spectroscopic elements are listed in Table 1.

### 3.2 Photometric Period

Goranskij (1990) published a light curve of V4641 Sgr obtained from 345 plates taken with the Crimean 0.40 m astrograph between about 1960 and 1990. After it was realized that V4641 Sgr was the optical counterpart of SAX J1819.2-2525, Goranskij reanalyzed the photographic data and made them publically available via the VSNet222http://www.kusastro.kyoto-u.ac.jp/vsnet/Mail/obs26000/msg00925.html. We searched these data for periodicities using both the technique of Stellingwerf (1978) and the Lomb-Scargle Periodogram (Lomb, 1976; Scargle, 1982). We excluded the uncertain measurements and also data from 1978 June, when the source had an optical flare. Fig. 1b shows the Lomb-Scargle periodogram in the period range of 0 to 5 days. Using a Monte Carlo procedure, we computed the power above which the peaks are significant at the $`3\sigma `$ level. We made 10,000 simulated light curves having the same observation times as the actual light curve and with magnitude points normally distributed about the mean magnitude. (By construction, the simulated light curves have the same variance as the actual light curve.) A Lomb-Scargle periodogram was computed for each simulated light curve and the power of the strongest peak was recorded. The $`3\sigma `$ threshold was then determined from the cumulative distribution of the peak powers. This threshold is shown by the dashed line in Fig. 1b. There are several peaks significant at the $`3\sigma `$ level. The peak with the most power is at a trial period of $`P=1.40865`$ days (0.7099 cycles day-1), which is half of the spectroscopic period. All of the other significant peaks are aliases of the main peak (e.g. $`1/11/1.40865=1/3.447`$, etc.) and their sidebands. Fig. 1c shows the $`\mathrm{\Theta }`$ statistic of Stellingwerf (1978) plotted as a function of the trial period in the range $`2.6P3.1`$ days (we used 20 phase bins of width 0.05 to compute the $`\mathrm{\Theta }`$ statistic—this combination yielded the smoothest minima). The minimum $`\mathrm{\Theta }`$ is at a period of $`P=2.81730`$ days. Yearly aliases of this period are also evident. Fig. 1d shows the $`\chi _\nu ^2`$ vs. $`P`$ curve from Fig. 1a plotted over the same period range as in Fig. 1c. The period of $`P=2.81730`$ days is the only trial period where there are minima in both the $`\mathrm{\Theta }`$ statistic and in the $`\chi _\nu ^2`$ vs. $`P`$ curve. Thus the photometric and spectroscopic periods of V4641 Sgr are nicely consistent with each other.

We computed the uncertainty in the photometric period by using another Monte Carlo procedure. We adopted a mean photometric error of 0.1 mag (Goranskij, private communication) and constructed 10,000 simulated data sets by assigning a new magnitude to each time using $`y_{\mathrm{new},i}=y_i+0.1\sigma `$, where $`\sigma `$ is a random Gaussian deviate with variance of unity and a zero mean. For each simulated data set we computed the Lomb-Scargle periodogram and the Stellingwerf $`\mathrm{\Theta }`$ statistic and recorded the resulting derived periods. The uncertainty in the photometric period was taken to be the standard deviation of the distribution of best periods from these simulated data. We find $`P=2.81730\pm 0.00001`$ days. Fig. 3 shows the light curve folded on this period and binned into 20 phase bins. The light curve resembles an ellipsoidal light curve with two maxima of roughly equal heights and two minima of unequal depths per orbital cycle. The amplitude is quite large; the range is about 0.5 mag from the deeper minimum to the maximum.

Goranskij (1990) gives the time of minimum light $`T_0`$(photo) as HJD 2,447,707.454 (no error is quoted). We measured $`T_0`$(photo) by folding the light curve on a dense grid of trial values of $`T_0`$(photo) and recording the $`\chi ^2`$ of an ellipsoidal model fit (the solid line in Fig. 3). The result is $`T_0`$(photo) = HJD 2,447,707.$`4865\pm 0.0038`$. The photometric phase of the time of maximum velocity $`T_0`$(spect) is $`0.75\pm 0.02`$, exactly as expected for an ellipsoidal light curve. (Note that the phase zero point used in Fig. 3 is $`T_0`$(photo)$`+0.5P`$ to conform with our ellipsoidal modelling code.)

## 4 Parameters for the Secondary Star

V4641 Sgr is relatively bright, and as a result one can easily obtain a spectrum with a relatively high signal-to-noise ratio (SNR) and high resolution. The eight NTT spectra from June 16 were Doppler shifted to zero velocity and combined to yield a “restframe” spectrum with a resolution of $`R=9000`$ and a SNR of about 50 per pixel. In the same way, we also made a restframe spectrum using the WHT data. The SNR of the WHT restframe spectrum is about 130 per pixel and is reasonably constant over much of the spectrum. We used synthetic spectra generated from solar metallicity Kurucz models333http://cfaku5.harvard.edu/ to determine three basic parameters for the secondary star, namely its effective temperature $`T_{\mathrm{eff}}`$, surface gravity $`\mathrm{log}g`$, and mean projected rotational velocity $`V_{\mathrm{rot}}\mathrm{sin}i`$. The IDL program SYNPLOT (Hubeny, Lanz, & Jeffrey, 1994) was used to compute the detailed model spectra. SYNPLOT can generate an optical spectrum with arbitrary spectral resolution, wavelength sampling, and rotational broadening. \[Currently SYNPLOT uses an analytic broadening kernel with a linear limb darkening coefficient of $`ϵ=0.6`$ (Gray, 1992)\]. For the NTT spectrum the resolution was taken to be $`0.4333`$ Å and the wavelength sampling was 0.1530 Å, and for the WHT spectra those quantities were 0.5802 Å and 0.2231 Å, respectively. We also have two synthetic spectra kindly provided by Peter Hauschildt computed using his PHOENIX code (Hauschildt & Baron, 1999). These models are fully line-blanketed with magnesium and iron computed in non-LTE (NLTE).

We used the technique outlined in Marsh, Robinson, & Wood (1994) to fit the model spectra to the observed ones. This fitting procedure has the advantage that it allows for the possibility of a (continuum) contribution from the accretion disk. The restframe spectra were normalized to their continuum fits. An input model spectrum characterized by the three parameters $`T_{\mathrm{eff}}`$, $`\mathrm{log}g`$, and $`V_{\mathrm{rot}}\mathrm{sin}i`$ is normalized to its continuum using the same procedure used for the observed spectra: it is scaled by a weight factor $`w`$, and subtracted from the observed spectrum. The scatter in the difference spectrum is measured by computing the reduced $`\chi ^2`$ from a low-order polynomial fit (after the interstellar lines were masked out of the fit). We varied $`T_{\mathrm{eff}}`$, $`\mathrm{log}g`$, and $`V_{\mathrm{rot}}\mathrm{sin}i`$ and looked for the “smoothest” difference spectrum. The Kurucz grid has models for gravities of $`\mathrm{log}g=3.0`$, 3.5, and 4.0 (cgs units) near the temperatures of interest (about 10000 K).

### 4.1 Effective Temperature and Spectral Classification

It quickly became clear using the NTT spectrum that the gravity of the secondary star in V4641 Sgr is quite close to $`\mathrm{log}g=3.5`$ since the $`\mathrm{log}g=3.0`$ models have Balmer lines that are much too narrow and the $`\mathrm{log}g=4.0`$ models have Balmer lines that are much too broad. We therefore fixed the gravity at $`\mathrm{log}g=3.5`$. The temperature is best constrained by the WHT spectrum. The He I line at 4024 Å is somewhat sensitive to the temperature, and we find $`T_{\mathrm{eff}}=10500`$ K. The statistical errors in the temperature and gravity are much smaller than the respective grid spacings (0.5 dex in $`\mathrm{log}g`$ and 250 K in $`T_{\mathrm{eff}}`$). No doubt there are systematic errors as well, but it is beyond our present ability to quantify these errors. We conservatively adopt $`\sigma _T=200`$ K and $`\sigma _{\mathrm{log}g}=0.1`$. The values of $`T_{\mathrm{eff}}`$ and $`\mathrm{log}g`$ we find for V4641 Sgr are not too different from those of the B9 III MK spectral type standard $`\gamma `$ Lyr for which Künzli et al. (1997) list $`T_{\mathrm{eff}}=9970\pm 540`$ and $`\mathrm{log}g=3.50`$; the gravity was estimated from a measurement of the H$`\beta `$ index (Balachandran et al., 1986). We fit synthetic spectra to the observed spectrum of $`\gamma `$ Lyr (obtained with the WHT) and find $`T_{\mathrm{eff}}=10000`$ K, in good agreement with the previous values, and $`\mathrm{log}g=3.0`$, somewhat lower than the previous values. Based on this approximate similarity, we assign a B9III classification to the secondary star in V4641 Sgr. We note that the spectral type of B9III derived here is somewhat earlier than the preliminary type of A2V we announced in IAU Circular 7440. The latter value was arrived at via a visual inspection of the Balmer lines in a low resolution spectrum. We believe the spectral type of B9III derived here is much more reliable since we used an impartial fitting procedure on a spectrum with higher resolution and a higher SNR.

### 4.2 Rotational Velocity of the Secondary

Two other interesting facts also came to light in this process. First, it was apparent that the rotational velocity of V4641 Sgr is quite high ($`0.5K_2`$) since the atomic lines are quite broadened by rotation. Second, some of the metal lines (most noteably Mg II near 4781 Å) are much stronger in V4641 Sgr than they are in the models. We therefore iterated the fits to determine the rotational velocity and the approximate abundances of Mg and other metals.

We started first with the determination of the rotational velocity. We computed model spectra for various values of $`V_{\mathrm{rot}}\mathrm{sin}i`$ (using $`T_{\mathrm{eff}}=10500`$ K and $`\mathrm{log}g=3.5`$) and computed the minimum reduced $`\chi ^2`$ of the polynomial fit to the difference spectra, masking out the few anomalous lines. It is well known that there is a systematic error in the measurement of $`V_{\mathrm{rot}}\mathrm{sin}i`$ caused by the non-spherical shape of the Roche lobe and by the variations in the temperature and limb darkening over the Roche surface. Therefore we made nightly restframe spectra using the two WHT spectra from August 12 and the five WHT spectra from August 13, thus insuring that the restframe spectra are composed of individual spectra with similarly distorted line profiles. Fig. 4a shows the reduced $`\chi ^2`$ as a function of the input $`V_{\mathrm{rot}}\mathrm{sin}i`$ for the two restframe spectra. The minimum reduced $`\chi ^2`$ for the August 13 spectrum (which has the larger SNR) is for $`V_{\mathrm{rot}}\mathrm{sin}i=118`$ km s-1. However, the curve is somewhat noisy, and there is another dip near $`V_{\mathrm{rot}}\mathrm{sin}i=113`$ km s-1. The corresponding curve for the August 12 restframe spectrum has its minimum at $`V_{\mathrm{rot}}\mathrm{sin}i=124`$ km s-1, and another broad minimum near $`V_{\mathrm{rot}}\mathrm{sin}i=130`$ km s-1.

It is somewhat difficult to choose the correct value of $`V_{\mathrm{rot}}\mathrm{sin}i`$ on the basis of these results. We repeated the exercise using phase-specific broadening kernels, computed using the ELC code (Orosz & Hauschildt, 2000). The ELC code uses specific intensities from model atmosphere computations rather than from blackbodies and a one- or two-parameter limb darkening law. Thus the distortions in the broadening kernel due to tidal distortions and to temperature and limb darkening variations are accounted for in detail. For the present problem we used the specific intensities from the Kurucz grid. We computed a broadening kernel for each of the two nights by averaging kernels made for each individual spectrum which went into the restframe average. We used a custom made IRAF routine to compute the convolution. Fig. 4b shows the resulting reduced $`\chi ^2`$ vs. $`V_{\mathrm{rot}}\mathrm{sin}i`$. The two curves are much “cleaner”, and both have well-defined minima near $`V_{\mathrm{rot}}\mathrm{sin}i=123`$ km s-1. Fig. 5 shows the numerical kernels used. In the case of the August 13 spectrum, the analytic kernel scaled by 118 km s-1 has almost the same full width at half maximum as the slightly asymmetric numerical kernel, scaled by 123 km s-1. Thus, to first order, we would expect a correction from 118 to 123 km s-1, which is what we find. Similarly, the analytic kernel scaled by 130 km s-1 approximately matches the full width at half maximum of the numerical kernel for August 12, scaled by 123 km s-1. Thus we adopt $`V_{\mathrm{rot}}\mathrm{sin}i=123`$ km s-1. The curves in Fig. 4b start to become steep at $`119`$ and 127 km s-1, so we adopt a $`1\sigma `$ error of 4 km s-1.

### 4.3 Abundance Anomalies

After determining the rotational velocity we then examined in more detail the metal lines that were unexpectedly strong. SYNPLOT can alter the abundance of any given element and compute an approximate line profile. The resulting line profile is reasonably accurate provided the abundance of the element in question is not altered by more than a factor of $`3`$ from its abundance in the input model atmosphere (I. Hubeny, private communication). Since we currently cannot compute a full grid of the fully line-blanketed Kurucz models with arbitrary abundance patterns, our results should be taken as only indicative. We computed models with altered abundances and compared them to the observed V4641 Sgr spectrum. In searching for elements with possible abundance anomalies, we were guided by the work of Israelian et al (1999) who found that the abundances of the $`\alpha `$-process elements nitrogen, oxygen, magnesium, silicon, and titanium in the F-star secondary of GRO J1655-40 are 6-10 times solar. We concentrated on these same elements in V4641 Sgr. As a check, we also computed models for $`\gamma `$ Lyr with altered abundances. The results are given in Table 2 and in Figs. 6, 7, and 9. Details for each element are provided below.

Taken at face value, the LTE Kurucz models indicate that the magnesium in V4641 Sgr is enhanced by about a factor of five to seven with respect to solar. However, we cannot fit the magnesium lines in $`\gamma `$ Lyr using the LTE Kurucz models (Fig.7), which is an indication that a NLTE treatment is needed. In Figure 8 we show the Kurucz models and the Hauschildt models for V4641 Sgr and $`\gamma `$ Lyr near the Mg II $`\lambda `$4481.22 line. In the case of V4641 Sgr ($`T_{\mathrm{eff}}=10500`$, $`\mathrm{log}g=3.5`$), the NLTE treatment makes almost no difference to the strength of the $`\lambda `$4481.22 line. In contrast, the NLTE treatment makes the $`\lambda `$4481.22 line much stronger in $`\gamma `$ Lyr ($`T_{\mathrm{eff}}=10000`$, $`\mathrm{log}g=3`$). We conclude the unusual strength of the $`\lambda `$4481.22 line in V4641 Sgr is due to an enhancement of magnesium of a factor of five to seven times solar.

The pair of lines at $`\lambda `$4128.04 and 4130.58 due to Si II indicates roughly solar abundance in V4641 Sgr; but as with magnesium, these Si II lines are not well modelled in $`\gamma `$ Lyr. NLTE effects probably are important for Si II as well.

Perhaps the most convincing case is for oxygen. The weak O I line at $`\lambda `$4368.2 observed in V4641 Sgr requires an oxygen enhancement of 3 times solar, whereas the solar abundance model for $`\gamma `$ Lyr fits the line profile in that star. For comparison, Balachandran et al. (1986) measured an oxygen abundance in $`\gamma `$ Lyr of about 1.3 times solar. The unresolved O I triplet near $`\lambda `$6454.4 appears in the VLT spectrum (resolution 3.2 Å FWHM, Fig. 9). Its strength suggests the abundance is even a bit more than three times solar (we did not attempt any further iteration, given the low resolution of the VLT spectrum).

The stellar Ca II line at $`\lambda `$3933.66 in V4641 Sgr was resolved from the interstellar line in the August 12 WHT spectrum. A calcium enhancement of twice solar is needed to fit the line. The Ca II line in $`\gamma `$ Lyr is fit by the solar abundance model.

There are no strong titanium lines that are sensitive to the abundance in the WHT spectrum. The VLT spectrum (resolution 4 Å FWHM) contains two Ti II lines that appeared to be much stronger than the corresponding model lines (Fig. 9). A titanium enhancement of a factor of about 10 can explain the profiles of the $`\lambda 5226.5`$ and $`\lambda `$5129.2 lines, but not the blend at $`\lambda `$5187.

There are no strong nitrogen lines in the WHT spectrum. There is a line due to N I in the VLT spectrum (resolution 3.2 Å FWHM) which is anomolous. A model with a nitrogen enhancement of a factor of 10 solar gives a reasonably good fit to the profile. We also searched the WHT and VLT spectra for suitable sulphur lines, but we did not find any.

To summarize the results of our exploratory abundance analysis, we find the following enhancements: nitrogen is $`10`$ times solar (low resolution spectrum only); oxygen is $`3`$ times solar (low and high resolution spectra); magnesium is enhanced by about five to seven times solar; silicon is roughly solar, but the LTE models may not be adequate; calcium is $`2`$ times solar (high resolution spectrum); and titanium is $`10`$ times solar (low resolution spectrum only). As we stated above, these results should be treated with extreme caution. We do not claim that the lines we examined are the most suitable ones for an abundance analysis, nor do we claim the solar metallicity Kurucz models are the most appropriate. Clearly we need better data (e.g. higher resolution spectra with lines from more than one ionization stage of the elements in question) and better models (e.g. NLTE treatment for all of the elements, etc.).

Finally, we repeated the procedure to find the rotational velocity, using a model with altered oxygen, magnesium, and calcium abundances. We found no change in the value of $`V_{\mathrm{rot}}\mathrm{sin}i`$. Fig. 10 shows the normalized WHT spectrum, the best (normalized) model, and the difference between the two. The agreement between the model and the data is quite good. The accretion disk is faint; it contributes less than 1% of the light. For comparison, we also show the difference between the model for $`\gamma `$ Lyr and the observed spectrum (spectrum e).

## 5 Astrophysical Parameters

Now that we have presented the observational facts concerning V4641 Sgr, we proceed to discuss the astrophysical implications of these observations.

### 5.1 Inclination Constraints

A straightforward upper limit to the inclination comes from the lack of X-ray eclipses. As we noted earlier, the source was active as a relatively weak X-ray source between 1999 February and September and was detected on numerous occasions by the PCA instrument on RXTE. C. Markwardt (private communication) informs us that the source was detected in X-rays several times near the inferior conjunction of the secondary star (i.e. a spectroscopic phase range of $`0.75\pm 0.04`$). Assuming the X-rays come from a relatively small region $`R0.01R_{}`$ centered on the compact object, one can show $`i70.7^{}`$ for $`Q=1.5`$.

In principle, one can determine the inclination by modelling the so-called “ellipsoidal” variations in the light curve of the B-star companion. The B-star fills its critical Roche lobe (there is ongoing mass transfer) and as such is quite distorted. As the distorted star moves about in its 2.8173 day orbit, the projected area on the sky changes, which in turn, gives rise to changes in the observed flux. The net result is that one observes a light curve with two maxima per cycle (corresponding to the quadrature phases when the star is seen “side-on”) and two minima per cycle (when the star is seen “end-on”). The amplitude of the light curve and the details of its shape is a function of the inclination.

In practice, however, there are difficulties which can introduce systematic errors into the ellipsoidal modelling. For example, the accretion disk can contribute a significant amount of flux, thereby “diluting” the ellipsoidal light curve observed from the secondary. From an empirical point of view, the black hole binaries with cool companions (G, K, or M-type companions) are subject to the biggest systematic uncertainties since their light curves change from one observing run to the next (McClintock & Remillard, 1986; Orosz et al., 1996; Haswell et al., 1993; Webb et al., 2000). On the other hand, the two black hole binaries with hot companions (GRO J1655-40 with an F6III companion and 4U 1543-47 with an A2V companion) have light curves which are quite stable, and it is apparent that the underlying ellipsoidal light curve of the secondary star dominates in these two cases. The light curves of GRO J1655-40 have been recently modelled and reasonably tight inclination constraints (to within about 5 degrees at 90% confidence) have been obtained (Greene, Bailyn, & Orosz 2000). Since V4641 Sgr contains a hot companion, it seems reasonable to expect that its light curve will be dominated by the ellipsoidal modulations from the secondary.

The only well-sampled light curve we have access to at the moment is the photometric light curve assembled by Goranskij (1990). Unfortunately, the individual measurements are relatively imprecise (the mean estimated error is 0.1 mag), and the folded and binned light curve (Fig. 3) has relatively large error bars. Nevertheless, the ellipsoidal modulation is apparent in the binned light curve. For the moment we will limit ourselves to some representative models and put off the more detailed modelling until precise CCD light curves in several filter bandpasses become available. We used the ELC code (Orosz & Hauschildt, 2000) to compute the models using specific intensities from the Kurucz grid. The dashed line in Fig. 3 is a Johnson $`B`$-band model for a lobe-filling star ($`T_{\mathrm{mean}}=10500`$ K) with an inclination of $`i=70^{}`$ and a mass ratio of $`Q=1.5`$. The amplitude of the model is too small ($`\chi ^2=35.196`$ for the 20 data points). The solid line in Fig. 3 is the same model, but with a faint accretion disk added (the outer radius is 98% of the Roche lobe radius). The depth of the minimum at phase 0.5 is increased by the partial eclipse of the B-star by the disk. In this case, the amplitude of the model is much closer to what is observed, and the model fit is improved ($`\chi ^2=25.39`$ for the 20 data points). The model disk contributes about 0.25% of the light in the $`B`$-band, and this would be compatible with the observed limits (Sec. 4). To find an approximate lower limit to the inclination, we adjusted the inclination of the model until the $`\chi ^2`$ of the fit increased by 9. The mass ratio and the radius of the disk were not adjusted. We find $`i60^{}`$.

We adopt an inclination range that is uniform in the interval $`60^{}i70.7^{}`$. A partial eclipse of the secondary star is needed to reproduce the relative depths of the minima. These conclusions (apart from the upper limit on $`i`$ which is a function of $`Q`$ only) should be treated with caution, given that the light curve is relatively uncertain.

### 5.2 Binary Parameters

For a Roche lobe filling star it is straightforward to show that

$$\frac{V_{\mathrm{rot}}\mathrm{sin}i}{K_2}=\frac{R_{L_1}(Q)}{a}\left(\frac{1+Q}{Q}\right),$$
(1)
where $`Q`$ is the mass ratio, $`R_{L_1}(Q)`$ is the sphere-equivalent radius of the Roche lobe, and $`V_{\mathrm{rot}}=2\pi R_{L_1}/P`$. We used the ELC code to numerically compute $`R_{L_1}(Q)`$ for $`V_{\mathrm{rot}}=123\pm 5`$ km s-1 and $`K_2=211.0\pm 3.1`$ km s-1; we find $`Q=1.50\pm 0.13`$ (90% confidence).

We used a simple Monte Carlo code to compute other interesting binary parameters (e.g. the component masses, secondary star radius, etc.) and their uncertainties using four input quantities, namely the period ($`P=2.81730\pm 0.00001`$ days, $`1\sigma `$), the velocity semiamplitude of the secondary star ($`K_2=211.0\pm 3.1`$ km s-1, $`1\sigma `$), the projected rotational velocity of the secondary star ($`V_{\mathrm{rot}}\mathrm{sin}i=123\pm 4`$ km s-1, $`1\sigma `$), and the inclination ($`60.0^{}i70.7^{}`$, uniform distribution). The nominal masses of the black hole and secondary star are $`9.6M_{}`$ and $`6.5M_{}`$, respectively; the precise masses and other key results are summarized in Table 3.

## 6 Discussion

The mass of the compact object in V4641 Sgr ($`9.6M_{}`$) is well above the maximum mass of a stable neutron star, which is usually taken to be $`3M_{}`$ (Chitre & Hartle, 1976). We therefore conclude that V4641 Sgr contains a black hole. V4641 Sgr is also interesting in several other respects. The B-star secondary is by far the most massive, the hottest, and the most luminous secondary of the ten dynamically confirmed transient black hole binaries. The total mass of the system, $`16.2M_{}`$, is also largest among the dynamically confirmed black hole binaries. The next most massive system is V404 Cyg whose total mass is about $`12.7M_{}`$ (Shahbaz et al., 1994). The other eight confirmed transient black hole binaries have total masses which are probably below $`10M_{}`$, although the uncertainties are large in many cases (Bailyn et al., 1998, and cited references).

The secondary star in V4641 Sgr appears to be in an unusual evolutionary state. To help illustrate this, we plot in Fig. 11ab the radius of the secondary and the luminosity of the secondary as a function of its assumed mass. Since the density of a Roche lobe filling star depends only on the orbital period of the binary to a good approximation, these curves are essentially independent of the inclination or the mass ratio of the binary. We also show the radius and luminosity of single stars, taken from the “Geneva” stellar evolution models (Schaller et al., 1992). For our adopted mass range of $`5.49M_28.14M_{}`$ (90% confidence), the secondary star is both smaller in size and significantly underluminous compared to single stars of similar masses. An underluminous secondary star usually indicates that the star went through a phase of thermally unstable mass transfer. A more detailed discussion of the evolutionary history of this star is beyond the scope of the current paper.

Two remarkable features of V4641 Sgr, which are discussed below, are a consequence of its great distance. For an assumed secondary star mass, we know precisely the corresponding radius. The corresponding luminosity is also precisely known since the temperature is well known from the spectrum. The distance then follows from the apparent $`V`$ magnitude and the extinction. Photometry of the source was obtained from CTIO 2000 July 23.1 (UT). We find $`B=14.00\pm 0.02`$ and $`V=13.68\pm 0.02`$. The intrinsic colors of this star are all close to zero, hence $`E(BV)=0.32`$. The star is of course variable, so for now we adopt an uncertainty on the $`B`$ and $`V`$ magnitudes of 0.1 mag. Fig. 11c shows the distance as a function of the assumed secondary star mass. Assuming $`T_{\mathrm{eff}}=10500\pm 200`$ K and $`5.49M_28.14M_{}`$ (90% confidence), we adopt $`7.40d12.31`$ kpc (90% confidence).

As we noted earlier, V4641 Sgr was marginally resolved in the first radio map made after the 12 crab X-ray flare. The proper motion of the jet is uncertain since it is not known when the jet was ejected. Hjellming et al. (2000) give three possibilities for the jet proper motion, based on known X-ray/radio correlations in other sources:

* $`\mu _{\mathrm{jet}}=0.22`$ arcseconds per day for an ejection at the start of the 4.5 crab X-ray flare, which implies $`V_{\mathrm{app}}9.47c`$ for $`d>7.4`$ kpc;
* $`\mu _{\mathrm{jet}}=0.36`$ arcseconds per day for an ejection at the start of the 12 crab X-ray flare (this is their preferred value), which implies $`V_{\mathrm{app}}15.4c`$ for $`d>7.4`$ kpc;
* $`\mu _{\mathrm{jet}}=1.1`$ arcseconds per day for an ejection at the end of the 12 crab flare, which implies $`V_{\mathrm{app}}47.7c`$ for $`d>7.4`$ kpc.

For comparison, all of the other Galactic sources with relativistic jets have apparent velocities of less than $`2c`$ (Mirabel & Rodriguez, 1999). Indeed, this source might be characterized as a galactic “microblazar”, since the jet velocity is similar to what is typcially seen in blazars (Jorstad et al., 2001). It is worthwhile to point out that there are other interpretations of the radio observations. For example, Chaty, et al. (2001) note that the long-lived radio remnant seen in Fig. 5c of Hjellming et al. (2000) did not show any movement between 1999 September 16 and 24. This raises the possibility that the radio emission was due to an interaction of ejecta with the surrounding interstellar medium at a distance of $`0.25d`$ AU from the source, where the distance $`d`$ is measured in parsecs. At a distance of 9.6 kpc, an angular size of 0.25 arcseconds corresponds to a physical size of 2400 AU. The corresponding light crossing time is 13.9 days, so the jet would have had to been ejected on $``$ September 2 at the latest in order to have interacted with matter at a distance of 2400 AU. Optical observers coordinated through the VSNet noted possible activity as early as August 8444http://www.kusastro.kyoto-u.ac.jp/vsnet/Mail/vsnet/msg01879.html, so a jet ejection as early as September 2 cannot be ruled out.

The apparent proper motion of the approaching jet depends on the intrinsic jet velocity $`\beta =v/c`$, jet angle $`\theta `$, and the distance $`d`$ (Mirabel & Rodriguez, 1999)

$$\mu _a=\frac{\beta \mathrm{sin}\theta }{1\beta \mathrm{cos}\theta }\frac{c}{d}.$$
(2)
Fig. 12 shows a contour plot of $`\mu _a`$ as a function of $`\theta `$ and $`\beta `$ for assumed distances of $`d=7.4`$ kpc (our lower limit at 90% confidence) and $`d=9.6`$ kpc (our preferred value). A proper motion of $`\mu _a>0.22`$ arcseconds per day requires $`\beta >0.9945`$ and $`\theta <12^{}`$ for $`d=7.4`$ kpc. The bulk Lorentz factor in this case is $`\mathrm{\Gamma }=(1\beta ^2)^{1/2}>9.5`$. The combination of $`d=9.6`$ kpc and $`\mu _a=0.36`$ arcseconds per day gives $`\beta >0.999`$ ($`\mathrm{\Gamma }>22.4`$) and $`\theta <6^{}`$. In any case, the bulk motions were extremely relativistic, and the jet in V4641 Sgr was nearly aligned with the line of sight. On the other hand, the large inclination indicates that the orbital plane is nearly perpendicular to the line of sight. If the jet is aligned with the spin vector of the black hole, then it seems likely that the spin vector of the black hole is misaligned (by a wide margin) with the orbital angular momentum vector. The inclination of the transient jet observed during the 1994 outburst of GRO J1655-40 was about $`85^{}`$ (Hjellming & Rupen, 1995), whereas the orbital inclination is close to $`70^{}`$ (Orosz & Bailyn, 1997; van der Hooft et al., 1998; Greene, Bailyn, & Orosz, 2000), indicating a misalignment of about $`15^{}`$. Thus if the inclination of the jet can be taken as the inclination of the black hole spin vector, then misaligned black hole spin vectors may be a common feature of these systems.

The peak isotropic X-ray luminosity observed during the large X-ray flare in the 2-10 keV band is $`L_{\mathrm{max}}=5\times 10^{36}(d/0.4\mathrm{kpc})^2`$ erg s-1 (Smith, Levine, & Morgan, 1999). Using our adopted distance we find $`\mathrm{log}L_{\mathrm{max}}=39.46\pm 0.23`$ (90% confidence) for the 2-10 keV band. The luminosity in the 10-100 keV band is comparable to the luminosity in the 2-10 keV band, so the bolometric luminosity will be a factor of $`2`$ higher. Thus, rather than the main outburst being substantially sub-Eddington, as initially thought (in ’t Zand et al., 2000), the outburst was actually super-Eddington, since $`\mathrm{log}L_{\mathrm{edd}}=39.10`$ for a mass of $`9.61M_{}`$. The bright soft X-ray transient 4U 1543-47 had a peak X-ray flux (1.5-3.8 keV) that was at least three times the Eddington limit (Orosz et al., 1998), so the super-Eddington peak X-ray flux observed for V4641 Sgr is by no means unique.

Fender & Kuulkers (2000) have recently shown that the peak X-ray flux $`P_X`$ observed in the outbursts of the black hole transients is positively correlated with the peak observed radio flux $`P_R`$. If all black hole X-ray transients produce a mildly relativistic jet (five sources have been resolved at radio wavelengths), then the existence of the $`P_X`$/$`P_R`$ correlation implies that on average the bulk Lorenz factors in the jets are relatively small (e.g. $`\mathrm{\Gamma }5`$) since very large values of $`\mathrm{\Gamma }`$ would result in many systems with radio emission beamed out of the line of sight and hence unobservable (Fender & Kuulkers, 2000). These authors point out that there is nothing unusual about V4641 Sgr in terms of its peak X-ray and peak radio flux. The radio emission observed in V4641 Sgr was strongly beamed, and the fact that V4641 Sgr sits on the $`P_X/P_R`$ correlation suggests that the X-rays may also have been beamed and that the peak outburst luminosity was not necessarily super-Eddington. X-ray beaming may also be responsible for the super-Eddington outburst of 4U 1543-47 and for the apparently large numbers (twelve or more) of the binary sources in M82 and Cen A with luminosities of up to $`10^{41}`$ ergs s-1 (Prestwich, 2000). These sources may be $`10M_{}`$ black holes with beamed X-ray emission rather than $`100M_{}`$ black holes emitting isotropically at the Eddington limit.

Finally, we note that we can make an independent estimate of the distance if we assume that the large observed radial velocity of the binary (Table 1) is due entirely to differential galactic rotation. For a source at a galactic longitude of $`\mathrm{}=6.77`$ and galactic latitude of $`b=4.79`$ we find $`d>7`$ kpc using the rotation curve given in Fich, Blitz, & Stark (1989) and the standard IAU rotation constants of $`R_0=8.5`$ kpc and $`\mathrm{\Theta }_0=220`$ km s-1. This lower limit is in good agreement with the lower limit we derived above. The space motion of the binary relative to its local standard of rest $`V_{\mathrm{corr}}`$ can be computed from the distance we derived above and the observed $`\gamma `$ velocity. However, $`V_{\mathrm{corr}}`$ is unfortunately rather uncertain owing to the binary’s close proximity to the galactic center. The distribution of $`V_{\mathrm{corr}}`$ values computed using the simple Monte Carlo code mentioned above is roughly uniform in the range $`80V_{\mathrm{corr}}80`$ km s-1.

Although our basic picture of V4641 Sgr is secure, there are many useful follow-up observations that can be done. The photometric light curve has relatively large errors; and light curves with greater statistical precision and for several bandpasses should be obtained. Our tentative conclusion is that the secondary star is partially eclipsed by the disk. It should be easy to confirm this with CCD light curves. The derived properties of the binary (e.g. component masses, etc.) are somewhat sensitive to the value of $`V_{\mathrm{rot}}\mathrm{sin}i`$. Thus another measurement of the rotational velocity of the B-star secondary should be obtained, preferably near the time of the inferior conjunction of the B-star, when the systematic bias introduced by the tidal distortion of the secondary star is minimized. Finally, since V4641 Sgr is so bright, one can easily obtain a high resolution, high signal-to-noise spectrum and do a more detailed abundance analysis. This was done for the black hole binary GRO J1655-40 by Israelian et al (1999) who reported a large overabundance of oxygen, magnesium, silicon, and sulphur, which they attributed to the capture of supernova ejecta from the progenitor of the present-day black hole. We find that the $`\alpha `$-process elements nitrogen, oxygen, calcium, and titanium appear to be overabundant with respect to solar in V4641 Sgr. However, we do not find the same levels of enrichment in V4641 Sgr that were found for GRO J1655-40. The amounts of the various elements synthesized in a supernova explosion depend on the mass of the progenitor star and on the energy of the explosion (Nomoto et al., 2000). Thus the differences in the levels of enrichment of the $`\alpha `$-process elements in GRO J1655-40 and V4641 Sgr may reflect the differences in the supernovae progenitors.

## 7 Summary

We have presented the results of our spectroscopic campaign on the fast X-ray transient and superluminal jet source V4641 Sgr. The spectroscopic period is $`P_{\mathrm{spect}}=2.81678\pm 0.00056`$ days and the radial velocity semiamplitude is $`K_2=211.0\pm 3.1`$ km s-1. The optical mass function is $`f(M)=2.74\pm 0.12M_{}`$. For the secondary star we measure $`T_{\mathrm{eff}}=10500\pm 200`$ K, $`\mathrm{log}g=3.5\pm 0.1`$, and $`V_{\mathrm{rot}}\mathrm{sin}i=123\pm 4`$ km s-1 (1$`\sigma `$ errors).

The photometric period measured from an archival photographic light curve is $`P_{\mathrm{photo}}=2.81730\pm 0.00001`$ days. The light curve folded on the photometric phase resembles an ellipsoidal light curve. Modelling of this light curve indicates a high inclination angle of $`i60^{}`$. The lack of X-ray eclipses implies an upper limit to the inclination of $`i70.7^{}`$. Using these inclination limits and the above spectroscopic parameters, we find a compact object mass in the range $`8.73M_111.70M_{}`$ (90% confidence). This mass range is well above the maximum mass of a stable neutron star and we conclude that V4641 Sgr contains a black hole.

The secondary is a late B-type star which has evolved off the main sequence. It is in a peculiar evolutionary state since both its radius and its luminosity are much smaller than the corresponding values for single stars with a similar mass.

Finally, we find a distance in the range $`7.40d12.31`$ kpc (90% confidence), which is at least a factor of $`15`$ larger than the initially assumed distance of $`500`$ pc. The peak X-ray luminosity was super-Eddington, and the apparent expansion velocity of the radio jet was $`9.5c`$.

We thank the unknown referee for several helpful comments which improved this paper, and Peter Hauschildt for providing the model spectra. We acknowledge the use of the PERIOD package developed by Vik Dhillon from the Starlink Software Collection, and the SYNPLOT spectrum synthesis programs written by Ivan Hubeny. It is a great pleasure to thank the numerous people at the four different observatories who made this project a success. P. Berlind and M. Calkins were the FLWO observers. M. Chadid and D. Hutsemekers were the support astronomers for the VLT observations. P. Leisy was the support astronomer for the NTT. Leisy and the rest of the NTT team headed by O. Hainaut went beyond the call of normal duty by reconfiguring EMMI on short notice. N. Walton was the support astronomer at the WHT. We acknowledge useful discussions with Ivan Hubeny, Vitaly Goranskij, Craig Markwardt, Peter Hauschildt, Frank Verbunt, Norbert Langer, Marten van Kerkwijk, and the late Bob Hjellming, who will be greatly missed.


format: text

text:
# Scaling of hadron masses and widths in thermal models for ultra-relativistic heavy-ion collisions

By means of a simple rescaling, modifications of hadron masses and widths are incorporated into the thermal analysis of particle ratios in ultra-relativistic heavy-ion collisions. We find that moderate, up to $`20\%`$, changes of hadron masses do not spoil the quality of the fits, which remain as good as those obtained without modifications. The fits with the modified masses yield, however, modified values of the optimal temperature and baryon chemical potential. In particular, with decreasing masses of all hadrons (except for pseudo-Goldstone bosons) the fitted values of the temperature and the baryon chemical potential are lowered, with the change approximately proportional to the scaling of masses. In addition, we find that the broadening of the hadron widths by less than a factor of two practically does not affect the fits.

Keywords:  in-medium properties of hadrons, ultra-relativistic heavy-ion collisions

The H. Niewodniczański Institute of Nuclear Physics,
ul. Radzikowskiego 152, PL-31342 Kraków, Poland

PACS: 25.75.Dw, 21.65.+f, 14.40.-n

Recent theoretical studies  show that the hadronic yields and ratios in ultra-relativistic heavy-ion collisions can be described well in the framework of simple thermal models. The thermodynamic parameters obtained in this approach define the so-called chemical freeze-out point, i.e., a stage in the evolution of the hadronic system when the “chemical composition” is established. Thermal-model fits to the SPS data show that the temperature at the chemical freeze-out, $`T_{chem}`$, as well as the baryon density, are large. One typically obtains $`T_{chem}`$ 170 MeV, which is close to the expected critical value for the deconfinement/hadronization phase transition. In this situation one may expect that hadron properties at the chemical freeze-out are strongly modified by the presence of the hadronic environment. Indeed, such modifications are predicted by model calculations , which inter alia helps to explain the low-mass dilepton enhancement observed in the CERES  and HELIOS  experiments. In this paper we incorporate possible modifications of hadron masses and widths into thermal analysis of particle ratios. We generalize the results of Refs.  where the problem was studied without refitting thermodynamical parameters.

In the first part we include only the mass modifications and calculate the particle densities from the ideal-gas expression222The use of Eq. (1) is valid when the in-medium hadrons can be regarded as good quasi-particles. A thermodynamically consistent approach has been constructed so far only for the lowest multiplets of hadrons . At SPS energies, however, it is crucial to include all hadrons with masses up to (at least) 1.8 GeV. For such a complicated system, a thermodynamically consistent approach is not at hand at the moment.

$$n_i=\frac{g_i}{2\pi ^2}\underset{0}{\overset{\mathrm{}}{}}\frac{p^2dp}{\mathrm{exp}\left[\left(E_i^{}\mu _{chem}^BB_i\mu _{chem}^SS_i\mu _{chem}^II_i\right)/T_{chem}\right]\pm 1},$$
(1)
where $`g_i`$ is the spin degeneracy factor of the $`i`$th hadron, $`B_i,S_i,I_i`$ are the baryon number, strangeness, and the third component of isospin, and $`E_i^{}=\sqrt{p^2+\left(m_i^{}\right)^2}`$ is the energy. The quantities $`\mu _{chem}^B,\mu _{chem}^S`$ and $`\mu _{chem}^I`$ are the chemical potentials enforcing the appropriate conservation laws. In standard thermal-model fits Eq. (1) is used with the vacuum masses, $`m_i^{}=m_i`$. The in-medium masses, $`m_i^{}`$, may depend on temperature and density in a complicated way. In order to explore possible different behavior of in-medium masses and, at the same time, keep simplicity, we do our calculations with the meson and baryon masses rescaled by the two universal parameters, $`x_M`$ and $`x_B`$, namely

$$m_M^{}=x_Mm_M,m_B^{}=x_Bm_B.$$
(2)
An exception from this rule are the masses of pseudo-Goldstone bosons ($`\pi ,K,\eta `$) which we keep constant.

Equation (1) is used to calculate the “primordial” density of stable hadrons and resonances at the chemical freeze-out. The final (observed) multiplicities receive contributions from the primordial stable hadrons, and from the secondary hadrons produced by decays of resonances after the freeze-out. We include all light-flavor hadrons listed in the newest review of particle physics , with a few exceptions for the cases where the properties of a listed particle are ambiguous or not known sufficiently well. Ref.  is also used to determine the branching ratios. We neglect the finite-size and excluded volume corrections: the former are negligible, whereas the latter do not affect the particle ratios.333 The last property is due to the fact that the overwhelming majority of hadrons is heavy and may be treated as classical particles. In this case we may replace the Fermi-Dirac (Bose-Einstein) distribution function in (1) by the Boltzmann distribution function. Thus, for equal excluded volumes of all hadrons the excluded volume corrections factorize and cancel out in the ratios. We have checked that the use of the classical statistics changes our results by a few percent only.

For given values of $`x_M`$ and $`x_B`$ we fit the temperature, $`T_{chem}`$, and the baryonic chemical potential, $`\mu _{chem}^B`$, by the minimization the expression $`\chi ^2=_k\left(R_k^{exp}R_k^{therm}\right)^2/\sigma _k^2`$ , where $`R_k^{exp}`$ is the $`k`$th measured ratio, $`\sigma _k`$ is the corresponding error, and $`R_k^{therm}`$ is the same ratio as determined from the thermal model. The strangeness chemical potential, $`\mu _{chem}^S`$, and the isospin chemical potential, $`\mu _{chem}^I`$, are determined from the requirements that: i) the initial strangeness of the system is zero, and ii) the ratio of the baryon number to the electric charge is the same as in the colliding nuclei.

In Fig. 1 we plot our results obtained for the experimental ratios for Pb + Pb collisions at SPS, as compiled in Ref. . In the case $`x_M=x_B=1`$ we find: $`T_{chem}=`$ 169 MeV, $`\mu _{chem}^B=`$ 250 MeV, $`\mu _{chem}^S=`$ 65 MeV, and $`\mu _{chem}^I=`$ -9 MeV. These values are in good agreement with those found in Ref. : $`T_{chem}=`$ 168 MeV, $`\mu _{chem}^B=`$ 266 MeV, $`\mu _{chem}^S=`$ 71 MeV, and $`\mu _{chem}^I=`$ -5 MeV. In Fig. 1 a) we give our values of $`\chi ^2`$. One can observe that a small decrease of the meson and baryon masses, $`x_M=x_B0.9`$, leads to a slightly better fit with the corresponding smaller values of the temperature and the baryon density, as shown in Figs. 1 b) and 1 c). To check the significance of this result we have also analyzed Si + Au collisions at AGS, and S + Au collisions at SPS. In these two cases $`\chi ^2`$ has a flat minimum at $`x_Mx_B1`$. We thus conclude that moderate modifications of hadron masses, say by 20%, do not spoil the quality of thermal fits. However, with the modified masses the thermodynamic parameters characterizing the fits do change. For example, if we rescale both meson and baryon masses (except for Goldstone bosons) in the same way, $`x=x_M=x_B`$, the temperature and the chemical potentials are to a very good approximation also rescaled by $`x`$. This follows from the fact that we study a system of equations which is invariant under rescaling of all quantities with the dimension of energy. If we allowed also for the changes of the masses of the Goldstone bosons, the thermodynamic parameters would scale exactly as $`T_{chem}(x)=xT_{chem}(x=0)`$ and $`\mu _{chem}(x)=x\mu _{chem}(x=0)`$. In this case $`\chi ^2`$ remains constant, independently of $`x`$. For fixed values of the Goldstone-boson masses, the scale invariance is broken, and $`\chi ^2`$ varies weakly with $`x`$, as shown in Fig. 1 a).

To account for finite in-medium widths, $`\mathrm{\Gamma }_i^{}`$, of the resonances we generalize Eq. (1) to the formula

$`n_i={\displaystyle \underset{M_0^2}{\overset{\mathrm{}}{}}}𝑑M^2{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}𝑑p{\displaystyle \frac{1}{\pi N}}{\displaystyle \frac{m_i^{}\mathrm{\Gamma }_i^{}}{(M^2(m_i^{})^2)^2+(m_i^{}\mathrm{\Gamma }_i^{})^2}}`$
$`\times {\displaystyle \frac{g_i}{2\pi ^2}}{\displaystyle \frac{p^2}{\mathrm{exp}\left[\left(\sqrt{M^2+p^2}\mu _{chem}^BB_i\mu _{chem}^SS_i\mu _{chem}^II_i\right)/T_{chem}\right]\pm 1}},`$ (3)
where $`N`$ is the normalization of the relativistic Breit-Wigner function, $`N=\frac{1}{2}+\frac{1}{\pi }\mathrm{arctan}\left[((m_i^{})^2M_0^2)/(m_i^{}\mathrm{\Gamma }_i^{})\right]1`$. The integral over $`M^2`$ is taken to start at the threshold $`M_0^2`$ corresponding to the dominant decay channel. In the limit $`\mathrm{\Gamma }_i^{}0`$ Eq. (3) obviously reduces to formula (1).

In order to analyze the effect of broadening of hadron widths we introduce the parameter $`y`$ in such a way that

$$\mathrm{\Gamma }_i^{}=y\mathrm{\Gamma }_i.$$
(4)
Here $`\mathrm{\Gamma }_i`$ are the vacuum widths, hence the case $`y=1`$ corresponds to the physical widths as measured in the vacuum, and the case $`y=0`$ represents the situation when the widths are neglected (our previous analysis based on Eq. (1)). In Fig. 2 we show the results of our fitting procedure. We observe that the inclusion of the vacuum widths does not change the value of $`\chi ^2`$, and the values of $`T_{chem}`$ and $`\rho _{chem}^B`$. An increase of the widths by a factor of 2 has also little effect. Only for larger modifications of the widths the fit gets worse.

In conclusion we state that the thermal analysis of particle ratios allows for scaling of hadron masses. Such changes, however, result in modifications of the thermodynamical parameters for which the fits are optimal. In particular, lowering of all the masses leads to a smaller values of $`T_{chem}`$ and $`\rho _{chem}^B`$. This is a desired effect, since $`T_{chem}`$170 MeV is large and may correspond to quark-gluon plasma rather than to a hadron gas. Our study of the modifications of the hadron widths shows that they have small impact on the ratios.


text:
# Localization properties of groups of eigenstates in chaotic systems

## I Introduction

The investigation of the quantum manifestations of classical chaos is at present a very active field of research . The relationship between quantum dynamics and classical invariants, mainly periodic orbits (PO) and their associated manifolds, is not satisfactorily understood yet . In this respect, an important achivement would be to disentangle the time evolution of a quantum system whose classical analog is chaotic, explaining this dynamics in terms of classical structures.

Fifteen years ago, a seminal step was done by Heller , who considered the evolution of a wave packet localized on an unstable PO. He found that the dynamics in the short time limit are controlled by the motion around the orbit and its Lyapunov exponent and is well understood in terms of semiclassical expressions .

For longer times the dynamical mixing induced by the chaotic behavior of the motion makes the situation more complicated, and one single orbit and the linearized dynamics around it is not enough to account for the evolution of the wave packet. Some time later, Tomsovic and Heller  showed that even for those long times, for which classical fine structure had developed on a scale much smaller than $`\mathrm{}`$, the semiclassical propagation of the packet can be carried out with good precission, by computing the corresponding correlation function, $`C_{\mathrm{scl}}(t)`$, as a sum of contributions of the homoclinic excursions of the PO. This procedure has nevertheless some drawbacks, since for example the number of orbits that is necessary to include in the calculation to obtain converged results grows dramatically with time.

A different approach, that maintains an interpretation of wave mechanics based on the more simple classical objects represented by POs has been presented elsewhere . In this theory  all quantum information of a bounded chaotic Hamiltonian system can be obtained using a number of POs that grows only linearly with the Heisenberg time.

Another point worth considering is the importance and/or convenience of using wave functions which are “dynamically adapted” in this type of studies . By “dynamically adapted”  we mean structures that contain information not only about the PO, but also about the linearized dynamics around the orbit. When considered in phase space, these structures appear no only localized over the PO, but are also spread over the corresponding manifolds.

The construction of functions of this kind can be carried out using different approaches. In Ref.  the necessary information is obtained by Fourier transforming the short–time exact quantum dynamics, so that wave functions averaged over the PO path are obtained. Kaplan and Heller  used a sum of coherent Gaussian wave packets centered along the PO to obtain the desired function. Finally, Vergini and Carlo  constructed semiclassically resonances along POs with minimum energy dispersion.

In spite of the differences existing among these functions, they are all suitable to efficiently investigate the quantum manifestations of the classical phase space complexity in classically chaotic systems.

The purpose of this paper is threefold. In the first place, we focus on the analysis of the topological characteristics (i.e. distribution of quantum probability) in phase space of this type of wave functions highly localized on an unstable PO. In the second place, special attention will be paid to the relationship existing between the quantum dynamics of this function and that of other classical POs of the system located in their vicinity, in order to explore the existence of dynamical connections and active interplay between them. Finally, another interesting aspect of these functions is their revivals, that take place in some cases for a surprisingly large number of cycles without showing any visible sign of dispersion . This fact, which is particularly surprising in view of the highly chaotic character of the problem that will be considered, is due to the special dynamical properties of these functions. Summarizing, it should be stressed that the most important result of our work is that, proceeding in this way, only a few short POs are necessary to adequately explain the quantum dynamics of this type of systems.

For our study we will use localized functions as defined in Ref. , applying them to the study of the dynamics in the Bunimovich stadium billiard , a system which is known to behave ergodically from the classical point of view.

The organization of the paper is as follows. In Sect. II we describe in detail how the method of Ref.  can be used to construct the desired wave functions, highly localized on a given unstable PO and their associated manifolds. Numerical results concerning the topological characteristics, both in configuration and phase space, of these functions and the corresponding time evolution will be presented and analyzed in Sects. III and IV, respectively. Finally, the main conclusions derived from our work are presented in Sect. V.

## II Constructing Localized Wave Functions

In this section we briefly describe the method  that will be used to construct states initially located on a given PO and the corresponding manifolds. The procedure has been previously applied to the analysis of the scarring effect of unstable POs in chaotic systems , as well as to the study of the quantum dynamics that takes place beyond the first recurrence of the corresponding motion . As it will be shown, these initial states can be used as an important tool to disentangle the complexity of the full quantum dynamics taking place in the area of influence of a given PO, making use of the associated classical mechanics as a guide in the process.

Initially, the method considers a wave function, $`|\varphi (0)`$, well localized in phase space in the vicinity of a given PO. The corresponding time evolution can be followed by means of the associated autocorrelation function,

$$C(t)=\varphi (0)|\varphi (t)=\varphi (0)|e^{i\widehat{H}t}|\varphi (0).$$
(1)
Recurrences in $`C(t)`$ are known to determine the low resolution structure of the corresponding spectra ,

$$I_T(E)=\frac{1}{2\pi }_{\mathrm{}}^{\mathrm{}}𝑑tW_T(t)C(t)e^{iEt},$$
(2)
and are the origin of the scarring effect, as first discussed in . $`W_T(t)`$ is a window function to take into account possible different resolutions in $`I_T(E)`$.

A widely spread choice for $`|\varphi (0)`$ is a minimun uncertainty harmonic oscillator coherent state ($`\mathrm{}`$ is set equal to one throughout this paper)

$`G_{𝐫^0,𝐏^0}(𝐫)`$ $`=`$ $`{\displaystyle \underset{j}{}}\left({\displaystyle \frac{1}{\pi \sigma _j^2}}\right)^{1/4}\mathrm{exp}\left[{\displaystyle \frac{(r_jr_j^0)^2}{2\sigma _j^2}}\right]`$ (4)
$`\times \mathrm{exp}\left[iP_j^0(r_jr_j^0)\right],`$
where ($`𝐫^0,𝐏^0`$) represents the coordinates and conjugate momenta of a suitable phase space point along the selected PO.

For bounded systems, the evolution of this packet can be followed quite conveniently by projection on a complete set of eigenfunctions, $`|n`$, of the Hamiltonian, $`\widehat{H}`$,

$$|\varphi (t)=e^{i\widehat{H}t}|\varphi (0)=\underset{n}{}|nn|\varphi (0)e^{iE_nt},$$
(5)
where $`\{E_n\}`$ is the corresponding set of eigenvalues. In our case the billiard eigenstates have been computed using the scaling method .

Regarding the smoothing function $`W_T(t)`$ one has many choices; i.e. “hat”, exponential or Gaussian functions, which render $`\mathrm{sin}x/x`$, Lorentzian or Gaussian–type line envelopes, respectively. In our case we use the Gaussian decay:

$$W_T(t)=\frac{e^{t^2/2T^2}}{T},$$
(6)
for which the spectrum takes the form

$$I_T(E)=\frac{1}{(2\pi )^{1/2}}\underset{n}{}|n|\varphi (0)|^2e^{T^2(EE_n)^2/2}$$
(7)
The center of the initial wave packet defined in Eq. (4) will follow for some time the classical PO on which it was launched, with the dispersion rate evolving according to the Lyapunov exponent of the orbit .

In a second step , a wave function highly scarred by the PO can be constructed very efficiently by considering the average of $`|\varphi (t)`$ along the dynamics generated by the PO during a given amount of time, usually taken of the order of the PO period. Let us remark the importance of this last requirement, since if the packet is allowed to evolve for a very long time, so that due to the chaotic character of the dynamics it has the opportunity to sample all the available phase space, just the (complicated) eigenfunctions of the system will be obtained. The corresponding expression for the localized wave function corresponding to an energy $`E_0`$ (usually taken as the center of one of the bands appearing in the low resolution spectrum generated by $`|\varphi (0)`$) reads as follows:

$`|\psi ^{E_0}`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑te^{iE_0t}W_T(t)|\varphi (t)`$ (8)
$`=`$ $`{\displaystyle \frac{1}{(2\pi )^{1/2}}}{\displaystyle \underset{n}{}}|nn|\varphi (0)e^{T^2(E_0E_n)^2/2}`$ (9)
Phase space representations for the wave functions described in this section can be constructed in a number of different ways. In this paper we will follow the procedure described in Ref., which relies on the use of normal derivatives of the wave functions evaluated at the boundary of the billiard. On this boundary, Birkhoff coordinates , ($`q,p`$), are used to define both classical and quantal Poincaré surfaces of section, such that $`q`$, is the arc length coordinate, and $`p=𝐩\widehat{𝐭}/|𝐩|`$ is the fraction of tangential momentum. The coherent states necessary to construct this representation (Husimi functions ) are then defined as:

$`G_{q^0p^0}(q)`$ $`=`$ $`\left({\displaystyle \frac{1}{\pi \sigma ^2}}\right)^{1/4}\mathrm{exp}\left[{\displaystyle \frac{1}{2\sigma ^2}}(qq^0)^2\right]`$ (11)
$`\times \mathrm{exp}[ip^0(qq^0)].`$
This expression corresponds to a boundary wave packet at point ($`q^0,p^0`$) in the surface of seccion; representing a bounce off a given boundary point with a specific tangential momentum. Then, for a wavefunction with normal derivative on the boundary $`\psi (q)`$ (extended periodically to the real line), the corresponding Husimi function is given by

$$H(q^0,p^0)=\left|G_{q^0p^0}|\psi \right|^2$$
(12)
## III Topological Characteristics of the Localized Wave Functions

The system that we have chosen to study is a particle of mass 1/2 moving in a desymmetrized stadium billiard with Newman boundary conditions on the symmetry axis (only even–even parity wave functions will then be considered). The radius $`r`$ is taken equal to unity and the enclosed area is 1+$`\pi `$/4. This system is known to exhibit a great degree of chaoticity, both from the classical and quantal points of view.

In this case, and due to the symmetry of the problem, symmetry adapted initial wave functions should be used; the corresponding expression being given by

$`x,y|\varphi (0)`$ $`=`$ $`G_{x^0,y^0,P_x^0,P_y^0}+G_{x^0,y^0,P_x^0,P_y^0}`$ (13)
$`+`$ $`G_{x^0,y^0,P_x^0,P_y^0}+G_{x^0,y^0,P_x^0,P_y^0}`$ (14)
$`+`$ $`c.c.`$ (15)
is obtained by imposing Newman boundary condition at the symmetry axis.

We consider the dynamics dominated by the horizontal PO, running along the $`x`$ axis with $`y=0`$. Accordingly, we use for our study a symmetry adapted Gaussian wave packet \[Eq. (15)\] with centers defined by the phase space point $`(x^0,y^0,P_x^0,P_y^0)=(1,0,k,0)`$ and width $`\sigma _x=\sigma _y=1.603/k^{1/2}`$, where $`k`$ is the usual wave number. The corresponding autocorrelation function for a packet with an energy at the center of $`E`$=3600, for which the period of the horizontal PO is equal to $`T_0`$=1/30, is shown in Fig. 1. As can be seen, two main peaks exist at $`t0.016`$ and 0.033 respectively, being the first one substantially lower than the second. After that, the correlation function becomes very complicated. The relative intensity of the two peaks can be easily explained if one takes into account that the initial coherent state that we are using is the sum of two packets, launched with opposite values of the momentum. The two packets first separate, one being scattered by the vertical axis and the other by the circular hard wall of the billiard, so that when they first return to the initial point, at a time roughly equal to $`T_0/2`$, they have picked up different phases. This gives rise to the first (small) recurrence peak in $`|C(t)|`$. Afterwards, the two components continue moving, colliding with the circular and vertical walls, respectively. In this way, when some time later, at $`tT_0`$, they return again close to the initial point, both arrive with approximately the same collisional phase, this giving rise to the second, more pronounced peak in the autocorrelation function.

Let us consider next the localization of the wave functions obtained with the method described in the previous section \[see Eq. (9)\] for different values of the smoothing time. We choose $`T`$=0.01, 0.03, 0.06, and 0.1, in order to cover a range of values smaller and greater than $`T_0`$. The results are shown in Fig. 2, where we present the corresponding squared wave functions, Husimi based quantum sufaces of sections, and the associated portions of the spectra generated from $`|\varphi (0)`$. Three comments are in order.

In the first place, the wave functions in configuration space (top tier) appear localized along the PO, presenting a number of nodes, $`m`$, in agrement with the value obtained from the Bohr–Sommerfeld quantization condition

$$k=\frac{2\pi }{L}\left(m+\frac{\nu }{4}\right)$$
(16)
when applied to the spectral peaks (bottom tier). For the horizontal PO, $`L=4`$ and the Maslov index $`\nu `$=3.

In the second place, and much more importantly, in the localized wave functions it is clearly observable the distorsion in the structure induced by the self focal point. A discussion of this effect can be found in Ref. . This is a very important point, and is a consequence of the fact that these functions have been defined in such a way–averaging over the motion along the orbit–that they are not only located in the right position of configuration space, but contain also information on the dynamics of the system in the vicinity of the PO. This is more clearly seen when examining the phase space picture of these functions, i.e. middle tier of Fig. 2. There, it is observed that the probability density do not just localize over the fixed point corresponding to the scarring orbit, but it is also spread significantly along their manifolds.

Finally, this structure over the manifolds changes as the value of $`T`$ is increased, giving the original packet the opportunity to explore more and more of the dynamics induced by the PO. Indeed, we see in Fig. 2 that for $`T`$=0.01, well before the first peak in $`|C(t)|`$, the quantum surface of section appear essentially localized on the fixed phase point corresponding to the horizontal PO, with just some probability sticking out of this region along the manifolds (also plotted in the figure). For $`T`$=0.03 this latter effect is more pronounced, and the pattern along the linearized part of the stable and unstable manifolds is fully developed. At the same time, the corresponding low resolution spectral envelopes $`I_T(E)`$ (continuous line in the bottom tier) define, with increasing precission, the clump of eigenstates contributing to this scar . Past this time, i.e. for $`T`$=0.06, well beyond the first (true) recurrence of the horizontal PO, the original packet has had the opportunity to enter into the truly nonlinear part of the periodic trajectory, exploring regions of phase space far apart from the fixed point, where the dynamical flux gets more complicated. However, we see that the averaged wave function is well confined on the manifolds, even following the sharp kinks that they present. This process continues for $`T`$=0.1, value at which the localized wave function appears covering practically all these invariant structures of phase space. In this case the resolution in the spectra is higher, and the corresponding spectral envelope define several bands, four in our case. These higher resolution bands  correspond to the interaction with other POs of longer periods, as will be discussed in the next section. The signatures of these longer POs are seen in the configuration space plots.

## IV Dynamics of the Localized Wave Functions

Another aspect of this problem which is worth considering is how the process described above is affected when the transition to the semiclassical limit is considered. For this purpose we have repeated the previous calculations, launching the packet from the same point as before but with different values of the energy. The results are shown in Fig. 3. In it we present the localized wave functions, Husimi based quantum Poincaré surfaces of section, and spectra for $`E_0`$=400 and $`T`$=0.1 \[column (a)\], and also for $`E_0`$=40267 and $`T`$=0.01 \[column (b)\]. These values of the smoothing times correspond to the periods of the horizontal PO at the selected energies, so that there are to be compared to those in Fig. 2 (b). As can be seen, the structure exhibited by these functions is essentially the same, if one takes into account the fact that all length features scale as $`k^{1/2}`$. This scaling is the origin of the better definition of the focalization effect existing in the configuration space wave functions at higher energies.

Using Gutzwiller trace formulae, similar smoothed wave function has recently been calculated , also showing that the invariant hyperbolic classical structure is contained in quantum mechanics. In this case no dynamical information was used, with the result that many more (of the order of 10000) states were needed in the summation in order to get the same kind of localization that we are showing in this paper.

In this section we investigate the dynamics of the localized wave functions described in the previous section, with the purpose of disentangling the complexity of the quantum dynamics in the problem that we are considering. As stated in the Introduction, our aim is to provide an explanation of the quantum dynamics of highly chaotic systems in terms of classical invariant structures.

In Fig. 4 we present the autocorrelation function calculated using the (non–stationary) localized wave functions \[Eq. (9)\] of Fig. 3 (a) and Fig. 2 (b). The results of part (a), corresponding to $`E_0`$=400 and $`T`$=0.1 , show a series of recurrences at time intervals approximately equal to 0.38. In each of them, practically all the initial probability is recovered, and the process shows no evidence of dispersion up to the times considered in our calculations.

This phenomenon has already been discussed in Ref. , where Gaussian wave packets “stretched” along the horizontal PO were used, and necessary conditions for the existence of revivals were given. As indicated before, in our case we use initial functions that are better adapted to the dynamics of the orbit, thus showing an improved revival behaviour. In this respect, it is important to note that, opossed to which it was found in , our functions show revivals, in this wavelength regime, for an ample range of energy values and also for other POs.

From the quantum mechanical point of view these recurrences are easy to interpret, since they are the result of the beating frequency originated from the two (almost exclusive) eigenstate contribution to the initial packet (see spectrum in the bottom part of Fig. 3). On the other hand, at the classical level, the phenomenon is more difficult to understand. To further investigate these recurrences from this point of view, we show in Fig. 5 a series of snapshots in phase space covering the period of time elapsed up to the first recurrence of $`|\psi ^{E_0}`$. Notice that due to the symmetry of the figure we have only plotted one half of the complete phase space of the billiard. To help interpreting this figure, the position of the fixed points and associated manifolds, corresponding to some short POs related to the horizontal one, have also been included. These orbits are presented in Fig. 6. For the sake of clarity, the labelling of the POs was not included in Fig. 5, but it can be found in Fig. 7.

In the first three images of Fig. 5 we see how the original probability density distribution elongates, stretching along the unstable manifolds of the horizontal PO and contracting along the stable ones. Afterwards, the process continues, populating the regions around the fixed points corresponding to orbits (3),(4) and (1),(2) , so that at $`t`$=0.127–0.169 a noticeably accumulation of probability around these points has taken place. Notice also the effect of the flux along the stable manifolds of these POs, which at $`t`$=0.127 is narrowing the elongated shape of the distribution along some specific lines, and at $`t`$=0.169 has split the packet in several pieces. At $`t`$=0.211 this process reaches its maximum, and the probability remaining in the initial region is practically zero. After this point, the evolution continues in a similar way with the probability returning gradually to the vicinity of the initial point, although obviously this time following a path going along the stable manifolds. This second half of the dynamics takes place following the same series of intermediate steps.

In Fig. 4 (b) the autocorrelation function for the initial state of Fig. 2 (b) is shown. In this case the recurrence pattern is much more complex, but it is nevertheless noticeable the existence of a fairly good number of recurrences in which the overlap with the initial distribution is greater than 50%. The difference with the example presented in the part (a) of the figure is clear. In this case, the initial wave function is formed by a substantially greater number of contributions \[see the associated spectrum at the bottom of Fig. 2 (b)\], and then the probability that all off them get in phase is much smaller.

The corresponding evolution in phase space is presented, up to the first big recurrence at $`t0.47`$, in Fig. 8. At this higher energy we see that the mechanism for the dispersion of the probability density is practically the same, although there exists certain differences worth commenting. In the first place, the time scale for the dispersion is obviously faster, i.e. compare the extent of the evolution in the first snapshot at $`t`$=0.059 with the third one at $`t`$=0.127 of Fig. 5. In the second place, the probability distribution appears in this case much more localized on the participating invariant classical structures for all values of time considered.

Finally, at certain values of $`t`$, the probability show peaks which localize on new POs, not important for the dynamical evolution at $`E_0`$=400 shown in Fig. 5. For example, at $`t=0.213`$ the probability density localize on orbits (5) and (6) of Fig. 6. This effect is not unexpected, since it is reasonable to admit that as we increase the energy, going towards the semiclassical limit, more POs are necessary to explain the quantum dynamics of our system, as discussed in Ref. .

## V Summary and Conclusions

The introduction of Gutzwiller’s trace formulae  in 1971 provided a method to understand the quantum mechanics of classically chaotic Hamiltonian systems in terms of POs. Its main drawback is the exponential growth in the number of orbits involved as the Heisenberg time, $`T_H`$, increases, and many efforts have been devoted to the development of resummation techniques that improve its convergence.

Our work considers an alternative approach, initiated with Refs. , based on the use of special functions well localized on a few POs and its manifolds, the number of which grows only linearly with $`T_H`$. To construct these localized functions we have used the method described in Ref. , based on the dynamics of a coherent wave packet averaged along a given PO and its neighborhood.

In this paper we have studied in detail some important characteristics and properties of these non–stationary scarred functions in the case of the horizontal PO in the desymmetrized Bunimovich stadium billiard, a paradigmatic example of chaos. In particular, we have considered how the localization properties, both in configuration and phase space, of these functions calculated at energies, $`E_0`$, close to the Bohr–Sommerfeld quantization conditions vary with the averaging time, $`T`$. Our results show that the probability density in phase space extends along the manifolds of the scarring PO, even when the associated dynamics have left the linear regime. In the second place, the transition to the semiclassical limit has also been investigated by simultaneously changing $`E_0`$ and $`T`$. As it takes place, all length features get better defined, for example the focal points corresponding to the hyperbolic structure around the fixed point. Finally, we have also studied the dynamics of these non–stationary wave functions. When the corresponding probability density in phase space is followed a dispersion along the manifolds is observed. In this process a building up of probability peaks takes place on a few POs related to the original one (the horizontal PO in our case). The recurrences in the corresponding autocorrelation function can be very important, and we have even shown one example in which they are particularly pronounced.

To conclude, we think that this work provides valuable information that can help in the disentanglement of the quantum mechanics of very chaotic Hamiltonian systems in terms of simple classical structures.

## Acknowledgments

This work was partially supported by PICT97 03–00050–01015, SECYT–ECOS (Argentina), and DGES (Spain) under contracts no. PB96–76, PB98–115, and BMF2000–437. DAW gratefully acknowledges support from CONICET (Argentina) and AECI (Spain).