text:
# Keck Mid-Infrared Imaging of the QSO 2237+0305

## 1. Introduction

QSO 2237+0305 (z=1.695, Huchra et al. 1985), referred to hereafter as “the Einstein Cross,” is the first gravitationally lensed QSO with confirmed variations due to lensing by the stars in the lens galaxy (z=0.0395), referred to as “microlensing” (Irwin et al. 1989). Optical variations of $`0.5`$ magnitudes due to microlensing are observed at least once per year (Racine 1992), which implies the optical emission region is smaller than an Einstein radius, $`R_\mathrm{E}=1.1\times 10^{17}M^{1/2}h^{1/2}\mathrm{cm}`$, where $`M`$ is the mass of the lensing star in solar masses and $`h=H_0/(100\mathrm{k}\mathrm{m}\mathrm{s}^1\mathrm{Mpc}^1)`$. The peak magnification and time scale of a high amplification event (caused by a caustic crossing the QSO) is strongly dependent on the size of the QSO. An event in 1988 has been used to limit the optical source size to about $`10^{15}`$ cm (Wambsganss, Schneider, & Paczyński 1990, Rauch & Blandford 1991, Wyithe et al. 2000).

Typically QSOs have a change in spectral slope near 1$`\mathrm{\mu m}`$ in the rest frame, which may imply that different physical components contribute longward and shortward of this wavelength (Neugebauer et al. 1987, Sanders et al. 1989, Barvainis 1993, Begelman 1994). One way to distinguish these components is reverberation mapping: observations of a few Seyfert galaxies indicate that the infrared emission region is more extended than the optical (Barvainis 1992, Nelson 1996). However, there exist no such direct constraints on the size of the infrared emission in QSOs; recent observations indicate this region may be rather compact, although the data are ambiguous (Neugebauer & Matthews 1999). Another technique for constraining the size of the emission region in QSOs is gravitational microlensing. Sources larger than an Einstein radius have smaller fluctuations than sources smaller than an Einstein radius (Kayser, Refsdal, & Stabell 1986). Despite its potential for probing the QSO central engine, microlensing has not yet been looked for in the mid-infrared ($`\lambda >3\mathrm{\mu m}`$). Currently the most popular model of the infrared (IR) spectrum in QSOs is thermal emission from hot dust. Evidence includes the lack of IR variability in most QSOs and the ubiquitous spectral dip around one micron which can be ascribed to dust sublimation at a temperature of $`T_s1800`$K. However, some radio quiet QSOs show phased IR/optical variability (Neugebauer & Matthews 1999), and dust absorption features are not generally seen in the mid-IR in bright AGN (Roche et al. 1991). Alternatively, the IR spectrum could be dominated by nonthermal synchrotron emission from the region near the black hole.

Microlensing variability in the mid-IR can place a limit on the size of the source, distinguishing between dust and synchrotron models. The minimum radius at which dust can exist near the central engine is

$$R_{\mathrm{dust}}=3\times 10^{18}\left(\frac{L}{2\times 10^{46}\frac{\mathrm{erg}}{\mathrm{s}}}\right)^{\frac{1}{2}}\left(\frac{T_s}{1800\mathrm{K}}\right)^2\left(\frac{ϵ}{0.1}\right)^{\frac{1}{2}}\mathrm{cm},$$
(1)
where $`L`$ is the QSO luminosity, $`ϵ`$ is the dust emissivity, and $`T_s`$ is the dust sublimation temperature. $`R_{\mathrm{dust}}`$ is larger than $`R_\mathrm{E}`$ so that if dust is responsible for the emission, microlensing fluctuations will be very small ($`0.05`$ magnitudes for $`R_{\mathrm{dust}}/R_\mathrm{E}=30`$, Refsdal and Stabell 1991).

On the other hand, if the infrared emission is due to synchrotron, then a turnover should be present at low frequencies due to self-absorption. Such turnovers are typically observed in quasars in the far infrared, but they are much steeper than the $`F_\nu \nu ^{5/2}`$ behavior expected from homogeneous self-absorbed synchrotron models (Hughes et al. 1993). Nevertheless, if such a turnover is due to self-absorption, then the requirement that the source be optically thick at the turnover frequency places a strong upper bound on the size of the emission region of

$$R_{\mathrm{synch}}10^{14}\left(\frac{F_\nu }{1\mathrm{m}\mathrm{J}\mathrm{y}}\right)^{1/2}\left(\frac{B}{1\mathrm{G}}\right)^{1/4}\left(\frac{\lambda }{100\mathrm{\mu m}}\right)^{5/4}\mathrm{cm},$$
(2)
where $`F_\nu `$ is the flux at the turnover wavelength $`\lambda `$, and $`B`$ is the magnetic field in the emission region. Unfortunately the Einstein Cross lies in a region of the sky which was not surveyed by IRAS, so the far infrared flux and turnover wavelength are unknown. The typical flux in the far-infrared for a QSO as bright as the Einstein Cross is $``$ 10 mJy, so the size of the synchrotron emission region is expected to be comparable to the optical size (where $`B`$ could be $`10^4`$ G). The limit on the size of a synchrotron emitting region is therefore much smaller than $`R_\mathrm{E}`$. If synchrotron emission is responsible for the IR, then microlensing variations in this region of the spectrum will be at least as large as in the optical.

A lack of microlensing might allow measurement of the amplification due to the gravitational lens galaxy. If the IR emission is due to dust, then it will not be microlensed, and, as emphasized by Kochanek (1991), would provide a direct measurement of the relative magnifications due to the lens galaxy rather than due to microlensing (although microlensing by clusters of stars can still affect a large emission region). This data is required for accurate modeling of this lens (Witt and Mao 1994). This has been attempted by using the CIII\] line, which originates from the extended broad line region and is therefore not microlensed (Yee & De Robertis 1992, Racine 1992, Fitte & Adam 1994, Saust 1994, Lewis et al. 1998); however, the CIII\] line is affected by extinction and continuum subtraction. To avoid these problems, the flux ratios were measured in the radio by Falco et al. (1996) with a signal-to-noise ratio (SNR) of only 2-4. One can do much better in the IR. Assuming a power-law mass distribution in the lensing galaxy ($`Mr^\beta `$), Wambsganss & Paczyński (1994) demonstrated that the total magnification of the lens could vary dramatically as $`\beta `$ varies. Their model predicts different flux ratios between the four images as $`\beta `$ varies, so measuring the flux ratios of an extended source might constrain the total amplification.

We present the first observations of the Einstein Cross near 10 microns ($``$3 microns rest frame). These observations were taken at Keck, which is the only telescope with a large enough collecting area to detect such a faint QSO. In addition, at 10 microns it is diffraction limited with about $`0\stackrel{}{\mathrm{.}}25`$ resolution, making it easy to distinguish between the four images which have a smallest separation of about 1$`\stackrel{}{\mathrm{.}}`$ In section 2, the observations are discussed. In section 3, the data reduction is discussed and the results are presented. In section 4 the implications for emission models and lens models are discussed. In section 5 the results are summarized and future directions are pointed out.

## 2. Observations

The Einstein Cross was observed with the Keck Long Wavelength Spectrometer (LWS) on three nights in 1999. The LWS is a 128$`\times `$ 128 Boeing Si:As array with a $`10\stackrel{}{\mathrm{.}}5`$ field of view at the Cassegrain focus of Keck I (Jones & Puetter 1993). The LWS was used in imaging mode with chop and nod which were necessary due to the faintness of the source and sky variability. The chop and nod were both fixed at $`10\mathrm{}`$, which is off of the LWS chip. The field was dithered by $`1\mathrm{}`$ in a 5-point pattern. The 8.9 and 11.7 $`\mathrm{\mu m}`$ filters were used, which have $`80`$% flux transmission from 8.4 to 9.2 $`\mathrm{\mu m}`$ and 11.2 to 12.2 $`\mathrm{\mu m}`$, respectively. Total observing times and on-source integration times are listed in Table 1. In addition to observing the QSO, the stars HR 5616, 4 Lac, and HR 8551 were observed (near the same airmass) at 11.7 microns, yielding 13,300 ADU/Jy/sec with a standard deviation of 5% for different objects and different nights. We observed stars HR 8551, $`\alpha `$ Boo, and HR 2443 at 8.9 microns, yielding 10,300 ADU/Jy/sec with a standard deviation of 6%.

Fig. 1: $`5\mathrm{}\times 5\mathrm{}`$ image of the Einstein Cross at 8.9 and 11.7 microns taken with LWS on July 28, 1999 and September 25 & 26 1999 (all data summed together; North is up). The images are labeled according to Yee (1988). bottom-right image.

## 3. Reduction

### 3.1. Basic Reduction

Successive chopped frames were subtracted and added to successive nod sets to get each frame. Frames were median clipped to get rid of cosmic rays and bad pixels. As the QSO is rather faint, it is not detected in any single frame, which typically consists of about five minutes of integration. To detect the QSO, the data were de-dithered and successive frames were added. Due to imperfect pointing, the de-dithering can cause some broadening of the point-spread-function (PSF). As the QSO is too faint to allow correction of pointing errors for individual frames, the stellar PSF was convolved with gaussians of various widths to see how this affects the results, as discussed below. The standard star PSF was constructed using an identical de-dithering procedure as that applied to the QSO images, but also was corrected for pointing errors. Scaled and broadened stellar PSFs were fitted and subtracted using the known relative QSO image positions (from Falco et al. 1996) while varying absolute position on the chip to minimize the residuals of the image minus the model using a linear least-squares fit. Using the relative positions measured by HST (Rix et al. 1992) changed the flux ratios by an amount much smaller than the error bars. The resulting fluxes are shown in Table 2. Figure 1 shows the summed image from both bands and both nights, smoothed with a 0.25$`\mathrm{}`$ boxcar to eliminate noise. The lens galaxy is undetected; neither is there any evidence of a fifth lensed QSO image.

For comparison to models, the fraction of the total Einstein Cross flux contained in each image was computed, hereafter referred to as the “flux ratio”; i.e. if all four images have the same flux, then the flux ratio is 1/4 for each image.

### 3.2. PSF

