Millisecond Pulsars and Black Holes in Globular Clusters

TL;DR

This study uses advanced N-body simulations to explore how millisecond pulsars form and distribute in globular clusters, revealing their relationship with black hole populations and providing insights into their observable properties.

Contribution

The paper introduces realistic N-body simulations to model MSP formation in GCs, highlighting the impact of dynamical interactions and black hole populations.

Findings

MSPs can be formed via electron-capture supernovae and accretion in GCs.

Number of MSPs correlates with cluster mass and anti-correlates with black hole retention.

MSPs are more centrally concentrated in clusters with fewer black holes.

Abstract

Over a hundred millisecond radio pulsars (MSPs) have been observed in globular clusters (GCs), motivating theoretical studies of the formation and evolution of these sources through stellar evolution coupled to stellar dynamics. Here we study MSPs in GCs using realistic $N$-body simulations with our Cluster Monte Carlo code. We show that neutron stars (NSs) formed in electron-capture supernovae (including both accretion-induced and merger-induced collapse of white dwarfs) can be spun up through mass transfer to form MSPs. Both NS formation and spin-up through accretion are greatly enhanced through dynamical interaction processes. We find that our models for average GCs at the present day with masses $\approx 2 \times 10^5\,M_\odot$ can produce up to $10-20$ MSPs, while a very massive GC model with mass $\approx 10^6\,M_\odot$ can produce close to $100$. We show that the number of MSPs is anti-correlated with the total number of stellar-mass black holes (BHs) retained in the host cluster. The radial distributions are also affected: MSPs are more concentrated towards the center in a host cluster with a smaller number of retained BHs. As a result, the number of MSPs in a GC could be used to place constraints on its BH population. Some intrinsic properties of MSP systems in our models (such as the magnetic fields and spin periods) are in good overall agreement with observations, while others (such as the distribution of binary companion types) less so, and we discuss the possible reasons for such discrepancies. Interestingly, our models also demonstrate the possibility of dynamically forming NS--NS and NS--BH binaries in GCs, although the predicted numbers are very small.

Figures

Tables

Table 1. Table 1: Cluster Model Properties

$Model$	$N$	$\frac{{\sigma }_{BH}}{{\sigma }_{NS}}$	$r_c$	$r_{hl}$	$M_{TOT}$	$N_{BH}$	$N_{NS}$	$N_{NS-DYN}$	$N_{PSR}$	$N_{MSP}$	$N_{sMSP}$	$N_{bMSP}$	$NS-BH$	$DNS$
	${10}^5$		$\mathrm{pc}$		${10}^5{M}_{\odot }$
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

Equations

$$\dot{P}=K\frac{B^2}{P},$$

$$B=B_0\exp \left(-\frac{T}{\tau }\right),$$

$$B=\frac{B_0}{1+\frac{\Delta M}{10^{-6}{M}_{\odot }}}\exp (-\frac{T-t_{\mathrm{acc}}}{\tau }),$$

$$\frac{B}{P^2}=0.17\times 10^{12}\mathrm{G}{\mathrm{s}}^{-2}$$

Full text

globular clusters: general — stars: neutron — pulsars: general — stars: kinematics and dynamics — methods: numerical

††software: CMC (Joshi et al., 2000, 2001; Fregeau et al., 2003; Fregeau & Rasio, 2007; Chatterjee et al., 2010, 2013b; Umbreit et al., 2012; Morscher et al., 2015; Rodriguez et al., 2016), Fewbody (Fregeau et al., 2004), BSE (Hurley et al., 2002), SSE (Hurley et al., 2000)

1 Introduction

