Millisecond Pulsars and Black Holes in Globular Clusters

We use our Cluster Monte Carlo code (CMC) to simulate a grid of models for this study. CMC is a parallelized Hénon-type Monte Carlo code (Hénon 1971a, 1971b) that has been developed and rigorously tested over many years (Joshi et al. 2000, 2001; Fregeau et al. 2003; Fregeau & Rasio 2007; Chatterjee et al. 2010, 2013a; Umbreit et al. 2012; Pattabiraman & Umbreit 2013; Rodriguez et al. 2018). CMC incorporates all the relevant physics for GC evolution, including two-body relaxation, three-body binary formation, strong three- and four-body interactions, and some post-Newtonian effects (Rodriguez et al. 2018). Updated versions of the SSE and BSE packages (Hurley et al. 2000, 2002) are used for single and binary stellar evolution in CMC and Fewbody (Fregeau et al. 2004; Fregeau & Rasio 2007) is used to directly integrate all three- and four-body gravitational encounters, with post-Newtonian effects (Antognini et al. 2014; Amaro-Seoane & Chen 2016).

In this study, we consider a set of 26 independent cluster models. In models 1–25 (which serve as our main grid of models with present-day properties typical of Milky Way GCs), we fix a number of initial cluster parameters: total star number N = 8 × 10⁵, binary fraction f_b = 5%, virial radius r_v = 1 pc, King concentration parameter W_o = 5, galactocentric distance r_g = 8 kpc, and metallicity Z = 0.001 (Morscher et al. 2015; Chatterjee et al. 2017). We use the initial mass function given in Kroupa (2001) ranging from 0.08 to 150 M_⊙ to sample the initial stellar masses. NS remnants from CCSNe at formation receive natal kicks drawn from a Maxwellian distribution with a standard deviation σ_NS = 265 km s⁻¹ (Hobbs et al. 2005).

Following Kremer et al. (2018), the only parameter varied in our models is the natal kick for BHs, which allows us to easily isolate and understand the effects of BH retention on the long-term GC evolution and the dynamical evolution of MSPs. Varying the BH natal kicks to alter the retained number of BHs at present times has a similar effect on the GC properties as varying more physically motivated initial conditions, such as the initial virial radius (Kremer et al. 2019). With all other parameters being fixed, all differences between models clearly originate from the differences in the fraction of BHs retained in the cluster upon formation (for more details see Kremer et al. 2018). All models were evolved for 12 Gyr.

The natal kicks for BHs are assumed to be independent of BH masses and are first drawn from the same Maxwellian distribution as the NSs. Their magnitudes are then multiplied by the ratio $\tfrac{{\sigma }_{\mathrm{BH}}}{{\sigma }_{\mathrm{NS}}}$, which is varied between models. We set different ratios from 0.005 to 1.0 for different models as shown in Table 1. The introduction of this simple parameter controlling the natal kicks received by the BHs gives us a key benefit for the purpose of this initial study: it allows us to control the BH retention fractions of our models without the needing to change any other initial parameters. As a result, the changes between the properties of these models can unambiguously be attributed to the difference in BH retention fractions.

Table 1.  Cluster Model Properties

Model	N	$\tfrac{{\sigma }_{\mathrm{BH}}}{{\sigma }_{\mathrm{NS}}}$	rc	rhl	MTOT	NBH	NNS	NNS−DYN	NPSR	NMSP	NsMSP	NbMSP	NS − BH	DNS
	105		pc		105 M⊙								9 Gyr < t < 12 Gyr
