DYNAMICAL FORMATION OF MILLISECOND PULSARS IN GLOBULAR CLUSTERS

It is well known that the formation rate per unit mass of low-mass X-ray binaries (LMXBs) is orders of magnitude greater in globular clusters (GCs) than in the Galactic field (Katz 1975; Clark 1975). The high formation rate of LMXBs in GCs is attributed to the frequent dynamical interactions in the dense stellar environment. As an additional formation channel, the binaries in a GC can also be a result of a standard evolutionary path identified for millisecond pulsar (MSP) formation in the Galactic field. This has stimulated many theoretical and observational studies to investigate the relative contribution of these two formation processes of compact binaries in the population of GCs (e.g., Fregeau 2008; Pooley et al. 2003; Pooley & Hut 2006).

With the superior sub-arcsecond spatial resolution of the Chandra X-Ray Observatory, remarkable progress has been made in the understanding of the formation processes of close binaries in GCs. For example, Pooley et al. (2003) found a positive correlation between the number of close X-ray binaries in GCs and the stellar encounter rate, Γ_c. Specifically, Pooley et al. (2003) found an approximately linear relationship between the number of LMXBs and Γ_c, indicating a dependence on the properties of GCs. A similar relationship has also been reported by Gendre et al. (2003), and, taken together with the results of Pooley et al. (2003), these studies provide evidence for the dynamical origin of LMXBs in GCs.

Since MSPs have long been proposed as the descendants of LMXBs, they are also expected to have a dynamical origin in GCs. Due to the existence of extensive pulsar surveys, 140 MSPs have been detected in 26 different clusters and a statistical study of their relationship to cluster parameters is desirable.⁶ However, previous studies were not successful in finding evidence for the dynamical origin of MSPs in the clusters due to the lack of a relation between the pulsar population and Γ_c in the GCs (e.g., Ransom 2008). This can be ascribed to the observational bias in the pulsar searches. As the distance of the GCs spans a rather wide range (cf. Harris 1996)⁷, the sensitivity of the observations can differ in the searches toward different clusters and hence induce selection effects in the observed sample (Ransom 2008). Therefore, the observed number of MSPs is not representative of an unbiased sample for the analysis.

In this paper, we present a method to alleviate the aforementioned problem and investigate the possible relationship between the number of MSPs and the cluster properties. In Section 2, an investigation of the cumulative luminosity distribution functions (CLFs) of MSPs in a number of selected GCs is carried out. We subsequently use the obtained results in a correlation analysis in Section 3 and discuss the physical implications of the possible correlation in Section 4.

Among the GCs, nine contain more than three MSPs for which their radio flux densities are reported. This enables one to create the CLFs of MSPs for these clusters and to estimate their logarithmic slopes by a power-law fit. The physical properties of these nine selected GCs are summarized in Table 1. As the observations for an individual cluster were conducted at different frequencies, we convert all flux densities to the values at 1.4 GHz by assuming a spectral index of −1.8, which is the mean value reported by Maron et al. (2000).

We model the CLFs with a form of N(>L_1.4 GHz) = N₀L^q_1.4 GHz, where N₀ represents the number of MSPs with the pseudo-radio luminosity >1 mJy kpc² at 1.4 GHz. Since no obvious turn-off is observed from the distribution of individual cluster, we have taken the entire sample into account. Following Hessels et al. (2007), we adopt the square root of N(>L_1.4 GHz) as the uncertainties in the analysis and fit log N(>L_1.4 GHz) versus log L_1.4 GHz with a linear regression analysis that take the uncertainties into account. The best-fit parameters are tabulated in Table 2. It can be seen that the best-fit values of q lie in a range ∼−0.6 to ∼−1.6 with the steepest slope inferred for the MSPs in M 3. In this analysis, we have also taken the uncertainties associated with the cluster distance determination into account. In reviewing different methods in determining the distances to GCs, it has been concluded that these various methods imply an uncertainty of ±6% in the cluster distance (Krauss & Chaboyer 2003; Chaboyer 2008). Adopting the uncertainty of distance in constructing the CLFs, this results in an additional error in the normalization beyond the statistical error from the fits.⁸ Both errors are combined in quadrature which are quoted in Table 2 and adopted in subsequent analysis in Section 3.

