Single Millisecond Pulsars from Dynamical Interaction Processes in Dense Star Clusters

We use the CMC code (Rodriguez et al. 2022 and references therein) to model GCs and all objects therein. CMC is a Hénon-style Monte Carlo code (Hénon 1971a, 1971b) that incorporates various relevant physics for cluster evolution, including two-body relaxation, tidal mass loss, strong few-body interactions, dynamical binary formation through collisions with giant stars and tidal capture interactions (Ye et al. 2022), and post-Newtonian effects (Rodriguez et al. 2018b). Stellar and binary star evolution are fully coupled to the dynamical interactions and are computed by the publicly available software COSMIC (Breivik et al. 2020), which is based on SSE (Hurley et al. 2000) and BSE (Hurley et al. 2002). Strong three- and four-body gravitational encounters are directly integrated using the Fewbody package (Fregeau et al. 2004; Fregeau & Rasio 2007), which includes post-Newtonian dynamics for black holes (BHs; Antognini et al. 2014; Amaro-Seoane & Chen 2016; Rodriguez et al. 2018a, 2018b).

CMC forms and evolves NSs and MSPs following Ye et al. (2019) and references therein. NSs are formed either in iron core-collapse supernovae (CCSNe) or electron-capture supernovae (ECSNe) depending on the mass and metallicity of the progenitor star. NSs formed in CCSNe receive large natal kicks drawn from a Maxwellian distribution with a standard deviation ${\sigma }_{\mathrm{CCSN}}\,=265\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (Hobbs et al. 2005). ECSNe include scenarios where the core of a ∼6–8 M_⊙ star collapses to an NS (Nomoto 1984, 1987), or the collapse of a massive WD to an NS triggered by the capture of electrons and the loss of electron pressure after mass accretion from or a merger with another WD (Saio & Nomoto 1985; Nomoto & Kondo 1991; Saio & Nomoto1998, 2004; Shen et al. 2012; Schwab et al. 2015; Schwab 2021). The critical mass for the collapse of an ONe core or WD is assumed to be 1.38 M_⊙. For more details about the different outcomes of WD–WD coalescences in CMC, see Kremer et al. (2021b). NSs formed in ECSNe receive small natal kicks drawn from a Maxwellian distribution with a standard deviation ${\sigma }_{\mathrm{ECSN}}=20\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (Kiel et al. 2008; Kiel & Hurley 2009). The small kicks allow these NSs to be preferentially retained in the clusters given the low cluster escape velocities. We set the maximum mass of an NS to be 2.5 M_⊙. Any object more massive than this will collapse into a BH.

All new NSs born in CCSNe and ECSNe are assigned initial magnetic fields in the range 10^11.5–10^13.8 G and spin periods in the range 30–1000 ms, similar to observed young radio pulsars (ATNF Catalog). Isolated young pulsars spin down through magnetic dipole radiation, and their surface magnetic fields decay exponentially over a 3 Gyr timescale. ¹⁰ On the other hand, during mass accretion in a binary, a pulsar's magnetic field decreases inversely proportional to the mass accreted as ${(1+{\rm{\Delta }}M/{10}^{-6}\,{M}_{\odot })}^{-1}$ according to the "magnetic field burying" scenario (e.g., Bhattacharya & van den Heuvel 1991, and references therein), and the pulsar is spun up from angular momentum transfer (Hurley et al. 2002; Kiel et al. 2008).

In addition to these standard assumptions, we have implemented in this study several updates into CMC for forming single MSPs as will be discussed in detail in Section 3.

For comparisons to pulsar and cluster observations, we focus in particular on one cluster, NGC 6752. NGC 6752 is a typical core-collapsed cluster in the Milky Way (Harris 1996; 2010 edition), thus our comparison to this particular cluster is a reasonable proxy for core-collapsed clusters in general. It hosts the second largest number of observed MSPs in core-collapsed clusters, ¹¹ with nine observed MSPs, where eight of them are single, and one is in a binary (GC Pulsar Catalog). To simulate this cluster, we adopt the initial conditions of the cluster model that has shown good agreement with the observations of NGC 6752 from Kremer et al. (2020b). Our simulations all have an initial number of stars N = 800,000, initial virial radius R_v = 0.5 pc, metallicity Z = 0.0002, and Galactocentric distance R_g = 8 kpc. The initial mass function follows a standard Kroupa broken power law (Kroupa 2001), and the stars are sampled between 0.08 and 150 M_⊙. The initial binary fraction is 5% across all primary masses, and the secondary masses are drawn from a uniform distribution in the range 0.1–1 of the primary mass (Duquennoy & Mayor 1991). We use a King profile (King 1966) with a concentration parameter W₀ = 5 for the initial mass density distribution.