1	8	0.005	2.88	4.00	2.10	383	444	15	1	1	0	1	0	0
2	8	0.01	1.85	4.45	2.12	366	405	12	2	2	1	1	0	0
3	8	0.02	2.72	4.04	2.12	332	414	20	1	1	0	1	0	0
4	8	0.03	1.85	3.92	2.14	339	431	14	1	1	0	1	0	0
5	8	0.04	1.54	4.17	2.15	328	396	11	0	0	0	0	0	0
6	8	0.05	2.31	4.39	2.11	338	424	20	2	2	0	2	0	0
7	8	0.06	2.04	3.49	2.15	293	446	14	1	1	0	1	0	0
8	8	0.07	1.53	3.73	2.18	277	437	12	1	1	0	1	0	0
9	8	0.08	1.44	3.18	2.20	237	414	14	1	1	0	1	0	0
10	8	0.09	1.53	3.10	2.22	210	431	13	1	1	0	1	0	0
11	8	0.1	1.20	3.02	2.24	189	489	21	1	1	0	1	1	0
12	8	0.11	0.68	2.60	2.28	136	445	18	2	2	0	2	0	0
13	8	0.12	0.77	2.31	2.29	114	452	20	1	0	0	0	0	0
14	8	0.13	0.51	1.79	2.31	90	440	21	1	1	0	1	0	0
15	8	0.14	0.54	1.83	2.31	66	464	18	0	0	0	0	0	0
16	8	0.15	0.48	1.82	2.32	67	473	20	1	1	0	1	0	0
17	8	0.16	0.18	1.53	2.33	39	462	14	2	2	1	1	2	0
18	8	0.17	0.36	1.44	2.35	31	436	19	5	4	0	4	0	0
19	8	0.18	0.24	1.51	2.35	21	443	27	3	1	0	1	0	0
20	8	0.19	0.14	1.43	2.29	11	453	32	1	0	0	0	1	0
21	8	0.2	0.25	1.48	2.32	5	493	48	4	1	0	1	3	2
22	8	0.4	0.20	1.95	2.29	1	505	178	8	7	1	6	2	10
23	8	0.6	0.34	2.07	2.34	1	481	173	9	8	3	5	0	7
24	8	0.8	0.22	2.05	2.31	1	522	187	17	13	1	12	7	15
25	8	1.0	0.24	2.03	2.28	0	524	184	10	7	2	5	0	4
26	35	1.0	0.18	1.41	11.18	20	3908	1211	132	83	14	69	24	55

Note. Columns 1–3: model number, initial number of stars, and BH natal kick scaling factor, respectively. Columns 4–13 show final model properties at 12 Gyr: projected core radius, projected half-light radius, total cluster mass, number of BHs, number of NSs, number of dynamically formed NSs, total number of pulsars including MSPs, total number of MSPs, number of single MSPs, and number of binary MSPs. The last two columns give the number of NS–BH binaries and DNSs that appear at any time in the model between 9 and 12 Gyr. N_NS−DYN is the number of NSs formed, or strongly affected, by dynamical processes, and includes all NSs formed through MIC or AIC, as well as all NSs that experienced a direct collision or merger with another star and had their properties reset (as explained in Section 2.2.1).

Download table as:  ASCIITypeset image

In addition, we simulate one model (model 26 in Table 1) meant to represent a very massive GC. This model has initial N = 3.5 × 10⁶, a binary fraction of 10%, and a BH natal kick $\tfrac{{\sigma }_{\mathrm{BH}}}{{\sigma }_{\mathrm{NS}}}=1.0$. Other initial conditions for this model are the same as in models 1–25. The final mass of the model is about 1.2 × 10⁶ M_⊙, which is close to the masses of 47 Tuc and Terzan 5. This model is an extreme case designed to form many MSPs as a result of the large initial N, higher binary fraction, and large BH kicks.

Throughout this paper, we refer to models with final numbers of retained BHs greater than 200 as "BH-rich" models, while those with <10 BHs at 12 Gyr are called "BH-poor" models; others in between these limits are referred to as "BH-intermediate" models.

For this work we have updated the prescriptions for the evolution of NSs from SSE and BSE in CMC. The version of BSE used in CMC is now nearly identical to COSMIC 2.0.0 (Breivik et al. 2019). These updates include changes to the magnetic field and spin-period evolution for single and binary NSs, and to the natal kick prescriptions for NSs formed in ECSNe (Kiel et al. 2008; Kiel & Hurley 2009).

2.2.1. Magnetic Field and Spin-period Evolution