Globular clusters (GCs) are known to be highly efficient at producing millisecond pulsars (MSPs). Since the discovery of radio MSPs in GCs in the 1980s (Lyne et al., 1987), multiple pulsar surveys have found 150 pulsars in 28 GCs111GC pulsar catalog: http://www.naic.edu/~pfreire/GCpsr.html (for reviews, see Camilo & Rasio, 2005; Ransom, 2008), including 38 in Terzan 5 and 25 in 47 Tuc. Although GCs make up only about $0.05\%$ of stars in the Milky Way, collectively, GCs contain more than one third of the total number of known MSPs in our Galaxy (Manchester et al., 2005)222ATNF pulsar catalog: http://www.atnf.csiro.au/research/pulsar/psrcat/.

GCs also contain many low-mass X-ray binaries (LMXBs; Clark, 1975) with neutron star (NS) accretors. The low surface magnetic fields ($\sim 10^{7}-10^{9}\,$G) and short spin periods ($\lesssim 30\,$ms) of MSPs suggest that they are “recycled” pulsars (Alpar et al., 1982) with LMXBs as their likely progenitors. Indeed some “transitional” MSPs, providing the link between LMXBs and MSPs, have recently been detected; these are observed to switch back and forth between phases of accretion-powered X-ray emission and rotation-powered radio emission. At present, there are three confirmed transitional MSPs, including one in a GC, M28 (Papitto et al., 2013), and two in the Galactic field (Archibald et al., 2009; Bassa et al., 2014; Roy et al., 2015). There are also a few additional candidates in GCs and in the field (Bahramian et al., 2018, and references therein). While the physics is far from being understood in detail, it is plausible that mass transfer onto old, slowly spinning NSs can “bury” their magnetic fields while at the same time spinning them up (e.g. Bhattacharya & van den Heuvel, 1991; Rappaport et al., 1995; Kiel et al., 2008; Tauris et al., 2012, and references therein). In addition, MSPs have low spin-down rates and thus long lifetimes ($\sim 10^{10}\,\rm{yr}$, compared to $\sim 10^{7}\,$yr for young pulsars), which makes them easier to observe in old stellar systems like Milky Way GCs.

The large numbers of MSPs and NS LMXBs suggest that a typical Galactic GC on average contain at least a few hundred NSs (e.g., Ivanova et al., 2008). However, observations show that the majority of NSs in the Galactic field are born with velocities $\gtrsim 200\,\rm{km\,s^{-1}}$ due to natal kicks associated with asymmetries in core collapse supernovae (CCSNe); see, e.g., Hobbs et al. (2005). Because of these large natal kicks, the majority of CCSN NSs born in GCs (where escape velocities are generally $<50-100\,\rm{km\,s^{-1}}$, even at early times when the clusters may have been more massive than at present) are ejected from the cluster at birth. This is seemingly at odds with the large numbers of NSs inferred to be present in GCs, which is often referred to as the NS “retention problem” (e.g. Pfahl et al., 2002b). However, the discovery of high-mass X-ray binaries with long orbital periods ($P_{orb}>30\,\rm{days}$) and low eccentricities ($e\lesssim 0.2$) (HMXBs; Pfahl et al., 2002a) suggests that some NSs must be born with very small natal kicks. Additionally, NSs born in massive binaries may be easier to retain in GCs, but Pfahl et al. (2002b) showed that the retention fraction of NSs formed through CCSNe in massive binaries is still at most a few percent, not enough to explain the large populations of NSs in GCs.

Later studies suggested that electron-capture supernovae (ECSNe) can solve the retention problem by producing many NSs with small natal kicks (Podsiadlowski et al., 2004; Ivanova et al., 2008). When electron capture occurs onto $\rm{Mg}^{24}$ and $\rm{Ne}^{20}$ in a degenerate ONeMg core, it triggers the core to collapse to a NS (Miyaji et al., 1980; Nomoto, 1984, 1987). The explosion energy of ECSNe is much lower than that of CCSNe for Fe cores. As a result, NSs formed in ECSNe could receive an order of magnitude smaller natal kicks compared to NSs formed in CCSNe (Podsiadlowski et al., 2004). Therefore a large fraction of ECSN NSs can be retained in GCs after formation (Kuranov & Postnov, 2006; Ivanova et al., 2008), in contrast to CCSN NSs, which are mostly lost.