As mentioned above, the pointing errors resulting from dithering can lead to a broadening of the PSF for long observations on faint sources, such as those presented here. A stellar PSF (corrected for pointing errors) was convolved with a gaussian, varying the FWHM until the best fit was obtained for each image. The typical FWHM varies between 3 and 5 pixels (see Table 2), but the flux ratio of each image on each night in each waveband varies by only 5% if the FWHM is changed from 3 to 5 pixels ($`0\stackrel{}{\mathrm{.}}25`$ to $`0\stackrel{}{\mathrm{.}}4`$). This demonstrates that the relative flux ratios are rather insensitive to the gaussian FWHM. The absolute fluxes, however, are quite sensitive to the assumed FWHM. Fortunately, the argument for microlensing depends only on the relative flux ratios.

### 3.3. Error bars

A Monte Carlo technique was used to compute the error bars on the fluxes. During these observations the LWS chip exhibited pattern noise which prevented a direct measurement of the error bars from photon counting statistics. Instead, the best-fit model of the four images was added back in at various points on the chip to create simulated images. The simulated images were run through the entire reduction procedure, measuring the position, gaussian FWHM, and fluxes of the four images. This was repeated $``$100 times to compute the standard deviation of the image fluxes, flux ratios, FWHM, and position. The error bars on the flux ratios are typically about a factor of 2 greater than what one would infer from shot-noise statistics only. Thus, only the Monte Carlo error bars are reported in this paper (Table 2; 1-sigma error bars are always quoted). The error bars on the absolute flux are quite large, mostly due to the uncertainty in the point spread function. The error bars on the measured FWHM are reported in Table 2; this is the FWHM gaussian which was convolved with the stellar PSF to approximate the QSO PSF.

### 3.4. Consistency of the data

The flux ratios measured on both nights and in both bands are consistent with being constant with time and with wavelength, as shown in Figure 2. Any disagreement between the measurements are consistent with random errors causing the differences; there is no need to invoke variability in the flux ratios with time or between the two bands, nor any need to find systematic errors to explain the differences.

Fig. 2: Fraction of the total flux of the Einstein Cross contained in each image. The data from September are plotted in bold; the 11.7 $`\mathrm{\mu m}`$ data with asterisks; and the 8.9 $`\mathrm{\mu m}`$ data with squares. The error bars reflect the statistical uncertainties

## 4. Discussion

### 4.1. Microlensing Constraints

Table 3 compares the mid-infrared flux ratios (averaged from both nights and both bands) to the flux ratios in other wave bands. The V band magnitudes are from the OGLE monitoring data on the dates closest to the Keck observations (available on their website; Woźniak et al. 2000); this data was corrected for extinction as described in the Appendix. The VLA radio data are from Falco et al. (1996) which were taken in 1995. The CIII\] ratios are from Racine (1992), data taken with CFHT in 1991. One can immediately draw the conclusion that the infrared is extended (and thus probably dust emission) from the following lines of evidence:

1) The flux ratios differ between the IR & V band by more than 6$`\sigma `$ for each image. Simultaneous with the Keck observations were strong optical microlensing events in images A and C at the time of the infrared observations (OGLE website; Woźniak et al. 2000). Since the optical is known to be microlensed, this indicates that the IR is extended on a scale comparable to or larger than an Einstein radius, causing dramatic differences in the flux ratios.

2) The image ratios in the IR are close to the ratios predicted by some lens models, lending weight to the idea that the IR is more extended and thus not affected by microlensing. Table 3 shows the predictions of the $`\lambda =1`$ model from Schmidt, Webster, & Lewis (1998, hereafter SWL), which were derived without using the flux ratios as a model constraint. For all published models $`\chi ^2=_{i=\{A,B,C,D\}}(f_i(obs)f_i(mod))^2/\sigma _i^2`$ was computed, where $`f_i`$ is fractional flux of each image and $`\sigma _i`$ is the measured error bar. The three best models are $`\lambda =1`$ of SWL with $`\chi ^2=7`$, model 4 of Kochanek (1991) with $`\chi ^2=7`$, and model 1 of Witt, Mao, & Schechter (1995) with $`\chi ^2=8`$. These $`\chi ^2`$ indicate that there may be deviations of the infrared magnifications from the predicted magnifications at the level of 10%. This may be due to inaccuracies in the modeling, an underestimation of the error bars, or microlensing fluctuations at the 10% level, which is to be expected for a source 15 times the size of an Einstein radius (Refsdal & Stabell 1991, assuming $`\sigma =0.5`$). All other published gravitational lens models for Einstein Cross have larger $`\chi ^2`$ (Kent & Falco 1988; Kochanek 1991; Rix, Schneider, & Bahcall 1992; Wambsganss & Paczyński 1994; Witt et al. 1995; Chae, Turnshek, & Khersonsky 1998).

3) The image ratios at 8.9 and 11.7 microns are consistent within the error bars, indicating that neither wavelength is affected by either extinction or differing microlensing amplification due to different sizes as a function of wavelength.

4) The radio emission is thought to be extended on a scale larger than an Einstein radius (possibly in a jet) and thus is unaffected by microlensing. Since the IR image ratios are consistent with the radio image ratios within the error bars, this implies that the IR may also be unaffected by microlensing. Note that the infrared SNR is much better than in the radio.

5) In the case of images C and D, the IR image ratios are not consistent with the CIII\] image ratios, but agree better than with the V band flux ratios. The deviations may be due to microlensing fluctuations of CIII\] if it is only a few times larger than an Einstein radius, or due to the systematic errors which affect the CIII\] measurements.

6) Image C dropped in optical flux by a factor of 1.4 between August 4 and September 26 (the two OGLE observations closest in time to the two Keck observations). Our observations are consistent with no variation in image C in either band between the two dates, although the possibility of a variation is not strongly ruled out.

7) Consider a model in which a point source and an extended source contribute to the infrared flux (which would be the case if both a dusty torus and synchrotron source contributed to the infrared flux). Then, the synchrotron source would be magnified by the same amount as the optical (due to microlensing), while the extended source would be magnified by the macrolensing magnification (assuming the SWL model holds). With these assumptions one can compute what the fraction of point source flux is relative to the extended source flux in order to give the infrared flux ratios. Carrying out this exercise, one finds that the point source contribution must be negative to explain the difference between images A and B or A and C. Simply put, the infrared flux ratio of B is larger than the model, while the point source (i.e. optical) is smaller than the model (see Table 3), so adding a point source to an extended source can only decrease the flux with respect to the model. Since the point source flux must be non-negative, it is most likely zero. Then, the disagreement between the infrared flux ratios and the model flux ratios is probably due to inaccuracies in the lens model, Poisson fluctuations due to a random number of lenses passing in front of a given image (since the source is not truly infinite), and/or errors in the measured fluxes.

Given these lines of evidence, one concludes that the mid-infrared emission region is larger than the Einstein radius of a typical star in the lens galaxy, about $`10^{17}`$ cm. The infrared may be extended by $`0\stackrel{}{\mathrm{.}}24`$ which corresponds to about 0.8 kpc in the source plane; however, this extension is possibly due to pointing errors as discussed above. Consequently, $`0.03\mathrm{pc}R_{IR}800\mathrm{p}\mathrm{c}`$. Comparing the optical and infrared data with simulated microlensing lightcurves with various source sizes, one can compute the probability distribution of the infrared source size, which will be done in future work (Wyithe & Agol, in preparation).

Fig. 3: Spectral energy distribution of the Einstein Cross (non-simultaneous; frequency in observed frame). The near-infrared and optical points are taken from Nadeau et al. (1991); the ultraviolet from Blanton, Turner, & Wambsganss (1998); the X-ray data (dashed line) is from Wambsganss et al. (1999); the radio points are from Falco et al. (1996); the mm points (diamonds) are 3 $`\sigma `$ upper limits from Richard Barvainis (priv. comm.); and the mid-infrared data are from this work. The solid line is the radio-quiet QSO composite from Elvis et al. (1994) scaled to match the optical flux. The dotted line is the dust shell model described in the text based on Barvainis (1987).

### 4.2. Spectral Energy Distribution

To study in more detail the properties of the Einstein Cross, an SED was constructed of all available data from the radio to X-ray from VLA, IRAM, CFHT, Keck, HST, and ROSAT. The lens galaxy has been removed (or does not contribute) for all the observations shown. The near-IR, optical, UV, and X-ray were corrected for extinction in the lens galaxy and the Milky Way as described in the Appendix. This brings the optical/UV slope into rough agreement with the average QSO slope. Figure 3 shows the SED in the observed frame with the extinction corrections, but without a correction for magnification. Figure 3 also shows the Elvis et al. (1994) radio-quiet QSO composite redshifted and scaled to the R band flux. Assuming a magnification of 16 (from the SWL model - this value is used throughout the rest of this paper), this implies a QSO luminosity of $`2\times 10^{46}`$ erg/s for the scaled composite. However, there appear to be deviations of the Einstein Cross from the composite QSO in the mm, infrared, UV, and X-ray. The mm upper limits are slightly lower than the average SED. The infrared has a narrower peak at higher frequency than the average QSO. A similar peak exists in 3C 273; however, 3C 273 peaks at 3 $`\mathrm{\mu m}`$ and is broader than the Einstein Cross (Barvainis 1987). The ultraviolet points are high, which may mean that this is a hot QSO or that the QSO has varied (or been magnified by microlensing) as the UV points were taken nearly a decade after the optical/near-IR points. The X-rays fall a factor of 3 below the average QSO, which is within the large scatter of optical to X-ray flux ratios of QSOs. The radio to optical flux ratio is identical to the QSO composite. The infrared discrepancy is of greatest relevance to this paper, which is discussed in turn.

### 4.3. Origin of Infrared Bump

