<figure><img src="image_1.png" /></figure>

# The MAVERIC Survey: Dynamical Origin of Radio Sources in Galactic Globular Clusters

Yue Zhao<fnref n="1" />, Craig O. Heinke<fnref n="1" />, Laura Shishkovsky<fnref n="2" />, Jay Strader<fnref n="2" />, Laura Chomiuk<fnref n="2" />, Thomas J. Maccarone<fnref n="3" />, Arash Bahramian<fnref n="4" />, Gregory R. Sivakoff<fnref n="1" />, James C. A. Miller-Jones<fnref n="4" />, and Evangelia Tremou<fnref n="5" />

Received 2020 December 8; revised 2021 April 9; accepted 2021 April 21; published 2021 June 17

## Abstract

We investigate potential correlations between radio source counts (after background corrections) of 22 Galactic globular clusters (GCs) from the MAVERIC survey and the stellar encounter rates ($\Gamma$) and masses ($M$) of the GCs. Applying a radio luminosity limit of $L_{lim} = 5.0 \times 10^{27}$ erg s$^{-1}$, we take a census of radio sources in the core and those within the half-light radius of each cluster. By following a maximum likelihood method and adopting a simplified linear model, we find an unambiguous dependence of core radio source counts on $\Gamma$ and/or $M$ at 90% confidence, but no clear dependence of source counts within the half-light radius on either $\Gamma$ or $M$. Five of the identified radio sources in GC cores above our adopted limit are millisecond pulsars or neutron star X-ray binaries, the dependence of which on $\Gamma$ is well known, but another is a published black hole (BH) X-ray binary candidate, and 10 others are not identified. Accounting for these verified cluster members increases the significance of the correlation with $M$ and/or $\Gamma$ (to 99% confidence) for fits to core and half-light region source counts, while excluding a dependence on $\Gamma$ alone at 90% (core) and 68% (half-light) confidence. This is consistent with published dynamical simulations of GC BH interactions that argue $\Gamma$ will be a poor predictor of the distribution of accreting BHs in GCs. Future multiwavelength follow-up to verify cluster membership will enable stronger constraints on the dependence of radio source classes on cluster properties, promising a new view on the dynamics of BHs in GCs.

*Unified Astronomy Thesaurus concepts: Globular star clusters (656); Low-mass x-ray binary stars (939); Compact objects (288); Neutron stars (1108); Black holes (162); Stellar dynamics (1596)*

## 1. Introduction

Since the launch of early X-ray missions like Uhuru and OSO-7, it has been known that Galactic globular clusters (GCs) host an overabundance of X-ray sources (Katz 1975). These X-ray sources are ascribed to close binaries of various kinds, whose origin is closely related to the dynamically active cores of GCs. Specifically, the very dense cores of GCs facilitate many close dynamical encounters, producing close binaries through multiple possible channels (Clark 1975; Fabian et al. 1975; Sutantyo 1975; Hills 1976). The most well-studied population comprises low-mass X-ray binaries (LMXBs), where neutron stars (NSs) or black holes (BHs) accrete matter from (usually) a (near) main-sequence donor star (Lewin & Joss 1983; Grindlay et al. 1984). LMXBs were discovered to dominate the bright X-ray source population, with typical X-ray luminosities $L_X \sim 10^{36}-10^{38}$ erg s$^{-1}$ (Giacconi et al. 1974; Canizares & Neighbours 1975; Clark et al. 1975).

With the advent of the Chandra X-ray Observatory's superior sensitivity and angular resolution, a plethora of faint ($L_X \lesssim 10^{33}$ erg s$^{-1}$) sources were detected in many GCs (e.g., Bassa et al. 2004; Heinke et al. 2005, 2006; Bassa et al. 2008; Lu et al. 2009; Grindlay et al. 2001a, 2001b; Pooley et al. 2002a, 2002b; Zhao et al. 2019). Some of these faint sources are thought to be quiescent LMXBs (qLMXBs), which are generally $\lesssim 10^4$ times fainter than actively accreting LMXBs (e.g., Campana et al. 1998; Rutledge et al. 2002; Heinke et al. 2003). Closely related to LMXBs are the millisecond pulsars (MSPs); these are fast-spinning radio pulsars spun up by accreted matter from donor stars during a prior LMXB phase (Bhattacharya & van den

Heuvel 1991). MSPs emit thermal or nonthermal X-rays and also contribute to the faint source population, especially below $10^{31}$ erg s$^{-1}$ (Saito et al. 1997; Bogdanov et al. 2006). A majority of the faint X-ray population between $10^{31} < L_X < 10^{32}$ erg s$^{-1}$ comprises cataclysmic variables (CVs), which are accreting white dwarfs in close orbit with normal stars (Hertz & Grindlay 1983; Cool et al. 1995; Pooley et al. 2002a; Cohn et al. 2010; Rivera Sandoval et al. 2018). Finally, chromospherically active binaries (ABs), which are tidally locked close binaries of normal or evolved stars, are very common in GCs below $10^{31}$ erg s$^{-1}$ (Bailyn et al. 1990; Dempsey et al. 1993; Grindlay et al. 2001a; Bassa et al. 2004; Heinke et al. 2005). The LMXB and MSP populations, along with the brighter CVs, are closely correlated with the frequency of dynamical encounters, quantified by the stellar encounter rate (Verbunt & Hut 1987; Johnston et al. 1992; Heinke et al. 2003; Pooley et al. 2003; Verbunt 2003; Pooley & Hut 2006; Hui et al. 2010; Bahramian et al. 2013; Heinke et al. 2020). However, other X-ray sources are not of dynamical origin; most ABs descend from primordial binaries (Bassa et al. 2004, 2008; Huang et al. 2010; Cheng et al. 2018), as do many CVs (Davies 1997; Ivanova et al. 2006; Kong et al. 2006; Belloni et al. 2019).

Binaries containing BHs may not be distributed in the same way as those containing NSs, as the number of BHs remaining in different GCs, and their level of mixture with other stars, may vary in complicated ways. Early theoretical calculations suggested that dynamical interactions involving BHs would quickly expel all, or nearly all, BHs from GCs (Sigurdsson & Hernquist 1993;

<footnote n="1">¹Department of Physics, University of Alberta, 4-183 CCIS, Edmonton, AB T6G 2E1, Canada</footnote>

<footnote n="2">²Center for Data Intensive and Time Domain Astronomy, Department of Physics and Astronomy, Michigan State University, East Lansing, MI 48824, USA</footnote>

<footnote n="3">³Department of Physics and Astronomy, Texas Tech University, Box 41051, Lubbock, TX 79409-1051, USA</footnote>

<footnote n="4">⁴International Centre for Radio Astronomy Research, Curtin University, GPO Box U1987, Perth, WA 6845, Australia</footnote>

<footnote n="5">⁵LESIA, Observatoire de Paris, CNRS, PSL, SU/UPD, Meudon, France</footnote>
Kulkarni et al. 1993). However, observations of candidate BH X-ray binaries in extragalactic (Maccarone et al. 2007) and Galactic GCs (Strader et al. 2012), along with GC simulations that left numerous BHs in clusters (Morscher et al. 2013; Sippel & Hurley 2013; Heggie & Giersz 2014; Morscher et al. 2015), changed the prevailing wisdom. There is now solid dynamical evidence of three (noninteracting) BH binaries in the GC NGC 3201 (Giesers et al. 2018, 2019). Current GC simulations predict large numbers of BHs in clusters with relatively little mass segregation (Peuten et al. 2016; Askar et al. 2018; Weatherford et al. 2018, 2020), though the number of detectable BHs in binary systems with other stars is not predicted to correlate with the total number of BHs in the clusters (Arca Sedda et al. 2018; Askar et al. 2018; Kremer et al. 2018).

