A Multimass Velocity Dispersion Model of 47 Tucanae Indicates No Evidence for an Intermediate-mass Black Hole

Our proper motion data come from two epochs of Hubble Space Telescope (HST) observations. HST program GO-12971 (PI: Richer) observed 47 Tuc over 10 orbits early in 2013, imaging the core in the F225W and F336W filters using Wide Field Camera 3. These data were used in conjunction with those of an earlier epoch; GO-9443 (PI: King) imaged the core in 2002 using the ACS High-Resolution Channel in the F475W filter. We match the stars between the F336W and F475W images, providing a baseline of ∼11 yr over which to calculate proper motions.

It is worth taking a moment to note the benefits of using such blue/UV filters. Historically, observations in the core of 47 Tuc have tended to be at longer wavelengths. In such filters, visibility in the central region is hampered both by crowding and by the dominating brightness of the many red giant stars located there (Sarajedini et al. 2007). These giants saturate a portion of the detector in any moderately deep image, bleeding outwards and reducing the ability to accurately measure their stellar neighbors. Our short-wavelength filters act to suppress the light from the cool red giants and allow for an unimpeded deep exposure right in the cluster's core. The shorter wavelengths also improve the diffraction limit of every star in the frame, reducing the degree of crowding. Figure 1 displays the higher core completeness of the bluer filters. For a detailed description of the observations and completeness correction process, see Goldsbury et al. (2016).

Zoom In Zoom Out Reset image size

As a result of this extremely clear view into the core, we achieve quality proper motion measurements of over 50,000 stars, including 12 within the central arcsecond, and 100 within the central 3''. Measurements this close to the cluster core have been extremely difficult in the past, and yet these are exactly the stars that probe the region of greatest potential IMBH influence. A dynamical estimate (Peebles 1972) sets the BH's radius of influence (r_I) to be where a circular orbital speed around the IMBH would equal the central cluster velocity dispersion:

$$\begin{eqnarray}&&{r}_{I}=\displaystyle \frac{{{GM}}_{\mathrm{IMBH}}}{{\sigma }_{c}^{2}},\end{eqnarray} \tag{ 1 }$$

where σ_c is the central velocity dispersion of the cluster. Adopting a distance to 47 Tuc of 4.69 kpc (Woodley et al. 2012) and measuring the central velocity dispersion to be σ_c ≈ 15 km s⁻¹ (value found in this study), we calculate the sphere of influence of a 2000 M_⊙ BH in 47 Tuc to have a radius of r_I ≈ 17. McLaughlin et al. (2006) report a line-of-sight velocity dispersion in the core somewhat lower than that found here, though scaling our velocities to their choice of distance (∼4 kpc) largely explains the difference. Additionally, their line-of-sight dispersion value is for the entire R < 25'' region rather than right in the core. Uncertainties in distance as well as discrepancies in σ_c measurements will affect r_I, but it remains clear that the central few arcseconds are where an IMBH has its most pronounced effect, so visibility in this region is critical.

Throughout this paper, we will be referring to different components and subpopulations of 47 Tuc with the following abbreviations in the text and quantity subscripts: IMBH, low-mass cluster objects (Cl), binaries (bin), white dwarfs (WDs), NSs, and sBHs.

When looking for a velocity dispersion signature in the core, it is important to separately account for populations of dynamically massive objects. These could be inherently massive objects like NSs or sBHs, or they could be binary systems of lighter objects, whose combined mass causes them to behave like heavy objects in the cluster. Relaxation effects cause more massive objects to be concentrated in the core while the less massive ones diffuse toward the outer reaches. This segregation will prove to be important in the context of velocity dispersion and will be discussed in later sections. For the purpose of our analysis, we make the distinction between high-mass (≳1.0 M_⊙) and low-mass (≲1.0 M_⊙) cluster objects. For the low-mass category, we group together the low-mass binaries, low-mass WDs, main-sequence stars, and giants as a single composite population. In the high-mass category, we individually model the populations of heavy binary systems, heavy WDs, NSs, and sBHs. These higher mass objects are more concentrated than the typical cluster object and have the potential to affect the velocity dispersion deep into the core, mimicking the effects of an IMBH (Zocchi et al. 2019).

2.2.1. Binary Population

We make an observationally motivated estimate of the binary population in 47 Tuc. On a color–magnitude diagram (CMD), the binary sequence lies above the main sequence (brighter and redder), and the extent of this offset depends on the mass ratio (q) of the binary pair. A system with a mass ratio of q = 1 has two identical stars and simply twice the flux in every filter, leading to a vertical rise of ∼0.75 mag and no shift in color. Systems with q < 1 will be brightened by a lesser amount and reddened to a small degree depending on the specific q value. While the UV filters mentioned in Section 2.1 provide excellent positional information and allow for high-quality proper motion determination, they unfortunately are not ideal for identifying binary systems. In the UV CMD, the binary population is largely blended with the broad scatter of the main sequence, likely due to multiple populations with differing compositions (Richer et al. 2013). Only those binaries with a mass ratio close to unity stand significantly apart from the main sequence, and even those have substantial contamination. Figure 2 shows a CMD of 47 Tuc from another data set using the redder filters F606W and F814W (Sarajedini et al. 2007, GO-10775). These optical/NIR data cover a similar field of view (central 100'') and provide a clearer distinction between the main sequence and the binary sequence, allowing us to distinguish systems down to q ≈ 0.5. The separation of the binaries from the main sequence occurs because the choice of filters can have a large impact on the morphology of the CMD. Both a small scatter and a shallower slope of the main sequence lead to a more distinct separation of binary systems from the main sequence. For the low-q binaries that remain hidden within the spread of the main sequence (q ≲ 0.5), the large majority of them have total binary masses below 1 M_⊙. Their masses will cause them to follow the distributions of similar-mass main-sequence stars and are thus left to be included as part of the general low-mass population.

Zoom In Zoom Out Reset image size