The Keck mid-IR fluxes (summed over all four images) give a spectral index of $`\alpha =0.5\pm 0.2`$ ($`F_\nu \nu ^\alpha `$), which is quite different from the average for PG QSOs, $`\alpha =1.4\pm 0.5`$ (Neugebauer et al. 1987). This unusual behavior may be due to calibration problems. However, we think our error bar estimates are accurate, and we are unable to identify any calibration errors which affect the Keck observations at the 50% level which would be required to bring the slope into agreement with that of an average QSO. The data were compared with ISO CAM observations (TDT 54500601; P.I. Sangeeta Malhotra) of the Einstein Cross from the ISO archive which show a similar infrared slope of $`\alpha =0.45`$. The ISO CAM fluxes at 9.1 and 11.7 microns agree with the LWS measured fluxes within the 1$`\sigma `$ error bars. The near-IR spectral index is dramatically different: Nadeau et al. (1991) report $`\alpha =3.7`$ between the K-band and L-band (3.3 micron band). This implies that there is a peak in $`\nu F_\nu `$ between 3.3 microns and 8.9 microns (1.2 microns and 3.3 microns in the rest frame); i.e. the peak infrared flux occurs somewhere near 2 microns. Typically, QSOs peak much further in the infrared, near 10 microns (see Figure 3).

This unusual infrared bump might be explained as a quasi-blackbody component. There are several candidate physical components which might contribute near these wavelengths: synchrotron, hot dust, a starburst, or the host galaxy. Synchrotron is ruled out by microlensing, so the other possibilities are considered next.

#### 4.3.1 Hot Dust

A blackbody component was fit at the dust sublimation temperature ($`2000`$ K), which agrees with the data remarkably well if the effective radius (area = $`\pi r_d^2`$) of the source is $`r_d=2\times 10^{18}`$ cm. Treating the dust emission as a blackbody assumes the emissivity is unity, appropriate for large grains; small grains would reduce the emissivity, requiring a larger area, but lowering the temperature to 1600 K (for $`ϵ=2\pi r/\lambda `$, where $`r1\mathrm{\mu m}`$ is the dust grain radius). The inferred source size agrees with the dust sublimation estimate and with the microlensing size constraint. However, this implies that the dust should be located at a constant radius so that when it is exposed to the QSO continuum it equilibrates at a constant temperature. This means that the dust is in a shell geometry; if there is dust at a range of radii, then it will radiate at lower temperatures than the sublimation temperature which will make the infrared SED wider and flatter than a blackbody, typically rising in flux towards longer wavelengths (Pier & Krolik 1992).

To quantify the constraint on geometry, the near-infrared bump was fit with the spectral model of Barvainis (1987). Barvainis’ model assumes that the dust is in a shell subtending $`\mathrm{\Omega }`$ steradians of the central source with an inner radius ($`r_1`$) set by the luminosity of the source ($`L`$) and the dust sublimation temperature ($`T_{\mathrm{max}}`$). The shell is assumed to be optically thin in the near-IR ($`\lambda 1\mathrm{\mu m}`$), with a UV optical depth of 3. Barvainis’ uniform dust model was adopted with $`a=0.05\mathrm{\mu m}`$, emissivity proportional to $`\nu ^{1.7}`$, and we additionally assume that the dust grains have a constant number density, that $`L=2\times 10^{46}`$ erg/s, and that $`T_{\mathrm{max}}=1800`$K, fixing the inner radius at 1.1 pc, consistent with the microlensing size constraint. The two free parameters left are $`r_2/r_1`$ (the ratio of the outer to inner radius of the shell) and $`\mathrm{\Omega }`$. An adequate fit was found for $`r_2/r_1=3`$ and $`\mathrm{\Omega }=7\pi /4`$. In this model, the dust outside $`r_2`$ is shielded, and so does not contribute to the reprocessed light. The grain number density at the inner edge is $`5\times 10^9`$ cm<sup>-3</sup>, so the implied dust mass is about $`4\mathrm{M}_{}`$ for a dust density of 1 g/cm<sup>3</sup>. The Toomre Q-parameter is quite small, $`10^9`$ (for a mass of $`2\times 10^9M_{}`$), which indicates that this region should be gravitationally unstable. A clumpy disk model should have a similar spectrum since the smaller filling factor can be compensated by an increase in dust number density. The conclusion is that the narrow infrared peak requires a narrow range of radii (factor of 3) so that reprocessed luminosity is absorbed by a region with a narrow range of dust temperature. This narrow range is not unprecedented as Barvainis (1992) found a narrow range of radii for the dust emission in Fairall 9 from reverberation $`r_10.3`$pc and $`r_21.3`$pc.

This, of course, is a rather simplistic model. X-rays may penetrate deeper, heating dust at larger radii. Transient heating by X-rays causes temperature flickering which broadens the spectrum; this is only important for dust grains with radii $`a0.01\mathrm{\mu m}`$ (Voit 1991). Infrared radiation transfer has been ignored since the shell is optically thin by construction; reality may not be so simple. A warped disk with a warp angle which scales logarithmically radius would produce a much wider infrared bump than observed.

#### 4.3.2 Host Galaxy or Starburst

Normal galaxies can be well described by a 3500 K blackbody in the near infrared spectral region (Schmitt et al. 1997), which peaks at a wavelength of about $`1\mathrm{\mu m}`$ (rest), while the Einstein Cross spectrum peaks closer to $`2\mathrm{\mu m}`$ (rest). Thus, to fit the infrared bump requires strong extinction as in a starburst galaxy: an extinction of $`A_V7`$ is needed to fit the infrared bump; most of the light is absorbed by dust. The absorbed luminosity is $`L_{\mathrm{abs}}=6\times 10^{46}`$ erg/s, which must be re-radiated in the far-infrared (FIR). The radio-FIR relation (Condon 1992) limits the FIR luminosity to $`6\times 10^{45}`$ erg/s, a factor of 10 times smaller than the absorbed luminosity. A test of the host galaxy hypothesis would be to look for an extended source in the IR images. Meurer et al. (1997) have found an empirical limit on the bolometric surface brightness of starburst galaxies of $`2\times 10^{11}\mathrm{L}_{}/\mathrm{kpc}^2`$, which would indicate that this starburst has a size of $``$ 3.5 kpc, or $`1\mathrm{}`$ (including the factor of $`2`$ stretch due to lensing). This is larger than the infrared upper lmit of $`0\stackrel{}{\mathrm{.}}25`$, ruling out a starburst as the origin of the infrared bump.

One can roughly estimate the host galaxy luminosity from the QSO luminosity as follows. Rauch & Blandford (1991) have argued that sub-Eddington accretion disk models are too large to be consistent with the optical size constraint from microlensing. We confirmed this by fitting the optical/UV SED with various thin accretion disk models from Hubeny et al. (2000); all disks are a factor of a few too big to be consistent with the microlensing constraint. A smaller emission region might be achieved by decreasing the mass of the black hole; however, once the accretion rate reaches Eddington, the luminosity will be roughly fixed at the Eddington rate. Thus, if one assumes that the observed (isotropic) luminosity equals the Eddington rate, then the black hole mass is about $`2\times 10^9\mathrm{M}_{}`$. Using the Magorrian et al. (1998) relation between black hole mass and bulge luminosity, the bulge luminosity is estimated to be about $`6\times 10^9\mathrm{L}_{}`$, with an uncertainty of about a factor of 10. This is more than 150 times smaller than the near-IR luminosity, yet another indication that the infrared bump may not be due to a quiescent host galaxy.

## 5. Conclusions

The observations presented here indicate the mid-infrared source size in the the Einstein Cross must be greater than $`10^{17}`$cm, ruling out synchrotron emission. If true for other radio-quiet QSOs, this implies that fits to the Big Blue Bump which include a synchrotron power law may be flawed (Sun & Malkan 1989, Laor 1990). Starburst and host galaxy contributions to the mid-infrared emission are ruled out based on the low radio flux, small size, and large infrared luminosity, leaving hot dust as the only viable model for the infrared emission in the Einstein Cross. To avoid sublimation, the dust must be at a distance of about 1 pc from the QSO; to explain the narrowness of the infrared bump it must be located within $``$ 3pc (according to a simple reprocessing model); and to explain its strength, the dust must intercept about half of the QSO continuum. Since the mid-IR is extended, the flux ratios should be nearly equal to those of the galaxy gravitational lens model, so future models should be constructed with the infrared flux ratios as model constraints. Improved microlensing statistics, magnifications, and time delays will result from better lens models. Accurate monitoring of the Einstein Cross at mid-infrared wavelengths with, say, NGST, may reveal flux variations at the few % level, which would allow a better estimate of the size of the dust emission region. The Einstein Cross mid-infrared spectrum is quite different from any other QSO; it should be observed with SIRTF to better constrain the dust emission model.

Data presented herein were obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. This work was supported by NSF grants AST-9529230 and AST-9970827 and California Space Institute grant CS-16-93. We would like to thank Greg Wirth and Randy Campbell for their expert assistance during our observing runs. We also thank Robert Antonucci for inspiring this project. We thank Robert Antonucci, Tim Heckman, Julian Krolik, Sangeeta Malhotra, Nancy Levenson, Henrique Schmitt, Mark Voit, and Rachel Webster for useful discussions. We thank Richard Barvainis for sharing his millimeter upper limits before publication. We thank Przemyslaw Woźniak and the OGLE team for making their data available on the web before publication. We have also used archival data from ISO, an ESA project with instruments funded by ESA Member States (especially the PI countries: France, Germany, the Netherlands and the United Kingdom) and with the participation of ISAS and NASA.

## Appendix A Appendix: Correction for Extinction

To correct for extinction in the lens galaxy, the correlation between $`gi`$ color (Gunn bands) and lens galaxy surface brightness found by Racine (1991) was used. Extrapolation to zero surface brightness then gives the color corrected for extinction in the lens galaxy, if one assumes that the dust column is proportional to surface brightness (this assumption may be invalid if the dust is patchy). To solve for the absolute extinction, we assume the extinction curve is described by $`R=3.1\pm 0.9`$, the average Milky Way extinction curve. $`A_V`$ for image A is 0.88$`\pm `$0.21, B is 0.84$`\pm `$0.20, C is 1.30$`\pm `$0.31, and D is 1.15$`\pm `$0.27. Using the Milky Way extinction $`E(BV)=0.071`$ measured by Schlegel, Finkbeiner, & Davis (1998), the Nadeau et al. spectrum was corrected for each image for galactic extinction using the $`R=3.1`$ extinction curve (Fitzpatrick 1998). It is reassuring that the corrected spectrum is similar to that of an average QSO (see Figure 3). Yee (1988) and Falco et al. (1999) correct for extinction assuming the Einstein Cross has a spectral index equal to the average value for QSOs; since individual QSOs do not behave as the average, we prefer the Racine technique.