The MAVERIC (Milky way ATCA and VLA Exploration of Radio sources In Clusters) survey involves deep radio imaging of 50 Galactic GCs (see Tremou et al. 2018; Shishkovsky et al. 2020, Sh20 hereafter, and V. Tudor et al. 2021, in preparation), dedicated to the discovery of potential BH LMXBs and other exotic radio sources. The key motivation is that accreting BHs are more radio luminous for a given X-ray luminosity than other systems (Maccarone 2005; Migliari & Fender 2006), making deep radio surveys an excellent method to find candidate BHs. Thanks to the superior sensitivity and subarcsecond resolving power of the Karl G. Jansky Very Large Array (VLA) and the Australia Telescope Compact Array (ATCA), it has led to fruitful revelations of faint radio sources that are strong candidate BH LMXBs (Strader et al. 2012; Chomik et al. 2013; Miller-Jones et al. 2015; Bahramian et al. 2017, 2020; Shishkovsky et al. 2018), and strong candidates for MSPs not yet detected via radio pulsations (Bahramian et al. 2018; Zhao et al. 2020; Urquhart et al. 2020). The 5σ detection limit of the observations reaches as low as $S_ν \sim 10 \ \mu Jy \text{ beam}^{-1}$ ($S_ν$ is the flux density), covering different kinds of radio-emitting objects in the fields of the GCs.

Of the detected radio sources, we expect a considerable fraction to be background sources (e.g., active galactic nuclei (AGNs)), while the sources that are actually associated with the clusters are mostly BH or neutron star (NS) LMXBs and MSPs. Radio emission has been observed in both NS (Migliari & Fender 2006) and BH (Gallo et al. 2014, 2018) LMXBs, characterized by flat ($-0.5 < \alpha < 0$)<fnref n="6" /> to inverted ($\alpha > 0$) radio spectra, while MSPs are generally steep-spectrum radio sources (Bates et al. 2013; $\alpha \approx -1.4$ with unit standard deviation). CVs and ABs may also be active radio emitters. For example, radio emission has been observed from both nonmagnetic (e.g., Coppejans et al. 2015, 2016) and magnetic CVs (e.g., Barrett et al. 2017), generally fainter than that from LMXBs, but visible at GC distances (typically a few kiloparsecs) during the peaks of outbursts (e.g., Mooley et al. 2017). ABs have also been observed in the radio, typically only visible at kiloparsec distances during flares (e.g., Osten et al. 2000). One abnormal radio-emitting binary observed in the GC M10 with a low measured mass function is possibly an RS CVn type of AB, or an unusual CV, or a face-on BH LMXB (Shishkovsky et al. 2018).

It is then intriguing to investigate if there exists a dependence, similar to that revealed by X-ray observations, of radio source populations on GC dynamical parameters. In this work, we compare radio source counts of different GCs at the same luminosity cutoff, searching for correlations between the number

of radio sources versus stellar encounter rate ($\Gamma$) and/or GC mass ($M$). The paper is organized as follows: in Section 2, we describe the MAVERIC observations involved in this work, presenting the method we used to fit the data; in Section 3 we present results and discussions based on the results; and in Section 4, we draw conclusions.

## 2. Observations and Method

We use the VLA catalog of 5σ radio point sources from the MAVERIC survey (Sh20). The catalog compiles radio point sources for 25 GCs, reporting source coordinates and radio flux densities in low- and high-frequency bands at 5 GHz ($S_{low}$) and 7.2 GHz ($S_{high}$),<fnref n="7" /> respectively. For each GC, we count the numbers of radio sources within the core radius ($r_c$) and the half-light radius ($r_h$). A radio source is selected if its 5 GHz radio luminosity ($L_{low}$) at the host cluster distance is higher than a limiting luminosity ($L_{lim}$). We set $L_{lim} = 5.0 \times 10^{27} \text{ erg s}^{-1}$ to include most relatively bright radio sources, excluding only the distant GCs M2, M3, and M54, as $L_{lim}$ converts to a limiting flux ($S_{lim}$) below the <5σ detection limit of the catalog for those clusters. The source counts and relevant GC parameters are summarized in Table 1.

To fit the data, we follow the maximum likelihood method described in (Verbunt et al. 2008, V08 hereafter). We briefly outline the method here.

The number of radio sources observed within the core or half-light radius of a GC follows a Poisson distribution:

$$P(N, \mu) = \frac{\mu^N}{N!} e^{-\mu}, \qquad (1)$$

where $\mu$ is the expected number of sources, and $N$ is the observed number of sources. This formula applies to both the number of actual cluster members ($N_c$) and the number of background sources ($N_b$), with $\mu_c$ members or $\mu_b$ background sources expected.<fnref n="8" /> $\mu_b$ is calculated for each cluster using the normalized source counts from Sh20 (see Table 4 of Sh20), given in differential form $S_ν^{2.5} dN/dS_ν$, while applying the $S_{lim}$ of each cluster as the lower bound of our integration. We assume that the expected number of cluster members, $\mu_c$, is determined by $\Gamma$ and cluster mass ($M$). Given the low numbers of sources (Table 1) at our designated luminosity cutoff ($L_{lim}$), we follow a simplified linear model as in V08, viz.

$$\mu_c = a\Gamma + bM, \qquad (2)$$

where $a$ and $b$ are positive coefficients, while for convenience, we renormalized $\Gamma$ and $M$, taken from Table 1, to fractions of the $\Gamma$ of M62 and the $M$ of M14. The form of the model is based on the assumption that one expects more LMXBs and MSPs in GCs with higher $\Gamma$ and/or $M$. The positive correlation with $\Gamma$ has been tested by using a census of X-ray sources from archival observations (e.g., Pooley & Hut 2006) and radio timing surveys (Hui et al. 2010).

The likelihood function is then the multiplication of the joint probability, $P(N_c, \mu_c)P(N_b, \mu_b)$, over all GCs in our sample,

<footnote n="6">⁶ $\alpha$ is the radio spectral index ($\alpha$) defined by $S_ν \propto \nu^\alpha$.</footnote>

<footnote n="7">⁷ Note that these frequencies are average central frequencies; the actual values may be slightly different between clusters.</footnote>