The evolution of NS magnetic fields and spin periods has long been a topic of debate and still remains uncertain (e.g., Faucher-Giguere & Kaspi 2006, and references therein). In our models, NS remnants, when formed, are assigned randomly sampled magnetic fields (range 10^11.5–10^13.8 G) and spin periods (range 30–1000 ms) to match observations of young pulsars. To model NS evolution, we follow the prescriptions described in Hurley et al. (2002) and Kiel et al. (2008).

As outlined in Kiel et al. (2008), we assume that the dominant spin-period evolution mechanism for single NSs is dipole radiation, and NSs are treated as solid spheres. The spin-down rate of single NSs is

$$\begin{eqnarray}&&\dot{P}=K\displaystyle \frac{{B}^{2}}{P},\end{eqnarray} \tag{ 1 }$$

where K = 9.87 × 10⁻⁴⁸ yr G⁻², B is the surface magnetic field, and P is the spin period.⁷

Additionally, magnetic fields of single NSs are assumed to decay exponentially,

$$\begin{eqnarray}&&B={B}_{0}\exp \left(-\displaystyle \frac{T}{\tau }\right),\end{eqnarray} \tag{ 2 }$$

on a timescale τ = 3 Gyr (Kiel et al. 2008). Here B₀ and T are the initial magnetic field and the NS's age, respectively. Faucher-Giguere & Kaspi (2006) showed that there is no significant magnetic field decay for single radio pulsars on a timescale of ∼100 Myr. Recent magneto-thermal models of NS magnetic field also show that the magnetic field evolution of single radio pulsars is compatible with no decay or weak decay (e.g., Popov et al. 2010; Viganò et al. 2013). The timescale we adopt here for magnetic field evolution of single radio pulsars is compatible with a very weak field decay. Further exploration of the effects of magnetic field decay on GC pulsars will be presented in future works.

For NSs in binaries, binary evolution is also taken into account. The evolution of NSs in detached binaries is the same as that for single NSs. On the other hand, during mass-transfer episodes, the magnetic fields and spin periods of NSs can change significantly on a short timescale. During accretion, the magnetic field is assumed to decay as

$$\begin{eqnarray}&&B=\displaystyle \frac{{B}_{0}}{1+\tfrac{{\rm{\Delta }}M}{{10}^{-6}{M}_{\odot }}}\exp \left(-\displaystyle \frac{T-{t}_{\mathrm{acc}}}{\tau }\right),\end{eqnarray} \tag{ 3 }$$

where t_acc is the duration of the NS accretion phase (during RLOF) and ΔM is the mass accreted.⁸ The NS spin period decreases accordingly through angular momentum transfer (see Hurley et al. 2002, Equation (54)). Wind mass loss, tidal evolution, magnetic braking, and supernova kicks have only small effects on the magnetic field and spin period evolution. Only stable mass transfer in a binary system can spin up NSs to MSPs in our models; we ignore the possibility of mass accretion by NSs during common-envelope phases (Hurley et al. 2002).

Occasionally, through dynamical or binary evolution, NSs merge with main-sequence (MS) stars, giants, or WDs. If the outcome of this merger is an NS (as opposed to a BH; see Hurley et al. 2002, Table 2), the magnetic field and spin period for this NS are reset by drawing from the same ranges of initial values as above. However, if an MSP is involved in such a merger, different initial magnetic fields and spin periods are assigned (relative to regular NSs) so that the new values still match those observed in MSPs. Specifically, for MSPs involved in mergers, a new magnetic field is drawn from the range 10⁸–10^8.8 G and a new spin period is drawn from the range 3–20 ms. In other words, we assume that an MSP involved in a merger remains an MSP after the merger, but allowing for a small random change in magnetic field and spin period.

In addition, for all MSPs we assume a lower limit for the NS magnetic fields to be 5 × 10⁷ G (Kiel et al. 2008). We do not set a lower limit for the spin periods. We use the standard expression

$$\begin{eqnarray}&&\displaystyle \frac{B}{{P}^{2}}=0.17\times {10}^{12}\,{\rm{G}}\ {{\rm{s}}}^{-2}\end{eqnarray} \tag{ 4 }$$

