1 - Astronomical Institute `Anton Pannekoek', University of Amsterdam, Kruislaan 403, 1098 SJ, Amsterdam, The Netherlands
2 - School of Physics and Astronomy, University of Southampton, Southampton, Hampshire SO17 1BJ, UK
3 - Department of Physics and Astronomy, Michigan State University, East Lansing MI, USA

Abstract
Optical and X-ray studies of six nearby galaxies show that the probability a globular cluster will be an X-ray source is consistent with being linearly proportional to its mass. We show that this result is consistent with some recent estimates of the velocity kick distributions for isolated radio pulsars - those which are the sum of two Maxwellians, with the slower distribution at about 100 km s^-1 - so long as a large fraction of the retained binaries are in binary systems with other massive stars. We confirm that over a large sample of galaxies, metallicity is clearly a factor in determining whether a globular cluster will contain an X-ray binary, and we estimate the transformations between color and metallicity for a large number of optical filter combinations. We also show that the core interaction rate is roughly linearly proportional to the stellar mass of a globular cluster for the Milky Way when one bins the clusters by mass.

Key words: stellar dynamics - binaries: close - globular clusters: general - galaxies: star clusters - X-rays: binaries - galaxies: elliptical and lenticular, cD

It has long been known that globular clusters contain much larger numbers of X-ray binaries per unit stellar mass than do field populations of galaxies (Katz 1975; Clark 1975). This overabundance has been ascribed to stellar interactions in globular clusters which allow neutron stars and/or black holes to enter new binary systems through either tidal capture (Clark 1975; Fabian et al. 1975) or exchange interactions (Hills 1976) long after the supernovae that produce them. In the Chandra era, it has been possible to extend this work to show that the clusters with the highest interaction rates are most likely to contain X-ray sources, and probably are most likely to contain accreting neutron stars (Pooley et al. 2003; Heinke et al. 2003; Gendre 2005).

While studies of Galactic globular clusters have thus been fruitful, the Milky Way's globular cluster system is rather small (containing only about 150 clusters), so it is not possible to study certain phenomena due to lack of statistical significance in a small sample of clusters. Furthermore, some of the Milky Way's globular clusters' parameters, such as metallicity and galactocentric radius, are strongly correlated with one another, making it difficult to isolate the causes of certain effects. As a result, extragalactic globular cluster systems can be invaluable in complementing the Galactic globular clusters for producing the best possible data on how globular cluster parameters affect X-ray binary production.

While ROSAT was able to observe globular cluster X-ray sources in M 31 (see e.g. Supper et al. 1997), it was not until the advent of the Chandra X-ray Observatory that there was sufficient angular resolution to resolve apart the sources in more distant galaxies, or to obtain sufficiently accurate positions that comparisons could be made with optical counterparts. Early on, by comparing the positions of X-ray sources with those of optically detected globular clusters, it was found that large fractions of the X-ray binaries in elliptical galaxies were in globular clusters (Sarazin et al. 2001; Angelini et al. 2001 - ALM). Larger samples of galaxies have shown a trend where the fraction of X-ray binaries in globular clusters in a galaxy seems to vary as a function of galaxy type, increasing from spiral galaxies to lenticular to normal ellipticals to cD galaxies (Maccarone et al. 2003; Juett 2005; Irwin 2005).

From the observations of elliptical galaxies, with their large globular cluster systems, it has been possible to find correlations between cluster properties and the probability that a cluster would contain an X-ray binary which are not statistically significant in the Milky Way, or which have ambiguous interpretations in the Milky Way because of the aforementioned correlations between different cluster parameters. Specifically, it has been shown clearly that metal rich clusters are far more likely to contain X-ray binaries than metal poor clusters (Kundu et al. 2002), confirming suggestive results from the Milky Way (Grindlay 1993; Bellazzini et al. 1995). These results have been found again in numerous subsequent papers (see e.g. Jordan et al. 2004; Minniti et al. 2004; Xu et al. 2005; Kim et al. 2005). It has also been shown that more luminous clusters are more likely to contain X-ray binaries than less luminous clusters (see e.g. ALM; KMZ). We refer the reader also to reviews by Verbunt & Lewin (2006) and Fabbiano (2006) for a broader overview of the literature.

Other results have been suggested by the data, but are not as clearly significant. While in the Milky Way, it is clear that denser globular clusters are more likely to contain large numbers of X-ray sources per unit stellar mass (Pooley et al. 2003; Heinke et al. 2003; Gendre 2005), the extragalactic clusters have core radii which have not been sufficiently well resolved spatially as to establish this effect clearly. KMZ found an anti-correlation between cluster half-light radius and LMXB hosting probability in NGC 4472, but noted that it was statistically insignificant. Jordan et al. (2004) estimated King model parameters of globular clusters in M 87. They then inferred stellar collision rates ( $\Gamma$ ) from the model fits, and found that the collision rates were strongly correlated with the likelihood that a cluster would contain an X-ray binary, but that this correlation was only marginally better than the correlation with cluster luminosity. The absence of strong signatures of cluster concentration should not be taken as evidence that cluster concentration is uncorrelated with X-ray binary production. Half-light radius is not expected to be very well correlated with core concentration, making the null result of KMZ unsurprising. Because the core radii of globular clusters in the Virgo Cluster are typically a small fraction of a pixel on the Advanced Camera for Surveys aboard the Hubble Space Telescope, neither is it surprising that Jordan et al. (2004) were unable to find a strong signature of cluster concentration.

In this paper, we analyze a sample of six galaxies for which there is good optical and X-ray data. We show unambigously that the metallicity effect is strong. We then isolate the effect of cluster mass, M, on probability that a cluster will contain an X-ray source, and find that the probability a globular cluster will contain an X-ray source scales as $M^{1.03\pm0.12}$ . We compare this value with the expected value from various kick velocity distributions, and find that this result is most consistent with double Gaussian distributions for pulsar kick velocities, where there exists a slow mode to the neutron star kick velocity distribution.

$\begin{figure} \par\renewcommand{\epsfsize}[2]{0.7 ...$	Figure 1: The globular clusters analyzed in this data set. Small dots represent optically detected globular clusters, while the larger open circles surround the clusters with X-ray sources.
Open with DEXTER

2 Observations

We use data from HST and Chandra for six galaxies: NGC 1399, NGC 3115, NGC 3379, NGC 4472, NGC 4594, and NGC 4649. The data analysis procedures used here are the same as those used in our previous work (see Maccarone et al. 2003, for details). We take the distances to these galaxies from Tonry et al. (2001). Partial results for all of these galaxies have been presented previously (ALM; Kundu et al. 2003, 2004; KMZ; Di Stefano et al. 2003), and a detailed analysis with full source catalogs will be presented in future work (Kundu et al., in prep.) In all, our sample includes 98 X-ray sources in 2276 globular clusters. These data points are plotted in Fig. 1.

Because the data are taken in different filter sets (B-I for NGC 1399 and V-I for the rest of the galaxies), we must develop an algorithm to convert from magnitudes to cluster masses and metallicities. The I-band magnitudes are used to estimate the cluster masses, while the colors are used to estimate the cluster metallicities. Following the procedure of Kundu & Whitmore (1998), we determine the color-to-metallicity conversions for all commonly used sets of photometric filters. We take all globular clusters from the Harris (1996) catalogue for which there are data in a given pair of filters, for which there is a spectroscopic metallicity measurement, and for which E(B-V)<0.4. We then de-redden the data, using the E(B-V) versus extinction relations from Cardelli et al. (1989), and fit [Fe/H] as a function of color and vice versa, and take the bisector of the two fits. This typically produces a color-metallicity relation with a scatter of about 0.2 dex in [Fe/H]. The fact that the scatter is typically so small (i.e. about the same size as what would be expected from measurement errors) indicates that linear color-metallicity relations are adequate for most purposes. The results are given in Table 1.

3 Analysis

Assuming that the probability that a globular cluster contains an LMXB can be described as a power law function of the cluster's mass times a power law function of its metallicity, we attempted to determine the exponents of the power laws. First we attempted to do the fitting using density estimation (Silverman 1986), applying Gaussian smoothing to the two dimensional data set with different smoothing factors. We began by making a grid of $100\times100$ in mass and metallicity, and tried smoothing factors from 3 to 10 (i.e. smearing with a Gaussian kernel with a full width half maximum ranging from 3 to 10 cells on the grid).

The alternative is to use binning. We then sorted the clusters by mass, and compared the total number of X-ray binaries in logarithmically spaced mass bins (varying the number of bins from 5 to 10) with the expression Y, defined as:

$\begin{figure} \par\renewcommand{\epsfsize}[2]{0.7 ...$	Figure 2: Metallicity corrected LMXB probability as a function of mass, fit using binning, with 7 bins. Fits with different numbers of bins produce results which are the same within statistical errors. Changing the metallicity correction between 0.20 and 0.30 has no effect on the mass exponent, as expected given the absence of correlations between cluster mass and metallicity.
Open with DEXTER

We then use minimization of the Cash statistic (Cash 1979) to estimate the mass exponent, and find that $P({\rm LMXB})\propto Y^{1.03\pm0.12}$ (see e.g. Fig. 2 for a fit using 7 bins), which means that $P({\rm LMXB})\propto M^{1.03\pm0.12}$ , since the exponent for Y is independent of the exponent for the metallicity term, meaning that all the dependence on Y is due to its mass dependence. We compare the results using different numbers of bins and find that the results are not strongly sensitive to the number of bins used (for numbers of bins in mass between 5 and 10, the best fitting value varies from 1.01 to 1.07, with a standard deviation of 0.11 or 0.12 for every fit). For a given number of bins, varying the metallicity exponent put into Y changes the fitted dependence of $P({\rm LMXB})$ on Y by only one part in one thousand - much less than the changes due to using different numbers of bins. We have also made linear-logarithmic and logarithmic-linear plots of the data and found that no reasonable fit can be made using such functions. Thus while a power law is obviously not a unique description of the data, it seems to be the most reasonable simple functional form for the data, and a reasonable choice for parameterizing the data. We note that for a smaller number of clusters in M 87, Jordan et al. (2004) found a metallicity exponent of $0.3\pm0.1$ and a mass exponent of $1.08\pm0.11$ , so our results are consistent with past work.

4 Discussion

4.1 Implications of the observational results for retention fractions

The number of X-ray binaries in a globular cluster is likely to be correlated most strongly with the stellar collision rate. The collision rate should be dominated by the collisions which take place in the core, and it is straightforward to compute analytically the collision rate in the core of a globular cluster:

The only galaxy where the core radii and central densities of a large number of clusters have been reliably measured (rather than inferred from model fits) is the Milky Way. The collision rates for the Milky Way are well correlated with the masses; when we bin the Milky Way's globular clusters by mass and fit a power law to the data. We find that the best fitting power law index varies between 0.9 and 1.3, depending on the number of bins used (see Fig. 3 for a characteristic fit). The statistical errors on the fits are always smaller than 0.1 in index. We thus adopt $1.1\pm0.2$ for the exponent here, but note that this is a systematic error, and really represents a hard bound on the possible values, rather than a $1\sigma$ error. This result is intermediate between the two cases where the spatial properties are independent of mass (i.e. where, among the King model parameters, only the central density varies, and the expectation value would be $\Gamma\propto M^{1.5}$ ) and where the concentration index (i.e. the ratio of tidal radius to core radius) is independent of mass (in which case the expectation value would be $\Gamma\propto M^{2/3})$ . The intermediate result is expected given that cluster concentration increases with increasing cluster mass (Djorgovski & Meylan 1994). It is possible to find other means of producing various exponents besides those presented above, if one allows for different relations between the fraction of the mass in the cluster core and the core radius with mass. The key point, though, is that empirically, the relation for Milky Way globular clusters is that $\Gamma$ scales approximately as M^1.1.

$\begin{figure} \par\renewcommand{\epsfsize}[2]{0.7 ...$	Figure 3: Binned collision rate versus mass for the Milky Way's globular clusters, using data from the Harris catalog with 8 bins.
Open with DEXTER

Jordan et al. (2004) found that the probability a cluster in M 87 will contain an X-ray binary is proportional to $\Gamma \rho_{\rm c}^{-0.4\pm0.1}$ , and suggested that the deviation from a purely linear relationship with $\Gamma$ might be due to destruction of binaries by dynamical interactions in the globular cluster. While clusters do undoubtedly destroy binary systems, this should be done preferentially for the longest period systems, perhaps by making them highly eccentric rather than by destroying them outright if they are Roche lobe overflowing systems(e.g. Hut & Paczynski 1984; Rasio & Heggie 1995; Maccarone 2005). Binary destruction should have relatively little impact on the short period systems which, due to their higher duty cycles (see e.g. Piro & Bildsten 2002; Bildsten & Deloye 2004; Ivanova & Kalogera 2006) should provide most of the bright X-ray binaries which can be seen out to Virgo Cluster distances (see Maccarone et al. 2005, for a discussion of the effects of destruction of X-ray binaries due to dynamical interactions on their orbital period distributions).

In light of our above finding that, empirically, for the Milky Way, $\Gamma$ scales with cluster mass to the power $1.1\pm0.2$ , it seems simpler just to consider that the estimates of the core radii of the globular clusters in M 87 are likely to be dominated by measurement uncertainties. In the absence of useful radial information, the derived central density of a globular cluster will scale linearly with the cluster's mass; the estimate of $\Gamma$ will then scale with the cluster mass to the power 3/2. We have shown that, for the Galactic globular cluster systems, $\Gamma$ scales linearly with the cluster mass; therefore, one will need to correct this derived value of $\Gamma$ by a term of order M^1/2 (or $\rho^{1/2}$ in the absence, again, or any useful radial information). The $\rho^{-0.4}$ term found by Jordan et al. (2004) is thus more likely this "correction'' term than a term indicating dynamical destruction of accreting binaries. The only other means of producing the $\Gamma$ -M correlation we show and the relation among $\Gamma,M,$ and $\rho_{\rm c}$ found by Jordan et al. (2004) would be a specific relationship between the fraction of mass in the core and the core radius for clusters, in which the cluster core density is uncorrelated with the cluster mass. This would, in principle, be possible, e.g. if the most massive clusters has progressively lower fraction of their mass in their core but is, in fact, not the case, given the strong observed correlation between cluster central density and cluster absolute magnitude (Djorgovski & Meylan 1994).

Finally, let us consider an additional important point. The collision rate $\Gamma$ is a collision rate of typical stars with one another - mostly main sequence stars with main sequence stars. The collision rate for neutron stars will be the collision rate itself, times the fraction of stars in the cluster core which are neutron stars. This fraction will come from three factors - stellar evolution, which determines what fraction of stellar mass ends up locked up in neutron stars, the retention fraction, which characterizes what fraction of those neutron stars remain in the cluster, and the effects of mass segregation, which determines how over-represented neutron stars are in the cores of globular clusters, compared to their representation in the clusters on the whole. We assume that the stellar evolution effects have no mass dependence, so any difference between the collision rate as a function of cluster mass and the LMXB hosting probability as a function of cluster mass is most naturally explained by the dependence of the retention fraction on cluster mass, and mass segregation variations as a function of mass.

As a result, we define a quantity, $\Phi$ , which should scale linearly with the expected formation rate of neutron star X-ray binaries in the core of a globular cluster. We add only a single assumption at this point, which is that mass segregation works sufficiently effectively on timescales of the ages of globular clusters that all the neutron stars are in the cluster cores. This introduces an additional factor of $M/M_{\rm core}$ to account for the "overdensity'' of neutron stars among the total stars in a cluster's core. Since $M_{\rm core}\propto r_{\rm c}^3\rho_{\rm c}$ , this yields:

$\begin{displaymath}\Phi=\Gamma M/M_{\rm core} = \rho^{1/2}M/r_{\rm c}. \end{displaymath}$

(4)

We also note that mass segregation of neutron stars is not sufficiently fast as to place all neutron stars in the cores of globular clusters; there are some millisecond pulsars with well-measured positions that place them outside the cores of their host clusters (see e.g. Camilo & Rasio 2005, for a reasonably up-to-date list of globular cluster pulsars, or Paulo Freire's continuously updated list at http://www.naic.edu/~pfreire/GCpsr.html). While, relaxing this assumption leads to a weaker dependence of the retention fraction on mass, the assumption is generally a good one. About half of all millisecond pulsars are within one core radius of the center of their host clusters, while the fraction of the total mass in the core is less than 10% even in 47 Tuc, which has the highest core fraction of any massive cluster.

$\begin{figure} \par\renewcommand{\epsfsize}[2]{0.7 ...$	Figure 4: Binned values of $\Phi$ versus mass for the Milky Way's globular clusters, using data from the Harris catalog with 6 bins.
Open with DEXTER

4.2 Scaling relations of the retention fraction

We have considered a variety of possible kick velocity distributions in order to determine whether any of them are consistent with the estimated retention fraction relation - $F_{\rm ret}\propto M^{0.5}$ . Essentially, two broad classes of velocity kick distributions have been proposed for radio pulsars: "zero-peaked'' distributions (e.g. Paczynski 1990; Hartman 1997; Hansen & Phinney 1997), where the most likely initial velocity for the neutron star is close to zero (although the distribution retains significant probability density out to hundreds of km s^-1), and distributions which are Gaussians peaked at the mean pulsar velocity (e.g. Lyne & Lorimer 1994) or Maxwellians, which have rather similar characteristics. Some recent work has modelled the pulsars' velocity distributions as double Gaussians or double Maxwellians (Arzoumanian et al. 2002; Brisken et al. 2003), typically with one peak at about 100 km s^-1 and containing about 20% of the neutron stars and the other peak at about 300-500 km s^-1 containing about 80% of the neutron stars. In essence, the two approaches have converged to rather similar distributions, since the low velocity mode "fills in'' the gap left at low velocities in the single mode approach. On the other hand, the most recent estimation of the kick velocity distribution argued for a return to a single-peaked Maxwellian distribution with a characteristic 1D velocity of about 265 km s^-1 (Hobbs et al. 2005), so there still exists controversy about the initial velocity distribution of radio pulsars.

4.2.1 Results of models tried

We have tried a variety of different models with different parameter values to determine whether they can produce retention fractions with negligible mass dependence. We present the methodology and a sampling of the results here. In the interests of brevity, we do not present every calculation we have done, but note that these are in the Master's degree thesis of the lead author of this paper albeit with a somewhat different interpretation section and can be obtained by contacting the authors of this paper.

To calculate reliable retention fractions we use 111 globular clusters from the Harris catalogue as templates. We used the core and tidal radii from the catalogue, converting arcminutes to parsecs using the distances provided in the catalogue. We also converted the V magnitude to mass assuming a constant mass-to-light ratio of 1.8, as found by Piatek et al. (1994). The 111 clusters are all clusters from the catalogue that are not core-collapsed and include the 3 parameters we need for fitting.

Table 2: Retention fraction fit results for a constant kick velocity distribution.

We use two different models to describe the early globular clusters. The first is the King model fitted to present-day average GC properties (from the Harris catalogue). All heavy stars are assumed to be in the center of the cluster at the end of their lifetimes and the retention fractions are calculated using the central escape velocity. The second model assumes that no mass segregation has taken place. Stars of all masses are spread equally through the cluster according to the cluster's mass density profile. The Plummer model is used because it puts the GCs with different central potentials on equal footing. GCs evolve to small core radii and large tidal radii as they age. This means that the present-day central potential has no meaning when applied to primordial GCs. Through dynamical friction, the high mass stars may fall to the center before their maximum age is reached or they may end somewhere between where they started and the center of the cluster. The escape velocity needed for the stellar remnant is estimated accordingly. In each cases, we simulate the results (i.e. retention or ejection) for 10 000 supernova explosions producing neutron stars.

We use the equations from Drukier (1996) to calculate retention fractions. The probability that a star is retained to the cluster after applying the supernova "kick'' is as follows:

$\begin{displaymath} P({\rm ret}\vert\hat{v}, \hat{v_{k}}) = \left\{ \begin{arr... ...t{v};\\ 0 & \hat{v_{k}} -1 \geq \hat{v} \end{array} \right\}, \end{displaymath}$

(5)

$\begin{figure} \par\renewcommand{\epsfsize}[2]{0.7 ...$	Figure 5: The retention fraction as a function of mass, fitted to the King-model data with a constant kick velocity distribution between 0 and 1000 km s^-1. Chi-square was minimized to obtain the fit.
Open with DEXTER

As a simple way to test the method we considered a constant kick velocity distribution. We tried 10 progressively larger kick velocity intervals. The intervals we used are [0, x] km, with x = 500, 600, ..., 1400. We use the 3 different models: King, P-center (Plummer model, all heavy stars in the core) and P-spread (Plummer model with the stars spread through the GC, no primordial mass segregation) in most cases, although since the different models result in qualitatively similar conclusions, we do not exhaustively present all the results. To each we fit a power law for the retention fraction as a function of mass (Fig. 5). The plots for all other kick velocity distributions are qualitatively similar, so for the remainder of the paper, we present tables showing the results rather than plots. We obtain the exponent and the average retention fraction. The results are in Table 2. The exponents are all approximately 0.50. This is because the escape velocity goes as the square-root of the cluster mass. The King model has a slightly higher exponent because more massive clusters tend to favor a slightly higher central potential and are more concentrated, resulting in a somewhat higher escape velocity. The average retention fraction is linear with x. This is as expected. If a constant probability distribution is spread out over an interval twice as large (0-1000 instead of 0-500), one expects half as many neutron stars to be retained.

Table 3: Retention fraction fit results for a linear kick velocity distribution with negative slope.

We have made similar studies of the effects of linear kick distributions (i.e. where the probability distribution varies linearly, with a maximum at v=0, and a probability of zero above some characteristic velocity) and Gaussian velocity kick distributions (which are taken to be Gaussian probability density functions with the characteristic velocity equal to the spread in the velocity distribution). The linear distributions are characterized in terms of a characteristic velocity x, and are presented for large and small velocities in Tables 3 and 4 respectively. The linear distributions, like the constant probability density distributions, generically produce a M^0.5 retention fraction mass dependence, unless the characteristic velocity is less than the escape velocity of the most massive clusters, in which case the dependence will be flatter. Both these functional forms are thus consistent with the observations of extragalactic X-ray binaries, but will most likely have trouble explaining the dearth of very low velocity pulsars, and simultaneously explaining the small, but non-neglible numbers of very high velocity pulsars (e.g. those with velocities of more than 800 km s^-1) while still producing retention fractions greater than $\sim$ 1%, which is needed to explain the large numbers of X-ray binaries and millisecond pulsars in globular clusters.

Table 4: Retention fraction fit results for a linear kick velocity distribution with very large negative slope, for a King model.

Table 5: Retention fraction fit results for a Gaussian kick velocity distribution. The exponent values with asterisks denote cases where the retention fractions are quite small and the estimations suffer from small number statistics.

Gaussian distributions produce a wide range of mass dependences, since for a non-zero Gaussian, there typically will exist a positive slope in the probability density distribution over the range of globular cluster escape velocities. The results of Gaussian distribution studies are shown in Table 5. It is clear that to reproduce the observations, a Gaussian distribution would need a characteristic velocity just less than 20 km s^-1. This is, on its face, strongly inconsistent with the proposed existing kick velocity distributions, but, upon consideration of the effects of binaries, it may become consistent with e.g. the results of Arzoumanian et al. (2002). The presence of a heavy binary companion at the time of the supernova explosion will dilute the velocity kick applied to the system, making retention more likely; this has already been suggested to be a big part of the solution to the retention fraction problem (e.g. Davies & Hansen 1998; Pfahl et al. 2002). The maximum possible dilution would occur if the envelope of the proto-neutron star has been completely stripped in the binary, while the companion star is very close to the maximum possible mass in order not to undergo a supernova explosion itself later. That is to say, the maximum dilution occurs for a supernova explosion of a 1.4 $M_\odot$ core with an 8 $M_\odot$ companion, and this dilution will be a factor of 6.7 [i.e. 1.4/(1.4+8.0)]. This would be enough to bring a 90 km s^-1 kick (the characteristic kick velocity of the slow mode found by Arzoumanian et al. 2002) down to about 13 km s^-1 (we do note that the Arzoumanian model is a Maxwellian, rather than a Gaussian, and a 90 km s^-1 Maxwellian will have slightly few low velocity objects than a 90 km s^-1 Gaussian, but that this should lead to only a minor difference). Thus we can say that the most extreme dilution factor is more than sufficient to produce a reasonable dependence of the retention fraction on mass. Additionally, heavy binaries will have lower initial velocities than the stars in the cluster as a whole, which should allow them to accomodate larger kicks (although in no case will this increase the kick velocity allowed by more than $\sqrt{2}$ ). Since mass segregation will bring the most massive stars to the centers of clusters, and these massive stars will have the highest cross-sections for interactions, it is likely that supernova progenitors will find themselves in binaries with other massive stars. Whether this actually happens to the extent required by the data must be addressed through numerical simulations and is beyond the scope of this paper.

Regardless, it seems that in no case will the data be consistent with a single Maxwellian or Gaussian with a post-dilution characteristic velocity greater than 20 km s^-1. This seems to rule out the Hobbs et al. (2005) suggestion that the pulsar velocity kick distribution can be a single Maxwellian at 265 km s^-1, unless the slower kicks in globular clusters are related to binary evolution (see e.g. Pfahl et al. 2002b; Dewi et al. 2005, for suggestions that binary evolution may affect pulsar kick velocities). It also indicates that the contribution of neutron stars from the fast mode in a distribution like that of Arzoumanian et al. (2002) to retention fractions should be negligible. These results are thus quite similar to the conclusions of Pfahl et al. (2002) who found that, while a single fast kick mode had trouble producing neutron star retention fractions large enough to match the observed numbers of millisecond pulsars in globular clusters like 47 Tuc, including also a slow velocity kick mode, and the effects of binaries could potentially solve the retention problem. With careful numerical simulations, the globular cluster X-ray binaries are likely to present the strongest constraints on the low velocity end of the pulsar kick velocity distribution, since field studies will always be complicated, e.g., by the intrinsic velocity dispersion in the field.

To obtain strong constraints on the actual asymmetric kick velocity applied to neutron stars at the times of their birth will require numerical simulations which can take account of both the distributions of the neutron star progenitors in the cluster at the time the supernovae occur, and what fraction of the neutron stars are in binaries, and what are the masses of those binary companions. The numerical simulations will also need to be made for ranges of cluster masses and initial concentration indices. At the present, the best studies of neutron star retention (e.g. Drukier 1996; Davies & Hansen 1998; Pfahl et al. 2002) have each considered several aspects of this problem; the simulation closest to meeting these requirements was that of Pfahl et al. (2002), in which all the relevant factors except for a range in cluster escape velocities were considered.

The globular cluster mass/low mass X-ray binary correlation: implications for kick velocity distributions from supernovae