In total, we have run 18 NGC 6752–like simulations with the same initial conditions while allowing the new implementations specific for forming single MSPs (Section 3) to be varied. All simulations are run for 13.8 Gyr. We compare all model surface brightness profiles and velocity dispersion profiles at about 12 Gyr to the observations in Figure 2 (the observed cluster age is about 11–14 Gyr; Buonanno et al. 1986; Gratton et al. 1997, 2003; Correnti et al. 2016; Souza et al. 2020; Bedin et al. 2023). We adopt 4.125 kpc as the distance to the cluster (Baumgardt & Vasiliev 2021). Note that Vasiliev & Baumgardt (2021) estimated that the semimajor axis of NGC 6752's orbit in the Milky Way is about 4.5 kpc with a small eccentricity of about 0.25, different from the Galactocentric distance we adopt here. This difference mostly affects the tidal boundary of the cluster and is not likely to affect the dynamical interactions in the cluster core which determine the number of single MSPs formed. As shown, all 18 models closely match the observed profiles without fine-tuning of the cluster's initial conditions.

Zoom In Zoom Out Reset image size

Perhaps the most straightforward way to form single MSPs is through binary-mediated dynamical encounters where the NS binaries that passed previously through a low-mass X-ray binary mass-transfer phase are disrupted dynamically. In addition, MSP binaries may also be disrupted through the evaporation of companion stars. Very low-mass (≲0.01 M_⊙) CO/ONe WD companions of some MSP binaries may be unstable and may evaporate, leaving behind single MSPs (Bildsten 2002). At the present day, ∼10 single MSPs which formed from these two channels are found per time step (models 1a and 1b in Table 2 by adding N_lmco to ${N}_{\mathrm{MSP}}^{s}$). This number alone roughly matches the number of observed single MSPs in NGC 6752, but the number of potentially observable MSPs is almost certainly much larger than the number we know now, given the sensitivity of radio telescopes and MSP selection effects (also see Section 5). In addition, there are similar amounts of binary MSPs in models 1a and 1b (Table 2 by subtracting N_lmco from ${N}_{\mathrm{MSP}}^{b}$), and the extreme ratio of single-to-binary MSPs in NGC 6752 and some other core-collapsed GCs (Figure 1 and Table 1) is not explained. Here we further discuss the dynamical disruption of MSP binaries.

The rate of dynamical disruption of NS binaries is determined by the density of the host cluster (often defined using the so-called Γ parameter; e.g., Verbunt & Freire 2014) and the sizes of the NS binaries themselves, which set their cross sections for encounters. A key parameter that determines the cross sections of these binaries is the stability criterion of binary mass transfer. The stability of mass transfer during Roche-lobe overflow in our simulations is determined by the critical mass ratio, q_crit = M_donor/M_accretor. Smaller q_crit allows for unstable mass transfer from lower-mass donors and a wider range of donor-to-accretor mass ratios, which essentially corresponds to higher instability during mass transfer. In turn, this leads to higher rates of the common-envelope phase (note that we assume there is no mass accretion during the common-envelope process) and the formation of more compact NS–WD binaries. Indeed, in our previous study (Ye et al. 2019), many MSPs were spun up at least partially by a WD binary companion. These compact NS–WD binaries have very small encounter cross sections and are difficult to disrupt dynamically. On the other hand, higher values for q_crit lead in general to more stable mass transfer, which allows NSs to be spun up during accretion from a giant-branch companion, leading to a population of relatively wider MSP–WD binaries with larger cross sections for subsequent dynamical encounters that may break apart the binaries.