(Ruderman & Sutherlandt 1975; Bhattacharya et al. 1992) for the pulsar death line. Note that the exact location of the death line is still under debate (see, e.g., Zhang et al. 2000; Zhang 2002; Zhou et al. 2017). Throughout this paper we define a pulsar as any NS with a magnetic field and spin period above the death line, an MSP as any pulsar with P ≲ 30 ms, and a young pulsar as any pulsar that is not an MSP.

2.2.2. Electron-capture Supernovae

After their formation at early times, most of the NSs in our models evolve like young pulsars in isolation until the end (12 Gyr). Their magnetic fields and spin rates slowly decrease until they die as pulsars. Dynamical interactions are essential for producing MSPs throughout the evolutionary histories of GCs (e.g., Ivanova et al. 2008). Mass segregation puts NSs close to the cluster centers in high stellar density regions. During frequent stellar encounters, single NSs can acquire companions through exchanges, and wide binaries can be hardened by repeated interactions. NSs in binaries can then be spun up to MSPs through mass transfer.

Figure 1 illustrates two examples of the dynamical interaction histories of MSPs. The progenitor NS of an MSP on the right of the figure was formed in EIC at 0.11 Gyr. It had a magnetic field of 6 × 10¹² G and a spin period of 320 ms at birth. The pulsar experienced its first encounter at 5.3 Gyr and acquired a WD companion through a binary–single exchange. A second binary–single exchange replaced the WD companion to another WD. The binary then experienced a few binary–single and binary–binary scatterings. During this time, there was no mass transfer, and the magnetic field and spin rates of the pulsar were decreasing (solid black line in Figure 2). A third binary–single exchange encounter gave the pulsar an MS star companion. The MS star spun up the pulsar via RLOF-driven mass transfer (the first red-dashed line in Figure 2). At this time it was not yet an MSP (not all mass transfer leads to MSPs; only extended periods of stable mass transfer can produce MSPs). The MS star later evolved to become a giant and the binary underwent a common-envelope phase (blue line in Figure 2), which circularized and shrank the binary orbit. The remnant of the MS star was a WD that continued to fill its Roche lobe and spun up the pulsar to a MSP (the second dashed red line in Figure 2). The final system has an MSP with a 5.9 × 10⁷ G magnetic field and a 1.56 ms spin period, and a companion with a mass of 0.0075 M_⊙ at 12 Gyr. Note that the companion mass is very small because it has been depleted through the extended period of mass transfer that spun up its pulsar companion, as in "black widow" type binary MSPs (e.g., Rasio et al. 2000).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The MSP on the left of the figure has a different evolutionary path. The NS was not formed until about 6.6 Gyr in AIC as a member of a double WD binary, which itself was formed through a series of dynamical interactions, including both binary–single scatterings and exchanges. After the formation of the NS there were no stellar encounters and the WD companion kept transferring mass to spin up the NS to a MSP. The MSP at 12 Gyr has a magnetic field of 7.5 × 10⁷ G and a spin period of 1.57 ms, again with a low-mass companion of mass about 0.0075 M_⊙.

Figures 2 and 3 show the evolution of the magnetic fields and spin periods of the same two MSPs as in Figure 1. The black dashed line is the death line (see Equation (4)), below which radio pulsars disappear because they can no longer support pair production.

Zoom In Zoom Out Reset image size

In Figure 2, for the first 6 Gyr the magnetic field of the pulsar decreased due to magnetic field decay, and the spin period increased slowly due to magnetic dipole radiation. The pulsar exchanged into a binary at 6.31 Gyr with an MS star companion, and was accreting material from the MS star for about 200 Myr, causing it to spin up from angular momentum transfer. During this time the pulsar temporarily went below the death line. The MS star then evolved to a giant star and went through a common-envelope phase, during which the pulsar was not accreting mass. The binary orbit was circularized and shrunk by the common envelope. At about 6.7 Gyr, the MS star became a WD and started mass-transferring, continuing to spin up the pulsar. The pulsar emerged out of the "graveyard" and was spun up to an MSP at late times.