<footnote n="8">⁸ Note that here $\mu_c$ and $\mu_b$ are shorthands for the expected numbers of members and background sources; in our application, they can represent the numbers in the core (indicated by a “core” subscript) or within the half-light region (indicated by a “half” subscript) as listed in Table 1.</footnote>
<table><caption>**Table 1**
Radio Source Counts in Different GCs</caption><thead><tr><th>GC Name</th><th>Distance (kpc)</th><th>r<sub>c</sub> (')</th><th>r<sub>h</sub> (')</th><th>S<sub>lm</sub> (μJy)</th><th>N<sub>core</sub></th><th>N<sub>half</sub></th><th>N<sub>m,core</sub></th><th>N<sub>m,half</sub></th><th>μ<sub>b,core</sub></th><th>μ<sub>b,half</sub></th><th>Γ</th><th>M (× 10<sup>5</sup>M<sub>☉</sub>)</th></tr></thead><tbody><tr><td>M62*</td><td>6.7</td><td>0.22</td><td>0.92</td><td>18.6</td><td>1</td><td>3</td><td>1</td><td>2</td><td>0.07</td><td>1.16</td><td>1670<sup>+710</sup><sub>-570</sub></td><td>6.90 ± 0.05</td></tr><tr><td>NGC 6440*</td><td>8.5</td><td>0.14</td><td>0.48</td><td>11.6</td><td>1</td><td>2</td><td>1</td><td>1</td><td>0.07</td><td>0.78</td><td>1400<sup>+630</sup><sub>-477</sub></td><td>5.33 ± 0.65</td></tr><tr><td>M28*</td><td>5.5</td><td>0.24</td><td>1.97</td><td>27.6</td><td>1</td><td>4</td><td>1</td><td>1</td><td>0.05</td><td>3.48</td><td>648<sup>+85</sup><sub>-91</sub></td><td>2.84 ± 0.15</td></tr><tr><td>M30</td><td>8.6</td><td>0.06</td><td>1.03</td><td>11.3</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0.01</td><td>3.58</td><td>324<sup>+124</sup><sub>-81</sub></td><td>1.39 ± 0.06</td></tr><tr><td>M92</td><td>8.9</td><td>0.26</td><td>1.02</td><td>10.6</td><td>0</td><td>3</td><td>0</td><td>0</td><td>0.23</td><td>3.51</td><td>270<sup>+30</sup><sub>-29</sub></td><td>3.12 ± 0.04</td></tr><tr><td>M19</td><td>8.2</td><td>0.43</td><td>1.32</td><td>12.4</td><td>1</td><td>3</td><td>0</td><td>0</td><td>0.40</td><td>3.74</td><td>200<sup>+67</sup><sub>-39</sub></td><td>6.57 ± 0.33</td></tr><tr><td>M5</td><td>7.7</td><td>0.44</td><td>1.77</td><td>14.1</td><td>0</td><td>2</td><td>0</td><td>0</td><td>0.42</td><td>6.73</td><td>164<sup>+39</sup><sub>-30</sub></td><td>3.68 ± 0.06</td></tr><tr><td>M9</td><td>7.8</td><td>0.45</td><td>0.96</td><td>13.7</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0.43</td><td>1.98</td><td>131<sup>+59</sup><sub>-42</sub></td><td>3.21 ± 0.25</td></tr><tr><td>M14</td><td>9.3</td><td>0.79</td><td>1.30</td><td>9.7</td><td>4</td><td>14</td><td>0</td><td>0</td><td>2.11</td><td>5.71</td><td>124<sup>+32</sup><sub>-30</sub></td><td>7.39 ± 0.37</td></tr><tr><td>NGC 6304</td><td>5.9</td><td>0.21</td><td>1.42</td><td>24.0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0.04</td><td>1.81</td><td>123<sup>+54</sup><sub>-22</sub></td><td>2.11 ± 0.18</td></tr><tr><td>NGC 6325</td><td>6.5</td><td>0.03</td><td>0.63</td><td>19.8</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0.001</td><td>0.54</td><td>118<sup>+45</sup><sub>-46</sub></td><td>0.76 ± 0.13</td></tr><tr><td>NGC 6544</td><td>3.0</td><td>0.05</td><td>1.21</td><td>92.9</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0.00</td><td>0.33</td><td>111<sup>+68</sup><sub>-36</sub></td><td>1.15 ± 0.11</td></tr><tr><td>M22</td><td>3.1</td><td>1.33</td><td>3.36</td><td>86.9</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0.40</td><td>2.54</td><td>78<sup>+31</sup><sub>-26</sub></td><td>4.09 ± 0.04</td></tr><tr><td>M13</td><td>7.6</td><td>0.62</td><td>1.69</td><td>14.5</td><td>1</td><td>5</td><td>1</td><td>0</td><td>0.83</td><td>6.13</td><td>69<sup>+18</sup><sub>-15</sub></td><td>4.69 ± 0.20</td></tr><tr><td>NGC 6760</td><td>7.4</td><td>0.34</td><td>1.27</td><td>15.3</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0.25</td><td>3.46</td><td>57<sup>+27</sup><sub>-19</sub></td><td>2.55 ± 0.30</td></tr><tr><td>NGC 6539*</td><td>7.8</td><td>0.38</td><td>1.70</td><td>13.7</td><td>2</td><td>5</td><td>1</td><td>1</td><td>0.04</td><td>6.20</td><td>42<sup>+29</sup><sub>-15</sub></td><td>2.61 ± 0.31</td></tr><tr><td>M10</td><td>4.4</td><td>0.77</td><td>1.95</td><td>43.1</td><td>0</td><td>3</td><td>0</td><td>0</td><td>0.25</td><td>1.59</td><td>31<sup>+5</sup><sub>-4</sub></td><td>1.89 ± 0.04</td></tr><tr><td>NGC 6712</td><td>8.0</td><td>0.76</td><td>1.33</td><td>13.1</td><td>2</td><td>2</td><td>1</td><td>1</td><td>1.24</td><td>3.80</td><td>31<sup>+5</sup><sub>-7</sub></td><td>1.19 ± 0.08</td></tr><tr><td>M4</td><td>1.8</td><td>1.16</td><td>4.33</td><td>258.0</td><td>0</td><td>2</td><td>0</td><td>0</td><td>0.11</td><td>1.49</td><td>27<sup>+12</sup><sub>-10</sub></td><td>0.93 ± 0.02</td></tr><tr><td>M12</td><td>5.2</td><td>0.79</td><td>1.77</td><td>30.9</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0.38</td><td>1.93</td><td>13<sup>+5</sup><sub>-4</sub></td><td>0.87 ± 0.03</td></tr><tr><td>M107</td><td>6.1</td><td>0.56</td><td>1.73</td><td>22.5</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0.28</td><td>2.69</td><td>6.8<sup>+2.3</sup><sub>-1.7</sub></td><td>0.81 ± 0.05</td></tr><tr><td>M55</td><td>5.7</td><td>1.80</td><td>2.83</td><td>25.7</td><td>2</td><td>5</td><td>0</td><td>0</td><td>2.91</td><td>7.19</td><td>3.2<sup>+1.4</sup><sub>-1.0</sub></td><td>1.88 ± 0.07</td></tr></tbody></table>

Note. The distances are from Tremou et al. (2018); the core radii ($r_c$) and half-light radii ($r_h$) are from Harris (1996, 2010 edition); the $\Gamma$'s are from Bahramian et al. (2013); the GC masses ($M$) are from Baumgardt & Hilker (2018, second version; https://people.smp.uq.edu.au/HolgerBaumgardt/globular/); and $N_{\text{core}}$ and $N_{\text{half}}$ are the observed numbers of radio sources within the core and the half-light region, with $L > 5 \times 10^{27} \text{ erg s}^{-1}$. $N_{\text{core}}$ and $N_{\text{half}}$ are numbers of confirmed members within the core radius and the half-light radius. GCs indicated with * have a significant radio source excess in the core over the expected background.

i.e.,