To correct the X-ray data for absorption, we adopt the dust-to-gas ratio given by Spitzer (1978): $`N_H=5.9\times 10^{21}E_{BV}\mathrm{mag}^1\mathrm{cm}^2`$. An average extinction of 1 for the lens galaxy and 0.2 for the Milky Way implies $`N_H=2.5\times 10^{21}\mathrm{cm}^2`$. The X-ray power law was assumed to be equal to that of an average QSO and the count rate of Wambsganss et al. (1999) was used to compute the unabsorbed X-ray flux in the ROSAT band.


format: text

text:
# Observation of resonance trapping in an open microwave cavity

\[

## Abstract

The coupling of a quantum mechanical system to open decay channels has been theoretically studied in numerous works, mainly in the context of nuclear physics but also in atomic, molecular and mesoscopic physics. Theory predicts that with increasing coupling strength to the channels the resonance widths of all states should first increase but finally decrease again for most of the states. In this letter, the first direct experimental verification of this effect, known as resonance trapping, is presented. In the experiment a microwave Sinai cavity with an attached waveguide with variable slit width was used.

PACS numbers: 05.45.-a, 03.65.Nk, 84.40.Az, 85.30.Vw

\]

Since more than ten years, interference phenomena in open quantum systems have been studied theoretically in the framework of different models. Common to all these studies is the appearance of different time scales as soon as the resonance states start to overlap see  and the recent papers  with references therein). Some of the states align with the decay channels and become short-lived while the remaining ones decouple to a great deal from the continuum and become long-lived (trapped). Due to this phenomenon, the number of relevant states will, in the short-time scale, be reduced while the system as a whole becomes dynamically stabilized. The phenomenologically introduced doorway states in nuclear physics provide an example for the alignment of the short-lived states with the channels . Calculations for microwave resonators showed that the trapped resonance states can be identified in the time-delay function and that short-lived collective modes are formed at large openings of the resonator . Resonance narrowing is inherent also in the Fano formalism . Similar effects have been found in the linewidths in a semiconductor microcavity with variable strength of normal-mode coupling . In spite of the many theoretical studies, the effect of resonance trapping has not yet been verified unambigously in an experiment. A theoretical study of neutron resonances in nuclei as a function of the interaction of a doorway state with narrow resonances  allowed only to draw the conclusion that resonance trapping is not in contradiction with experimental data. For a clear experimental demonstration of the trapping effect, the coupling strength to the decay channels should be tunable, which was not possible in all above mentioned experiments.

The mechanism of resonance trapping can be illustrated best on the basis of a schematical model. In an open quantum system the resonance states are allowed to decay, i. e. their energies are complex, $`_R=E_R\frac{i}{2}\mathrm{\Gamma }_R`$. The Hamilton operator is non-hermitian,

$$=_0i\alpha VV^{}.$$
(1)
Here $`_0`$ describes the $`N`$ discrete states of the closed quantum system coupled to $`K`$ decay channels by the $`N\times K`$ matrix $`V`$. $`_0`$ and $`VV^{}`$ are hermitian and $`\alpha `$ is a real parameter for the total coupling strength between the closed system and the channels. The complex eigenvalues of $``$ give the energy positions $`E_R`$ and widths $`\mathrm{\Gamma }_R`$ of the resonance states. Studies on the basis of this model were performed by different groups with different assumptions on the properties of $`_0`$ and $`V`$, see e.g. . The Hamiltonian $`_0`$ can be given by, e.g., a Gaussian orthogonal ensemble or a Poissonian othogonal ensemble and the elements of $`V`$ may follow from different assumptions on their statistical distribution. In all cases, the results show clearly the formation of different time scales due to the anti-hermitian part of $``$. At a certain critical value of the coupling parameter $`\alpha `$, the widths of $`NK`$ states start to decrease with increasing $`\alpha `$ and approach zero with $`\alpha \mathrm{}`$ while the widths of $`K`$ states increase.

The widths of the states of realistic systems show a more complicated behaviour than described in the statistical model. The widths of the long-lived states of molecules saturate with increasing coupling strength to the continuum, but do not approach zero . A saturation of the widths of the trapped states occurs, however, also in the calculations with the schematical model when it is improved by considering different coupling strengths for the different decay channels . This improvement, as well as the introduction of a complex coupling parameter instead of the real $`\alpha `$, are justified since they follow from formally rewriting the Schrödinger equation in the function space of both discrete and scattering states . Similar results are obtained for nuclei  and for open microwave cavities . On the basis of this result it is possible to clarify another problem discussed for resonance states in molecules at high level density : The fundamental quantum mechanical relation between the average width of the states and the average lifetime is violated when all resonance states are considered in the averaging procedure. The mechanism of resonance trapping makes, however, the averaging procedure meaningful only for either the long-lived or the short-lived states. Performing the average over the trapped states only, the relation between the average width and the average lifetime is recovered.

While $`\mathrm{\Gamma }_R0`$ with $`\alpha \mathrm{}`$ is a necessary condition for trapped states in the schematical model with the Hamiltonian (1), the calculations for realistic systems show that the widths of the trapped states may decrease, stop to increase and even slowly increase with increasing coupling strength to the continuum. Only the ratio between the widths of the trapped states and the widths of the short-lived states decreases with increasing coupling strength according to these studies. This behaviour may lead to such an interesting effect as population trapping resulting from the coupling of two states of an atom by means of a strong laser field . The interplay between the direct interaction of two states and the real and imaginary parts of their interaction via the continuum causes unexpected and sometimes contra-intuitive effects. In any case, the local processes between individual resonances are of central importance for the properties of the system.

Let us sketch the main ingredients of the trapping mechanism on the basis of calculations for two neighboured resonances. For $`N=2`$, $`K=1`$, and $`\alpha `$ replaced by $`\alpha e^{i\beta }`$, the Hamiltonian (1) reads

$$H_2=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)2i\alpha e^{i\beta }\left(\begin{array}{cc}\mathrm{cos}^2\phi & \mathrm{cos}\phi \mathrm{sin}\phi \\ \mathrm{cos}\phi \mathrm{sin}\phi & \mathrm{sin}^2\phi \end{array}\right)$$
(2)
where $`\beta `$ describes the ratio between the real and imaginary parts of the coupling matrix, and $`\phi `$ gives the relative coupling strength of the two states. In Fig. 1, the widths $`\mathrm{\Gamma }_R/2`$ and eigenvalues of the two resonances are shown as a function of $`\alpha `$ for $`\beta =\pi /18`$ and $`\phi =\pi /5`$. For small $`\alpha `$, the widths of both states increase with $`\alpha `$. Thereafter the two states attract each other in energy and their widths bifurcate: The width of one state starts to decrease with increasing $`\alpha `$ while the width of the other one increases more strongly. At still larger $`\alpha `$, the broad resonance gets shifted towards lower energies due to $`\beta 0`$.

The goal of this paper is to present an experimental verification of the effect of resonance trapping in a microwave cavity by considering the local interactions between individual resonances. The coupling of the discrete states of the cavity to an attached waveguide makes the system open.

The measuring technique is described elsewhere . Here we note only that for flat cavities the electromagnetic spectrum is equivalent to the quantum mechanical spectrum of the corresponding system, as long as one does not surpass the frequency $`\nu _{\mathrm{max}}=c/2h`$, where $`h`$ is the height of the cavity. The quantum mechanical energy $`E`$ corresponds to the square of the wavenumber $`k`$. For the measurement we used a quartered Sinai billiard with an attached waveguide (see Fig. 2). The actual form of the cavity is of little importance since resonance trapping is expected to take place in chaotic as well as in regular billiards . The frequency range of the waveguide is 8.2 to 12.5 GHz. This corresponds to energies between $`E=2.95\mathrm{cm}^2`$ and $`6.85\mathrm{cm}^2`$. Between the two thresholds at $`E=1.83\mathrm{cm}^2`$ and $`7.33\mathrm{cm}^2`$ only one mode can propagate through the waveguide. The microwaves are coupled into the system through an antenna at the end of the waveguide. They enter the billiard through a slit, the opening of which can be varied. The coupling matrix elements $`V`$, Eq. (1), are related to the width of the opening. The exact expression for this relation is not known for our situation with a limited number of resonances. We performed 102 measurements for different openings $`d`$ from $`3`$ to $`23.2`$ mm in steps of 0.2 mm. For all openings we measured the relative amplitude and relative phase of the reflected wave as a function of $`E`$ and expressed the values through the complex reflection coefficient $`R(E)`$. Due to the finite length of the waveguide with reflecting end there are also broad channel resonances in the measured $`R(E)`$. This, together with the fact that a part of the incoming flux gets absorbed in the walls of the resonator, complicates the data analysis somewhat.

In Fig. 3.a, a part of the measured spectrum $`|R|`$ for $`d=14`$ mm is shown. Due to the wall absorption we have $`|R|1`$ and the resonances show up as dips in $`|R|`$. Around $`E=5.75\mathrm{cm}^2`$ a pair of closely lying resonances can be seen. It will be evident later (see Fig. 4) that resonance trapping takes place between these two resonances.

The method used in this letter to obtain the energy positions $`E_R`$ and widths $`\mathrm{\Gamma }_R`$ of the resonances is as follows: In Fig. 3.b we plot the real and imaginary part of $`R(E)`$ for $`d=14`$ mm in the region $`5.54\mathrm{cm}^2E5.77\mathrm{cm}^2`$. This plot is known as the Argand diagram . In the figure, the narrow resonances show up as small circles superimposed on larger circles caused by the broader structures. In order to analyse the Argand diagram, we propose a centered time-delay analysis (CTDA). From the