Figure 3 is similar to Figure 2 but for the other MSP on the left side of Figure 1. There were no dynamical interactions between the time the NS formed and 12 Gyr. In this case the NS–WD binary stayed intact after AIC and the WD continued through RLOF and spun the pulsar up to a MSP. A few studies have suggested that NSs formed in AIC are unlikely to be spun up to MSPs in the same binary (e.g., Verbunt & Hut 1987), because MSPs formed in this way cannot explain the overabundance of LMXBs and MSPs in GCs. However, this contrasts the results from our models. In the progenitor binaries of AIC NS, the companions may fill the Roche lobe and transfer mass onto the ONe WD, which likely leads to common-envelope evolution in BSE. The common-envelope evolution circularizes and shrinks the binary orbits, producing tight (semimajor axis about 0.001 au) WD–WD or NS–WD binaries. Because the binaries are so tight, the same WD companions can continue to fill the Roche lobe and spin up the NSs to MSPs.

Consider all NSs still in our model clusters at 12 Gyr. These represent about 8% of all NSs ever born in the N = 8 × 10⁵ models, and about 15% in the N = 3.5 × 10⁶ model. Of all these retained NSs at 12 Gyr, in the N = 8 × 10⁵ models (models 1–25), about 5% came from CCSN NSs, while about 95% came from ECSNe. For the large-N model (model 26), about 15% NSs came from CCSNe, and about 85% came from ECSNe. For models with fewer retained BHs, a larger fraction of NSs retained at 12 Gyr were formed through dynamically influenced channels (mainly AIC and MIC): about 0.5% of the NSs that are retained in BH-rich (models 1–10) and BH-intermediate models (models 11–19) and about 5% in BH-poor models (models 20–26) are formed through these channels. In all the models, most of the retained NSs at 12 Gyr were formed through EIC.

Now we turn to the sample of NSs in our models that would be potentially detectable as radio pulsars at 12 Gyr (all NSs above the death line defined by Equation (4)). In total, there are 208 pulsars in all our models, including 142 MSPs (i.e., pulsars with P ≲ 30 ms) and 66 young pulsars. Note that these numbers do not take into account beaming (and the beaming fraction is thought to be much smaller for young pulsars; see, e.g., Lorimer 2008, and references therein). Because CCSN NSs are such a small fraction of all NSs retained at late times, this class of NS does not contribute significantly to the pulsar population observable in old GCs. In contrast, 26% of all MSPs in our models were formed through EIC, 68% through AIC, and only 5% through CCSNe and 1% through MIC. Therefore, ECSN NSs are the principal source for MSPs in GCs. MIC NSs do not contribute much in general to MSPs. This is because there are only a few MIC NSs in the models (the total number is comparable to CCSN NSs in BH-poor models, but almost zero in the other models), and most of them formed at late times (they do not have enough time to be spun up).

Most of the MSPs in our models are in binaries, with a binary fraction of about 80%, and most of the companions are low-mass stars (only 9 companions are MS stars). There are also a few DNSs and NS–BH binaries in our models (see Section 3.4). Single MSPs result mostly from dynamical encounters in which the MSPs are exchanged out of the binaries.

We find a clear anti-correlation between the number of retained BHs and the number of MSPs in our models. Figure 4 shows the number of BHs and MSPs between 9 and 12 Gyr in models 1–25. We include multiple points from the same model sampled at different times, so a single model can provide different numbers of BHs and MSPs in the figure. For example, there are four models at N_BH = 1, showing a range of N_MSP the models can have for the same BH number. For cluster models with only a few retained BHs, there can be as many as 16 MSPs in the cluster, while for cluster models with more than 200 BHs, there are at most 2 MSPs. In the large-N model (model 26) there are 83 MSPs and 20 BHs at 12 Gyr, which we do not include in the figure, simply so that we can isolate the BH–MSP relation for models with similar N. It is worth noting that if the AIC MSPs are excluded from Figure 4, the anti-correlation between the numbers of BHs and MSPs still holds.

Zoom In Zoom Out Reset image size