$$ \mathcal{L} = \prod_i P(N_{c,i}, \mu_{c,i}) P(N_{b,i}, \mu_{b,i}), \quad (3) $$

where $i$ indexes the GCs. The best-fitting model is given by a combination of $a$ and $b$ that maximizes $\mathcal{L}$.

We set up a grid of $a$ and $b$ values ranging from 0 to 4, with a spacing of 0.02 in both $a$ and $b$, and for each pair of $a$ and $b$, we generate 1000 random Poisson realizations. Each realization draws a random integer as per a Poisson distribution given $\mu_c$, which is assigned to $N_c$, $N_b$ is then calculated by subtracting $N_c$ from the observed number of radio sources as listed in the sixth and the seventh column of Table 1. To use the Poisson probabilities, we keep $N_c$ below the observed number of sources ($N_{\text{core}}$ or $N_{\text{half}}$), avoiding negative $N_b$ values; this is done by setting the probability $P(N_b, \mu_b) = 0$ when $N_b$ is negative while keeping the total number of realizations (1000) for each pair of $a$ and $b$ unchanged. Zero probability leads to zero likelihood, so equivalently these realizations are excluded from the maximization of the combined likelihood ($\mathcal{L}$).

We also perform a somewhat more constrained fit by keeping $N_c$ to be at least the number of confirmed members in the cores ($N_{\text{m,core}}$, the eighth column in Table 1) or half-light regions ($N_{\text{m,half}}$, the ninth column of Table 1). The members include confirmed core MSPs and LMXBs listed in Table 2 and PSR J1701−3006A in M62, which is outside the core (Lynch et al. 2012). Similarly, realizations with $N_c < N_{\text{m,core}}$ or $N_{\text{m,half}}$ are excluded from the fit by setting $P(N_c, \mu_c) = 0$. For simplicity, this fit is referred to as the “constrained fit”—to differentiate it from the “unconstrained fit” where we only

apply upper bounds on $N_c$ ($N_c$ less than or equal to the observed number of sources).

To plot the result and compute confidence contours, we follow the definition of $Z$ in V08:

$$ Z \equiv -2[\log(\mathcal{L}) - \log(\mathcal{L})_{\text{max}}]. \quad (4) $$

Confidence contours are calculated assuming that $Z$ follows a $\chi^2$ distribution with one degree of freedom, so the best fit ($\mathcal{L} = \mathcal{L}_{\text{max}}$) corresponds to $Z=0$, while the 68%, 90%, and 99% confidence intervals are where $\Delta Z = 1.00, 2.71, \text{and } 6.63$, respectively.

### 3. Results and Discussion

The resulting distributions of $Z$ for unconstrained and constrained fits are presented in Figures 1 and 2, overplotted with confidence contours. The unconstrained fit to source counts in the core gives the best-fitting $a = 0.52$ and $b = 0.58$, excluding noncorrelation ($a = b = 0$) at the 90% confidence level; whereas the best fit to source counts within the half-light region suggests no correlation: $a = b = 0$.

The constrained fit to source counts in the core consistently favors a correlation (best fit: $a = 0.10$ and $b = 1.09$), now at the 99% confidence level. Furthermore, the entire $a$-axis lies outside the 90% confidence contour, indicating a dependence of source counts on GC mass ($b \neq 0$) at 90% confidence. We note that the constrained fit to source counts within the half-light region (best fit: $a = 0.86$, and $b = 0.56$) also excludes noncorrelation at 99% confidence, in contrast to the unconstrained fit, while it rules out the $\Gamma$-only dependence ($b = 0$) at a lower (68%) confidence level.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** 68% (solid red), 90% (dashed green), and 99% (dashed-dotted gray) contours in the *ab* plane for fits to radio source counts in cluster cores (left) and half-light regions (right). The best-fitting values of *a* and *b* are indicated with a black cross in each panel. The color bar presents the intensity of Z. For source counts in the core, *a* ≠ 0 or *b* ≠ 0 at the 90% level, while *a* = *b* = 0 is consistent with the data for sources within the half-light regions.</figcaption></figure>

<table><caption>**Table 2**
Radio Sources in the Core with $L_{\text{low}} > 5 \times 10^{27} \text{ erg s}^{-1}$</caption><thead><tr><th>Source ID</th><th>L<sub>low</sub><br>(× 10<sup>27</sup> erg s<sup>-1</sup>)</th><th>S<sub>low</sub>/rms<br>(σ)</th><th>α<br>(S<sub>ν</sub> ∝ ν<sup>α</sup>)</th><th>Notes</th><th>References</th></tr></thead><tbody><tr><td>M9-VLA24</td><td>5.4 ± 0.6</td><td>8.6</td><td>-0.9<sup>+0.6</sup><sub>-0.6</sub></td><td>...</td><td>...</td></tr><tr><td>M13-VLA21</td><td>5.0 ± 0.7</td><td>7.3</td><td>-1.3<sup>+0.7</sup><sub>-0.8</sub></td><td>PSR B1639 + 36A</td><td>K91, W20</td></tr><tr><td>M14-VLA8</td><td>24.2 ± 1.0</td><td>26.1</td><td>-1.9<sup>+0.3</sup><sub>-0.3</sub></td><td>...</td><td>...</td></tr><tr><td>M14-VLA11</td><td>20.0 ± 1.0</td><td>21.4</td><td>-1.2<sup>+0.3</sup><sub>-0.3</sub></td><td>...</td><td>...</td></tr><tr><td>M14-VLA15</td><td>18.6 ± 1.0</td><td>20.0</td><td>-0.2<sup>+0.2</sup><sub>-0.2</sub></td><td>...</td><td>...</td></tr><tr><td>M14-VLA45</td><td>6.3 ± 1.0</td><td>6.8</td><td>&lt;0.0</td><td>...</td><td>...</td></tr><tr><td>M19-VLA34</td><td>5.1 ± 0.9</td><td>5.5</td><td>&lt;0.6</td><td>...</td><td>...</td></tr><tr><td>M28-VLA3</td><td>24.3 ± 0.5</td><td>53.7</td><td>-2.2<sup>+0.1</sup><sub>-0.1</sub></td><td>PSR B1821-24A</td><td>L87, F88, C04, R04, B11, J13</td></tr><tr><td>M55-VLA6</td><td>8.7 ± 0.4</td><td>19.6</td><td>-1.2<sup>+0.3</sup><sub>-0.3</sub></td><td>...</td><td>...</td></tr><tr><td>M55-VLA15</td><td>5.1 ± 0.5</td><td>11.3</td><td>-2.3<sup>+0.7</sup><sub>-0.7</sub></td><td>...</td><td>...</td></tr><tr><td>M62-VLA1</td><td>6.0 ± 1.0</td><td>7.0</td><td>-0.4<sup>+0.6</sup><sub>-0.5</sub></td><td>A BH candidate</td><td>C13</td></tr><tr><td>NGC 6440-VLA6</td><td>33.2 ± 1.3</td><td>28.4</td><td>-2.0<sup>+0.2</sup><sub>-0.2</sub></td><td>PSR B1745-20</td><td>L96, F08</td></tr><tr><td>NGC 6539-VLA11</td><td>9.4 ± 0.7</td><td>9.4</td><td>-0.7<sup>+0.7</sup><sub>-0.8</sub></td><td>PSR B1802-07</td><td>D93, T93, T99</td></tr><tr><td>NGC 6539-VLA12</td><td>5.2 ± 0.6</td><td>8.5</td><td>-0.4<sup>+0.7</sup><sub>-0.7</sub></td><td>...</td><td>...</td></tr><tr><td>NGC 6712-VLA7</td><td>33.8 ± 0.9</td><td>38.4</td><td>0.1<sup>+0.1</sup><sub>-0.1</sub></td><td>4U 1850-087</td><td>S76, A93, H96, S06</td></tr><tr><td>NGC 6712-VLA9</td><td>25.0 ± 0.9</td><td>28.3</td><td>-0.6<sup>+0.1</sup><sub>-0.1</sub></td><td>...</td><td>...</td></tr></tbody></table>

