How black holes shape globular clusters: Modeling NGC 3201
Kyle Kremer Claire S. Ye Sourav Chatterjee Carl L. Rodriguez Frederic A. Rasio

Modeling NGC 3201
Kremer et al.




0000-0002-4086-3180]Kyle Kremer
 Department of Physics & Astronomy, Northwestern University, Evanston, IL 60202, USA
 Center for Interdisciplinary Exploration & Research in Astrophysics (CIERA), Evanston, IL 60202, USA

 Department of Physics & Astronomy, Northwestern University, Evanston, IL 60202, USA
 Center for Interdisciplinary Exploration & Research in Astrophysics (CIERA), Evanston, IL 60202, USA
0000-0002-3680-2684]Sourav Chatterjee
 Department of Physics & Astronomy, Northwestern University, Evanston, IL 60202, USA
 Center for Interdisciplinary Exploration & Research in Astrophysics (CIERA), Evanston, IL 60202, USA
Tata Institute of Fundamental Research, Homi Bhabha Road, Mumbai 400005, India

MIT-Kavli Institute for Astrophysics and Space Research, Cambridge, MA 02139, USA
0000-0002-7132-418X]Frederic A. Rasio
 Department of Physics & Astronomy, Northwestern University, Evanston, IL 60202, USA
 Center for Interdisciplinary Exploration & Research in Astrophysics (CIERA), Evanston, IL 60202, USA


Numerical simulations have shown that black holes (BHs) can strongly influence the evolution and present-day observational properties of globular clusters (GCs). Using a Monte Carlo code, we construct GC models that match the Milky Way (MW) cluster NGC 3201, the first cluster in which a stellar-mass BH was identified through radial-velocity measurements. We predict that NGC 3201 contains ≳ 200 stellar-mass BHs. Furthermore, we explore the dynamical formation of main sequence–BH binaries and demonstrate that systems similar to the observed BH binary in NGC 3201 are produced naturally. Additionally, our models predict the existence of bright blue-straggler–BH binaries unique to core-collapsed clusters, which otherwise retain few BHs.







§ INTRODUCTION



In recent years, an increasing number of stellar-mass black hole (BH) candidates have been observed within globular clusters (GCs) as members of binary systems with luminous companions. The first such systems were observed in mass-transferring configurations and identified through radio and/or X-ray observations. To date, mass-transferring binaries with BH-candidate accretors have been identified in four Milky Way (MW) GCs: M10 <cit.>, M22 <cit.>, 47 Tuc <cit.>, and M62 <cit.>.  Additionally, multiple candidates have been identified in extragalactic GCs <cit.>.

Most recently, a stellar-mass BH candidate was identified in the MW GC NGC 3201 through radial-velocity observations in a detached binary with a main-sequence (MS) companion <cit.>. This observation marks the first detection of a stellar-mass BH candidate through radial-velocity methods.

The effect of BHs on the evolution of GCs is of high current interest, especially since the recent discoveries of merging binary BHs (BBHs) by LIGO <cit.>. Several recent studies have demonstrated that dynamical processing in GCs can be a dominant formation channel for the BBHs observed by LIGO <cit.>. Furthermore, numerical simulations have shown that the retention fraction of BHs in GCs can significantly influence the way the host cluster evolves <cit.>, which in turn affects the overall BH dynamics and BBH formation. 

The observation of stellar-mass BH candidates in GCs has sparked extensive theoretical study of the dynamical formation of BH–non-BH binaries in GCs. Traditionally, it has been argued that finding even a few BH–binary candidates in GCs indicates a much larger population of undetectable BHs retained in those clusters <cit.>. However, recent analyses <cit.> have shown that while BH–non-BH binaries can readily form in GCs in both detached and mass-transferring configurations, the number of BH–non-BH binaries retained within a cluster at late times is independent of the total number of BHs retained.

In this letter, we use our Cluster Monte Carlo code () to construct a GC model that matches the observed properties of NGC 3201, with the goal of exploring the formation of detached BH–non-BH binaries in such a cluster. We demonstrate a large population of BHs (≳ 200) is necessary to produce a GC similar to NGC 3201. Furthermore, we show that our GC models that are most similar to NGC 3201 contain BH–MS binaries with orbital parameters similar to those of the system recently observed in NGC 3201.

In Section <ref>, we describe our method for modeling GCs and introduce the grid of models used for this analysis. In Section <ref>, we display various observational properties of our models and identify the models which best resemble the observed properties of NGC 3201. In Section <ref>, we discuss the properties and dynamical history of BH–MS binaries found in our models at late times and compare these systems to the BH–MS recently identified in NGC 3201. In Section <ref>, we discuss the formation of blue straggler–BH binaries found in our core-collapsed GC models which retain few BHs. We discuss our results and conclude in Section <ref>.