For a more refined characterization of the binaries, we look at the binary fraction in three different q ranges. We find that for the ranges 0.5 < q < 0.7, 0.7 < q < 0.9, and q > 0.9, the binary fractions equal 0.0265, 0.0098, and 0.0039, respectively. These are displayed in Figure 2, and the specific values are used in computing the binary mass distribution in Section 3.3.

The isochrones used to calculate stellar masses were built from MESA stellar evolution models (Modules for Experiments in Stellar Astrophysics; Paxton et al. 2011, 2013, 2015) with a metallicity z = 0.003 and age 11 Gyr. The fluxes in each filter were determined using Phoenix atmospheres⁴ by Baraffe et al. (2015) and Allard (2016).

2.2.2. Stellar Remnants

In addition to the concentrated binary systems, 47 Tuc harbors a collection of dark stellar remnants in the form of heavy WDs, NSs, and sBHs. Most of these objects are virtually impossible to directly detect, but we can infer their parameters based on an initial mass function (IMF) for the cluster.

In Figure 3, we have plotted the observed stellar-mass function of the visual/NIR data as blue points. Overlaid as a dashed red line is a Kroupa (2001) IMF normalized to the heaviest observed objects. The Kroupa IMF follows a broken power law, N ∝ M^−α, where α = 0.3, 1.3, and 2.3 in the mass ranges M < 0.08, 0.08 < M < 0.5, and M > 0.5 M_⊙, respectively. The deficit of low-mass stars between the data and IMF is presumably due to the preferential loss of low-mass stars from our field of view (segregation) and from the cluster as a whole (evaporation). Baumgardt & Sollima (2017) verified this strong depletion of low-mass stars, and the present-day mass function agrees well with an evolved Kroupa IMF (Baumgardt & Hilker 2018). We make the assumption that the highest mass stars still observed in the cluster (∼0.85 M_⊙) remain at their natal abundance and use this IMF to estimate how many progenitor objects formed at different masses. This is, in several ways, a conservative approach to estimating the number of stellar remnants. The Kroupa IMF is among the steepest mass functions, and assuming some fraction of the 0.85 M_⊙ stars have been lost would only serve to raise the IMF and increase the remnant estimates. By keeping our estimates conservative the final predicted IMBH mass is likely to be an overestimate, as there is potentially more remnant mass that we have not accounted for.

Zoom In Zoom Out Reset image size

As mentioned previously, we are interested in tallying objects that have a mass above 1.0 M_⊙ and would therefore be distributed more centrally than main-sequence stars. The red region of Figure 3 marks progenitor stars that would form WDs of <1.0 M_⊙. These WDs would relax into a distribution comparable to the main-sequence stars. The blue region marks progenitors that produced WDs of mass 1.0 M_⊙ < M < 1.4 M_⊙ (Cummings et al. 2016). The green region indicates the more massive progenitors that form NSs, and the purple region marks the progenitors that produce sBHs of up to 10 M_⊙ (Spera & Mapelli 2017). The number of objects in each colored mass range is determined by integrating the mass function. We cut off the IMF beyond M = 50 M_⊙ because the BHs above 10 M_⊙ tend to be completely lost during cluster evolution (Morscher et al. 2015).

We have to consider that natal kicks imparted upon formation as well as gravitational encounters within the cluster will cause the loss of NSs and sBHs. The retention fraction in globular clusters is not a well-constrained quantity. Simulations tend to produce values anywhere in the range of about 5%–50% (Mackey et al. 2008; Moody & Sigurdsson 2009; Morscher et al. 2015; Baumgardt & Sollima 2017). We adopt a value of 8.5% retention for NSs and sBHs. This choice will be explained in Section 4.3, where we also explore a range of retention values and their effect on the velocity dispersion profile. Building the mass distributions of these heavy remnant objects is explained in the context of the velocity dispersion model in Section 3.4.

In order to determine whether or not our observations are consistent with a central IMBH, we create a theoretical velocity dispersion model that includes the contributions of all the non-IMBH cluster objects as well as an IMBH mass parameter that affects the velocity dispersion signal in the core. We use a 3D version of the King density template (King 1962) as a starting point in determining the shape of the velocity dispersion profile. The density distribution,

$$\begin{eqnarray}&&\rho (r)=K\left[\displaystyle \frac{1}{{\left({r}^{2}+{a}^{2}\right)}^{3/2}}-\displaystyle \frac{1}{{\left({r}_{{\rm{t}}}^{2}+{a}^{2}\right)}^{3/2}}\right],(r\leqslant {r}_{{\rm{t}}}),\end{eqnarray} \tag{ 2 }$$

describes how the density (ρ) falls off with 3D radius (r) as we move away from the core. The density drops to zero at a prescribed tidal radius (r_t), which is set by the Galactic potential. We adopt a fixed tidal radius of r_t = 42' (Harris 1996, 2010 edition) for all populations discussed in this paper. The profile scales with K, a constant that encodes the total mass contained in the distribution, and the parameter a describes an effective core radius. For radii interior to a, the distribution is largely flat.

The density profile and the associated enclosed mass as a function of radius, M(r), can be used in the isotropic Jeans equation,

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dr}}(\rho {\sigma }^{2})=-\rho \displaystyle \frac{d{\rm{\Phi }}}{{dr}}=-\rho \displaystyle \frac{G}{{r}^{2}}M(r),\end{eqnarray} \tag{ 3 }$$

to isolate the square of the velocity dispersion σ²(r). Here, Φ and G are the gravitational potential and constant, respectively. We use the nonrotating and isotropic version of the Jeans equation as the degree of rotation and anisotropy is very low in the core region of 47 Tuc where our proper motion data occur (Watkins et al. 2015; Bellini et al. 2017; Heyl et al. 2017). The squared velocity dispersion then takes the form

$$\begin{eqnarray}&&{\sigma }^{2}(r)=-\displaystyle \frac{G}{\rho (r)}{\int }_{0}^{r}\displaystyle \frac{\rho (r)M(r)}{{r}^{2}}{dr}.\end{eqnarray} \tag{ 4 }$$