**Note.** rms noises are adapted from Sh20. α is the spectral index from Sh20. References: A93: Anderson et al. (1993), B11: Bogdanov et al. (2011), C04: Cognard & Backer (2004), C13: Chomik et al. (2013), D93: D'Amico et al. (1993), F88: Foster et al. (1988), F08: Freire et al. (2008), H96: Homer et al. (1996), J13: Johnson et al. (2013), K91: Kulkarni et al. (1991), L87: Lyne et al. (1987), L96: Lyne et al. (1996), R04: Rutledge et al. (2004), S76: Swank et al. (1976), S06: Sidoli et al. (2006), T93: Thorsett et al. (1993), T99: Thorsett & Chakrabarty (1999), W20: Wang et al. (2020).

For both the constrained and unconstrained fits, we generate 1000 random data sets based on the best-fitting *a* and *b* and calculate the corresponding combined likelihood *L* for each data set. We note that the maximum likelihoods corresponding to the best fits are greater than those for all random data sets; the best fits are therefore appropriate for modeling the radio source counts.

Looking at Table 1, we note that a few GCs (indicated with an *) show significant core source excesses over the predicted numbers of background sources (6th column versus 10th column). Statistically, the probability of observing 1 while expecting 0.04–0.07 background sources (Table 1) is around

4%–7% (Equation (1)), so the excess is not very likely given by background sources. In Figure 3, we plot the model predicted counts using the best-fitting parameters versus the observed counts in the core for each GC. We also overplot error bars that indicate uncertainties propagated from the fitting parameters (*a* and/or *b*) and from the corresponding GC parameters (Γ and/or *M*).

Source counts within the half-light region have no clear dependence on either Γ or *M* (with the best-fitting *a* and *b* = 0), when $N_c$ is only limited by $N_{\text{half}}$. However, when confirmed members are specified, the result excludes nondependence ($a = b = 0$) at 99% confidence. Of the six GCs that have confirmed
<figure><img src="image_3.png" /><figcaption>**Figure 2.** The marker and the contours are the same as those of Figure 1, but this fit required $N_c$ to be equal to or greater than the number of confirmed members in the core (left) or within the half-light region (right).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Observed core source counts (black empty squares) and observed half-light region source counts (red empty diamonds) vs. predicted source counts in the core (black filled squares) and in the half-light region (red filled diamonds) based on the constrained fits plotted for all 22 GCs. The error bars on the observed numbers indicate 84% confidence upper and lower limits (Gehrels 1986); for GCs with zero counts, we only plot the upper bar. The uncertainties in predicted numbers are propagated from uncertainties in $\Gamma$ and $M$, and the 68% confidence ranges of a and b. The black dashed and red dotted lines in each panel mark the number of expected background sources in the core ($\mu_{b,core}$) and in the half-light region ($\mu_{b,half}$), respectively.</figcaption></figure>

members, three (M62, NGC 6440, and M28) are large in both $\Gamma$ and $M$ (Table 1); specifying confirmed members for these GCs favors a positive correlation between $\mu_c$ and $\Gamma$ and/or $M$, although the observed source counts ($N_{half}$) are generally consistent with the predicted numbers of background sources ($\mu_{b,half}$) within the 90% confidence level derived for Poisson statistics (Gehrels upper and lower limits; see Gehrels 1986).

One exception is M14—the observed half-light region source count (14) is consistent with the background estimate (5.71) at a higher, 98%, confidence level, considering the Gehrels upper limit. We note that M14 has the lowest $S_{lim}$ of ≈9.2 µJy, and for all GCs but M14, $S_{lim}$ is above the very first flux bin of Sh20's radio source counts (between 7.90 µJy and 10.61 µJy; see Table 4 of Sh20). For M14, part of this bin goes below the 5σ flux limit of the observation (≈9 µJy). Including this flux bin in the integration will incorporate faint sources, between 9 and 10.61 µJy, and so gives a somewhat higher background estimate (7.82) for M14. However, the first flux bin also misses faint sources between 7.9 and 9 µJy. It is hence clear that the expected number of background sources ($\mu_{b,half}$) might have been underestimated for M14 even when the first

flux bin is considered. We run our constrained fits without M14 and note that the best fit to core radio source counts consistently excludes $a=b=0$ at 99% confidence, with the maximum likelihood at $a=0.08$ and $b=1.03$, and the best fit to the half-light region counts also excludes $a=b=0$ at 99% confidence, giving $a=0.93$ and $b=0.46$.

The 22 GCs listed in Table 1 host a total of 16 sources in their cores. We further investigate all these sources by matching their radio positions to those measured by timing observations (Table 2, and references therein). We find four known MSPs, three of which are in M28, NGC 6440, and NGC 6539, GCs that have core source excesses over the background (Table 1), and one in M13, which is marginally above $L_{lim}$ (Table 2). Two more sources are known cluster members—a BH candidate in M62 (Chomiuk et al. 2013) and an NS LMXB in NGC 6712 (Swank et al. 1976).