§ METHOD





 §.§ Globular cluster models


In order to model the long-term evolution of GCs, we use  <cit.>.  incorporates all physics relevant to the evolution of BHs and BH binaries in GCs, including two-body relaxation <cit.>, single and binary star evolution <cit.>, three-body binary formation <cit.>, galactic tides <cit.>, and small-N gravitational encounters calculated using  <cit.>, updated to incorporate post-Newtonian effects in all three- and four-body encounters <cit.>.



 
 Here, we present a grid of 16 GC models aimed at representing a cluster like NGC 3201. Fixed initial cluster parameters include total particle number (N=8 × 10^5), virial radius (r_v = 1 pc), King concentration parameter (W_0=5), binary fraction (f_b = 5%), and metallicity (Z=0.001).


The initial stellar masses are sampled from the initial mass function (IMF) given in <cit.> in the range 0.08–150 M_⊙. 

The initial number of binary stars is chosen based on the specified f_b and N. Primary masses for all binaries are sampled from the same IMF as single stars, and secondary star masses are assigned assuming a flat distribution in mass ratios. As in previous papers <cit.>, initial orbital periods, P, for all binaries are drawn from a distribution of the form dn/dlog P∝ P^0 and the initial eccentricities are thermal.
 
To treat stellar remnant formation, we use a modified prescription from that implemented in  and  by using the results of <cit.> and <cit.>. Upon formation, all neutron star (NS) remnants receive natal kicks drawn from a Maxwellian distribution with σ_NS = 265 km s^-1 <cit.>. For BH natal kicks, we assume kick magnitudes independent of the BH masses and drawn from a Maxwellian with dispersion width, σ_BH. The ratio σ_BH/σ_NS is varied between models from σ_BH/σ_NS = 0.005 up to σ_BH/σ_NS = 1.0 (full NS kicks), as listed in column 2 of Table <ref>. Effectively, we use BH natal kicks <cit.> as a simple and convenient way to vary the number of BHs retained in our GC models. BH natal kicks are intended to serve as a proxy for more physical mechanisms that determine the numbers of BHs retained in GCs over long timescales.
 


c|c||c|cc|c|c

0pt
Cluster properties for all model GCs

	Model    
    σ_BH/σ_NS    
    M_tot   
    r_c   
    r_h   
    N_BH   
    N_BH-MS

        
        
    (10^5 M_⊙)   
    2c(pc)   
       
    10 Gyr <t< 12 Gyr

1     0.005     2.11     1.27     4.25     404     2 

2     0.01     2.1     1.97     4.46     424     2 

3     0.02     2.11     3.04     4.91     397     1 

4     0.03     2.09     3.77     4.91     389     0 

5     0.04     2.14     1.83     4.84     379     10 

6     0.05     2.12     3.13     4.12     383     0 

7     0.06     2.2     1.88     3.55     319     4 

8     0.07     2.16     1.64     4.43     285     1 

9     0.08     2.21     1.63     3.26     286     5 

10     0.09     2.22     2.23     2.8     249     2 

11     0.1     2.24     0.9     2.58     192     1 

12     0.2     2.31     0.11     1.64     7     9 

13     0.4     2.22     0.22     1.82     2     3 

14     0.6     2.33     0.21     1.78     4     0 

15     0.8     2.36     0.32     1.71     2     1 

16     1.0     2.3     0.21     1.56     3     2 