measured points in a small region around an energy $`E`$, we define a local circle segment. We calculate the center $`C(E)`$ of this segment and define the angle $`\theta ^c(E)`$ as the angle of the complex number $`R(E)C(E)`$. The point $`C`$ moves with $`E`$ in such a way that $`\theta ^c(E)`$ increases by $`2\pi `$ in the neighbourhood of each resonance, even if the corresponding circle does not go around the origin. We define the corresponding ”time delay” as $`\tau ^c(E)=d\theta ^c(E)/dE`$. For overlapping resonances, $`\tau ^c`$ is just the sum of the contributions $`\tau _w^R`$ from the individual resonances. In the approximation of Breit-Wigner shaped resonances we have,

$$\tau _w^R(E)=\mathrm{\Gamma }_R/\left((EE_R)^2+\mathrm{\Gamma }_R^2/4\right).$$
(3)
In Fig. 3.c, the $`\tau ^c`$ obtained directly from the measurement are shown for $`d=14`$ mm (full line). Further, we show the $`\tau ^c`$ (dashed line) obtained from a fit to $`_R\tau _w^R(E)`$, Eq. (3). The measured data is split into smaller energy intervals by searching for sufficiently small minima in $`\tau ^c`$. The fitted curve agrees well with the measured values.

The usual fit of the complex $`R(E)`$ to a sum of Lorentzians was not possible because of the broad structures in $`R(E)`$ caused by the channel resonances. By means of the CTDA, however, any broad structures are automatically removed. It was also possible to do an automatic evaluation of the measured data, which was mandatory in view of the large amount of data. To test the CTDA method we performed an analysis of spectra created theoretically from $`H^{\mathrm{eff}}`$, Eq. (1), with the wall absorption simulated by adding a constant to the diagonal of the imaginary part of $`H^{\mathrm{eff}}`$. The comparison between the theoretical values and those extracted from the $`R(E)`$ showed a good agreement. A more detailed study of the CTDA method together with a comparison to possible other methods, e.g. the filter diagonalization method , will be published elsewhere.

The $`E_R`$ and $`\mathrm{\Gamma }_R/2`$ for four resonances obtained from the measured data with the CTDA method are presented in Fig. 4. The motion of the $`E_R`$ and $`\mathrm{\Gamma }_R/2`$ with $`d`$ is clearly observable. At least three states start to decrease in width with increasing $`d`$, i. e. there are some evident cases of resonance trapping. Note that for each opening $`d`$ of the slit the data have been fitted independently. The smoothness of the curves is thus a measure for the reliability by which the $`E_R`$ and $`\mathrm{\Gamma }_R`$ are extracted from the measured data. Fig. 4 even shows hierarchical trapping: one state first traps its closest neighbour only to be trapped at stronger coupling by another state.

We followed the motion of all the 127 resonance poles $`E_Ri/2\mathrm{\Gamma }_R`$ in the complex plane as a function of increasing $`d`$. (Due to the small number of resonances, no comparison to statistical predictions is made.) The widths of 47 resonances decrease with $`d`$ while those of 11 resonances stop to increase (within the accuracy of the CTDA). These two groups can unambigously be identified as trapped resonances, i. e. at least $`46\%`$ of the resonances get trapped. For 16 of the resonances increasing further in width, it was possible to identify other resonances getting trapped by them. By this they aquired an extra contribution to their widths. In total, we have found that at least $`58\%`$ of the resonances are strongly affected by trapping. These resonances are in the whole energy interval considered.

The estimation $`58\%`$ of the resonances getting affected by resonance trapping is the lowest limit. According to theoretical studies of the diagonal elements of the effective Hamilton operator Eq. (1) (i.e. without taking resonance trapping into account) and studies for isolated resonances , one would expect a more or less uniform increase in width of all resonances. The experiment shows, however, that the widths of the remaining $`42\%`$ increase non-uniformly in $`d`$. This indicates that some of these resonances gain (or loose) width at the cost (or in favour) of other ones. Thus practically all of the resonances are affected by the mechanism of resonance trapping.

In conclusion, we have demonstrated experimentally that resonance trapping takes place in a microwave cavity coupled to a waveguide. This is, to the best of our knowledge, the first unambigous experimental verification of the effect. This demonstration was possible by tracing the motion of the resonance poles as a function of the opening of the slit starting close to the real axis and following them into the region of overlapping resonances. Our experimental proof of resonance trapping does not depend on any model assumptions.

Theoretical studies have shown that the phenomenon of resonance trapping does neither depend on the number of resonances nor on the special shape of the microwave cavity. The effect appears also in various open many-body quantum systems. In any case, the decoupling of some states from the decay channels takes place. However, the experimental verification of resonance trapping is important not only for an understanding of the properties of open quantum systems with overlapping resonances. More important is, maybe, the necessity of a good knowledge of the detailed properties of open quantum systems for the design of mesoscopic systems. By tuning the coupling strength to the decay channels, the properties of the system can be controlled. Further experimental and theoretical studies of this interesting topic, including the conductance, have to be performed.

Acknowledgment: The participation of Y.-H. Kim, R. Schäfer, and H. Schanze in performing the experiment is gratefully acknowledged. E. P. thanks the SFB 383 ’Unordnung in Festkörpern auf mesoskopischen Skalen’ for taking the costs of a three-month stay in Marburg. P. Thomas and F. Jahnke are thanked for stimulating discussions.


text:
# First search for gravitational wave bursts with a network of detectors

## Abstract

We report the initial results from a search for bursts of gravitational radiation by a network of five cryogenic resonant detectors during 1997 and 1998. This is the first significant search with more than two detectors observing simultaneously. No gravitational wave burst was detected. The false alarm rate was lower than 1 per $`10^4`$ years when three or more detectors were operating simultaneously. The typical threshold was $`H4\times 10^{21}Hz^1`$ on the Fourier component at $`10^3Hz`$ of the gravitational wave strain amplitude. New upper limits for amplitude and rate of g.w. bursts have been set.

preprint: HEP/123-qed

The direct detection of gravitational waves will be a watershed event for both the physics of gravitation and the investigation of compact astronomical objects. A variety of astrophysical events are expected to produce gravitational waves of short duration ($`1`$ sec), or gw bursts, such as the gravitational collapse of stars or the final few orbits and the subsequent coalescence of a close binary system of neutron stars or black holes . Due to the inherent weakness of such signals, and the difficulty in distinguishing them from a myriad of noise sources, the direct detection of a gw burst will require coincident detection by multiple detectors with uncorrelated noise. Searches for gw bursts over periods of observation of $`13`$ months have been performed in the past by pairs of cryogenic resonant bar detectors , setting upper limits on the incoming rate. A few days of observation have been reported for simultaneous operation of three cryogenic bar detectors  and, with much less sensitivity, of a pair of short-arm interferometric detectors . Upper limits on gw signals from coalescing binaries has been recently reported also by a single interferometric detector for $`25`$ hours of observation .

In the last few years, the increase of the number of cryogenic resonant detectors in simultaneous operation has greatly improved the prospects of obtaining a confident detection of gw bursts. There are now five operational cryogenic bar detectors: ALLEGRO (Baton Rouge, Louisiana, USA) , AURIGA (Legnaro, Italy) , EXPLORER (CERN) , NAUTILUS (Frascati, Italy)  and NIOBE (Perth, Australia) . The groups operating these detectors agreed in $`1997`$ to start a global search for short ($`1ms`$) gw bursts under common protocols, by establishing the International Gravitational Event Collaboration (IGEC) .

All these detectors use the same principles of operation. The gw excites the first longitudinal mode of the cylindrical bar, which is cooled to cryogenic temperatures to reduce the thermal noise and is isolated from seismic and acoustic disturbances. To measure the strain of the bar, a secondary mechanical resonator tuned to the cited mode is mounted on one bar face and a sensor measures the displacement between the secondary resonator and the bar face. The resulting noise of the detectors in terms of strain at the input is $`510\times 10^{22}/\sqrt{\mathrm{H}z}`$ in a bandwidth of $`1Hz`$ surrounding the two coupled-mode frequencies. Some of the important physical parameters of the five detectors are shown in Table 1. The detector response is optimal for a gw incoming perpendicular to the bar axis and polarized along it . The axes of all the bar detectors are aligned to within a few degrees of one another, so that the chance of coincidence detection is maximized. This makes the amplitude acceptance of the detectors for the Galactic Center direction greater than $`0.7`$ for about $`60\%`$ of the time .

Each detector output is processed by filters optimized for short gw bursts, giving the estimate for the Fourier component $`H(\omega )`$ of the strain amplitude $`h(t)`$ in the detection bandwidth of $`1Hz`$ around the mode frequencies listed in Table 1. More specifically, $`h(t)`$ is the gw amplitude multiplied by the antenna pattern of the detector . With the exception of the ALLEGRO detector, the noise of the detectors was typically not stationary over long observation times and was affected by some unmodeled noise sources, whose correlation with common environmental noise sources was found to be weak . Fig. 1 shows the variability of the noise of each detector during 1997-1998, in terms of the Fourier component of the gw corresponding to unity signal-to-noise ratio, $`H_{rms}`$. The detectors had quite similar noise levels, since the typical values of $`H_{rms}`$ were all within a factor of 3.

We point out that this search for bursts is suitable for any transient gw which shows a nearly flat Fourier transform $`H(\omega )`$ of its amplitude $`h(t)`$ at the two resonant frequencies of each detector. The metric perturbation $`h(t)`$ can either be a millisecond pulse, a signal made by a few millisecond cycles or a signal sweeping in frequency through the detector resonances. The IGEC search is therefore sensitive to different kinds of gw sources such as a stellar gravitational collapse , the last stable orbits of an spiralling $`NS`$ or $`BH`$ binary, its merging and its final ringdown . The computation of $`h`$ from the measured Fourier component $`H`$ requires a model for the signal shape. A conventionally chosen shape is a pulse lasting $`10^3s`$ ; in this case, $`H`$ should be multiplied by $`10^3Hz`$ to get the corresponding strain amplitude, $`h`$.

This letter reports the results of the first coincidence search for gw bursts performed by the IGEC observatory. The observations covered most of $`19971998`$, including $`625.0`$ days with at least one detector in operation, $`260.4`$ days with at least two detectors in simultaneous operation, $`89.7`$ days with three detectors, and $`15.5`$ days with four. This is the first search with significant observation time with more than two detectors. The duration of simultaneous operation would have been greater if it had been possible to operate these instruments with higher duty factors, which were typically $`50\%`$ during this period with the exception of ALLEGRO. More details on the observatory, its data exchange protocol and the exchanged data set can be found in Ref. .