The high stellar densities in GC cores lead to frequent dynamical encounters and high formation rates of NSs, MSPs and LMXBS (Clark, 1975; Hut et al., 1992; Pooley et al., 2003; Hui et al., 2010; Bahramian et al., 2013). For example, dynamical interactions can enhance the WD–WD merger rate (e.g. Shara & Hurley, 2002), which can lead to a higher merger-induced collapse rate of WDs to NSs (Saio & Nomoto, 1985). NS binaries are also created at an increased rate since single NSs can acquire companions through exchange interactions (Sigurdsson & Phinney, 1993, 1995; Rasio et al., 2000; Ivanova et al., 2008). Subsequent stellar evolution of the companion and hardening through repeated close flybys can then trigger Roche-lobe overflow (RLOF) and mass transfer, possibly turning the NS into a MSP. More than half of the observed MSPs in GCs are known to be in binaries.

Several studies have shown that the number of MSPs in GCs is correlated with the stellar encounter rate, $\propto\rho_{0}^{2}r_{c}^{3}/\sigma_{0}$, where $\rho_{0}$ is the central luminosity density, $r_{c}$ is the core radius and $\sigma_{0}$ is the central velocity dispersion (Verbunt & Hut, 1987; Hui et al., 2010; Bahramian et al., 2013). In clusters with larger encounter rates, NSs undergo more dynamical interactions, resulting in more NSs in binaries and binaries with shorter orbital periods. For example, while double NSs (DNSs) are very rare, the only confirmed GC DNS is in M15, which is a core-collapsed cluster with extremely high central density. It has been suggested that the frequent stellar encounters in M15 led to the formation of this DNS (Anderson et al., 1990; Prince et al., 1991; Phinney & Sigurdsson, 1991; Deich & Kulkarni, 1996; Jacoby et al., 2006). Additionally, BHs could have a strong influence on the formation of NS binaries by altering the evolution of their host GCs. Once formed, BHs quickly mass-segregate to the cluster core through dynamical friction. Several recent studies have shown that the long-term retention fraction of BHs in GCs strongly affects the GC core densities and the ability of other compact objects to participate in the core dynamics (Mackey et al., 2008; Morscher et al., 2015; Chatterjee et al., 2017; Arca Sedda et al., 2018; Kremer et al., 2018a; Fragione et al., 2018b).

Very few previous works have attempted to study the formation and evolution of MSPs in the context of fully realistic $N$-body simulations of GCs. Ivanova et al. (2008) modeled the formation and evolution of NSs in GCs using the population synthesis code StarTrack (Belczynski et al., 2002b, 2008) to follow single and binary star evolution in a fixed cluster background, and the small $N$-body integrator Fewbody (Fregeau et al., 2004) to compute dynamical interactions in the cluster core. They included ECSNe for NS formation and showed that the low natal kicks associated with these NSs were crucial for matching the observed numbers of MSPs in GCs. Here, we build upon this previous study by performing full, self-consistent $N$-body simulations for the cluster dynamics, and also incorporating the effects of NS magnetic field and spin period evolution during mass-transfer.

In Section 2, we describe how we model MSPs in GCs. We describe the set of models used for this study and our methods and assumptions for tracking the magnetic fields and spin periods of NSs. We also provide some examples of the complicated dynamical evolution of NSs in our models. In Section 3 we provide an overview of results. We explore the expected anti-correlation between the numbers of MSPs and BHs in GCs. We also compare our results broadly with pulsar observations. We discuss selection effects and summarize our findings in Section 4.

2 METHODS