Column 2 shows the scaling of BH natal kicks, σ_BH/σ_NS, used for each model. Columns 3-6 show properties of each model at t=12 Gyr. Note that all models form ∼ 1500 BHs initially. Column 6 shows the number of BHs that are retained at t=12 Gyr. Column 7 shows the number of different detached BH–MS binaries that appear in snapshots in the range 10 Gyr <t< 12 Gyr, spanning the approximate uncertainty in the age of NGC 3201.














 §.§ `Observing' model clusters



It is not straightforward to compare GC models to observational data for a specific GC such as NGC 3201.

Each of our GC models contains ∼100-1000 snapshots in time spaced ∼10-100 Myr apart 
spanning the full evolution of the cluster, from formation to present day. From these snapshots, we create two-dimensional projections for each model at the corresponding snapshot times, assuming spherical symmetry. The half-light radius, r_hl, is estimated by finding the projected radius which contains half of the cluster's total light. To estimate the core radius, r_c, we use the method described in <cit.> and <cit.>: we fit an analytic King model approximation <cit.> to the cumulative stellar luminosity at a given projected radius including stars within a projected distance of r_hl.

The methods used to calculate the surface brightness profile and velocity dispersion profile are described in detail in the Appendix.



§ RESULTS





Figure <ref> shows the evolution of r_c and r_hl for the 16 models in Table <ref>. Black lines show evolution of models 1-11, all of which retain ≳ 200 BHs at late times. Blue lines show models 12-16 which contain ≲ 10 BHs at late times (see Table <ref>). Horizontal orange bars show r_c and r_hl of MW GCs with absolute V-band magnitude and heliocentric distance similar to NGC 3201 (in the ranges of -8.95 to -5.95 magnitude and 6-12 kpc, respectively; data taken from <cit.>). As Figure <ref> shows, these MW GCs span a wide range in r_c and r_hl. As a limiting case, this set of MW GCs includes the core-collapsed cluster NGC 6397 (the green bar), which is shown to align with models containing very few BHs at present.


The red star marks the observed r_c and r_hl for NGC 3201, which are consistent with models that retain large numbers of BHs.



The sustained heating and expansion shown by the red curves in Figure <ref> are primarily due to the frequent super-elastic interactions in the centrally-concentrated BH system: energy is deposited onto the stellar background when a recoiled BH sinks back via dynamical friction
<cit.>.

Figure <ref> shows the surface brightness profile for each of the 16 models at t=12 Gyr (solid black and blue lines) compared to the observed surface brightness profile of NGC 3201 <cit.>. As in Figure <ref>, black lines show models 1-11 of Table <ref> and blue lines show models 12-16. 

As Figure <ref> shows, models that retain many BHs have surface brightness profiles closest to matching NGC 3201. The surface brightness profiles of models retaining few BHs exhibit cusps at low r, representative of so-called “core-collapsed” MW GCs.

Henceforth, we use the term “core-collapsed” to simply denote those models (12-16) that remain much more centrally concentrated than their BH-heated counterparts (models 1-11), throughout their evolution, and have surface brightness profiles with prominent central cusps at late times as shown in Figure <ref>.

The surface brightness profiles of all models converge with NGC 3201 for r ≳ 100 arcsec. This is simply because all models have similar total mass at late times (see column 3 of Table <ref>), which is a consequence of the choice of initial particle number (N=8 × 10^5 for all models), and the same galactocentric distance (and hence tidal radius).



In order to identify a single best-fit model which accurately matches NGC 3201, we calculate the surface brightness profile, velocity dispersion profile, r_c, and r_hl for all snapshots at late times for each of the 16 models and compare to the observational data of NGC 3201. Note that several different estimates for the current age of NGC 3201 exist, ranging from ∼10 Gyr <cit.> up to ∼12 Gyr <cit.>. To reflect this uncertainty in age, we simply consider all cluster snapshots in the range  10 Gyr ≤ t ≤ 12 Gyr (the maximum evolution time considered in our grid of models).

We identify model 5 at t=11.2 Gyr as our best-fit model. For this model, r_c=1.85 pc and r_hl=4.42 pc, which are in exact agreement with the r_c and r_hl values for NGC 3201 given in <cit.>. This particular model retains 389 BHs, including 4 detached BH–MS binaries. The assumptions that go into choosing our best-fit model are discussed in greater detail in the Appendix. 




 §.§ BH–MS binaries



As shown in <cit.> and <cit.>, the formation of BH–non-BH binaries in a GC is a self-regulated process limited by a complex competition between the number of BHs and the number density of non-BH–non-BH binaries in the region of the GC where BHs mix with non-BHs. This competition ensures that the number of BH–non-BH binaries found in a GC at late times remains independent of the total number of BHs retained in the cluster.

All GC models considered here produce at most a few BH–MS binaries at late times (column 7 of Table <ref>), regardless of the total number of BHs retained at late times (column 6), in agreement with the results of <cit.> and <cit.>

In order to explore the formation of BH–MS binaries similar to the binary system recently detected in NGC 3201 <cit.>, we search for detached BH–MS binaries in models 1-11 (models that retain large numbers of BHs at late times and have structural parameters similar to those of NGC 3201; see Figure <ref>).

In total, we find 24 different detached BH–MS binaries in models 1-11 in cluster snapshots between 10 Gyr < t < 12 Gyr. The top panel of Figure <ref> shows the MS companion mass, M_MS, vs semi-major axis, a, for each of these binaries. Filled circles in Figure <ref> mark the 4 BH–MS binaries found in model 5 at t=11.2 Gyr (our best-fit model described in the previous section). Open circles show BH–MS binaries found in all remaining snapshots in the range 10 Gyr < t < 12 Gyr for models 1-11. The horizontal dashed line marks the turn-off mass (M_to=0.83 M_⊙) for clusters of this particular age and metallicity. The black “x” marks the quoted values of the MS mass (0.81± 0.05 M_⊙) and minimum a (1.03 ± 0.03 au) for the BH–MS binary detected in NGC 3201 <cit.>.

The bottom panel of Figure <ref> shows the eccentricity, e, vs a for all BH–MS found in models 1-11. Clearly, the majority of the detached BH–MS binaries identified at late times have moderate to high eccentricities, a consequence of the dynamical encounters experienced by these binaries.

Indeed, dynamical encounters play a critical role in the formation and evolution of these binaries. Of the 24 BH–MS binaries identified in these models, all are formed dynamically (none are primordial binaries). The median number of encounters experienced by the BHs in these binaries is 9. The median number of dynamical encounters experienced by the BH–MS binaries after formation, via an exchange encounter, is 2. Note that these BH–MS binaries are susceptible to destruction as the result of exchange encounters with other cluster objects. The median lifetime for the BH–MS binaries shown in Figure <ref> is ∼ 400 Myr.

The median BH mass for these 24 binaries is 9 M_⊙, consistent with the minimum mass estimated for the NGC 3201 BH-candidate  <cit.>. A low BH-mass is anticipated based on our understanding of the evolution of BHs in GCs. The most massive BHs retained after formation will be ejected at early times through dynamical encounters. Only the least massive BHs will remain at late times <cit.>. Thus the detected BH–MS binary has properties that would be commonly produced dynamically in a GC similar to NGC 3201. 

Note that, in addition to the 24 detached BH–MS binaries, we also find 14 BH–MS binaries in mass-transferring configurations at late times, similar to those studied in <cit.>.





 §.§ Blue straggler–BH binaries



Although models 12-16 are core-collapsed at late times and are therefore poor representations of NGC 3201, we show orbital parameters of all detached  BH–MS binaries found in these models in Figure <ref>. Interestingly, for these models, 


a large fraction of the MS star masses lie above the turn-off mass. These systems are blue stragglers (BSs) formed through collisions earlier in the evolution of their host cluster. Three of these BS companions have mass M_MS≳ 2   M_⊙, which is more than twice the turn-off mass. These are created via multiple collisions.




Among our models, the BH–BS binaries are unique to the core-collapsed clusters retaining few BHs. This is consistent with our understanding of GC evolution. As the number of retained BHs decreases, the central density increases, ultimately leading to core collapse. Additionally, in the absence of large numbers of retained BHs, the high-mass MS stars and MS binaries increasingly dominate the cluster's core, leading to an increased rate of dynamical interactions between multiple MS stars. Most of these interactions are resonant in nature and can often result in collisions, efficiently creating BSs <cit.>. Being more massive than the average cluster object, the BSs readily interact and form binaries with the remaining BHs, like the systems seen in Figure <ref>. Such BH–BS binaries could be detectable through radial-velocity methods similar to those used to identify the BH candidate in NGC 3201. We do note that BSs have been observed in non-core-collapsed MW GCs, including NGC 3201 <cit.>. We do find BSs in our BH-rich models, however none of these BSs are found with BH binary companions. We reserve a more detailed examination of BS formation in GCs for a future study.

Comparison of Figures <ref> and <ref> show that the BH–MS binaries found in NGC 3201-like GCs and those found in core-collapsed clusters occupy significantly different regions of the parameter space. In addition to the presence of BH–BS binaries in core-collapsed models, the BH–MS binaries found in NGC 3201-like GCs tend to have higher eccentricities and are found over a wider range in a (no BH–MS binaries with a ≳ 10 au are found in core-collapsed clusters). This is a direct consequence of higher central velocity dispersions in the more compact core-collapsed clusters which allows only very hard binaries to survive.




§ CONCLUSION AND DISCUSSION



Using our Cluster Monte Carlo code, , we have demonstrated a method to model the MW GC NGC 3201. We showed that by varying the magnitude of BH natal kicks,  the number of retained BHs at late times changes substantially, which in turn, determines the observational features of our GC models at late times. In particular, we showed that in order to produce a GC model that matches the observational properties (including surface brightness profile, velocity dispersion profile, core radius, and half-light radius) of NGC 3201 at present, the model must retain ≳ 200 BHs.

Additionally, we demonstrated that GC models retaining large numbers of BHs at late times can harbor detached BH–MS binaries similar in properties to the system recently detected through radial-velocity measurements in NGC 3201.

We also showed that although GCs which retain few BHs are poor representations of NGC 3201, they may serve as good representations of core-collapsed MW GCs. 

We demonstrated that such core-collapsed GCs may contain bright blue straggler–BH binaries, that could be detectable using radial-velocity methods similar to those used to identify the BH–MS system in NGC 3201.

Previous analyses have modeled other GCs with known stellar-mass BH candidates. For example, <cit.> (Monte Carlo) and <cit.> (N-body) explored models of M22, and demonstrated similar results to those of this letter; namely, that retained BHs provide an energy source capable of producing models that match GCs with large observed core radii at late times and that a fraction of the BHs remaining at late times are likely to be found in binaries with luminous companions.

Additionally, previous analyses have shown that, in the absence of BHs, binary hardening may provide sufficient energy to roughly balance the energy diffusion from two-body relaxation, after contraction of the cluster's core. This allows the cluster to attain a quasi-steady state even when the binary fraction is only a few percent <cit.>.

<cit.> showed how clusters that have entered this “binary-burning” phase exhibit cuspy surface brightness profiles typical of the observed core-collapsed MW GCs and similar to GC models presented here that retain few BHs. The binary fraction for NGC 3201 is estimated to be low <cit.>, consistent with our choice of initial binary fraction. This suggests that large numbers of BHs must be retained to produce a cluster similar to NGC 3201 since models with too few retained BHs would lead to surface brightness profiles vastly different from that of NGC 3201. See the Appendix for further discussion of the possible effects of binaries.

In this study we used BH natal kicks as a proxy for adjusting the number of BHs retained in GCs at late times. In reality, other factors likely determine whether or not a particular GC retains large numbers of BHs at present. In particular, initial cluster parameters such as virial radius, cluster mass, binary fraction, concentration, and galactocentric distance determine how dynamically evolved a GC is at present. This in turn may determine how efficiently BHs are dynamically ejected from GCs via recoil from dynamical encounters. We will explore effects of these other initial parameters in future work. 




While we consider only BH–MS binaries in this paper, populations of other classes of compact objects may also correlate with GC evolution, and thus, the retained  BH population. In particular, the number of neutron stars (NSs) in GCs that can be observed as X-ray binaries (NS–XRBs) or as millisecond pulsars may correlate with the number of BHs retained at late times. We note that no pulsars or NS–XRBs have been observed in NGC 3201,  unlike, for example, the core-collapsed cluster NGC 6397. We also find that NS-XRBs and millisecond pulsars rarely form in our models retaining many BHs (NGC 3201-like models), but instead we find these objects form readily in our core-collapsed models retaining few BHs.









We thank the referee for helpful comments and suggestions.

This work was supported by NASA ATP Grant NNX14AP92G 
and NSF Grant AST-1716762. K.K. acknowledges support by the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE-1324585.
S.C. acknowledges support from
CIERA, the National Aeronautics and Space Administration
through a Chandra Award Number TM5-16004X/NAS8-
03060 issued by the Chandra X-ray Observatory Center
(operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics
Space Administration under contract NAS8-03060), 
and Hubble Space Telescope Archival research 
grant HST-AR-14555.001-A (from the Space Telescope 
Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555).


[Abbott et al.(2016a)]Abbott2016a Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016a, ApJL, 818, L22
[Abbott et al.(2016b)]Abbott2016b Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016b, PhRvL, 116, 061102
[Abbott et al.(2016c)]Abbott2016c Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016c, PhRvL, 116, 241103
[Abbott et al.(2016d)]Abbott2016d Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016d, PhRvL, 116, 061102
[Abbott et al.(2016e)]Abbott2016e Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016e, PhRvL 118 221101
[Abbott et al.(2017)]Abbott2017  Abbott, B. P., Abbott, R., Abbott, T. D., et al., PhRvL 118, 221101
[Alessandrini et al.(2016)] Alessandrini2016 Alessandrini, E., Lanzoni, B., Ferraro, F. R., Miocchi, P., &
Vesperini, E. 2016, ApJ, 833, 252
[Amaro-Seoane & Chen et al.(2016)]AmaroSeoane2016 Amaro-Seoane, P. & Chen, X. 2016, MNRAS, 458, 3075
[Antognini et al.(2014)] Antognini2014 Antognini, J. M., Shappee, B. J., Thompson, T. A., & Amaro-Seoane, P. 2014, MNRAS, 439, 1079
[Banerjee et al.(2010)] Banerjee2010 Banerjee, S. and Baumgardt, H. & Kroupa, P. 2010, MNRAS, 402, 371
[Belczynski et al. (2002)] Belczynski2002 Belczynski, K., Kalogera, V., & Bulik, T. 2002, ApJ, 572, 407
[Chatterjee et al.(2010)] Chatterjee2010 Chatterjee, S., Fregeau, J. M., Umbreit, S., & Rasio, F. A. 2010, ApJ, 719, 915
[Chatterjee et al.(2013a)] Chatterjee2013a Chatterjee, S., Umbreit, S., Fregeau, J. M., & Rasio, F. A. 2013, MNRAS, 429, 2881
[Chatterjee et al.(2013b)] Chatterjee2013b Chatterjee, S., Rasio, F. A., Sills, A., & Glebbeek, E. 2013, ApJ, 777, 106
[Chatterjee et al.(2017)] Chatterjee2017 Chatterjee, S., Rodriguez, C. L., & Rasio, F. A. 2017, ApJ, 834, 68
[Chomiuk et al.(2013)] Chomiuk2013 Chomiuk, L., Strader, J., Maccarone, T. J., Miller-Jones, J. C. A., Heinke, C., Noyola, E., Seth, A. C., & Ransom, S. 2013, ApJ, 777, 69
[Côté et al.(1994)] Cote1994 Côté, P., Welch, D. L., Fischer, P., Da Costa, G. S., Tamblyn, P et al. 1994, ApJSS, 90, 83
[Côté et al.(1995)] Cote1995 Côté, P., Welch, D. L., Fischer, P., & Gebhardt, K. 1995, ApJ, 454, 788
[Covino & Ortolani(1997)] Covino1997 Covino, S. & Ortolani, S. 1997, A&A, 318, 40
[Forbes & Bridges(2010)]Forbes2010 Forbes, D. A. & Bridges, T. 2010, MNRAS, 404, 3
[Fregeau et al.(2003)] Fregeau2003 Fregeau, J. M., Gu ̈rka, M. A., Joshi, K. J., & Rasio, F. A. 2003, ApJ, 593, 772
[Fregeau et al.(2004)] Fregeau2004 Fregeau, J.M., Cheung, P., Portegies Zwart, S.F., & Rasio, F. A. 2004,MNRAS, 352, 1 
[Fregeau & Rasio(2007)] Fregeau2007 Fregeau, J. M., & Rasio, F. A. 2007, ApJ, 658, 1047
[Fryer & Kalogera (2001)] Fryer2001 Fryer, C. L. & Kalogera, V. 2001, ApJ, 554, 548
[Giesers et al.(2018)] Giesers2018 Giesers, B., Dreizler, S., Husser, T.-O., Kamann, S., Escudé, G. A., et al. 2018, MNRAS, 475, L15
[Harris1996(2010 edition)] Harris1996 Harris, W.E. 1996, AJ, 112, 1487
[Heggie & Hut(2003)] Heggie2003 Heggie D. & Hut P. 2003, The Gravitational Million-Body Problem: A Multidisciplinary Approach to Star Cluster Dynamics. Cambridge University Press, Cambridge
[Heggie & Giersz(2014)] Heggie2014 Heggie, D. C. & Giersz, M. 2014, MNRAS, 439, 2459
[Hénon(1971a)] Henon1971a Hénon, M. 1971, AP&SS, 14, 151
[Hénon(1971b)] Henon1971b Hénon, M. 1971, AP&SS, 13, 284
[Hobbs et al.(2005)] Hobbs2005 Hobbs, G., Lorimer, D. R., Lyne, A. G., Kramer, M. 2005, MRNAS, 360, 974
[Hurley et al.(2000)] Hurley2000 Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543
[Hurley et al.(2002)] Hurley2002 Hurley, J. R., Tout, C. A. & Pols, O. R. 2002, MNRAS, 329, 897
[Irwin et al.(2010)] Irwin2010 Irwin, J. A., Brink, T. G., Bregman, J. N., & Roberts, T. P. 2010, ApJL, 712, L1
[Joshi et al.(2000)] Joshi2000 Joshi, K. J., Rasio, F. A., & Portegies Zwart, S. 2000, ApJ, 540, 969
[Joshi et al.(2001)] Joshi2001 Joshi, K. J., Nave, C. P., & Rasio, F. A. 2001, ApJ, 550, 691
[Kamann et al.(2018)] Kamann2017 Kamann, S., Husser, T.-O., Dreizler, S., Emsellem, E., Weilbacher, P. M. et al. 2018, MNRAS, 473, 4
[King(1962)] King1962 King, I. 1962, AJ, 67, 471
[Kremer et al.(2018)] Kremer2018 Kremer, K., Chatterjee, S., Rodriguez, C. L., & Rasio, F. A. 2018, ApJ, 852, 29
[Kroupa(2001)] Kroupa2001 Kroupa, P. 2001, MNRAS, 322, 23
[Maccarone et al.(2007)] Maccarone2007 Maccarone, T. J., Kundu, A., Zepf, S. E., & Rhode, K. L. 2007, Nature, 445, 183
[Mackey et al(2007)]Mackey2007 Mackey, A. D., Wilkinson, M. I., Davies, M. B., & Gilmore, G. F. 2007, MNRAS, 379, L40
[Mackey et al.(2008)] Mackey2008 Mackey, A. D., Wilkinson, M. I., Davies, M. B., & Gilmore, G. F. 2008, MNRAS, 386, 65
[Mandel(2016)] Mandel2016 Mandel, I. 2016, MNRAS, 456, 578
[Miller-Jones et al.(2014)] Miller-Jones2014 Miller-Jones, J., Maccarone, T., Chomiuk, L., Strader, J., Bogdanov, S., Sivakoff, G., & Heinke, C. 2014a, A new black hole candidate in the globular cluster 47 Tucanae, ATNF Proposal
[Morscher et al.(2013)] Morscher2013 Morscher, M., Pattabiraman, B., Rodriguez, C., Rasio, F. A., & Umbreit, S. 2015, ApJ, 800, 9
[Morscher et al.(2015)] Morscher2015 Morscher, M., Umbreit, S., Farr, W. M., & Rasio, F. A. 2013, ApJL, 763, L15
[Pattabiraman et al.(2013)] Pattabiraman2013 Pattabiraman, B., Umbreit, S., Liao, W.-k., et al. 2013, ApJS, 204, 15
[Pryor & Meylan(1993)] Pryor1993 Pryor C. & Meylan G., 1993, in Djorgovski S. G., Meylan G., eds, Astronomical Society of the Pacific Conference Series Vol. 50, Structure and Dynamics of Globular Clusters, p. 357
[Repetto et al.(2012)] Repetto2012 Repetto, S., Davies, M. B., & Sigurdsson, S. 2012, MNRAS, 425, 2799
[Rodriguez et al.(2015)] Rodriguez2015 Rodriguez, C. L., Morscher, M., Pattabiraman, B., et al. 2015, PhRvL, 115, 051101
[Rodriguez et al.(2016)] Rodriguez2016 Rodriguez, C. L., Chatterjee, S., & Rasio, F. A. 2016, PhRvD, 93, 84029
[Rodriguez et al(2018)] Rodriguez2018 Rodriguez, C. L., Amaro-Seoane, P., Chatterjee, S., & Rasio, F. A. 2018, preprint (arXiv: 1712.04937)
[Sills et al.(2013)] Sills2013 Sills, A., Glebbeek, E., Chatterjee, S., & Rasio, F. A. 2013, ApJ, 777, 105
[Simunovic & Puzia(2014)] Simunovic2014 Simunovic, M. & Puzia, T. H. 2014, ApJ, 782, 49
[Sippel & Hurley(2012)] Sippel2012 Sippel, A. C. & Hurley, J. R. 2012, MNRAS, 430, L30
[Strader et al.(2012)] Strader2012 Strader, J., Chomiuk, L., Maccarone, T. J., Miller-Jones, J. C. A., & Seth, A. C. 2012, Nature, 490, 71
[Strader(2014)] Strader2014 Strader, J. 2014, A black hole in the Galactic Globular Cluster, Chandra Proposal, M10
[Trager et al.(1995)]Trager1995 Trager, S.C., King, I.R., and Djorgovski, S. 1995, AJ 109, 218.
[Umbreit et al.(2012)] Umbreit2012 Umbreit, S., Fregeau, J. M., Chatterjee, S., & Rasio, F. A. 2012, ApJ, 750, 31
[Usher et al.(2017)]Usher2017 Usher, C., Pastorello, N., Bellstedt, S., Alabi, A., Cerulo, P. et al. 2017, MNRAS, 468, 4
[Vesperini & Chernoff(1994)] Vesperini1994 Vesperini E. & Chernoff D. F. 1994, ApJ, 431, 231
[Ziosi et al.(2014)] Ziosi2014 Ziosi, B. M., Mapelli, M., Branchesi, M., & Tormen, G. 2014, MNRAS, 441, 3703
[Zocchi et al.(2012)] Zocchi2012 Zocchi, A., Bertin, G., & Varri, A. L. 2012, A&A, 539, A65






§ CALCULATING SURFACE BRIGHTNESS AND VELOCITY DISPERSION PROFILES
 




Figure <ref> shows the surface brightness profile (top panel) and velocity dispersion profile (bottom two panels) of our best-fit model: model 5 in Table <ref> at t=11.2 Gyr.  Here we describe our method for calculating surface brightness profiles and velocity dispersion profiles for our models at different points in time.

To calculate the surface brightness profile of a model GC, we produce two-dimensional projections of each time snapshot as described in Section <ref>, then divide these two-dimensional projections into 50 equally spaced radial bins. To calculate the surface brightness within each bin, we remove all stars with luminosity L_⋆ > 15   L_⊙.  We compare to the observed surface brightness profile from <cit.>. Note that the central density (inside the core radius) of our best-fit model in physical units is 4.4 × 10^3 pc^-3.

To calculate the velocity dispersion profile, we use the method of <cit.> to find the velocity dispersion within a set of radial bins. We compare to two observational profiles (shown here with 2σ error bars): <cit.>, which calculates the velocity dispersion from observations of giants using radial velocity data from <cit.>, and the  survey <cit.>, which uses observations of main sequence stars and giants to calculate the dispersion. To compare to <cit.> (middle panel of Figure <ref>), we include only giants and group into radial bins of 25 stars each. For comparison to MUSE, we include all main-sequence and giant stars with M ≥ 0.4 M_⊙ and choose radial bins such that each bin contains at least 100 stars and covered an annulus of log(r/arcsec) ≥ 0.2, as in <cit.>.

Because the positions and velocities of stars in our models can change significantly from snapshot to snapshot, and as a consequence of the low number of stars per bin, the calculated value of σ_v will fluctuate between snapshots, particularly at small r. Nevertheless, as Figure <ref> shows, the velocity dispersion profile of our best-fit model falls within the 2σ error bars of the observations.

Uncertainty in several observational features of NGC 3201 at present, in particular the cluster's age and distance from Earth, give us freedom to explore ranges in these parameters in order to find a model that accurately matches the observational data of NGC 3201, including the surface brightness and velocity dispersion profiles, as well as r_c and r_hl (calculated as described in Section <ref>).

The uncertainty in the age of NGC 3201 <cit.> allows us to treat all cluster snapshots with ages in range 10 Gyr <t< 12 Gyr as equally valid representations of NGC 3201 at present. Generally, the total cluster mass and total number of retained BHs change by only ∼ 5-10% within this time range, however, as shown in Figure <ref>, r_c and r_hl can fluctuate significantly between successive snapshots. As a result, by exploring different time snapshots we can identify models with r_c and r_hl most similar to NGC 3201.

Because the total cluster mass changes by only small amounts in the time range 10–12 Gyr, the outer parts of the surface brightness and velocity dispersion profiles (r ≳ 10^2 arcsec), which are most sensitive to the total cluster mass, do not change significantly between snapshots. However, in the innermost parts (r ≲ 10^2 arcsec), both the surface brightness and velocity dispersion profiles fluctuate between snapshots because of the increased error due to low N in these regions. 

Secondly, the uncertainty in the distance to NGC 3201 <cit.>, allows the surface brightness profiles for our models to be shifted by small amounts horizontally in either direction which alter the fit to the observational data. For our best-fit model, we adopt a distance of 5.7 kpc, within 2σ of the published distance to the cluster.

Figure <ref> shows how the calculated surface brightness profile for model 5 (our best-fit model) changes based on our choice of cluster age and distance from Earth. Here, the blue plots show surface brightness profiles for all snapshots with 10 Gyr <t< 12 Gyr assuming a distance of 4.3 kpc. Orange plots show profiles for all snapshots assuming a distance of 5 kpc. Red shows the same for a distance of 5.7 kpc, which is distance adopted for our best-fit model, which is shown in the figure as the black curve.

As the figure shows, at small r, changes in snapshot time are the dominant effect in the variations in the surface brightness profile, while at large r, the choice of distance is the dominant effect.





§ BINARY FRACTION


As Figure <ref> demonstrates, with all other initial 
cluster parameters fixed, the scaling of BH natal kicks, and, as a consequence, the number of BHs retained, has a substantial effect on the observational characteristics of a GC model at late times. However, it is unclear from the grid of models used in this study whether variation of other initial parameters may also have similar effects. In particular, previous studies have shown that binaries may provide a sufficient energy source to prevent deep core-collapse, and in absence of sufficient energy production from BH dynamics, energy can be extracted via binary-hardening leading to a quasi-stedy core the size of which depends on the binary fraction <cit.>.

For the grid of models in this analysis, we choose an initial binary fraction of 5%, consistent with the observed binary fraction of NGC 3201 <cit.>. However, we also ran a single model identical to model 16 in Table <ref> (with σ_BH=σ_NS), but with an initial binary fraction of 50%. At t=12 Gyr, this model retains few BHs (N_BH=3) and is core-collapsed. This suggests that even with a binary fraction significantly higher than that suggested by observations, BHs are likely necessary to reproduce the observational features of NGC 3201.




§ VIRIAL RADIUS


In this analysis, we assume an initial virial radius of 1 pc, but we also ran a second grid of models identical to those listed in Table <ref>, but with r_v=2 pc. Initially larger clusters are less dynamically evolved at a fixed present-day age. Thus, in general, these models retain ∼ 10% more BHs than their r_v=1 pc counterparts, but the trend shown in Figure <ref> still holds: models which retain few BHs at late times are core-collapsed, while models with many BHs have observational features similar to NGC 3201. Taking advantage of the uncertainty in distance and age for NGC 3201, as described in Section <ref>, we can similarly identify a best-fit model from the r_v=2pc grid, and draw similar conclusions to those drawn in this analysis. In particular, an NGC 3201-like model must retain several hundred BHs, and such a model can naturally produce BH–MS binaries similar to the BH candidate recently observed in NGC 3201 composed of a low-mass BH and a MS star close in mass to the turn-off.