To demonstrate the effects of q_crit on the orbital separations of MSP binaries, we have run three NS isolated binary populations using COSMIC (Breivik et al. 2020). The initial mass function, metallicity, and the initial binary properties are sampled the same way as in the cluster simulations (Section 2). All three binary simulations have the same initial conditions, including the number of binaries, and we vary only the prescriptions for q_crit. We select the population of interest at the time of their formation. Figure 3 compares the orbital period distributions of MSP binaries for three different q_crit prescriptions. The black, coral, and blue histograms contain 3936, 261, and 2992 MSP binaries at a Hubble time, respectively. In general, the critical mass ratios from Belczynski et al. (2008; blue distribution) are larger for stars at the giant branch or the asymptotic giant branch (when the core mass of the giant is ≲90% of its total mass), which lead to a much smaller fraction of tight MSP binaries. In the dense, dynamical environments of GCs, larger fractions of wide binaries may lead to more single MSPs since wide binaries have larger cross sections and more frequent dynamical encounters; thus, the original MSP binaries may be more easily disrupted and leave behind single MSPs.

Zoom In Zoom Out Reset image size

To test the role of this mechanism in clusters, we run two cluster simulations using the critical mass ratio prescriptions from Belczynski et al. (2008; models 9a and 9b in Table 2). All other simulations use the prescriptions from Hurley et al. (2002) and Hjellming & Webbink (1987), as was adopted in our previous paper (Ye et al. 2019). Model 9a has about four times more single MSPs than model 1a (not including single MSPs from WD–WD coalescences, although model 9b only has about the same number of single MSPs as model 6b), but is still not able to reproduce the observed ratio of single-to-binary MSPs. This suggests that the critical mass ratio for unstable mass transfer may play a role in forming single MSPs in GCs, even though the evolution of binaries in core-collapsed clusters is likely strongly influenced by dynamical encounters.

Second, NSs could be spun up through disruptions of main-sequence stars and accretion of the debris, producing single MSPs directly (Kremer et al. 2022). This phenomenon would most frequently occur in core-collapsed GCs. CMC self-consistently tracks close encounters between NSs and main-sequence stars and computes the tidal disruption process.

The accretion of debris following the tidal disruption of a main-sequence star may spin up an NS to millisecond periods (e.g., Davies et al. 1992; Lee et al. 1996; Kremer et al. 2022). NS–main-sequence star TDEs can occur during close encounters with pericenter distance

$$\begin{eqnarray}&&{r}_{p}\lt {\left(\displaystyle \frac{{M}_{\mathrm{NS}}}{{M}_{* }}\right)}^{1/3}{R}_{* },\end{eqnarray} \tag{ 1 }$$

where M_NS is the initial mass of the NS, and M_* and R_* are the initial mass and radius of the main-sequence star, respectively. For this study, we adopt the prescriptions in Kremer et al. (2022) for spinning up NSs through TDEs.

The mass accretion rate is calculated as

$$\begin{eqnarray}&&\begin{array}{l}{\dot{M}}_{\mathrm{acc}}\approx \left(\displaystyle \frac{{M}_{d}}{{t}_{v}}\right){\left(\displaystyle \frac{{R}_{\mathrm{acc}}}{{R}_{d}}\right)}^{s}\\ \times {\left[1+3(1-C)\left(\displaystyle \frac{t}{{t}_{v}}\right)\right]}^{-\tfrac{1+3(1+2s/3)(1-C)}{3(1-C)}},\end{array}\end{eqnarray} \tag{ 2 }$$

where M_d , R_d , and t_v are the initial disk mass, disk radius, and viscous accretion time, respectively. We set M_d = 0.9 M_* following the hydrodynamic simulations of Kremer et al. (2022), where it is shown that, in general, ∼90% of the initial stellar mass is bound to the NS following these TDEs. This is a reasonable assumption for encounters in GCs where the relative velocity of a pair of objects is in general much less than the stellar escape velocity. We assume that R_d = 2 r_p , t_v = 1 day, and R_acc equals the radius of the NS (Kremer et al. 2022). The exponent s ∈ [0, 1] is a free parameter (e.g., Blandford & Begelman 1999) defining the amount of material transported from the outer edge of the disk to the NS. This parameterization accounts for the fact that in these highly super-Eddington disks, an uncertain fraction of disk material is expected to be launched from the disk via winds (e.g., Metzger et al. 2008; Kremer et al. 2022). We adopt C = 2s/(2s + 1), equivalent to the assumption that the disk wind outflow produces no net torque on the disk (e.g., Kumar et al. 2008).

The accretion rate of the spin angular momentum is estimated as