2.1 Models

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, b) which 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; Umbreit et al., 2012; Pattabiraman et al., 2013; Chatterjee et al., 2013a; 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\times 10^{5}$, binary fraction $f_{b}=5\%$, virial radius $r_{v}=1\,\rm{pc}$, King concentration parameter $W_{o}=5$, galactocentric distance $r_{g}=8\,\rm{kpc}$ and metallicity $\rm{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_{\odot}$ to sample the initial stellar masses. NS remnants from CCSNe at formation receive natal kicks drawn from a Maxwellian distribution with a standard deviation $\sigma_{NS}=265\,\rm{km\,s^{-1}}$ (Hobbs et al., 2005).

Following Kremer et al. (2018a), 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., 2018b). 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., 2018a). All models were evolved for $12\,\rm{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 $\frac{\sigma_{BH}}{\sigma_{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 necessity 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.

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\times 10^{6}$, binary fraction $10\%$ and BH natal kick $\frac{\sigma_{BH}}{\sigma_{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\times 10^{6}\,M_{\odot}$, 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.

\startlongtable

2.2 Simulating NS evolution

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}\,\rm{G}$) and spin periods (range $30-1000\,\rm{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

[TABLE]

where $K=9.87\times 10^{-48}\,\rm{yr/G^{2}}$, $B$ is the surface magnetic field and $P$ is the spin period333Note that there is a typo in the value of $K$ in Kiel et al. (2008); our adopted value is the correct one..

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

[TABLE]

on a timescale $\tau=3$ Gyr (Kiel et al., 2008). Here $B_{0}$, 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 $\sim 100\,\rm{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 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

[TABLE]

where $t_{acc}$ is the duration of the NS accretion phase (during RLOF) and $\Delta M$ is the mass accreted444We use the standard value $10^{-6}\,M_{\odot}$ for the threshold accreted mass, as adopted in COSMIC. This ensures that the first factor acts as a simple switch, lowering the NS magnetic field as soon as even a small amount of mass is accreted. Changing the value adopted here, as long as it is very small, produces negligible changes in the overall results.. 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 a 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 a 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^{8}-10^{8.8}\,\rm{G}$ and a new spin period is drawn from the range $3-20\,\rm{ms}$. In other words, we assume that a MSP involved in a merger remains a 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\times 10^{7}\,\rm{G}$ (Kiel et al., 2008). We do not set a lower limit for the spin periods. We use the standard expression

[TABLE]

(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 magnetic field and spin period above the death line, a MSP as any pulsar with $P\lesssim 30\,$ms, and a young pulsar as any pulsar that is not a MSP.

2.2.2 Electron-capture Supernovae

In our models, NSs formed in ECSNe are the dominant type of retained NSs in the cluster, and in all NS-LMXBs and MSPs. We assume that an ECSN happens whenever an ONeMg WD reaches a critical mass $M=1.38\,M_{\odot}$ so that electron capture is triggered on $\rm{Mg}^{24}$ and $\rm{Ne}^{20}$ and the WD undergoes collapse from the sudden lack of electron pressure support (Miyaji et al., 1980; Nomoto, 1984, 1987). We give small natal kicks to ECSN NSs, drawn from a Maxwellian distribution with a dispersion $\sigma_{ECSN}=20\,\rm{km\,s^{-1}}$ (Kiel et al., 2008).

In CMC we have three different evolutionary paths that lead to ECSNe (Ivanova et al., 2008). The first one is evolution-induced collapse (EIC) of a single star with initial mass in the range $6-8\,M_{\odot}$ (Nomoto, 1984, 1987; Kiel et al., 2008). For main-sequence (MS) and giant stars, if their initial masses are smaller than the carbon ignition mass and larger than the critical NS formation mass, they can form NSs in ECSNe. Helium stars with masses between $1.6\,M_{\odot}$ and $2.25\,M_{\odot}$ (Hurley et al., 2000) also go through ECSNe. The second path is accretion-induced collapse (AIC) (Nomoto & Kondo, 1991; Saio & Nomoto, 2004). If an ONeMg WD accretes ONe or CO material from its companion during RLOF, the WD is assumed to collapse to a NS when its mass is larger than the ECSN critical mass and smaller than the maximum NS mass set by BSE, above which it will become a BH. The last path is the merger-induced collapse (MIC) of the product of a merger or collision between two WDs (Saio & Nomoto, 1985); these can be either ONe WDs or CO WDs. All three paths generally produce NSs in binaries, which can often lead to subsequent RLOF and the production of LMXBs and MSPs.

It is important to note that the mass range for ECSNe progenitors is under debate (see, e.g., Poelarends et al., 2017, and references therein). However, changing the mass range for EIC in our models has only a small effect on the number of retained NSs and MSPs at late times. This is because, by assumption, the ECSN kicks are not mass dependent. In general, since the mass range for helium cores to collapse in ECSNe is narrow, small changes to the EIC mass range will not affect the overall number of EIC NSs and MSPs significantly. Thus we simply adopt the SSE prescription as described above.

2.3 Influence of Dynamics on MSP Formation

After their formation at early times, most of the NSs in our models evolve like young pulsars in isolation until the end ($12\,\rm{Gyr}$). Their magnetic fields and spin rates slowly decrease until they die as pulsars. Dynamical interactions are essential in 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 a MSP on the right of the figure was formed in EIC at $0.11\,\rm{Gyr}$. It had a magnetic field of $6\times 10^{12}\,\rm{G}$ and a spin period of $320\,\rm{ms}$ at birth. The pulsar experienced its first encounter at $5.3\,\rm{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 Fig. 2). A third binary-single exchange encounter gave the pulsar a MS star companion. The MS star spun up the pulsar via RLOF-driven mass transfer (the first red-dashed line in Fig. 2). At this time it was not yet a 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 Fig. 2), which circularized and shrank the binary orbit. The remnant of the MS star was a WD which continued to fill its Roche lobe and spun up the pulsar to a MSP (The second red-dashed line in Fig. 2). The final system has a MSP with a $5.9\times 10^{7}\,\rm{G}$ magnetic field and a $1.57\,\rm{ms}$ spin period, and a companion with mass $0.0075\,M_{\odot}$ at $12\,\rm{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).

The MSP on the left of the figure has a different evolutionary path. The NS was not formed until about $6.6\,\rm{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\,\rm{Gyr}$ has a magnetic field of $7.5\times 10^{7}\,\rm{G}$ and a spin period of $1.56\,\rm{ms}$, again with a low-mass companion of mass about $0.0075\,M_{\odot}$.

Figures 2 and 3 show the evolution of 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 as they can no longer support pair production.

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\,\rm{Gyr}$ with a 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\,\rm{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 a 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\,\rm{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 is in contrast to 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 (semi-major axis about $0.001\,\rm{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.

3 RESULTS

3.1 Overview of Statistics

Consider all NSs still in our model clusters at 12 Gyr. These represent about $8\%$ of all NSs ever born in the $N=8\times 10^{5}$ models, and about $15\%$ in the $N=3.5\times 10^{6}$ model. Of all these retained NSs at 12 Gyr, in the $N=8\times 10^{5}$ 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\lesssim 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 does 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 Sec. 3.4). Single MSPs result mostly from dynamical encounters in which the MSPs are exchanged out of the binaries.

3.2 BH–MSP Anti-correlation

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\,\rm{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 4 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\,\rm{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.

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 (Fragione et al., 2018b; Arca Sedda 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 towards the cluster centers, where the average stellar densities within the median 2D-projected radii (Fig. 5) for the BH-rich, BH-intermediate and BH-poor models are about $6.6\times 10^{3}\,\rm{pc^{-3}}$, about $8.9\times 10^{4}\,\rm{pc^{-3}}$ and about $1.0\times 10^{6}\,\rm{pc^{-3}}$, 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, a larger chance to become 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.

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. However, we intend to perform a more detailed study of the effects of changing metallicity in these models in future work.

3.3 Spin Periods and Spin Period Derivatives

We plot all the pulsars in our models on top of all observed pulsars in the GC pulsar catalog555http://www.naic.edu/~pfreire/GCpsr.html in Figure 6. The intrinsic spin period derivatives of the model pulsars are derived from their magnetic fields using $\dot{P}=K\frac{B^{2}}{P}$ (see Sec. 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 Eq. (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).

There are both MSPs ($P\lesssim 30\,$ms) and young pulsars in our models (Fig. 6), as are observed in GCs. Most of the pulsars observed in Galactic GCs are MSPs, as expected since 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 a 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 (Fig. 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). However, this does not happen at a significant rate in our models: only 3 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.

3.4 DNSs and NS–BH Binaries

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 M15666There is also a candidate DNS in NGC 6544; see Lynch et al. (2012).. M15C has a spin period of $30.5\,\rm{ms}$ and a derived magnetic field of $1.2\times 10^{10}\,\rm{G}$. Both the pulsar and its NS companion have masses $1.35\,M_{\odot}$. This DNS system has an 8-hour 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\,\rm{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 Fig. 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) estimate that there are at most $\sim 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\,\rm{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\,\rm{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 a MSP of $1.32\,M_{\odot}$ and a NS of $1.24\,M_{\odot}$ 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\,\rm{min}$ (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.

4 Discussion and Summary

It is well understood that pulsar observations are limited by strong selection effects. The large distances of most GCs from Earth, Doppler shifts in short-period binaries, and dispersion of signals by the interstellar medium all make it difficult to detect radio pulsars. The luminosity function of GC pulsars, $d\log\,N\approx-d\log\,L$ (Hessels et al., 2007; Hui et al., 2010), suggests that some low-luminosity pulsars may remain unobserved.

Given these observational biases, the true number of pulsars in any given cluster remains quite uncertain. Many studies have attempted to constrain the true population of radio pulsars in Milky Way GCs (Kulkarni et al., 1990; Hui et al., 2010; Bagchi et al., 2011; Chennamangalam et al., 2013; Turk & Lorimer, 2013). For example, the predicted numbers of potentially observable pulsars in Terzan 5 and 47 Tuc are about $150$ and $80$, respectively, with large error bars ($\sim 50\%$) (Bagchi et al., 2011; Chennamangalam et al., 2013). Assuming the beaming fraction to be about $50\%$, there are around $70$ potentially observable pulsars in our model 26 (Table 1), a very promising result since this model has a mass at 12 Gyr close to that of 47 Tuc (e.g., Giersz & Heggie, 2011).

There are $5$ conspicuously sub-millisecond pulsars (sub-MSPs) in our models, as shown in Figure 6. The progenitor NSs were exchanged into binaries with eccentric orbits and with MS star or giant star companions in dynamical interactions, and they were spun up all the way to these very short spin periods through mass transfer during the evolution of their companion stars. However, there is no sub-MSP observed in any GC or in the Galactic field. Since we do not assume a lower limit for the NS spin period, sub-MSPs can appear in our models through continued mass transfer. The absence of detected sub-MSPs indicates that some physical mechanism not included in our models, such as gravitational radiation (see e.g. Bildsten, 1998; Chakrabarty et al., 2003), may set a fundamental limit on how fast a NS can be spinning. Additionally, many of our model MSPs are still accreting (at very low mass-transfer rates; $\dot{M}\lesssim 10^{-11}\,M_{\odot}\rm{yr^{-1}}$) from extremely low-mass ($\lesssim 0.01\,M_{\odot}$) companions at late times. These binaries could potentially be observed as either radio pulsars or LMXBs. The companions may eventually be completely evaporated (as observed in “black widow” systems), a process that we do not currently include in our treatment of binary evolution. In the absence of a sufficiently detailed treatment of these systems in BSE, we have assumed all of them to be binary MSPs in our current models.

About $20\%$ of all MSPs in our models are single. Of the remaining $\sim 80\%$ found in binaries, about $75\%$ of MSPs have WD companions and about $5\%$ have MS star companions. Observationally, there are 17 black-widow systems and 12 red-back systems in the GC pulsar catalog. Assuming that the ratio between the numbers of black widows and red backs is not affected by selection effects, the fraction of MSP binaries with H-rich companions (red backs) in our models seems too low, and the fraction of MSP binaries with low-mass WD companions (black widows) seems too high, compared to observations. As this is our very first attempt at studying the formation of MSPs within our CMC code, for simplicity, we did not explore the various uncertainties associated with the treatment of binary evolution in BSE, and instead we focused on general trends (for example the anti-correlation between MSPs and BHs). We do get good overall agreement with observations, in the sense that our models produce reasonable numbers of all observed type of systems (single vs binary MSPs, very low-mass companions vs MS star companions, slow vs fast pulsars). More closely matching observations would require a more sophisticated treatment of the binary evolution physics, which is beyond the scope of this analysis, but will be considered in future works.

To summarize, in this study, we have simulated pulsars in GCs using our CMC code. We updated BSE to incorporate ECSNe for NS formation via three channels: EIC for single stars, and AIC and MIC for binary stars. For ECSN NSs, we apply low natal kicks with $\sigma_{ECSN}=20\,\rm{km\,s^{-1}}$. They are the major source for MSPs in our GC models. We also incorporate magnetic field and spin period evolution of NSs in BSE to model realistic pulsars. In our models, magnetic field evolution follows simple prescriptions for magnetic field decay or lowering the field through mass accretion. Spin periods then evolve according to the change in magnetic field and angular momentum of the NSs.

Dynamical interactions are clearly the key to understanding the formation of MSPs in GCs. Stellar dynamics can enhance the formation of NSs through AIC and MIC, which are likely to be retained in GCs, and enhance mass transfer in binaries with NS accretors. Furthermore, binaries can be harden or get to higher eccentricities through interactions, which leads to mass accretion and spinning up of NSs.

We have studied the pulsar population in 26 GC models. We found that the number of MSPs is anti-correlated with the number of BHs retained in GCs. This results from the dynamical coupling of BHs and NSs in our models, where most of the pulsars are formed dynamically, including the few young pulsars that are active at late times. Additionally, NS–BH and NS–NS binaries are also more readily formed in models with few BHs retained. The large specific abundances of MSPs in GCs relative to the Galactic field can be explained naturally by dynamical formation. This also provides a way to estimate the number of BHs in GCs given the number of observed pulsars. Furthermore, the spin periods and spin period derivatives of our model pulsars agree reasonably with observations, showing that CMC is able to model realistic pulsar populations in clusters. Our code can also produce a realistically large number of pulsars in a large-$N$ model, opening the door to future detailed modeling of very massive GCs such as 47 Tuc (Giersz & Heggie, 2011).

We thank Craig Heinke, Paulo Freire, Ron Taam, and the anonymous referee for useful discussions and comments on the manuscript. Our work was supported by NASA ATP Grant NNX14AP92G and NSF Grant AST-1716762 at Northwestern University. Computations were supported in part through the resources and staff contributions provided for the Quest high-performance computing facility at Northwestern. C.S.Y. acknowledges support from the NSF GK-12 Fellowship Program under Grant DGE-0948017. K.K. acknowledges support from the NSF Graduate Research Fellowship Program under Grant DGE-1324585. S.C. acknowledges support from CIERA, and from NASA through Chandra Award TM5-16004X/NAS8-03060, issued by the Chandra X-ray Observatory Center (operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-03060). C.L.R. is supported by a Pappalardo Postdoctoral Fellowship at MIT. F.A.R. acknowledges support from NSF Grant PHY-1607611 while at the Aspen Center for Physics.