By using the best-fit CLFs for extrapolation, we estimate the number of MSPs in each cluster with L_1.4 GHz>0.5 mJy kpc². Almost all the MSPs in GCs considered in this study have their radio luminosities above this threshold level. Since the errors of both q and N₀ are considered in the extrapolation, the uncertainties of these estimates are larger than that of N₀. These population estimates are also given in Table 2. For those clusters which have been searched deep enough that the sensitivity is close to this threshold, e.g., 47 Tuc, the population estimates based on the CLFs are close to the observed number of uncovered MSPs. On the other hand, for those clusters where pulsation searches have not yet reached this sensitivity level, our results provide predictions for their MSP populations when searches toward these clusters become sufficiently deep.

Apart from estimating the CLFs for these individual GCs, we have also constructed the luminosity function by combining all the cluster MSPs (i.e., 76 pulsars in total) used in this study. The combined CLFs of the selected cluster MSP population are displayed in Figure 1. In the combined distribution, we have observed that there is a turn-off for the pseudo-luminosities smaller than ∼1.5 mJy kpc². Hessels et al. (2007) have also observed the same behavior in an independent analysis. In order to compare our result with Hessels et al. (2007), we followed their procedure and used a minimum luminosity cutoff of L_1.4 GHz = 1.5 mJy kpc² for the fitting procedure. The results of this analysis are summarized in Table 3. The logarithmic slope inferred from this collective sample is q = −0.83 ± 0.05, which is consistent with the value deduced by Hessels et al. (2007; i.e., q = −0.77 ± 0.03) within 1σ error.

Zoom In Zoom Out Reset image size

Hessels et al. (2007) have further compared the luminosity distribution of the MSP population in GCs with that of other pulsar populations. They suggested that the distribution of cluster MSPs is marginally consistent with that of the MSPs in the Galactic field reported by Cordes & Chernoff (1997) and Lyne et al. (1998). However, the samples adopted in these investigations contain not more than 22 pulsars. This relative small sample size may introduce bias in the analysis and inference. Thanks to the extensive pulsar searches in the recent years, the whole pulsar population has been significantly increased. Hessels et al. (2007) have also compared their results with a more updated pulsar sample reported by Lorimer et al. (2006) which used 1008 pulsars in their study. Lorimer et al. (2006) have found a distribution of dlog N ∼ −0.8dlog L for their sample, which is very similar to that inferred from the cluster population as reported by Hessels et al. (2007) and this paper. Nevertheless, the sample used by Lorimer et al. (2006) consists of non-recycled canonical pulsars. It is more instructive to compare the CLF of the cluster MSP population and the MSPs in the Galactic field, both of which have undergone the recycling processes in binary systems. Therefore, we construct the luminosity function for the MSP population in the Galactic field with all the available data in the ATNF catalog (Manchester et al. 2005)⁹.

To do this, we specifically selected all pulsars in the Galactic field with P< 20 ms. In total, there are 51 MSPs with flux densities reported in the catalog, which enable us to construct the CLF. For comparison, the CLF of the Galactic population is over-plotted along with the cluster population in Figure 1. To be consistent with the analysis of the cluster population, we used a minimum luminosity cutoff of L_1.4 GHz = 1.5 mJy kpc² in the fitting. The inferred slope of the Galactic population is q = −0.48 ± 0.04, which is found to be flatter than the value deduced from the cluster population.

To further investigate the difference between two populations, we separate each population into two sub-samples, namely the binary and the solitary MSPs. The results are summarized in Table 3.¹⁰ By studying these sub-populations in GCs separately, we deduce the slopes to be q = −0.73 ± 0.08 and q = −0.89 ± 0.11 for the binary and the solitary MSPs, respectively. These values are consistent with those inferred by the Hessels et al. (2007) within 1σ uncertainties. In the Galactic field, we also examine the luminosity functions separately for the binary and solitary MSPs. The slopes inferred for the binary and solitary populations are q = −0.49 ± 0.05 and q = −0.24 ± 0.11, respectively, which are rather different from those inferred from the cluster population. For completeness, we have also repeated the above analysis with all the MSPs (i.e., without excluding those with their pseudo-luminosities smaller than 1.5 mJy kpc²). The results are also tabulated in Table 3 for the sake of comparison.