This, of course, is the velocity dispersion profile for the 3D radius r. Once σ²(r) is obtained, it is integrated against the density along the line of sight to produce the projected squared velocity dispersion,

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{p}^{2}(R) & = & \displaystyle \frac{{\displaystyle \int }_{R}^{\infty }{\sigma }^{2}(r)\rho (r)r{\left({r}^{2}-{R}^{2}\right)}^{-1/2}{dr}}{{\displaystyle \int }_{R}^{\infty }\rho (r)r{\left({r}^{2}-{R}^{2}\right)}^{-1/2}{dr}}\\ & = & {M}_{{\rm{T}}}F(R),\end{array}\end{eqnarray} \tag{ 5 }$$

where R is the on-sky projected radius from the center. In the second line, we have simply parceled the equation into two parts: a function that holds the R dependence and a scaling term out front. F(R) will have units of (velocity² mass⁻¹) and describes how the velocity dispersion falls off with R. The M_T term is the total mass of the distribution in question and scales the magnitude of the velocity dispersion profile. It arises from M(r) in Equation (4), which involves an integral over the density profile and a coefficient that can be written in terms of M_T. We adopt this format for ease in writing separate component contributions later on in Section 3.4.

The operations between Equations (3) and (5) are fully generic manipulations of the Jeans equation and make no specifications about the density profile. In our full density distribution, we combine different King templates with individual total masses and core radii to capture the unique dynamical contributions of the various cluster populations. For an analogous procedure using Gaussian distributions and with detailed formalism, see Emsellem et al. (1994) and Cappellari (2008). Before we apply the Jeans equation directly to our cluster model, there are several nuances to consider which we describe in the following subsections.

The parameter a in Equation (2) acts as a characteristic core radius. However, when looking at different mass groups, there is not a single global value for a because objects of different masses will be distributed with different core radii. It makes more sense to describe populations of different masses using individual a values. To correct for the inevitability that our tracer stars (those for which we have proper motions) are not perfectly representative of the true mass distribution, we make a distinction between two different core radius parameters, a_Cl and a_PM.

In the right-hand side of Equation (3), the factor −dΦ/dr = −GM(r)/r² generates the force acting on each particle. The M(r) in this part of the equation must therefore reflect the true underlying mass distribution of the cluster, regardless of whether or not those objects are detected in the photometry. Calculating this M(r) term thus requires using a core radius parameter that encompasses all masses, even objects below our detection threshold. This is our a_Cl. It is an average core radius for all components of the cluster and is left as a free fitting parameter given its role of encompassing both detected and undetected mass.

Conversely, the ρ terms of Equation (3) describe the distribution of tracer particles whose motions we are observing. These tracers need not be the same objects that are generating the full potential. We therefore calculate the remaining ρ(r) terms (in Equations (3)–(5)) using the measured core radius of the stars for which we have proper motions, a_PM. This a_PM is a fixed quantity with a value measured by conducting a likelihood fit of the tracer particle positions to the King density template of Equation (3). In doing so, we find a_PM = 360.

Because our tracer stars (even with completeness corrections) are likely somewhat biased toward higher mass stars, we expect a difference between a_Cl and a_PM. This formulation allows a_Cl and a_PM the freedom to differ.

While a_PM is calculated by fitting stellar positions, the free parameter a_Cl is fit using the velocity information of the tracer stars. This is described in more detail in later sections.

The binary and stellar remnant populations that we discuss in the next subsections will also have their own distribution parameters. Due to the differing concentrations, the functional form of the velocity dispersion in Equation (5) now becomes $F(R| a)$, where a is specific to the mass distribution in question.

If there is no IMBH in the mass distribution, we expect the velocity dispersion to flatten off in the core as $M(r)\to 0$. However, a central BH point mass causes M(r) to remain positive for any r > 0 and produces a rising slope into the cluster core. A concentrated population of objects in the cluster can also cause the velocity dispersion profile to rise farther into the core before flattening out, mimicking the central effects of a BH. When testing for an IMBH, it is therefore important to take careful consideration of populations more centrally concentrated than the typical main-sequence star, which may cause the velocity dispersion to rise deeper into the core.

For this reason, we model the binary population separately from the general lower mass cluster stars. However, the binary systems are not a single homogeneous population in terms of their mass. For a given primary mass (i.e., position along the main sequence), systems with high mass ratios are more massive than those with low q-values, and for a given mass ratio, binaries at the top of the main sequence are more massive than those below them. The higher mass systems will have correspondingly small distribution parameters. We therefore break the binary contribution into components to account for binary subpopulations with different distributions. The expression for the squared velocity dispersion of the binary stars becomes

$$\begin{eqnarray}&&{\sigma }_{p}^{2}(R)=\displaystyle \sum _{j,\alpha }{M}_{\mathrm{bin}}^{j,\alpha }F(R| {a}_{\mathrm{bin}}^{j,\alpha }).\end{eqnarray} \tag{ 6 }$$

The index j signifies different mass bins down the main sequence in the range of 0.85–0.55 M_⊙. The index α specifies one of the three mass ratio ranges seen in Figure 2: lower (0.5 < q < 0.7), moderate (0.7 < q < 0.9), or high (q > 0.9). The mass of binaries in each j-bin and α-range (${M}_{\mathrm{bin}}^{j,\alpha }$) is calculated as

$$\begin{eqnarray}&&{M}_{\mathrm{bin}}^{j,\alpha }={f}^{\alpha }\left(\displaystyle \frac{{A}_{\mathrm{MS}}}{{A}_{\mathrm{bin}}}\right)\left(1+\langle q{\rangle }^{\alpha }\right)\left(\displaystyle \frac{{M}_{\mathrm{MS}}^{j}}{{M}_{\mathrm{obs}}}\right){M}_{\mathrm{Cl}}.\end{eqnarray} \tag{ 7 }$$