These sources are all confirmed cluster members, with 4 hosted by GCs that have a core radio source excess, namely M28, M62, NGC 6440, and NGC 6539. Three of these are GCs with high mass and $\Gamma$ (Table 1), in agreement with the positive correlation shown by the fit.
NGC 6539 possesses a clear core source excess (2 observed versus 0.04 expected) while having a low Γ of 42 and an intermediate $M (=2.6 \times 10^5 M_\odot)$ versus median mass of $2.8 \times 10^5 M_\odot$). There are two core sources, VLA11 and VLA12. VLA11 is positionally consistent with the known MSP, PSR B1745−20 (D'Amico et al. 1993), while VLA12 has an inverted-to-steep index ($\alpha = -0.4 \pm 0.7$, Table 2), which can overlap with multiple scenarios. The probability of finding a background source in the core is ≈8%, not very likely but not ruled out. Considering the low Γ, it is possible that VLA12 is a primordial cluster member. For example, short (~hours) flares observed in some short-period RS CVn ABs can be detected in the radio at kpc distances (e.g., Osten et al. 2000).

MSPs in cluster cores are expected to scale directly with $\Gamma$, as they are produced by LMXBs that scale with the stellar encounter rate (e.g., Heinke et al. 2003; Pooley & Hut 2006; Hui et al. 2010; Bahramian et al. 2013). The presence of four confirmed MSPs among the cluster radio sources indicates that at least some cluster radio sources must be dynamically formed. However, our analysis including known cluster members requires a contribution from cluster mass, not just stellar encounter rate. It is therefore of great interest that the relative numbers of radio-emitting BH LMXBs are not predicted to scale directly with relative stellar encounter rates $\Gamma$ (which are calculated for visible stars, not the BH subsystems). The numbers of these BH binaries are instead predicted to depend in complicated ways on the total number of BHs in each cluster and on the interactions of these BH populations with other cluster stars (Arca Sedda et al. 2018; Weatherford et al. 2018). Thus, our work gives tentative support to the idea that a portion of the cluster radio sources are BH LMXBs, distributed in a complicated way among clusters, although we cannot rule out a contribution by other kinds of sources.

Another interesting test involves comparing the radio spectral indices of our sample with the known distribution of spectral indices of MSPs. In Figure 4, we present normalized histograms of spectral indices ($\alpha$) for pulsars from the ATNF catalog (Manchester et al. 2005, and references therein), sources from the MAVERIC catalog that are outside of the half-light radii and thus are mostly background AGNs, and for radio sources in Table 2 (excluding sources with unconstrained $\alpha$ values). We first compare equal-size samples of MSPs ($P \le 20$ ms) and normal pulsars from the ATNF catalog, and note that there is no significant difference between the distribution of their spectral indices ($p$-value of 0.25 for a two-sample Anderson–Darling test). We can therefore randomly draw pulsars from the whole ATNF catalog to form a sample the same size as that of the AGNs. The core sources seem to follow a bimodal distribution, containing a group of steep ($-2.5 \lesssim \alpha \lesssim -1.2$) sources, and a group of relatively flatter ($-0.8 \lesssim \alpha \lesssim 0.2$) sources. The latter significantly deviates from the observed pulsar distribution and exceeds the AGN distribution around $-0.5 \lesssim \alpha \lesssim -0.2$, which could be partially contributed by LMXBs. We performed a two-sample Anderson–Darling test (Scholz & Stephens 1987) comparing the core source distribution with that of the ATNF pulsars and that of AGNs. A test comparing the core sources with ATNF pulsars rejects the null hypothesis that the core sources and ATNF pulsars are drawn from the same population at a level more significant than 0.1%. A test comparing core sources with the AGN sample rejects the null hypothesis at 3% significance.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Top: A histogram showing the distribution of the spectral index ($\alpha$) for core radio sources with constrained $\alpha$ values as listed in Table 2 (orange). Bottom: Histograms of $\alpha$ for pulsars from the ATNF catalog (hatched green), AGNs (hatched red), and pulsar and AGN combined (solid black contour). The counts in each bin are normalized by dividing the total counts of each sample and the bin width. The core sources divide into steep-index and flatter-index groups; the flatter group is more numerous than expected from a pulsar-only distribution. Despite the fact that AGNs can partly account for the flatter group, there is an overabundance of core sources around $-0.5 \lesssim \alpha \lesssim -0.2$.</figcaption></figure>

Finally, a test comparing the core sources with the AGN and pulsar samples combined rejects the null hypothesis at a level more significant than 0.1%. In fact, the flatter group contains six sources, of which one is an NS LMXB, one is the known BH candidate in M62, and one is the known MSP in NGC 6539; the other three sources could be either background AGNs or cluster members.

Our analyses have shown that additional membership information results in tighter constraints on the parameters for both source counts in the core and those in the half-light region. In this regard, future careful multiwavelength follow-up will therefore be necessary to separate background sources from cluster members, allowing investigation of radial trends in
the radio source properties. For example, proper motion analysis of potential optical counterparts, or with radio VLBI (Tetarenko et al. 2016), can be effective in excluding background sources. Detailed investigation at other wavelengths is beyond the scope of this paper but will be addressed in future work.

There are several potential factors that could affect our final results: uncertainties in $\Gamma$ or mass of clusters, differences in sensitivity to radio sources across the cluster cores, and unusual cluster histories. Uncertainties in $\Gamma$ can be up to 70% of the central values (Bahramian et al. 2013) but are not included in our aforementioned fits. To address this, we fit *a* and *b* simultaneously to the same observed counts but with $\Gamma$ and *M* randomly drawn from the error ranges for each cluster, and observe the distribution of the best-fitting values and the corresponding confidence contours. We ran a total of 50 random constrained fits and find the point $a=b=0$ is outside of the 99% confidence contour in all of the fits; running 50 random constrained fits for source counts in the half-light regions also excludes $a=b=0$ at 99% confidence, so the correlation holds for both core and half-light region source counts, regardless of uncertainties in $\Gamma$ and *M*.

Another factor to consider is the increase in noise in the radio images with increasing off-axis distance. To first order, the rms noise scales as the reciprocal of the primary beam sensitivity. As a result, sources at the edges of the beam might have a flux above the average $5\sigma$ level but might have been missed because they are below the local $5\sigma$ level. This effect is expected to be minor in our analyses of the core sources, as most of the clusters have core sizes much smaller than the half-power primary beam size ($\approx 8'2$ at 5.5 GHz). The largest core, in the cluster M55, is $1'8$ in radius. If we made the assumption that all of our core sources (Table 2) are located $1'8$ from the beam center, we would expect the sensitivity to drop to $\approx 88\%$ of the center value (assuming the primary beam follows a Gaussian sensitivity curve); this corresponds to an increase in noise by a factor of $\approx 1.1$, so the detection procedure might miss sources between $5\sigma$ and $5.5\sigma$. Only one (of 16) of our core sources, M19-VLA34, falls in this flux range; we thus estimate that $<6\%$ of the core sources above our stated flux limit might have been missed by our observations. Each GC has one to three core sources, corresponding to 0.1–0.2 missing sources. Therefore, the effect of noise variation is minor in our analyses.

A final concern is whether the unusual histories of these GCs could mean that the current mass and stellar interaction rate do not represent their values during the period in which X-ray binaries or MSPs were formed. The prototype for this argument is the GC NGC 6712, which was shown to have lost >99% of its mass, leaving preferentially high-mass stars and binaries in the core (de Marchi et al. 1999; Andreuzzi et al. 2001). This cluster was probably much more massive and dense when the two detected radio sources (one a known NS LMXB) were formed, explaining their presence in an apparently low-mass, low-density cluster (see Ferraro et al. 2000). However, running our fits excluding NGC 6712 gives best-fitting $a = 0.32$ and $b = 0.78$ for the core source counts, and $a = 0.97, b = 0.42$ for the half-light region counts, while both fits consistently exclude $a=b=0$ at 99% confidence (when lower limits on $N_c$ are included). Although similar cluster mass loss is likely to have affected other clusters (see, e.g., Moreno et al. 2014), this effect has not removed the detectable correlation of stellar

interactions on X-ray sources (Pooley & Hut 2006; Bahramian et al. 2013), and therefore we suspect its effects on radio sources are likely also modest.

### 4. Conclusion

We investigate linear correlations of radio source counts with the encounter rate ($\Gamma$) and mass (*M*) for a total of 22 GCs from the MAVERIC survey. Including information about confirmed cluster members in our analysis, the fit to source counts in the core rules out no correlation with *M* or $\Gamma$ at 99% confidence and excludes a dependence on $\Gamma$ alone at 90% confidence. The fit to source counts in the half-light region also rules out no correlation at 99% confidence, but just excludes a $\Gamma$-only dependence at 68% confidence.

The histogram of spectral indices of our radio sources is inconsistent with that of pulsars and/or AGN alone, indicating another component with flatter spectra is present. Our findings are intriguingly consistent with the expectations of detectable radio sources—some will be MSPs (produced dynamically, and thus correlated with $\Gamma$), and some may be BH LMXBs, the numbers of which are unlikely to be directly correlated with $\Gamma$ (Arca Sedda et al. 2018; Kremer et al. 2018).

Our analyses also suggest that membership information can lead to significantly tighter constraints on the parameters. Future secure determinations of the nature of individual systems, obtained with careful follow-up, will therefore allow us to test predictions of the distribution of BH binaries in globular clusters in more detail. More detailed analyses will also use the full luminosity range of detected radio sources in each cluster and (with sufficient follow-up efforts) will study their radial distributions.

We thank the anonymous referee for helpful comments. C.O.H. and G.R.S. are supported by NSERC Discovery Grants RGPIN-2016-04602 and RGPIN-2016-06569, respectively. J.S. is supported by NSF grants AST-1308124 and AST-1514763, and a Packard Fellowship. J.C.A.M.-J is the recipient of an Australian Research Council Future Fellowship (FT140101082), funded by the Australian government. The MAVERIC source catalog used in this work is derived from observations made by the Karl G. Jansky Very Large Array under the National Radio Astronomy Observatory. The National Radio Astronomy Observatory is a facility of the National Science Foundation (NSF) operated under cooperative agreement by Associated Universities, Inc. (AUI).

### ORCID iDs

Yue Zhao https://orcid.org/0000-0002-9547-8677

Craig O. Heinke https://orcid.org/0000-0003-3944-6109

Laura Shishkovsky https://orcid.org/0000-0003-0286-7858

Jay Strader https://orcid.org/0000-0002-1468-9668

Laura Chomiuk https://orcid.org/0000-0002-8400-3705

Arash Bahramian https://orcid.org/0000-0003-2506-6041

Gregory R. Sivakoff https://orcid.org/0000-0001-6682-916X

James C. A. Miller-Jones https://orcid.org/0000-0003-3124-2814

Evangelia Tremou https://orcid.org/0000-0002-4039-6703

### References

Anderson, S. F., Margon, B., Deutsch, E. W., & Downes, R. A. 1993, AJ, 106, 1049

Andreuzzi, G., De Marchi, G., Ferraro, F. R., et al. 2001, A&A, 372, 851
Arca Sedda, M., Askar, A., & Giersz, M. 2018, MNRAS, 479, 4652

Askar, A., Arca Sedda, M., & Giersz, M. 2018, MNRAS, 478, 1844

Bahramian, A., Heinke, C. O., Sivakoff, G. R., & Gladstone, J. C. 2013, ApJ, 766, 136

Bahramian, A., Heinke, C. O., Tudor, V., et al. 2017, MNRAS, 467, 2199

Bahramian, A., Strader, J., Chomiuk, L., et al. 2018, ApJ, 864, 28

Bahramian, A., Strader, J., Miller-Jones, J. C. A., et al. 2020, ApJ, 901, 57

Bailyn, C. D., Grindlay, J. E., & Garcia, M. R. 1990, ApJL, 357, L35

Barrett, P. E., Dieck, C., Beasley, A. J., Singh, K. P., & Mason, P. A. 2017, AJ, 154, 252

Bassa, C., Pooley, D., Homer, L., et al. 2004, ApJ, 609, 755

Bassa, C. G., Pooley, D., Verbunt, F., et al. 2008, A&A, 488, 921

Bates, S. D., Lorimer, D. R., & Verbiest, J. P. W. 2013, MNRAS, 431, 1352

Baumgardt, H., & Hilker, M. 2018, MNRAS, 478, 1520

Belloni, D., Giersz, M., Rivera Sandoval, L. E., Askar, A., & Ciecielåg, P. 2019, MNRAS, 483, 315

Bhattacharya, D., & van den Heuvel, E. P. J. 1991, PhR, 203, 1

Bogdanov, S., Grindlay, J. E., Heinke, C. O., et al. 2006, ApJ, 646, 1104

Bogdanov, S., van den Berg, M., Servillat, M., et al. 2011, ApJ, 730, 81

Campana, S., Stella, L., Mereghetti, S., et al. 1998, ApJL, 499, L65

Canizares, C. R., & Neighbours, J. E. 1975, ApJL, 199, L97

Cheng, Z., Li, Z., Xu, X., & Li, X. 2018, ApJ, 858, 33

Chomiuk, L., Strader, J., Maccarone, T. J., et al. 2013, ApJ, 777, 69

Clark, G. W. 1975, ApJL, 199, L143

Clark, G. W., Markert, T. H., & Li, F. K. 1975, ApJL, 199, L93

Cognard, I., & Backer, D. C. 2004, ApJL, 612, L125

Cohn, H. N., Lugger, P. M., Couch, S. M., et al. 2010, ApJ, 722, 20

Cool, A. M., Grindlay, J. E., Cohn, H. N., Lugger, P. M., & Slavin, S. D. 1995, ApJ, 439, 695

Coppejans, D. L., Körding, E. G., Miller-Jones, J. C. A., et al. 2015, MNRAS, 451, 3801

Coppejans, D. L., Körding, E. G., Miller-Jones, J. C. A., et al. 2016, MNRAS, 463, 2229

D'Amico, N., Bailes, M., Lyne, A. G., et al. 1993, MNRAS, 260, L7

Davies, M. B. 1997, MNRAS, 288, 117

de Marchi, G., Leibundgut, B., Paresce, F., & Pulone, L. 1999, A&A, 343, L9

Dempsey, R. C., Linsky, J. L., Fleming, T. A., & Schmitt, J. H. M. M. 1993, ApJS, 86, 599

Fabian, A. C., Pringle, J. E., & Rees, M. J. 1975, MNRAS, 172, 15

Ferraro, F. R., Paltrinieri, B., Paresce, F., & De Marchi, G. 2000, ApJL, 542, L29

Foster, R. S., Backer, D. C., Taylor, J. H., & Goss, W. M. 1988, ApJL, 326, L13

Freire, P. C. C., Ransom, S. M., Bégin, S., et al. 2008, ApJ, 675, 670

Gallo, E., Degenaar, N., & van den Eijnden, J. 2018, MNRAS, 478, L132

Gallo, E., Miller-Jones, J. C. A., Russell, D. M., et al. 2014, MNRAS, 445, 290

Gehrels, N. 1986, ApJ, 303, 336

Giacconi, R., Murray, S., Gursky, H., et al. 1974, ApJS, 27, 37

Giesers, B., Dreizler, S., Husser, T.-O., et al. 2018, MNRAS, 475, L15

Giesers, B., Kamann, S., Dreizler, S., et al. 2019, A&A, 632, A3

Grindlay, J. E., Heinke, C., Edmonds, P. D., & Murray, S. S. 2001a, Sci, 292, 2290

Grindlay, J. E., Heinke, C. O., Edmonds, P. D., Murray, S. S., & Cool, A. M. 2001b, ApJL, 563, L53

Grindlay, J. E., Hertz, P., Steiner, J. E., Murray, S. S., & Lightman, A. P. 1984, ApJL, 282, L13

Harris, W. E. 1996, AJ, 112, 1487

Heggie, D. C., & Giersz, M. 2014, MNRAS, 439, 2459

Heinke, C. O., Grindlay, J. E., Edmonds, P. D., et al. 2003, ApJ, 598, L16

Heinke, C. O., Grindlay, J. E., Edmonds, P. D., et al. 2005, ApJ, 625, 796

Heinke, C. O., Ivanov, M. G., Koch, E. W., et al. 2020, MNRAS, 492, 5684

Heinke, C. O., Wijnands, R., Cohn, H. N., et al. 2006, ApJ, 651, 1098

Hertz, P., & Grindlay, J. E. 1983, ApJL, 267, L83

Hills, J. G. 1976, MNRAS, 175, 1P

Homer, L., Charles, P. A., Naylor, T., et al. 1996, MNRAS, 282, L37

Huang, R. H. H., Becker, W., Edmonds, P. D., et al. 2010, A&A, 513, A16

Hui, C. Y., Cheng, K. S., & Taam, R. E. 2010, ApJ, 714, 1149

Ivanova, N., Heinke, C. O., Rasio, F. A., et al. 2006, MNRAS, 372, 1043

Johnson, T. J., Guillemot, L., Kerr, M., et al. 2013, ApJ, 778, 106

Johnston, H. M., Kulkarni, S. R., & Phinney, E. S. 1992, in X-Ray Binaries and the Formation of Binary and Millisecond Radio Pulsars, ed. E. P. J. van den Heuvel (Dordrecht: Kluwer), 349

Katz, J. I. 1975, Natur, 253, 698

Kong, A. K. H., Bassa, C., Pooley, D., et al. 2006, ApJ, 647, 1065

Kremer, K., Chatterjee, S., Rodriguez, C. L., & Rasio, F. A. 2018, ApJ, 852, 29

Kulkarni, S. R., Anderson, S. B., Prince, T. A., & Wolszczan, A. 1991, Natur, 349, 47

Kulkarni, S. R., Hut, P., & McMillan, S. 1993, Natur, 364, 421

Lewin, W. H. G., & Joss, P. C. 1983, in Accretion-Driven Stellar X-ray Sources, ed. W. H. G. Lewin & E. P. J. van den Heuvel (Cambridge: Cambridge Univ. Press), 41

Lu, T.-N., Kong, A. K. H., Bassa, C., et al. 2009, ApJ, 705, 175

Lynch, R. S., Freire, P. C. C., Ransom, S. M., & Jacoby, B. A. 2012, ApJ, 745, 109

Lyne, A. G., Brinklow, A., Middleditch, J., Kulkarni, S. R., & Backer, D. C. 1987, Natur, 328, 399

Lyne, A. G., Manchester, R. N., & D'Amico, N. 1996, ApJL, 460, L41

Maccarone, T. J. 2005, MNRAS, 360, L30

Maccarone, T. J., Kundu, A., Zepf, S. E., & Rhode, K. L. 2007, Natur, 445, 183

Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005, AJ, 129, 1993

Migliari, S., & Fender, R. P. 2006, MNRAS, 366, 79

Miller-Jones, J. C. A., Strader, J., Heinke, C. O., et al. 2015, MNRAS, 453, 3918

Mooley, K. P., Miller-Jones, J. C. A., Fender, R. P., et al. 2017, MNRAS, 467, L31

Moreno, E., Pichardo, B., & Velázquez, H. 2014, ApJ, 793, 110

Morscher, M., Pattabiraman, B., Rodriguez, C., Rasio, F. A., & Umbreit, S. 2015, ApJ, 800, 9

Morscher, M., Umbreit, S., Farr, W. M., & Rasio, F. A. 2013, ApJL, 763, L15

Osten, R. A., Brown, A., Ayres, T. R., et al. 2000, ApJ, 544, 953

Peuten, M., Zocchi, A., Gieles, M., Gualandris, A., & Hénault-Brunet, V. 2016, MNRAS, 462, 2333

Pooley, D., & Hut, P. 2006, ApJL, 646, L143

Pooley, D., Lewin, W. H. G., Anderson, S. F., et al. 2003, ApJL, 591, L131

Pooley, D., Lewin, W. H. G., Homer, L., et al. 2002a, ApJ, 569, 405

Pooley, D., Lewin, W. H. G., Verbunt, F., et al. 2002b, ApJ, 573, 184

Rivera Sandoval, L. E., van den Berg, M., Heinke, C. O., et al. 2018, MNRAS, 475, 4841

Rutledge, R. E., Bildsten, L., Brown, E. F., Pavlov, G. G., & Zavlin, V. E. 2002, ApJ, 578, 405

Rutledge, R. E., Fox, D. W., Kulkarni, S. R., et al. 2004, ApJ, 613, 522

Saito, Y., Kawai, N., Kamae, T., et al. 1997, ApJL, 477, L37

Scholz, F. W., & Stephens, M. A. 1987, JASA, 82, 918

Shishkovsky, L., Strader, J., Chomiuk, L., et al. 2018, ApJ, 855, 55

Shishkovsky, L., Strader, J., Chomiuk, L., et al. 2020, ApJ, 903, 73

Sidoli, L., Paizis, A., Bazzano, A., & Mereghetti, S. 2006, A&A, 460, 229

Sigurdsson, S., & Hernquist, L. 1993, Natur, 364, 423

Sippel, A. C., & Hurley, J. R. 2013, MNRAS, 430, L30

Strader, J., Chomiuk, L., Maccarone, T. J., Miller-Jones, J. C. A., & Seth, A. C. 2012, Natur, 490, 71

Sutantyo, W. 1975, A&A, 44, 227

Swank, J. H., Becker, R. H., Pravdo, S. H., Saba, J. R., & Serlemitsos, P. J. 1976, IAU Circ., 3010, 1

Tetarenko, B. E., Bahramian, A., Arnason, R. M., et al. 2016, ApJ, 825, 10

Thorsett, S. E., Arzoumanian, Z., McKinnon, M. M., & Taylor, J. H. 1993, ApJL, 405, L29

Thorsett, S. E., & Chakrabarty, D. 1999, ApJ, 512, 288

Tremou, E., Strader, J., Chomiuk, L., et al. 2018, ApJ, 862, 16

Urquhart, R., Bahramian, A., Strader, J., et al. 2020, ApJ, 904, 57

Verbunt, F. 2003, in ASP Conf. Ser. 296, New Horizons in Globular Cluster Astronomy, ed. G. Piotto et al. (San Francisco, CA: ASP), 245

Verbunt, F., & Hut, P. 1987, in IAU Symp. Vol. 125, The Origin and Evolution of Neutron Stars, ed. D. J. Helfand & J. H. Huang (Cambridge: Cambridge Univ. Press), 187

Verbunt, F., Pooley, D., & Bassa, C. 2008, in IAU Symp. 246, Dynamical Evolution of Dense Stellar Systems, ed. E. Vesperini, M. Giersz, & A. Sills (Cambridge: Cambridge Univ. Press), 301

Wang, L., Peng, B., Stappers, B. W., et al. 2020, ApJ, 892, 43

Weatherford, N. C., Chatterjee, S., Kremer, K., & Rasio, F. A. 2020, ApJ, 898, 162

Weatherford, N. C., Chatterjee, S., Rodriguez, C. L., & Rasio, F. A. 2018, ApJ, 864, 13

Zhao, Y., Heinke, C. O., Cohn, H. N., Lugger, P. M., & Cool, A. M. 2019, MNRAS, 483, 4560

Zhao, Y., Heinke, C. O., Tudor, V., et al. 2020, MNRAS, 493, 6033