We have attempted to search for correlations between the observed number of MSPs in each GC with cluster parameters. As the most promising correlation is expected with the two-body encounter rate, we begin our analysis with this parameter. The two-body encounter rate in a cluster can be estimated as Γ_c ∝ ρ²₀r³_cσ⁻¹₀, where ρ₀ is the central luminosity density, r_c is the core radius, and σ₀ is the velocity dispersion at the cluster center. In adopting the central luminosity density as an estimate of the stellar density at the center, there is an underlying assumption that the average luminosities of the stars in the cluster centers are approximately equal to ∼1 L_☉. Figure 2 shows the relation between the observed populations with Γ_c. However, no obvious correlations can be identified with this observed sample. These two quantities are correlated at a confidence level of ∼75% only. This conclusion is similar to that reported by Ransom (2008). Owing to the limited amount of telescope time, many clusters that host only a single pulsar have not been searched to the same sensitivity level as those of the specifically selected targets, such as 47 Tuc. Therefore, the selection effect biases the observed numbers of MSPs in different clusters.

Zoom In Zoom Out Reset image size

In order to alleviate this problem, we suggest the use of the CLFs of the investigated clusters. With the best fits of the CLFs (see Section 2), we are able to estimate the number of MSP in these GCs above a given luminosity threshold and thus obtain an unbiased sample. Specifically, we take the best-fit values of N₀ in these GCs to estimate the numbers of the MSPs in these clusters with their pseudo-luminosities above >1 mJy kpc² and examine whether it is related to different physical quantities of the clusters. In this analysis, the possible correlation between N₀ with two-body encounter rate Γ_c, metallicity [Fe/H], cluster mass M_GC, velocity dispersion σ₀, and escape velocity v_escape at the cluster center are explored. All these quantities are speculated to have influence on the binary formation and hence the MSP population in a cluster.

While the two-body encounter rate Γ_c is related to the binary population resulting from dynamical interactions, the metallicity [Fe/H] of a cluster can have a profound influence on the evolution of LMXBs (see Ivanova 2006 and the discussion in Section 4). On the other hand, if stellar encounters were not the major channel of the binary formation, one would expect the binary population to be correlated with the cluster mass M_GC. Assuming a constant mass-to-light ratio, M_GC can be estimated from the absolute visual magnitude M_V: $M_{\rm GC}\propto 10^{-0.4M_{V}}$. We have also tested the correlation with σ₀ and v_escape which may possibly be related to the retention of the neutron stars in a cluster.

Without a priori knowledge of the distributions of the tested quantities, a nonparametric correlation analysis is adopted. The computed Spearman rank correlation coefficients between N₀ and the various quantities are tabulated in Table 4. Among all the tested quantities, the strongest correlation is found between N₀ and Γ_c. The corresponding Spearman correlation is 0.78 with a chance correlation probability of 0.0125. The plot of N₀–Γ_c is displayed in Figure 3(a). The correlation between N₀ and [Fe/H] with a Spearman correlation = 0.72 has also been found to be significant with a chance correlation probability = 0.0298, which is plotted in Figure 3(b). By taking the errors of N₀ as the weight in the linear regression analysis, the logarithmic slopes of the N₀–Γ_c and N₀–[Fe/H] relations are found to be 0.69 ± 0.11 and 0.72 ± 0.11, respectively. For the other tested quantities, there are marginal correlations of N₀ versus v_escape and σ₀ at a confidence level ≳89%, though it is not sufficiently significant to secure the relations. It is not surprising to note that the rank correlation coefficients are the same for these two quantities, as Gnedin et al. (2002) have found that the ratio of v_escape to σ₀ has a narrow range between ∼3 and 5. Among all the tested quantities, the weakest correlation is found for the N₀–M_GC relation which has a chance correlation probability over 60%.

Zoom In Zoom Out Reset image size

Table 4. Spearman Rank Correlation Coefficient of N₀ Versus Various Physical Properties of GCs

Parameter	Spearman Rank Correlation Coefficient	Probabilitya
With N(L1.4 GHz>1 mJy kpc2)
Γc ∝ ρ20r3cσ−10	0.7833	0.9875
[Fe/H]	0.7167	0.9702
$M_{{\rm GC}}\propto 10^{-0.4M_{V}}$	−0.2000	0.3941
σ0	0.5667	0.8884
vescape	0.5667	0.8884
With N(L1.4 GHz>0.5 mJy kpc2)
Γc ∝ ρ20r3cσ−10	0.8333	0.9947
[Fe/H]	0.6167	0.9231
$M_{{\rm GC}}\propto 10^{-0.4M_{V}}$	−0.1333	0.2676
σ0	0.5167	0.8456
vescape	0.5167	0.8456
With N(L1.4 GHz>2 mJy kpc2)
Γc ∝ ρ20r3cσ−10	0.7628	0.9832
[Fe/H]	0.7798	0.9868
$M_{{\rm GC}}\propto 10^{-0.4M_{V}}$	−0.2034	0.4003
σ0	0.6272	0.9294
vescape	0.6272	0.9294