The first item, f^α, is the observed binary fraction for a given mass ratio range (α). These values are determined by counting stars in a section of the CMD as seen in Figure 2 and Table 2. We assume that these binary fraction values hold true for the entire main sequence that is within our sensitivity limit. The A factors are the proportion of the population that is visible within the 100'' field of view according to the distribution parameter for that population. Therefore, the term A_MS/A_bin modifies the observed binary fraction to describe the entire cluster, including objects beyond the image field. This is necessary, as the binary fraction is expected to drop away from the core because the binary systems are more concentrated than the main-sequence stars. The factor $(1+\langle q{\rangle }^{\alpha })$ accounts for the mass of the secondary star in the binary system where $\langle q{\rangle }^{\alpha }$ is the midpoint of the mass ratio range α. At the end of Equation (7), the quantity ${M}_{{\rm{MS}}}^{j}$ is the observed mass in main-sequence bin j, obtained by counting stars and consulting the isochrone. It is divided by M_obs, the total mass of the observed main sequence (again by consulting the isochrone), to determine the mass proportion of the entire cluster our selection j corresponds to. Finally, multiplying by M_Cl scales this proportion to the fitted dynamical cluster mass.

Table 2.  Model Inputs for Binaries and Stellar Remnants

Object Type	Subgroup Parameters/Type	Primary Mass (Mp) Range	$\langle {M}_{p}\rangle$	Object Mass	Retention Fraction	Population Mass	a
		(M⊙)	(${M}_{\odot })$	(${M}_{\odot })$		(${M}_{\odot }$)	('')
Binaries	q = [0.5, 0.7]	0.80–0.86	0.83	1.33	1	7102	18.61
	($\langle q\rangle =0.60$)	0.74–0.80	0.77	1.23	1	6537	20.07
	(f = 0.00265)	0.67–0.74	0.71	1.13	1	5979	21.72
		0.61–0.67	0.64	1.03	1	5110	23.71
		0.55–0.61	0.58	0.93	1	5011	26.12
	q = [0.7, 0.9]	0.80–0.86	0.83	1.50	1	2824	16.64
	($\langle q\rangle =0.80$)	0.74–0.80	0.77	1.38	1	2592	17.95
	(f = 0.00098)	0.67–0.74	0.71	1.27	1	2363	19.42
		0.61–0.67	0.64	1.16	1	2013	21.20
		0.55–0.61	0.58	1.05	1	1965	23.36
	q = [0.9, 1.0]	0.80–0.86	0.83	1.62	1	1193	15.42
	($\langle q\rangle =0.95$)	0.74–0.80	0.77	1.50	1	1094	16.63
	(f = 0.00039)	0.67–0.74	0.71	1.38	1	995	18.00
		0.61–0.67	0.64	1.26	1	845	19.65
		0.55–0.61	0.58	1.13	1	822	21.65
Remnants	WD	⋯	⋯	1.2	1	24000	20.52
	NS	⋯	⋯	1.4	0.085	23800	17.72
	sBH	⋯	⋯	10	0.085	19000	2.74

Note. The object mass within a binary subgroup is calculated as $\langle M\rangle (1-\langle q\rangle )$. There is a total mass of ∼46,500 M_⊙ in binaries and ∼27,600 M_⊙ in stellar remnants (taking into account the 8.5% retention of NSs and sBHs). Population masses for binaries are calculated as described in Equation (7). Population masses for stellar remnants are calculated by assuming an object mass, determining the distribution parameter (a) from that mass, and then consulting the predicted progenitor counts of Figure 3 that have been corrected to include the entire cluster instead of the limited field of view. Population mass values presented are before retention is considered.

Download table as:  ASCIITypeset image

Because of the shape and the width of the main sequence and giant branch in Figure 2, the binary stars cannot be identified for all masses, and thus, their distributions cannot be measured directly. We make use of the results of Goldsbury et al. (2013) to determine the ${a}_{\mathrm{bin}}^{j,\alpha }$ parameter of Equation (6). In order to quantify the mass segregation in 54 globular clusters (including 47 Tuc), Goldsbury et al. (2013) measured the distributions of various stellar groups, determining a power-law relation between the objects' mass and their distribution parameter. They fit a mass segregation power law of R_c = A(M/M_⊙)^B, finding best-fit values of A = 24.4 ± 0.6 and B = −0.95 ± 0.05 for 47 Tuc. In their paper, R_c works out to be equivalent to our a parameter. We use this power law to determine ${a}_{\mathrm{bin}}^{j,\alpha }$ given that the mass of the binary system is the mass of the primary times the factor $(1+\langle q{\rangle }^{\alpha })$.

The dark stellar remnants are included in a similar but simpler manner. Each of the heavy WD, NS, and sBH populations are given a single distribution specified by their total mass and core radius parameter. Having no direct measurements of the dark remnant populations, we are required to make some simplifying assumptions. The heavy WDs are assumed to have a mean mass of 1.2 M_⊙, the midpoint of the considered mass range (1 M_⊙ to the Chandrasekhar limit 1.4 M_⊙). The NSs are all conservatively estimated to have masses of 1.4 M_⊙, and the sBHs are all assumed to have masses of 10 M_⊙ (Morscher et al. 2015; Cummings et al. 2016; Spera & Mapelli 2017). The populations' overall masses are simply given by the predicted number of objects multiplied by their mass. The number of objects in each of Figure 3's colored mass ranges is determined by integrating the mass function. This value is then corrected to include everything beyond the image field of view. This correction process is analogous to the treatment used for binaries within Equation (7). Using estimates of the remnant masses, we can infer their distributions and include them in our velocity dispersion model using the Goldsbury et al. (2013) mass segregation power law as was done for the binary systems. However, the sBHs and NSs will suffer additional losses, largely due to natal kicks from their progenitor supernovae. We assume a moderate value of 8.5% retention for NSs and sBHs. The projected squared velocity dispersion contribution from each dark stellar remnant populations becomes

$$\begin{eqnarray}&&{\sigma }_{p,k}^{2}(R)={M}_{k}F(R| {a}_{k}),\end{eqnarray} \tag{ 8 }$$