The analysis of the data can be divided into two parts: a generation of candidate event lists for each of the individual detectors, and a time coincidence analysis using the lists. This approach, though not optimal, has the advantage of being easily implemented and provides for a satisfactory effectiveness.

Each IGEC group extracted the candidates for gw bursts, or events, by applying a threshold to the filtered output of the detector. The events were described by their Fourier magnitude $`H`$, their arrival time, the detector noise at that time and other auxiliary information. To limit the expected rate of accidental coincidences, each detector threshold was adaptively set to obtain a maximum event rate of $`100`$ / day, with typical values in the range $`H_{det}26\times 10^{21}Hz^1`$ corresponding to magnitude signal-to-noise ratio SNR $`35`$. Single spurious excitations are vetoed against disturbances detected by environmental sensors. The AURIGA detector checked each event against the expected waveform template by means of a $`\chi ^2`$ test . The lists of the events exchanged within IGEC by each detector also include declarations of the off- and on- times for the detectors.

All searches for coincident events used a time window of 1.0 second. This choice limits the false dismissal probability to less than a few per cent while it ensures a very low false alarm probability when at least three detectors are observing simultaneously. No three- and four-fold coincidence was detected, and therefore we did not identify candidates for gravitational wave detection in the $`89.7`$ days of three-fold observation. The detector thresholds were typically $`3\times 10^{21}Hz^1`$ for the most sensitive three-fold configuration (ALLEGRO-AURIGA-NAUTILUS) and $`5\times 10^{21}Hz^1`$ for the others. To give examples of detectable signals, these thresholds would correspond to respectively $`0.04`$ and $`0.11M_{}`$ converted to isotropic radiation in the optimal polarization at the distance of the Galactic Center ($`10kPc`$), assuming a g.w. burst of $`1ms`$ duration . For comparison, the signal expected from the last stable orbits of an optimally oriented $`NS`$ coalescing binary at $`10kPc`$ with $`2\times 1.4M_{}`$, would give $`H(\omega )=34\times 10^{21}Hz^1`$ at the detector resonant frequencies. The number and amplitude of the two-fold coincidences found in the $`260.4`$ days of two-fold observation are in agreement with the estimated accidental background .

The estimation of the false alarm rate is a crucial element in any gw search. It allows for the interpretation of any observed coincidences as well as the evaluation of the potential of the observatory. Since the events arrival times of each detector are randomly distributed with a non stationary rate, the expected background of accidental coincidences can be computed by two methods: i) by modeling the event times as Poisson point process and using the measured rates of events for each detector, and ii) by counting the coincidences after performing even time shifts of the event times of one detector with respect to the others .

In the first approach, the expected rate of accidental coincidences is

$$\lambda =N\frac{(\mathrm{\Delta }t)^{N1}}{T_{obs}^N}\underset{i=1}{\overset{N}{}}n_i,$$
(1)
where $`N`$ is the number of detectors simultaneously operating, $`T_{obs}`$ their common observation time, $`\mathrm{\Delta }t=1s`$ the maximum time separation for a coincidence, $`n_i`$ the number of events of the $`ith`$ detector during $`T_{obs}`$. This equation holds even if the event rates of detectors are not stationary as long as they are uncorrelated among different detectors.

The second approach is more empirical. In the case of the two-fold coincidence searches, the time shift results are in agreement with those predicted through Eq. 1 , and demonstrate that the event rates of different detectors are uncorrelated.

The capabilities of the IGEC observatory with respect to the false alarm probability are shown in Fig. 2 for a few sample configurations of the observatory. The accidental rate is calculated as a function of a signal amplitude threshold $`H_{thr}`$ at the detectors by applying Eq. 1 to the number of events of the detectors whose amplitude is $`H_{thr}`$. The typical time variability of the instantaneous accidental rate $`\lambda `$ has been calculated by means of a Monte Carlo simulation based on the measured non-stationary behavior of event rates on single detectors. This variability is about one order of magnitude with respect to the mean and is mainly determined by the non stationary performances of the detectors. The estimated mean background of two-fold coincidences is still fairly high, unless $`H_{thr}`$ is raised well above the data exchange threshold $`H_{det}`$. On the other hand, a three-fold or four-fold coincidence search keeps a high statistical significance even for $`H_{thr}H_{det}`$, since the expected accidental rates are low enough: respectively less than 1 false alarm per $`10^4`$ or $`10^6`$ years of observation at $`H_{thr}4\times 10^{21}Hz^1`$. These fall rapidly as $`H_{thr}`$ increases. In fact, the IGEC accidental background noise would remain negligible even after centuries of observation time.

The $`260`$ days of observation with two or more detectors in simultaneous operation improved by about a factor of three the previously set upper limit on the rate of gw bursts incident on the Earth . Assuming the emission is described by a stationary Poisson point process and using the same procedure as in Ref. , the limiting rate set for a gw burst emission from isotropically distributed sources is $`4year^1`$ for $`H_{gw}10^{20}Hz^1`$ (Fig. 2) with 95% confidence. A more complete analysis is in progress.

The IGEC observatory can also be used to set an upper limit on the amplitude of gw bursts corresponding to astronomical events, such as supernovae or gamma ray bursts. For time windows of the order of the hour or larger, each detector is likely to show accidental events and therefore this upper limit benefits from a multiple coincidence search among the operating detectors.

A sample of the upper limits on the amplitude of gw bursts occurring within a time span of $`1`$ hour, is shown in Fig. 3 for a few weeks of 1998, when up to four detectors were operating. These limits apply to the component of the radiation emitted with optimal polarization from a source optimally oriented with respect to the detectors. We are $`95\%`$ confident that there was no radiation above this level hour by hour. In all of 1998, the limits set by the IGEC observatory were better than $`H_{gw}=6`$ and $`4\times 10^{21}Hz^1`$ for $`94\%`$ and $`21\%`$ of the year respectively. To specialize these upper limits for a specific source direction, each detector response should be divided by its antenna pattern . The corresponding observation times of the Galactic Center by IGEC within the same limits have been respectively $`44\%`$ and $`7.5\%`$ of 1998. For a source at the Galactic Center emitting isotropically a $`1ms`$ burst , the above upper limits correspond to about $`0.16`$ and $`0.07M_{}`$ converted in the optimal polarization.

Finally, we remark that the IGEC observatory is capable of monitoring the strongest galactic sources with a very low false alarm probability when at least three detectors are simultaneously operating. All the groups involved are actively working for upgrading the current detector performances and therefore we expect in the near future to extend the observation range to the Local Group of galaxies, which means an increase of a factor of 10 of the observed mass.

The ALLEGRO group was supported by the National Science Foundation grant 9970742 and LSU, the NIOBE group by the Australian Research Council, and the Italian groups by the Istituto Nazionale di Fisica Nucleare and in part by a grant COFIN’97 from Ministero dell’Università e della Ricerca Scientifica e Tecnologica.


text:
# Mass Loss from Globular Clusters

## 1. Introduction

Among the reasons for studying this problem are the following:

1. There is a growing amount of direct observational evidence showing that globular clusters lose mass (cf. Sec.2 and the paper by Meylan in these Proceedings, and several of the poster papers).
2. Preferential loss of stars of low mass causes the mass function to evolve from its primordial form (cf. Sec.2).
3. It affects the statistical distribution of cluster properties in the galaxy, causing the death of clusters of low mass or large radius (Sec.2).
4. It has much in common with the problem of disruption of satellite galaxies, in which there has been much interest recently (e.g. Kroupa 1997, Ibata et al. 1997).
5. Decayed globular clusters may be a substantial contribution to the stellar halo of the Galaxy (cf. Ashman & Zepf 1998).
6. On the theoretical side, understanding mass loss turns out to be a vital issue in establishing the right boundary conditions to use in simplified models, such as Fokker-Planck models (Takahashi & Portegies Zwart 1998, 1999; Spurzem & Takahashi 2000).
7. Also on theory, mass loss turns out to be a key to understanding how to scale $`N`$-body models with $`N`$ (cf. Sec.3).

The emphasis in this review is on theoretical problems, but it begins with a summary of the observational position.

## 2. Observational Issues

### 2.1. Tidal tails and extensions

Two surveys (Grillmair et al. 1995; Meylan et al. 1999, Leon, Meylan & Combes 2000, and Meylan, these Proceedings) and some data on individual objects (Zaggia et al. 1997 for M55, Testa et al. 2000 for M92) have shown the presence of stars beyond the conventional tidal radius. In the case of the surveys the shapes tend to resemble those of tidal tails observed in simulations (cf. Fig.1), while the studies of M55 and M92 revealed roughly spherical extensions.

These studies raise a number of theoretical questions. In Grillmair et al. the cluster with the strongest tidal tail is M2 (=NGC7089), and modelling based on this data (Johnston et al. 1998) suggests that its lifetime is of order 1-3Gyr. Theoretical estimates, however, imply that this object is losing mass at one of the very slowest rates of all Galactic globular clusters, on a time of order 10–20 Hubble times (Gnedin & Ostriker 1997, Aguilar et al. 1988). The direction of the tail is not, at first sight, consistent with the most plausible cause (shocking during a recent passage through the Galactic disk). Meylan et al. themselves draw attention to two clusters in their sample with very similar orbits (and hence strength of tidal effects) but tails of very different strength. Furthermore the dependence on radius is flatter than in any of their simulations. The presence of a nearly spherical extra-tidal extension, as in the case of M92, raises the question of what the real value of the tidal radius is.

### 2.2. Kinematic Evidence

By studying the radial velocities of large numbers of stars in M15, Drukier et al (1998) found a slight rise in the velocity dispersion at large radii (but still well within the tidal radius). This was interpreted at first as evidence of tidal heating, but it was then realised, using $`N`$-body modelling, that this effect can occur even in clusters on circular galactic orbits. (This can also be seen in Giersz & Heggie 1997.) In this case the tide is steady (in a rotating frame), and its effect cannot be described as “heating” (Sec.3).