$$\begin{eqnarray}&&{\dot{J}}_{\mathrm{acc}}\approx {\dot{M}}_{\mathrm{acc}}\left(\sqrt{{{GM}}_{\mathrm{NS}}{R}_{\mathrm{acc}}}+\displaystyle \frac{1}{2}{M}_{\mathrm{acc}}\sqrt{\displaystyle \frac{{{GR}}_{\mathrm{acc}}}{{M}_{\mathrm{NS}}}}\right),\end{eqnarray} \tag{ 3 }$$

where M_acc is the total mass accreted onto the NS.

We analytically integrate Equations (2)–(3) (Equations (3)–(4) in Kremer et al. 2022) for 10⁶ s to compute the total mass and angular momentum accreted onto the NS. During the TDE, the NS is spun up accordingly through angular momentum transfer (Hurley et al. 2002, Equation (54)), and its magnetic field is assumed to decay following the same "magnetic field burying" scenario as in Ye et al. (2019, their Equation (3)) for mass accretion in a binary.

Figure 4 shows the properties of all NSs in the simulations that tidally disrupted main-sequence stars. ¹² The amounts of mass accreted onto the NSs depend on the mass inflow rate and can differ by ∼three orders of magnitude between high (s = 0.2) and low (s = 0.8) efficiency (upper right panel; Figure 3 in Kremer et al. 2022). In the optimistic case where the inflow rate is high, almost all NSs that went through main-sequence star TDEs were spun up to MSPs (lower panels). This can boost the number of single MSPs in a typical core-collapsed cluster by a factor of ≳10 (Table 2). On the other hand, low inflow rates increase the numbers by a factor of ≲3. In addition, 11 NSs in all relevant simulations exceeded the maximum NS mass after accreting from the debris of the disrupted main-sequence stars and collapsed to become BHs.

Zoom In Zoom Out Reset image size

We want to point out that ∼10% of the NS–main-sequence star TDEs include more than two objects in the close encounters. The outcomes of these few-body interactions are uncertain and require careful studies with hydrodynamics simulations. Here we allow these NSs to be spun up to MSPs to maximize the single MSPs formed in a model.

The average rate of NS–main-sequence star TDEs is about 60 Gpc⁻³ yr⁻¹ at ≳9 Gyr (up to redshift z ∼ 0.5) for typical core-collapsed clusters assuming a GC number density of 2.31 Mpc⁻³. Figure 5 illustrates the time of collisions or TDEs between an NS and a main-sequence star in all simulations. If a fraction of the energy is converted to luminosity during accretion, we can roughly estimate the luminosity to be ${L}_{\mathrm{TDE}}\sim \eta \dot{{M}_{\mathrm{acc}}}{c}^{2}$, which is ∼2 × 10⁴³ erg s⁻¹ assuming η = 0.01 and ${\dot{M}}_{\mathrm{acc}}\sim 0.04\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ (Kremer et al. 2022).

Zoom In Zoom Out Reset image size

In reality, the outcomes of these TDEs are uncertain and additional effects may limit the total mass growth by the NS. For instance, accretion feedback (e.g., Armitage & Livio 2000; Papish et al. 2015) or feedback due to nuclear burning on the NS surface (e.g., Hansen & van Horn 1975; Taam 1985; Bildsten 1998) may unbind material initially bound to the NS and inhibit growth, provided the feedback energy can couple mechanically with the bound material. Our intent here is to explore whether these TDEs may account for the populations of observed single MSPs in the most optimal scenario where these feedback effects play a negligible role.

Finally, we also include WD–main-sequence star TDEs in CMC. Following the treatments in Kremer et al. (2019), we assume an immediate collision between a WD and a main-sequence star if their pericenter distance during the first passage is smaller than the tidal disruption radius of the main-sequence star, ${({M}_{\mathrm{WD}}/{M}_{* })}^{1/3}{R}_{* }$. The mass of the collision product is the total mass of the WD and the main-sequence star, so the product may evolve to be a more massive WD (Hurley et al. 2002), increasing the WD interaction rates and the probability of subsequent double WD collisions and NS formation.