where k is an index that refers to the dark stellar remnant in question (WD, NS, sBH). The full and final expression for the squared velocity dispersion comes together as

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{p}^{2}(R) & = & {M}_{\mathrm{IMBH}}F(R)\\ & & +\ {M}_{\mathrm{Cl}}F(R| {a}_{\mathrm{Cl}})\\ & & +\ \displaystyle \sum _{j,\alpha }{M}_{\mathrm{bin}}^{j,\alpha }F(R| {a}_{\mathrm{bin}}^{j,\alpha })\\ & & +\ \displaystyle \sum _{k}{M}_{k}F(R| {a}_{k}).\end{array}\end{eqnarray} \tag{ 9 }$$

In each case, the F function encapsulates the radial dependence of the velocity dispersion for the component in question and is scaled by the total mass M term of that component. The IMBH term is slightly different in that the enclosed mass function (M(r) of Equations (3) and (4)) is simply the constant M_IMBH because we have considered it a point mass centered on r = 0. It therefore is influenced only by the a_PM value within ρ(r). It has no distribution parameter equivalent to the a_Cl, a_bin, or a_k found in the other terms.

We have chosen to split the heavy objects by their type rather than their mass group simply for ease of estimation. One could also imagine lumping all objects of similar mass regardless of whether they are binary, WD, NS, or sBH. We found it simpler to estimate object counts by considering direct observational evidence of the binary sequence and estimating progenitor counts for the remnants.

We have now built a σ_p²(R) distribution from 47 Tuc's different mass components. Using measured binary fractions and inferred distributions, we account for the individual contributions of binary subpopulations that have differing distributions and total masses. An IMF normalized to observed cluster stars allows us to estimate how many dark stellar remnants remain in the cluster. The fit parameters M_Cl and a_Cl describe the gross shape of the velocity dispersion curve (its vertical scale and the radius at which it falls off), while the concentrated binaries and remnants bring up the velocity dispersion value toward the core. Any central rise observed in the data that is not due to the concentrated populations is taken up by the IMBH fit parameter M_IMBH.

Putting all of this together produces a radial profile for the projected velocity dispersion that depends on three parameters: M_IMBH, M_Cl, and a_Cl. The rest of the inputs are measured or inferred quantities. The projected velocity dispersion profile, ${\sigma }_{p}(R| {M}_{\mathrm{IMBH}},{M}_{\mathrm{Cl}},{a}_{\mathrm{Cl}})$, is now in a form that can be fit to the data and the optimal values for the three fitting parameters found.

With a theoretical velocity dispersion model now built, we are in a position to find the set of the three parameters (M_IMBH, M_Cl, and a_Cl) that best fit the observations. This is most simply and commonly done by radially binning the data, computing the velocity dispersion in each bin, then finding a least-squares fit of the model to the data.

We must, however, be cautious of trusting the results of a binned data set. Any binning scheme has some degree of arbitrariness that may influence the result in unpredictable ways. At the very least, introducing bins will smear out the radius-dependent signal over a larger range of radii. At worst, the choice of binning system may unintentionally amplify or suppress a particular result.

As a test of the influence of bin choice, we set up a logarithmically spaced binning scheme then steadily change the number (and thus the size) of the bins. Figure 4 shows how the binned-fit parameters fluctuate as we slowly scale the bin sizes. We find that the jitter the parameters experience is generally within the typical uncertainty of each parameter. Despite the agreement between binning schemes, we still decide to use an unbinned fitting method for two reasons: an unbinned system is more inclusive of every piece of information available to us, and it also removes any human element of bin choice. The mean values of each of the Figure 2 parameters are slightly offset from the unbinned results reported later in Section 4. We suspect this is due to the smearing-out of the data when employing bins, though the different methods do agree within their respective uncertainties. We use the following unbinned likelihood maximization technique to find the best-fit velocity dispersion model.

Zoom In Zoom Out Reset image size

Under the assumption of isotropy in the core for 47 Tuc, we assume the observed velocity components (v_x, v_y) of each star are samples from Gaussian-distributed populations centered on zero. The widths of these Gaussian populations are estimated by the model velocity dispersion value at each star's projected distance from the core, broadened by the uncertainty in the velocity measurements. We find the set of M_IMBH, M_Cl, and a_Cl parameters for σ_p(R) that maximize the global log-likelihood function,

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{ln}\,L & = & \displaystyle \sum _{i}-\displaystyle \frac{{v}_{x,i}^{2}}{2({\sigma }_{p,i}^{2}+{\delta }_{x,i}^{2})}-\displaystyle \frac{{v}_{y,i}^{2}}{2({\sigma }_{p,i}^{2}+{\delta }_{y,i}^{2})}\\ & & -\ \displaystyle \frac{1}{2}\mathrm{ln}({\sigma }_{p,i}^{2}+{\delta }_{x,i}^{2})-\displaystyle \frac{1}{2}\mathrm{ln}({\sigma }_{p,i}^{2}+{\delta }_{y,i}^{2}).\end{array}\end{eqnarray} \tag{ 10 }$$

Here, ${\sigma }_{p,i}={\sigma }_{p}({R}_{i}| {M}_{\mathrm{IMBH}},{M}_{\mathrm{Cl}},{a}_{\mathrm{Cl}})$ is the model velocity dispersion value at the projected radius of star i, and δ_x,i and δ_y,i are the measurement uncertainties in the velocity components of star i that arise from our photometric data. This likelihood estimator is equivalent to Equation (8) of Walker et al. (2006), but has been adapted for 2D proper motions and a continuous σ_p²(R) rather than 1D radial velocities and binned velocity dispersion values. Furthermore, we ignore the constant offset introduced to the $\mathrm{ln}L$ value by Gaussian normalization. For a full derivation and explanation of the entire unbinned likelihood technique, see Goldsbury et al. (2016). The specific parameters that went into our model are shown in Table 2, and the resulting velocity dispersion profiles are plotted in Figure 6.