In two individual clusters, stars are observed with radial velocity exceeding the escape velocity (Gunn & Griffin 1979; Meylan, Dubath & Mayor 1991). As with extratidal extensions, the significance of these observations depends on the reliability of the models which yield the value of the escape speed. If taken at face value, however, such observations are usually interpreted as ejecta from three-body interactions (e.g. Davies 1992).

An alternative interpretation is suggested by recent $`N`$-body simulations (Fig.2). They show that the number of stars inside the tidal radius with energies above the energy of escape (“potential escapers”) is surprisingly large. Furthermore, the maximum mass of this population, $`M_{\mathrm{pe}}`$, decreases surprisingly slowly with increasing $`N`$: $`M_{\mathrm{pe}}/M(0)0.1(N/3\times 10^4)^{0.23}`$, where $`M(0)`$ is the initial cluster mass. Therefore even for a real cluster the proportion could be as large as $``$5%, and a considerably larger proportion of the current cluster mass. These results apply to models on circular galactic orbits, however.

### 2.3. Evidence from Mass Functions

At one time there seemed very clear evidence that the slope of the mass function in observed clusters is correlated with galactocentric radius and $`z`$, consistent with an interpretation in terms of preferential loss of low-mass stars (Capaccioli et al. 1993). The best that can be said now (Piotto & Zoccali 1999) is that, while there seems to be a trend, assessing its statistical significance is hard.

One cluster in which the effect of mass loss may be most pronounced is NGC 6712, which definitely appears to have a flatter mass function than any other comparable cluster (de Marchi et al. 1999), and has received considerable attention. Takahashi & Portegies Zwart (1999) gave an example of an $`N`$-body model which, after losing 99% of its mass, resembles this cluster. This is not implausible for a cluster whose present mass is estimated at $`10^5M_{}`$, and a tentative reconstruction of its orbit (Dauphole et al. 1996) suggests a very small galactocentric radius and therefore heavy tidal mass loss. A large mass loss rate for this object was also estimated by Gnedin & Ostriker (1997), though there is nothing about the appearance of this cluster to suggest that it is close to death.

### 2.4. Evolution of the Galactic Globular Cluster system

There have been many studies, especially in recent years, of the effects of mass loss and other dynamical processes on the statistical properties of the system of Galactic globular clusters. A few (e.g. Aguilar et al. 1988, Gnedin & Ostriker 1997) have given estimates for all individual clusters about which enough is known. In particular, Gnedin & Ostriker (following Fall & Rees 1977) plot a survival diagram of half-mass radius and cluster mass to show that the disruptive effects of relaxation, tidal shocking and dynamical friction roughly account for the present-day distribution of the Galactic globular clusters. Even at a statistical level the agreement is imperfect – there are too many low-mass clusters – which suggests that either initial conditions also have a role, or the dynamical processes are not yet sufficiently understood quantitatively.

## 3. Theoretical Issues

Several processes which lead to mass loss from star clusters have been suggested. In roughly the historical order in which they have come to prominence, they are

(i) two-body relaxation, which was the main process considered until the 1960s;

(ii) tidal effects;

(iii) mass loss from stellar evolution, which began to be allowed for routinely in this subject only after Applegate (1986); and

(iv) massive dark halo objects, a relatively speculative mechanism about which no more will be said here.

### 3.1. Isolated Systems and Systems with a Tidal Cutoff

Though highly idealised, the simplest problem of all would seem to be mass loss from equal-mass systems with no external field. In fact this problem has been more-or-less abandoned without being completely understood. The essential issue is whether the time scale of mass loss follows the relaxation time scale $`t_rt_{cr}N/\mathrm{ln}N`$ or the “strong encounter” time scale $`t_r^{}t_{cr}N`$, where $`t_{cr}`$ is the crossing time. It is not inconsistent to suppose that internal processes like core collapse proceed on the time scale $`t_r`$ while the time scale of escape, which may depend on energetic encounters, is more like $`t_r^{}`$. Long ago Hénon (1960, 1969) gave formulae for the escape rate involving essentially $`t_r^{}`$, and this prediction is confirmed by $`N`$-body results (Giersz & Heggie 1994). Only slightly more recently, however, Spitzer & Shapiro (1972) gave a theory for isolated systems which gives escape on a time scale $`t_r`$.

These issues are not irrelevant when we turn to tidally limited systems, at least if the tide is represented as a cutoff. There is no reason why Hénon’s beautiful but little-known general formula for the escape rate (Hénon 1960), i.e.

$$\dot{N}=\frac{256}{3}\sqrt{2}\pi ^4G^2m^2_0^{r_{max}}r^2𝑑r\frac{f(E)f(E^{})(E+E^{}\varphi )^{3/2}}{E^2}𝑑E𝑑E^{},$$

which leads to an escape time scale $`t_r^{}`$, is inapplicable, and yet virtually all theorists would adopt a time scale based on the usual relaxation time $`t_r`$. One reason for this is that so much work is based on the Fokker-Planck equation, which has no other time scale. Among theorists one exception is Wielen (e.g. 1988), who argues on the basis of $`N`$-body data that strong encounters greatly dominate for clusters with a few hundred stars. Since the Coulomb logarithm varies so slowly with $`N`$ he concludes that they still dominate for systems of the size of globular clusters.

### 3.2. Steady Tides

Even in a steady tide, e.g. a cluster in a circular galactic orbit, it is useful to distinguish three effects: (i) the finite boundary, usually taken at the radius of the Lagrangian point (cf. Spitzer 1987); (ii) the reduced escape energy; and (iii) “induced” mass loss, which is described further below.

#### Definition of Escape

Several escape criteria are in common use, and all are wrong. Unless the potential of the cluster changes, a star may have an energy above the escape energy and yet remain within the tidal radius, $`r_t`$ indefinitely (Fukushige & Heggie 2000 and references therein). Within the usual tidal approximation it is possible for a star to recede to an arbitrary distance and eventually return to the cluster (Ross, Mennim & Heggie 1997), and so a simple distance or apocentre criterion fails. Ross et al. give a rigorous criterion, but in practice what is usually done in $`N`$-body simulations is to measure the mass within $`r_t`$, and then the escape criterion for simplified models (gas or Fokker-Planck), though based on physical ideas, is adjusted for reasonable agreement (e.g. Portegies Zwart & Takahashi 1999).

#### Induced Mass Loss

This non-standard term is used here to describe the following fact. If a cluster in a steady tide loses mass for any reason (e.g. mass loss of stellar evolution), then the tidal radius decreases and the escape threshold changes, so that still more mass is lost as the cluster readjusts. This additional mass loss is referred to as “induced mass loss”. If the cluster has little mass near the tidal radius or near the escape energy then induced mass loss is small. In King models it is a few percent of the initial mass loss. In Woolley models, however, in which the velocity distribution is a Maxwellian truncated at the escape speed, it is a much large factor, and in fact models with a Woolley concentration parameter $`W_0<4`$ are unstable to induced mass loss.

Consider a cluster in which the tide is treated as a cutoff, the cluster loses mass by stellar evolution, and relaxation is very slow. (Weinberg (1993) gives results of such models.) In this situation even an initial King-like distribution evolves to one more like a Woolley model. It seems likely that the above instability is responsible for the rapid loss of mass which such models appear to exhibit (Chernoff & Weinberg 1990, Fukushige & Heggie 1995) at the end of their life. Clusters in which mass loss is dominated by relaxation do not show this behaviour. The distinction between these two kinds of evolution is sometimes referred to as the “ski-jump” problem (Fig.3), and was highlighted by Portegies Zwart et al. (1998).

#### The Rate of Mass Loss

The traditional picture (in the absence of mass loss by stellar evolution) is that stars gain sufficient energy to escape on the relaxation time scale, $`t_r`$, while escape itself takes place on a crossing time scale. As the latter is much shorter when $`N`$ is large, it follows that the time scale of mass loss is proportional to $`t_r`$. $`N`$-body simulations give results that are sufficiently different from this prediction (Fig.4) that extrapolation to, say, $`N10^6`$, which would be required for scaling these models to a typical globular cluster, would yield results differing by a factor of order two.

Much work has been done on this so-called “scaling” problem since its discovery. The problem is reduced somewhat if self-consistent initial conditions are used, i.e. a model which is in equilibrium within the tidal field, unlike a standard King model (cf. Fukushige & Heggie 2000). The main problem, however, is the time scale on which stars escape, and the resulting population of “potential escapers” (Fig.2). It has been known for a long time (King, pers. comm.) that stars with energy above the escape energy may remain in a cluster for many crossing times before, if ever, finding a way out, but it is only recently that this aspect has been quantified (Fukushige & Heggie 2000) as a function of energy. These authors also determined the fraction of stars with energy above the escape energy which would never escape (if the potential remained fixed). How this new understanding resolves the scaling problem is a story which has been taken up by Baumgardt (these proceedings, and Baumgardt 2000).

### 3.3. Unsteady Tides: Bulge Shocking

Unsteady tides usually considered are of two kinds: (i) bulge shocking, (ii) disk shocking. The dynamical problems they pose are similar, and the remainder of this review concentrates heavily on bulge shocking. It happens only to clusters on elliptic orbits, however, whereas disk shocking would affect any cluster except one lying strictly in the galactic plane.

#### Theory of Bulge Shocking

This mechanism is often referred to as bulge heating, and indeed considerable emphasis is placed on its effect on the energy of stars and clusters. The effect of the heating, however, is to accelerate the expansion of the halo of a cluster across the tidal boundary, and so this mechanism is a legitimate subject of interest in this review.

Spitzer (1987) treated the problem impulsively, with qualitative discussion of the so-called adiabatic correction (for the fact that the orbital period of a star within the cluster need not be much greater than the time scale of variation in the tide). The argument was quantified and elaborated by Aguilar, Hut & Ostriker (1988), using an exponential form of adiabatic correction due to Spitzer in another context. The next interesting development was Weinberg’s realisation (Weinberg 1994a,b,c) that adiabatic invariants do not protect resonant combinations of frequencies in multi-dimensional systems (such as three-dimensional motion in a fixed cluster potential). (Some of the dynamical issues are ably expounded also in Henrard 1982 and Sridhar & Touma 1996.) The effect of this is a significant change from an exponential adiabatic correction to a power law, as has been checked (with some empirical improvements) in simulations (Gnedin & Ostriker 1999).