This anti-correlation was anticipated, based on our understanding of BH populations in GCs. Figure 5 shows how the BHs dynamically influence the NSs, including the MSPs in the clusters. As long as many are present, BHs dominate the cluster cores because of mass segregation and they prevent the NSs from concentrating in the high-density central region. Furthermore, GCs with more retained BHs have lower core densities due to the heating of the cores from BH interactions (Arca Sedda et al. 2018; Fragione et al. 2018b; Kremer et al. 2018). The upper panel of Figure 5 shows the distribution of the 2D-projected radii of all NSs (step histograms) and only MSPs (filled histograms) in models 1–25. For models with fewer than 10 BHs, the NSs are located closer to the GC centers; while for models with a large number of BHs, the projected radii for most of the NSs are about an order of magnitude larger. This also affects the number of encounters the NSs can have during the cluster evolution. The stellar densities are higher toward the cluster centers, where the average stellar densities within the median 2D-projected radii (Figure 5) for the BH-rich, BH-intermediate and BH-poor models are about 6.6 × 10³ pc⁻³, about 8.9 × 10⁴ pc⁻³ and about 1.0 × 10⁶ pc⁻³, respectively. As a result, NSs located closer to the centers go through more dynamical interactions. Therefore, NSs in the BH-poor models are more likely to acquire companions and accrete mass, and thus have a higher chance of becoming MSPs at late times. The numbers of encounters of all NSs and MSPs are shown in the lower panel of Figure 5. This trend can also be seen in the upper panel, where NSs in the BH-poor models are scattered more and have a wider radial distribution, in contrast to NSs in the BH-rich models.

Zoom In Zoom Out Reset image size

The number of BHs can also be constrained by other measurable quantities such as the core radius of a cluster, which is strongly correlated with the number of retained BHs (Table 1). In our models, retaining more BHs leads to larger core radii, lower central (non-BH) stellar densities, and lower rates for all dynamical interactions. Thus, the basic properties of the MSP population (number and radial profile) and the core radius correlate similarly with the number of BHs present in a cluster. There might also be intrinsic differences between the properties of the MSPs in denser versus less dense clusters (such as systematic differences in spin periods or binary companion types). However, we do not have sufficient coverage of the statistics in either the observed sample (also affected strongly by selection bias) or in our current set of models (too few MSPs in low-density clusters; 11 total in models 1–10) to study this at present.

Metallicity has a minor effect on the number of retained NSs in clusters and is unlikely to affect this anti-correlation between the number of retained BHs and the number of MSPs in a cluster (de Menezes et al. 2019). However, we intend to perform a more detailed study of the effects of changing metallicity in these models in future work.

We plot all the pulsars in our models on top of all observed pulsars in the GC pulsar catalog⁹ in Figure 6. The intrinsic spin period derivatives of the model pulsars are derived from their magnetic fields using $\dot{P}=K\tfrac{{B}^{2}}{P}$ (see Section 2.2.1). In the figure, the blue dots show the spin periods and spin period derivatives of the observed GC pulsars and the orange dots show the intrinsic P and $\dot{P}$ values for model pulsars. We calculate the maximum accelerations the pulsars can have in their cluster potential using Equation (2.5) of Phinney 1992. The upper and lower limits for the "observable" $\dot{P}$ values of model pulsars with accelerations are shown as green circles. We also include a few pulsars in 47 Tuc with derived intrinsic $\dot{P}$ values, and in Terzan 5 with inferred magnetic fields, in Figure 6; these are shown by red stars and triangles, respectively (Freire et al. 2017; Prager et al. 2017).

Zoom In Zoom Out Reset image size

There are both MSPs ($P\lesssim 30\,$ ms) and young pulsars in our models (Figure 6), as are observed in GCs. Most of the pulsars observed in Galactic GCs are MSPs, as expected because MSPs have long lifetimes and can exist in GCs for many Gyr. In contrast, young pulsars have relatively short lifetimes, and those formed at early times in GCs are no longer there. However, through dynamical interactions such as collisions between an MS star and a WD, young pulsars can be formed at the present time in old GCs. Almost all young pulsars in our models were formed by collisions at late times (Figure 6), including newborn NSs formed in WD–MS star or WD–WD collisions, and old NSs that partially accreted during collisions with MS or giant stars. It has also been suggested that (apparently) young pulsars in GCs could be formed by partial recycling of NSs in mass-transferring binaries (Verbunt & Freire 2014; de Menezes et al. 2019). However, this does not happen at a significant rate in our models: only three young pulsars in our models were formed by partially spinning up a NS through accretion in a binary.