We should note that this model is not restricted from fitting a negative M_IMBH parameter and clarify what a negative value indicates. The model cannot measure mass directly, only the amount of dispersion present in the stellar velocities. A negative IMBH mass, while clearly an unphysical phenomenon if taken at face value, simply indicates that the chosen components of the model (e.g., Table 2) produce a higher velocity dispersion in the core than is seen in the data. The fitting process adopts a negative IMBH mass parameter in an attempt to bring down the velocity dispersion in the core to improve the overall fit. This indicates that some imperfect assumptions and estimations were made in building the model and the negative IMBH mass is a sign that the model places more mass in the core than can be explained by the observed velocities.

Having built a velocity dispersion model that includes measured binary fractions, inferred dark stellar remnants, and an NS and sBH retention fraction of 8.5%, we obtain the following best-fit parameter values: M_IMBH = 40 ± 1650 M_⊙, M_Cl = (1.39 ± 0.07) × 10⁶ M_⊙, and a_Cl = 448 ± 16. This provides us with an M_IMBH/M_Cl ratio of 0.00% ± 0.12%. Note that the binaries and remnants themselves sum to ∼5.3% of M_Cl (see Table 2) and are not included in the M_Cl parameter. Therefore, the entire mass of 47 Tuc comes to 1.47 × 10⁶ M_⊙. The uncertainties in the fit parameters were determined through a bootstrapping method in which we randomly sampled the stars (with replacement) to obtain a new data set to which we fit the velocity dispersion model. The spreads in the distributions of the resulting fit parameters provide our estimates of uncertainty (shown in Figure 5).

Zoom In Zoom Out Reset image size

The velocity dispersion model described by these best-fit parameters is shown in Figure 6, plotted over binned velocity dispersion data. The binned data are included for visualization purposes only as the fit was made using the unbinned likelihood maximization described in the previous section. For concentrated heavy populations, the binaries and sBHs play dominant roles in the central velocity dispersion rise. The binary population contains a large overall mass and their distribution causes the velocity dispersion to continue rising into ∼10''. The sBHs, while comprising a lesser total mass, are much more concentrated and continue this velocity dispersion rise into ∼2''. The combined effect leaves no room for any IMBH that would drive the velocity dispersion curve even higher.

Zoom In Zoom Out Reset image size

The results reported here are somewhat influenced by our choice of cluster distance. We have adopted a distance of d = 4.69 kpc from Woodley et al. (2012), who fit spectral energy distributions of 47 Tuc's WDs. The choice of distance scales the observed velocities, raising or lowering the velocity dispersion profile. This has a direct effect on the fitted cluster mass (set by the height of the profile) but has a less straightforward influence on the IMBH mass parameter, due to our inclusion of dark stellar remnants in the model. The remnants were calculated independent of any velocity measurement, so they do not scale with distance. Therefore, if we were to adopt a lower cluster distance (producing lower stellar velocities), the cluster mass would scale downward while the remnant population remains fixed. This would make the contribution of the remnants proportionally more important in the core and further reduce the need for an IMBH influence.

We verified this effect by considering the recent parallax distance of 4.45 kpc for 47 Tuc from Chen et al. (2018). Fitting the model with this lower distance dropped the cluster mass down to M_Cl = 1.20 × 10⁶ M_⊙ and the fitted IMBH mass to a negative value, M_IMBH = −560 M_⊙ (at 8.5% sBH and NS retention). Our calculated stellar remnants already accounted for enough concentrated mass in the core. A smaller cluster distance only serves to push the fitted IMBH mass lower.

Our fitted cluster mass is relatively large. Giersz & Heggie (2011) and Lane et al. (2010) found 0.9 and 1.1 × 10⁶ M_⊙, respectively. The discrepancy could be influenced by the choice of distance as discussed above or potentially be due to our data limitations. Our proper motion data only extend to 120'', a small fraction of the tidal radius, r_t ≈ 42' (Harris 1996, 2010 edition). This means we are really only fitting the central region of the cluster. Our cluster mass estimate involves integrating the entire density distribution and therefore assumes that the density beyond 120'' follows a perfect King template out to the tidal radius. Deviation in the tails may contribute to errors in the cluster mass.

We would like to test our velocity dispersion model to ensure it is accurately representing the data it is being fit to. Being able to find a set of best-fit model parameters does not necessarily mean one has modeled the data well. For example, any set of curved data will have a straight line that fits it best, but that does not mean a straight line is a particularly good description of those data. To determine the quality of our model, we employ a system of statistical resampling in order to assess the efficacy of our modeling process.

We generate new v_x and v_y velocity components for each star. These velocities are drawn randomly from Gaussian distributions whose widths are determined by the best-fit model velocity dispersion value at the radius of each data point and broadened by the uncertainty of the velocity measurement. We keep the original positions for each star to ensure the radial distribution remains unchanged, maintaining a consistent radial sampling in the core where there are few stars. The velocity dispersion model is fit to this new set of artificial data, and the fit parameters as well as a goodness-of-fit measure, the log-likelihood (ln L), are recorded. Repeating this process many times builds distributions of the fit parameters and ln L values. The artificial data sets will clearly produce a Gaussian spread around the input model, but the pertinent question is whether or not the true data produce a similar fit. Figure 7 shows the distributions of the generated fit parameters and ln L values. As can be seen, the orange lines/dots indicating the true values of the data fall very near the middle of each distribution. Such close agreement of the M_IMBH, M_Cl, and a_Cl parameters tell us our model suffers no systematic bias in recovering input parameters. The close agreement of the ln L value indicates our model is a good morphological description of the data. We can conclude that artificial velocity data generated from our model closely resemble the cluster's velocity data.

Zoom In Zoom Out Reset image size

We explored a range of sBH retention rates to determine the effect it can have on the final fitted IMBH mass (see Table 3 and Figure 8). The negative IMBH models (dark curves of Figure 8) can be instructive in signifying inappropriate model assumptions, but they are ultimately unphysical systems. Our adoption of an 8.5% retention fraction arises from choosing the highest likelihood model that produces nonnegative IMBH mass parameter (Table 3). There is a somewhat degenerate trade-off between the influence of tightly distributed massive objects versus a central IMBH. This trade-off happens in the 2''–10'' region where the sBH distribution can produce a slope in the velocity dispersion. Interior to ∼2'', nothing can produce a slope except for a point mass, so the IMBH term always dominates in this region.