Note. ^aThe probability that the correlation coefficient is different from zero.

Download table as:  ASCIITypeset image

As this choice of luminosity threshold is arbitrary, we further check the robustness of the correlation analysis results by repeating the investigation with different thresholds. We have repeated the analysis by adopting N(L_1.4 GHz>0.5) in Table 2, which provides the estimates for the number of MSPs with L_1.4 GHz>0.5 mJy kpc². Almost all the MSPs in GCs considered in this investigation have their radio luminosities above this threshold. With these new values, the correlations of the MSP number versus Γ_c, [Fe/H], M_GC, and v_escape (or σ₀) are found at the confidence levels of 99.47%, 92.31%, 26.76%, and 84.56%, respectively.

For a further test of the robustness, by using the best-fit CLFs in Table 2, we have also repeated the analysis for a minimum luminosity cutoff of L_1.4 GHz>2 mJy kpc². In this case, the correlations with Γ_c, [Fe/H], M_GC, and v_escape (or σ₀) are confident at the levels of 98.32%, 98.68%, 40.03%, and 92.94%, respectively. Therefore, the degrees of correlation for the tested quantities are found to be insensitive to the choice of the threshold. We conclude that the correlation between Γ_c and the MSP number is the most robust among all the tested quantities, which have a confidence level >98% regardless of the chosen threshold.

The CLFs of nine GCs, each containing a population of MSPs, have been examined. Upon comparison of the MSP population in GCs with that in the Galactic field, it has been found that the slopes of the CLFs inferred in these two populations significantly differ. It is natural to speculate that the CLF is somehow related to the magnetic field and spin of the MSPs. Wang et al. (2005) have compared the distributions of the spin period and the dipolar surface magnetic field for both cluster and disk populations (cf. Figures 2 and 3 in their paper). Despite the broader distribution for the disk population, their mean values are not dissimilar in both populations and therefore cannot solely explain the difference of CLFs.

Apart from the radio luminosity functions, the X-ray emission properties of the MSPs in the GCs are also found to be very different from those in the Galactic field. While the MSPs in the Galactic field generally require a hot polar cap component plus a nonthermal power-law tail to model their X-ray spectra (cf. Zavlin 2006), the X-rays from a majority of the MSPs in GCs are purely thermal in nature (see Hui et al. 2009 and the references therein for a recent review). Cheng & Taam (2003) suggest the absence of nonthermal X-ray from the cluster MSPs can be possibly related to the complicated magnetic field structure. Since the stellar interaction in GCs is much more frequent than that in the Galactic field, MSPs in the GCs can possibly change their companion several times throughout their lives. Since the orientation of the binary after each exchange can differ, the direction of the angular momentum accreted during the mass transfer phase subsequent to each exchange can vary possibly affecting the magnetic field configuration at the neutron star surface. Such an evolution could lead to a much more complicated magnetic field structure for the MSPs in the GCs than in the case of the Galactic field. In such a complicated magnetic field, Ruderman & Cheng (1988) have argued that high-energy curvature photons will be emitted and subsequently converted into pairs to quench the accelerating region. This provides an explanation for the absence of nonthermal emission in the cluster MSPs. For the same reason, the complicated magnetic field structure can also possibly alter the coherent radio emission and result in a different radio luminosity of the cluster MSPs in comparison with the disk population.

Adopting the best-fit normalization inferred from the CLFs of individual cluster as an unbiased estimate of the number of MSPs, we have further examined the relationships between the pulsar population and the physical properties in GCs. We have found the positive correlations of N₀ versus Γ_c as well as N₀ versus [Fe/H] at a relatively high confidence level. A marginal positive correlation between N₀ and v_escape is also suggested. Although a high escape speed implies the presence of a deeper gravitational potential well and hence a higher neutron star retention, this correlation is not sufficiently significant to warrant such an interpretation. Hence, we do not discuss this relation any further and focus on the physical implications of the N₀–Γ_c and N₀–[Fe/H] relations.