One application of these results is their incorporation into Fokker-Planck models. This has been accomplished by Weinberg (1994c) and by Gnedin & Ostriker (1997), and applied to the evolution of clusters by Murali & Weinberg (1997) and Gnedin, Lee & Ostriker (1999). Some results must be treated with caution, however, as the formulation of Murali & Weinberg results in bulge heating even in the case of circular orbits, where the tide is steady and cannot change the energy of a star.

Such work can give a first answer to the question whether it is tidal shocking or two-body relaxation which dominates the mass loss of galactic clusters. One answer, based on the work of Gnedin & Ostriker (1997), is given in Fig.5. When we come to study simulations in the next sub-subsection, this is a question on which we shall concentrate.

#### Simulations of Bulge Shocking

Much can be learned from the motion of test particles in a fixed cluster-like potential orbiting in a galaxy-like potential. Keenan & Innanen (1975) used a King-like cluster model and both point-mass and axisymmetric galaxy potentials, while Keenan (1981ab) used point masses for both. This case is the “elliptic restricted three-body problem”, which has experienced a modest renaissance recently in other contexts (e.g. Benest 1998 and references therein).

Relaxation can be added to this simple model by incorporating diffusion (Oh, Lin & Aarseth 1992), and in this way Oh & Lin (1992) concluded that the presence of a tidal field could suppress escape. Perhaps this result, which at face value contradicts all other wisdom on this question, may be explained by the increase in the escape time scale in a tidal field, as already discussed.

Self-consistent calculations of tidal effects on clusters on elliptical orbits were carried out by Chernoff & Weinberg (1991), though the initial models considerably exceeded their tidal radii, and therefore perhaps relate better to young clusters than to the old Galactic globular clusters. Much the same may be said of the later studies by Johnston, Sigurdsson and Hernquist (1999), though they took care to quantify the effects of relaxation. This was also done in the simulations by Combes et al. (1999), though the choice of initial conditions is too patchy to arrive at general conclusions on the relative importance of bulge heating and relaxation.

An independent attack on this question was begun by Baumgardt (1997), who computed some small $`N`$-body models of clusters on elliptical orbits about a point mass galaxy with initial conditions similar to those of Murali & Weinberg (1997), i.e. the cluster is started at apocentre with a limiting radius equal to the tidal radius at pericentre.

One attraction of these initial conditions is that they correspond to one of the conclusions of Oh & Lin (1992), which was that the limiting radius of clusters is determined at pericentre and maintained all along the orbit. Another attraction is that they allow a systematic exploration of the relative importance of tides and relaxation. Collaborative work with H. Baumgardt is currently under way at Edinburgh on this problem (Fig.6).

One of the interesting questions which even this research will leave unanswered is the behaviour of clusters in axisymmetric potentials. The essential difference is that the perigalactic distance varies considerably from one bulge passage to the next. It may be that mass loss through bulge shocking is much less regular and more infrequent than is suggested by Fig.6.

## 4. Summary and Conclusions

Observation of tidal tails of Galactic globular clusters now seem well established and consistent, but the detailed modelling and understanding of various individual objects remains almost unexplored. Kinematic data also suggests that tidal effects may impose an observable signature, while it is suggested in this paper that stars apparently observed above the escape energy may be genuine cluster members trapped for long periods in the combined field of the cluster and galaxy. It seems certain that mass loss has altered both the stellar mass function of each cluster and the distribution of cluster properties throughout the galaxy. The expected correlations among the stellar mass functions are a little elusive, however, and it is only in extreme cases that the effects of mass loss seem clear cut.

There are still loose ends to be tied up in the theory of escape from a cluster with a tidal cutoff: is the escape time scale proportional to the relaxation time, with the Coulomb logarithm, or is it dominated by single energetic scatterings, as in the formula of Hénon? The “induced” loss of mass may help to understand the distinction between “skiing” and “jumping” models, but a quantitative understanding of this distinction is still missing. When the tide is represented by a proper external field, the effect of the large population of “potential escapers” complicates the time scale for mass loss, and a theory of this problem has recently been developed by Baumgardt. In the more realistic case of an elliptic galactic orbit, an interesting focus for current research is the question of which mass loss process dominates: relaxation, or bulge shocking? An interesting future development will be to understand the effect of the varying perigalactic distance in axisymmetric potentials.

##### Acknowledgments.

I thank G.A. Drukier, G. Meylan and R. Wielen for additional information and comment on their results, and H. Baumgardt for many discussions on theoretical matters. $`N`$-body calculations were performed on hardware funded under SERC/PPARC grants GR/H93941 and GR/J79461, and supplied by the GRAPE team at the University of Tokyo, and partly using NBODY4 by S.J. Aarseth. Work with H. Baumgardt is supported by PPARC under grant 1998/00044. I also thank the organisers of Star2000 for their financial support.

## References

Aguilar L., Hut P., Ostriker J.P., 1988, ApJ, 335, 720

Applegate J.H., 1986, ApJ, 301, 132

Ashman K.M., Zepf S.E., 1998, Globular Cluster Systems (Cambridge: Cambridge University Press)

Baumgardt H., 1997, PhD thesis, University of Heidelberg

Baumgardt H., 2000, preprint

Benest D., 1998, A&A, 332, 1147

Capaccioli M., Piotto G., Stiavelli M., 1993, MNRAS, 261, 819

Chernoff D.F., Weinberg M., 1990, ApJ, 351, 121

Chernoff D.F., Weinberg M., 1991, in ASP Conf. Ser. Vol. 13, The Formation and Evolution of Star Clusters, ed. K. Janes (San Francisco: ASP), 373

Combes F., Leon S., Meylan G., 1999, A&A, 352, 149

Dauphole B., Geffert M., Colin J., Ducourant C., Odenkirchen M., Tucholke H.-J., 1996, A&A, 313, 119

Davies M.B., 1992, PhD thesis, Harvard University

de Marchi G., Leibundgut B., Paresce F., Pulone L., 1999, A&A, 343, L9

Drukier G., Slavin S.D., Cohn H.N., Lugger P.M., Berrington R.C., Murphy B.W., Seitzer P.O., 1998, AJ, 115, 708

Fall S.M., Rees M.J., 1977, MNRAS, 181, 37P

Fukushige T., Heggie D.C., 1995, MNRAS, 276, 206

Fukushige T., Heggie D.C., 2000, MNRAS, submitted

Giersz M., Heggie D.C., 1994, MNRAS, 268, 257

Giersz M., Heggie D.C., 1997, MNRAS, 286, 709

Gnedin O.Y., Lee H.M., Ostriker J.P., 1999, ApJ, 522, 935

Gnedin O.Y., Ostriker J.P., 1997, ApJ, 474, 223

Gnedin O.Y., Ostriker J.P., 1999, ApJ, 513, 626

Grillmair C.J., Freeman K.C., Irwin M., Quinn P.J., 1995, AJ, 109, 2553

Gunn J.E., Griffin R.F., 1979, AJ, 84, 752

Heggie D.C., Giersz M., Spurzem R., Takahashi K., 1998, in Highlights of Astronomy, vol.11A, ed J. Andersen (Kluwer, Dordrecht), p.591

Heggie D.C., Ramamani N., 1995, MNRAS, 272, 317

Hénon M., 1960, Ann. d’Astr., 23, 668

Hénon M., 1969, A&A, 2, 151

Henrard J., 1982, in NATO ASI Ser. C, 82, Applications of Modern Dynamics to Celestial Mechanics and Astrodynamics, ed. V. Szebehely (Dordrecht: Reidel), 153

Ibata R.A., Wyse R.F.G., Gilmore G., Irwin M.J., Suntzeff N.B., 1997, AJ, 113, 634

Johnston K.V., Sigurdsson S., Hernquist L., 1999, MNRAS, 302, 771

Keenan D.W., 1981ab, A&A, 95, 334 and 340

Keenan D.W., Innanen K.A., 1975, AJ, 80, 290

Kroupa P., 1997, New Astron., 2, 139

Leon S., Meylan G., Combes F., 2000, astro-ph/0006100

Meylan G., Dubath P., Mayor M., 1991, ApJ, 383, 587

Meylan G., Leon S., Combes F., 1999, astro-ph/9912041

Murali C., Weinberg M., 1997, MNRAS, 288, 749

Oh K.S., Lin D.N.C., Aarseth S.J., 1992, ApJ, 386, 506

Oh K.S., Lin D.N.C., 1992, ApJ, 386, 519

Piotto G., Zoccali M., 1999, A&A, 345, 485

Portegies Zwart S.F., Hut P., Makino J., McMillan S.L.W., 1998, A&A, 337, 363

Portegies Zwart S.F., Takahashi K., 1999, CeMDA, 73, 179

Ross D.J., Mennim A., Heggie D.C., 1997, MNRAS, 284, 811

Spitzer L., Jr, 1987, Dynamical Evolution of Globular Clusters. Princeton U.P., Princeton.

Spitzer L., Jr, Shapiro S.L., 1972, ApJ, 173, 529

Spurzem R., Takahashi K., 2000, preprint

Sridhar S., Touma J., 1996, MNRAS, 279, 1263

Takahashi K., Portegies Zwart, 1998, ApJ, 503L, 49

Takahashi K., Portegies Zwart, 1999, astro-ph 9903366

Testa V., Zaggia S.R., Andreon S., Longo G., Scaramella R., Djorgovski S.G., de Carvalho R., 2000, A&A, 356, 127

Weinberg M.D., 1993, in ASP Conf. Ser. Vol. 48, The Globular Cluster-Galaxy Connection, ed. G.H. Smith, J.P. Brodie (San Francisco: ASP), 689

Weinberg M.D., 1994abc, AJ, 108, 1398, 1403 and 1414

Wielen R., 1988, in IAU Symp.126., The Harlow-Shapley Symposium on Globular Cluster Systems in Galaxies, ed. J.E. Grindlay, A.G.D. Philip (Dordrecht: D. Reidel), p.393

Zaggia S.R., Piotto G., Capaccioli M., 1997, A&A, 327, 1004