A third way to form single MSPs is through the merger-induced collapse of WDs. To understand this channel fully, we have implemented tidal capture interactions between two WDs into CMC. During single–single or few-body interactions, if the pericenter distance of the approaching WDs is smaller than twice the sum of their radii (see the Appendix for more details), we assume an immediate collision of the two WDs. This is a reasonable assumption since tidally captured double WD binaries are very compact and will merge in a short timescale (e.g., two solar-mass WDs captured at the maximum pericenter distance and with eccentricity e = 0.99 will merge in about 10,000 yr, less than the typical timescale between dynamical encounters in a typical GC). ¹³

In this study, we focus on WD coalescences that can lead to the formation of NSs. These include both direct coalescences when two WDs physically collide (i.e., the pericenter distance of the encounter is smaller than the sum of the radii of the two WDs) and coalescences that follow more distant encounters leading to tidal capture and gravitational-wave inspiral. In addition to these dynamically mediated scenarios, WD coalescences can also occur following unstable mass transfer in a gravitational-wave inspiralling WD–WD binary (e.g., Webbink 1984; Nelemans et al. 2001; Marsh et al. 2004; Kremer et al. 2021b). We assume unstable mass transfer ensues for WD–WD binaries with mass ratio q > 0.628 (e.g., Breivik et al. 2020). Such mass-transfer-mediated binary mergers are qualitatively similar to the physical collisions and tidal capture encounters facilitated by close flybys, since in both cases the coalescence occurs on the dynamical timescale of the WDs involved. Thus, we assume that all of these scenarios result in the same outcome. ¹⁴

Traditionally, coalescences of massive (i.e., super-Chandrasekhar) WD pairs have been connected to Type Ia SNe (e.g., Webbink 1984). However, other studies have noted that in some cases, mechanisms such as runaway electron capture (e.g., Nomoto & Iben 1985) or off-center carbon ignition leading to the formation of a low-mass iron core (e.g., Schwab et al. 2016) may allow a thermonuclear explosion to be avoided. In the latter cases, collapse to an NS is the likely result. Here we assume all super-Chandrasekhar coalescences where the components are CO WDs and/or ONe WDs lead to collapse (e.g., also Kremer et al. 2023a), thus maximizing the formation rate of single NSs through this channel. The outcomes of various combinations of super-Chandrasekhar mass WD coalescences adopted in our simulations are shown in Table 3. We also list the assumed outcomes for cases where a WD accretes to the Chandrasekhar limit via stable Roche-lobe overflow in a binary (assumed for mass ratios q < 0.628), which is a distinct process from dynamical-timescale coalescences.

Following the dynamical coalescence of a massive WD pair (e.g., Dan et al. 2014), subsequent phases of viscous and thermal evolution (e.g., Schwab et al. 2012; Shen et al. 2012) that can last as long as 10⁵ yr and potentially include phases of dusty wind mass loss (e.g., Schwab et al. 2016) determine whether an NS forms and, if so, its final properties (e.g., mass, spin period, and magnetic field; Schwab 2021). Simulating in detail these various phases is beyond the scope of CMC, so in order to bracket some of the uncertainties inherent to these phases, we adopt two limiting assumptions in our simulations. As an upper limit (model set b in Table 2), we assume that the total mass lost during these phases is negligible so that the total mass of the WDs is conserved up through collapse. Accounting for the fractional loss of mass through gravitational-wave emission, we then assume the final masses of the NSs born following the collapses are 90% of the initial WD pair. In practice, this choice serves to maximize the rate of single NS formation.

As an alternative lower-limit scenario, the other half of the simulations (model set a) instead follow the default prescriptions from Breivik et al. (2020). In these simulations, the more massive WD accretes a total mass ${\rm{\Delta }}M={\dot{M}}_{\mathrm{edd}}{\tau }_{\dot{{\rm{M}}}}$ from the donor WD. Here ${\dot{M}}_{\mathrm{edd}}$ is the Eddington accretion limit and ${\tau }_{\dot{{\rm{M}}}}=\sqrt{{\tau }_{\mathrm{KH}}{\tau }_{\mathrm{dyn}}}$ is the characteristic mass-transfer timescale where τ_KH and τ_dyn are the Kelvin–Helmholtz and dynamical timescale of the donor WD, respectively. For WDs that are ∼1 M_⊙, ${\dot{M}}_{\mathrm{edd}}\approx {10}^{-5}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ and ${\tau }_{\dot{{\rm{M}}}}\approx 2000\,$ yr, and thus a very small amount of mass ΔM ≈ 10⁻³ M_⊙ is accreted (Breivik et al. 2020; Kremer et al. 2021b). The masses of all NSs formed in merger-induced collapses in these simulations (set a) are fixed to 0.9 × 1.38 M_⊙, the same as for NSs formed in accretion-induced collapses.