Zoom In Zoom Out Reset image size

Table 3.  Exploring Retention Fraction

Retention	MIMBH	MCl	aCl	$\tfrac{{M}_{\mathrm{IMBH}}}{{M}_{\mathrm{Cl}}}$	ln L
(%)	(M⊙)	(M⊙)	('')
0	4125	1346143	42.8	0.0031	−310214.24
5	1713	1374272	43.9	0.0012	−310213.41
8	278	1392757	44.7	0.0002	−310213.01
8.5	40	1395974	44.8	0.0000	−310212.95
10	−672	1405850	45.2	−0.0005	−310212.78
15	−3020	1441640	46.6	−0.0021	−310212.40
20	−5317	1482612	48.1	−0.0036	−310212.32
25	−7545	1530150	49.8	−0.0049	−310212.63

Note. Enforcing various retention fractions for NSs and sBHs provides different fit results. Typical uncertainties for these parameters are very similar to those reported in Section 4.1 (those results are the 8.5% retention case shown here). Retention below ∼8.5% is needed to find a positive fit mass for an IMBH. A retention of essentially 0% is required to obtain the Kızıltan et al. (2017) mass fraction of 0.30%. It appears that retention above ∼8.5% starts to become incompatible with the core velocity dispersion seen in this 47 Tuc data. Corresponding velocity dispersion profiles are shown in Figure 8.

Download table as:  ASCIITypeset image

We observe that in order to find the M_IMBH/M_Cl fraction concluded by Kızıltan et al. (2017), we would need to assume that essentially all 10 M_⊙ sBHs are ejected from the cluster. The only way we can see evidence of an IMBH at all is if we assume an sBH and NS retention fraction very close to zero. Pushing the retention any higher than ∼8.5% creates a larger velocity dispersion in the core than is observed, requiring a negative IMBH mass fit. Interestingly, from the trials shown in Table 3 and Figure 8, the ln L values' quadratic shape peaks near 18.4% retention. This may indicate a model that overestimates the core contribution from remnant/binary populations (requiring negative M_IMBH) but fits the numerous data points at intermediate radii very well. We clearly see the trade-off between sBH populations and an IMBH that can produce similar dynamical effects.

Our use of UV imaging data has allowed us to glean velocity information from stars right into the very center of the globular cluster 47 Tucanae. Historically, crowding has been a limiting factor in IMBH velocity dispersion searches as it severely reduces the numbers of stars measurable at a very small projected radius. Resolution and accuracy in this region are of key importance in trying to estimate the mass of a central IMBH through this method. With visibility right in to the cluster core, we are able to probe the region of the cluster where an IMBH's influence would be most pronounced.

Beginning with a 3D King density template, we used the isotropic Jeans equation to build a velocity dispersion profile with the added influence of a central IMBH. The velocity dispersion profile was built with contributions from several subcomponents: a central IMBH point mass, a collective group of low-mass cluster objects (≲1.0 M_⊙), and individual contributions from the concentrated populations of the binary systems, heavy WDs, NSs, and sBHs.

We employed an unbinned likelihood analysis when fitting our model. This technique weighs the combined likelihoods of each individual star's velocities being drawn from a given velocity dispersion model. It prevents the radial smearing-out that occurs when data are binned together and removes any arbitrary human factor in choosing a binning scheme.

To determine how good a fit our velocity dispersion model is to the data, we randomly resampled the velocity components of our stars based on the best-fit velocity dispersion curve. The log-likelihood of the real data's best-fit model is in very good agreement with the distribution of resampled log-likelihood values, indicating our model does a good job of describing the data.

Our analysis produces a best-fit value of M_IMBH = 40 ± 1650 M_⊙ for a central IMBH in the core of 47 Tucanae, giving an M_IMBH/M_Cl ratio of 0.00% ± 0.12%. We find no evidence to suspect an IMBH is present in the core of 47 Tuc unless one also accepts a very low sBH and NS retention fraction. Our result is consistent with zero and in contention with the results of Kızıltan et al. (2017), who found an M_IMBH/M_Cl ratio of $0.30{ \% }_{-0.12 \% }^{+0.20 \% }$. We emphasize that the proper characterization of the binary population and sBHs can have a marked effect on the IMBH estimate and their omission will likely bias that estimate high.

The velocity dispersion profile investigated in this study does not support an IMBH detection for 47 Tuc unless one also posits a very low remnant retention fraction. The concentrated populations of binary stars and dark stellar remnants alone appear to be enough to explain the velocity dispersion in the core.

Though the main thrust of this paper is the analysis of 47 Tuc's stellar velocity dispersion, we can also briefly look at complementary avenues of analysis. In the current state of the field, a firm case for an IMBH detection (or absence) is likely to be drawn from a suite of findings rather than any single result.

There is numerical evidence that the presence of an IMBH may cause a weak central cusp in the cluster's density/surface brightness profile (Baumgardt et al. 2005; Noyola & Baumgardt 2011); however, other works have shown that an observed cusp may also be caused by a state change in the degree of core collapse (Trenti et al. 2010; Vesperini & Trenti 2010) and is not sufficient on its own to infer an IMBH. Our UV data show no evidence for a central cusp in the density distribution. Figure 9 shows the cumulative distribution of objects (corrected for incompleteness) as a function of projected distance from the core. This distribution is compared against an R² relation, which would indicate a perfectly flat surface density profile. The mild deviation from R² in the center is very small and is not significant, as it is at the low end of the logarithmic scale and comprises very few objects (≲10). P-values from the Kolmogorov–Smirnov test shown in the bottom panel of Figure 9 indicate that we cannot reject the null hypothesis (i.e., that the data are drawn from an R² distribution) to any reasonable significance until we look beyond ∼20'', where the cluster density begins to fall off.