Due to the different selection effects in the pulsar search surveys, it is not feasible to directly use the detected MSP populations in GCs for a statistical analysis. Instead, we alleviate the problem by taking N₀ as the estimator for the number of pulsars with pseudo-radio luminosities at 1.4 GHz larger than 1 mJy kpc². With this consideration, we have found a correlation between N₀ and Γ_c at a confidence level >98%. We have further found that the strength of this correlation is robust and independent of the choice of the luminosity cutoff by repeating the analysis with different thresholds. This provides evidence for the dynamical formation of MSPs in GCs. For a competing scenario that the MSPs have a binary origin similar to the Galactic field, one should expect the number of MSPs to scale with the cluster mass, M_GC, instead of Γ_c. However, we do not find any convincing relationship between N₀ and M_GC (see Table 4). The absence of correlation with M_GC provides additional support for the dynamical formation scenario. Taken together with the difference in the X-ray luminosity functions of LMXBs in the field and in GCs (see Voss et al. 2009; Kim et al. 2009), it is likely that the MSPs have different origins/evolutions in GCs relative to the Galactic field.

We note that the logarithmic slope of the power-law fit in the N₀–Γ_c relationship (i.e., 0.69 ± 0.11) is not dissimilar to that of the number of X-ray sources versus Γ_c (0.74 ± 0.36; Pooley et al. 2003). This dependence on the two-body encounter rate suggests a possible relationship between the MSP population and close X-ray binaries in GCs. Apart from the whole X-ray binary population, Pooley et al. (2003) and Gendre et al. (2003) have also examined the relationship for the individual class of LMXBs which has a logarithmic slope of 0.97 ± 0.5. Although the large uncertainty of this slope resulting from the limited sample of LMXBs precludes a definitive conclusion concerning the link between LMXBs and MSPs, it is consistent with such an interpretation.

Theoretical arguments (Verbunt & Hut 1987) suggest that the number of LMXBs is linearly proportional to the stellar encounter rate of the cluster; however, direct comparison of their relationship with the current two-body encounter rate may be misleading. As the MSPs are long-lived and are produced by the previous generations of LMXBs, they can have a different formation rate from the LMXB population currently observed. This point is important since the relaxation time at the cluster core is generally longer than the lifetime of LMXBs (cf. Harris 1996)⁷. Therefore, the continuous mass segregation at the cluster center can result in an evolution of the stellar collision frequency and hence a varying formation rate of compact binary systems. Nevertheless, the combination of X-ray and Hubble Space Telescope observations of Cen A (see Jordán et al. 2007) indicates that GCs with LMXBs are characterized by higher stellar encounter rates than those devoid of LMXBs.

In addition to the N₀–Γ_c relation, we have also found a positive correlation between N₀ and the metallicity of the GCs. It has been noted that observational evidence suggests that bright LMXBs are preferably formed in metal-rich clusters in our Milky Way as well as other galaxies (e.g., Bellazzini et al. 1995; Maccarone et al. 2004; Jordán et al. 2004). Ivanova (2006) proposes that the absence of the outer convective zone in metal-poor main-sequence donor stars in the mass range of 0.85–1.25 M_☉ in comparison to their metal-rich counterparts can be responsible, since the absence of magnetic braking in such stars precludes orbital shrinkage, thereby, significantly reducing the binary parameter space for the production of bright LMXBs. For the conventional scenario that LMXBs are the progenitors of MSPs, the positive correlation between N₀ and [Fe/H] is not unexpected, since the MSP number should scale with that of their progenitors.

While the stellar encounter rate has been widely accepted as a parameter to indicate which clusters are likely to host a large MSP population, our study suggests that the metallicity can also be an important parameter. To explore this hypothesis, we suggest that pulsar searches be carried out toward metal-rich GCs, such as Liller 1 which has the highest metallicity ([Fe/H] = 0.22) among all 150 GCs in the Milky Way (cf. Harris 1996). Furthermore, its two-body encounter rate is estimated to be comparable with that of 47 Tuc. Therefore, according to these parameters, it is very likely to host a considerable number of MSPs. With a dedicated search, this hidden population may be revealed.

K.S.C. was supported by a GRF grant of Hong Kong Government under HKU700908P. R.E.T. was supported in part by NSF grant AST-0703950 at Northwestern University and by the Theoretical Institute for Advanced Research in Astrophysics (TIARA) operated under the Academia Sinica Institute of Astronomy and Astrophysics in Taipei, Taiwan.