The number of WD–WD coalescences that collapse to NSs in each simulation is shown in Table 2. ¹⁵ Model set b differs from set a by (1) the mass of the NSs formed in merger-induced collapse, with model set b having more massive NSs; and (2) the amount of mass accreted onto the primary WD during dynamical mass transfer, with model set b conserving the mass and having a larger number of NSs formed from dynamical mass transfer between two WDs. Note that even without additional dynamical interactions such as WD–WD tidal captures and NS spun up through TDEs (Section 3.2), the number of single MSPs is boosted simply by assuming that the dynamical mass transfer in double WD binaries and the subsequent collapses of the super-Chandrasekhar-mass merger products conserve mass (e.g., model 1b versus 1a, where there are more mergers in 1b).

Overall, model set b (with the updated mass prescription) has more combined WD–WD coalescences that produce NSs than model set a (with the default mass prescription). This is largely because the conservation of mass during dynamical mass transfer in double WD binaries leads to many more NS remnants from mergers, despite model set b having fewer NS remnants from collisions. The decreased amount of collisions is probably caused by mergers reducing some of the WDs available for collisions and the larger number of low-mass BHs leading to larger WD radial distributions through energy release during BH dynamical encounters (e.g., Kremer et al. 2020a). Figure 6 illustrates the radial distributions of the compact objects within 1 pc from the cluster centers. The WDs in model set a (using the default mass prescription) are slightly more concentrated than those in model set b (using the updated mass prescription).

Zoom In Zoom Out Reset image size

The masses of the component WDs involved in the coalescences are plotted in Figure 7. Overall, CO WD–CO WD interactions account for ∼70% of mergers and ∼60% of the collisions. The rest of the collisions or mergers include at least one ONe WD. The overdensity of WDs at around 1 M_⊙ comes from the zero-age main-sequence (ZAMS) mass to WD mass relation we adopt for low metallicities (Hurley et al. 2000; Figure 1 in Kremer et al. 2021b). These WDs are the remnants of ZAMS stars of about 2.5–3.5 M_⊙.

Zoom In Zoom Out Reset image size

We also show the distributions of component mass, mass ratio, and coalescence time of the WDs in Figure 8. In addition to the peaks at about 1 M_⊙ as mentioned above, the peak at roughly 0.6 M_⊙ arises from the initial mass function of ZAMS stars of about 1 M_⊙ (also see Kremer et al. 2021b). Meanwhile, the mass segregation of massive WDs and previous collisions/mergers contribute to the peak at around 1.3 M_⊙.

Zoom In Zoom Out Reset image size

Most mergers and collisions have mass ratios larger than 0.6. The peak at around 1 for collisions can be naturally explained by mass segregation and more frequent close encounters of massive WDs in the cluster cores (see also Figure 6). Furthermore, the critical mass ratio for unstable mass transfer in double WD binaries (0.628 in all simulations) limits the mergers to a narrower mass ratio range than the collisions.

Lastly, almost all WD–WD collisions and about 60% of WD–WD mergers occur at late times (≳6 Gyr) of the host clusters' evolution. At early times, BHs dominated the cluster cores because of mass segregation, while many WDs have yet to form. Frequent dynamical interactions eject many BHs from the clusters in a few billion years (e.g., Kremer et al. 2020a) and the most abundant WDs come to dominate the cluster cores and engage in dynamical encounters at late times (see also Figure 6). Almost all WD–WD binaries that merged at late times are dynamically assembled. In contrast, more than 90% of the mergers at early times are from primordial binaries.

By assuming that all these WD–WD collisions and mergers produce MSPs directly (for alternatives, see Section 5), these two channels can also significantly increase the number of single MSPs in GCs (Table 2). The number of single MSPs formed in these channels is about half of those formed in TDEs (e.g., model 3a) or roughly comparable if WD–WD tidal capture is included or if masses are conserved during dynamical mass transfer in double WD binaries and merger-induced collapse of heavy WDs (e.g., models 6a and 3b). Thus the WD channels alone can also boost the number of single MSPs by a factor of ∼5–10.