Zoom In Zoom Out Reset image size

In addition, we find no sign of a small population of high-velocity stars. If the IMBH were to become part of a binary system, three-body interactions should periodically eject high-velocity stars into the cluster on largely radial orbits. Our data do not show any stars with velocities beyond ∼60 km s⁻¹, while N-body simulations by H. Baumgardt, which are filtered to match our observational limitations, show a small population of stars with velocities up to ∼80 km s⁻¹ when they include a central IMBH.

Though we find conflicting results with the Kızıltan et al. (2017) study, contrary findings are nothing new in this field. Of particular note in recent years are the contradictory results of of Lützgendorf et al. (2015) and Lanzoni et al. (2013) on NGC 6388, van der Marel & Anderson (2010) and Noyola et al. (2008) on ω Centauri, and Perera et al. (2017) and Gieles et al. (2018) on NGC 6624. There is clearly more work to be done in characterizing and evaluating the shortcomings and strengths of different IMBH detection methods, and we could benefit from more studies that evaluate the accuracy and biases of these techniques. For example, de Vita et al. (2017) provided insight into how using integrated-light spectroscopy might affect results versus using velocity measurements of individual objects. Simulations can be extremely useful with their ability to explicitly incorporate a central IMBH of known mass to determine its observable effects, but the expensive computation limitations require simplified physics and assumptions about initial conditions and short-timescale dynamics. As has been shown in this study, some of the short-timescale phenomena such as stellar evolution/remnant retention and binary formation can play a very key role in IMBH determination.

Consideration of the sBH retention fraction led us to an unanticipated additional result of this study. Even employing a conservative IMF to estimate the number of sBHs, we still find that a retention of ≳8.5% becomes incompatible with our velocity dispersion observations. Above this fraction, our model fits a negative IMBH mass in an attempt to drop the velocity dispersion values in the core. This clearly unphysical result indicates that the model's inferred massive populations provide more concentrated mass in the core than can be explained by the the observed velocity dispersion. For 47 Tuc, it seems that we can place the constraint on the sBH retention to be ≲8.5%.

For the sake of clarity and for the guidance of future research, we would like to acknowledge a few of the areas where this study could benefit from improvement.

Including additional data sets for regions beyond the central ∼100'' would likely provide tighter constraints on the M_Cl and a_Cl parameters. Data tight in the core are important to constrain M_IMBH, but the other parameters are sensitive to the entire profile beyond our field of view. Data sets would have to be included in a careful manner to appropriately account for their differing spatial and wavelength coverages, as well as their individual completeness limitations.

Line-of-sight velocities in the core could also help to boost star counts where there are few, though these types of data have their own distinct biases and limitations to consider.

Our determination of binary counts and masses relied heavily on the use of isochrones. Stellar evolution modeling has come a long way but is not without its uncertainties. The effects of using different isochrones and decontaminating the binary sequence are topics that have bearing on our results but merit full studies in their own right.

There are a few aspects specific to our model that were approximated to some degree or perhaps not expanded to their fullest potential. The sums of the King templates that we used for the density distributions in the Jeans equation can provide an internally consistent model (i.e., a distribution function f(E) that satisfies Jeans' theorem) for a spherical cluster with isotropic velocity dispersion through Eddington's formula (Binney & Tremaine 2008) for a vanishingly small BH mass if the values of a do not differ by too large of a factor. Because there is no evidence for anisotropy or rotation in the core region (Watkins et al. 2015; Bellini et al. 2017; Heyl et al. 2017), an isotropic distribution function f(E) provides a good approximation. Rather than constructing f(E) explicitly, we compare the values of the velocity dispersion as calculated from the model using the Jeans equation with those of the data. In fitting the data, we assume a normal distribution for the velocities near the center of the cluster (Equation (10)) when fitting for the observed velocity dispersion; this final assumption is not consistent with the King template and Jeans theorem, but it does provide an excellent approximation, given the observational data. Furthermore, with the inclusion of the central BH, we can construct an isotropic distribution function using Eddington's formula everywhere except in the immediate vicinity of the BH (about 10'' for the largest BH mass considered), so our isotropic models do not satisfy Jeans theorem there. A more complete model could use the observed isotropy and the lack of a density cusp to provide an additional constraint on the BH mass. In the absence of an IMBH, we find that the concentration of the retained sBHs is sufficient to render our assumed density profile and isotropy inconsistent with Jeans theorem within about 7'' of the center. For future work, a more detailed model could be constrained with data at larger radii and with the lack of anisotropy or a density cusp at small radii. One might also choose to set up the model to simultaneously fit the cluster's density and surface brightness profile along with the velocity dispersion for a built-in self-consistency check.

Some of the details comprising our chosen model are worth consideration as well. For instance, we rely on the Goldsbury et al. (2016) mass segregation power law to infer population concentrations for unseen objects based on their mass. Goldsbury et al. (2016) only tested this relation up to ∼0.8 M_⊙, so we have no guarantee that objects above this mass regime will follow the same relation.

Additionally, one might get a better picture of the dark stellar remnant populations by doing a careful IMF evolution. Such a process may provide more accurate masses and counts for these dark objects.

Where possible, we attempted to make conservative estimates of the concentrated mass components. The fact that our model predicts no IMBH despite these conservative assumptions degrades any support for the idea of an IMBH in 47 Tuc according to this analysis.

We have shown here that proper motion velocity dispersion analysis can prove to be a useful tool in the search for IMBHs, even in the crowded environments of a globular cluster core. With the right choice of photometric filters, we can overcome the crowding issues and glean a multitude of data in the region most sensitive to an IMBH's influence. The inclusion of concentrated populations in our mass model is sufficient to explain the observed velocity dispersion in 47 Tuc's core. That being said, the velocity dispersion signature of an IMBH and a tight population of sBHs is very similar. A trade-off in mass between them would not be very evident in the velocity data. Independent estimates of the sBH retention and progenitor IMF are likely the best ways to disentangle this degeneracy.