Figure 6 shows that the spin periods and spin period derivatives of our model pulsars are in reasonable agreement with observations, suggesting that our current very simple treatment for the evolution of magnetic fields and spin periods of NSs in CMC is sufficient for this first attempt at a detailed comparison.

DNSs and NS–BH binaries provide important probes for general relativity and their mergers can now be detected as gravitational-wave sources and can power short gamma-ray bursts visible throughout the universe (e.g., Belczynski et al. 2002a; Clausen et al. 2014). Such systems may form in a variety of astrophysical environments, including through isolated massive binary evolution in galactic fields (e.g., Belczynski et al. 2002b), as well as through dynamical processes in dense star clusters such as galactic nuclei (e.g., Fragione et al. 2018a) and GCs (e.g., Sigurdsson & Phinney 1995).

So far, only one DNS system has been identified in a GC, PSR 2127 + 11 C in M15.¹⁰ M15C has a spin period of 30.5 ms and a derived magnetic field of 1.2 × 10¹⁰ G. Both the pulsar and its NS companion have masses 1.35 M_⊙. This DNS system has an 8 hr orbital period and a highly eccentric orbit with e = 0.68 (Anderson et al. 1990; Prince et al. 1991; Deich & Kulkarni 1996; Jacoby et al. 2006). The projected radius of M15C from the center of M15 is 2.7 pc, which, combined with its eccentricity, suggests that it was most likely formed by an exchange interaction in the cluster core with recoil to its current location (Phinney & Sigurdsson 1991).

Primordial DNSs must be very rare in GCs because of natal kicks for NSs that break up the progenitor binaries or eject them out of the clusters. However, one way of producing these systems in GCs is through dynamical interactions. M15 is "core-collapsed", i.e., it has an extremely high stellar density near its center (Harris 1996, 2010 edition), which provides a good environment for DNS formation. NS–BH binaries should be even more rare than DNSs in GCs, even taking dynamics into account. Indeed, in clusters with many BHs, mass segregation and the strong heating by BH interactions prevent the NSs from interacting with the BHs (see Figure 5), while in clusters with few BHs, NSs have few interactions with BHs compared to other stars. It is therefore not surprising that no NS–BH binary has ever been detected. Based on simple models of the dynamical evolution of BH binaries in fixed GC backgrounds (described by multi-mass King models), Clausen et al. (2014) estimated that there are at most ∼10 BH–MSP binaries in the entire Milky Way GC system. We plan to further explore the formation and merger rates of NS–BH binaries in GCs in future studies using CMC.

In our current models we found that several NS–BH binaries and DNSs formed dynamically. The last two columns of Table 1 show the total numbers of these systems that appeared between 9 and 12 Gyr in our models. Two of the DNSs contain MSPs, while 32 DNSs and 6 NS–BH binaries contain young pulsars. Most of the NS–BH binaries and DNSs in our models appeared for a short time (<100 Myr) and were subsequently disrupted by dynamical encounters. However, there are also a few that survived for longer times. For example, in the large-N model there is a long-surviving DNS formed through a binary–binary interaction involving an MSP of 1.32 M_⊙ and an NS of 1.24 M_⊙ at about 11.5 Gyr. In the final model at 12 Gyr the binary has nearly completed its inspiral, reaching an orbital period of just 8 minutes (the remaining time to merger is about 0.03 Myr).

In total, we identify 40 unique NS–BH binaries and 93 unique DNSs in our models. Of these, 11 NS–BH binaries and 22 DNSs have gravitational-radiation inspiral times less than a Hubble time, and 9 NS–BH binaries and 15 DNSs merged during the model clusters' evolution time. If such binaries merge in the local universe, they may be observed as gravitational-wave sources by detectors such as LIGO/Virgo (Abbott et al. 2017) and they may produce short-hard gamma-ray bursts (Meszaros 2006; Lee & Ramirez-Ruiz 2007). A more detailed study of the merger rates and properties of NS–BH binaries and DNSs in GCs will be presented in a follow-up paper.